3
# pp_setversion $VERSION; # haven't worked out why it breaks my system (CS)
4
pp_beginwrap; # required for overload to work
6
# pp_def functions go into the PDL::Complex namespace
7
# to avoid clashing with PDL::FFTW funcs of the same name that go
8
# into the PDL namespace
9
# it should be of no effect to the user of the module but you
11
pp_bless('PDL::Complex');
13
pp_addpm {At => Top}, <<'EOD';
18
use vars qw($sep $sep2);
22
pp_addpm {At => Top}, <<'EOD';
25
PDL::Complex - handle complex numbers
34
This module features a growing number of functions manipulating complex
35
numbers. These are usually represented as a pair C<[ real imag ]> or
36
C<[ angle phase ]>. If not explicitly mentioned, the functions can work
37
inplace (not yet implemented!!!) and require rectangular form.
39
While there is a procedural interface available (C<< $a/$b*$c <=> Cmul
40
(Cdiv $a, $b), $c) >>), you can also opt to cast your pdl's into the
41
C<PDL::Complex> datatype, which works just like your normal piddles, but
42
with all the normal perl operators overloaded.
44
The latter means that C<sin($a) + $b/$c> will be evaluated using the
45
normal rules of complex numbers, while other pdl functions (like C<max>)
46
just treat the piddle as a real-valued piddle with a lowest dimension of
47
size 2, so C<max> will return the maximum of all real and imaginary parts,
48
not the "highest" (for some definition)
50
=head1 TIPS, TRICKS & CAVEATS
56
C<i> is a constant exported by this module, which represents
57
C<-1**0.5>, i.e. the imaginary unit. it can be used to quickly and
58
conviniently write complex constants like this: C<4+3*i>.
62
Use C<r2C(real-values)> to convert from real to complex, as in C<$r
63
= Cpow $cplx, r2C 2>. The overloaded operators automatically do that for
64
you, all the other functions, do not. So C<Croots 1, 5> will return all
65
the fifths roots of 1+1*i (due to threading).
69
use C<cplx(real-valued-piddle)> to cast from normal piddles into the
70
complex datatype. Use C<real(complex-valued-piddle)> to cast back. This
71
requires a copy, though.
75
This module has received some testing by Vanuxem Gr�gory
76
(g.vanuxem at wanadoo dot fr). Please report any other errors you
81
=head1 EXAMPLE WALK-THROUGH
83
The complex constant five is equal to C<pdl(1,0)>:
88
Now calculate the three cubic roots of of five:
90
pdl> p $r = Croots $x, 3
91
[1.70998 +0i -0.854988 +1.48088i -0.854988 -1.48088i]
93
Check that these really are the roots:
96
[5 +0i 5 -1.22465e-15i 5 -7.65714e-15i]
98
Duh! Could be better. Now try by multiplying C<$r> three times with itself:
101
[5 +0i 5 -4.72647e-15i 5 -7.53694e-15i]
103
Well... maybe C<Cpow> (which is used by the C<**> operator) isn't as
104
bad as I thought. Now multiply by C<i> and negate, which is just a very
105
expensive way of swapping real and imaginary parts.
108
[0 -1.70998i 1.48088 +0.854988i -1.48088 +0.854988i]
110
Now plot the magnitude of (part of) the complex sine. First generate the
113
pdl> $sin = i * zeroes(50)->xlinvals(2,4) + zeroes(50)->xlinvals(0,7)
115
Now plot the imaginary part, the real part and the magnitude of the sine
116
into the same diagram:
118
pdl> use PDL::Graphics::Gnuplot
119
pdl> gplot( with => 'lines',
120
PDL::cat(im ( sin $sin ),
124
An ASCII version of this plot looks like this:
126
30 ++-----+------+------+------+------+------+------+------+------+-----++
127
+ + + + + + + + + + +
145
5 ++ $$$############ * # ++
146
|*****$$$### ### * # |
160
+ + + + + + + ### + ### + + +
161
-15 ++-----+------+------+------+------+------+-----###-----+------+-----++
162
0 5 10 15 20 25 30 35 40 45 50
167
BEGIN { $i = bless pdl 0,1 }
168
sub i () { $i->copy };
171
for (qw(Ctan Catan re im i cplx real)) {
172
pp_add_exported '', $_;
180
# define M_PI 3.1415926535897932384626433832795029
183
# define M_2PI (2. * M_PI)
186
#if __GLIBC__ > 1 && (defined __USE_MISC || defined __USE_XOPEN || defined __USE_ISOC9X)
187
# define CABS(r,i) hypot (r, i)
190
CABS (double r, double i)
206
return r * sqrt (1 + t*t);
210
#if __GLIBC__ >= 2 && __GLIBC_MINOR__ >= 1 && defined __USE_GNU
211
# define SINCOS(x,s,c) sincos ((x), &(s), &(c))
213
# define SINCOS(x,s,c) \
219
#define CSQRT(type,ar,ai,cr,ci) \
220
type mag = CABS ((ar), (ai)); \
227
t = sqrt (0.5 * (mag + (ar))); \
229
(ci) = 0.5 * (ai) / t; \
233
t = sqrt (0.5 * (mag - (ar))); \
238
(cr) = 0.5 * (ai) / t; \
243
#define CLOG(ar,ai,cr,ci) \
244
(cr) = log (CABS ((ar), (ai))); \
245
(ci) = atan2 ((ai), (ar));
251
=head2 cplx real-valued-pdl
253
Cast a real-valued piddle to the complex datatype. The first dimension of
254
the piddle must be of size 2. After this the usual (complex) arithmetic
255
operators are applied to this pdl, rather than the normal elementwise pdl
256
operators. Dataflow to the complex parent works. Use C<sever> on the result
257
if you don't want this.
260
=head2 complex real-valued-pdl
262
Cast a real-valued piddle to the complex datatype I<without> dataflow
263
and I<inplace>. Achieved by merely reblessing a piddle. The first dimension of
264
the piddle must be of size 2.
266
=head2 real cplx-valued-pdl
268
Cast a complex valued pdl back to the "normal" pdl datatype. Afterwards
269
the normal elementwise pdl operators are used in operations. Dataflow
270
to the real parent works. Use C<sever> on the result if you don't want this.
276
return $_[0] if UNIVERSAL::isa($_[0],'PDL::Complex'); # NOOP if just piddle
277
croak "first dimsize must be 2" unless $_[0]->dims > 0 && $_[0]->dim(0) == 2;
278
bless $_[0]->slice('');
282
return $_[0] if UNIVERSAL::isa($_[0],'PDL::Complex'); # NOOP if just piddle
283
croak "first dimsize must be 2" unless $_[0]->dims > 0 && $_[0]->dim(0) == 2;
288
*PDL::complex = \&complex;
291
return $_[0] unless UNIVERSAL::isa($_[0],'PDL::Complex'); # NOOP unless complex
292
bless $_[0]->slice(''), 'PDL';
298
Pars => 'r(); [o]c(m=2)',
299
Doc => 'convert real to complex, assuming an imaginary part of zero',
302
*PDL::r2C = \&PDL::Complex::r2C;
303
sub PDL::Complex::r2C($) {
304
return $_[0] if UNIVERSAL::isa($_[0],'PDL::Complex');
305
my $r = __PACKAGE__->initialize;
306
&PDL::Complex::_r2C_int($_[0], $r);
317
Pars => 'r(); [o]c(m=2)',
318
Doc => 'convert imaginary to complex, assuming a real part of zero',
319
PMCode => '*PDL::i2C = \&PDL::Complex::i2C; sub PDL::Complex::i2C($) { my $r = __PACKAGE__->initialize; &PDL::Complex::_i2C_int($_[0], $r); $r }',
327
Pars => 'r(m=2); float+ [o]p(m=2)',
329
Doc => 'convert complex numbers in rectangular form to polar (mod,arg) form. Works inplace',
331
$GENERIC() x = $r(m=>0);
332
$GENERIC() y = $r(m=>1);
333
$p(m=>0) = CABS (x, y);
334
$p(m=>1) = atan2 (y, x);
339
Pars => 'r(m=2); [o]p(m=2)',
341
GenericTypes => [F,D],
342
Doc => 'convert complex numbers in polar (mod,arg) form to rectangular form. Works inplace',
344
$GENERIC() m = $r(m=>0);
345
$GENERIC() a = $r(m=>1);
354
pp_def 'Cadd', # this is here for a) completeness and b) not having to mess with PDL::Ops
355
Pars => 'a(m=2); b(m=2); [o]c(m=2)',
358
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
359
$GENERIC() br = $b(m=>0), bi = $b(m=>1);
365
pp_def 'Csub', # this is here for a) completeness and b) not having to mess with PDL::Ops
366
Pars => 'a(m=2); b(m=2); [o]c(m=2)',
369
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
370
$GENERIC() br = $b(m=>0), bi = $b(m=>1);
377
Pars => 'a(m=2); b(m=2); [o]c(m=2)',
378
Doc => 'complex multiplication',
380
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
381
$GENERIC() br = $b(m=>0), bi = $b(m=>1);
382
$c(m=>0) = ar*br - ai*bi;
383
$c(m=>1) = ar*bi + ai*br;
388
Pars => 'a(m=2,n); [o]c(m=2)',
389
Doc => 'Project via product to N-1 dimension',
392
$GENERIC() br, bi, cr, ci,tmp;
395
for (iter=1; iter < $SIZE(n);iter++)
397
br = $a(m=>0,n=>iter);
398
bi = $a(m=>1,n=>iter);
409
Pars => 'a(m=2); b(); [o]c(m=2)',
410
Doc => 'mixed complex/real multiplication',
412
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
413
$c(m=>0) = ar * $b();
414
$c(m=>1) = ai * $b();
419
Pars => 'a(m=2); b(m=2); [o]c(m=2)',
420
GenericTypes => [F,D],
421
Doc => 'complex division',
423
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
424
$GENERIC() br = $b(m=>0), bi = $b(m=>1);
426
if (fabs (br) > fabs (bi))
428
$GENERIC() tt = bi / br;
429
$GENERIC() dn = br + tt * bi;
430
$c(m=>0) = (ar + tt * ai) / dn;
431
$c(m=>1) = (ai - tt * ar) / dn;
435
$GENERIC() tt = br / bi;
436
$GENERIC() dn = br * tt + bi;
437
$c(m=>0) = (ar * tt + ai) / dn;
438
$c(m=>1) = (ai * tt - ar) / dn;
444
Pars => 'a(m=2); b(m=2); [o]c()',
445
GenericTypes => [F,D],
446
Doc => 'Complex comparison oeprator (spaceship). It orders by real first, then by imaginary. Hm, but it is mathematical nonsense! Complex numbers cannot be ordered.',
450
a = $a(m=>0), b = $b(m=>0);
452
$c() = (a > b) * 2 - 1;
455
a = $a(m=>1), b = $b(m=>1);
463
Pars => 'a(m=2); [o]c(m=2)',
465
Doc => 'complex conjugation. Works inplace',
468
$c(m=>1) = -$a(m=>1);
473
Pars => 'a(m=2); [o]c()',
474
GenericTypes => [F,D],
475
Doc => 'complex C<abs()> (also known as I<modulus>)',
476
PMCode => q^sub PDL::Complex::Cabs($) {
479
&PDL::Complex::_Cabs_int($pdl, $abs);
483
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
484
$c() = CABS (ar, ai);
489
Pars => 'a(m=2); [o]c()',
490
Doc => 'complex squared C<abs()> (also known I<squared modulus>)',
491
PMCode => q^sub PDL::Complex::Cabs2($) {
493
my $abs2 = PDL->null;
494
&PDL::Complex::_Cabs2_int($pdl, $abs2);
498
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
499
$c() = ar*ar + ai*ai;
504
Pars => 'a(m=2); [o]c()',
505
GenericTypes => [F,D],
506
Doc => 'complex argument function ("angle")',
507
PMCode => q^sub PDL::Complex::Carg($) {
510
&PDL::Complex::_Carg_int($pdl, $arg);
514
$c() = atan2 ($a(m=>1), $a(m=>0));
519
Pars => 'a(m=2); [o]c(m=2)',
521
GenericTypes => [F,D],
522
Doc => ' sin (a) = 1/(2*i) * (exp (a*i) - exp (-a*i)). Works inplace',
524
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
528
$c(m=>0) = s * cosh (ai);
529
$c(m=>1) = c * sinh (ai);
534
Pars => 'a(m=2); [o]c(m=2)',
536
GenericTypes => [F,D],
537
Doc => ' cos (a) = 1/2 * (exp (a*i) + exp (-a*i)). Works inplace',
539
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
543
$c(m=>0) = c * cosh (ai);
544
$c(m=>1) = - s * sinh (ai);
550
=head2 Ctan a [not inplace]
552
tan (a) = -i * (exp (a*i) - exp (-a*i)) / (exp (a*i) + exp (-a*i))
556
sub Ctan($) { Csin($_[0]) / Ccos($_[0]) }
562
Pars => 'a(m=2); [o]c(m=2)',
564
GenericTypes => [F,D],
565
Doc => 'exp (a) = exp (real (a)) * (cos (imag (a)) + i * sin (imag (a))). Works inplace',
567
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
568
$GENERIC() ex = exp (ar);
578
Pars => 'a(m=2); [o]c(m=2)',
580
GenericTypes => [F,D],
581
Doc => 'log (a) = log (cabs (a)) + i * carg (a). Works inplace',
583
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
585
CLOG (ar, ai, $c(m=>0), $c(m=>1));
590
Pars => 'a(m=2); b(m=2); [o]c(m=2)',
592
GenericTypes => [F,D],
593
Doc => 'complex C<pow()> (C<**>-operator)',
595
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
596
$GENERIC() br = $b(m=>0), bi = $b(m=>1);
598
double logr, logi, x, y;
601
if(ar == 0 && ai == 0){
602
if(br == 0 && bi == 0) {
612
CLOG (ar, ai, logr, logi);
613
x = exp (logr*br - logi*bi);
614
y = logr*bi + logi*br;
619
if(ai == 0 && bi == 0) $c(m=>1) = 0;
620
else $c(m=>1) = x * s;
626
Pars => 'a(m=2); [o]c(m=2)',
628
GenericTypes => [F,D],
629
Doc => 'Works inplace',
631
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
633
CSQRT ($GENERIC(), ar, ai, $c(m=>0), $c(m=>1));
638
Pars => 'a(m=2); [o]c(m=2)',
640
GenericTypes => [F,D],
641
Doc => 'Works inplace',
643
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
645
$GENERIC() t1 = sqrt ((ar+1)*(ar+1) + ai*ai);
646
$GENERIC() t2 = sqrt ((ar-1)*(ar-1) + ai*ai);
647
$GENERIC() alpha = (t1+t2)*0.5;
648
$GENERIC() beta = (t1-t2)*0.5;
650
if (alpha < 1) alpha = 1;
651
if (beta > 1) beta = 1;
652
else if (beta < -1) beta = -1;
654
$c(m=>0) = atan2 (beta, sqrt (1-beta*beta));
655
$c(m=>1) = - log (alpha + sqrt (alpha*alpha-1));
656
if (ai > 0 || (ai == 0 && ar < -1))
657
$c(m=>1) = - $c(m=>1);
662
Pars => 'a(m=2); [o]c(m=2)',
664
GenericTypes => [F,D],
665
Doc => 'Works inplace',
667
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
669
$GENERIC() t1 = sqrt ((ar+1)*(ar+1) + ai*ai);
670
$GENERIC() t2 = sqrt ((ar-1)*(ar-1) + ai*ai);
671
$GENERIC() alpha = (t1+t2)*0.5;
672
$GENERIC() beta = (t1-t2)*0.5;
674
if (alpha < 1) alpha = 1;
675
if (beta > 1) beta = 1;
676
else if (beta < -1) beta = -1;
678
$c(m=>0) = atan2 (sqrt (1-beta*beta), beta);
679
$c(m=>1) = log (alpha + sqrt (alpha*alpha-1));
680
if (ai > 0 || (ai == 0 && ar < -1))
681
$c(m=>1) = - $c(m=>1);
687
=head2 Catan cplx [not inplace]
689
Return the complex C<atan()>.
695
Cmul Clog(Cdiv (PDL::Complex::i+$z, PDL::Complex::i-$z)), pdl(0, 0.5);
701
Pars => 'a(m=2); [o]c(m=2)',
703
GenericTypes => [F,D],
704
Doc => ' sinh (a) = (exp (a) - exp (-a)) / 2. Works inplace',
706
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
710
$c(m=>0) = sinh (ar) * c;
711
$c(m=>1) = cosh (ar) * s;
716
Pars => 'a(m=2); [o]c(m=2)',
718
GenericTypes => [F,D],
719
Doc => ' cosh (a) = (exp (a) + exp (-a)) / 2. Works inplace',
721
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
725
$c(m=>0) = cosh (ar) * c;
726
$c(m=>1) = sinh (ar) * s;
731
Pars => 'a(m=2); [o]c(m=2)',
733
GenericTypes => [F,D],
734
Doc => 'Works inplace',
736
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
741
den = cosh (2*ar) + c;
743
$c(m=>0) = sinh (2*ar) / den;
749
Pars => 'a(m=2); [o]c(m=2)',
751
GenericTypes => [F,D],
752
Doc => 'Works inplace',
754
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
755
$GENERIC() yr = (ar-ai) * (ar+ai) + 1;
756
$GENERIC() yi = 2*ar*ai;
758
CSQRT ($GENERIC(), yr, yi, yr, yi)
763
CLOG (yr, yi, $c(m=>0), $c(m=>1));
768
Pars => 'a(m=2); [o]c(m=2)',
770
GenericTypes => [F,D],
771
Doc => 'Works inplace',
773
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
775
$GENERIC() yr = (ar-ai) * (ar+ai) - 1;
776
$GENERIC() yi = 2*ar*ai;
778
CSQRT ($GENERIC(), yr, yi, yr, yi)
783
CLOG (yr, yi, $c(m=>0), $c(m=>1));
788
Pars => 'a(m=2); [o]c(m=2)',
790
GenericTypes => [F,D],
791
Doc => 'Works inplace',
793
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
796
double num = i2 + (1+ar) * (1+ar);
797
double den = i2 + (1-ar) * (1-ar);
799
$c(m=>0) = 0.25 * (log(num) - log(den));
800
$c(m=>1) = 0.5 * atan2 (2*ai, 1 - ar*ar - i2);
805
Pars => 'a(m=2); [o]c(m=2)',
807
GenericTypes => [F,D],
808
Doc => 'compute the projection of a complex number to the riemann sphere. Works inplace',
810
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
812
double den = ar*ar + ai*ai + 1;
814
$c(m=>0) = 2*ar / den;
815
$c(m=>1) = 2*ai / den;
820
Pars => 'a(m=2); [o]c(m=2,n)',
821
OtherPars => 'int n => n',
822
GenericTypes => [F,D],
823
Doc => 'Compute the C<n> roots of C<a>. C<n> must be a positive integer. The result will always be a complex type!',
824
PMCode => q^sub PDL::Complex::Croots($$) {
827
&PDL::Complex::_Croots_int($pdl, $r, $n);
832
double ar = $a(m=>0), ai = $a(m=>1),
833
n1 = 1 / (double)$COMP(n),
834
rr = pow (CABS (ar, ai), n1), /* do not optimize the sqrt out of this expr! */
835
at = atan2 (ai, ar) * n1,
851
=head2 re cplx, im cplx
853
Return the real or imaginary part of the complex number(s) given. These
854
are slicing operators, so data flow works. The real and imaginary parts
855
are returned as piddles (ref eq PDL).
859
sub re($) { bless $_[0]->slice("(0)"), 'PDL'; }
860
sub im($) { bless $_[0]->slice("(1)"), 'PDL'; }
862
*PDL::Complex::re = \&re;
863
*PDL::Complex::im = \&im;
867
pp_def 'rCpolynomial',
868
Pars => 'coeffs(n); x(c=2,m); [o]out(c=2,m)',
869
Doc => 'evaluate the polynomial with (real) coefficients C<coeffs> at the (complex) position(s) C<x>. C<coeffs[0]> is the constant term.',
870
GenericTypes => [F,D],
880
or += $coeffs() * xr;
881
oi += $coeffs() * xi;
884
xr = Xr * $x(c=>0) - xi * $x(c=>1);
885
xi = xi * $x(c=>0) + Xr * $x(c=>1);
896
pp_addpm {At => Bot}, <<'EOD';
898
# overload must be here, so that all the functions can be seen
900
# undocumented compatibility functions
901
sub Catan2($$) { Catan Cdiv $_[1], $_[0] }
902
sub atan2($$) { Catan Cdiv $_[1], $_[0] }
907
if (/(\S+)\+(\w+)/) {
908
$sub = eval 'sub { '.$2.' $_[0], ref $_[1] eq __PACKAGE__ ? $_[1] : r2C $_[1] }';
909
} elsif (/(\S+)\-(\w+)/) {
910
$sub = eval 'sub { my $b = ref $_[1] eq __PACKAGE__ ? $_[1] : r2C $_[1];
911
$_[2] ? '.$2.' $b, $_[0] : '.$2.' $_[0], $b }';
915
if($1 eq "atan2" || $1 eq "<=>") { return ($1, $sub) }
916
($1, $sub, "$1=", $sub);
920
my ($op, $func) = ($_[0] =~ /(.+)@(\w+)/);
921
*$op = \&$func if $op =~ /\w+/; # create an alias
922
($op, eval 'sub { '.$func.' $_[0] }');
926
($_[0], eval 'sub { my $b = ref $_[1] eq __PACKAGE__ ? $_[1] : r2C $_[1];
927
($_[2] ? $b <=> $_[0] : $_[0] <=> $b) '.$_[0].' 0 }');
931
# Bless a null PDL into the supplied 1st arg package
932
# If 1st arg is a ref, get the package from it
933
bless PDL->null, ref($_[0]) ? ref($_[0]) : $_[0];
937
(map _gen_biop($_), qw(++Cadd --Csub *+Cmul /-Cdiv **-Cpow atan2-Catan2 <=>-Ccmp)),
938
(map _gen_unop($_), qw(sin@Csin cos@Ccos exp@Cexp abs@Cabs log@Clog sqrt@Csqrt abs@Cabs)),
939
(map _gen_cpop($_), qw(< <= == != >= >)),
940
'++' => sub { $_[0] += 1 },
941
'--' => sub { $_[0] -= 1 },
942
'""' => \&PDL::Complex::string
945
# overwrite PDL's overloading to honour subclass methods in + - * /
948
# This strange usage of BEGINs is to ensure the
949
# warning messages get disabled and enabled in the
950
# proper order. Without the BEGIN's the 'use overload'
951
# would be called first.
952
BEGIN {$warningFlag = $^W; # Temporarily disable warnings caused by
953
$^W = 0; # redefining PDL's subs
960
&& (ref $_[1] ne 'PDL')
961
&& defined ($foo = overload::Method($_[1],'+')))
962
{ &$foo($_[1], $_[0], !$_[2])}
963
else { PDL::plus (@_)}
969
&& (ref $_[1] ne 'PDL')
970
&& defined ($foo = overload::Method($_[1],'*')))
971
{ &$foo($_[1], $_[0], !$_[2])}
972
else { PDL::mult (@_)}
978
&& (ref $_[1] ne 'PDL')
979
&& defined ($foo = overload::Method($_[1],'-')))
980
{ &$foo($_[1], $_[0], !$_[2])}
981
else { PDL::minus (@_)}
987
&& (ref $_[1] ne 'PDL')
988
&& defined ($foo = overload::Method($_[1],'/')))
989
{ &$foo($_[1], $_[0], !$_[2])}
990
else { PDL::divide (@_)}
994
# Used in overriding standard PDL +, -, *, / ops in the complex subclass.
1004
BEGIN{ $^W = $warningFlag;} # Put Back Warnings
1010
our $floatformat = "%4.4g"; # Default print format for long numbers
1011
our $doubleformat = "%6.6g";
1013
$PDL::Complex::_STRINGIZING = 0;
1015
sub PDL::Complex::string {
1016
my($self,$format1,$format2)=@_;
1017
my @dims = $self->dims;
1018
return PDL::string($self) if ($dims[0] != 2);
1020
if($PDL::Complex::_STRINGIZING) {
1021
return "ALREADY_STRINGIZING_NO_LOOPS";
1023
local $PDL::Complex::_STRINGIZING = 1;
1024
my $ndims = $self->getndims;
1025
if($self->nelem > $PDL::toolongtoprint) {
1026
return "TOO LONG TO PRINT";
1029
PDL::Core::string($self,$format1);
1031
return "Null" if $self->isnull;
1032
return "Empty" if $self->isempty; # Empty piddle
1033
local $sep = $PDL::use_commas ? ", " : " ";
1034
local $sep2 = $PDL::use_commas ? ", " : "";
1036
return str1D($self,$format1,$format2);
1039
return strND($self,$format1,$format2,0);
1046
my $tmp = $x->mv(0,1)->clump(0,2)->mv(1,0)->sumover;
1047
return $tmp->squeeze;
1052
PDL::Ufunc::sumover($m->xchg(0,1));
1057
my($self,$format1,$format2,$level)=@_;
1058
my @dims = $self->dims;
1061
return str2D($self,$format1,$format2,$level);
1064
my $secbas = join '',map {":,"} @dims[0..$#dims-1];
1065
my $ret="\n"." "x$level ."["; my $j;
1066
for ($j=0; $j<$dims[$#dims]; $j++) {
1067
my $sec = $secbas . "($j)";
1069
$ret .= strND($self->slice($sec),$format1,$format2, $level+1);
1070
chop $ret; $ret .= $sep2;
1072
chop $ret if $PDL::use_commas;
1073
$ret .= "\n" ." "x$level ."]\n";
1079
# String 1D array in nice format
1082
my($self,$format1,$format2)=@_;
1083
barf "Not 1D" if $self->getndims() > 2;
1084
my $x = PDL::Core::listref_c($self);
1085
my ($ret,$dformat,$t, $i);
1087
my $dtype = $self->get_datatype();
1088
$dformat = $PDL::Complex::floatformat if $dtype == $PDL_F;
1089
$dformat = $PDL::Complex::doubleformat if $dtype == $PDL_D;
1091
$ret = "[" if $self->getndims() > 1;
1092
my $badflag = $self->badflag();
1093
for($i=0; $i<=$#$x; $i++){
1095
if ( $badflag and $t eq "BAD" ) {
1097
} elsif ($format1) {
1098
$t = sprintf $format1,$t;
1100
if ($dformat && length($t)>7) { # Try smaller
1101
$t = sprintf $dformat,$t;
1105
$i<$#$x ? $t."i$sep" : $t."i"
1106
: substr($$x[$i+1],0,1) eq "-" ? "$t " : $t." +";
1108
$ret.="]" if $self->getndims() > 1;
1114
my($self,$format1,$format2,$level)=@_;
1115
my @dims = $self->dims();
1116
barf "Not 2D" if scalar(@dims)!=3;
1117
my $x = PDL::Core::listref_c($self);
1118
my ($i, $f, $t, $len1, $len2, $ret);
1120
my $dtype = $self->get_datatype();
1121
my $badflag = $self->badflag();
1125
if (!defined $format1 || !defined $format2 ||
1126
$format1 eq '' || $format2 eq '') {
1130
for ($i=0; $i<=$#$x; $i++) {
1131
if ( $$x[$i] eq "BAD" ) {
1135
$f = length($$x[$i]);
1138
$len2 = $f if $f > $len2;
1141
$len1 = $f if $f > $len1;
1145
for ($i=0; $i<=$#$x; $i++) {
1146
$f = length($$x[$i]);
1148
$len2 = $f if $f > $len2;
1151
$len1 = $f if $f > $len1;
1156
$format1 = '%'.$len1.'s';
1157
$format2 = '%'.$len2.'s';
1160
if ($dtype == $PDL_F) {
1161
$format1 = $PDL::Complex::floatformat;
1163
} elsif ($dtype == $PDL_D) {
1164
$format1 = $PDL::Complex::doubleformat;
1171
if ($dtype == $PDL_F) {
1172
$format2 = $PDL::Complex::floatformat;
1174
} elsif ($dtype == $PDL_D) {
1175
$format2 = $PDL::Complex::doubleformat;
1178
$findmax = 0 unless $findmax;
1187
for($i=0; $i<=$#$x; $i++){
1189
if ( $$x[$i] eq 'BAD' ){
1193
$f = $findmax ? length(sprintf $format2,$$x[$i]) :
1194
length(sprintf $format1,$$x[$i]);
1197
$len2 = $f if $f > $len2;
1200
$len1 = $f if $f > $len1;
1204
for ($i=0; $i<=$#$x; $i++) {
1206
$f = length(sprintf $format2,$$x[$i]);
1207
$len2 = $f if $f > $len2;
1210
$f = length(sprintf $format1,$$x[$i]);
1211
$len1 = $f if $f > $len1;
1219
$ret = "\n" . ' 'x$level . "[\n";
1221
my $level = $level+1;
1222
$ret .= ' 'x$level .'[';
1225
for ($i=0; $i<=$#$x; $i++) {
1228
if ( $badflag and $$x[$i] eq 'BAD' ){
1230
#($findmax && $$x[$i - 1 ] eq 'BAD') ||
1231
#(!$findmax && $$x[$i +1 ] eq 'BAD')){
1235
$f = sprintf $format2, $$x[$i];
1236
if (substr($$x[$i],0,1) eq '-'){
1240
$f =~ s/(\s*)(.*)/+$2i/;
1243
$t = $len2-length($f);
1246
if ( $badflag and $$x[$i] eq 'BAD' ){
1250
$f = sprintf $format1, $$x[$i];
1251
$t = $len1-length($f);
1255
$f = ' 'x$t.$f if $t>0;
1258
if (($i+1)%($dims[1]*2)) {
1259
$ret.=$sep if $findmax;
1261
else{ # End of output line
1263
if ($i==$#$x) { # very last number
1267
$ret.= $sep2."\n" . ' 'x$level .'[';
1272
$ret .= ' 'x$level."]\n";
1280
Copyright (C) 2000 Marc Lehmann <pcg@goof.com>.
1281
All rights reserved. There is no warranty. You are allowed
1282
to redistribute this software / documentation as described
1283
in the file COPYING in the PDL distribution.
added
removed
.pc/fix_pod2man_errors.patch/Basic/Complex/complex.pd
