diff options
Diffstat (limited to 'contrib/perl5/lib/Math')
| -rw-r--r-- | contrib/perl5/lib/Math/BigFloat.pm | 327 | ||||
| -rw-r--r-- | contrib/perl5/lib/Math/BigInt.pm | 415 | ||||
| -rw-r--r-- | contrib/perl5/lib/Math/Complex.pm | 1775 | ||||
| -rw-r--r-- | contrib/perl5/lib/Math/Trig.pm | 419 |
4 files changed, 2936 insertions, 0 deletions
diff --git a/contrib/perl5/lib/Math/BigFloat.pm b/contrib/perl5/lib/Math/BigFloat.pm new file mode 100644 index 0000000..576f341 --- /dev/null +++ b/contrib/perl5/lib/Math/BigFloat.pm @@ -0,0 +1,327 @@ +package Math::BigFloat; + +use Math::BigInt; + +use Exporter; # just for use to be happy +@ISA = (Exporter); + +use overload +'+' => sub {new Math::BigFloat &fadd}, +'-' => sub {new Math::BigFloat + $_[2]? fsub($_[1],${$_[0]}) : fsub(${$_[0]},$_[1])}, +'<=>' => sub {new Math::BigFloat + $_[2]? fcmp($_[1],${$_[0]}) : fcmp(${$_[0]},$_[1])}, +'cmp' => sub {new Math::BigFloat + $_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])}, +'*' => sub {new Math::BigFloat &fmul}, +'/' => sub {new Math::BigFloat + $_[2]? scalar fdiv($_[1],${$_[0]}) : + scalar fdiv(${$_[0]},$_[1])}, +'neg' => sub {new Math::BigFloat &fneg}, +'abs' => sub {new Math::BigFloat &fabs}, + +qw( +"" stringify +0+ numify) # Order of arguments unsignificant +; + +sub new { + my ($class) = shift; + my ($foo) = fnorm(shift); + panic("Not a number initialized to Math::BigFloat") if $foo eq "NaN"; + bless \$foo, $class; +} +sub numify { 0 + "${$_[0]}" } # Not needed, additional overhead + # comparing to direct compilation based on + # stringify +sub stringify { + my $n = ${$_[0]}; + + my $minus = ($n =~ s/^([+-])// && $1 eq '-'); + $n =~ s/E//; + + $n =~ s/([-+]\d+)$//; + + my $e = $1; + my $ln = length($n); + + if ($e > 0) { + $n .= "0" x $e . '.'; + } elsif (abs($e) < $ln) { + substr($n, $ln + $e, 0) = '.'; + } else { + $n = '.' . ("0" x (abs($e) - $ln)) . $n; + } + $n = "-$n" if $minus; + + # 1 while $n =~ s/(.*\d)(\d\d\d)/$1,$2/; + + return $n; +} + +$div_scale = 40; + +# Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'. + +$rnd_mode = 'even'; + +sub fadd; sub fsub; sub fmul; sub fdiv; +sub fneg; sub fabs; sub fcmp; +sub fround; sub ffround; +sub fnorm; sub fsqrt; + +# Convert a number to canonical string form. +# Takes something that looks like a number and converts it to +# the form /^[+-]\d+E[+-]\d+$/. +sub fnorm { #(string) return fnum_str + local($_) = @_; + s/\s+//g; # strip white space + if (/^([+-]?)(\d*)(\.(\d*))?([Ee]([+-]?\d+))?$/ && "$2$4" ne '') { + &norm(($1 ? "$1$2$4" : "+$2$4"),(($4 ne '') ? $6-length($4) : $6)); + } else { + 'NaN'; + } +} + +# normalize number -- for internal use +sub norm { #(mantissa, exponent) return fnum_str + local($_, $exp) = @_; + if ($_ eq 'NaN') { + 'NaN'; + } else { + s/^([+-])0+/$1/; # strip leading zeros + if (length($_) == 1) { + '+0E+0'; + } else { + $exp += length($1) if (s/(0+)$//); # strip trailing zeros + sprintf("%sE%+ld", $_, $exp); + } + } +} + +# negation +sub fneg { #(fnum_str) return fnum_str + local($_) = fnorm($_[$[]); + vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0E+0'; # flip sign + s/^H/N/; + $_; +} + +# absolute value +sub fabs { #(fnum_str) return fnum_str + local($_) = fnorm($_[$[]); + s/^-/+/; # mash sign + $_; +} + +# multiplication +sub fmul { #(fnum_str, fnum_str) return fnum_str + local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1])); + if ($x eq 'NaN' || $y eq 'NaN') { + 'NaN'; + } else { + local($xm,$xe) = split('E',$x); + local($ym,$ye) = split('E',$y); + &norm(Math::BigInt::bmul($xm,$ym),$xe+$ye); + } +} + +# addition +sub fadd { #(fnum_str, fnum_str) return fnum_str + local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1])); + if ($x eq 'NaN' || $y eq 'NaN') { + 'NaN'; + } else { + local($xm,$xe) = split('E',$x); + local($ym,$ye) = split('E',$y); + ($xm,$xe,$ym,$ye) = ($ym,$ye,$xm,$xe) if ($xe < $ye); + &norm(Math::BigInt::badd($ym,$xm.('0' x ($xe-$ye))),$ye); + } +} + +# subtraction +sub fsub { #(fnum_str, fnum_str) return fnum_str + fadd($_[$[],fneg($_[$[+1])); +} + +# division +# args are dividend, divisor, scale (optional) +# result has at most max(scale, length(dividend), length(divisor)) digits +sub fdiv #(fnum_str, fnum_str[,scale]) return fnum_str +{ + local($x,$y,$scale) = (fnorm($_[$[]),fnorm($_[$[+1]),$_[$[+2]); + if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0E+0') { + 'NaN'; + } else { + local($xm,$xe) = split('E',$x); + local($ym,$ye) = split('E',$y); + $scale = $div_scale if (!$scale); + $scale = length($xm)-1 if (length($xm)-1 > $scale); + $scale = length($ym)-1 if (length($ym)-1 > $scale); + $scale = $scale + length($ym) - length($xm); + &norm(&round(Math::BigInt::bdiv($xm.('0' x $scale),$ym),$ym), + $xe-$ye-$scale); + } +} + +# round int $q based on fraction $r/$base using $rnd_mode +sub round { #(int_str, int_str, int_str) return int_str + local($q,$r,$base) = @_; + if ($q eq 'NaN' || $r eq 'NaN') { + 'NaN'; + } elsif ($rnd_mode eq 'trunc') { + $q; # just truncate + } else { + local($cmp) = Math::BigInt::bcmp(Math::BigInt::bmul($r,'+2'),$base); + if ( $cmp < 0 || + ($cmp == 0 && + ( $rnd_mode eq 'zero' || + ($rnd_mode eq '-inf' && (substr($q,$[,1) eq '+')) || + ($rnd_mode eq '+inf' && (substr($q,$[,1) eq '-')) || + ($rnd_mode eq 'even' && $q =~ /[24680]$/) || + ($rnd_mode eq 'odd' && $q =~ /[13579]$/) )) ) { + $q; # round down + } else { + Math::BigInt::badd($q, ((substr($q,$[,1) eq '-') ? '-1' : '+1')); + # round up + } + } +} + +# round the mantissa of $x to $scale digits +sub fround { #(fnum_str, scale) return fnum_str + local($x,$scale) = (fnorm($_[$[]),$_[$[+1]); + if ($x eq 'NaN' || $scale <= 0) { + $x; + } else { + local($xm,$xe) = split('E',$x); + if (length($xm)-1 <= $scale) { + $x; + } else { + &norm(&round(substr($xm,$[,$scale+1), + "+0".substr($xm,$[+$scale+1,1),"+10"), + $xe+length($xm)-$scale-1); + } + } +} + +# round $x at the 10 to the $scale digit place +sub ffround { #(fnum_str, scale) return fnum_str + local($x,$scale) = (fnorm($_[$[]),$_[$[+1]); + if ($x eq 'NaN') { + 'NaN'; + } else { + local($xm,$xe) = split('E',$x); + if ($xe >= $scale) { + $x; + } else { + $xe = length($xm)+$xe-$scale; + if ($xe < 1) { + '+0E+0'; + } elsif ($xe == 1) { + &norm(&round('+0',"+0".substr($xm,$[+1,1),"+10"), $scale); + } else { + &norm(&round(substr($xm,$[,$xe), + "+0".substr($xm,$[+$xe,1),"+10"), $scale); + } + } + } +} + +# compare 2 values returns one of undef, <0, =0, >0 +# returns undef if either or both input value are not numbers +sub fcmp #(fnum_str, fnum_str) return cond_code +{ + local($x, $y) = (fnorm($_[$[]),fnorm($_[$[+1])); + if ($x eq "NaN" || $y eq "NaN") { + undef; + } else { + ord($y) <=> ord($x) + || + ( local($xm,$xe,$ym,$ye) = split('E', $x."E$y"), + (($xe <=> $ye) * (substr($x,$[,1).'1') + || Math::BigInt::cmp($xm,$ym)) + ); + } +} + +# square root by Newtons method. +sub fsqrt { #(fnum_str[, scale]) return fnum_str + local($x, $scale) = (fnorm($_[$[]), $_[$[+1]); + if ($x eq 'NaN' || $x =~ /^-/) { + 'NaN'; + } elsif ($x eq '+0E+0') { + '+0E+0'; + } else { + local($xm, $xe) = split('E',$x); + $scale = $div_scale if (!$scale); + $scale = length($xm)-1 if ($scale < length($xm)-1); + local($gs, $guess) = (1, sprintf("1E%+d", (length($xm)+$xe-1)/2)); + while ($gs < 2*$scale) { + $guess = fmul(fadd($guess,fdiv($x,$guess,$gs*2)),".5"); + $gs *= 2; + } + new Math::BigFloat &fround($guess, $scale); + } +} + +1; +__END__ + +=head1 NAME + +Math::BigFloat - Arbitrary length float math package + +=head1 SYNOPSIS + + use Math::BigFloat; + $f = Math::BigFloat->new($string); + + $f->fadd(NSTR) return NSTR addition + $f->fsub(NSTR) return NSTR subtraction + $f->fmul(NSTR) return NSTR multiplication + $f->fdiv(NSTR[,SCALE]) returns NSTR division to SCALE places + $f->fneg() return NSTR negation + $f->fabs() return NSTR absolute value + $f->fcmp(NSTR) return CODE compare undef,<0,=0,>0 + $f->fround(SCALE) return NSTR round to SCALE digits + $f->ffround(SCALE) return NSTR round at SCALEth place + $f->fnorm() return (NSTR) normalize + $f->fsqrt([SCALE]) return NSTR sqrt to SCALE places + +=head1 DESCRIPTION + +All basic math operations are overloaded if you declare your big +floats as + + $float = new Math::BigFloat "2.123123123123123123123123123123123"; + +=over 2 + +=item number format + +canonical strings have the form /[+-]\d+E[+-]\d+/ . Input values can +have inbedded whitespace. + +=item Error returns 'NaN' + +An input parameter was "Not a Number" or divide by zero or sqrt of +negative number. + +=item Division is computed to + +C<max($div_scale,length(dividend)+length(divisor))> digits by default. +Also used for default sqrt scale. + +=back + +=head1 BUGS + +The current version of this module is a preliminary version of the +real thing that is currently (as of perl5.002) under development. + +=head1 AUTHOR + +Mark Biggar + +=cut diff --git a/contrib/perl5/lib/Math/BigInt.pm b/contrib/perl5/lib/Math/BigInt.pm new file mode 100644 index 0000000..ef4af61 --- /dev/null +++ b/contrib/perl5/lib/Math/BigInt.pm @@ -0,0 +1,415 @@ +package Math::BigInt; + +use overload +'+' => sub {new Math::BigInt &badd}, +'-' => sub {new Math::BigInt + $_[2]? bsub($_[1],${$_[0]}) : bsub(${$_[0]},$_[1])}, +'<=>' => sub {new Math::BigInt + $_[2]? bcmp($_[1],${$_[0]}) : bcmp(${$_[0]},$_[1])}, +'cmp' => sub {new Math::BigInt + $_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])}, +'*' => sub {new Math::BigInt &bmul}, +'/' => sub {new Math::BigInt + $_[2]? scalar bdiv($_[1],${$_[0]}) : + scalar bdiv(${$_[0]},$_[1])}, +'%' => sub {new Math::BigInt + $_[2]? bmod($_[1],${$_[0]}) : bmod(${$_[0]},$_[1])}, +'**' => sub {new Math::BigInt + $_[2]? bpow($_[1],${$_[0]}) : bpow(${$_[0]},$_[1])}, +'neg' => sub {new Math::BigInt &bneg}, +'abs' => sub {new Math::BigInt &babs}, + +qw( +"" stringify +0+ numify) # Order of arguments unsignificant +; + +$NaNOK=1; + +sub new { + my($class) = shift; + my($foo) = bnorm(shift); + die "Not a number initialized to Math::BigInt" if !$NaNOK && $foo eq "NaN"; + bless \$foo, $class; +} +sub stringify { "${$_[0]}" } +sub numify { 0 + "${$_[0]}" } # Not needed, additional overhead + # comparing to direct compilation based on + # stringify +sub import { + shift; + return unless @_; + die "unknown import: @_" unless @_ == 1 and $_[0] eq ':constant'; + overload::constant integer => sub {Math::BigInt->new(shift)}; +} + +$zero = 0; + + +# normalize string form of number. Strip leading zeros. Strip any +# white space and add a sign, if missing. +# Strings that are not numbers result the value 'NaN'. + +sub bnorm { #(num_str) return num_str + local($_) = @_; + s/\s+//g; # strip white space + if (s/^([+-]?)0*(\d+)$/$1$2/) { # test if number + substr($_,$[,0) = '+' unless $1; # Add missing sign + s/^-0/+0/; + $_; + } else { + 'NaN'; + } +} + +# Convert a number from string format to internal base 100000 format. +# Assumes normalized value as input. +sub internal { #(num_str) return int_num_array + local($d) = @_; + ($is,$il) = (substr($d,$[,1),length($d)-2); + substr($d,$[,1) = ''; + ($is, reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $d))); +} + +# Convert a number from internal base 100000 format to string format. +# This routine scribbles all over input array. +sub external { #(int_num_array) return num_str + $es = shift; + grep($_ > 9999 || ($_ = substr('0000'.$_,-5)), @_); # zero pad + &bnorm(join('', $es, reverse(@_))); # reverse concat and normalize +} + +# Negate input value. +sub bneg { #(num_str) return num_str + local($_) = &bnorm(@_); + return $_ if $_ eq '+0' or $_ eq 'NaN'; + vec($_,0,8) ^= ord('+') ^ ord('-'); + $_; +} + +# Returns the absolute value of the input. +sub babs { #(num_str) return num_str + &abs(&bnorm(@_)); +} + +sub abs { # post-normalized abs for internal use + local($_) = @_; + s/^-/+/; + $_; +} + +# Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort) +sub bcmp { #(num_str, num_str) return cond_code + local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1])); + if ($x eq 'NaN') { + undef; + } elsif ($y eq 'NaN') { + undef; + } else { + &cmp($x,$y) <=> 0; + } +} + +sub cmp { # post-normalized compare for internal use + local($cx, $cy) = @_; + + return 0 if ($cx eq $cy); + + local($sx, $sy) = (substr($cx, 0, 1), substr($cy, 0, 1)); + local($ld); + + if ($sx eq '+') { + return 1 if ($sy eq '-' || $cy eq '+0'); + $ld = length($cx) - length($cy); + return $ld if ($ld); + return $cx cmp $cy; + } else { # $sx eq '-' + return -1 if ($sy eq '+'); + $ld = length($cy) - length($cx); + return $ld if ($ld); + return $cy cmp $cx; + } +} + +sub badd { #(num_str, num_str) return num_str + local(*x, *y); ($x, $y) = (&bnorm($_[$[]),&bnorm($_[$[+1])); + if ($x eq 'NaN') { + 'NaN'; + } elsif ($y eq 'NaN') { + 'NaN'; + } else { + @x = &internal($x); # convert to internal form + @y = &internal($y); + local($sx, $sy) = (shift @x, shift @y); # get signs + if ($sx eq $sy) { + &external($sx, &add(*x, *y)); # if same sign add + } else { + ($x, $y) = (&abs($x),&abs($y)); # make abs + if (&cmp($y,$x) > 0) { + &external($sy, &sub(*y, *x)); + } else { + &external($sx, &sub(*x, *y)); + } + } + } +} + +sub bsub { #(num_str, num_str) return num_str + &badd($_[$[],&bneg($_[$[+1])); +} + +# GCD -- Euclids algorithm Knuth Vol 2 pg 296 +sub bgcd { #(num_str, num_str) return num_str + local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1])); + if ($x eq 'NaN' || $y eq 'NaN') { + 'NaN'; + } else { + ($x, $y) = ($y,&bmod($x,$y)) while $y ne '+0'; + $x; + } +} + +# routine to add two base 1e5 numbers +# stolen from Knuth Vol 2 Algorithm A pg 231 +# there are separate routines to add and sub as per Kunth pg 233 +sub add { #(int_num_array, int_num_array) return int_num_array + local(*x, *y) = @_; + $car = 0; + for $x (@x) { + last unless @y || $car; + $x -= 1e5 if $car = (($x += (@y ? shift(@y) : 0) + $car) >= 1e5) ? 1 : 0; + } + for $y (@y) { + last unless $car; + $y -= 1e5 if $car = (($y += $car) >= 1e5) ? 1 : 0; + } + (@x, @y, $car); +} + +# subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y +sub sub { #(int_num_array, int_num_array) return int_num_array + local(*sx, *sy) = @_; + $bar = 0; + for $sx (@sx) { + last unless @sy || $bar; + $sx += 1e5 if $bar = (($sx -= (@sy ? shift(@sy) : 0) + $bar) < 0); + } + @sx; +} + +# multiply two numbers -- stolen from Knuth Vol 2 pg 233 +sub bmul { #(num_str, num_str) return num_str + local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1])); + if ($x eq 'NaN') { + 'NaN'; + } elsif ($y eq 'NaN') { + 'NaN'; + } else { + @x = &internal($x); + @y = &internal($y); + &external(&mul(*x,*y)); + } +} + +# multiply two numbers in internal representation +# destroys the arguments, supposes that two arguments are different +sub mul { #(*int_num_array, *int_num_array) return int_num_array + local(*x, *y) = (shift, shift); + local($signr) = (shift @x ne shift @y) ? '-' : '+'; + @prod = (); + for $x (@x) { + ($car, $cty) = (0, $[); + for $y (@y) { + $prod = $x * $y + ($prod[$cty] || 0) + $car; + $prod[$cty++] = + $prod - ($car = int($prod * 1e-5)) * 1e5; + } + $prod[$cty] += $car if $car; + $x = shift @prod; + } + ($signr, @x, @prod); +} + +# modulus +sub bmod { #(num_str, num_str) return num_str + (&bdiv(@_))[$[+1]; +} + +sub bdiv { #(dividend: num_str, divisor: num_str) return num_str + local (*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1])); + return wantarray ? ('NaN','NaN') : 'NaN' + if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0'); + return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0); + @x = &internal($x); @y = &internal($y); + $srem = $y[$[]; + $sr = (shift @x ne shift @y) ? '-' : '+'; + $car = $bar = $prd = 0; + if (($dd = int(1e5/($y[$#y]+1))) != 1) { + for $x (@x) { + $x = $x * $dd + $car; + $x -= ($car = int($x * 1e-5)) * 1e5; + } + push(@x, $car); $car = 0; + for $y (@y) { + $y = $y * $dd + $car; + $y -= ($car = int($y * 1e-5)) * 1e5; + } + } + else { + push(@x, 0); + } + @q = (); ($v2,$v1) = @y[-2,-1]; + while ($#x > $#y) { + ($u2,$u1,$u0) = @x[-3..-1]; + $q = (($u0 == $v1) ? 99999 : int(($u0*1e5+$u1)/$v1)); + --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*1e5+$u2); + if ($q) { + ($car, $bar) = (0,0); + for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) { + $prd = $q * $y[$y] + $car; + $prd -= ($car = int($prd * 1e-5)) * 1e5; + $x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0)); + } + if ($x[$#x] < $car + $bar) { + $car = 0; --$q; + for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) { + $x[$x] -= 1e5 + if ($car = (($x[$x] += $y[$y] + $car) > 1e5)); + } + } + } + pop(@x); unshift(@q, $q); + } + if (wantarray) { + @d = (); + if ($dd != 1) { + $car = 0; + for $x (reverse @x) { + $prd = $car * 1e5 + $x; + $car = $prd - ($tmp = int($prd / $dd)) * $dd; + unshift(@d, $tmp); + } + } + else { + @d = @x; + } + (&external($sr, @q), &external($srem, @d, $zero)); + } else { + &external($sr, @q); + } +} + +# compute power of two numbers -- stolen from Knuth Vol 2 pg 233 +sub bpow { #(num_str, num_str) return num_str + local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1])); + if ($x eq 'NaN') { + 'NaN'; + } elsif ($y eq 'NaN') { + 'NaN'; + } elsif ($x eq '+1') { + '+1'; + } elsif ($x eq '-1') { + &bmod($x,2) ? '-1': '+1'; + } elsif ($y =~ /^-/) { + 'NaN'; + } elsif ($x eq '+0' && $y eq '+0') { + 'NaN'; + } else { + @x = &internal($x); + local(@pow2)=@x; + local(@pow)=&internal("+1"); + local($y1,$res,@tmp1,@tmp2)=(1); # need tmp to send to mul + while ($y ne '+0') { + ($y,$res)=&bdiv($y,2); + if ($res ne '+0') {@tmp=@pow2; @pow=&mul(*pow,*tmp);} + if ($y ne '+0') {@tmp=@pow2;@pow2=&mul(*pow2,*tmp);} + } + &external(@pow); + } +} + +1; +__END__ + +=head1 NAME + +Math::BigInt - Arbitrary size integer math package + +=head1 SYNOPSIS + + use Math::BigInt; + $i = Math::BigInt->new($string); + + $i->bneg return BINT negation + $i->babs return BINT absolute value + $i->bcmp(BINT) return CODE compare numbers (undef,<0,=0,>0) + $i->badd(BINT) return BINT addition + $i->bsub(BINT) return BINT subtraction + $i->bmul(BINT) return BINT multiplication + $i->bdiv(BINT) return (BINT,BINT) division (quo,rem) just quo if scalar + $i->bmod(BINT) return BINT modulus + $i->bgcd(BINT) return BINT greatest common divisor + $i->bnorm return BINT normalization + +=head1 DESCRIPTION + +All basic math operations are overloaded if you declare your big +integers as + + $i = new Math::BigInt '123 456 789 123 456 789'; + + +=over 2 + +=item Canonical notation + +Big integer value are strings of the form C</^[+-]\d+$/> with leading +zeros suppressed. + +=item Input + +Input values to these routines may be strings of the form +C</^\s*[+-]?[\d\s]+$/>. + +=item Output + +Output values always always in canonical form + +=back + +Actual math is done in an internal format consisting of an array +whose first element is the sign (/^[+-]$/) and whose remaining +elements are base 100000 digits with the least significant digit first. +The string 'NaN' is used to represent the result when input arguments +are not numbers, as well as the result of dividing by zero. + +=head1 EXAMPLES + + '+0' canonical zero value + ' -123 123 123' canonical value '-123123123' + '1 23 456 7890' canonical value '+1234567890' + + +=head1 Autocreating constants + +After C<use Math::BigInt ':constant'> all the integer decimal constants +in the given scope are converted to C<Math::BigInt>. This conversion +happens at compile time. + +In particular + + perl -MMath::BigInt=:constant -e 'print 2**100' + +print the integer value of C<2**100>. Note that without convertion of +constants the expression 2**100 will be calculatted as floating point number. + +=head1 BUGS + +The current version of this module is a preliminary version of the +real thing that is currently (as of perl5.002) under development. + +=head1 AUTHOR + +Mark Biggar, overloaded interface by Ilya Zakharevich. + +=cut diff --git a/contrib/perl5/lib/Math/Complex.pm b/contrib/perl5/lib/Math/Complex.pm new file mode 100644 index 0000000..e711c14 --- /dev/null +++ b/contrib/perl5/lib/Math/Complex.pm @@ -0,0 +1,1775 @@ +# +# Complex numbers and associated mathematical functions +# -- Raphael Manfredi Since Sep 1996 +# -- Jarkko Hietaniemi Since Mar 1997 +# -- Daniel S. Lewart Since Sep 1997 +# + +require Exporter; +package Math::Complex; + +use strict; + +use vars qw($VERSION @ISA @EXPORT %EXPORT_TAGS); + +my ( $i, $ip2, %logn ); + +$VERSION = sprintf("%s", q$Id: Complex.pm,v 1.25 1998/02/05 16:07:37 jhi Exp $ =~ /(\d+\.\d+)/); + +@ISA = qw(Exporter); + +my @trig = qw( + pi + tan + csc cosec sec cot cotan + asin acos atan + acsc acosec asec acot acotan + sinh cosh tanh + csch cosech sech coth cotanh + asinh acosh atanh + acsch acosech asech acoth acotanh + ); + +@EXPORT = (qw( + i Re Im rho theta arg + sqrt log ln + log10 logn cbrt root + cplx cplxe + ), + @trig); + +%EXPORT_TAGS = ( + 'trig' => [@trig], +); + +use overload + '+' => \&plus, + '-' => \&minus, + '*' => \&multiply, + '/' => \÷, + '**' => \&power, + '<=>' => \&spaceship, + 'neg' => \&negate, + '~' => \&conjugate, + 'abs' => \&abs, + 'sqrt' => \&sqrt, + 'exp' => \&exp, + 'log' => \&log, + 'sin' => \&sin, + 'cos' => \&cos, + 'tan' => \&tan, + 'atan2' => \&atan2, + qw("" stringify); + +# +# Package "privates" +# + +my $package = 'Math::Complex'; # Package name +my $display = 'cartesian'; # Default display format +my $eps = 1e-14; # Epsilon + +# +# Object attributes (internal): +# cartesian [real, imaginary] -- cartesian form +# polar [rho, theta] -- polar form +# c_dirty cartesian form not up-to-date +# p_dirty polar form not up-to-date +# display display format (package's global when not set) +# + +# Die on bad *make() arguments. + +sub _cannot_make { + die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n"; +} + +# +# ->make +# +# Create a new complex number (cartesian form) +# +sub make { + my $self = bless {}, shift; + my ($re, $im) = @_; + my $rre = ref $re; + if ( $rre ) { + if ( $rre eq ref $self ) { + $re = Re($re); + } else { + _cannot_make("real part", $rre); + } + } + my $rim = ref $im; + if ( $rim ) { + if ( $rim eq ref $self ) { + $im = Im($im); + } else { + _cannot_make("imaginary part", $rim); + } + } + $self->{'cartesian'} = [ $re, $im ]; + $self->{c_dirty} = 0; + $self->{p_dirty} = 1; + $self->display_format('cartesian'); + return $self; +} + +# +# ->emake +# +# Create a new complex number (exponential form) +# +sub emake { + my $self = bless {}, shift; + my ($rho, $theta) = @_; + my $rrh = ref $rho; + if ( $rrh ) { + if ( $rrh eq ref $self ) { + $rho = rho($rho); + } else { + _cannot_make("rho", $rrh); + } + } + my $rth = ref $theta; + if ( $rth ) { + if ( $rth eq ref $self ) { + $theta = theta($theta); + } else { + _cannot_make("theta", $rth); + } + } + if ($rho < 0) { + $rho = -$rho; + $theta = ($theta <= 0) ? $theta + pi() : $theta - pi(); + } + $self->{'polar'} = [$rho, $theta]; + $self->{p_dirty} = 0; + $self->{c_dirty} = 1; + $self->display_format('polar'); + return $self; +} + +sub new { &make } # For backward compatibility only. + +# +# cplx +# +# Creates a complex number from a (re, im) tuple. +# This avoids the burden of writing Math::Complex->make(re, im). +# +sub cplx { + my ($re, $im) = @_; + return $package->make($re, defined $im ? $im : 0); +} + +# +# cplxe +# +# Creates a complex number from a (rho, theta) tuple. +# This avoids the burden of writing Math::Complex->emake(rho, theta). +# +sub cplxe { + my ($rho, $theta) = @_; + return $package->emake($rho, defined $theta ? $theta : 0); +} + +# +# pi +# +# The number defined as pi = 180 degrees +# +use constant pi => 4 * CORE::atan2(1, 1); + +# +# pit2 +# +# The full circle +# +use constant pit2 => 2 * pi; + +# +# pip2 +# +# The quarter circle +# +use constant pip2 => pi / 2; + +# +# deg1 +# +# One degree in radians, used in stringify_polar. +# + +use constant deg1 => pi / 180; + +# +# uplog10 +# +# Used in log10(). +# +use constant uplog10 => 1 / CORE::log(10); + +# +# i +# +# The number defined as i*i = -1; +# +sub i () { + return $i if ($i); + $i = bless {}; + $i->{'cartesian'} = [0, 1]; + $i->{'polar'} = [1, pip2]; + $i->{c_dirty} = 0; + $i->{p_dirty} = 0; + return $i; +} + +# +# Attribute access/set routines +# + +sub cartesian {$_[0]->{c_dirty} ? + $_[0]->update_cartesian : $_[0]->{'cartesian'}} +sub polar {$_[0]->{p_dirty} ? + $_[0]->update_polar : $_[0]->{'polar'}} + +sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] } +sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] } + +# +# ->update_cartesian +# +# Recompute and return the cartesian form, given accurate polar form. +# +sub update_cartesian { + my $self = shift; + my ($r, $t) = @{$self->{'polar'}}; + $self->{c_dirty} = 0; + return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)]; +} + +# +# +# ->update_polar +# +# Recompute and return the polar form, given accurate cartesian form. +# +sub update_polar { + my $self = shift; + my ($x, $y) = @{$self->{'cartesian'}}; + $self->{p_dirty} = 0; + return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0; + return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y), CORE::atan2($y, $x)]; +} + +# +# (plus) +# +# Computes z1+z2. +# +sub plus { + my ($z1, $z2, $regular) = @_; + my ($re1, $im1) = @{$z1->cartesian}; + $z2 = cplx($z2) unless ref $z2; + my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); + unless (defined $regular) { + $z1->set_cartesian([$re1 + $re2, $im1 + $im2]); + return $z1; + } + return (ref $z1)->make($re1 + $re2, $im1 + $im2); +} + +# +# (minus) +# +# Computes z1-z2. +# +sub minus { + my ($z1, $z2, $inverted) = @_; + my ($re1, $im1) = @{$z1->cartesian}; + $z2 = cplx($z2) unless ref $z2; + my ($re2, $im2) = @{$z2->cartesian}; + unless (defined $inverted) { + $z1->set_cartesian([$re1 - $re2, $im1 - $im2]); + return $z1; + } + return $inverted ? + (ref $z1)->make($re2 - $re1, $im2 - $im1) : + (ref $z1)->make($re1 - $re2, $im1 - $im2); + +} + +# +# (multiply) +# +# Computes z1*z2. +# +sub multiply { + my ($z1, $z2, $regular) = @_; + if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) { + # if both polar better use polar to avoid rounding errors + my ($r1, $t1) = @{$z1->polar}; + my ($r2, $t2) = @{$z2->polar}; + my $t = $t1 + $t2; + if ($t > pi()) { $t -= pit2 } + elsif ($t <= -pi()) { $t += pit2 } + unless (defined $regular) { + $z1->set_polar([$r1 * $r2, $t]); + return $z1; + } + return (ref $z1)->emake($r1 * $r2, $t); + } else { + my ($x1, $y1) = @{$z1->cartesian}; + if (ref $z2) { + my ($x2, $y2) = @{$z2->cartesian}; + return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2); + } else { + return (ref $z1)->make($x1*$z2, $y1*$z2); + } + } +} + +# +# _divbyzero +# +# Die on division by zero. +# +sub _divbyzero { + my $mess = "$_[0]: Division by zero.\n"; + + if (defined $_[1]) { + $mess .= "(Because in the definition of $_[0], the divisor "; + $mess .= "$_[1] " unless ($_[1] eq '0'); + $mess .= "is 0)\n"; + } + + my @up = caller(1); + + $mess .= "Died at $up[1] line $up[2].\n"; + + die $mess; +} + +# +# (divide) +# +# Computes z1/z2. +# +sub divide { + my ($z1, $z2, $inverted) = @_; + if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) { + # if both polar better use polar to avoid rounding errors + my ($r1, $t1) = @{$z1->polar}; + my ($r2, $t2) = @{$z2->polar}; + my $t; + if ($inverted) { + _divbyzero "$z2/0" if ($r1 == 0); + $t = $t2 - $t1; + if ($t > pi()) { $t -= pit2 } + elsif ($t <= -pi()) { $t += pit2 } + return (ref $z1)->emake($r2 / $r1, $t); + } else { + _divbyzero "$z1/0" if ($r2 == 0); + $t = $t1 - $t2; + if ($t > pi()) { $t -= pit2 } + elsif ($t <= -pi()) { $t += pit2 } + return (ref $z1)->emake($r1 / $r2, $t); + } + } else { + my ($d, $x2, $y2); + if ($inverted) { + ($x2, $y2) = @{$z1->cartesian}; + $d = $x2*$x2 + $y2*$y2; + _divbyzero "$z2/0" if $d == 0; + return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d); + } else { + my ($x1, $y1) = @{$z1->cartesian}; + if (ref $z2) { + ($x2, $y2) = @{$z2->cartesian}; + $d = $x2*$x2 + $y2*$y2; + _divbyzero "$z1/0" if $d == 0; + my $u = ($x1*$x2 + $y1*$y2)/$d; + my $v = ($y1*$x2 - $x1*$y2)/$d; + return (ref $z1)->make($u, $v); + } else { + _divbyzero "$z1/0" if $z2 == 0; + return (ref $z1)->make($x1/$z2, $y1/$z2); + } + } + } +} + +# +# _zerotozero +# +# Die on zero raised to the zeroth. +# +sub _zerotozero { + my $mess = "The zero raised to the zeroth power is not defined.\n"; + + my @up = caller(1); + + $mess .= "Died at $up[1] line $up[2].\n"; + + die $mess; +} + +# +# (power) +# +# Computes z1**z2 = exp(z2 * log z1)). +# +sub power { + my ($z1, $z2, $inverted) = @_; + my $z1z = $z1 == 0; + my $z2z = $z2 == 0; + _zerotozero if ($z1z and $z2z); + if ($inverted) { + return 0 if ($z2z); + return 1 if ($z1z or $z2 == 1); + } else { + return 0 if ($z1z); + return 1 if ($z2z or $z1 == 1); + } + my $w = $inverted ? CORE::exp($z1 * CORE::log($z2)) : CORE::exp($z2 * CORE::log($z1)); + # If both arguments cartesian, return cartesian, else polar. + return $z1->{c_dirty} == 0 && + (not ref $z2 or $z2->{c_dirty} == 0) ? + cplx(@{$w->cartesian}) : $w; +} + +# +# (spaceship) +# +# Computes z1 <=> z2. +# Sorts on the real part first, then on the imaginary part. Thus 2-4i > 3+8i. +# +sub spaceship { + my ($z1, $z2, $inverted) = @_; + my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0); + my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); + my $sgn = $inverted ? -1 : 1; + return $sgn * ($re1 <=> $re2) if $re1 != $re2; + return $sgn * ($im1 <=> $im2); +} + +# +# (negate) +# +# Computes -z. +# +sub negate { + my ($z) = @_; + if ($z->{c_dirty}) { + my ($r, $t) = @{$z->polar}; + $t = ($t <= 0) ? $t + pi : $t - pi; + return (ref $z)->emake($r, $t); + } + my ($re, $im) = @{$z->cartesian}; + return (ref $z)->make(-$re, -$im); +} + +# +# (conjugate) +# +# Compute complex's conjugate. +# +sub conjugate { + my ($z) = @_; + if ($z->{c_dirty}) { + my ($r, $t) = @{$z->polar}; + return (ref $z)->emake($r, -$t); + } + my ($re, $im) = @{$z->cartesian}; + return (ref $z)->make($re, -$im); +} + +# +# (abs) +# +# Compute or set complex's norm (rho). +# +sub abs { + my ($z, $rho) = @_; + return $z unless ref $z; + if (defined $rho) { + $z->{'polar'} = [ $rho, ${$z->polar}[1] ]; + $z->{p_dirty} = 0; + $z->{c_dirty} = 1; + return $rho; + } else { + return ${$z->polar}[0]; + } +} + +sub _theta { + my $theta = $_[0]; + + if ($$theta > pi()) { $$theta -= pit2 } + elsif ($$theta <= -pi()) { $$theta += pit2 } +} + +# +# arg +# +# Compute or set complex's argument (theta). +# +sub arg { + my ($z, $theta) = @_; + return $z unless ref $z; + if (defined $theta) { + _theta(\$theta); + $z->{'polar'} = [ ${$z->polar}[0], $theta ]; + $z->{p_dirty} = 0; + $z->{c_dirty} = 1; + } else { + $theta = ${$z->polar}[1]; + _theta(\$theta); + } + return $theta; +} + +# +# (sqrt) +# +# Compute sqrt(z). +# +# It is quite tempting to use wantarray here so that in list context +# sqrt() would return the two solutions. This, however, would +# break things like +# +# print "sqrt(z) = ", sqrt($z), "\n"; +# +# The two values would be printed side by side without no intervening +# whitespace, quite confusing. +# Therefore if you want the two solutions use the root(). +# +sub sqrt { + my ($z) = @_; + my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0); + return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re) if $im == 0; + my ($r, $t) = @{$z->polar}; + return (ref $z)->emake(CORE::sqrt($r), $t/2); +} + +# +# cbrt +# +# Compute cbrt(z) (cubic root). +# +# Why are we not returning three values? The same answer as for sqrt(). +# +sub cbrt { + my ($z) = @_; + return $z < 0 ? -CORE::exp(CORE::log(-$z)/3) : ($z > 0 ? CORE::exp(CORE::log($z)/3): 0) + unless ref $z; + my ($r, $t) = @{$z->polar}; + return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3); +} + +# +# _rootbad +# +# Die on bad root. +# +sub _rootbad { + my $mess = "Root $_[0] not defined, root must be positive integer.\n"; + + my @up = caller(1); + + $mess .= "Died at $up[1] line $up[2].\n"; + + die $mess; +} + +# +# root +# +# Computes all nth root for z, returning an array whose size is n. +# `n' must be a positive integer. +# +# The roots are given by (for k = 0..n-1): +# +# z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n)) +# +sub root { + my ($z, $n) = @_; + _rootbad($n) if ($n < 1 or int($n) != $n); + my ($r, $t) = ref $z ? @{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi); + my @root; + my $k; + my $theta_inc = pit2 / $n; + my $rho = $r ** (1/$n); + my $theta; + my $cartesian = ref $z && $z->{c_dirty} == 0; + for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) { + my $w = cplxe($rho, $theta); + # Yes, $cartesian is loop invariant. + push @root, $cartesian ? cplx(@{$w->cartesian}) : $w; + } + return @root; +} + +# +# Re +# +# Return or set Re(z). +# +sub Re { + my ($z, $Re) = @_; + return $z unless ref $z; + if (defined $Re) { + $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ]; + $z->{c_dirty} = 0; + $z->{p_dirty} = 1; + } else { + return ${$z->cartesian}[0]; + } +} + +# +# Im +# +# Return or set Im(z). +# +sub Im { + my ($z, $Im) = @_; + return $z unless ref $z; + if (defined $Im) { + $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ]; + $z->{c_dirty} = 0; + $z->{p_dirty} = 1; + } else { + return ${$z->cartesian}[1]; + } +} + +# +# rho +# +# Return or set rho(w). +# +sub rho { + Math::Complex::abs(@_); +} + +# +# theta +# +# Return or set theta(w). +# +sub theta { + Math::Complex::arg(@_); +} + +# +# (exp) +# +# Computes exp(z). +# +sub exp { + my ($z) = @_; + my ($x, $y) = @{$z->cartesian}; + return (ref $z)->emake(CORE::exp($x), $y); +} + +# +# _logofzero +# +# Die on logarithm of zero. +# +sub _logofzero { + my $mess = "$_[0]: Logarithm of zero.\n"; + + if (defined $_[1]) { + $mess .= "(Because in the definition of $_[0], the argument "; + $mess .= "$_[1] " unless ($_[1] eq '0'); + $mess .= "is 0)\n"; + } + + my @up = caller(1); + + $mess .= "Died at $up[1] line $up[2].\n"; + + die $mess; +} + +# +# (log) +# +# Compute log(z). +# +sub log { + my ($z) = @_; + unless (ref $z) { + _logofzero("log") if $z == 0; + return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi); + } + my ($r, $t) = @{$z->polar}; + _logofzero("log") if $r == 0; + if ($t > pi()) { $t -= pit2 } + elsif ($t <= -pi()) { $t += pit2 } + return (ref $z)->make(CORE::log($r), $t); +} + +# +# ln +# +# Alias for log(). +# +sub ln { Math::Complex::log(@_) } + +# +# log10 +# +# Compute log10(z). +# + +sub log10 { + return Math::Complex::log($_[0]) * uplog10; +} + +# +# logn +# +# Compute logn(z,n) = log(z) / log(n) +# +sub logn { + my ($z, $n) = @_; + $z = cplx($z, 0) unless ref $z; + my $logn = $logn{$n}; + $logn = $logn{$n} = CORE::log($n) unless defined $logn; # Cache log(n) + return CORE::log($z) / $logn; +} + +# +# (cos) +# +# Compute cos(z) = (exp(iz) + exp(-iz))/2. +# +sub cos { + my ($z) = @_; + my ($x, $y) = @{$z->cartesian}; + my $ey = CORE::exp($y); + my $ey_1 = 1 / $ey; + return (ref $z)->make(CORE::cos($x) * ($ey + $ey_1)/2, + CORE::sin($x) * ($ey_1 - $ey)/2); +} + +# +# (sin) +# +# Compute sin(z) = (exp(iz) - exp(-iz))/2. +# +sub sin { + my ($z) = @_; + my ($x, $y) = @{$z->cartesian}; + my $ey = CORE::exp($y); + my $ey_1 = 1 / $ey; + return (ref $z)->make(CORE::sin($x) * ($ey + $ey_1)/2, + CORE::cos($x) * ($ey - $ey_1)/2); +} + +# +# tan +# +# Compute tan(z) = sin(z) / cos(z). +# +sub tan { + my ($z) = @_; + my $cz = CORE::cos($z); + _divbyzero "tan($z)", "cos($z)" if (CORE::abs($cz) < $eps); + return CORE::sin($z) / $cz; +} + +# +# sec +# +# Computes the secant sec(z) = 1 / cos(z). +# +sub sec { + my ($z) = @_; + my $cz = CORE::cos($z); + _divbyzero "sec($z)", "cos($z)" if ($cz == 0); + return 1 / $cz; +} + +# +# csc +# +# Computes the cosecant csc(z) = 1 / sin(z). +# +sub csc { + my ($z) = @_; + my $sz = CORE::sin($z); + _divbyzero "csc($z)", "sin($z)" if ($sz == 0); + return 1 / $sz; +} + +# +# cosec +# +# Alias for csc(). +# +sub cosec { Math::Complex::csc(@_) } + +# +# cot +# +# Computes cot(z) = cos(z) / sin(z). +# +sub cot { + my ($z) = @_; + my $sz = CORE::sin($z); + _divbyzero "cot($z)", "sin($z)" if ($sz == 0); + return CORE::cos($z) / $sz; +} + +# +# cotan +# +# Alias for cot(). +# +sub cotan { Math::Complex::cot(@_) } + +# +# acos +# +# Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)). +# +sub acos { + my $z = $_[0]; + return CORE::atan2(CORE::sqrt(1-$z*$z), $z) if (! ref $z) && CORE::abs($z) <= 1; + my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0); + my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y); + my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y); + my $alpha = ($t1 + $t2)/2; + my $beta = ($t1 - $t2)/2; + $alpha = 1 if $alpha < 1; + if ($beta > 1) { $beta = 1 } + elsif ($beta < -1) { $beta = -1 } + my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta); + my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1)); + $v = -$v if $y > 0 || ($y == 0 && $x < -1); + return $package->make($u, $v); +} + +# +# asin +# +# Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)). +# +sub asin { + my $z = $_[0]; + return CORE::atan2($z, CORE::sqrt(1-$z*$z)) if (! ref $z) && CORE::abs($z) <= 1; + my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0); + my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y); + my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y); + my $alpha = ($t1 + $t2)/2; + my $beta = ($t1 - $t2)/2; + $alpha = 1 if $alpha < 1; + if ($beta > 1) { $beta = 1 } + elsif ($beta < -1) { $beta = -1 } + my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta)); + my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1)); + $v = -$v if $y > 0 || ($y == 0 && $x < -1); + return $package->make($u, $v); +} + +# +# atan +# +# Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)). +# +sub atan { + my ($z) = @_; + return CORE::atan2($z, 1) unless ref $z; + _divbyzero "atan(i)" if ( $z == i); + _divbyzero "atan(-i)" if (-$z == i); + my $log = CORE::log((i + $z) / (i - $z)); + $ip2 = 0.5 * i unless defined $ip2; + return $ip2 * $log; +} + +# +# asec +# +# Computes the arc secant asec(z) = acos(1 / z). +# +sub asec { + my ($z) = @_; + _divbyzero "asec($z)", $z if ($z == 0); + return acos(1 / $z); +} + +# +# acsc +# +# Computes the arc cosecant acsc(z) = asin(1 / z). +# +sub acsc { + my ($z) = @_; + _divbyzero "acsc($z)", $z if ($z == 0); + return asin(1 / $z); +} + +# +# acosec +# +# Alias for acsc(). +# +sub acosec { Math::Complex::acsc(@_) } + +# +# acot +# +# Computes the arc cotangent acot(z) = atan(1 / z) +# +sub acot { + my ($z) = @_; + _divbyzero "acot(0)" if (CORE::abs($z) < $eps); + return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z) unless ref $z; + _divbyzero "acot(i)" if (CORE::abs($z - i) < $eps); + _logofzero "acot(-i)" if (CORE::abs($z + i) < $eps); + return atan(1 / $z); +} + +# +# acotan +# +# Alias for acot(). +# +sub acotan { Math::Complex::acot(@_) } + +# +# cosh +# +# Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2. +# +sub cosh { + my ($z) = @_; + my $ex; + unless (ref $z) { + $ex = CORE::exp($z); + return ($ex + 1/$ex)/2; + } + my ($x, $y) = @{$z->cartesian}; + $ex = CORE::exp($x); + my $ex_1 = 1 / $ex; + return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2, + CORE::sin($y) * ($ex - $ex_1)/2); +} + +# +# sinh +# +# Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2. +# +sub sinh { + my ($z) = @_; + my $ex; + unless (ref $z) { + $ex = CORE::exp($z); + return ($ex - 1/$ex)/2; + } + my ($x, $y) = @{$z->cartesian}; + $ex = CORE::exp($x); + my $ex_1 = 1 / $ex; + return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2, + CORE::sin($y) * ($ex + $ex_1)/2); +} + +# +# tanh +# +# Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z). +# +sub tanh { + my ($z) = @_; + my $cz = cosh($z); + _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0); + return sinh($z) / $cz; +} + +# +# sech +# +# Computes the hyperbolic secant sech(z) = 1 / cosh(z). +# +sub sech { + my ($z) = @_; + my $cz = cosh($z); + _divbyzero "sech($z)", "cosh($z)" if ($cz == 0); + return 1 / $cz; +} + +# +# csch +# +# Computes the hyperbolic cosecant csch(z) = 1 / sinh(z). +# +sub csch { + my ($z) = @_; + my $sz = sinh($z); + _divbyzero "csch($z)", "sinh($z)" if ($sz == 0); + return 1 / $sz; +} + +# +# cosech +# +# Alias for csch(). +# +sub cosech { Math::Complex::csch(@_) } + +# +# coth +# +# Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z). +# +sub coth { + my ($z) = @_; + my $sz = sinh($z); + _divbyzero "coth($z)", "sinh($z)" if ($sz == 0); + return cosh($z) / $sz; +} + +# +# cotanh +# +# Alias for coth(). +# +sub cotanh { Math::Complex::coth(@_) } + +# +# acosh +# +# Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)). +# +sub acosh { + my ($z) = @_; + unless (ref $z) { + return CORE::log($z + CORE::sqrt($z*$z-1)) if $z >= 1; + $z = cplx($z, 0); + } + my ($re, $im) = @{$z->cartesian}; + if ($im == 0) { + return cplx(CORE::log($re + CORE::sqrt($re*$re - 1)), 0) if $re >= 1; + return cplx(0, CORE::atan2(CORE::sqrt(1-$re*$re), $re)) if CORE::abs($re) <= 1; + } + return CORE::log($z + CORE::sqrt($z*$z - 1)); +} + +# +# asinh +# +# Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1)) +# +sub asinh { + my ($z) = @_; + return CORE::log($z + CORE::sqrt($z*$z + 1)); +} + +# +# atanh +# +# Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)). +# +sub atanh { + my ($z) = @_; + unless (ref $z) { + return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1; + $z = cplx($z, 0); + } + _divbyzero 'atanh(1)', "1 - $z" if ($z == 1); + _logofzero 'atanh(-1)' if ($z == -1); + return 0.5 * CORE::log((1 + $z) / (1 - $z)); +} + +# +# asech +# +# Computes the hyperbolic arc secant asech(z) = acosh(1 / z). +# +sub asech { + my ($z) = @_; + _divbyzero 'asech(0)', $z if ($z == 0); + return acosh(1 / $z); +} + +# +# acsch +# +# Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z). +# +sub acsch { + my ($z) = @_; + _divbyzero 'acsch(0)', $z if ($z == 0); + return asinh(1 / $z); +} + +# +# acosech +# +# Alias for acosh(). +# +sub acosech { Math::Complex::acsch(@_) } + +# +# acoth +# +# Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)). +# +sub acoth { + my ($z) = @_; + _divbyzero 'acoth(0)' if (CORE::abs($z) < $eps); + unless (ref $z) { + return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1; + $z = cplx($z, 0); + } + _divbyzero 'acoth(1)', "$z - 1" if (CORE::abs($z - 1) < $eps); + _logofzero 'acoth(-1)', "1 / $z" if (CORE::abs($z + 1) < $eps); + return CORE::log((1 + $z) / ($z - 1)) / 2; +} + +# +# acotanh +# +# Alias for acot(). +# +sub acotanh { Math::Complex::acoth(@_) } + +# +# (atan2) +# +# Compute atan(z1/z2). +# +sub atan2 { + my ($z1, $z2, $inverted) = @_; + my ($re1, $im1, $re2, $im2); + if ($inverted) { + ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); + ($re2, $im2) = @{$z1->cartesian}; + } else { + ($re1, $im1) = @{$z1->cartesian}; + ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); + } + if ($im2 == 0) { + return cplx(CORE::atan2($re1, $re2), 0) if $im1 == 0; + return cplx(($im1<=>0) * pip2, 0) if $re2 == 0; + } + my $w = atan($z1/$z2); + my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0); + $u += pi if $re2 < 0; + $u -= pit2 if $u > pi; + return cplx($u, $v); +} + +# +# display_format +# ->display_format +# +# Set (fetch if no argument) display format for all complex numbers that +# don't happen to have overridden it via ->display_format +# +# When called as a method, this actually sets the display format for +# the current object. +# +# Valid object formats are 'c' and 'p' for cartesian and polar. The first +# letter is used actually, so the type can be fully spelled out for clarity. +# +sub display_format { + my $self = shift; + my $format = undef; + + if (ref $self) { # Called as a method + $format = shift; + } else { # Regular procedure call + $format = $self; + undef $self; + } + + if (defined $self) { + return defined $self->{display} ? $self->{display} : $display + unless defined $format; + return $self->{display} = $format; + } + + return $display unless defined $format; + return $display = $format; +} + +# +# (stringify) +# +# Show nicely formatted complex number under its cartesian or polar form, +# depending on the current display format: +# +# . If a specific display format has been recorded for this object, use it. +# . Otherwise, use the generic current default for all complex numbers, +# which is a package global variable. +# +sub stringify { + my ($z) = shift; + my $format; + + $format = $display; + $format = $z->{display} if defined $z->{display}; + + return $z->stringify_polar if $format =~ /^p/i; + return $z->stringify_cartesian; +} + +# +# ->stringify_cartesian +# +# Stringify as a cartesian representation 'a+bi'. +# +sub stringify_cartesian { + my $z = shift; + my ($x, $y) = @{$z->cartesian}; + my ($re, $im); + + $x = int($x + ($x < 0 ? -1 : 1) * $eps) + if int(CORE::abs($x)) != int(CORE::abs($x) + $eps); + $y = int($y + ($y < 0 ? -1 : 1) * $eps) + if int(CORE::abs($y)) != int(CORE::abs($y) + $eps); + + $re = "$x" if CORE::abs($x) >= $eps; + if ($y == 1) { $im = 'i' } + elsif ($y == -1) { $im = '-i' } + elsif (CORE::abs($y) >= $eps) { $im = $y . "i" } + + my $str = ''; + $str = $re if defined $re; + $str .= "+$im" if defined $im; + $str =~ s/\+-/-/; + $str =~ s/^\+//; + $str =~ s/([-+])1i/$1i/; # Not redundant with the above 1/-1 tests. + $str = '0' unless $str; + + return $str; +} + + +# Helper for stringify_polar, a Greatest Common Divisor with a memory. + +sub _gcd { + my ($a, $b) = @_; + + use integer; + + # Loops forever if given negative inputs. + + if ($b and $a > $b) { return gcd($a % $b, $b) } + elsif ($a and $b > $a) { return gcd($b % $a, $a) } + else { return $a ? $a : $b } +} + +my %gcd; + +sub gcd { + my ($a, $b) = @_; + + my $id = "$a $b"; + + unless (exists $gcd{$id}) { + $gcd{$id} = _gcd($a, $b); + $gcd{"$b $a"} = $gcd{$id}; + } + + return $gcd{$id}; +} + +# +# ->stringify_polar +# +# Stringify as a polar representation '[r,t]'. +# +sub stringify_polar { + my $z = shift; + my ($r, $t) = @{$z->polar}; + my $theta; + + return '[0,0]' if $r <= $eps; + + my $nt = $t / pit2; + $nt = ($nt - int($nt)) * pit2; + $nt += pit2 if $nt < 0; # Range [0, 2pi] + + if (CORE::abs($nt) <= $eps) { $theta = 0 } + elsif (CORE::abs(pi-$nt) <= $eps) { $theta = 'pi' } + + if (defined $theta) { + $r = int($r + ($r < 0 ? -1 : 1) * $eps) + if int(CORE::abs($r)) != int(CORE::abs($r) + $eps); + $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps) + if ($theta ne 'pi' and + int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps)); + return "\[$r,$theta\]"; + } + + # + # Okay, number is not a real. Try to identify pi/n and friends... + # + + $nt -= pit2 if $nt > pi; + + if (CORE::abs($nt) >= deg1) { + my ($n, $k, $kpi); + + for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) { + $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5); + if (CORE::abs($kpi/$n - $nt) <= $eps) { + $n = CORE::abs($n); + my $gcd = gcd($k, $n); + if ($gcd > 1) { + $k /= $gcd; + $n /= $gcd; + } + next if $n > 360; + $theta = ($nt < 0 ? '-':''). + ($k == 1 ? 'pi':"${k}pi"); + $theta .= '/'.$n if $n > 1; + last; + } + } + } + + $theta = $nt unless defined $theta; + + $r = int($r + ($r < 0 ? -1 : 1) * $eps) + if int(CORE::abs($r)) != int(CORE::abs($r) + $eps); + $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps) + if ($theta !~ m(^-?\d*pi/\d+$) and + int(CORE::abs($theta)) != int(CORE::abs($theta) + $eps)); + + return "\[$r,$theta\]"; +} + +1; +__END__ + +=head1 NAME + +Math::Complex - complex numbers and associated mathematical functions + +=head1 SYNOPSIS + + use Math::Complex; + + $z = Math::Complex->make(5, 6); + $t = 4 - 3*i + $z; + $j = cplxe(1, 2*pi/3); + +=head1 DESCRIPTION + +This package lets you create and manipulate complex numbers. By default, +I<Perl> limits itself to real numbers, but an extra C<use> statement brings +full complex support, along with a full set of mathematical functions +typically associated with and/or extended to complex numbers. + +If you wonder what complex numbers are, they were invented to be able to solve +the following equation: + + x*x = -1 + +and by definition, the solution is noted I<i> (engineers use I<j> instead since +I<i> usually denotes an intensity, but the name does not matter). The number +I<i> is a pure I<imaginary> number. + +The arithmetics with pure imaginary numbers works just like you would expect +it with real numbers... you just have to remember that + + i*i = -1 + +so you have: + + 5i + 7i = i * (5 + 7) = 12i + 4i - 3i = i * (4 - 3) = i + 4i * 2i = -8 + 6i / 2i = 3 + 1 / i = -i + +Complex numbers are numbers that have both a real part and an imaginary +part, and are usually noted: + + a + bi + +where C<a> is the I<real> part and C<b> is the I<imaginary> part. The +arithmetic with complex numbers is straightforward. You have to +keep track of the real and the imaginary parts, but otherwise the +rules used for real numbers just apply: + + (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i + (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i + +A graphical representation of complex numbers is possible in a plane +(also called the I<complex plane>, but it's really a 2D plane). +The number + + z = a + bi + +is the point whose coordinates are (a, b). Actually, it would +be the vector originating from (0, 0) to (a, b). It follows that the addition +of two complex numbers is a vectorial addition. + +Since there is a bijection between a point in the 2D plane and a complex +number (i.e. the mapping is unique and reciprocal), a complex number +can also be uniquely identified with polar coordinates: + + [rho, theta] + +where C<rho> is the distance to the origin, and C<theta> the angle between +the vector and the I<x> axis. There is a notation for this using the +exponential form, which is: + + rho * exp(i * theta) + +where I<i> is the famous imaginary number introduced above. Conversion +between this form and the cartesian form C<a + bi> is immediate: + + a = rho * cos(theta) + b = rho * sin(theta) + +which is also expressed by this formula: + + z = rho * exp(i * theta) = rho * (cos theta + i * sin theta) + +In other words, it's the projection of the vector onto the I<x> and I<y> +axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta> +the I<argument> of the complex number. The I<norm> of C<z> will be +noted C<abs(z)>. + +The polar notation (also known as the trigonometric +representation) is much more handy for performing multiplications and +divisions of complex numbers, whilst the cartesian notation is better +suited for additions and subtractions. Real numbers are on the I<x> +axis, and therefore I<theta> is zero or I<pi>. + +All the common operations that can be performed on a real number have +been defined to work on complex numbers as well, and are merely +I<extensions> of the operations defined on real numbers. This means +they keep their natural meaning when there is no imaginary part, provided +the number is within their definition set. + +For instance, the C<sqrt> routine which computes the square root of +its argument is only defined for non-negative real numbers and yields a +non-negative real number (it is an application from B<R+> to B<R+>). +If we allow it to return a complex number, then it can be extended to +negative real numbers to become an application from B<R> to B<C> (the +set of complex numbers): + + sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i + +It can also be extended to be an application from B<C> to B<C>, +whilst its restriction to B<R> behaves as defined above by using +the following definition: + + sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2) + +Indeed, a negative real number can be noted C<[x,pi]> (the modulus +I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative +number) and the above definition states that + + sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i + +which is exactly what we had defined for negative real numbers above. +The C<sqrt> returns only one of the solutions: if you want the both, +use the C<root> function. + +All the common mathematical functions defined on real numbers that +are extended to complex numbers share that same property of working +I<as usual> when the imaginary part is zero (otherwise, it would not +be called an extension, would it?). + +A I<new> operation possible on a complex number that is +the identity for real numbers is called the I<conjugate>, and is noted +with an horizontal bar above the number, or C<~z> here. + + z = a + bi + ~z = a - bi + +Simple... Now look: + + z * ~z = (a + bi) * (a - bi) = a*a + b*b + +We saw that the norm of C<z> was noted C<abs(z)> and was defined as the +distance to the origin, also known as: + + rho = abs(z) = sqrt(a*a + b*b) + +so + + z * ~z = abs(z) ** 2 + +If z is a pure real number (i.e. C<b == 0>), then the above yields: + + a * a = abs(a) ** 2 + +which is true (C<abs> has the regular meaning for real number, i.e. stands +for the absolute value). This example explains why the norm of C<z> is +noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet +is the regular C<abs> we know when the complex number actually has no +imaginary part... This justifies I<a posteriori> our use of the C<abs> +notation for the norm. + +=head1 OPERATIONS + +Given the following notations: + + z1 = a + bi = r1 * exp(i * t1) + z2 = c + di = r2 * exp(i * t2) + z = <any complex or real number> + +the following (overloaded) operations are supported on complex numbers: + + z1 + z2 = (a + c) + i(b + d) + z1 - z2 = (a - c) + i(b - d) + z1 * z2 = (r1 * r2) * exp(i * (t1 + t2)) + z1 / z2 = (r1 / r2) * exp(i * (t1 - t2)) + z1 ** z2 = exp(z2 * log z1) + ~z = a - bi + abs(z) = r1 = sqrt(a*a + b*b) + sqrt(z) = sqrt(r1) * exp(i * t/2) + exp(z) = exp(a) * exp(i * b) + log(z) = log(r1) + i*t + sin(z) = 1/2i (exp(i * z1) - exp(-i * z)) + cos(z) = 1/2 (exp(i * z1) + exp(-i * z)) + atan2(z1, z2) = atan(z1/z2) + +The following extra operations are supported on both real and complex +numbers: + + Re(z) = a + Im(z) = b + arg(z) = t + abs(z) = r + + cbrt(z) = z ** (1/3) + log10(z) = log(z) / log(10) + logn(z, n) = log(z) / log(n) + + tan(z) = sin(z) / cos(z) + + csc(z) = 1 / sin(z) + sec(z) = 1 / cos(z) + cot(z) = 1 / tan(z) + + asin(z) = -i * log(i*z + sqrt(1-z*z)) + acos(z) = -i * log(z + i*sqrt(1-z*z)) + atan(z) = i/2 * log((i+z) / (i-z)) + + acsc(z) = asin(1 / z) + asec(z) = acos(1 / z) + acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i)) + + sinh(z) = 1/2 (exp(z) - exp(-z)) + cosh(z) = 1/2 (exp(z) + exp(-z)) + tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z)) + + csch(z) = 1 / sinh(z) + sech(z) = 1 / cosh(z) + coth(z) = 1 / tanh(z) + + asinh(z) = log(z + sqrt(z*z+1)) + acosh(z) = log(z + sqrt(z*z-1)) + atanh(z) = 1/2 * log((1+z) / (1-z)) + + acsch(z) = asinh(1 / z) + asech(z) = acosh(1 / z) + acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1)) + +I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, +I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>, +I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>, +I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>, +C<rho>, and C<theta> can be used also also mutators. The C<cbrt> +returns only one of the solutions: if you want all three, use the +C<root> function. + +The I<root> function is available to compute all the I<n> +roots of some complex, where I<n> is a strictly positive integer. +There are exactly I<n> such roots, returned as a list. Getting the +number mathematicians call C<j> such that: + + 1 + j + j*j = 0; + +is a simple matter of writing: + + $j = ((root(1, 3))[1]; + +The I<k>th root for C<z = [r,t]> is given by: + + (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n) + +The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In +order to ensure its restriction to real numbers is conform to what you +would expect, the comparison is run on the real part of the complex +number first, and imaginary parts are compared only when the real +parts match. + +=head1 CREATION + +To create a complex number, use either: + + $z = Math::Complex->make(3, 4); + $z = cplx(3, 4); + +if you know the cartesian form of the number, or + + $z = 3 + 4*i; + +if you like. To create a number using the polar form, use either: + + $z = Math::Complex->emake(5, pi/3); + $x = cplxe(5, pi/3); + +instead. The first argument is the modulus, the second is the angle +(in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a +notation for complex numbers in the polar form). + +It is possible to write: + + $x = cplxe(-3, pi/4); + +but that will be silently converted into C<[3,-3pi/4]>, since the modulus +must be non-negative (it represents the distance to the origin in the complex +plane). + +It is also possible to have a complex number as either argument of +either the C<make> or C<emake>: the appropriate component of +the argument will be used. + + $z1 = cplx(-2, 1); + $z2 = cplx($z1, 4); + +=head1 STRINGIFICATION + +When printed, a complex number is usually shown under its cartesian +form I<a+bi>, but there are legitimate cases where the polar format +I<[r,t]> is more appropriate. + +By calling the routine C<Math::Complex::display_format> and supplying either +C<"polar"> or C<"cartesian">, you override the default display format, +which is C<"cartesian">. Not supplying any argument returns the current +setting. + +This default can be overridden on a per-number basis by calling the +C<display_format> method instead. As before, not supplying any argument +returns the current display format for this number. Otherwise whatever you +specify will be the new display format for I<this> particular number. + +For instance: + + use Math::Complex; + + Math::Complex::display_format('polar'); + $j = ((root(1, 3))[1]; + print "j = $j\n"; # Prints "j = [1,2pi/3] + $j->display_format('cartesian'); + print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i" + +The polar format attempts to emphasize arguments like I<k*pi/n> +(where I<n> is a positive integer and I<k> an integer within [-9,+9]). + +=head1 USAGE + +Thanks to overloading, the handling of arithmetics with complex numbers +is simple and almost transparent. + +Here are some examples: + + use Math::Complex; + + $j = cplxe(1, 2*pi/3); # $j ** 3 == 1 + print "j = $j, j**3 = ", $j ** 3, "\n"; + print "1 + j + j**2 = ", 1 + $j + $j**2, "\n"; + + $z = -16 + 0*i; # Force it to be a complex + print "sqrt($z) = ", sqrt($z), "\n"; + + $k = exp(i * 2*pi/3); + print "$j - $k = ", $j - $k, "\n"; + + $z->Re(3); # Re, Im, arg, abs, + $j->arg(2); # (the last two aka rho, theta) + # can be used also as mutators. + +=head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO + +The division (/) and the following functions + + log ln log10 logn + tan sec csc cot + atan asec acsc acot + tanh sech csch coth + atanh asech acsch acoth + +cannot be computed for all arguments because that would mean dividing +by zero or taking logarithm of zero. These situations cause fatal +runtime errors looking like this + + cot(0): Division by zero. + (Because in the definition of cot(0), the divisor sin(0) is 0) + Died at ... + +or + + atanh(-1): Logarithm of zero. + Died at... + +For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>, +C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the the +logarithmic functions and the C<atanh>, C<acoth>, the argument cannot +be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be +C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be +C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument +cannot be C<-i> (the negative imaginary unit). For the C<tan>, +C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k> +is any integer. + +Note that because we are operating on approximations of real numbers, +these errors can happen when merely `too close' to the singularities +listed above. For example C<tan(2*atan2(1,1)+1e-15)> will die of +division by zero. + +=head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS + +The C<make> and C<emake> accept both real and complex arguments. +When they cannot recognize the arguments they will die with error +messages like the following + + Math::Complex::make: Cannot take real part of ... + Math::Complex::make: Cannot take real part of ... + Math::Complex::emake: Cannot take rho of ... + Math::Complex::emake: Cannot take theta of ... + +=head1 BUGS + +Saying C<use Math::Complex;> exports many mathematical routines in the +caller environment and even overrides some (C<sqrt>, C<log>). +This is construed as a feature by the Authors, actually... ;-) + +All routines expect to be given real or complex numbers. Don't attempt to +use BigFloat, since Perl has currently no rule to disambiguate a '+' +operation (for instance) between two overloaded entities. + +In Cray UNICOS there is some strange numerical instability that results +in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware. +The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex. +Whatever it is, it does not manifest itself anywhere else where Perl runs. + +=head1 AUTHORS + +Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>> and +Jarkko Hietaniemi <F<jhi@iki.fi>>. + +Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>. + +=cut + +1; + +# eof diff --git a/contrib/perl5/lib/Math/Trig.pm b/contrib/perl5/lib/Math/Trig.pm new file mode 100644 index 0000000..b7b5d5d --- /dev/null +++ b/contrib/perl5/lib/Math/Trig.pm @@ -0,0 +1,419 @@ +# +# Trigonometric functions, mostly inherited from Math::Complex. +# -- Jarkko Hietaniemi, since April 1997 +# -- Raphael Manfredi, September 1996 (indirectly: because of Math::Complex) +# + +require Exporter; +package Math::Trig; + +use strict; + +use Math::Complex qw(:trig); + +use vars qw($VERSION $PACKAGE + @ISA + @EXPORT @EXPORT_OK %EXPORT_TAGS); + +@ISA = qw(Exporter); + +$VERSION = 1.00; + +my @angcnv = qw(rad2deg rad2grad + deg2rad deg2grad + grad2rad grad2deg); + +@EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}}, + @angcnv); + +my @rdlcnv = qw(cartesian_to_cylindrical + cartesian_to_spherical + cylindrical_to_cartesian + cylindrical_to_spherical + spherical_to_cartesian + spherical_to_cylindrical); + +@EXPORT_OK = (@rdlcnv, 'great_circle_distance'); + +%EXPORT_TAGS = ('radial' => [ @rdlcnv ]); + +use constant pi2 => 2 * pi; +use constant pip2 => pi / 2; +use constant DR => pi2/360; +use constant RD => 360/pi2; +use constant DG => 400/360; +use constant GD => 360/400; +use constant RG => 400/pi2; +use constant GR => pi2/400; + +# +# Truncating remainder. +# + +sub remt ($$) { + # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available. + $_[0] - $_[1] * int($_[0] / $_[1]); +} + +# +# Angle conversions. +# + +sub rad2deg ($) { remt(RD * $_[0], 360) } + +sub deg2rad ($) { remt(DR * $_[0], pi2) } + +sub grad2deg ($) { remt(GD * $_[0], 360) } + +sub deg2grad ($) { remt(DG * $_[0], 400) } + +sub rad2grad ($) { remt(RG * $_[0], 400) } + +sub grad2rad ($) { remt(GR * $_[0], pi2) } + +sub cartesian_to_spherical { + my ( $x, $y, $z ) = @_; + + my $rho = sqrt( $x * $x + $y * $y + $z * $z ); + + return ( $rho, + atan2( $y, $x ), + $rho ? acos( $z / $rho ) : 0 ); +} + +sub spherical_to_cartesian { + my ( $rho, $theta, $phi ) = @_; + + return ( $rho * cos( $theta ) * sin( $phi ), + $rho * sin( $theta ) * sin( $phi ), + $rho * cos( $phi ) ); +} + +sub spherical_to_cylindrical { + my ( $x, $y, $z ) = spherical_to_cartesian( @_ ); + + return ( sqrt( $x * $x + $y * $y ), $_[1], $z ); +} + +sub cartesian_to_cylindrical { + my ( $x, $y, $z ) = @_; + + return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z ); +} + +sub cylindrical_to_cartesian { + my ( $rho, $theta, $z ) = @_; + + return ( $rho * cos( $theta ), $rho * sin( $theta ), $z ); +} + +sub cylindrical_to_spherical { + return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) ); +} + +sub great_circle_distance { + my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_; + + $rho = 1 unless defined $rho; # Default to the unit sphere. + + my $lat0 = pip2 - $phi0; + my $lat1 = pip2 - $phi1; + + return $rho * + acos(cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) + + sin( $lat0 ) * sin( $lat1 ) ); +} + +=pod + +=head1 NAME + +Math::Trig - trigonometric functions + +=head1 SYNOPSIS + + use Math::Trig; + + $x = tan(0.9); + $y = acos(3.7); + $z = asin(2.4); + + $halfpi = pi/2; + + $rad = deg2rad(120); + +=head1 DESCRIPTION + +C<Math::Trig> defines many trigonometric functions not defined by the +core Perl which defines only the C<sin()> and C<cos()>. The constant +B<pi> is also defined as are a few convenience functions for angle +conversions. + +=head1 TRIGONOMETRIC FUNCTIONS + +The tangent + +=over 4 + +=item B<tan> + +=back + +The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot +are aliases) + +B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan> + +The arcus (also known as the inverse) functions of the sine, cosine, +and tangent + +B<asin>, B<acos>, B<atan> + +The principal value of the arc tangent of y/x + +B<atan2>(y, x) + +The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc +and acotan/acot are aliases) + +B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan> + +The hyperbolic sine, cosine, and tangent + +B<sinh>, B<cosh>, B<tanh> + +The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch +and cotanh/coth are aliases) + +B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh> + +The arcus (also known as the inverse) functions of the hyperbolic +sine, cosine, and tangent + +B<asinh>, B<acosh>, B<atanh> + +The arcus cofunctions of the hyperbolic sine, cosine, and tangent +(acsch/acosech and acoth/acotanh are aliases) + +B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh> + +The trigonometric constant B<pi> is also defined. + +$pi2 = 2 * B<pi>; + +=head2 ERRORS DUE TO DIVISION BY ZERO + +The following functions + + acoth + acsc + acsch + asec + asech + atanh + cot + coth + csc + csch + sec + sech + tan + tanh + +cannot be computed for all arguments because that would mean dividing +by zero or taking logarithm of zero. These situations cause fatal +runtime errors looking like this + + cot(0): Division by zero. + (Because in the definition of cot(0), the divisor sin(0) is 0) + Died at ... + +or + + atanh(-1): Logarithm of zero. + Died at... + +For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>, +C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the +C<atanh>, C<acoth>, the argument cannot be C<1> (one). For the +C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one). For the +C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k * +pi>, where I<k> is any integer. + +=head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS + +Please note that some of the trigonometric functions can break out +from the B<real axis> into the B<complex plane>. For example +C<asin(2)> has no definition for plain real numbers but it has +definition for complex numbers. + +In Perl terms this means that supplying the usual Perl numbers (also +known as scalars, please see L<perldata>) as input for the +trigonometric functions might produce as output results that no more +are simple real numbers: instead they are complex numbers. + +The C<Math::Trig> handles this by using the C<Math::Complex> package +which knows how to handle complex numbers, please see L<Math::Complex> +for more information. In practice you need not to worry about getting +complex numbers as results because the C<Math::Complex> takes care of +details like for example how to display complex numbers. For example: + + print asin(2), "\n"; + +should produce something like this (take or leave few last decimals): + + 1.5707963267949-1.31695789692482i + +That is, a complex number with the real part of approximately C<1.571> +and the imaginary part of approximately C<-1.317>. + +=head1 PLANE ANGLE CONVERSIONS + +(Plane, 2-dimensional) angles may be converted with the following functions. + + $radians = deg2rad($degrees); + $radians = grad2rad($gradians); + + $degrees = rad2deg($radians); + $degrees = grad2deg($gradians); + + $gradians = deg2grad($degrees); + $gradians = rad2grad($radians); + +The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians. + +=head1 RADIAL COORDINATE CONVERSIONS + +B<Radial coordinate systems> are the B<spherical> and the B<cylindrical> +systems, explained shortly in more detail. + +You can import radial coordinate conversion functions by using the +C<:radial> tag: + + use Math::Trig ':radial'; + + ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z); + ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z); + ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z); + ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z); + ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi); + ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi); + +B<All angles are in radians>. + +=head2 COORDINATE SYSTEMS + +B<Cartesian> coordinates are the usual rectangular I<(x, y, +z)>-coordinates. + +Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional +coordinates which define a point in three-dimensional space. They are +based on a sphere surface. The radius of the sphere is B<rho>, also +known as the I<radial> coordinate. The angle in the I<xy>-plane +(around the I<z>-axis) is B<theta>, also known as the I<azimuthal> +coordinate. The angle from the I<z>-axis is B<phi>, also known as the +I<polar> coordinate. The `North Pole' is therefore I<0, 0, rho>, and +the `Bay of Guinea' (think of the missing big chunk of Africa) I<0, +pi/2, rho>. + +B<Beware>: some texts define I<theta> and I<phi> the other way round, +some texts define the I<phi> to start from the horizontal plane, some +texts use I<r> in place of I<rho>. + +Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional +coordinates which define a point in three-dimensional space. They are +based on a cylinder surface. The radius of the cylinder is B<rho>, +also known as the I<radial> coordinate. The angle in the I<xy>-plane +(around the I<z>-axis) is B<theta>, also known as the I<azimuthal> +coordinate. The third coordinate is the I<z>, pointing up from the +B<theta>-plane. + +=head2 3-D ANGLE CONVERSIONS + +Conversions to and from spherical and cylindrical coordinates are +available. Please notice that the conversions are not necessarily +reversible because of the equalities like I<pi> angles being equal to +I<-pi> angles. + +=over 4 + +=item cartesian_to_cylindrical + + ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z); + +=item cartesian_to_spherical + + ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z); + +=item cylindrical_to_cartesian + + ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z); + +=item cylindrical_to_spherical + + ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z); + +Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>. + +=item spherical_to_cartesian + + ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi); + +=item spherical_to_cylindrical + + ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi); + +Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>. + +=back + +=head1 GREAT CIRCLE DISTANCES + +You can compute spherical distances, called B<great circle distances>, +by importing the C<great_circle_distance> function: + + use Math::Trig 'great_circle_distance' + + $distance = great_circle_distance($theta0, $phi0, $theta1, $phi, [, $rho]); + +The I<great circle distance> is the shortest distance between two +points on a sphere. The distance is in C<$rho> units. The C<$rho> is +optional, it defaults to 1 (the unit sphere), therefore the distance +defaults to radians. + +=head1 EXAMPLES + +To calculate the distance between London (51.3N 0.5W) and Tokyo (35.7N +139.8E) in kilometers: + + use Math::Trig qw(great_circle_distance deg2rad); + + # Notice the 90 - latitude: phi zero is at the North Pole. + @L = (deg2rad(-0.5), deg2rad(90 - 51.3)); + @T = (deg2rad(139.8),deg2rad(90 - 35.7)); + + $km = great_circle_distance(@L, @T, 6378); + +The answer may be off by up to 0.3% because of the irregular (slightly +aspherical) form of the Earth. + +=head1 BUGS + +Saying C<use Math::Trig;> exports many mathematical routines in the +caller environment and even overrides some (C<sin>, C<cos>). This is +construed as a feature by the Authors, actually... ;-) + +The code is not optimized for speed, especially because we use +C<Math::Complex> and thus go quite near complex numbers while doing +the computations even when the arguments are not. This, however, +cannot be completely avoided if we want things like C<asin(2)> to give +an answer instead of giving a fatal runtime error. + +=head1 AUTHORS + +Jarkko Hietaniemi <F<jhi@iki.fi>> and +Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>>. + +=cut + +# eof |
