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