3
# pp_setversion $VERSION; # haven't worked out why it breaks my system (CS)
4
pp_beginwrap; # required for overload to work
6
pp_addpm {At => Top}, <<'EOD';
9
PDL::Complex - handle complex numbers
18
This module features a growing number of functions manipulating complex
19
numbers. These are usually represented as a pair C<[ real imag ]> or
20
C<[ angle phase ]>. If not explicitly mentioned, the functions can work
21
inplace (not yet implemented!!!) and require rectangular form.
23
While there is a procedural interface available (C<$a/$b*$c <=> Cmul
24
(Cdiv $a, $b), $c)>), you can also opt to cast your pdl's into the
25
C<PDL::Complex> datatype, which works just like your normal piddles, but
26
with all the normal perl operators overloaded.
28
The latter means that C<sin($a) + $b/$c> will be evaluated using the
29
normal rules of complex numbers, while other pdl functions (like C<max>)
30
just treat the piddle as a real-valued piddle with a lowest dimension of
31
size 2, so C<max> will return the maximum of all real and imaginary parts,
32
not the "highest" (for some definition)
34
=head1 TIPS, TRICKS & CAVEATS
40
C<i> is a constant exported by this module, which represents
41
C<-1**0.5>, i.e. the imaginary unit. it can be used to quickly and
42
conviniently write complex constants like this: C<4+3*i>.
46
Use C<r2C(real-values)> to convert from real to complex, as in C<$r
47
= Cpow $cplx, r2C 2>. The overloaded operators automatically do that for
48
you, all the other functions, do not. So C<Croots 1, 5> will return all
49
the fifths roots of 1+1*i (due to threading).
53
use C<cplx(real-valued-piddle)> to cast from normal piddles intot he
54
complex datatype. Use C<real(complex-valued-piddle)> to cast back. This
55
requires a copy, though.
59
BEWARE: This module has not been extensively tested. Watch out and
64
=head1 EXAMPLE WALK-THROUGH
66
The complex constant five is equal to C<pdl(1,0)>:
71
Now calculate the three roots of of five:
73
perldl> p $r = Croots $x, 3
77
[-0.85498797 1.4808826]
78
[-0.85498797 -1.4808826]
81
Check that these really are the roots of unity:
91
Duh! Could be better. Now try by multiplying C<$r> three times with itself:
101
Well... maybe C<Cpow> (which is used by the C<**> operator) isn't as
102
bad as I thought. Now multiply by C<i> and negate, which is just a very
103
expensive way of swapping real and imaginary parts.
109
[ 1.4808826 -0.85498797]
110
[ -1.4808826 -0.85498797]
113
Now plot the magnitude of (part of) the complex sine. First generate the
116
perldl> $sin = i * zeroes(50)->xlinvals(2,4)
117
+ zeroes(50)->xlinvals(0,7)
119
Now plot the imaginary part, the real part and the magnitude of the sine
120
into the same diagram:
122
perldl> line im sin $sin; hold
123
perldl> line re sin $sin
124
perldl> line abs sin $sin
126
Sorry, but I didn't yet try to reproduce the diagram in this
127
text. Just run the commands yourself, making sure that you have loaded
128
C<PDL::Complex> (and C<PDL::Graphics::PGPLOT>).
135
for (qw(Ctan Catan re im i cplx real)) {
136
pp_add_exported '', $_;
144
# define M_PI 3.1415926535897932384626433832795029
147
# define M_2PI (2. * M_PI)
150
#if __GLIBC__ > 1 && (defined __USE_MISC || defined __USE_XOPEN || defined __USE_ISOC9X)
151
# define CABS(r,i) hypot (r, i)
154
CABS (double r, double i)
170
return r * sqrt (1 + t*t);
174
#if __GLIBC__ >= 2 && __GLIBC_MINOR__ >= 1 && defined __USE_GNU
175
# define SINCOS(x,s,c) sincos ((x), &(s), &(c))
177
# define SINCOS(x,s,c) \
183
#define CSQRT(type,ar,ai,cr,ci) \
184
type mag = CABS ((ar), (ai)); \
191
t = sqrt (0.5 * (mag + (ar))); \
193
(ci) = 0.5 * (ai) / t; \
197
t = sqrt (0.5 * (mag - (ar))); \
203
(cr) = 0.5 * (ai) / t; \
206
#define CLOG(ar,ai,cr,ci) \
207
(cr) = log (CABS ((ar), (ai))); \
208
(ci) = atan2 ((ai), (ar));
214
=head2 cplx real-valued-pdl
216
Cast a real-valued piddle to the complex datatype. The first dimension of
217
the piddle must be of size 2. After this the usual (complex) arithmetic
218
operators are applied to this pdl, rather than the normal elementwise pdl
219
operators. Dataflow to the complex parent works. Use C<sever> on the result
220
if you don't want this.
223
=head2 real cplx-valued-pdl
225
Cast a complex vlaued pdl back to the "normal" pdl datatype. Afterwards
226
the normal elementwise pdl operators are used in operations. Dataflow
227
to the real parent works. Use C<sever> on the result if you don't want this.
232
bless $_[0]->slice('');
238
bless $_[0]->slice(''), 'PDL';
244
Pars => 'r(); [o]c(m=2)',
245
Doc => 'convert real to complex, assuming an imaginary part of zero',
246
PMCode => 'sub PDL::r2C($) { my $r = __PACKAGE__->initialize; &PDL::_r2C_int($_[0], $r); $r }',
254
Pars => 'r(); [o]c(m=2)',
255
Doc => 'convert imaginary to complex, assuming a real part of zero',
256
PMCode => 'sub PDL::i2C($) { my $r = __PACKAGE__->initialize; &PDL::_i2C_int($_[0], $r); $r }',
264
Pars => 'r(m=2); float+ [o]p(m=2)',
265
Doc => 'convert complex numbers in rectangular form to polar (mod,arg) form',
267
$GENERIC() x = $r(m=>0);
268
$GENERIC() y = $r(m=>1);
269
$p(m=>0) = CABS (x, y);
270
$p(m=>1) = atan2 (y, x);
275
Pars => 'r(m=2); [o]p(m=2)',
276
GenericTypes => [F,D],
277
Doc => 'convert complex numbers in polar (mod,arg) form to rectangular form',
279
$GENERIC() m = $r(m=>0);
280
$GENERIC() a = $r(m=>1);
289
pp_def 'Cadd', # this is here for a) completeness and b) not having to mess with PDL::Ops
290
Pars => 'a(m=2); b(m=2); [o]c(m=2)',
293
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
294
$GENERIC() br = $b(m=>0), bi = $b(m=>1);
300
pp_def 'Csub', # this is here for a) completeness and b) not having to mess with PDL::Ops
301
Pars => 'a(m=2); b(m=2); [o]c(m=2)',
304
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
305
$GENERIC() br = $b(m=>0), bi = $b(m=>1);
312
Pars => 'a(m=2); b(m=2); [o]c(m=2)',
313
Doc => 'complex multiplication',
315
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
316
$GENERIC() br = $b(m=>0), bi = $b(m=>1);
317
$c(m=>0) = ar*br - ai*bi;
318
$c(m=>1) = ar*bi + ai*br;
323
Pars => 'a(m=2); b(); [o]c(m=2)',
324
Doc => 'mixed complex/real multiplication',
326
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
327
$c(m=>0) = ar * $b();
328
$c(m=>1) = ai * $b();
333
Pars => 'a(m=2); b(m=2); [o]c(m=2)',
334
GenericTypes => [F,D],
335
Doc => 'complex division',
337
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
338
$GENERIC() br = $b(m=>0), bi = $b(m=>1);
340
if (fabs (br) > fabs (bi))
342
$GENERIC() tt = bi / br;
343
$GENERIC() dn = br + tt * bi;
344
$c(m=>0) = (ar + tt * ai) / dn;
345
$c(m=>1) = (ai - tt * ar) / dn;
349
$GENERIC() tt = br / bi;
350
$GENERIC() dn = br * tt + bi;
351
$c(m=>0) = (ar * tt + ai) / dn;
352
$c(m=>1) = (ai * tt - ar) / dn;
358
Pars => 'a(m=2); b(m=2); [o]c()',
359
GenericTypes => [F,D],
360
Doc => 'Complex comparison oeprator (spaceship). It orders by real first, then by imaginary.',
364
a = $a(m=>0), b = $b(m=>0);
366
$c() = (a > b) * 2 - 1;
369
a = $a(m=>1), b = $b(m=>1);
377
Pars => 'a(m=2); [o]c(m=2)',
378
Doc => 'complex conjugation',
381
$c(m=>1) = -$a(m=>1);
386
Pars => 'a(m=2); [o]c()',
387
GenericTypes => [F,D],
388
Doc => 'complex C<abs()> (also known as I<modulus>)',
390
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
391
$c() = CABS (ar, ai);
396
Pars => 'a(m=2); [o]c()',
397
Doc => 'complex squared C<abs()> (also known I<squared modulus>)',
399
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
400
$c() = ar*ar + ai*ai;
405
Pars => 'a(m=2); [o]c()',
406
GenericTypes => [F,D],
407
Doc => 'complex argument function ("angle")',
409
$c() = atan2 ($a(m=>1), $a(m=>0));
414
Pars => 'a(m=2); [o]c(m=2)',
415
GenericTypes => [F,D],
416
Doc => ' sin (a) = 1/(2*i) * (exp (a*i) - exp (-a*i))',
418
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
422
$c(m=>0) = s * cosh (ai);
423
$c(m=>1) = c * sinh (ai);
428
Pars => 'a(m=2); [o]c(m=2)',
429
GenericTypes => [F,D],
430
Doc => ' cos (a) = 1/2 * (exp (a*i) + exp (-a*i))',
432
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
436
$c(m=>0) = c * cosh (ai);
437
$c(m=>1) = - s * sinh (ai);
443
=head2 Ctan a [not inplace]
445
tan (a) = -i * (exp (a*i) - exp (-a*i)) / (exp (a*i) + exp (-a*i))
449
sub Ctan($) { Csin($_[0]) / Ccos($_[0]) }
454
Pars => 'a(m=2); [o]c(m=2)',
455
GenericTypes => [F,D],
456
Doc => 'exp (a) = exp (real (a)) * (cos (imag (a)) + i * sin (imag (a)))',
458
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
459
$GENERIC() ex = exp (ar);
469
Pars => 'a(m=2); [o]c(m=2)',
470
GenericTypes => [F,D],
471
Doc => 'log (a) = log (cabs (a)) + i * carg (a)',
473
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
475
CLOG (ar, ai, $c(m=>0), $c(m=>1));
480
Pars => 'a(m=2); b(m=2); [o]c(m=2)',
481
GenericTypes => [F,D],
482
Doc => 'complex C<pow()> (C<**>-operator)',
484
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
485
$GENERIC() br = $b(m=>0), bi = $b(m=>1);
487
double logr, logi, x, y;
490
CLOG (ar, ai, logr, logi);
491
x = exp (logr*br - logi*bi);
492
y = logr*bi + logi*br;
501
Pars => 'a(m=2); [o]c(m=2)',
502
GenericTypes => [F,D],
505
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
507
CSQRT ($GENERIC(), ar, ai, $c(m=>0), $c(m=>1));
512
Pars => 'a(m=2); [o]c(m=2)',
513
GenericTypes => [F,D],
516
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
518
$GENERIC() t1 = sqrt ((ar+1)*(ar+1) + ai*ai);
519
$GENERIC() t2 = sqrt ((ar-1)*(ar-1) + ai*ai);
520
$GENERIC() alpha = (t1+t2)*0.5;
521
$GENERIC() beta = (t1-t2)*0.5;
523
if (alpha < 1) alpha = 1;
524
if (beta > 1) beta = 1;
525
else if (beta < -1) beta = -1;
527
$c(m=>0) = atan2 (beta, sqrt (1-beta*beta));
528
$c(m=>1) = - log (alpha + sqrt (alpha*alpha-1));
529
if (ai > 0 || (ai == 0 && ar < -1))
530
$c(m=>1) = - $c(m=>1);
535
Pars => 'a(m=2); [o]c(m=2)',
536
GenericTypes => [F,D],
539
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
541
$GENERIC() t1 = sqrt ((ar+1)*(ar+1) + ai*ai);
542
$GENERIC() t2 = sqrt ((ar-1)*(ar-1) + ai*ai);
543
$GENERIC() alpha = (t1+t2)*0.5;
544
$GENERIC() beta = (t1-t2)*0.5;
546
if (alpha < 1) alpha = 1;
547
if (beta > 1) beta = 1;
548
else if (beta < -1) beta = -1;
550
$c(m=>0) = atan2 (sqrt (1-beta*beta), beta);
551
$c(m=>1) = log (alpha + sqrt (alpha*alpha-1));
552
if (ai > 0 || (ai == 0 && ar < -1))
553
$c(m=>1) = - $c(m=>1);
559
=head2 Catan cplx [not inplace]
561
Return the complex C<atan()>.
567
Cmul pdl(0, 0.5), Clog Cdiv PDL::Complex::i+$z, PDL::Complex::i-$z;
573
Pars => 'a(m=2); [o]c(m=2)',
574
GenericTypes => [F,D],
575
Doc => ' sinh (a) = (exp (a) - exp (-a)) / 2',
577
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
581
$c(m=>0) = sinh (ar) * c;
582
$c(m=>0) = cosh (ar) * s;
587
Pars => 'a(m=2); [o]c(m=2)',
588
GenericTypes => [F,D],
589
Doc => ' cosh (a) = (exp (a) + exp (-a)) / 2',
591
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
595
$c(m=>0) = cosh (ar) * c;
596
$c(m=>0) = sinh (ar) * s;
601
Pars => 'a(m=2); [o]c(m=2)',
602
GenericTypes => [F,D],
605
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
610
den = cosh (2*ar) + c;
612
$c(m=>0) = sinh (2*ar) / den;
618
Pars => 'a(m=2); [o]c(m=2)',
619
GenericTypes => [F,D],
622
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
624
$GENERIC() yr = (ar-ai) * (ar+ai) + 1;
625
$GENERIC() yi = 2*ar*ai;
627
CSQRT ($GENERIC(), yr, yi, yr, yi)
632
CLOG (yr, yi, $c(m=>0), $c(m=>1));
637
Pars => 'a(m=2); [o]c(m=2)',
638
GenericTypes => [F,D],
641
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
643
$GENERIC() yr = (ar-ai) * (ar+ai) - 1;
644
$GENERIC() yi = 2*ar*ai;
646
CSQRT ($GENERIC(), yr, yi, yr, yi)
651
CLOG (yr, yi, $c(m=>0), $c(m=>1));
656
Pars => 'a(m=2); [o]c(m=2)',
657
GenericTypes => [F,D],
660
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
663
double num = i2 + (1+ar) * (1+ar);
664
double den = i2 + (1-ar) * (1-ar);
666
$c(m=>0) = 0.25 * (log(num) - log(den));
667
$c(m=>1) = 0.5 * atan2 (2*ai, 1 - ar*ar - i2);
672
Pars => 'a(m=2); [o]c(m=2)',
673
GenericTypes => [F,D],
674
Doc => 'compute the projection of a complex number to the riemann sphere',
676
$GENERIC() ar = $a(m=>0), ai = $a(m=>1);
678
double den = ar*ar + ai*ai + 1;
680
$c(m=>0) = 2*ar / den;
681
$c(m=>1) = 2*ai / den;
686
Pars => 'a(m=2); [o]c(m=2,n)',
687
OtherPars => 'int n => n',
688
GenericTypes => [F,D],
689
Doc => 'Compute the C<n> roots of C<a>. C<n> must be a positive integer. The result will always be a complex type!',
690
PMCode => q^sub PDL::Croots($$) {
693
&PDL::_Croots_int($pdl, $r, $n);
697
double ar = $a(m=>0), ai = $a(m=>1),
698
n1 = 1 / (double)$COMP(n),
699
rr = pow (CABS (ar, ai), n1), /* do not optimize the sqrt out of this expr! */
718
=head2 re cplx, im cplx
720
Return the real or imaginary part of the complex number(s) given. These
721
are slicing operators, so data flow works. The real and imaginary parts
722
are returned as piddles (ref eq PDL).
726
sub re($) { bless $_[0]->slice("(0)"), 'PDL' }
727
sub im($) { bless $_[0]->slice("(1)"), 'PDL' }
731
pp_def 'rCpolynomial',
732
Pars => 'coeffs(n); x(c=2,m); [o]out(c=2,m)',
733
Doc => 'evaluate the polynomial with (real) coefficients C<coeffs> at the (complex) position(s) C<x>. C<coeffs[0]> is the constant term.',
734
GenericTypes => [F,D],
744
or += $coeffs() * xr;
745
oi += $coeffs() * xi;
748
xr = Xr * $x(c=>0) - xi * $x(c=>1);
749
xi = xi * $x(c=>0) + Xr * $x(c=>1);
760
pp_addpm {At => Bot}, <<'EOD';
762
# overload must be here, so that all the functions can be seen
764
# undocumented compatibility functions
765
sub Catan2($$) { Catan Cdiv $_[1], $_[0] }
766
sub atan2($$) { Catan Cdiv $_[1], $_[0] }
771
if (/(\S+)\+(\w+)/) {
772
$sub = eval 'sub { '.$2.' $_[0], ref $_[1] eq __PACKAGE__ ? $_[1] : r2C $_[1] }';
773
} elsif (/(\S+)\-(\w+)/) {
774
$sub = eval 'sub { my $b = ref $_[1] eq __PACKAGE__ ? $_[1] : r2C $_[1];
775
$_[2] ? '.$2.' $b, $_[0] : '.$2.' $_[0], $b }';
779
return ($1, $sub) if $1 eq "atan2";
780
($1, $sub, "$1=", $sub);
784
my ($op, $func) = ($_[0] =~ /(.+)@(\w+)/);
785
*$op = \&$func if $op =~ /\w+/; # create an alias
786
($op, eval 'sub { '.$func.' $_[0] }');
790
($_[0], eval 'sub { my $b = ref $_[1] eq __PACKAGE__ ? $_[1] : r2C $_[1];
791
($_[2] ? $b <=> $_[0] : $_[0] <=> $b) '.$_[0].' 0 }');
795
bless PDL->null, $_[0];
798
BEGIN { $i = bless pdl 0,1 }
799
sub i () { $i->copy };
802
(map _gen_biop($_), qw(++Cadd --Csub *+Cmul /-Cdiv **-Cpow atan2-Catan2 <=>-Ccmp)),
803
(map _gen_unop($_), qw(sin@Csin cos@Ccos exp@Cexp abs@Cabs log@Clog sqrt@Csqrt abs@Cabs)),
804
(map _gen_cpop($_), qw(< <= == != => >)),
805
'++' => sub { $_[0] += 1 },
806
'--' => sub { $_[0] -= 1 },
809
# overwrite PDL's overloading to honour subclass methods in + - * /
812
# This strange usage of BEGINs is to ensure the
813
# warning messages get disabled and enabled in the
814
# proper order. Without the BEGIN's the 'use overload'
815
# would be called first.
816
BEGIN {$warningFlag = $^W; # Temporarily disable warnings caused by
817
$^W = 0; # redefining PDL's subs
824
&& (ref $_[1] ne 'PDL')
825
&& defined ($foo = overload::Method($_[1],'+')))
826
{ &$foo($_[1], $_[0], !$_[2])}
827
else { PDL::plus (@_)}
833
&& (ref $_[1] ne 'PDL')
834
&& defined ($foo = overload::Method($_[1],'*')))
835
{ &$foo($_[1], $_[0], !$_[2])}
836
else { PDL::mult (@_)}
842
&& (ref $_[1] ne 'PDL')
843
&& defined ($foo = overload::Method($_[1],'-')))
844
{ &$foo($_[1], $_[0], !$_[2])}
845
else { PDL::minus (@_)}
851
&& (ref $_[1] ne 'PDL')
852
&& defined ($foo = overload::Method($_[1],'/')))
853
{ &$foo($_[1], $_[0], !$_[2])}
854
else { PDL::divide (@_)}
858
# Used in overriding standard PDL +, -, *, / ops in the complex subclass.
868
BEGIN{ $^W = $warningFlag;} # Put Back Warnings
873
Copyright (C) 2000 Marc Lehmann <pcg@goof.com>.
874
All rights reserved. There is no warranty. You are allowed
875
to redistribute this software / documentation as described
876
in the file COPYING in the PDL distribution.
added
removed
Basic/Complex/complex.pd
