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