lib/Math/BigFloat.pm

   1 package Math::BigFloat;
   2
   3 use Math::BigInt;
   4
   5 use Exporter;  # just for use to be happy
   6 @ISA = (Exporter);
   7
   8 %OVERLOAD = (
   9                                 # Anonymous subroutines:
  10 '+'     =>      sub {new BigFloat &fadd},
  11 '-'     =>      sub {new BigFloat
  12                        $_[2]? fsub($_[1],${$_[0]}) : fsub(${$_[0]},$_[1])},
  13 '<=>'   =>      sub {new BigFloat
  14                        $_[2]? fcmp($_[1],${$_[0]}) : fcmp(${$_[0]},$_[1])},
  15 'cmp'   =>      sub {new BigFloat
  16                        $_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])},
  17 '*'     =>      sub {new BigFloat &fmul},
  18 '/'     =>      sub {new BigFloat
  19                        $_[2]? scalar fdiv($_[1],${$_[0]}) :
  20                          scalar fdiv(${$_[0]},$_[1])},
  21 'neg'   =>      sub {new BigFloat &fneg},
  22 'abs'   =>      sub {new BigFloat &fabs},
  23
  24 qw(
  25 ""      stringify
  26 0+      numify)                 # Order of arguments unsignificant
  27 );
  28
  29 sub new {
  30   my $foo = fnorm($_[1]);
  31   panic("Not a number initialized to BigFloat") if $foo eq "NaN";
  32   bless \$foo;
  33 }
  34 sub numify { 0 + "${$_[0]}" }   # Not needed, additional overhead
  35                                 # comparing to direct compilation based on
  36                                 # stringify
  37 sub stringify {
  38     my $n = ${$_[0]};
  39
  40     $n =~ s/^\+//;
  41     $n =~ s/E//;
  42
  43     $n =~ s/([-+]\d+)$//;
  44
  45     my $e = $1;
  46     my $ln = length($n);
  47
  48     if ($e > 0) {
  49         $n .= "0" x $e . '.';
  50     } elsif (abs($e) < $ln) {
  51         substr($n, $ln + $e, 0) = '.';
  52     } else {
  53         $n = '.' . ("0" x (abs($e) - $ln)) . $n;
  54     }
  55
  56     # 1 while $n =~ s/(.*\d)(\d\d\d)/$1,$2/;
  57
  58     return $n;
  59 }
  60
  61 # Arbitrary length float math package
  62 #
  63 # by Mark Biggar
  64 #
  65 # number format
  66 #   canonical strings have the form /[+-]\d+E[+-]\d+/
  67 #   Input values can have inbedded whitespace
  68 # Error returns
  69 #   'NaN'           An input parameter was "Not a Number" or
  70 #                       divide by zero or sqrt of negative number
  71 # Division is computed to
  72 #   max($div_scale,length(dividend)+length(divisor))
  73 #   digits by default.
  74 # Also used for default sqrt scale
  75
  76 $div_scale = 40;
  77
  78 # Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
  79
  80 $rnd_mode = 'even';
  81
  82 sub fadd; sub fsub; sub fmul; sub fdiv;
  83 sub fneg; sub fabs; sub fcmp;
  84 sub fround; sub ffround;
  85 sub fnorm; sub fsqrt;
  86
  87 #   bigfloat routines
  88 #
  89 #   fadd(NSTR, NSTR) return NSTR            addition
  90 #   fsub(NSTR, NSTR) return NSTR            subtraction
  91 #   fmul(NSTR, NSTR) return NSTR            multiplication
  92 #   fdiv(NSTR, NSTR[,SCALE]) returns NSTR   division to SCALE places
  93 #   fneg(NSTR) return NSTR                  negation
  94 #   fabs(NSTR) return NSTR                  absolute value
  95 #   fcmp(NSTR,NSTR) return CODE             compare undef,<0,=0,>0
  96 #   fround(NSTR, SCALE) return NSTR         round to SCALE digits
  97 #   ffround(NSTR, SCALE) return NSTR        round at SCALEth place
  98 #   fnorm(NSTR) return (NSTR)               normalize
  99 #   fsqrt(NSTR[, SCALE]) return NSTR        sqrt to SCALE places
 100
 101 \f
 102 # Convert a number to canonical string form.
 103 #   Takes something that looks like a number and converts it to
 104 #   the form /^[+-]\d+E[+-]\d+$/.
 105 sub fnorm { #(string) return fnum_str
 106     local($_) = @_;
 107     s/\s+//g;                               # strip white space
 108     if (/^([+-]?)(\d*)(\.(\d*))?([Ee]([+-]?\d+))?$/ && "$2$4" ne '') {
 109         &norm(($1 ? "$1$2$4" : "+$2$4"),(($4 ne '') ? $6-length($4) : $6));
 110     } else {
 111         'NaN';
 112     }
 113 }
 114
 115 # normalize number -- for internal use
 116 sub norm { #(mantissa, exponent) return fnum_str
 117     local($_, $exp) = @_;
 118     if ($_ eq 'NaN') {
 119         'NaN';
 120     } else {
 121         s/^([+-])0+/$1/;                        # strip leading zeros
 122         if (length($_) == 1) {
 123             '+0E+0';
 124         } else {
 125             $exp += length($1) if (s/(0+)$//);  # strip trailing zeros
 126             sprintf("%sE%+ld", $_, $exp);
 127         }
 128     }
 129 }
 130
 131 # negation
 132 sub fneg { #(fnum_str) return fnum_str
 133     local($_) = fnorm($_[$[]);
 134     vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0E+0'; # flip sign
 135     s/^H/N/;
 136     $_;
 137 }
 138
 139 # absolute value
 140 sub fabs { #(fnum_str) return fnum_str
 141     local($_) = fnorm($_[$[]);
 142     s/^-/+/;                                   # mash sign
 143     $_;
 144 }
 145
 146 # multiplication
 147 sub fmul { #(fnum_str, fnum_str) return fnum_str
 148     local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1]));
 149     if ($x eq 'NaN' || $y eq 'NaN') {
 150         'NaN';
 151     } else {
 152         local($xm,$xe) = split('E',$x);
 153         local($ym,$ye) = split('E',$y);
 154         &norm(Math::BigInt::bmul($xm,$ym),$xe+$ye);
 155     }
 156 }
 157 \f
 158 # addition
 159 sub fadd { #(fnum_str, fnum_str) return fnum_str
 160     local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1]));
 161     if ($x eq 'NaN' || $y eq 'NaN') {
 162         'NaN';
 163     } else {
 164         local($xm,$xe) = split('E',$x);
 165         local($ym,$ye) = split('E',$y);
 166         ($xm,$xe,$ym,$ye) = ($ym,$ye,$xm,$xe) if ($xe < $ye);
 167         &norm(Math::BigInt::badd($ym,$xm.('0' x ($xe-$ye))),$ye);
 168     }
 169 }
 170
 171 # subtraction
 172 sub fsub { #(fnum_str, fnum_str) return fnum_str
 173     fadd($_[$[],fneg($_[$[+1]));
 174 }
 175
 176 # division
 177 #   args are dividend, divisor, scale (optional)
 178 #   result has at most max(scale, length(dividend), length(divisor)) digits
 179 sub fdiv #(fnum_str, fnum_str[,scale]) return fnum_str
 180 {
 181     local($x,$y,$scale) = (fnorm($_[$[]),fnorm($_[$[+1]),$_[$[+2]);
 182     if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0E+0') {
 183         'NaN';
 184     } else {
 185         local($xm,$xe) = split('E',$x);
 186         local($ym,$ye) = split('E',$y);
 187         $scale = $div_scale if (!$scale);
 188         $scale = length($xm)-1 if (length($xm)-1 > $scale);
 189         $scale = length($ym)-1 if (length($ym)-1 > $scale);
 190         $scale = $scale + length($ym) - length($xm);
 191         &norm(&round(Math::BigInt::bdiv($xm.('0' x $scale),$ym),$ym),
 192             $xe-$ye-$scale);
 193     }
 194 }
 195 \f
 196 # round int $q based on fraction $r/$base using $rnd_mode
 197 sub round { #(int_str, int_str, int_str) return int_str
 198     local($q,$r,$base) = @_;
 199     if ($q eq 'NaN' || $r eq 'NaN') {
 200         'NaN';
 201     } elsif ($rnd_mode eq 'trunc') {
 202         $q;                         # just truncate
 203     } else {
 204         local($cmp) = Math::BigInt::bcmp(Math::BigInt::bmul($r,'+2'),$base);
 205         if ( $cmp < 0 ||
 206                  ($cmp == 0 &&
 207                   ( $rnd_mode eq 'zero'                             ||
 208                    ($rnd_mode eq '-inf' && (substr($q,$[,1) eq '+')) ||
 209                    ($rnd_mode eq '+inf' && (substr($q,$[,1) eq '-')) ||
 210                    ($rnd_mode eq 'even' && $q =~ /[24680]$/)        ||
 211                    ($rnd_mode eq 'odd'  && $q =~ /[13579]$/)        )) ) {
 212             $q;                     # round down
 213         } else {
 214             Math::BigInt::badd($q, ((substr($q,$[,1) eq '-') ? '-1' : '+1'));
 215                                     # round up
 216         }
 217     }
 218 }
 219
 220 # round the mantissa of $x to $scale digits
 221 sub fround { #(fnum_str, scale) return fnum_str
 222     local($x,$scale) = (fnorm($_[$[]),$_[$[+1]);
 223     if ($x eq 'NaN' || $scale <= 0) {
 224         $x;
 225     } else {
 226         local($xm,$xe) = split('E',$x);
 227         if (length($xm)-1 <= $scale) {
 228             $x;
 229         } else {
 230             &norm(&round(substr($xm,$[,$scale+1),
 231                          "+0".substr($xm,$[+$scale+1,1),"+10"),
 232                   $xe+length($xm)-$scale-1);
 233         }
 234     }
 235 }
 236 \f
 237 # round $x at the 10 to the $scale digit place
 238 sub ffround { #(fnum_str, scale) return fnum_str
 239     local($x,$scale) = (fnorm($_[$[]),$_[$[+1]);
 240     if ($x eq 'NaN') {
 241         'NaN';
 242     } else {
 243         local($xm,$xe) = split('E',$x);
 244         if ($xe >= $scale) {
 245             $x;
 246         } else {
 247             $xe = length($xm)+$xe-$scale;
 248             if ($xe < 1) {
 249                 '+0E+0';
 250             } elsif ($xe == 1) {
 251                 &norm(&round('+0',"+0".substr($xm,$[+1,1),"+10"), $scale);
 252             } else {
 253                 &norm(&round(substr($xm,$[,$xe),
 254                       "+0".substr($xm,$[+$xe,1),"+10"), $scale);
 255             }
 256         }
 257     }
 258 }
 259
 260 # compare 2 values returns one of undef, <0, =0, >0
 261 #   returns undef if either or both input value are not numbers
 262 sub fcmp #(fnum_str, fnum_str) return cond_code
 263 {
 264     local($x, $y) = (fnorm($_[$[]),fnorm($_[$[+1]));
 265     if ($x eq "NaN" || $y eq "NaN") {
 266         undef;
 267     } else {
 268         ord($y) <=> ord($x)
 269         ||
 270         (  local($xm,$xe,$ym,$ye) = split('E', $x."E$y"),
 271              (($xe <=> $ye) * (substr($x,$[,1).'1')
 272              || Math::BigInt::cmp($xm,$ym))
 273         );
 274     }
 275 }
 276 \f
 277 # square root by Newtons method.
 278 sub fsqrt { #(fnum_str[, scale]) return fnum_str
 279     local($x, $scale) = (fnorm($_[$[]), $_[$[+1]);
 280     if ($x eq 'NaN' || $x =~ /^-/) {
 281         'NaN';
 282     } elsif ($x eq '+0E+0') {
 283         '+0E+0';
 284     } else {
 285         local($xm, $xe) = split('E',$x);
 286         $scale = $div_scale if (!$scale);
 287         $scale = length($xm)-1 if ($scale < length($xm)-1);
 288         local($gs, $guess) = (1, sprintf("1E%+d", (length($xm)+$xe-1)/2));
 289         while ($gs < 2*$scale) {
 290             $guess = fmul(fadd($guess,fdiv($x,$guess,$gs*2)),".5");
 291             $gs *= 2;
 292         }
 293         new BigFloat &fround($guess, $scale);
 294     }
 295 }
 296
 297 1;