home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / perl_mlb.zip / Math / Complex.pm < prev    next >
Text File  |  1997-11-25  |  35KB  |  1,593 lines
  1. #
  2. # Complex numbers and associated mathematical functions
  3. # -- Raphael Manfredi    September 1996
  4. # -- Jarkko Hietaniemi    March-October 1997
  5. # -- Daniel S. Lewart    September-October 1997
  6. #
  7.  
  8. require Exporter;
  9. package Math::Complex;
  10.  
  11. $VERSION = 1.05;
  12.  
  13. # $Id: Complex.pm,v 1.2 1997/10/15 10:08:39 jhi Exp $
  14.  
  15. use strict;
  16.  
  17. use vars qw($VERSION @ISA
  18.         @EXPORT %EXPORT_TAGS
  19.         $package $display
  20.         $i $ip2 $logn %logn);
  21.  
  22. @ISA = qw(Exporter);
  23.  
  24. my @trig = qw(
  25.           pi
  26.           tan
  27.           csc cosec sec cot cotan
  28.           asin acos atan
  29.           acsc acosec asec acot acotan
  30.           sinh cosh tanh
  31.           csch cosech sech coth cotanh
  32.           asinh acosh atanh
  33.           acsch acosech asech acoth acotanh
  34.          );
  35.  
  36. @EXPORT = (qw(
  37.          i Re Im arg
  38.          sqrt log ln
  39.          log10 logn cbrt root
  40.          cplx cplxe
  41.          ),
  42.        @trig);
  43.  
  44. %EXPORT_TAGS = (
  45.     'trig' => [@trig],
  46. );
  47.  
  48. use overload
  49.     '+'    => \&plus,
  50.     '-'    => \&minus,
  51.     '*'    => \&multiply,
  52.     '/'    => \÷,
  53.     '**'    => \&power,
  54.     '<=>'    => \&spaceship,
  55.     'neg'    => \&negate,
  56.     '~'    => \&conjugate,
  57.     'abs'    => \&abs,
  58.     'sqrt'    => \&sqrt,
  59.     'exp'    => \&exp,
  60.     'log'    => \&log,
  61.     'sin'    => \&sin,
  62.     'cos'    => \&cos,
  63.     'tan'    => \&tan,
  64.     'atan2'    => \&atan2,
  65.     qw("" stringify);
  66.  
  67. #
  68. # Package globals
  69. #
  70.  
  71. $package = 'Math::Complex';        # Package name
  72. $display = 'cartesian';            # Default display format
  73.  
  74. #
  75. # Object attributes (internal):
  76. #    cartesian    [real, imaginary] -- cartesian form
  77. #    polar        [rho, theta] -- polar form
  78. #    c_dirty        cartesian form not up-to-date
  79. #    p_dirty        polar form not up-to-date
  80. #    display        display format (package's global when not set)
  81. #
  82.  
  83. #
  84. # ->make
  85. #
  86. # Create a new complex number (cartesian form)
  87. #
  88. sub make {
  89.     my $self = bless {}, shift;
  90.     my ($re, $im) = @_;
  91.     $self->{'cartesian'} = [$re, $im];
  92.     $self->{c_dirty} = 0;
  93.     $self->{p_dirty} = 1;
  94.     return $self;
  95. }
  96.  
  97. #
  98. # ->emake
  99. #
  100. # Create a new complex number (exponential form)
  101. #
  102. sub emake {
  103.     my $self = bless {}, shift;
  104.     my ($rho, $theta) = @_;
  105.     if ($rho < 0) {
  106.         $rho   = -$rho;
  107.         $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
  108.     }
  109.     $self->{'polar'} = [$rho, $theta];
  110.     $self->{p_dirty} = 0;
  111.     $self->{c_dirty} = 1;
  112.     return $self;
  113. }
  114.  
  115. sub new { &make }        # For backward compatibility only.
  116.  
  117. #
  118. # cplx
  119. #
  120. # Creates a complex number from a (re, im) tuple.
  121. # This avoids the burden of writing Math::Complex->make(re, im).
  122. #
  123. sub cplx {
  124.     my ($re, $im) = @_;
  125.     return $package->make($re, defined $im ? $im : 0);
  126. }
  127.  
  128. #
  129. # cplxe
  130. #
  131. # Creates a complex number from a (rho, theta) tuple.
  132. # This avoids the burden of writing Math::Complex->emake(rho, theta).
  133. #
  134. sub cplxe {
  135.     my ($rho, $theta) = @_;
  136.     return $package->emake($rho, defined $theta ? $theta : 0);
  137. }
  138.  
  139. #
  140. # pi
  141. #
  142. # The number defined as pi = 180 degrees
  143. #
  144. use constant pi => 4 * atan2(1, 1);
  145.  
  146. #
  147. # pit2
  148. #
  149. # The full circle
  150. #
  151. use constant pit2 => 2 * pi;
  152.  
  153. #
  154. # pip2
  155. #
  156. # The quarter circle
  157. #
  158. use constant pip2 => pi / 2;
  159.  
  160. #
  161. # uplog10
  162. #
  163. # Used in log10().
  164. #
  165. use constant uplog10 => 1 / log(10);
  166.  
  167. #
  168. # i
  169. #
  170. # The number defined as i*i = -1;
  171. #
  172. sub i () {
  173.         return $i if ($i);
  174.     $i = bless {};
  175.     $i->{'cartesian'} = [0, 1];
  176.     $i->{'polar'}     = [1, pip2];
  177.     $i->{c_dirty} = 0;
  178.     $i->{p_dirty} = 0;
  179.     return $i;
  180. }
  181.  
  182. #
  183. # Attribute access/set routines
  184. #
  185.  
  186. sub cartesian {$_[0]->{c_dirty} ?
  187.            $_[0]->update_cartesian : $_[0]->{'cartesian'}}
  188. sub polar     {$_[0]->{p_dirty} ?
  189.            $_[0]->update_polar : $_[0]->{'polar'}}
  190.  
  191. sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
  192. sub set_polar     { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
  193.  
  194. #
  195. # ->update_cartesian
  196. #
  197. # Recompute and return the cartesian form, given accurate polar form.
  198. #
  199. sub update_cartesian {
  200.     my $self = shift;
  201.     my ($r, $t) = @{$self->{'polar'}};
  202.     $self->{c_dirty} = 0;
  203.     return $self->{'cartesian'} = [$r * cos $t, $r * sin $t];
  204. }
  205.  
  206. #
  207. #
  208. # ->update_polar
  209. #
  210. # Recompute and return the polar form, given accurate cartesian form.
  211. #
  212. sub update_polar {
  213.     my $self = shift;
  214.     my ($x, $y) = @{$self->{'cartesian'}};
  215.     $self->{p_dirty} = 0;
  216.     return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
  217.     return $self->{'polar'} = [sqrt($x*$x + $y*$y), atan2($y, $x)];
  218. }
  219.  
  220. #
  221. # (plus)
  222. #
  223. # Computes z1+z2.
  224. #
  225. sub plus {
  226.     my ($z1, $z2, $regular) = @_;
  227.     my ($re1, $im1) = @{$z1->cartesian};
  228.     $z2 = cplx($z2) unless ref $z2;
  229.     my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
  230.     unless (defined $regular) {
  231.         $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
  232.         return $z1;
  233.     }
  234.     return (ref $z1)->make($re1 + $re2, $im1 + $im2);
  235. }
  236.  
  237. #
  238. # (minus)
  239. #
  240. # Computes z1-z2.
  241. #
  242. sub minus {
  243.     my ($z1, $z2, $inverted) = @_;
  244.     my ($re1, $im1) = @{$z1->cartesian};
  245.     $z2 = cplx($z2) unless ref $z2;
  246.     my ($re2, $im2) = @{$z2->cartesian};
  247.     unless (defined $inverted) {
  248.         $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
  249.         return $z1;
  250.     }
  251.     return $inverted ?
  252.         (ref $z1)->make($re2 - $re1, $im2 - $im1) :
  253.         (ref $z1)->make($re1 - $re2, $im1 - $im2);
  254.  
  255. }
  256.  
  257. #
  258. # (multiply)
  259. #
  260. # Computes z1*z2.
  261. #
  262. sub multiply {
  263.         my ($z1, $z2, $regular) = @_;
  264.     if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
  265.         # if both polar better use polar to avoid rounding errors
  266.         my ($r1, $t1) = @{$z1->polar};
  267.         my ($r2, $t2) = @{$z2->polar};
  268.         my $t = $t1 + $t2;
  269.         if    ($t >   pi()) { $t -= pit2 }
  270.         elsif ($t <= -pi()) { $t += pit2 }
  271.         unless (defined $regular) {
  272.         $z1->set_polar([$r1 * $r2, $t]);
  273.         return $z1;
  274.         }
  275.         return (ref $z1)->emake($r1 * $r2, $t);
  276.     } else {
  277.         my ($x1, $y1) = @{$z1->cartesian};
  278.         if (ref $z2) {
  279.         my ($x2, $y2) = @{$z2->cartesian};
  280.         return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
  281.         } else {
  282.         return (ref $z1)->make($x1*$z2, $y1*$z2);
  283.         }
  284.     }
  285. }
  286.  
  287. #
  288. # _divbyzero
  289. #
  290. # Die on division by zero.
  291. #
  292. sub _divbyzero {
  293.     my $mess = "$_[0]: Division by zero.\n";
  294.  
  295.     if (defined $_[1]) {
  296.     $mess .= "(Because in the definition of $_[0], the divisor ";
  297.     $mess .= "$_[1] " unless ($_[1] eq '0');
  298.     $mess .= "is 0)\n";
  299.     }
  300.  
  301.     my @up = caller(1);
  302.  
  303.     $mess .= "Died at $up[1] line $up[2].\n";
  304.  
  305.     die $mess;
  306. }
  307.  
  308. #
  309. # (divide)
  310. #
  311. # Computes z1/z2.
  312. #
  313. sub divide {
  314.     my ($z1, $z2, $inverted) = @_;
  315.     if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
  316.         # if both polar better use polar to avoid rounding errors
  317.         my ($r1, $t1) = @{$z1->polar};
  318.         my ($r2, $t2) = @{$z2->polar};
  319.         my $t;
  320.         if ($inverted) {
  321.         _divbyzero "$z2/0" if ($r1 == 0);
  322.         $t = $t2 - $t1;
  323.         if    ($t >   pi()) { $t -= pit2 }
  324.         elsif ($t <= -pi()) { $t += pit2 }
  325.         return (ref $z1)->emake($r2 / $r1, $t);
  326.         } else {
  327.         _divbyzero "$z1/0" if ($r2 == 0);
  328.         $t = $t1 - $t2;
  329.         if    ($t >   pi()) { $t -= pit2 }
  330.         elsif ($t <= -pi()) { $t += pit2 }
  331.         return (ref $z1)->emake($r1 / $r2, $t);
  332.         }
  333.     } else {
  334.         my ($d, $x2, $y2);
  335.         if ($inverted) {
  336.         ($x2, $y2) = @{$z1->cartesian};
  337.         $d = $x2*$x2 + $y2*$y2;
  338.         _divbyzero "$z2/0" if $d == 0;
  339.         return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
  340.         } else {
  341.         my ($x1, $y1) = @{$z1->cartesian};
  342.         if (ref $z2) {
  343.             ($x2, $y2) = @{$z2->cartesian};
  344.             $d = $x2*$x2 + $y2*$y2;
  345.             _divbyzero "$z1/0" if $d == 0;
  346.             my $u = ($x1*$x2 + $y1*$y2)/$d;
  347.             my $v = ($y1*$x2 - $x1*$y2)/$d;
  348.             return (ref $z1)->make($u, $v);
  349.         } else {
  350.             _divbyzero "$z1/0" if $z2 == 0;
  351.             return (ref $z1)->make($x1/$z2, $y1/$z2);
  352.         }
  353.         }
  354.     }
  355. }
  356.  
  357. #
  358. # _zerotozero
  359. #
  360. # Die on zero raised to the zeroth.
  361. #
  362. sub _zerotozero {
  363.     my $mess = "The zero raised to the zeroth power is not defined.\n";
  364.  
  365.     my @up = caller(1);
  366.  
  367.     $mess .= "Died at $up[1] line $up[2].\n";
  368.  
  369.     die $mess;
  370. }
  371.  
  372. #
  373. # (power)
  374. #
  375. # Computes z1**z2 = exp(z2 * log z1)).
  376. #
  377. sub power {
  378.     my ($z1, $z2, $inverted) = @_;
  379.     my $z1z = $z1 == 0;
  380.     my $z2z = $z2 == 0;
  381.     _zerotozero if ($z1z and $z2z);
  382.     if ($inverted) {
  383.         return 0 if ($z2z);
  384.         return 1 if ($z1z or $z2 == 1);
  385.     } else {
  386.         return 0 if ($z1z);
  387.         return 1 if ($z2z or $z1 == 1);
  388.     }
  389.     return $inverted ? exp($z1 * log $z2) : exp($z2 * log $z1);
  390. }
  391.  
  392. #
  393. # (spaceship)
  394. #
  395. # Computes z1 <=> z2.
  396. # Sorts on the real part first, then on the imaginary part. Thus 2-4i > 3+8i.
  397. #
  398. sub spaceship {
  399.     my ($z1, $z2, $inverted) = @_;
  400.     my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
  401.     my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
  402.     my $sgn = $inverted ? -1 : 1;
  403.     return $sgn * ($re1 <=> $re2) if $re1 != $re2;
  404.     return $sgn * ($im1 <=> $im2);
  405. }
  406.  
  407. #
  408. # (negate)
  409. #
  410. # Computes -z.
  411. #
  412. sub negate {
  413.     my ($z) = @_;
  414.     if ($z->{c_dirty}) {
  415.         my ($r, $t) = @{$z->polar};
  416.         $t = ($t <= 0) ? $t + pi : $t - pi;
  417.         return (ref $z)->emake($r, $t);
  418.     }
  419.     my ($re, $im) = @{$z->cartesian};
  420.     return (ref $z)->make(-$re, -$im);
  421. }
  422.  
  423. #
  424. # (conjugate)
  425. #
  426. # Compute complex's conjugate.
  427. #
  428. sub conjugate {
  429.     my ($z) = @_;
  430.     if ($z->{c_dirty}) {
  431.         my ($r, $t) = @{$z->polar};
  432.         return (ref $z)->emake($r, -$t);
  433.     }
  434.     my ($re, $im) = @{$z->cartesian};
  435.     return (ref $z)->make($re, -$im);
  436. }
  437.  
  438. #
  439. # (abs)
  440. #
  441. # Compute complex's norm (rho).
  442. #
  443. sub abs {
  444.     my ($z) = @_;
  445.     my ($r, $t) = @{$z->polar};
  446.     return $r;
  447. }
  448.  
  449. #
  450. # arg
  451. #
  452. # Compute complex's argument (theta).
  453. #
  454. sub arg {
  455.     my ($z) = @_;
  456.     return ($z < 0 ? pi : 0) unless ref $z;
  457.     my ($r, $t) = @{$z->polar};
  458.     if    ($t >   pi()) { $t -= pit2 }
  459.     elsif ($t <= -pi()) { $t += pit2 }
  460.     return $t;
  461. }
  462.  
  463. #
  464. # (sqrt)
  465. #
  466. # Compute sqrt(z).
  467. #
  468. sub sqrt {
  469.     my ($z) = @_;
  470.     return $z >= 0 ? sqrt($z) : cplx(0, sqrt(-$z)) unless ref $z;
  471.     my ($re, $im) = @{$z->cartesian};
  472.     return cplx($re < 0 ? (0, sqrt(-$re)) : (sqrt($re), 0)) if $im == 0;
  473.     my ($r, $t) = @{$z->polar};
  474.     return (ref $z)->emake(sqrt($r), $t/2);
  475. }
  476.  
  477. #
  478. # cbrt
  479. #
  480. # Compute cbrt(z) (cubic root).
  481. #
  482. sub cbrt {
  483.     my ($z) = @_;
  484.     return $z < 0 ? -exp(log(-$z)/3) : ($z > 0 ? exp(log($z)/3): 0)
  485.         unless ref $z;
  486.     my ($r, $t) = @{$z->polar};
  487.     return (ref $z)->emake(exp(log($r)/3), $t/3);
  488. }
  489.  
  490. #
  491. # _rootbad
  492. #
  493. # Die on bad root.
  494. #
  495. sub _rootbad {
  496.     my $mess = "Root $_[0] not defined, root must be positive integer.\n";
  497.  
  498.     my @up = caller(1);
  499.  
  500.     $mess .= "Died at $up[1] line $up[2].\n";
  501.  
  502.     die $mess;
  503. }
  504.  
  505. #
  506. # root
  507. #
  508. # Computes all nth root for z, returning an array whose size is n.
  509. # `n' must be a positive integer.
  510. #
  511. # The roots are given by (for k = 0..n-1):
  512. #
  513. # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
  514. #
  515. sub root {
  516.     my ($z, $n) = @_;
  517.     _rootbad($n) if ($n < 1 or int($n) != $n);
  518.     my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi);
  519.     my @root;
  520.     my $k;
  521.     my $theta_inc = pit2 / $n;
  522.     my $rho = $r ** (1/$n);
  523.     my $theta;
  524.     my $complex = ref($z) || $package;
  525.     for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
  526.         push(@root, $complex->emake($rho, $theta));
  527.     }
  528.     return @root;
  529. }
  530.  
  531. #
  532. # Re
  533. #
  534. # Return Re(z).
  535. #
  536. sub Re {
  537.     my ($z) = @_;
  538.     return $z unless ref $z;
  539.     my ($re, $im) = @{$z->cartesian};
  540.     return $re;
  541. }
  542.  
  543. #
  544. # Im
  545. #
  546. # Return Im(z).
  547. #
  548. sub Im {
  549.     my ($z) = @_;
  550.     return 0 unless ref $z;
  551.     my ($re, $im) = @{$z->cartesian};
  552.     return $im;
  553. }
  554.  
  555. #
  556. # (exp)
  557. #
  558. # Computes exp(z).
  559. #
  560. sub exp {
  561.     my ($z) = @_;
  562.     my ($x, $y) = @{$z->cartesian};
  563.     return (ref $z)->emake(exp($x), $y);
  564. }
  565.  
  566. #
  567. # _logofzero
  568. #
  569. # Die on logarithm of zero.
  570. #
  571. sub _logofzero {
  572.     my $mess = "$_[0]: Logarithm of zero.\n";
  573.  
  574.     if (defined $_[1]) {
  575.     $mess .= "(Because in the definition of $_[0], the argument ";
  576.     $mess .= "$_[1] " unless ($_[1] eq '0');
  577.     $mess .= "is 0)\n";
  578.     }
  579.  
  580.     my @up = caller(1);
  581.  
  582.     $mess .= "Died at $up[1] line $up[2].\n";
  583.  
  584.     die $mess;
  585. }
  586.  
  587. #
  588. # (log)
  589. #
  590. # Compute log(z).
  591. #
  592. sub log {
  593.     my ($z) = @_;
  594.     unless (ref $z) {
  595.         _logofzero("log") if $z == 0;
  596.         return $z > 0 ? log($z) : cplx(log(-$z), pi);
  597.     }
  598.     my ($r, $t) = @{$z->polar};
  599.     _logofzero("log") if $r == 0;
  600.     if    ($t >   pi()) { $t -= pit2 }
  601.     elsif ($t <= -pi()) { $t += pit2 }
  602.     return (ref $z)->make(log($r), $t);
  603. }
  604.  
  605. #
  606. # ln
  607. #
  608. # Alias for log().
  609. #
  610. sub ln { Math::Complex::log(@_) }
  611.  
  612. #
  613. # log10
  614. #
  615. # Compute log10(z).
  616. #
  617.  
  618. sub log10 {
  619.     return Math::Complex::log($_[0]) * uplog10;
  620. }
  621.  
  622. #
  623. # logn
  624. #
  625. # Compute logn(z,n) = log(z) / log(n)
  626. #
  627. sub logn {
  628.     my ($z, $n) = @_;
  629.     $z = cplx($z, 0) unless ref $z;
  630.     my $logn = $logn{$n};
  631.     $logn = $logn{$n} = log($n) unless defined $logn;    # Cache log(n)
  632.     return log($z) / $logn;
  633. }
  634.  
  635. #
  636. # (cos)
  637. #
  638. # Compute cos(z) = (exp(iz) + exp(-iz))/2.
  639. #
  640. sub cos {
  641.     my ($z) = @_;
  642.     my ($x, $y) = @{$z->cartesian};
  643.     my $ey = exp($y);
  644.     my $ey_1 = 1 / $ey;
  645.     return (ref $z)->make(cos($x) * ($ey + $ey_1)/2,
  646.                   sin($x) * ($ey_1 - $ey)/2);
  647. }
  648.  
  649. #
  650. # (sin)
  651. #
  652. # Compute sin(z) = (exp(iz) - exp(-iz))/2.
  653. #
  654. sub sin {
  655.     my ($z) = @_;
  656.     my ($x, $y) = @{$z->cartesian};
  657.     my $ey = exp($y);
  658.     my $ey_1 = 1 / $ey;
  659.     return (ref $z)->make(sin($x) * ($ey + $ey_1)/2,
  660.                   cos($x) * ($ey - $ey_1)/2);
  661. }
  662.  
  663. #
  664. # tan
  665. #
  666. # Compute tan(z) = sin(z) / cos(z).
  667. #
  668. sub tan {
  669.     my ($z) = @_;
  670.     my $cz = cos($z);
  671.     _divbyzero "tan($z)", "cos($z)" if ($cz == 0);
  672.     return sin($z) / $cz;
  673. }
  674.  
  675. #
  676. # sec
  677. #
  678. # Computes the secant sec(z) = 1 / cos(z).
  679. #
  680. sub sec {
  681.     my ($z) = @_;
  682.     my $cz = cos($z);
  683.     _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
  684.     return 1 / $cz;
  685. }
  686.  
  687. #
  688. # csc
  689. #
  690. # Computes the cosecant csc(z) = 1 / sin(z).
  691. #
  692. sub csc {
  693.     my ($z) = @_;
  694.     my $sz = sin($z);
  695.     _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
  696.     return 1 / $sz;
  697. }
  698.  
  699. #
  700. # cosec
  701. #
  702. # Alias for csc().
  703. #
  704. sub cosec { Math::Complex::csc(@_) }
  705.  
  706. #
  707. # cot
  708. #
  709. # Computes cot(z) = cos(z) / sin(z).
  710. #
  711. sub cot {
  712.     my ($z) = @_;
  713.     my $sz = sin($z);
  714.     _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
  715.     return cos($z) / $sz;
  716. }
  717.  
  718. #
  719. # cotan
  720. #
  721. # Alias for cot().
  722. #
  723. sub cotan { Math::Complex::cot(@_) }
  724.  
  725. #
  726. # acos
  727. #
  728. # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
  729. #
  730. sub acos {
  731.     my $z = $_[0];
  732.     return atan2(sqrt(1-$z*$z), $z) if (! ref $z) && abs($z) <= 1;
  733.     my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
  734.     my $t1 = sqrt(($x+1)*($x+1) + $y*$y);
  735.     my $t2 = sqrt(($x-1)*($x-1) + $y*$y);
  736.     my $alpha = ($t1 + $t2)/2;
  737.     my $beta  = ($t1 - $t2)/2;
  738.     $alpha = 1 if $alpha < 1;
  739.     if    ($beta >  1) { $beta =  1 }
  740.     elsif ($beta < -1) { $beta = -1 }
  741.     my $u = atan2(sqrt(1-$beta*$beta), $beta);
  742.     my $v = log($alpha + sqrt($alpha*$alpha-1));
  743.     $v = -$v if $y > 0 || ($y == 0 && $x < -1);
  744.     return $package->make($u, $v);
  745. }
  746.  
  747. #
  748. # asin
  749. #
  750. # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
  751. #
  752. sub asin {
  753.     my $z = $_[0];
  754.     return atan2($z, sqrt(1-$z*$z)) if (! ref $z) && abs($z) <= 1;
  755.     my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
  756.     my $t1 = sqrt(($x+1)*($x+1) + $y*$y);
  757.     my $t2 = sqrt(($x-1)*($x-1) + $y*$y);
  758.     my $alpha = ($t1 + $t2)/2;
  759.     my $beta  = ($t1 - $t2)/2;
  760.     $alpha = 1 if $alpha < 1;
  761.     if    ($beta >  1) { $beta =  1 }
  762.     elsif ($beta < -1) { $beta = -1 }
  763.     my $u =  atan2($beta, sqrt(1-$beta*$beta));
  764.     my $v = -log($alpha + sqrt($alpha*$alpha-1));
  765.     $v = -$v if $y > 0 || ($y == 0 && $x < -1);
  766.     return $package->make($u, $v);
  767. }
  768.  
  769. #
  770. # atan
  771. #
  772. # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
  773. #
  774. sub atan {
  775.     my ($z) = @_;
  776.     return atan2($z, 1) unless ref $z;
  777.     _divbyzero "atan(i)"  if ( $z == i);
  778.     _divbyzero "atan(-i)" if (-$z == i);
  779.     my $log = log((i + $z) / (i - $z));
  780.     $ip2 = 0.5 * i unless defined $ip2;
  781.     return $ip2 * $log;
  782. }
  783.  
  784. #
  785. # asec
  786. #
  787. # Computes the arc secant asec(z) = acos(1 / z).
  788. #
  789. sub asec {
  790.     my ($z) = @_;
  791.     _divbyzero "asec($z)", $z if ($z == 0);
  792.     return acos(1 / $z);
  793. }
  794.  
  795. #
  796. # acsc
  797. #
  798. # Computes the arc cosecant acsc(z) = asin(1 / z).
  799. #
  800. sub acsc {
  801.     my ($z) = @_;
  802.     _divbyzero "acsc($z)", $z if ($z == 0);
  803.     return asin(1 / $z);
  804. }
  805.  
  806. #
  807. # acosec
  808. #
  809. # Alias for acsc().
  810. #
  811. sub acosec { Math::Complex::acsc(@_) }
  812.  
  813. #
  814. # acot
  815. #
  816. # Computes the arc cotangent acot(z) = atan(1 / z)
  817. #
  818. sub acot {
  819.     my ($z) = @_;
  820.     return ($z >= 0) ? atan2(1, $z) : atan2(-1, -$z) unless ref $z;
  821.     _divbyzero "acot(i)", if ( $z == i);
  822.     _divbyzero "acot(-i)" if (-$z == i);
  823.     return atan(1 / $z);
  824. }
  825.  
  826. #
  827. # acotan
  828. #
  829. # Alias for acot().
  830. #
  831. sub acotan { Math::Complex::acot(@_) }
  832.  
  833. #
  834. # cosh
  835. #
  836. # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
  837. #
  838. sub cosh {
  839.     my ($z) = @_;
  840.     my $ex;
  841.     unless (ref $z) {
  842.         $ex = exp($z);
  843.         return ($ex + 1/$ex)/2;
  844.     }
  845.     my ($x, $y) = @{$z->cartesian};
  846.     $ex = exp($x);
  847.     my $ex_1 = 1 / $ex;
  848.     return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
  849.                   sin($y) * ($ex - $ex_1)/2);
  850. }
  851.  
  852. #
  853. # sinh
  854. #
  855. # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
  856. #
  857. sub sinh {
  858.     my ($z) = @_;
  859.     my $ex;
  860.     unless (ref $z) {
  861.         $ex = exp($z);
  862.         return ($ex - 1/$ex)/2;
  863.     }
  864.     my ($x, $y) = @{$z->cartesian};
  865.     $ex = exp($x);
  866.     my $ex_1 = 1 / $ex;
  867.     return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
  868.                   sin($y) * ($ex + $ex_1)/2);
  869. }
  870.  
  871. #
  872. # tanh
  873. #
  874. # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
  875. #
  876. sub tanh {
  877.     my ($z) = @_;
  878.     my $cz = cosh($z);
  879.     _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
  880.     return sinh($z) / $cz;
  881. }
  882.  
  883. #
  884. # sech
  885. #
  886. # Computes the hyperbolic secant sech(z) = 1 / cosh(z).
  887. #
  888. sub sech {
  889.     my ($z) = @_;
  890.     my $cz = cosh($z);
  891.     _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
  892.     return 1 / $cz;
  893. }
  894.  
  895. #
  896. # csch
  897. #
  898. # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
  899. #
  900. sub csch {
  901.     my ($z) = @_;
  902.     my $sz = sinh($z);
  903.     _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
  904.     return 1 / $sz;
  905. }
  906.  
  907. #
  908. # cosech
  909. #
  910. # Alias for csch().
  911. #
  912. sub cosech { Math::Complex::csch(@_) }
  913.  
  914. #
  915. # coth
  916. #
  917. # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
  918. #
  919. sub coth {
  920.     my ($z) = @_;
  921.     my $sz = sinh($z);
  922.     _divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
  923.     return cosh($z) / $sz;
  924. }
  925.  
  926. #
  927. # cotanh
  928. #
  929. # Alias for coth().
  930. #
  931. sub cotanh { Math::Complex::coth(@_) }
  932.  
  933. #
  934. # acosh
  935. #
  936. # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
  937. #
  938. sub acosh {
  939.     my ($z) = @_;
  940.     unless (ref $z) {
  941.         return log($z + sqrt($z*$z-1)) if $z >= 1;
  942.         $z = cplx($z, 0);
  943.     }
  944.     my ($re, $im) = @{$z->cartesian};
  945.     if ($im == 0) {
  946.         return cplx(log($re + sqrt($re*$re - 1)), 0) if $re >= 1;
  947.         return cplx(0, atan2(sqrt(1-$re*$re), $re)) if abs($re) <= 1;
  948.     }
  949.     return log($z + sqrt($z*$z - 1));
  950. }
  951.  
  952. #
  953. # asinh
  954. #
  955. # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
  956. #
  957. sub asinh {
  958.     my ($z) = @_;
  959.     return log($z + sqrt($z*$z + 1));
  960. }
  961.  
  962. #
  963. # atanh
  964. #
  965. # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
  966. #
  967. sub atanh {
  968.     my ($z) = @_;
  969.     unless (ref $z) {
  970.         return log((1 + $z)/(1 - $z))/2 if abs($z) < 1;
  971.         $z = cplx($z, 0);
  972.     }
  973.     _divbyzero 'atanh(1)',  "1 - $z" if ($z ==  1);
  974.     _logofzero 'atanh(-1)'           if ($z == -1);
  975.     return 0.5 * log((1 + $z) / (1 - $z));
  976. }
  977.  
  978. #
  979. # asech
  980. #
  981. # Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
  982. #
  983. sub asech {
  984.     my ($z) = @_;
  985.     _divbyzero 'asech(0)', $z if ($z == 0);
  986.     return acosh(1 / $z);
  987. }
  988.  
  989. #
  990. # acsch
  991. #
  992. # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
  993. #
  994. sub acsch {
  995.     my ($z) = @_;
  996.     _divbyzero 'acsch(0)', $z if ($z == 0);
  997.     return asinh(1 / $z);
  998. }
  999.  
  1000. #
  1001. # acosech
  1002. #
  1003. # Alias for acosh().
  1004. #
  1005. sub acosech { Math::Complex::acsch(@_) }
  1006.  
  1007. #
  1008. # acoth
  1009. #
  1010. # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
  1011. #
  1012. sub acoth {
  1013.     my ($z) = @_;
  1014.     unless (ref $z) {
  1015.         return log(($z + 1)/($z - 1))/2 if abs($z) > 1;
  1016.         $z = cplx($z, 0);
  1017.     }
  1018.     _divbyzero 'acoth(1)', "$z - 1" if ($z ==  1);
  1019.     _logofzero 'acoth(-1)'          if ($z == -1);
  1020.     return log((1 + $z) / ($z - 1)) / 2;
  1021. }
  1022.  
  1023. #
  1024. # acotanh
  1025. #
  1026. # Alias for acot().
  1027. #
  1028. sub acotanh { Math::Complex::acoth(@_) }
  1029.  
  1030. #
  1031. # (atan2)
  1032. #
  1033. # Compute atan(z1/z2).
  1034. #
  1035. sub atan2 {
  1036.     my ($z1, $z2, $inverted) = @_;
  1037.     my ($re1, $im1, $re2, $im2);
  1038.     if ($inverted) {
  1039.         ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
  1040.         ($re2, $im2) = @{$z1->cartesian};
  1041.     } else {
  1042.         ($re1, $im1) = @{$z1->cartesian};
  1043.         ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
  1044.     }
  1045.     if ($im2 == 0) {
  1046.         return cplx(atan2($re1, $re2), 0) if $im1 == 0;
  1047.         return cplx(($im1<=>0) * pip2, 0) if $re2 == 0;
  1048.     }
  1049.     my $w = atan($z1/$z2);
  1050.     my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
  1051.     $u += pi   if $re2 < 0;
  1052.     $u -= pit2 if $u > pi;
  1053.     return cplx($u, $v);
  1054. }
  1055.  
  1056. #
  1057. # display_format
  1058. # ->display_format
  1059. #
  1060. # Set (fetch if no argument) display format for all complex numbers that
  1061. # don't happen to have overridden it via ->display_format
  1062. #
  1063. # When called as a method, this actually sets the display format for
  1064. # the current object.
  1065. #
  1066. # Valid object formats are 'c' and 'p' for cartesian and polar. The first
  1067. # letter is used actually, so the type can be fully spelled out for clarity.
  1068. #
  1069. sub display_format {
  1070.     my $self = shift;
  1071.     my $format = undef;
  1072.  
  1073.     if (ref $self) {            # Called as a method
  1074.         $format = shift;
  1075.     } else {                # Regular procedure call
  1076.         $format = $self;
  1077.         undef $self;
  1078.     }
  1079.  
  1080.     if (defined $self) {
  1081.         return defined $self->{display} ? $self->{display} : $display
  1082.             unless defined $format;
  1083.         return $self->{display} = $format;
  1084.     }
  1085.  
  1086.     return $display unless defined $format;
  1087.     return $display = $format;
  1088. }
  1089.  
  1090. #
  1091. # (stringify)
  1092. #
  1093. # Show nicely formatted complex number under its cartesian or polar form,
  1094. # depending on the current display format:
  1095. #
  1096. # . If a specific display format has been recorded for this object, use it.
  1097. # . Otherwise, use the generic current default for all complex numbers,
  1098. #   which is a package global variable.
  1099. #
  1100. sub stringify {
  1101.     my ($z) = shift;
  1102.     my $format;
  1103.  
  1104.     $format = $display;
  1105.     $format = $z->{display} if defined $z->{display};
  1106.  
  1107.     return $z->stringify_polar if $format =~ /^p/i;
  1108.     return $z->stringify_cartesian;
  1109. }
  1110.  
  1111. #
  1112. # ->stringify_cartesian
  1113. #
  1114. # Stringify as a cartesian representation 'a+bi'.
  1115. #
  1116. sub stringify_cartesian {
  1117.     my $z  = shift;
  1118.     my ($x, $y) = @{$z->cartesian};
  1119.     my ($re, $im);
  1120.     my $eps = 1e-14;
  1121.  
  1122.     $x = int($x + ($x < 0 ? -1 : 1) * $eps)
  1123.         if int(abs($x)) != int(abs($x) + $eps);
  1124.     $y = int($y + ($y < 0 ? -1 : 1) * $eps)
  1125.         if int(abs($y)) != int(abs($y) + $eps);
  1126.  
  1127.     $re = "$x" if abs($x) >= $eps;
  1128.         if ($y == 1)                           { $im = 'i' }
  1129.         elsif ($y == -1)                       { $im = '-i' }
  1130.         elsif (abs($y) >= $eps)                { $im = $y . "i" }
  1131.  
  1132.     my $str = '';
  1133.     $str = $re if defined $re;
  1134.     $str .= "+$im" if defined $im;
  1135.     $str =~ s/\+-/-/;
  1136.     $str =~ s/^\+//;
  1137.     $str = '0' unless $str;
  1138.  
  1139.     return $str;
  1140. }
  1141.  
  1142. #
  1143. # ->stringify_polar
  1144. #
  1145. # Stringify as a polar representation '[r,t]'.
  1146. #
  1147. sub stringify_polar {
  1148.     my $z  = shift;
  1149.     my ($r, $t) = @{$z->polar};
  1150.     my $theta;
  1151.     my $eps = 1e-14;
  1152.  
  1153.     return '[0,0]' if $r <= $eps;
  1154.  
  1155.     my $nt = $t / pit2;
  1156.     $nt = ($nt - int($nt)) * pit2;
  1157.     $nt += pit2 if $nt < 0;            # Range [0, 2pi]
  1158.  
  1159.     if (abs($nt) <= $eps)        { $theta = 0 }
  1160.     elsif (abs(pi-$nt) <= $eps)    { $theta = 'pi' }
  1161.  
  1162.     if (defined $theta) {
  1163.         $r = int($r + ($r < 0 ? -1 : 1) * $eps)
  1164.             if int(abs($r)) != int(abs($r) + $eps);
  1165.         $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
  1166.             if ($theta ne 'pi' and
  1167.                 int(abs($theta)) != int(abs($theta) + $eps));
  1168.         return "\[$r,$theta\]";
  1169.     }
  1170.  
  1171.     #
  1172.     # Okay, number is not a real. Try to identify pi/n and friends...
  1173.     #
  1174.  
  1175.     $nt -= pit2 if $nt > pi;
  1176.     my ($n, $k, $kpi);
  1177.  
  1178.     for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
  1179.         $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
  1180.         if (abs($kpi/$n - $nt) <= $eps) {
  1181.             $theta = ($nt < 0 ? '-':'').
  1182.                  ($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
  1183.             last;
  1184.         }
  1185.     }
  1186.  
  1187.     $theta = $nt unless defined $theta;
  1188.  
  1189.     $r = int($r + ($r < 0 ? -1 : 1) * $eps)
  1190.         if int(abs($r)) != int(abs($r) + $eps);
  1191.     $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
  1192.         if ($theta !~ m(^-?\d*pi/\d+$) and
  1193.             int(abs($theta)) != int(abs($theta) + $eps));
  1194.  
  1195.     return "\[$r,$theta\]";
  1196. }
  1197.  
  1198. 1;
  1199. __END__
  1200.  
  1201. =head1 NAME
  1202.  
  1203. Math::Complex - complex numbers and associated mathematical functions
  1204.  
  1205. =head1 SYNOPSIS
  1206.  
  1207.     use Math::Complex;
  1208.  
  1209.     $z = Math::Complex->make(5, 6);
  1210.     $t = 4 - 3*i + $z;
  1211.     $j = cplxe(1, 2*pi/3);
  1212.  
  1213. =head1 DESCRIPTION
  1214.  
  1215. This package lets you create and manipulate complex numbers. By default,
  1216. I<Perl> limits itself to real numbers, but an extra C<use> statement brings
  1217. full complex support, along with a full set of mathematical functions
  1218. typically associated with and/or extended to complex numbers.
  1219.  
  1220. If you wonder what complex numbers are, they were invented to be able to solve
  1221. the following equation:
  1222.  
  1223.     x*x = -1
  1224.  
  1225. and by definition, the solution is noted I<i> (engineers use I<j> instead since
  1226. I<i> usually denotes an intensity, but the name does not matter). The number
  1227. I<i> is a pure I<imaginary> number.
  1228.  
  1229. The arithmetics with pure imaginary numbers works just like you would expect
  1230. it with real numbers... you just have to remember that
  1231.  
  1232.     i*i = -1
  1233.  
  1234. so you have:
  1235.  
  1236.     5i + 7i = i * (5 + 7) = 12i
  1237.     4i - 3i = i * (4 - 3) = i
  1238.     4i * 2i = -8
  1239.     6i / 2i = 3
  1240.     1 / i = -i
  1241.  
  1242. Complex numbers are numbers that have both a real part and an imaginary
  1243. part, and are usually noted:
  1244.  
  1245.     a + bi
  1246.  
  1247. where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
  1248. arithmetic with complex numbers is straightforward. You have to
  1249. keep track of the real and the imaginary parts, but otherwise the
  1250. rules used for real numbers just apply:
  1251.  
  1252.     (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
  1253.     (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
  1254.  
  1255. A graphical representation of complex numbers is possible in a plane
  1256. (also called the I<complex plane>, but it's really a 2D plane).
  1257. The number
  1258.  
  1259.     z = a + bi
  1260.  
  1261. is the point whose coordinates are (a, b). Actually, it would
  1262. be the vector originating from (0, 0) to (a, b). It follows that the addition
  1263. of two complex numbers is a vectorial addition.
  1264.  
  1265. Since there is a bijection between a point in the 2D plane and a complex
  1266. number (i.e. the mapping is unique and reciprocal), a complex number
  1267. can also be uniquely identified with polar coordinates:
  1268.  
  1269.     [rho, theta]
  1270.  
  1271. where C<rho> is the distance to the origin, and C<theta> the angle between
  1272. the vector and the I<x> axis. There is a notation for this using the
  1273. exponential form, which is:
  1274.  
  1275.     rho * exp(i * theta)
  1276.  
  1277. where I<i> is the famous imaginary number introduced above. Conversion
  1278. between this form and the cartesian form C<a + bi> is immediate:
  1279.  
  1280.     a = rho * cos(theta)
  1281.     b = rho * sin(theta)
  1282.  
  1283. which is also expressed by this formula:
  1284.  
  1285.     z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
  1286.  
  1287. In other words, it's the projection of the vector onto the I<x> and I<y>
  1288. axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
  1289. the I<argument> of the complex number. The I<norm> of C<z> will be
  1290. noted C<abs(z)>.
  1291.  
  1292. The polar notation (also known as the trigonometric
  1293. representation) is much more handy for performing multiplications and
  1294. divisions of complex numbers, whilst the cartesian notation is better
  1295. suited for additions and subtractions. Real numbers are on the I<x>
  1296. axis, and therefore I<theta> is zero or I<pi>.
  1297.  
  1298. All the common operations that can be performed on a real number have
  1299. been defined to work on complex numbers as well, and are merely
  1300. I<extensions> of the operations defined on real numbers. This means
  1301. they keep their natural meaning when there is no imaginary part, provided
  1302. the number is within their definition set.
  1303.  
  1304. For instance, the C<sqrt> routine which computes the square root of
  1305. its argument is only defined for non-negative real numbers and yields a
  1306. non-negative real number (it is an application from B<R+> to B<R+>).
  1307. If we allow it to return a complex number, then it can be extended to
  1308. negative real numbers to become an application from B<R> to B<C> (the
  1309. set of complex numbers):
  1310.  
  1311.     sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
  1312.  
  1313. It can also be extended to be an application from B<C> to B<C>,
  1314. whilst its restriction to B<R> behaves as defined above by using
  1315. the following definition:
  1316.  
  1317.     sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
  1318.  
  1319. Indeed, a negative real number can be noted C<[x,pi]> (the modulus
  1320. I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
  1321. number) and the above definition states that
  1322.  
  1323.     sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
  1324.  
  1325. which is exactly what we had defined for negative real numbers above.
  1326.  
  1327. All the common mathematical functions defined on real numbers that
  1328. are extended to complex numbers share that same property of working
  1329. I<as usual> when the imaginary part is zero (otherwise, it would not
  1330. be called an extension, would it?).
  1331.  
  1332. A I<new> operation possible on a complex number that is
  1333. the identity for real numbers is called the I<conjugate>, and is noted
  1334. with an horizontal bar above the number, or C<~z> here.
  1335.  
  1336.      z = a + bi
  1337.     ~z = a - bi
  1338.  
  1339. Simple... Now look:
  1340.  
  1341.     z * ~z = (a + bi) * (a - bi) = a*a + b*b
  1342.  
  1343. We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
  1344. distance to the origin, also known as:
  1345.  
  1346.     rho = abs(z) = sqrt(a*a + b*b)
  1347.  
  1348. so
  1349.  
  1350.     z * ~z = abs(z) ** 2
  1351.  
  1352. If z is a pure real number (i.e. C<b == 0>), then the above yields:
  1353.  
  1354.     a * a = abs(a) ** 2
  1355.  
  1356. which is true (C<abs> has the regular meaning for real number, i.e. stands
  1357. for the absolute value). This example explains why the norm of C<z> is
  1358. noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
  1359. is the regular C<abs> we know when the complex number actually has no
  1360. imaginary part... This justifies I<a posteriori> our use of the C<abs>
  1361. notation for the norm.
  1362.  
  1363. =head1 OPERATIONS
  1364.  
  1365. Given the following notations:
  1366.  
  1367.     z1 = a + bi = r1 * exp(i * t1)
  1368.     z2 = c + di = r2 * exp(i * t2)
  1369.     z = <any complex or real number>
  1370.  
  1371. the following (overloaded) operations are supported on complex numbers:
  1372.  
  1373.     z1 + z2 = (a + c) + i(b + d)
  1374.     z1 - z2 = (a - c) + i(b - d)
  1375.     z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
  1376.     z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
  1377.     z1 ** z2 = exp(z2 * log z1)
  1378.     ~z1 = a - bi
  1379.     abs(z1) = r1 = sqrt(a*a + b*b)
  1380.     sqrt(z1) = sqrt(r1) * exp(i * t1/2)
  1381.     exp(z1) = exp(a) * exp(i * b)
  1382.     log(z1) = log(r1) + i*t1
  1383.     sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1))
  1384.     cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1))
  1385.     atan2(z1, z2) = atan(z1/z2)
  1386.  
  1387. The following extra operations are supported on both real and complex
  1388. numbers:
  1389.  
  1390.     Re(z) = a
  1391.     Im(z) = b
  1392.     arg(z) = t
  1393.  
  1394.     cbrt(z) = z ** (1/3)
  1395.     log10(z) = log(z) / log(10)
  1396.     logn(z, n) = log(z) / log(n)
  1397.  
  1398.     tan(z) = sin(z) / cos(z)
  1399.  
  1400.     csc(z) = 1 / sin(z)
  1401.     sec(z) = 1 / cos(z)
  1402.     cot(z) = 1 / tan(z)
  1403.  
  1404.     asin(z) = -i * log(i*z + sqrt(1-z*z))
  1405.     acos(z) = -i * log(z + i*sqrt(1-z*z))
  1406.     atan(z) = i/2 * log((i+z) / (i-z))
  1407.  
  1408.     acsc(z) = asin(1 / z)
  1409.     asec(z) = acos(1 / z)
  1410.     acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
  1411.  
  1412.     sinh(z) = 1/2 (exp(z) - exp(-z))
  1413.     cosh(z) = 1/2 (exp(z) + exp(-z))
  1414.     tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
  1415.  
  1416.     csch(z) = 1 / sinh(z)
  1417.     sech(z) = 1 / cosh(z)
  1418.     coth(z) = 1 / tanh(z)
  1419.  
  1420.     asinh(z) = log(z + sqrt(z*z+1))
  1421.     acosh(z) = log(z + sqrt(z*z-1))
  1422.     atanh(z) = 1/2 * log((1+z) / (1-z))
  1423.  
  1424.     acsch(z) = asinh(1 / z)
  1425.     asech(z) = acosh(1 / z)
  1426.     acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
  1427.  
  1428. I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, I<coth>,
  1429. I<acosech>, I<acotanh>, have aliases I<ln>, I<cosec>, I<cotan>,
  1430. I<acosec>, I<acotan>, I<cosech>, I<cotanh>, I<acosech>, I<acotanh>,
  1431. respectively.
  1432.  
  1433. The I<root> function is available to compute all the I<n>
  1434. roots of some complex, where I<n> is a strictly positive integer.
  1435. There are exactly I<n> such roots, returned as a list. Getting the
  1436. number mathematicians call C<j> such that:
  1437.  
  1438.     1 + j + j*j = 0;
  1439.  
  1440. is a simple matter of writing:
  1441.  
  1442.     $j = ((root(1, 3))[1];
  1443.  
  1444. The I<k>th root for C<z = [r,t]> is given by:
  1445.  
  1446.     (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
  1447.  
  1448. The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
  1449. order to ensure its restriction to real numbers is conform to what you
  1450. would expect, the comparison is run on the real part of the complex
  1451. number first, and imaginary parts are compared only when the real
  1452. parts match.
  1453.  
  1454. =head1 CREATION
  1455.  
  1456. To create a complex number, use either:
  1457.  
  1458.     $z = Math::Complex->make(3, 4);
  1459.     $z = cplx(3, 4);
  1460.  
  1461. if you know the cartesian form of the number, or
  1462.  
  1463.     $z = 3 + 4*i;
  1464.  
  1465. if you like. To create a number using the polar form, use either:
  1466.  
  1467.     $z = Math::Complex->emake(5, pi/3);
  1468.     $x = cplxe(5, pi/3);
  1469.  
  1470. instead. The first argument is the modulus, the second is the angle
  1471. (in radians, the full circle is 2*pi).  (Mnemonic: C<e> is used as a
  1472. notation for complex numbers in the polar form).
  1473.  
  1474. It is possible to write:
  1475.  
  1476.     $x = cplxe(-3, pi/4);
  1477.  
  1478. but that will be silently converted into C<[3,-3pi/4]>, since the modulus
  1479. must be non-negative (it represents the distance to the origin in the complex
  1480. plane).
  1481.  
  1482. =head1 STRINGIFICATION
  1483.  
  1484. When printed, a complex number is usually shown under its cartesian
  1485. form I<a+bi>, but there are legitimate cases where the polar format
  1486. I<[r,t]> is more appropriate.
  1487.  
  1488. By calling the routine C<Math::Complex::display_format> and supplying either
  1489. C<"polar"> or C<"cartesian">, you override the default display format,
  1490. which is C<"cartesian">. Not supplying any argument returns the current
  1491. setting.
  1492.  
  1493. This default can be overridden on a per-number basis by calling the
  1494. C<display_format> method instead. As before, not supplying any argument
  1495. returns the current display format for this number. Otherwise whatever you
  1496. specify will be the new display format for I<this> particular number.
  1497.  
  1498. For instance:
  1499.  
  1500.     use Math::Complex;
  1501.  
  1502.     Math::Complex::display_format('polar');
  1503.     $j = ((root(1, 3))[1];
  1504.     print "j = $j\n";        # Prints "j = [1,2pi/3]
  1505.     $j->display_format('cartesian');
  1506.     print "j = $j\n";        # Prints "j = -0.5+0.866025403784439i"
  1507.  
  1508. The polar format attempts to emphasize arguments like I<k*pi/n>
  1509. (where I<n> is a positive integer and I<k> an integer within [-9,+9]).
  1510.  
  1511. =head1 USAGE
  1512.  
  1513. Thanks to overloading, the handling of arithmetics with complex numbers
  1514. is simple and almost transparent.
  1515.  
  1516. Here are some examples:
  1517.  
  1518.     use Math::Complex;
  1519.  
  1520.     $j = cplxe(1, 2*pi/3);    # $j ** 3 == 1
  1521.     print "j = $j, j**3 = ", $j ** 3, "\n";
  1522.     print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
  1523.  
  1524.     $z = -16 + 0*i;            # Force it to be a complex
  1525.     print "sqrt($z) = ", sqrt($z), "\n";
  1526.  
  1527.     $k = exp(i * 2*pi/3);
  1528.     print "$j - $k = ", $j - $k, "\n";
  1529.  
  1530. =head1 ERRORS DUE TO DIVISION BY ZERO
  1531.  
  1532. The division (/) and the following functions
  1533.  
  1534.     tan
  1535.     sec
  1536.     csc
  1537.     cot
  1538.     asec
  1539.     acsc
  1540.     atan
  1541.     acot
  1542.     tanh
  1543.     sech
  1544.     csch
  1545.     coth
  1546.     atanh
  1547.     asech
  1548.     acsch
  1549.     acoth
  1550.  
  1551. cannot be computed for all arguments because that would mean dividing
  1552. by zero or taking logarithm of zero. These situations cause fatal
  1553. runtime errors looking like this
  1554.  
  1555.     cot(0): Division by zero.
  1556.     (Because in the definition of cot(0), the divisor sin(0) is 0)
  1557.     Died at ...
  1558.  
  1559. or
  1560.  
  1561.     atanh(-1): Logarithm of zero.
  1562.     Died at...
  1563.  
  1564. For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
  1565. C<asech>, C<acsch>, the argument cannot be C<0> (zero).  For the
  1566. C<atanh>, C<acoth>, the argument cannot be C<1> (one).  For the
  1567. C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one).  For the
  1568. C<atan>, C<acot>, the argument cannot be C<i> (the imaginary unit).
  1569. For the C<atan>, C<acoth>, the argument cannot be C<-i> (the negative
  1570. imaginary unit).  For the C<tan>, C<sec>, C<tanh>, C<sech>, the
  1571. argument cannot be I<pi/2 + k * pi>, where I<k> is any integer.
  1572.  
  1573. =head1 BUGS
  1574.  
  1575. Saying C<use Math::Complex;> exports many mathematical routines in the
  1576. caller environment and even overrides some (C<sqrt>, C<log>).
  1577. This is construed as a feature by the Authors, actually... ;-)
  1578.  
  1579. All routines expect to be given real or complex numbers. Don't attempt to
  1580. use BigFloat, since Perl has currently no rule to disambiguate a '+'
  1581. operation (for instance) between two overloaded entities.
  1582.  
  1583. =head1 AUTHORS
  1584.  
  1585. Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>> and
  1586. Jarkko Hietaniemi <F<jhi@iki.fi>>.
  1587.  
  1588. Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.
  1589.  
  1590. =cut
  1591.  
  1592. # eof
  1593.