diff options
Diffstat (limited to 'contrib/perl5/lib/Math')
| -rw-r--r-- | contrib/perl5/lib/Math/BigFloat.pm | 398 | ||||
| -rw-r--r-- | contrib/perl5/lib/Math/BigInt.pm | 519 | ||||
| -rw-r--r-- | contrib/perl5/lib/Math/Complex.pm | 1889 | ||||
| -rw-r--r-- | contrib/perl5/lib/Math/Trig.pm | 456 |
4 files changed, 0 insertions, 3262 deletions
diff --git a/contrib/perl5/lib/Math/BigFloat.pm b/contrib/perl5/lib/Math/BigFloat.pm deleted file mode 100644 index 1eefac2..0000000 --- a/contrib/perl5/lib/Math/BigFloat.pm +++ /dev/null @@ -1,398 +0,0 @@ -package Math::BigFloat; - -use Math::BigInt; - -use Exporter; # just for use to be happy -@ISA = (Exporter); -$VERSION = '0.02'; - -use overload -'+' => sub {new Math::BigFloat &fadd}, -'-' => sub {new Math::BigFloat - $_[2]? fsub($_[1],${$_[0]}) : fsub(${$_[0]},$_[1])}, -'<=>' => sub {$_[2]? fcmp($_[1],${$_[0]}) : fcmp(${$_[0]},$_[1])}, -'cmp' => sub {$_[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])}, -'%' => sub {new Math::BigFloat - $_[2]? scalar fmod($_[1],${$_[0]}) : - scalar fmod(${$_[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); - 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 ( defined $e ) - { - 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 - no warnings; # $4 and $5 below might legitimately be undefined - 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) = @_; - $exp = 0 unless defined $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), - Math::BigInt::babs($ym)), - $xe-$ye-$scale); - } -} - -# modular division -# args are dividend, divisor -sub fmod #(fnum_str, fnum_str) return fnum_str -{ - local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1])); - 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); - if ( $xe < $ye ) - { - $ym .= ('0' x ($ye-$xe)); - } - else - { - $xm .= ('0' x ($xe-$ye)); - } - &norm(Math::BigInt::bmod($xm,$ym)); - } -} -# 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"."0" x length(substr($xm,$[+$scale+1))), - $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) { - # The first substr preserves the sign, passing a non- - # normalized "-0" to &round when rounding -0.006 (for - # example), purely so &round won't lose the sign. - &norm(&round(substr($xm,$[,1).'0', - "+0".substr($xm,$[+1), - "+1"."0" x length(substr($xm,$[+1))), $scale); - } else { - &norm(&round(substr($xm,$[,$xe), - "+0".substr($xm,$[+$xe), - "+1"."0" x length(substr($xm,$[+$xe))), $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 { - local($xm,$xe,$ym,$ye) = split('E', $x."E$y"); - if ($xm eq '+0' || $ym eq '+0') { - return $xm <=> $ym; - } - if ( $xe < $ye ) # adjust the exponents to be equal - { - $ym .= '0' x ($ye - $xe); - $ye = $xe; - } - elsif ( $ye < $xe ) # same here - { - $xm .= '0' x ($xe - $ye); - $xe = $ye; - } - return 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->fmod(NSTR) returns NSTR modular remainder - $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 embedded 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($Math::BigFloat::div_scale,length(dividend)+length(divisor))> -digits by default. -Also used for default sqrt scale. - -=item Rounding is performed - -according to the value of -C<$Math::BigFloat::rnd_mode>: - - trunc truncate the value - zero round towards 0 - +inf round towards +infinity (round up) - -inf round towards -infinity (round down) - even round to the nearest, .5 to the even digit - odd round to the nearest, .5 to the odd digit - -The default is C<even> rounding. - -=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. - -The printf subroutine does not use the value of -C<$Math::BigFloat::rnd_mode> when rounding values for printing. -Consequently, the way to print rounded values is -to specify the number of digits both as an -argument to C<ffround> and in the C<%f> printf string, -as follows: - - printf "%.3f\n", $bigfloat->ffround(-3); - -=head1 AUTHOR - -Mark Biggar -Patches by John Peacock Apr 2001 -=cut diff --git a/contrib/perl5/lib/Math/BigInt.pm b/contrib/perl5/lib/Math/BigInt.pm deleted file mode 100644 index 066577d..0000000 --- a/contrib/perl5/lib/Math/BigInt.pm +++ /dev/null @@ -1,519 +0,0 @@ -package Math::BigInt; -$VERSION='0.01'; - -use overload -'+' => sub {new Math::BigInt &badd}, -'-' => sub {new Math::BigInt - $_[2]? bsub($_[1],${$_[0]}) : bsub(${$_[0]},$_[1])}, -'<=>' => sub {$_[2]? bcmp($_[1],${$_[0]}) : bcmp(${$_[0]},$_[1])}, -'cmp' => sub {$_[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}, -'<<' => sub {new Math::BigInt - $_[2]? blsft($_[1],${$_[0]}) : blsft(${$_[0]},$_[1])}, -'>>' => sub {new Math::BigInt - $_[2]? brsft($_[1],${$_[0]}) : brsft(${$_[0]},$_[1])}, -'&' => sub {new Math::BigInt &band}, -'|' => sub {new Math::BigInt &bior}, -'^' => sub {new Math::BigInt &bxor}, -'~' => sub {new Math::BigInt &bnot}, - -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; - -# overcome a floating point problem on certain osnames (posix-bc, os390) -BEGIN { - my $x = 100000.0; - my $use_mult = int($x*1e-5)*1e5 == $x ? 1 : 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; - if ($use_mult) { - $prod[$cty++] = - $prod - ($car = int($prod * 1e-5)) * 1e5; - } - else { - $prod[$cty++] = - $prod - ($car = int($prod / 1e5)) * 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; - if ($use_mult) { - $x -= ($car = int($x * 1e-5)) * 1e5; - } - else { - $x -= ($car = int($x / 1e5)) * 1e5; - } - } - push(@x, $car); $car = 0; - for $y (@y) { - $y = $y * $dd + $car; - if ($use_mult) { - $y -= ($car = int($y * 1e-5)) * 1e5; - } - else { - $y -= ($car = int($y / 1e5)) * 1e5; - } - } - } - else { - push(@x, 0); - } - @q = (); ($v2,$v1) = @y[-2,-1]; - $v2 = 0 unless $v2; - while ($#x > $#y) { - ($u2,$u1,$u0) = @x[-3..-1]; - $u2 = 0 unless $u2; - $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; - if ($use_mult) { - $prd -= ($car = int($prd * 1e-5)) * 1e5; - } - else { - $prd -= ($car = int($prd / 1e5)) * 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); - } -} - -# compute x << y, y >= 0 -sub blsft { #(num_str, num_str) return num_str - &bmul($_[$[], &bpow(2, $_[$[+1])); -} - -# compute x >> y, y >= 0 -sub brsft { #(num_str, num_str) return num_str - &bdiv($_[$[], &bpow(2, $_[$[+1])); -} - -# compute x & y -sub band { #(num_str, num_str) return num_str - local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1); - if ($x eq 'NaN' || $y eq 'NaN') { - 'NaN'; - } else { - while ($x ne '+0' && $y ne '+0') { - ($x, $xr) = &bdiv($x, 0x10000); - ($y, $yr) = &bdiv($y, 0x10000); - $r = &badd(&bmul(int $xr & $yr, $m), $r); - $m = &bmul($m, 0x10000); - } - $r; - } -} - -# compute x | y -sub bior { #(num_str, num_str) return num_str - local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1); - if ($x eq 'NaN' || $y eq 'NaN') { - 'NaN'; - } else { - while ($x ne '+0' || $y ne '+0') { - ($x, $xr) = &bdiv($x, 0x10000); - ($y, $yr) = &bdiv($y, 0x10000); - $r = &badd(&bmul(int $xr | $yr, $m), $r); - $m = &bmul($m, 0x10000); - } - $r; - } -} - -# compute x ^ y -sub bxor { #(num_str, num_str) return num_str - local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1); - if ($x eq 'NaN' || $y eq 'NaN') { - 'NaN'; - } else { - while ($x ne '+0' || $y ne '+0') { - ($x, $xr) = &bdiv($x, 0x10000); - ($y, $yr) = &bdiv($y, 0x10000); - $r = &badd(&bmul(int $xr ^ $yr, $m), $r); - $m = &bmul($m, 0x10000); - } - $r; - } -} - -# represent ~x as twos-complement number -sub bnot { #(num_str) return num_str - &bsub(-1,$_[$[]); -} - -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 - $i->blsft(BINT) return BINT left shift - $i->brsft(BINT) return (BINT,BINT) right shift (quo,rem) just quo if scalar - $i->band(BINT) return BINT bit-wise and - $i->bior(BINT) return BINT bit-wise inclusive or - $i->bxor(BINT) return BINT bit-wise exclusive or - $i->bnot return BINT bit-wise not - -=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 conversion of -constants the expression 2**100 will be calculated 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 deleted file mode 100644 index 9812513..0000000 --- a/contrib/perl5/lib/Math/Complex.pm +++ /dev/null @@ -1,1889 +0,0 @@ -# -# Complex numbers and associated mathematical functions -# -- Raphael Manfredi Since Sep 1996 -# -- Jarkko Hietaniemi Since Mar 1997 -# -- Daniel S. Lewart Since Sep 1997 -# - -package Math::Complex; - -our($VERSION, @ISA, @EXPORT, %EXPORT_TAGS, $Inf); - -$VERSION = 1.31; - -BEGIN { - unless ($^O eq 'unicosmk') { - my $e = $!; - # We do want an arithmetic overflow, Inf INF inf Infinity:. - undef $Inf unless eval <<'EOE' and $Inf =~ /^inf(?:inity)?$/i; - local $SIG{FPE} = sub {die}; - my $t = CORE::exp 30; - $Inf = CORE::exp $t; -EOE - if (!defined $Inf) { # Try a different method - undef $Inf unless eval <<'EOE' and $Inf =~ /^inf(?:inity)?$/i; - local $SIG{FPE} = sub {die}; - my $t = 1; - $Inf = $t + "1e99999999999999999999999999999999"; -EOE - } - $! = $e; # Clear ERANGE. - } - $Inf = "Inf" if !defined $Inf || !($Inf > 0); # Desperation. -} - -use strict; - -my $i; -my %LOGN; - -require Exporter; - -@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, - '==' => \&numeq, - '<=>' => \&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 %DISPLAY_FORMAT = ('style' => 'cartesian', - 'polar_pretty_print' => 1); -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 -# -sub pi () { 4 * CORE::atan2(1, 1) } - -# -# pit2 -# -# The full circle -# -sub pit2 () { 2 * pi } - -# -# pip2 -# -# The quarter circle -# -sub pip2 () { pi / 2 } - -# -# deg1 -# -# One degree in radians, used in stringify_polar. -# - -sub deg1 () { pi / 180 } - -# -# uplog10 -# -# Used in log10(). -# -sub 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; -} - -# -# ip2 -# -# Half of i. -# -sub ip2 () { i / 2 } - -# -# 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); - } - } - } -} - -# -# (power) -# -# Computes z1**z2 = exp(z2 * log z1)). -# -sub power { - my ($z1, $z2, $inverted) = @_; - if ($inverted) { - return 1 if $z1 == 0 || $z2 == 1; - return 0 if $z2 == 0 && Re($z1) > 0; - } else { - return 1 if $z2 == 0 || $z1 == 1; - return 0 if $z1 == 0 && Re($z2) > 0; - } - my $w = $inverted ? &exp($z1 * &log($z2)) - : &exp($z2 * &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); -} - -# -# (numeq) -# -# Computes z1 == z2. -# -# (Required in addition to spaceship() because of NaNs.) -sub numeq { - my ($z1, $z2, $inverted) = @_; - my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0); - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); - return $re1 == $re2 && $im1 == $im2 ? 1 : 0; -} - -# -# (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) = @_; - unless (ref $z) { - if (@_ == 2) { - $_[0] = $_[1]; - } else { - return CORE::abs($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 0 if $r == 0; - return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3); -} - -# -# _rootbad -# -# Die on bad root. -# -sub _rootbad { - my $mess = "Root $_[0] illegal, root rank 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 0 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 &log($z) / $logn; -} - -# -# (cos) -# -# Compute cos(z) = (exp(iz) + exp(-iz))/2. -# -sub cos { - my ($z) = @_; - return CORE::cos($z) unless ref $z; - my ($x, $y) = @{$z->cartesian}; - my $ey = CORE::exp($y); - my $sx = CORE::sin($x); - my $cx = CORE::cos($x); - my $ey_1 = $ey ? 1 / $ey : $Inf; - return (ref $z)->make($cx * ($ey + $ey_1)/2, - $sx * ($ey_1 - $ey)/2); -} - -# -# (sin) -# -# Compute sin(z) = (exp(iz) - exp(-iz))/2. -# -sub sin { - my ($z) = @_; - return CORE::sin($z) unless ref $z; - my ($x, $y) = @{$z->cartesian}; - my $ey = CORE::exp($y); - my $sx = CORE::sin($x); - my $cx = CORE::cos($x); - my $ey_1 = $ey ? 1 / $ey : $Inf; - return (ref $z)->make($sx * ($ey + $ey_1)/2, - $cx * ($ey - $ey_1)/2); -} - -# -# tan -# -# Compute tan(z) = sin(z) / cos(z). -# -sub tan { - my ($z) = @_; - my $cz = &cos($z); - _divbyzero "tan($z)", "cos($z)" if $cz == 0; - return &sin($z) / $cz; -} - -# -# sec -# -# Computes the secant sec(z) = 1 / cos(z). -# -sub sec { - my ($z) = @_; - my $cz = &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 = &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 = &sin($z); - _divbyzero "cot($z)", "sin($z)" if ($sz == 0); - return &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; - $z = cplx($z, 0) unless ref $z; - my ($x, $y) = @{$z->cartesian}; - return 0 if $x == 1 && $y == 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 (ref $z)->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; - $z = cplx($z, 0) unless ref $z; - my ($x, $y) = @{$z->cartesian}; - return 0 if $x == 0 && $y == 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 (ref $z)->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; - my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0); - return 0 if $x == 0 && $y == 0; - _divbyzero "atan(i)" if ( $z == i); - _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test... - my $log = &log((i + $z) / (i - $z)); - 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 $z == 0; - return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z) - unless ref $z; - _divbyzero "acot(i)" if ($z - i == 0); - _logofzero "acot(-i)" if ($z + i == 0); - 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 ? ($ex + 1/$ex)/2 : $Inf; - } - my ($x, $y) = @{$z->cartesian}; - $ex = CORE::exp($x); - my $ex_1 = $ex ? 1 / $ex : $Inf; - 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) { - return 0 if $z == 0; - $ex = CORE::exp($z); - return $ex ? ($ex - 1/$ex)/2 : "-$Inf"; - } - my ($x, $y) = @{$z->cartesian}; - my $cy = CORE::cos($y); - my $sy = CORE::sin($y); - $ex = CORE::exp($x); - my $ex_1 = $ex ? 1 / $ex : $Inf; - 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) { - $z = cplx($z, 0); - } - my ($re, $im) = @{$z->cartesian}; - if ($im == 0) { - return CORE::log($re + CORE::sqrt($re*$re - 1)) - if $re >= 1; - return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re)) - if CORE::abs($re) < 1; - } - my $t = &sqrt($z * $z - 1) + $z; - # Try Taylor if looking bad (this usually means that - # $z was large negative, therefore the sqrt is really - # close to abs(z), summing that with z...) - $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7) - if $t == 0; - my $u = &log($t); - $u->Im(-$u->Im) if $re < 0 && $im == 0; - return $re < 0 ? -$u : $u; -} - -# -# asinh -# -# Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z+1)) -# -sub asinh { - my ($z) = @_; - unless (ref $z) { - my $t = $z + CORE::sqrt($z*$z + 1); - return CORE::log($t) if $t; - } - my $t = &sqrt($z * $z + 1) + $z; - # Try Taylor if looking bad (this usually means that - # $z was large negative, therefore the sqrt is really - # close to abs(z), summing that with z...) - $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7) - if $t == 0; - return &log($t); -} - -# -# 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 (1 - $z == 0); - _logofzero 'atanh(-1)' if (1 + $z == 0); - return 0.5 * &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 ($z == 0); - 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 ($z - 1 == 0); - _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0); - return &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 CORE::atan2($re1, $re2) if $im1 == 0; - return ($im1<=>0) * pip2 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 (get if no argument) the display format for all complex numbers that -# don't happen to have overridden it via ->display_format -# -# When called as an object 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 %display_format = %DISPLAY_FORMAT; - - if (ref $self) { # Called as an object method - if (exists $self->{display_format}) { - my %obj = %{$self->{display_format}}; - @display_format{keys %obj} = values %obj; - } - } - if (@_ == 1) { - $display_format{style} = shift; - } else { - my %new = @_; - @display_format{keys %new} = values %new; - } - - if (ref $self) { # Called as an object method - $self->{display_format} = { %display_format }; - return - wantarray ? - %{$self->{display_format}} : - $self->{display_format}->{style}; - } - - # Called as a class method - %DISPLAY_FORMAT = %display_format; - return - wantarray ? - %DISPLAY_FORMAT : - $DISPLAY_FORMAT{style}; -} - -# -# (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 $style = $z->display_format; - - $style = $DISPLAY_FORMAT{style} unless defined $style; - - return $z->stringify_polar if $style =~ /^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); - - my %format = $z->display_format; - my $format = $format{format}; - - if ($x) { - if ($x =~ /^NaN[QS]?$/i) { - $re = $x; - } else { - if ($x =~ /^-?$Inf$/oi) { - $re = $x; - } else { - $re = defined $format ? sprintf($format, $x) : $x; - } - } - } else { - undef $re; - } - - if ($y) { - if ($y =~ /^(NaN[QS]?)$/i) { - $im = $y; - } else { - if ($y =~ /^-?$Inf$/oi) { - $im = $y; - } else { - $im = - defined $format ? - sprintf($format, $y) : - ($y == 1 ? "" : ($y == -1 ? "-" : $y)); - } - } - $im .= "i"; - } else { - undef $im; - } - - my $str = $re; - - if (defined $im) { - if ($y < 0) { - $str .= $im; - } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i) { - $str .= "+" if defined $re; - $str .= $im; - } - } elsif (!defined $re) { - $str = "0"; - } - - return $str; -} - - -# -# ->stringify_polar -# -# Stringify as a polar representation '[r,t]'. -# -sub stringify_polar { - my $z = shift; - my ($r, $t) = @{$z->polar}; - my $theta; - - my %format = $z->display_format; - my $format = $format{format}; - - if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?$Inf$/oi) { - $theta = $t; - } elsif ($t == pi) { - $theta = "pi"; - } elsif ($r == 0 || $t == 0) { - $theta = defined $format ? sprintf($format, $t) : $t; - } - - return "[$r,$theta]" if defined $theta; - - # - # Try to identify pi/n and friends. - # - - $t -= int(CORE::abs($t) / pit2) * pit2; - - if ($format{polar_pretty_print} && $t) { - my ($a, $b); - for $a (2..9) { - $b = $t * $a / pi; - if ($b =~ /^-?\d+$/) { - $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1; - $theta = "${b}pi/$a"; - last; - } - } - } - - if (defined $format) { - $r = sprintf($format, $r); - $theta = sprintf($format, $theta) unless defined $theta; - } else { - $theta = $t unless defined $theta; - } - - return "[$r,$theta]"; -} - -1; -__END__ - -=pod - -=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 -style I<a+bi>, but there are legitimate cases where the polar style -I<[r,t]> is more appropriate. - -By calling the class method C<Math::Complex::display_format> and -supplying either C<"polar"> or C<"cartesian"> as an argument, you -override the default display style, which is C<"cartesian">. Not -supplying any argument returns the current settings. - -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 style for this number. Otherwise whatever you -specify will be the new display style 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 style attempts to emphasize arguments like I<k*pi/n> -(where I<n> is a positive integer and I<k> an integer within [-9, +9]), -this is called I<polar pretty-printing>. - -=head2 CHANGED IN PERL 5.6 - -The C<display_format> class method and the corresponding -C<display_format> object method can now be called using -a parameter hash instead of just a one parameter. - -The old display format style, which can have values C<"cartesian"> or -C<"polar">, can be changed using the C<"style"> parameter. - - $j->display_format(style => "polar"); - -The one parameter calling convention also still works. - - $j->display_format("polar"); - -There are two new display parameters. - -The first one is C<"format">, which is a sprintf()-style format string -to be used for both numeric parts of the complex number(s). The is -somewhat system-dependent but most often it corresponds to C<"%.15g">. -You can revert to the default by setting the C<format> to C<undef>. - - # the $j from the above example - - $j->display_format('format' => '%.5f'); - print "j = $j\n"; # Prints "j = -0.50000+0.86603i" - $j->display_format('format' => undef); - print "j = $j\n"; # Prints "j = -0.5+0.86603i" - -Notice that this affects also the return values of the -C<display_format> methods: in list context the whole parameter hash -will be returned, as opposed to only the style parameter value. -This is a potential incompatibility with earlier versions if you -have been calling the C<display_format> method in list context. - -The second new display parameter is C<"polar_pretty_print">, which can -be set to true or false, the default being true. See the previous -section for what this means. - -=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. - -=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@pobox.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 deleted file mode 100644 index b28f150..0000000 --- a/contrib/perl5/lib/Math/Trig.pm +++ /dev/null @@ -1,456 +0,0 @@ -# -# 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 5.005_64; -use strict; - -use Math::Complex qw(:trig); - -our($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 ]); - -sub pi2 () { 2 * pi } -sub pip2 () { pi / 2 } - -sub DR () { pi2/360 } -sub RD () { 360/pi2 } -sub DG () { 400/360 } -sub GD () { 360/400 } -sub RG () { 400/pi2 } -sub 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 rad2rad($) { remt($_[0], pi2) } - -sub deg2deg($) { remt($_[0], 360) } - -sub grad2grad($) { remt($_[0], 400) } - -sub rad2deg ($;$) { my $d = RD * $_[0]; $_[1] ? $d : deg2deg($d) } - -sub deg2rad ($;$) { my $d = DR * $_[0]; $_[1] ? $d : rad2rad($d) } - -sub grad2deg ($;$) { my $d = GD * $_[0]; $_[1] ? $d : deg2deg($d) } - -sub deg2grad ($;$) { my $d = DG * $_[0]; $_[1] ? $d : grad2grad($d) } - -sub rad2grad ($;$) { my $d = RG * $_[0]; $_[1] ? $d : grad2grad($d) } - -sub grad2rad ($;$) { my $d = GR * $_[0]; $_[1] ? $d : rad2rad($d) } - -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. -The result is by default wrapped to be inside the [0, {2pi,360,400}[ circle. -If you don't want this, supply a true second argument: - - $zillions_of_radians = deg2rad($zillions_of_degrees, 1); - $negative_degrees = rad2deg($negative_radians, 1); - -You can also do the wrapping explicitly by rad2rad(), deg2deg(), and -grad2grad(). - -=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>. In geographical terms I<phi> is latitude (northward -positive, southward negative) and I<theta> is longitude (eastward -positive, westward negative). - -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, $phi1, [, $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. - -If you think geographically the I<theta> are longitudes: zero at the -Greenwhich meridian, eastward positive, westward negative--and the -I<phi> are latitudes: zero at the North Pole, northward positive, -southward negative. B<NOTE>: this formula thinks in mathematics, not -geographically: the I<phi> zero is at the North Pole, not at the -Equator on the west coast of Africa (Bay of Guinea). You need to -subtract your geographical coordinates from I<pi/2> (also known as 90 -degrees). - - $distance = great_circle_distance($lon0, pi/2 - $lat0, - $lon1, pi/2 - $lat1, $rho); - -=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 few percentages because of the irregular -(slightly aspherical) form of the Earth. The used formula - - lat0 = 90 degrees - phi0 - lat1 = 90 degrees - phi1 - d = R * arccos(cos(lat0) * cos(lat1) * cos(lon1 - lon01) + - sin(lat0) * sin(lat1)) - -is also somewhat unreliable for small distances (for locations -separated less than about five degrees) because it uses arc cosine -which is rather ill-conditioned for values close to zero. - -=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@pobox.com>>. - -=cut - -# eof |
