projects / p5sagit/p5-mst-13.2.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | history | HEAD
A little better comments.
[p5sagit/p5-mst-13.2.git] / lib / Math / BigFloat.pm
CommitLineData
58cde26e 1#!/usr/bin/perl -w
2
58cde26e 3# The following hash values are internally used:
4# _e: exponent (BigInt)
5# _m: mantissa (absolute BigInt)
6# sign: +,-,"NaN" if not a number
7# _a: accuracy
8# _p: precision
0716bf9b 9# _f: flags, used to signal MBI not to touch our private parts
58cde26e 10
a0d0e21e 11package Math::BigFloat;
12
394e6ffb 13$VERSION = '1.26';
58cde26e 14require 5.005;
15use Exporter;
0716bf9b 16use Math::BigInt qw/objectify/;
58cde26e 17@ISA = qw( Exporter Math::BigInt);
394e6ffb 18#@EXPORT_OK = qw(
19# bcmp
20# badd bmul bdiv bmod bnorm bsub
21# bgcd blcm bround bfround
22# bpow bnan bzero bfloor bceil
23# bacmp bstr binc bdec binf
24# is_odd is_even is_nan is_inf is_positive is_negative
25# is_zero is_one sign
26# );
27
58cde26e 28use strict;
027dc388 29use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode/;
58cde26e 30my $class = "Math::BigFloat";
a0d0e21e 31
a5f75d66 32use overload
bd05a461 33'<=>' => sub { $_[2] ?
34 ref($_[0])->bcmp($_[1],$_[0]) :
35 ref($_[0])->bcmp($_[0],$_[1])},
0716bf9b 36'int' => sub { $_[0]->as_number() }, # 'trunc' to bigint
a5f75d66 37;
a0d0e21e 38
0716bf9b 39##############################################################################
40# global constants, flags and accessory
41
42use constant MB_NEVER_ROUND => 0x0001;
43
58cde26e 44# are NaNs ok?
45my $NaNOK=1;
58cde26e 46# constant for easier life
47my $nan = 'NaN';
58cde26e 48
ee15d750 49# class constants, use Class->constant_name() to access
50$round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
51$accuracy = undef;
52$precision = undef;
53$div_scale = 40;
58cde26e 54
027dc388 55##############################################################################
56# the old code had $rnd_mode, so we need to support it, too
57
58$rnd_mode = 'even';
59sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
60sub FETCH { return $round_mode; }
61sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
62
63BEGIN { tie $rnd_mode, 'Math::BigFloat'; }
64
65##############################################################################
66
574bacfe 67# in case we call SUPER::->foo() and this wants to call modify()
68# sub modify () { 0; }
69
58cde26e 70{
ee15d750 71 # valid method aliases for AUTOLOAD
58cde26e 72 my %methods = map { $_ => 1 }
73 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
394e6ffb 74 fint facmp fcmp fzero fnan finf finc fdec
75 fceil ffloor frsft flsft fone
ee15d750 76 /;
77 # valid method's that need to be hand-ed up (for AUTOLOAD)
78 my %hand_ups = map { $_ => 1 }
79 qw / is_nan is_inf is_negative is_positive
394e6ffb 80 accuracy precision div_scale round_mode fneg fabs babs fnot
58cde26e 81 /;
82
ee15d750 83 sub method_alias { return exists $methods{$_[0]||''}; }
84 sub method_hand_up { return exists $hand_ups{$_[0]||''}; }
a0d0e21e 85}
0e8b9368 86
58cde26e 87##############################################################################
88# constructors
a0d0e21e 89
58cde26e 90sub new
91 {
92 # create a new BigFloat object from a string or another bigfloat object.
93 # _e: exponent
94 # _m: mantissa
95 # sign => sign (+/-), or "NaN"
a0d0e21e 96
58cde26e 97 my $class = shift;
98
99 my $wanted = shift; # avoid numify call by not using || here
100 return $class->bzero() if !defined $wanted; # default to 0
101 return $wanted->copy() if ref($wanted) eq $class;
a0d0e21e 102
58cde26e 103 my $round = shift; $round = 0 if !defined $round; # no rounding as default
104 my $self = {}; bless $self, $class;
b22b3e31 105 # shortcut for bigints and its subclasses
0716bf9b 106 if ((ref($wanted)) && (ref($wanted) ne $class))
58cde26e 107 {
0716bf9b 108 $self->{_m} = $wanted->as_number(); # get us a bigint copy
027dc388 109 $self->{_e} = Math::BigInt->bzero();
58cde26e 110 $self->{_m}->babs();
111 $self->{sign} = $wanted->sign();
0716bf9b 112 return $self->bnorm();
58cde26e 113 }
114 # got string
115 # handle '+inf', '-inf' first
ee15d750 116 if ($wanted =~ /^[+-]?inf$/)
58cde26e 117 {
027dc388 118 $self->{_e} = Math::BigInt->bzero();
119 $self->{_m} = Math::BigInt->bzero();
58cde26e 120 $self->{sign} = $wanted;
ee15d750 121 $self->{sign} = '+inf' if $self->{sign} eq 'inf';
0716bf9b 122 return $self->bnorm();
58cde26e 123 }
124 #print "new string '$wanted'\n";
125 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
126 if (!ref $mis)
127 {
128 die "$wanted is not a number initialized to $class" if !$NaNOK;
027dc388 129 $self->{_e} = Math::BigInt->bzero();
130 $self->{_m} = Math::BigInt->bzero();
58cde26e 131 $self->{sign} = $nan;
132 }
133 else
134 {
135 # make integer from mantissa by adjusting exp, then convert to bigint
136 $self->{_e} = Math::BigInt->new("$$es$$ev"); # exponent
027dc388 137 $self->{_m} = Math::BigInt->new("$$miv$$mfv"); # create mantissa
58cde26e 138 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
027dc388 139 $self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0;
140 $self->{sign} = $$mis;
58cde26e 141 }
0716bf9b 142 #print "$wanted => $self->{sign} $self->{value}\n";
58cde26e 143 $self->bnorm(); # first normalize
144 # if any of the globals is set, round to them and thus store them insid $self
ee15d750 145 $self->round($accuracy,$precision,$class->round_mode)
58cde26e 146 if defined $accuracy || defined $precision;
147 return $self;
148 }
a0d0e21e 149
58cde26e 150sub bnan
151 {
152 # create a bigfloat 'NaN', if given a BigFloat, set it to 'NaN'
153 my $self = shift;
154 $self = $class if !defined $self;
155 if (!ref($self))
288d023a 156 {
58cde26e 157 my $c = $self; $self = {}; bless $self, $c;
a0d0e21e 158 }
574bacfe 159 $self->{_m} = Math::BigInt->bzero();
160 $self->{_e} = Math::BigInt->bzero();
58cde26e 161 $self->{sign} = $nan;
394e6ffb 162 ($self->{_a},$self->{_p}) = @_ if @_ > 0;
58cde26e 163 return $self;
164 }
a0d0e21e 165
58cde26e 166sub binf
167 {
168 # create a bigfloat '+-inf', if given a BigFloat, set it to '+-inf'
169 my $self = shift;
170 my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
a0d0e21e 171
58cde26e 172 $self = $class if !defined $self;
173 if (!ref($self))
174 {
175 my $c = $self; $self = {}; bless $self, $c;
176 }
574bacfe 177 $self->{_m} = Math::BigInt->bzero();
178 $self->{_e} = Math::BigInt->bzero();
58cde26e 179 $self->{sign} = $sign.'inf';
394e6ffb 180 ($self->{_a},$self->{_p}) = @_ if @_ > 0;
58cde26e 181 return $self;
182 }
a0d0e21e 183
574bacfe 184sub bone
185 {
186 # create a bigfloat '+-1', if given a BigFloat, set it to '+-1'
187 my $self = shift;
188 my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
189
190 $self = $class if !defined $self;
191 if (!ref($self))
192 {
193 my $c = $self; $self = {}; bless $self, $c;
194 }
195 $self->{_m} = Math::BigInt->bone();
196 $self->{_e} = Math::BigInt->bzero();
197 $self->{sign} = $sign;
394e6ffb 198 ($self->{_a},$self->{_p}) = @_ if @_ > 0;
574bacfe 199 return $self;
200 }
201
58cde26e 202sub bzero
203 {
204 # create a bigfloat '+0', if given a BigFloat, set it to 0
205 my $self = shift;
206 $self = $class if !defined $self;
207 if (!ref($self))
208 {
209 my $c = $self; $self = {}; bless $self, $c;
210 }
574bacfe 211 $self->{_m} = Math::BigInt->bzero();
212 $self->{_e} = Math::BigInt->bone();
58cde26e 213 $self->{sign} = '+';
394e6ffb 214 ($self->{_a},$self->{_p}) = @_ if @_ > 0;
58cde26e 215 return $self;
216 }
217
218##############################################################################
219# string conversation
220
221sub bstr
222 {
223 # (ref to BFLOAT or num_str ) return num_str
224 # Convert number from internal format to (non-scientific) string format.
225 # internal format is always normalized (no leading zeros, "-0" => "+0")
ee15d750 226 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
227 #my $x = shift; my $class = ref($x) || $x;
228 #$x = $class->new(shift) unless ref($x);
58cde26e 229
574bacfe 230 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
231 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
232 if ($x->{sign} !~ /^[+-]$/)
58cde26e 233 {
574bacfe 234 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
235 return 'inf'; # +inf
58cde26e 236 }
237
574bacfe 238 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
239
240 my $not_zero = !$x->is_zero();
241 if ($not_zero)
58cde26e 242 {
574bacfe 243 $es = $x->{_m}->bstr();
244 $len = CORE::length($es);
245 if (!$x->{_e}->is_zero())
246# {
247# $es = $x->{sign}.$es if $x->{sign} eq '-';
248# }
249# else
58cde26e 250 {
574bacfe 251 if ($x->{_e}->sign() eq '-')
252 {
253 $dot = '';
254 if ($x->{_e} <= -$len)
255 {
256 # print "style: 0.xxxx\n";
257 my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
258 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
259 }
260 else
261 {
262 # print "insert '.' at $x->{_e} in '$es'\n";
263 substr($es,$x->{_e},0) = '.'; $cad = $x->{_e};
264 }
265 }
266 else
267 {
268 # expand with zeros
269 $es .= '0' x $x->{_e}; $len += $x->{_e}; $cad = 0;
270 }
82cf049f 271 }
574bacfe 272 } # if not zero
273 $es = $x->{sign}.$es if $x->{sign} eq '-';
274 # if set accuracy or precision, pad with zeros
275 if ((defined $x->{_a}) && ($not_zero))
276 {
277 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
278 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
279 $zeros = $x->{_a} - $len if $cad != $len;
280 #print "acc padd $x->{_a} $zeros (len $len cad $cad)\n";
281 $es .= $dot.'0' x $zeros if $zeros > 0;
82cf049f 282 }
574bacfe 283 elsif ($x->{_p} || 0 < 0)
58cde26e 284 {
574bacfe 285 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
286 my $zeros = -$x->{_p} + $cad;
287 #print "pre padd $x->{_p} $zeros (len $len cad $cad)\n";
288 $es .= $dot.'0' x $zeros if $zeros > 0;
58cde26e 289 }
58cde26e 290 return $es;
82cf049f 291 }
f216259d 292
58cde26e 293sub bsstr
294 {
295 # (ref to BFLOAT or num_str ) return num_str
296 # Convert number from internal format to scientific string format.
297 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
ee15d750 298 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
299 #my $x = shift; my $class = ref($x) || $x;
300 #$x = $class->new(shift) unless ref($x);
a0d0e21e 301
574bacfe 302 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
303 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
304 if ($x->{sign} !~ /^[+-]$/)
305 {
306 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
307 return 'inf'; # +inf
308 }
58cde26e 309 my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
310 my $sep = 'e'.$sign;
311 return $x->{_m}->bstr().$sep.$x->{_e}->bstr();
312 }
313
314sub numify
315 {
316 # Make a number from a BigFloat object
574bacfe 317 # simple return string and let Perl's atoi()/atof() handle the rest
ee15d750 318 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 319 return $x->bsstr();
320 }
a0d0e21e 321
58cde26e 322##############################################################################
323# public stuff (usually prefixed with "b")
324
574bacfe 325# tels 2001-08-04
326# todo: this must be overwritten and return NaN for non-integer values
327# band(), bior(), bxor(), too
58cde26e 328#sub bnot
329# {
330# $class->SUPER::bnot($class,@_);
331# }
332
333sub bcmp
334 {
335 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
336 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
337 my ($self,$x,$y) = objectify(2,@_);
58cde26e 338
0716bf9b 339 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
340 {
341 # handle +-inf and NaN
342 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
343 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
344 return +1 if $x->{sign} eq '+inf';
345 return -1 if $x->{sign} eq '-inf';
346 return -1 if $y->{sign} eq '+inf';
347 return +1 if $y->{sign} eq '-inf';
348 }
349
350 # check sign for speed first
574bacfe 351 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
58cde26e 352 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
353
574bacfe 354 # shortcut
355 my $xz = $x->is_zero();
356 my $yz = $y->is_zero();
357 return 0 if $xz && $yz; # 0 <=> 0
358 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
359 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
58cde26e 360
361 # adjust so that exponents are equal
bd05a461 362 my $lxm = $x->{_m}->length();
363 my $lym = $y->{_m}->length();
364 my $lx = $lxm + $x->{_e};
365 my $ly = $lym + $y->{_e};
58cde26e 366 # print "x $x y $y lx $lx ly $ly\n";
367 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
368 # print "$l $x->{sign}\n";
bd05a461 369 return $l <=> 0 if $l != 0;
58cde26e 370
bd05a461 371 # lengths (corrected by exponent) are equal
372 # so make mantissa euqal length by padding with zero (shift left)
373 my $diff = $lxm - $lym;
374 my $xm = $x->{_m}; # not yet copy it
375 my $ym = $y->{_m};
376 if ($diff > 0)
377 {
378 $ym = $y->{_m}->copy()->blsft($diff,10);
379 }
380 elsif ($diff < 0)
381 {
382 $xm = $x->{_m}->copy()->blsft(-$diff,10);
383 }
384 my $rc = $xm->bcmp($ym);
58cde26e 385 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
bd05a461 386 return $rc <=> 0;
58cde26e 387 }
388
389sub bacmp
390 {
391 # Compares 2 values, ignoring their signs.
392 # Returns one of undef, <0, =0, >0. (suitable for sort)
393 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
394 my ($self,$x,$y) = objectify(2,@_);
ee15d750 395
396 # handle +-inf and NaN's
abcfbf51 397 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
ee15d750 398 {
399 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
400 return 0 if ($x->is_inf() && $y->is_inf());
401 return 1 if ($x->is_inf() && !$y->is_inf());
402 return -1 if (!$x->is_inf() && $y->is_inf());
403 }
404
405 # shortcut
406 my $xz = $x->is_zero();
407 my $yz = $y->is_zero();
408 return 0 if $xz && $yz; # 0 <=> 0
409 return -1 if $xz && !$yz; # 0 <=> +y
410 return 1 if $yz && !$xz; # +x <=> 0
411
412 # adjust so that exponents are equal
413 my $lxm = $x->{_m}->length();
414 my $lym = $y->{_m}->length();
415 my $lx = $lxm + $x->{_e};
416 my $ly = $lym + $y->{_e};
417 # print "x $x y $y lx $lx ly $ly\n";
394e6ffb 418 my $l = $lx - $ly;
ee15d750 419 # print "$l $x->{sign}\n";
420 return $l <=> 0 if $l != 0;
58cde26e 421
ee15d750 422 # lengths (corrected by exponent) are equal
394e6ffb 423 # so make mantissa equal-length by padding with zero (shift left)
ee15d750 424 my $diff = $lxm - $lym;
425 my $xm = $x->{_m}; # not yet copy it
426 my $ym = $y->{_m};
427 if ($diff > 0)
428 {
429 $ym = $y->{_m}->copy()->blsft($diff,10);
430 }
431 elsif ($diff < 0)
432 {
433 $xm = $x->{_m}->copy()->blsft(-$diff,10);
434 }
435 my $rc = $xm->bcmp($ym);
ee15d750 436 return $rc <=> 0;
58cde26e 437 }
a0d0e21e 438
58cde26e 439sub badd
440 {
441 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
442 # return result as BFLOAT
58cde26e 443 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
444
574bacfe 445 # inf and NaN handling
446 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
447 {
448 # NaN first
449 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
450 # inf handline
451 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
452 {
453 # + and + => +, - and - => -, + and - => 0, - and + => 0
454 return $x->bzero() if $x->{sign} ne $y->{sign};
455 return $x;
456 }
457 # +-inf + something => +inf
458 # something +-inf => +-inf
459 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
460 return $x;
461 }
462
58cde26e 463 # speed: no add for 0+y or x+0
464 return $x if $y->is_zero(); # x+0
465 if ($x->is_zero()) # 0+y
466 {
467 # make copy, clobbering up x (modify in place!)
468 $x->{_e} = $y->{_e}->copy();
469 $x->{_m} = $y->{_m}->copy();
470 $x->{sign} = $y->{sign} || $nan;
471 return $x->round($a,$p,$r,$y);
a0d0e21e 472 }
58cde26e 473
474 # take lower of the two e's and adapt m1 to it to match m2
475 my $e = $y->{_e}; $e = Math::BigInt::bzero() if !defined $e; # if no BFLOAT
476 $e = $e - $x->{_e};
477 my $add = $y->{_m}->copy();
478 if ($e < 0)
479 {
0716bf9b 480 # print "e < 0\n";
58cde26e 481 #print "\$x->{_m}: $x->{_m} ";
482 #print "\$x->{_e}: $x->{_e}\n";
483 my $e1 = $e->copy()->babs();
484 $x->{_m} *= (10 ** $e1);
485 $x->{_e} += $e; # need the sign of e
486 #$x->{_m} += $y->{_m};
487 #print "\$x->{_m}: $x->{_m} ";
488 #print "\$x->{_e}: $x->{_e}\n";
489 }
490 elsif ($e > 0)
491 {
0716bf9b 492 # print "e > 0\n";
58cde26e 493 #print "\$x->{_m}: $x->{_m} \$y->{_m}: $y->{_m} \$e: $e ",ref($e),"\n";
494 $add *= (10 ** $e);
495 #$x->{_m} += $y->{_m} * (10 ** $e);
496 #print "\$x->{_m}: $x->{_m}\n";
497 }
498 # else: both e are same, so leave them
499 #print "badd $x->{sign}$x->{_m} + $y->{sign}$add\n";
500 # fiddle with signs
501 $x->{_m}->{sign} = $x->{sign};
502 $add->{sign} = $y->{sign};
503 # finally do add/sub
504 $x->{_m} += $add;
505 # re-adjust signs
506 $x->{sign} = $x->{_m}->{sign};
507 $x->{_m}->{sign} = '+';
0716bf9b 508 #$x->bnorm(); # delete trailing zeros
58cde26e 509 return $x->round($a,$p,$r,$y);
510 }
511
512sub bsub
513 {
0716bf9b 514 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
58cde26e 515 # subtract second arg from first, modify first
e745a66c 516 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
a0d0e21e 517
e745a66c 518 if (!$y->is_zero()) # don't need to do anything if $y is 0
519 {
520 $y->{sign} =~ tr/+\-/-+/; # does nothing for NaN
521 $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
522 $y->{sign} =~ tr/+\-/-+/; # refix $y (does nothing for NaN)
523 }
524 $x; # already rounded by badd()
58cde26e 525 }
526
527sub binc
528 {
529 # increment arg by one
ee15d750 530 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
e745a66c 531
532 if ($x->{_e}->sign() eq '-')
533 {
534 return $x->badd($self->bone(),$a,$p,$r); # digits after dot
535 }
536
537 if (!$x->{_e}->is_zero())
538 {
539 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
540 $x->{_e}->bzero();
541 }
542 # now $x->{_e} == 0
543 if ($x->{sign} eq '+')
544 {
545 $x->{_m}->binc();
546 return $x->bnorm()->bround($a,$p,$r);
547 }
548 elsif ($x->{sign} eq '-')
549 {
550 $x->{_m}->bdec();
551 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
552 return $x->bnorm()->bround($a,$p,$r);
553 }
554 # inf, nan handling etc
555 $x->badd($self->__one(),$a,$p,$r); # does round
58cde26e 556 }
557
558sub bdec
559 {
560 # decrement arg by one
ee15d750 561 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
e745a66c 562
563 if ($x->{_e}->sign() eq '-')
564 {
565 return $x->badd($self->bone('-'),$a,$p,$r); # digits after dot
566 }
567
568 if (!$x->{_e}->is_zero())
569 {
570 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
571 $x->{_e}->bzero();
572 }
573 # now $x->{_e} == 0
574 my $zero = $x->is_zero();
575 # <= 0
576 if (($x->{sign} eq '-') || $zero)
577 {
578 $x->{_m}->binc();
579 $x->{sign} = '-' if $zero; # 0 => 1 => -1
580 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
581 return $x->bnorm()->round($a,$p,$r);
582 }
583 # > 0
584 elsif ($x->{sign} eq '+')
585 {
586 $x->{_m}->bdec();
587 return $x->bnorm()->round($a,$p,$r);
588 }
589 # inf, nan handling etc
590 $x->badd($self->bone('-'),$a,$p,$r); # does round
58cde26e 591 }
592
593sub blcm
594 {
ee15d750 595 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
58cde26e 596 # does not modify arguments, but returns new object
597 # Lowest Common Multiplicator
58cde26e 598
599 my ($self,@arg) = objectify(0,@_);
600 my $x = $self->new(shift @arg);
601 while (@arg) { $x = _lcm($x,shift @arg); }
602 $x;
603 }
604
605sub bgcd
606 {
ee15d750 607 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
58cde26e 608 # does not modify arguments, but returns new object
609 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
58cde26e 610
611 my ($self,@arg) = objectify(0,@_);
612 my $x = $self->new(shift @arg);
613 while (@arg) { $x = _gcd($x,shift @arg); }
614 $x;
615 }
616
617sub is_zero
618 {
ee15d750 619 # return true if arg (BFLOAT or num_str) is zero (array '+', '0')
620 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
574bacfe 621
622 return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
623 return 0;
58cde26e 624 }
625
626sub is_one
627 {
ee15d750 628 # return true if arg (BFLOAT or num_str) is +1 (array '+', '1')
58cde26e 629 # or -1 if signis given
ee15d750 630 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
631
632 my $sign = shift || ''; $sign = '+' if $sign ne '-';
633 return 1
634 if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
635 return 0;
58cde26e 636 }
637
638sub is_odd
639 {
ee15d750 640 # return true if arg (BFLOAT or num_str) is odd or false if even
641 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
0716bf9b 642
643 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
ee15d750 644 return 1 if ($x->{_e}->is_zero() && $x->{_m}->is_odd());
645 return 0;
58cde26e 646 }
647
648sub is_even
649 {
b22b3e31 650 # return true if arg (BINT or num_str) is even or false if odd
ee15d750 651 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
0716bf9b 652
653 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
654 return 1 if $x->{_m}->is_zero(); # 0e1 is even
ee15d750 655 return 1 if ($x->{_e}->is_zero() && $x->{_m}->is_even()); # 123.45 is never
656 return 0;
58cde26e 657 }
658
659sub bmul
660 {
661 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
662 # (BINT or num_str, BINT or num_str) return BINT
663 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
58cde26e 664
0716bf9b 665 # print "mbf bmul $x->{_m}e$x->{_e} $y->{_m}e$y->{_e}\n";
58cde26e 666 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
667
574bacfe 668 # handle result = 0
669 return $x->bzero() if $x->is_zero() || $y->is_zero();
670 # inf handling
671 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
672 {
673 # result will always be +-inf:
674 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
675 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
676 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
677 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
678 return $x->binf('-');
679 }
680
58cde26e 681 # aEb * cEd = (a*c)E(b+d)
394e6ffb 682 $x->{_m}->bmul($y->{_m});
683 $x->{_e}->badd($y->{_e});
58cde26e 684 # adjust sign:
685 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
394e6ffb 686 return $x->bnorm()->round($a,$p,$r,$y);
58cde26e 687 }
688
689sub bdiv
690 {
691 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
692 # (BFLOAT,BFLOAT) (quo,rem) or BINT (only rem)
693 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
694
ee15d750 695
574bacfe 696 # x / +-inf => 0, reminder x
697 return wantarray ? ($x->bzero(),$x->copy()) : $x->bzero()
698 if $y->{sign} =~ /^[+-]inf$/;
699
700 # NaN if x == NaN or y == NaN or x==y==0
58cde26e 701 return wantarray ? ($x->bnan(),bnan()) : $x->bnan()
574bacfe 702 if (($x->is_nan() || $y->is_nan()) ||
703 ($x->is_zero() && $y->is_zero()));
704
705 # 5 / 0 => +inf, -6 / 0 => -inf
706 return wantarray
707 ? ($x->binf($x->{sign}),$self->bnan()) : $x->binf($x->{sign})
708 if ($x->{sign} =~ /^[+-]$/ && $y->is_zero());
0716bf9b 709
394e6ffb 710 # x== 0 or y == 1 or y == -1
711 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
0716bf9b 712
58cde26e 713 # we need to limit the accuracy to protect against overflow
574bacfe 714 my $fallback = 0;
ee15d750 715 my $scale = 0;
ee15d750 716 my @params = $x->_find_round_parameters($a,$p,$r,$y);
717
718 # no rounding at all, so must use fallback
719 if (scalar @params == 1)
58cde26e 720 {
0716bf9b 721 # simulate old behaviour
ee15d750 722 $scale = $self->div_scale()+1; # at least one more for proper round
723 $params[1] = $self->div_scale(); # and round to it as accuracy
724 $params[3] = $r; # round mode by caller or undef
725 $fallback = 1; # to clear a/p afterwards
726 }
727 else
728 {
729 # the 4 below is empirical, and there might be cases where it is not
730 # enough...
731 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
a0d0e21e 732 }
0716bf9b 733 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
58cde26e 734 $scale = $lx if $lx > $scale;
58cde26e 735 $scale = $ly if $ly > $scale;
0716bf9b 736 my $diff = $ly - $lx;
737 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
a0d0e21e 738
58cde26e 739 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
a0d0e21e 740
58cde26e 741 # check for / +-1 ( +/- 1E0)
394e6ffb 742 if (!$y->is_one())
58cde26e 743 {
394e6ffb 744 # promote BigInts and it's subclasses (except when already a BigFloat)
745 $y = $self->new($y) unless $y->isa('Math::BigFloat');
746
747 # calculate the result to $scale digits and then round it
748 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
749 $x->{_m}->blsft($scale,10);
750 $x->{_m}->bdiv( $y->{_m} ); # a/c
751 $x->{_e}->bsub( $y->{_e} ); # b-d
752 $x->{_e}->bsub($scale); # correct for 10**scale
753 $x->bnorm(); # remove trailing 0's
a0d0e21e 754 }
a5f75d66 755
ee15d750 756 # shortcut to not run trough _find_round_parameters again
757 if (defined $params[1])
758 {
759 $x->bround($params[1],undef,$params[3]); # then round accordingly
760 }
761 else
762 {
763 $x->bfround($params[2],$params[3]); # then round accordingly
764 }
574bacfe 765 if ($fallback)
766 {
767 # clear a/p after round, since user did not request it
ee15d750 768 $x->{_a} = undef; $x->{_p} = undef;
574bacfe 769 }
0716bf9b 770
58cde26e 771 if (wantarray)
772 {
394e6ffb 773 my $rem;
774 if (!$y->is_one())
775 {
776 $rem = $x->copy();
777 $rem->bmod($y,$params[1],$params[2],$params[3]);
778 }
779 else
780 {
781 $rem = $self->bzero();
782 }
574bacfe 783 if ($fallback)
784 {
785 # clear a/p after round, since user did not request it
ee15d750 786 $rem->{_a} = undef; $rem->{_p} = undef;
574bacfe 787 }
0716bf9b 788 return ($x,$rem);
58cde26e 789 }
790 return $x;
791 }
a0d0e21e 792
58cde26e 793sub bmod
794 {
795 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
796 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
a0d0e21e 797
abcfbf51 798 return $x->bnan() if ($x->is_nan() || $y->is_nan() || $y->is_zero());
58cde26e 799 return $x->bzero() if $y->is_one();
800
801 # XXX tels: not done yet
802 return $x->round($a,$p,$r,$y);
803 }
804
805sub bsqrt
806 {
0716bf9b 807 # calculate square root; this should probably
808 # use a different test to see whether the accuracy we want is...
ee15d750 809 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 810
0716bf9b 811 return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
812 return $x if $x->{sign} eq '+inf'; # +inf
394e6ffb 813 return $x if $x->is_zero() || $x->is_one();
58cde26e 814
ee15d750 815 # we need to limit the accuracy to protect against overflow (ignore $p)
816 my ($scale) = $x->_scale_a($self->accuracy(),$self->round_mode,$a,$r);
574bacfe 817 my $fallback = 0;
0716bf9b 818 if (!defined $scale)
819 {
820 # simulate old behaviour
ee15d750 821 $scale = $self->div_scale()+1; # one more for proper riund
822 $a = $self->div_scale(); # and round to it
823 $fallback = 1; # to clear a/p afterwards
0716bf9b 824 }
394e6ffb 825 my $xas = $x->as_number();
826 my $gs = $xas->copy()->bsqrt(); # some guess
827 if (($x->{_e}->{sign} ne '-') # guess can't be accurate if there are
828 # digits after the dot
829 && ($xas->bcmp($gs * $gs) == 0)) # guess hit the nail on the head?
830 {
831 # exact result
832 $x->{_m} = $gs;
833 # leave alone if _e is already right
834 $x->{_e} = Math::BigInt->bzero();
835 return $x->bnorm()->round($a,$p,$r)
836 }
837 $gs = $self->new( $gs );
838
0716bf9b 839 my $lx = $x->{_m}->length();
840 $scale = $lx if $scale < $lx;
394e6ffb 841 my $e = $self->new("1E-$scale"); # make test variable
58cde26e 842 return $x->bnan() if $e->sign() eq 'NaN';
843
58cde26e 844 # start with some reasonable guess
394e6ffb 845# $lx = $lx+$x->{_e};
846# $lx = $lx / 2;
847# $lx = 1 if $lx < 1;
848 # my $gs = Math::BigFloat->new("1E$lx");
849
850# print "first guess: $gs (x $x) scale $scale\n";
851# # use BigInt:sqrt as reasonabe guess
852# print "second guess: $gs (x $x) scale $scale\n";
853
58cde26e 854 my $diff = $e;
855 my $y = $x->copy();
394e6ffb 856 my $two = $self->new(2);
ee15d750 857 # promote BigInts and it's subclasses (except when already a BigFloat)
858 $y = $self->new($y) unless $y->isa('Math::BigFloat');
ee15d750 859 my $rem;
394e6ffb 860# my $steps = 0;
58cde26e 861 while ($diff >= $e)
862 {
394e6ffb 863 # return $x->bnan() if $gs->is_zero();
864
865 $x = $y->copy()->bdiv($gs,$scale)->badd($gs)->bdiv($two,$scale);
58cde26e 866 $diff = $x->copy()->bsub($gs)->babs();
58cde26e 867 $gs = $x->copy();
394e6ffb 868# $steps++;
a0d0e21e 869 }
394e6ffb 870# print "steps $steps\n";
0716bf9b 871 $x->round($a,$p,$r);
574bacfe 872 if ($fallback)
873 {
874 # clear a/p after round, since user did not request it
ee15d750 875 $x->{_a} = undef; $x->{_p} = undef;
574bacfe 876 }
877 $x;
58cde26e 878 }
879
880sub bpow
881 {
882 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
883 # compute power of two numbers, second arg is used as integer
884 # modifies first argument
885
886 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
887
0716bf9b 888 return $x if $x->{sign} =~ /^[+-]inf$/;
58cde26e 889 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
574bacfe 890 return $x->bone() if $y->is_zero();
58cde26e 891 return $x if $x->is_one() || $y->is_one();
ee15d750 892 my $y1 = $y->as_number(); # make bigint (trunc)
394e6ffb 893 # if ($x == -1)
894 if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
58cde26e 895 {
896 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
0716bf9b 897 return $y1->is_odd() ? $x : $x->babs(1);
288d023a 898 }
58cde26e 899 return $x if $x->is_zero() && $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
574bacfe 900 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
901 return $x->binf() if $x->is_zero() && $y->{sign} eq '-';
58cde26e 902
903 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
904 $y1->babs();
905 $x->{_m}->bpow($y1);
906 $x->{_e}->bmul($y1);
907 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
908 $x->bnorm();
909 if ($y->{sign} eq '-')
910 {
911 # modify $x in place!
0716bf9b 912 my $z = $x->copy(); $x->bzero()->binc();
58cde26e 913 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
a0d0e21e 914 }
58cde26e 915 return $x->round($a,$p,$r,$y);
916 }
917
918###############################################################################
919# rounding functions
920
921sub bfround
922 {
923 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
924 # $n == 0 means round to integer
925 # expects and returns normalized numbers!
ee15d750 926 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
a0d0e21e 927
58cde26e 928 return $x if $x->modify('bfround');
929
ee15d750 930 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
58cde26e 931 return $x if !defined $scale; # no-op
932
574bacfe 933 # never round a 0, +-inf, NaN
934 return $x if $x->{sign} !~ /^[+-]$/ || $x->is_zero();
58cde26e 935 # print "MBF bfround $x to scale $scale mode $mode\n";
58cde26e 936
ee15d750 937 # don't round if x already has lower precision
938 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
939
940 $x->{_p} = $scale; # remember round in any case
941 $x->{_a} = undef; # and clear A
58cde26e 942 if ($scale < 0)
943 {
944 # print "bfround scale $scale e $x->{_e}\n";
945 # round right from the '.'
946 return $x if $x->{_e} >= 0; # nothing to round
947 $scale = -$scale; # positive for simplicity
948 my $len = $x->{_m}->length(); # length of mantissa
949 my $dad = -$x->{_e}; # digits after dot
950 my $zad = 0; # zeros after dot
951 $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
ee15d750 952 #print "scale $scale dad $dad zad $zad len $len\n";
58cde26e 953
954 # number bsstr len zad dad
955 # 0.123 123e-3 3 0 3
956 # 0.0123 123e-4 3 1 4
957 # 0.001 1e-3 1 2 3
958 # 1.23 123e-2 3 0 2
959 # 1.2345 12345e-4 5 0 4
960
961 # do not round after/right of the $dad
962 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
963
ee15d750 964 # round to zero if rounding inside the $zad, but not for last zero like:
965 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
966 return $x->bzero() if $scale < $zad;
967 if ($scale == $zad) # for 0.006, scale -3 and trunc
58cde26e 968 {
ee15d750 969 $scale = -$len-1;
58cde26e 970 }
971 else
972 {
973 # adjust round-point to be inside mantissa
974 if ($zad != 0)
975 {
976 $scale = $scale-$zad;
977 }
978 else
979 {
980 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
981 $scale = $dbd+$scale;
982 }
983 }
984 # print "round to $x->{_m} to $scale\n";
a0d0e21e 985 }
58cde26e 986 else
987 {
988 # 123 => 100 means length(123) = 3 - $scale (2) => 1
a5f75d66 989
58cde26e 990 # calculate digits before dot
991 my $dbt = $x->{_m}->length(); $dbt += $x->{_e} if $x->{_e}->sign() eq '-';
ee15d750 992 # if not enough digits before dot, round to zero
993 return $x->bzero() if ($scale > $dbt) && ($dbt < 0);
994 # scale always >= 0 here
995 if ($dbt == 0)
58cde26e 996 {
997 # 0.49->bfround(1): scale == 1, dbt == 0: => 0.0
998 # 0.51->bfround(0): scale == 0, dbt == 0: => 1.0
999 # 0.5->bfround(0): scale == 0, dbt == 0: => 0
1000 # 0.05->bfround(0): scale == 0, dbt == 0: => 0
1001 # print "$scale $dbt $x->{_m}\n";
1002 $scale = -$x->{_m}->length();
1003 }
1004 elsif ($dbt > 0)
1005 {
1006 # correct by subtracting scale
1007 $scale = $dbt - $scale;
1008 }
1009 else
1010 {
1011 $scale = $x->{_m}->length() - $scale;
1012 }
a0d0e21e 1013 }
574bacfe 1014 # print "using $scale for $x->{_m} with '$mode'\n";
1015 # pass sign to bround for rounding modes '+inf' and '-inf'
58cde26e 1016 $x->{_m}->{sign} = $x->{sign};
1017 $x->{_m}->bround($scale,$mode);
1018 $x->{_m}->{sign} = '+'; # fix sign back
1019 $x->bnorm();
1020 }
1021
1022sub bround
1023 {
1024 # accuracy: preserve $N digits, and overwrite the rest with 0's
ee15d750 1025 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1026
1027 die ('bround() needs positive accuracy') if ($_[0] || 0) < 0;
58cde26e 1028
ee15d750 1029 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1030 return $x if !defined $scale; # no-op
1031
58cde26e 1032 return $x if $x->modify('bround');
ee15d750 1033
1034 # scale is now either $x->{_a}, $accuracy, or the user parameter
1035 # test whether $x already has lower accuracy, do nothing in this case
1036 # but do round if the accuracy is the same, since a math operation might
1037 # want to round a number with A=5 to 5 digits afterwards again
1038 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
58cde26e 1039
1040 # print "bround $scale $mode\n";
1041 # 0 => return all digits, scale < 0 makes no sense
1042 return $x if ($scale <= 0);
574bacfe 1043 # never round a 0, +-inf, NaN
1044 return $x if $x->{sign} !~ /^[+-]$/ || $x->is_zero();
58cde26e 1045
1046 # if $e longer than $m, we have 0.0000xxxyyy style number, and must
1047 # subtract the delta from scale, to simulate keeping the zeros
1048 # -5 +5 => 1; -10 +5 => -4
1049 my $delta = $x->{_e} + $x->{_m}->length() + 1;
58cde26e 1050
1051 # if we should keep more digits than the mantissa has, do nothing
1052 return $x if $x->{_m}->length() <= $scale;
f216259d 1053
58cde26e 1054 # pass sign to bround for '+inf' and '-inf' rounding modes
1055 $x->{_m}->{sign} = $x->{sign};
1056 $x->{_m}->bround($scale,$mode); # round mantissa
1057 $x->{_m}->{sign} = '+'; # fix sign back
ee15d750 1058 $x->{_a} = $scale; # remember rounding
1059 $x->{_p} = undef; # and clear P
574bacfe 1060 $x->bnorm(); # del trailing zeros gen. by bround()
58cde26e 1061 }
1062
1063sub bfloor
1064 {
1065 # return integer less or equal then $x
ee15d750 1066 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 1067
1068 return $x if $x->modify('bfloor');
1069
1070 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1071
1072 # if $x has digits after dot
1073 if ($x->{_e}->{sign} eq '-')
1074 {
1075 $x->{_m}->brsft(-$x->{_e},10);
1076 $x->{_e}->bzero();
1077 $x-- if $x->{sign} eq '-';
f216259d 1078 }
58cde26e 1079 return $x->round($a,$p,$r);
1080 }
288d023a 1081
58cde26e 1082sub bceil
1083 {
1084 # return integer greater or equal then $x
ee15d750 1085 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 1086
1087 return $x if $x->modify('bceil');
1088 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1089
1090 # if $x has digits after dot
1091 if ($x->{_e}->{sign} eq '-')
1092 {
1093 $x->{_m}->brsft(-$x->{_e},10);
1094 $x->{_e}->bzero();
1095 $x++ if $x->{sign} eq '+';
a0d0e21e 1096 }
58cde26e 1097 return $x->round($a,$p,$r);
1098 }
1099
394e6ffb 1100sub brsft
1101 {
1102 # shift right by $y (divide by power of 2)
1103 my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1104
1105 return $x if $x->modify('brsft');
1106 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1107
1108 $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1109 $x->bdiv($n ** $y,$a,$p,$r,$y);
1110 }
1111
1112sub blsft
1113 {
1114 # shift right by $y (divide by power of 2)
1115 my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1116
1117 return $x if $x->modify('brsft');
1118 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1119
1120 $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1121 $x->bmul($n ** $y,$a,$p,$r,$y);
1122 }
1123
58cde26e 1124###############################################################################
a5f75d66 1125
58cde26e 1126sub DESTROY
1127 {
ee15d750 1128 # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
58cde26e 1129 }
1130
1131sub AUTOLOAD
1132 {
1133 # make fxxx and bxxx work
1134 # my $self = $_[0];
1135 my $name = $AUTOLOAD;
1136
1137 $name =~ s/.*:://; # split package
1138 #print "$name\n";
ee15d750 1139 no strict 'refs';
1140 if (!method_alias($name))
58cde26e 1141 {
ee15d750 1142 if (!defined $name)
1143 {
1144 # delayed load of Carp and avoid recursion
1145 require Carp;
1146 Carp::croak ("Can't call a method without name");
1147 }
ee15d750 1148 if (!method_hand_up($name))
1149 {
1150 # delayed load of Carp and avoid recursion
1151 require Carp;
1152 Carp::croak ("Can't call $class\-\>$name, not a valid method");
1153 }
1154 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1155 $name =~ s/^f/b/;
1156 return &{'Math::BigInt'."::$name"}(@_);
a0d0e21e 1157 }
58cde26e 1158 my $bname = $name; $bname =~ s/^f/b/;
1159 *{$class."\:\:$name"} = \&$bname;
1160 &$bname; # uses @_
1161 }
1162
1163sub exponent
1164 {
1165 # return a copy of the exponent
ee15d750 1166 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 1167
ee15d750 1168 if ($x->{sign} !~ /^[+-]$/)
1169 {
1170 my $s = $x->{sign}; $s =~ s/^[+-]//;
1171 return $self->new($s); # -inf, +inf => +inf
1172 }
1173 return $x->{_e}->copy();
58cde26e 1174 }
1175
1176sub mantissa
1177 {
1178 # return a copy of the mantissa
ee15d750 1179 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 1180
ee15d750 1181 if ($x->{sign} !~ /^[+-]$/)
1182 {
1183 my $s = $x->{sign}; $s =~ s/^[+]//;
1184 return $self->new($s); # -inf, +inf => +inf
1185 }
1186 my $m = $x->{_m}->copy(); # faster than going via bstr()
1187 $m->bneg() if $x->{sign} eq '-';
58cde26e 1188
1189 return $m;
1190 }
1191
1192sub parts
1193 {
1194 # return a copy of both the exponent and the mantissa
ee15d750 1195 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 1196
ee15d750 1197 if ($x->{sign} !~ /^[+-]$/)
1198 {
1199 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
1200 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
1201 }
1202 my $m = $x->{_m}->copy(); # faster than going via bstr()
1203 $m->bneg() if $x->{sign} eq '-';
1204 return ($m,$x->{_e}->copy());
58cde26e 1205 }
1206
1207##############################################################################
1208# private stuff (internal use only)
1209
58cde26e 1210sub import
1211 {
1212 my $self = shift;
58cde26e 1213 for ( my $i = 0; $i < @_ ; $i++ )
1214 {
1215 if ( $_[$i] eq ':constant' )
1216 {
1217 # this rest causes overlord er load to step in
1218 # print "overload @_\n";
1219 overload::constant float => sub { $self->new(shift); };
1220 splice @_, $i, 1; last;
1221 }
1222 }
1223 # any non :constant stuff is handled by our parent, Exporter
1224 # even if @_ is empty, to give it a chance
abcfbf51 1225 $self->SUPER::import(@_); # for subclasses
1226 $self->export_to_level(1,$self,@_); # need this, too
58cde26e 1227 }
1228
1229sub bnorm
1230 {
1231 # adjust m and e so that m is smallest possible
1232 # round number according to accuracy and precision settings
ee15d750 1233 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 1234
0716bf9b 1235 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
58cde26e 1236
1237 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
1238 if ($zeros != 0)
1239 {
1240 $x->{_m}->brsft($zeros,10); $x->{_e} += $zeros;
1241 }
ee15d750 1242 # for something like 0Ey, set y to 1, and -0 => +0
1243 $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
1244 # this is to prevent automatically rounding when MBI's globals are set
0716bf9b 1245 $x->{_m}->{_f} = MB_NEVER_ROUND;
1246 $x->{_e}->{_f} = MB_NEVER_ROUND;
ee15d750 1247 # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
1248 $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
1249 $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
394e6ffb 1250 return $x; # MBI bnorm is no-op, so dont call it
58cde26e 1251 }
1252
1253##############################################################################
1254# internal calculation routines
1255
1256sub as_number
1257 {
394e6ffb 1258 # return copy as a bigint representation of this BigFloat number
1259 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 1260
1261 my $z;
1262 if ($x->{_e}->is_zero())
1263 {
1264 $z = $x->{_m}->copy();
1265 $z->{sign} = $x->{sign};
1266 return $z;
1267 }
0716bf9b 1268 $z = $x->{_m}->copy();
58cde26e 1269 if ($x->{_e} < 0)
1270 {
0716bf9b 1271 $z->brsft(-$x->{_e},10);
1272 }
1273 else
1274 {
1275 $z->blsft($x->{_e},10);
58cde26e 1276 }
58cde26e 1277 $z->{sign} = $x->{sign};
1278 return $z;
1279 }
1280
1281sub length
1282 {
ee15d750 1283 my $x = shift;
1284 my $class = ref($x) || $x;
1285 $x = $class->new(shift) unless ref($x);
58cde26e 1286
ee15d750 1287 return 1 if $x->{_m}->is_zero();
58cde26e 1288 my $len = $x->{_m}->length();
1289 $len += $x->{_e} if $x->{_e}->sign() eq '+';
1290 if (wantarray())
1291 {
1292 my $t = Math::BigInt::bzero();
1293 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1294 return ($len,$t);
1295 }
1296 return $len;
1297 }
a0d0e21e 1298
12991;
a5f75d66 1300__END__
1301
1302=head1 NAME
1303
58cde26e 1304Math::BigFloat - Arbitrary size floating point math package
a5f75d66 1305
1306=head1 SYNOPSIS
1307
a2008d6d 1308 use Math::BigFloat;
58cde26e 1309
1310 # Number creation
1311 $x = Math::BigInt->new($str); # defaults to 0
1312 $nan = Math::BigInt->bnan(); # create a NotANumber
1313 $zero = Math::BigInt->bzero();# create a "+0"
1314
1315 # Testing
1316 $x->is_zero(); # return whether arg is zero or not
0716bf9b 1317 $x->is_nan(); # return whether arg is NaN or not
1318 $x->is_one(); # true if arg is +1
1319 $x->is_one('-'); # true if arg is -1
1320 $x->is_odd(); # true if odd, false for even
1321 $x->is_even(); # true if even, false for odd
1322 $x->is_positive(); # true if >= 0
1323 $x->is_negative(); # true if < 0
1324 $x->is_inf(sign) # true if +inf or -inf (sign default '+')
58cde26e 1325 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
1326 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
1327 $x->sign(); # return the sign, either +,- or NaN
1328
1329 # The following all modify their first argument:
1330
1331 # set
1332 $x->bzero(); # set $i to 0
1333 $x->bnan(); # set $i to NaN
1334
1335 $x->bneg(); # negation
1336 $x->babs(); # absolute value
1337 $x->bnorm(); # normalize (no-op)
1338 $x->bnot(); # two's complement (bit wise not)
1339 $x->binc(); # increment x by 1
1340 $x->bdec(); # decrement x by 1
1341
1342 $x->badd($y); # addition (add $y to $x)
1343 $x->bsub($y); # subtraction (subtract $y from $x)
1344 $x->bmul($y); # multiplication (multiply $x by $y)
1345 $x->bdiv($y); # divide, set $i to quotient
1346 # return (quo,rem) or quo if scalar
1347
1348 $x->bmod($y); # modulus
1349 $x->bpow($y); # power of arguments (a**b)
1350 $x->blsft($y); # left shift
1351 $x->brsft($y); # right shift
1352 # return (quo,rem) or quo if scalar
1353
1354 $x->band($y); # bit-wise and
1355 $x->bior($y); # bit-wise inclusive or
1356 $x->bxor($y); # bit-wise exclusive or
1357 $x->bnot(); # bit-wise not (two's complement)
1358
1359 $x->bround($N); # accuracy: preserver $N digits
1360 $x->bfround($N); # precision: round to the $Nth digit
1361
1362 # The following do not modify their arguments:
1363
1364 bgcd(@values); # greatest common divisor
1365 blcm(@values); # lowest common multiplicator
1366
1367 $x->bstr(); # return string
1368 $x->bsstr(); # return string in scientific notation
1369
1370 $x->exponent(); # return exponent as BigInt
1371 $x->mantissa(); # return mantissa as BigInt
1372 $x->parts(); # return (mantissa,exponent) as BigInt
1373
1374 $x->length(); # number of digits (w/o sign and '.')
1375 ($l,$f) = $x->length(); # number of digits, and length of fraction
a5f75d66 1376
1377=head1 DESCRIPTION
1378
58cde26e 1379All operators (inlcuding basic math operations) are overloaded if you
1380declare your big floating point numbers as
a5f75d66 1381
58cde26e 1382 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1383
1384Operations with overloaded operators preserve the arguments, which is
1385exactly what you expect.
1386
1387=head2 Canonical notation
1388
1389Input to these routines are either BigFloat objects, or strings of the
1390following four forms:
a5f75d66 1391
1392=over 2
1393
58cde26e 1394=item *
1395
1396C</^[+-]\d+$/>
a5f75d66 1397
58cde26e 1398=item *
a5f75d66 1399
58cde26e 1400C</^[+-]\d+\.\d*$/>
a5f75d66 1401
58cde26e 1402=item *
a5f75d66 1403
58cde26e 1404C</^[+-]\d+E[+-]?\d+$/>
a5f75d66 1405
58cde26e 1406=item *
a5f75d66 1407
58cde26e 1408C</^[+-]\d*\.\d+E[+-]?\d+$/>
5d7098d5 1409
58cde26e 1410=back
1411
1412all with optional leading and trailing zeros and/or spaces. Additonally,
1413numbers are allowed to have an underscore between any two digits.
1414
1415Empty strings as well as other illegal numbers results in 'NaN'.
1416
1417bnorm() on a BigFloat object is now effectively a no-op, since the numbers
1418are always stored in normalized form. On a string, it creates a BigFloat
1419object.
1420
1421=head2 Output
1422
1423Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1424
1425The string output will always have leading and trailing zeros stripped and drop
1426a plus sign. C<bstr()> will give you always the form with a decimal point,
1427while C<bsstr()> (for scientific) gives you the scientific notation.
1428
1429 Input bstr() bsstr()
1430 '-0' '0' '0E1'
1431 ' -123 123 123' '-123123123' '-123123123E0'
1432 '00.0123' '0.0123' '123E-4'
1433 '123.45E-2' '1.2345' '12345E-4'
1434 '10E+3' '10000' '1E4'
1435
1436Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1437C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1438return either undef, <0, 0 or >0 and are suited for sort.
1439
1440Actual math is done by using BigInts to represent the mantissa and exponent.
1441The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
1442represent the result when input arguments are not numbers, as well as
1443the result of dividing by zero.
1444
1445=head2 C<mantissa()>, C<exponent()> and C<parts()>
1446
1447C<mantissa()> and C<exponent()> return the said parts of the BigFloat
1448as BigInts such that:
1449
1450 $m = $x->mantissa();
1451 $e = $x->exponent();
1452 $y = $m * ( 10 ** $e );
1453 print "ok\n" if $x == $y;
1454
1455C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
1456
1457A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
1458
1459Currently the mantissa is reduced as much as possible, favouring higher
1460exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
1461This might change in the future, so do not depend on it.
1462
1463=head2 Accuracy vs. Precision
1464
1465See also: L<Rounding|Rounding>.
1466
027dc388 1467Math::BigFloat supports both precision and accuracy. For a full documentation,
1468examples and tips on these topics please see the large section in
1469L<Math::BigInt>.
5d7098d5 1470
58cde26e 1471Since things like sqrt(2) or 1/3 must presented with a limited precision lest
1472a operation consumes all resources, each operation produces no more than
1473C<Math::BigFloat::precision()> digits.
1474
1475In case the result of one operation has more precision than specified,
1476it is rounded. The rounding mode taken is either the default mode, or the one
1477supplied to the operation after the I<scale>:
1478
1479 $x = Math::BigFloat->new(2);
1480 Math::BigFloat::precision(5); # 5 digits max
1481 $y = $x->copy()->bdiv(3); # will give 0.66666
1482 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1483 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
1484 Math::BigFloat::round_mode('zero');
1485 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1486
1487=head2 Rounding
1488
1489=over 2
1490
5dc6f178 1491=item ffround ( +$scale )
58cde26e 1492
0716bf9b 1493Rounds to the $scale'th place left from the '.', counting from the dot.
1494The first digit is numbered 1.
58cde26e 1495
5dc6f178 1496=item ffround ( -$scale )
58cde26e 1497
0716bf9b 1498Rounds to the $scale'th place right from the '.', counting from the dot.
58cde26e 1499
5dc6f178 1500=item ffround ( 0 )
1501
0716bf9b 1502Rounds to an integer.
5dc6f178 1503
1504=item fround ( +$scale )
1505
0716bf9b 1506Preserves accuracy to $scale digits from the left (aka significant digits)
1507and pads the rest with zeros. If the number is between 1 and -1, the
1508significant digits count from the first non-zero after the '.'
5dc6f178 1509
1510=item fround ( -$scale ) and fround ( 0 )
1511
0716bf9b 1512These are effetively no-ops.
5d7098d5 1513
a5f75d66 1514=back
1515
0716bf9b 1516All rounding functions take as a second parameter a rounding mode from one of
1517the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
58cde26e 1518
1519The default rounding mode is 'even'. By using
ee15d750 1520C<< Math::BigFloat::round_mode($round_mode); >> you can get and set the default
1521mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
0716bf9b 1522no longer supported.
b22b3e31 1523The second parameter to the round functions then overrides the default
0716bf9b 1524temporarily.
58cde26e 1525
1526The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
1527'trunc' as rounding mode to make it equivalent to:
1528
1529 $x = 2.5;
1530 $y = int($x) + 2;
1531
1532You can override this by passing the desired rounding mode as parameter to
1533C<as_number()>:
1534
1535 $x = Math::BigFloat->new(2.5);
1536 $y = $x->as_number('odd'); # $y = 3
1537
1538=head1 EXAMPLES
1539
58cde26e 1540 # not ready yet
58cde26e 1541
1542=head1 Autocreating constants
1543
1544After C<use Math::BigFloat ':constant'> all the floating point constants
1545in the given scope are converted to C<Math::BigFloat>. This conversion
1546happens at compile time.
1547
1548In particular
1549
1550 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
1551
1552prints the value of C<2E-100>. Note that without conversion of
1553constants the expression 2E-100 will be calculated as normal floating point
1554number.
1555
a5f75d66 1556=head1 BUGS
1557
58cde26e 1558=over 2
1559
1560=item *
1561
1562The following does not work yet:
1563
1564 $m = $x->mantissa();
1565 $e = $x->exponent();
1566 $y = $m * ( 10 ** $e );
1567 print "ok\n" if $x == $y;
1568
1569=item *
1570
1571There is no fmod() function yet.
1572
1573=back
1574
1575=head1 CAVEAT
1576
1577=over 1
1578
1579=item stringify, bstr()
1580
1581Both stringify and bstr() now drop the leading '+'. The old code would return
1582'+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
1583reasoning and details.
1584
1585=item bdiv
1586
1587The following will probably not do what you expect:
1588
1589 print $c->bdiv(123.456),"\n";
1590
1591It prints both quotient and reminder since print works in list context. Also,
1592bdiv() will modify $c, so be carefull. You probably want to use
1593
1594 print $c / 123.456,"\n";
1595 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
1596
1597instead.
1598
1599=item Modifying and =
1600
1601Beware of:
1602
1603 $x = Math::BigFloat->new(5);
1604 $y = $x;
1605
1606It will not do what you think, e.g. making a copy of $x. Instead it just makes
1607a second reference to the B<same> object and stores it in $y. Thus anything
1608that modifies $x will modify $y, and vice versa.
1609
1610 $x->bmul(2);
1611 print "$x, $y\n"; # prints '10, 10'
1612
1613If you want a true copy of $x, use:
1614
1615 $y = $x->copy();
1616
1617See also the documentation in L<overload> regarding C<=>.
1618
1619=item bpow
1620
1621C<bpow()> now modifies the first argument, unlike the old code which left
1622it alone and only returned the result. This is to be consistent with
1623C<badd()> etc. The first will modify $x, the second one won't:
1624
1625 print bpow($x,$i),"\n"; # modify $x
1626 print $x->bpow($i),"\n"; # ditto
1627 print $x ** $i,"\n"; # leave $x alone
1628
1629=back
1630
1631=head1 LICENSE
a5f75d66 1632
58cde26e 1633This program is free software; you may redistribute it and/or modify it under
1634the same terms as Perl itself.
5d7098d5 1635
58cde26e 1636=head1 AUTHORS
5d7098d5 1637
58cde26e 1638Mark Biggar, overloaded interface by Ilya Zakharevich.
1639Completely rewritten by Tels http://bloodgate.com in 2001.
a5f75d66 1640
a5f75d66 1641=cut
Unnamed repository; edit this file to name it for gitweb.