This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Math::Big* test tweaks to work better with core:
[perl5.git] / lib / Math / BigFloat.pm
CommitLineData
58cde26e
JH
1#!/usr/bin/perl -w
2
58cde26e
JH
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
JH
10# _cow: Copy-On-Write (NRY)
11
a0d0e21e
LW
12package Math::BigFloat;
13
ee15d750 14$VERSION = '1.23';
58cde26e
JH
15require 5.005;
16use Exporter;
0716bf9b 17use Math::BigInt qw/objectify/;
58cde26e
JH
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
JH
27 is_zero is_one sign
28 );
a0d0e21e 29
58cde26e
JH
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
JH
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
JH
42##############################################################################
43# global constants, flags and accessory
44
45use constant MB_NEVER_ROUND => 0x0001;
46
58cde26e
JH
47# are NaNs ok?
48my $NaNOK=1;
58cde26e
JH
49# constant for easier life
50my $nan = 'NaN';
58cde26e 51
ee15d750
JH
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
JH
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
JH
63 my %methods = map { $_ => 1 }
64 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
ee15d750
JH
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
JH
72 /;
73
ee15d750
JH
74 sub method_alias { return exists $methods{$_[0]||''}; }
75 sub method_hand_up { return exists $hand_ups{$_[0]||''}; }
a0d0e21e 76}
0e8b9368 77
58cde26e
JH
78##############################################################################
79# constructors
a0d0e21e 80
58cde26e
JH
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
JH
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
JH
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
JH
100 $self->{_e} = Math::BigInt->new(0);
101 $self->{_m}->babs();
102 $self->{sign} = $wanted->sign();
0716bf9b 103 return $self->bnorm();
58cde26e
JH
104 }
105 # got string
106 # handle '+inf', '-inf' first
ee15d750 107 if ($wanted =~ /^[+-]?inf$/)
58cde26e
JH
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
JH
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
JH
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
JH
137 if defined $accuracy || defined $precision;
138 return $self;
139 }
a0d0e21e 140
58cde26e
JH
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
JH
150 $self->{_m} = Math::BigInt->bzero();
151 $self->{_e} = Math::BigInt->bzero();
58cde26e 152 $self->{sign} = $nan;
58cde26e
JH
153 return $self;
154 }
a0d0e21e 155
58cde26e
JH
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
JH
162 $self = $class if !defined $self;
163 if (!ref($self))
164 {
165 my $c = $self; $self = {}; bless $self, $c;
166 }
574bacfe
JH
167 $self->{_m} = Math::BigInt->bzero();
168 $self->{_e} = Math::BigInt->bzero();
58cde26e 169 $self->{sign} = $sign.'inf';
58cde26e
JH
170 return $self;
171 }
a0d0e21e 172
574bacfe
JH
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
JH
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
JH
199 $self->{_m} = Math::BigInt->bzero();
200 $self->{_e} = Math::BigInt->bone();
58cde26e 201 $self->{sign} = '+';
58cde26e
JH
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
JH
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
JH
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
JH
221 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
222 return 'inf'; # +inf
58cde26e
JH
223 }
224
574bacfe
JH
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
JH
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
JH
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
JH
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
JH
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
JH
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
JH
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
JH
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
JH
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
JH
306 return $x->bsstr();
307 }
a0d0e21e 308
58cde26e
JH
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
JH
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
JH
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
JH
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
JH
349 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
350
574bacfe
JH
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
JH
357
358 # adjust so that exponents are equal
bd05a461
JH
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
JH
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
JH
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
JH
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
JH
392
393 # handle +-inf and NaN's
394 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]/)
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
JH
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
JH
451sub badd
452 {
453 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
454 # return result as BFLOAT
58cde26e
JH
455 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
456
574bacfe
JH
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
JH
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
JH
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
JH
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
JH
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
JH
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
JH
527 # subtract second arg from first, modify first
528 my ($self,$x,$y) = objectify(2,@_);
a0d0e21e 529
58cde26e
JH
530 $x->badd($y->bneg()); # badd does not leave internal zeros
531 $y->bneg(); # refix y, assumes no one reads $y in between
0716bf9b 532 return $x; # badd() already normalized and rounded
58cde26e
JH
533 }
534
535sub binc
536 {
537 # increment arg by one
ee15d750
JH
538 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
539 $x->badd($self->bone())->round($a,$p,$r);
58cde26e
JH
540 }
541
542sub bdec
543 {
544 # decrement arg by one
ee15d750
JH
545 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
546 $x->badd($self->bone('-'))->round($a,$p,$r);
58cde26e
JH
547 }
548
549sub blcm
550 {
ee15d750 551 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
58cde26e
JH
552 # does not modify arguments, but returns new object
553 # Lowest Common Multiplicator
58cde26e
JH
554
555 my ($self,@arg) = objectify(0,@_);
556 my $x = $self->new(shift @arg);
557 while (@arg) { $x = _lcm($x,shift @arg); }
558 $x;
559 }
560
561sub bgcd
562 {
ee15d750 563 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
58cde26e
JH
564 # does not modify arguments, but returns new object
565 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
58cde26e
JH
566
567 my ($self,@arg) = objectify(0,@_);
568 my $x = $self->new(shift @arg);
569 while (@arg) { $x = _gcd($x,shift @arg); }
570 $x;
571 }
572
573sub is_zero
574 {
ee15d750
JH
575 # return true if arg (BFLOAT or num_str) is zero (array '+', '0')
576 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
574bacfe
JH
577
578 return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
579 return 0;
58cde26e
JH
580 }
581
582sub is_one
583 {
ee15d750 584 # return true if arg (BFLOAT or num_str) is +1 (array '+', '1')
58cde26e 585 # or -1 if signis given
ee15d750
JH
586 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
587
588 my $sign = shift || ''; $sign = '+' if $sign ne '-';
589 return 1
590 if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
591 return 0;
58cde26e
JH
592 }
593
594sub is_odd
595 {
ee15d750
JH
596 # return true if arg (BFLOAT or num_str) is odd or false if even
597 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
0716bf9b
JH
598
599 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
ee15d750
JH
600 return 1 if ($x->{_e}->is_zero() && $x->{_m}->is_odd());
601 return 0;
58cde26e
JH
602 }
603
604sub is_even
605 {
b22b3e31 606 # return true if arg (BINT or num_str) is even or false if odd
ee15d750 607 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
0716bf9b
JH
608
609 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
610 return 1 if $x->{_m}->is_zero(); # 0e1 is even
ee15d750
JH
611 return 1 if ($x->{_e}->is_zero() && $x->{_m}->is_even()); # 123.45 is never
612 return 0;
58cde26e
JH
613 }
614
615sub bmul
616 {
617 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
618 # (BINT or num_str, BINT or num_str) return BINT
619 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
58cde26e 620
0716bf9b 621 # print "mbf bmul $x->{_m}e$x->{_e} $y->{_m}e$y->{_e}\n";
58cde26e
JH
622 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
623
574bacfe
JH
624 # handle result = 0
625 return $x->bzero() if $x->is_zero() || $y->is_zero();
626 # inf handling
627 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
628 {
629 # result will always be +-inf:
630 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
631 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
632 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
633 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
634 return $x->binf('-');
635 }
636
58cde26e
JH
637 # aEb * cEd = (a*c)E(b+d)
638 $x->{_m} = $x->{_m} * $y->{_m};
639 #print "m: $x->{_m}\n";
640 $x->{_e} = $x->{_e} + $y->{_e};
641 #print "e: $x->{_m}\n";
642 # adjust sign:
643 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
644 #print "s: $x->{sign}\n";
0716bf9b 645 $x->bnorm();
58cde26e
JH
646 return $x->round($a,$p,$r,$y);
647 }
648
649sub bdiv
650 {
651 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
652 # (BFLOAT,BFLOAT) (quo,rem) or BINT (only rem)
653 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
654
ee15d750 655
574bacfe
JH
656 # x / +-inf => 0, reminder x
657 return wantarray ? ($x->bzero(),$x->copy()) : $x->bzero()
658 if $y->{sign} =~ /^[+-]inf$/;
659
660 # NaN if x == NaN or y == NaN or x==y==0
58cde26e 661 return wantarray ? ($x->bnan(),bnan()) : $x->bnan()
574bacfe
JH
662 if (($x->is_nan() || $y->is_nan()) ||
663 ($x->is_zero() && $y->is_zero()));
664
665 # 5 / 0 => +inf, -6 / 0 => -inf
666 return wantarray
667 ? ($x->binf($x->{sign}),$self->bnan()) : $x->binf($x->{sign})
668 if ($x->{sign} =~ /^[+-]$/ && $y->is_zero());
0716bf9b 669
ee15d750
JH
670 # promote BigInts and it's subclasses (except when already a BigFloat)
671 $y = $self->new($y) unless $y->isa('Math::BigFloat');
672
673 # old, broken way
674 # $y = $class->new($y) if ref($y) ne $self; # promote bigints
0716bf9b
JH
675
676 # print "mbf bdiv $x ",ref($x)," ",$y," ",ref($y),"\n";
58cde26e 677 # we need to limit the accuracy to protect against overflow
ee15d750 678
574bacfe 679 my $fallback = 0;
ee15d750
JH
680 my $scale = 0;
681# print "s=$scale a=",$a||'undef'," p=",$p||'undef'," r=",$r||'undef',"\n";
682 my @params = $x->_find_round_parameters($a,$p,$r,$y);
683
684 # no rounding at all, so must use fallback
685 if (scalar @params == 1)
58cde26e 686 {
0716bf9b 687 # simulate old behaviour
ee15d750
JH
688 $scale = $self->div_scale()+1; # at least one more for proper round
689 $params[1] = $self->div_scale(); # and round to it as accuracy
690 $params[3] = $r; # round mode by caller or undef
691 $fallback = 1; # to clear a/p afterwards
692 }
693 else
694 {
695 # the 4 below is empirical, and there might be cases where it is not
696 # enough...
697 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
a0d0e21e 698 }
ee15d750 699 # print "s=$scale a=",$params[1]||'undef'," p=",$params[2]||'undef'," f=$fallback\n";
0716bf9b 700 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
58cde26e 701 $scale = $lx if $lx > $scale;
58cde26e 702 $scale = $ly if $ly > $scale;
ee15d750 703# print "scale $scale $lx $ly\n";
0716bf9b
JH
704 my $diff = $ly - $lx;
705 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
a0d0e21e 706
58cde26e
JH
707 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
708
709 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
a0d0e21e 710
58cde26e
JH
711 # check for / +-1 ( +/- 1E0)
712 if ($y->is_one())
713 {
ee15d750 714 return wantarray ? ($x,$self->bzero()) : $x;
a0d0e21e 715 }
a5f75d66 716
ee15d750 717 # calculate the result to $scale digits and then round it
58cde26e 718 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
ee15d750 719 #$scale = 82;
58cde26e 720 #print "self: $self x: $x ref(x) ", ref($x)," m: $x->{_m}\n";
58cde26e 721 $x->{_m}->blsft($scale,10);
0716bf9b 722 #print "m: $x->{_m} $y->{_m}\n";
58cde26e
JH
723 $x->{_m}->bdiv( $y->{_m} ); # a/c
724 #print "m: $x->{_m}\n";
ee15d750 725 #print "e: $x->{_e} $y->{_e} ",$scale,"\n";
58cde26e
JH
726 $x->{_e}->bsub($y->{_e}); # b-d
727 #print "e: $x->{_e}\n";
728 $x->{_e}->bsub($scale); # correct for 10**scale
0716bf9b
JH
729 #print "after div: m: $x->{_m} e: $x->{_e}\n";
730 $x->bnorm(); # remove trailing 0's
ee15d750
JH
731 #print "after norm: m: $x->{_m} e: $x->{_e}\n";
732
733 # shortcut to not run trough _find_round_parameters again
734 if (defined $params[1])
735 {
736 $x->bround($params[1],undef,$params[3]); # then round accordingly
737 }
738 else
739 {
740 $x->bfround($params[2],$params[3]); # then round accordingly
741 }
574bacfe
JH
742 if ($fallback)
743 {
744 # clear a/p after round, since user did not request it
ee15d750 745 $x->{_a} = undef; $x->{_p} = undef;
574bacfe 746 }
0716bf9b 747
58cde26e
JH
748 if (wantarray)
749 {
750 my $rem = $x->copy();
ee15d750 751 $rem->bmod($y,$params[1],$params[2],$params[3]);
574bacfe
JH
752 if ($fallback)
753 {
754 # clear a/p after round, since user did not request it
ee15d750 755 $rem->{_a} = undef; $rem->{_p} = undef;
574bacfe 756 }
0716bf9b 757 return ($x,$rem);
58cde26e
JH
758 }
759 return $x;
760 }
a0d0e21e 761
58cde26e
JH
762sub bmod
763 {
764 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
765 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
a0d0e21e 766
58cde26e
JH
767 return $x->bnan() if ($x->{sign} eq $nan || $y->is_nan() || $y->is_zero());
768 return $x->bzero() if $y->is_one();
769
770 # XXX tels: not done yet
771 return $x->round($a,$p,$r,$y);
772 }
773
774sub bsqrt
775 {
0716bf9b
JH
776 # calculate square root; this should probably
777 # use a different test to see whether the accuracy we want is...
ee15d750 778 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 779
0716bf9b
JH
780 return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
781 return $x if $x->{sign} eq '+inf'; # +inf
58cde26e
JH
782 return $x if $x->is_zero() || $x == 1;
783
ee15d750
JH
784 # we need to limit the accuracy to protect against overflow (ignore $p)
785 my ($scale) = $x->_scale_a($self->accuracy(),$self->round_mode,$a,$r);
574bacfe 786 my $fallback = 0;
0716bf9b
JH
787 if (!defined $scale)
788 {
789 # simulate old behaviour
ee15d750
JH
790 $scale = $self->div_scale()+1; # one more for proper riund
791 $a = $self->div_scale(); # and round to it
792 $fallback = 1; # to clear a/p afterwards
0716bf9b
JH
793 }
794 my $lx = $x->{_m}->length();
795 $scale = $lx if $scale < $lx;
58cde26e
JH
796 my $e = Math::BigFloat->new("1E-$scale"); # make test variable
797 return $x->bnan() if $e->sign() eq 'NaN';
798
58cde26e 799 # start with some reasonable guess
0716bf9b 800 #$x *= 10 ** ($len - $org->{_e}); $x /= 2; # !?!?
574bacfe 801 $lx = $lx+$x->{_e};
0716bf9b
JH
802 $lx = 1 if $lx < 1;
803 my $gs = Math::BigFloat->new('1'. ('0' x $lx));
804
ee15d750 805# print "first guess: $gs (x $x) scale $scale\n";
58cde26e
JH
806
807 my $diff = $e;
808 my $y = $x->copy();
809 my $two = Math::BigFloat->new(2);
ee15d750
JH
810 # promote BigInts and it's subclasses (except when already a BigFloat)
811 $y = $self->new($y) unless $y->isa('Math::BigFloat');
812 # old, broken way
813 # $x = Math::BigFloat->new($x) if ref($x) ne $class; # promote BigInts
814 my $rem;
58cde26e
JH
815 # $scale = 2;
816 while ($diff >= $e)
817 {
58cde26e 818 return $x->bnan() if $gs->is_zero();
ee15d750
JH
819 $rem = $y->copy(); $rem->bdiv($gs,$scale);
820 #print "y $y gs $gs ($gs->{_a}) rem (y/gs)\n $rem\n";
821 $x = ($rem + $gs);
822 #print "x $x rem $rem gs $gs gsa: $gs->{_a}\n";
823 $x->bdiv($two,$scale);
824 #print "x $x (/2)\n";
58cde26e 825 $diff = $x->copy()->bsub($gs)->babs();
58cde26e 826 $gs = $x->copy();
a0d0e21e 827 }
ee15d750 828# print "before $x $x->{_a} ",$a||'a undef'," ",$p||'p undef',"\n";
0716bf9b 829 $x->round($a,$p,$r);
ee15d750 830# print "after $x $x->{_a} ",$a||'a undef'," ",$p||'p undef',"\n";
574bacfe
JH
831 if ($fallback)
832 {
833 # clear a/p after round, since user did not request it
ee15d750 834 $x->{_a} = undef; $x->{_p} = undef;
574bacfe
JH
835 }
836 $x;
58cde26e
JH
837 }
838
839sub bpow
840 {
841 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
842 # compute power of two numbers, second arg is used as integer
843 # modifies first argument
844
845 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
846
0716bf9b 847 return $x if $x->{sign} =~ /^[+-]inf$/;
58cde26e 848 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
574bacfe 849 return $x->bone() if $y->is_zero();
58cde26e 850 return $x if $x->is_one() || $y->is_one();
ee15d750 851 my $y1 = $y->as_number(); # make bigint (trunc)
58cde26e
JH
852 if ($x == -1)
853 {
854 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
0716bf9b 855 return $y1->is_odd() ? $x : $x->babs(1);
288d023a 856 }
58cde26e 857 return $x if $x->is_zero() && $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
574bacfe
JH
858 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
859 return $x->binf() if $x->is_zero() && $y->{sign} eq '-';
58cde26e
JH
860
861 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
862 $y1->babs();
863 $x->{_m}->bpow($y1);
864 $x->{_e}->bmul($y1);
865 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
866 $x->bnorm();
867 if ($y->{sign} eq '-')
868 {
869 # modify $x in place!
0716bf9b 870 my $z = $x->copy(); $x->bzero()->binc();
58cde26e 871 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
a0d0e21e 872 }
58cde26e
JH
873 return $x->round($a,$p,$r,$y);
874 }
875
876###############################################################################
877# rounding functions
878
879sub bfround
880 {
881 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
882 # $n == 0 means round to integer
883 # expects and returns normalized numbers!
ee15d750 884 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
a0d0e21e 885
58cde26e
JH
886 return $x if $x->modify('bfround');
887
ee15d750 888 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
58cde26e
JH
889 return $x if !defined $scale; # no-op
890
574bacfe
JH
891 # never round a 0, +-inf, NaN
892 return $x if $x->{sign} !~ /^[+-]$/ || $x->is_zero();
58cde26e 893 # print "MBF bfround $x to scale $scale mode $mode\n";
58cde26e 894
ee15d750
JH
895 # don't round if x already has lower precision
896 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
897
898 $x->{_p} = $scale; # remember round in any case
899 $x->{_a} = undef; # and clear A
58cde26e
JH
900 if ($scale < 0)
901 {
902 # print "bfround scale $scale e $x->{_e}\n";
903 # round right from the '.'
904 return $x if $x->{_e} >= 0; # nothing to round
905 $scale = -$scale; # positive for simplicity
906 my $len = $x->{_m}->length(); # length of mantissa
907 my $dad = -$x->{_e}; # digits after dot
908 my $zad = 0; # zeros after dot
909 $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
ee15d750 910 #print "scale $scale dad $dad zad $zad len $len\n";
58cde26e
JH
911
912 # number bsstr len zad dad
913 # 0.123 123e-3 3 0 3
914 # 0.0123 123e-4 3 1 4
915 # 0.001 1e-3 1 2 3
916 # 1.23 123e-2 3 0 2
917 # 1.2345 12345e-4 5 0 4
918
919 # do not round after/right of the $dad
920 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
921
ee15d750
JH
922 # round to zero if rounding inside the $zad, but not for last zero like:
923 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
924 return $x->bzero() if $scale < $zad;
925 if ($scale == $zad) # for 0.006, scale -3 and trunc
58cde26e 926 {
ee15d750 927 $scale = -$len-1;
58cde26e
JH
928 }
929 else
930 {
931 # adjust round-point to be inside mantissa
932 if ($zad != 0)
933 {
934 $scale = $scale-$zad;
935 }
936 else
937 {
938 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
939 $scale = $dbd+$scale;
940 }
941 }
942 # print "round to $x->{_m} to $scale\n";
a0d0e21e 943 }
58cde26e
JH
944 else
945 {
946 # 123 => 100 means length(123) = 3 - $scale (2) => 1
a5f75d66 947
58cde26e
JH
948 # calculate digits before dot
949 my $dbt = $x->{_m}->length(); $dbt += $x->{_e} if $x->{_e}->sign() eq '-';
ee15d750
JH
950 # if not enough digits before dot, round to zero
951 return $x->bzero() if ($scale > $dbt) && ($dbt < 0);
952 # scale always >= 0 here
953 if ($dbt == 0)
58cde26e
JH
954 {
955 # 0.49->bfround(1): scale == 1, dbt == 0: => 0.0
956 # 0.51->bfround(0): scale == 0, dbt == 0: => 1.0
957 # 0.5->bfround(0): scale == 0, dbt == 0: => 0
958 # 0.05->bfround(0): scale == 0, dbt == 0: => 0
959 # print "$scale $dbt $x->{_m}\n";
960 $scale = -$x->{_m}->length();
961 }
962 elsif ($dbt > 0)
963 {
964 # correct by subtracting scale
965 $scale = $dbt - $scale;
966 }
967 else
968 {
969 $scale = $x->{_m}->length() - $scale;
970 }
a0d0e21e 971 }
574bacfe
JH
972 # print "using $scale for $x->{_m} with '$mode'\n";
973 # pass sign to bround for rounding modes '+inf' and '-inf'
58cde26e
JH
974 $x->{_m}->{sign} = $x->{sign};
975 $x->{_m}->bround($scale,$mode);
976 $x->{_m}->{sign} = '+'; # fix sign back
977 $x->bnorm();
978 }
979
980sub bround
981 {
982 # accuracy: preserve $N digits, and overwrite the rest with 0's
ee15d750
JH
983 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
984
985 die ('bround() needs positive accuracy') if ($_[0] || 0) < 0;
58cde26e 986
ee15d750
JH
987 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
988 return $x if !defined $scale; # no-op
989
58cde26e 990 return $x if $x->modify('bround');
ee15d750
JH
991
992 # scale is now either $x->{_a}, $accuracy, or the user parameter
993 # test whether $x already has lower accuracy, do nothing in this case
994 # but do round if the accuracy is the same, since a math operation might
995 # want to round a number with A=5 to 5 digits afterwards again
996 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
58cde26e
JH
997
998 # print "bround $scale $mode\n";
999 # 0 => return all digits, scale < 0 makes no sense
1000 return $x if ($scale <= 0);
574bacfe
JH
1001 # never round a 0, +-inf, NaN
1002 return $x if $x->{sign} !~ /^[+-]$/ || $x->is_zero();
58cde26e
JH
1003
1004 # if $e longer than $m, we have 0.0000xxxyyy style number, and must
1005 # subtract the delta from scale, to simulate keeping the zeros
1006 # -5 +5 => 1; -10 +5 => -4
1007 my $delta = $x->{_e} + $x->{_m}->length() + 1;
58cde26e
JH
1008
1009 # if we should keep more digits than the mantissa has, do nothing
1010 return $x if $x->{_m}->length() <= $scale;
f216259d 1011
58cde26e
JH
1012 # pass sign to bround for '+inf' and '-inf' rounding modes
1013 $x->{_m}->{sign} = $x->{sign};
1014 $x->{_m}->bround($scale,$mode); # round mantissa
1015 $x->{_m}->{sign} = '+'; # fix sign back
ee15d750
JH
1016 $x->{_a} = $scale; # remember rounding
1017 $x->{_p} = undef; # and clear P
574bacfe 1018 $x->bnorm(); # del trailing zeros gen. by bround()
58cde26e
JH
1019 }
1020
1021sub bfloor
1022 {
1023 # return integer less or equal then $x
ee15d750 1024 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e
JH
1025
1026 return $x if $x->modify('bfloor');
1027
1028 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1029
1030 # if $x has digits after dot
1031 if ($x->{_e}->{sign} eq '-')
1032 {
1033 $x->{_m}->brsft(-$x->{_e},10);
1034 $x->{_e}->bzero();
1035 $x-- if $x->{sign} eq '-';
f216259d 1036 }
58cde26e
JH
1037 return $x->round($a,$p,$r);
1038 }
288d023a 1039
58cde26e
JH
1040sub bceil
1041 {
1042 # return integer greater or equal then $x
ee15d750 1043 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e
JH
1044
1045 return $x if $x->modify('bceil');
1046 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1047
1048 # if $x has digits after dot
1049 if ($x->{_e}->{sign} eq '-')
1050 {
1051 $x->{_m}->brsft(-$x->{_e},10);
1052 $x->{_e}->bzero();
1053 $x++ if $x->{sign} eq '+';
a0d0e21e 1054 }
58cde26e
JH
1055 return $x->round($a,$p,$r);
1056 }
1057
1058###############################################################################
a5f75d66 1059
58cde26e
JH
1060sub DESTROY
1061 {
ee15d750 1062 # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
58cde26e
JH
1063 }
1064
1065sub AUTOLOAD
1066 {
1067 # make fxxx and bxxx work
1068 # my $self = $_[0];
1069 my $name = $AUTOLOAD;
1070
1071 $name =~ s/.*:://; # split package
1072 #print "$name\n";
ee15d750
JH
1073 no strict 'refs';
1074 if (!method_alias($name))
58cde26e 1075 {
ee15d750
JH
1076 if (!defined $name)
1077 {
1078 # delayed load of Carp and avoid recursion
1079 require Carp;
1080 Carp::croak ("Can't call a method without name");
1081 }
1082 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1083 if (!method_hand_up($name))
1084 {
1085 # delayed load of Carp and avoid recursion
1086 require Carp;
1087 Carp::croak ("Can't call $class\-\>$name, not a valid method");
1088 }
1089 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1090 $name =~ s/^f/b/;
1091 return &{'Math::BigInt'."::$name"}(@_);
a0d0e21e 1092 }
58cde26e
JH
1093 my $bname = $name; $bname =~ s/^f/b/;
1094 *{$class."\:\:$name"} = \&$bname;
1095 &$bname; # uses @_
1096 }
1097
1098sub exponent
1099 {
1100 # return a copy of the exponent
ee15d750 1101 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 1102
ee15d750
JH
1103 if ($x->{sign} !~ /^[+-]$/)
1104 {
1105 my $s = $x->{sign}; $s =~ s/^[+-]//;
1106 return $self->new($s); # -inf, +inf => +inf
1107 }
1108 return $x->{_e}->copy();
58cde26e
JH
1109 }
1110
1111sub mantissa
1112 {
1113 # return a copy of the mantissa
ee15d750 1114 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 1115
ee15d750
JH
1116 if ($x->{sign} !~ /^[+-]$/)
1117 {
1118 my $s = $x->{sign}; $s =~ s/^[+]//;
1119 return $self->new($s); # -inf, +inf => +inf
1120 }
1121 my $m = $x->{_m}->copy(); # faster than going via bstr()
1122 $m->bneg() if $x->{sign} eq '-';
58cde26e
JH
1123
1124 return $m;
1125 }
1126
1127sub parts
1128 {
1129 # return a copy of both the exponent and the mantissa
ee15d750 1130 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 1131
ee15d750
JH
1132 if ($x->{sign} !~ /^[+-]$/)
1133 {
1134 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
1135 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
1136 }
1137 my $m = $x->{_m}->copy(); # faster than going via bstr()
1138 $m->bneg() if $x->{sign} eq '-';
1139 return ($m,$x->{_e}->copy());
58cde26e
JH
1140 }
1141
1142##############################################################################
1143# private stuff (internal use only)
1144
58cde26e
JH
1145sub import
1146 {
1147 my $self = shift;
58cde26e
JH
1148 for ( my $i = 0; $i < @_ ; $i++ )
1149 {
1150 if ( $_[$i] eq ':constant' )
1151 {
1152 # this rest causes overlord er load to step in
1153 # print "overload @_\n";
1154 overload::constant float => sub { $self->new(shift); };
1155 splice @_, $i, 1; last;
1156 }
1157 }
1158 # any non :constant stuff is handled by our parent, Exporter
1159 # even if @_ is empty, to give it a chance
1160 #$self->SUPER::import(@_); # does not work (would call MBI)
1161 $self->export_to_level(1,$self,@_); # need this instead
1162 }
1163
1164sub bnorm
1165 {
1166 # adjust m and e so that m is smallest possible
1167 # round number according to accuracy and precision settings
ee15d750 1168 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 1169
0716bf9b 1170 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
58cde26e
JH
1171
1172 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
1173 if ($zeros != 0)
1174 {
1175 $x->{_m}->brsft($zeros,10); $x->{_e} += $zeros;
1176 }
ee15d750
JH
1177 # for something like 0Ey, set y to 1, and -0 => +0
1178 $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
1179 # this is to prevent automatically rounding when MBI's globals are set
0716bf9b
JH
1180 $x->{_m}->{_f} = MB_NEVER_ROUND;
1181 $x->{_e}->{_f} = MB_NEVER_ROUND;
ee15d750
JH
1182 # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
1183 $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
1184 $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
0716bf9b 1185 return $x; # MBI bnorm is no-op
58cde26e
JH
1186 }
1187
1188##############################################################################
1189# internal calculation routines
1190
1191sub as_number
1192 {
1193 # return a bigint representation of this BigFloat number
ee15d750 1194 my $x = shift; my $class = ref($x) || $x; $x = $class->new(shift) unless ref($x);
58cde26e
JH
1195
1196 my $z;
1197 if ($x->{_e}->is_zero())
1198 {
1199 $z = $x->{_m}->copy();
1200 $z->{sign} = $x->{sign};
1201 return $z;
1202 }
0716bf9b 1203 $z = $x->{_m}->copy();
58cde26e
JH
1204 if ($x->{_e} < 0)
1205 {
0716bf9b
JH
1206 $z->brsft(-$x->{_e},10);
1207 }
1208 else
1209 {
1210 $z->blsft($x->{_e},10);
58cde26e 1211 }
58cde26e
JH
1212 $z->{sign} = $x->{sign};
1213 return $z;
1214 }
1215
1216sub length
1217 {
ee15d750
JH
1218 my $x = shift;
1219 my $class = ref($x) || $x;
1220 $x = $class->new(shift) unless ref($x);
58cde26e 1221
ee15d750 1222 return 1 if $x->{_m}->is_zero();
58cde26e
JH
1223 my $len = $x->{_m}->length();
1224 $len += $x->{_e} if $x->{_e}->sign() eq '+';
1225 if (wantarray())
1226 {
1227 my $t = Math::BigInt::bzero();
1228 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1229 return ($len,$t);
1230 }
1231 return $len;
1232 }
a0d0e21e
LW
1233
12341;
a5f75d66
AD
1235__END__
1236
1237=head1 NAME
1238
58cde26e 1239Math::BigFloat - Arbitrary size floating point math package
a5f75d66
AD
1240
1241=head1 SYNOPSIS
1242
a2008d6d 1243 use Math::BigFloat;
58cde26e
JH
1244
1245 # Number creation
1246 $x = Math::BigInt->new($str); # defaults to 0
1247 $nan = Math::BigInt->bnan(); # create a NotANumber
1248 $zero = Math::BigInt->bzero();# create a "+0"
1249
1250 # Testing
1251 $x->is_zero(); # return whether arg is zero or not
0716bf9b
JH
1252 $x->is_nan(); # return whether arg is NaN or not
1253 $x->is_one(); # true if arg is +1
1254 $x->is_one('-'); # true if arg is -1
1255 $x->is_odd(); # true if odd, false for even
1256 $x->is_even(); # true if even, false for odd
1257 $x->is_positive(); # true if >= 0
1258 $x->is_negative(); # true if < 0
1259 $x->is_inf(sign) # true if +inf or -inf (sign default '+')
58cde26e
JH
1260 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
1261 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
1262 $x->sign(); # return the sign, either +,- or NaN
1263
1264 # The following all modify their first argument:
1265
1266 # set
1267 $x->bzero(); # set $i to 0
1268 $x->bnan(); # set $i to NaN
1269
1270 $x->bneg(); # negation
1271 $x->babs(); # absolute value
1272 $x->bnorm(); # normalize (no-op)
1273 $x->bnot(); # two's complement (bit wise not)
1274 $x->binc(); # increment x by 1
1275 $x->bdec(); # decrement x by 1
1276
1277 $x->badd($y); # addition (add $y to $x)
1278 $x->bsub($y); # subtraction (subtract $y from $x)
1279 $x->bmul($y); # multiplication (multiply $x by $y)
1280 $x->bdiv($y); # divide, set $i to quotient
1281 # return (quo,rem) or quo if scalar
1282
1283 $x->bmod($y); # modulus
1284 $x->bpow($y); # power of arguments (a**b)
1285 $x->blsft($y); # left shift
1286 $x->brsft($y); # right shift
1287 # return (quo,rem) or quo if scalar
1288
1289 $x->band($y); # bit-wise and
1290 $x->bior($y); # bit-wise inclusive or
1291 $x->bxor($y); # bit-wise exclusive or
1292 $x->bnot(); # bit-wise not (two's complement)
1293
1294 $x->bround($N); # accuracy: preserver $N digits
1295 $x->bfround($N); # precision: round to the $Nth digit
1296
1297 # The following do not modify their arguments:
1298
1299 bgcd(@values); # greatest common divisor
1300 blcm(@values); # lowest common multiplicator
1301
1302 $x->bstr(); # return string
1303 $x->bsstr(); # return string in scientific notation
1304
1305 $x->exponent(); # return exponent as BigInt
1306 $x->mantissa(); # return mantissa as BigInt
1307 $x->parts(); # return (mantissa,exponent) as BigInt
1308
1309 $x->length(); # number of digits (w/o sign and '.')
1310 ($l,$f) = $x->length(); # number of digits, and length of fraction
a5f75d66
AD
1311
1312=head1 DESCRIPTION
1313
58cde26e
JH
1314All operators (inlcuding basic math operations) are overloaded if you
1315declare your big floating point numbers as
a5f75d66 1316
58cde26e
JH
1317 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1318
1319Operations with overloaded operators preserve the arguments, which is
1320exactly what you expect.
1321
1322=head2 Canonical notation
1323
1324Input to these routines are either BigFloat objects, or strings of the
1325following four forms:
a5f75d66
AD
1326
1327=over 2
1328
58cde26e
JH
1329=item *
1330
1331C</^[+-]\d+$/>
a5f75d66 1332
58cde26e 1333=item *
a5f75d66 1334
58cde26e 1335C</^[+-]\d+\.\d*$/>
a5f75d66 1336
58cde26e 1337=item *
a5f75d66 1338
58cde26e 1339C</^[+-]\d+E[+-]?\d+$/>
a5f75d66 1340
58cde26e 1341=item *
a5f75d66 1342
58cde26e 1343C</^[+-]\d*\.\d+E[+-]?\d+$/>
5d7098d5 1344
58cde26e
JH
1345=back
1346
1347all with optional leading and trailing zeros and/or spaces. Additonally,
1348numbers are allowed to have an underscore between any two digits.
1349
1350Empty strings as well as other illegal numbers results in 'NaN'.
1351
1352bnorm() on a BigFloat object is now effectively a no-op, since the numbers
1353are always stored in normalized form. On a string, it creates a BigFloat
1354object.
1355
1356=head2 Output
1357
1358Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1359
1360The string output will always have leading and trailing zeros stripped and drop
1361a plus sign. C<bstr()> will give you always the form with a decimal point,
1362while C<bsstr()> (for scientific) gives you the scientific notation.
1363
1364 Input bstr() bsstr()
1365 '-0' '0' '0E1'
1366 ' -123 123 123' '-123123123' '-123123123E0'
1367 '00.0123' '0.0123' '123E-4'
1368 '123.45E-2' '1.2345' '12345E-4'
1369 '10E+3' '10000' '1E4'
1370
1371Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1372C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1373return either undef, <0, 0 or >0 and are suited for sort.
1374
1375Actual math is done by using BigInts to represent the mantissa and exponent.
1376The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
1377represent the result when input arguments are not numbers, as well as
1378the result of dividing by zero.
1379
1380=head2 C<mantissa()>, C<exponent()> and C<parts()>
1381
1382C<mantissa()> and C<exponent()> return the said parts of the BigFloat
1383as BigInts such that:
1384
1385 $m = $x->mantissa();
1386 $e = $x->exponent();
1387 $y = $m * ( 10 ** $e );
1388 print "ok\n" if $x == $y;
1389
1390C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
1391
1392A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
1393
1394Currently the mantissa is reduced as much as possible, favouring higher
1395exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
1396This might change in the future, so do not depend on it.
1397
1398=head2 Accuracy vs. Precision
1399
1400See also: L<Rounding|Rounding>.
1401
1402Math::BigFloat supports both precision and accuracy. (here should follow
1403a short description of both).
5d7098d5 1404
58cde26e
JH
1405Precision: digits after the '.', laber, schwad
1406Accuracy: Significant digits blah blah
5d7098d5 1407
58cde26e
JH
1408Since things like sqrt(2) or 1/3 must presented with a limited precision lest
1409a operation consumes all resources, each operation produces no more than
1410C<Math::BigFloat::precision()> digits.
1411
1412In case the result of one operation has more precision than specified,
1413it is rounded. The rounding mode taken is either the default mode, or the one
1414supplied to the operation after the I<scale>:
1415
1416 $x = Math::BigFloat->new(2);
1417 Math::BigFloat::precision(5); # 5 digits max
1418 $y = $x->copy()->bdiv(3); # will give 0.66666
1419 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1420 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
1421 Math::BigFloat::round_mode('zero');
1422 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1423
1424=head2 Rounding
1425
1426=over 2
1427
5dc6f178 1428=item ffround ( +$scale )
58cde26e 1429
0716bf9b
JH
1430Rounds to the $scale'th place left from the '.', counting from the dot.
1431The first digit is numbered 1.
58cde26e 1432
5dc6f178 1433=item ffround ( -$scale )
58cde26e 1434
0716bf9b 1435Rounds to the $scale'th place right from the '.', counting from the dot.
58cde26e 1436
5dc6f178
JH
1437=item ffround ( 0 )
1438
0716bf9b 1439Rounds to an integer.
5dc6f178
JH
1440
1441=item fround ( +$scale )
1442
0716bf9b
JH
1443Preserves accuracy to $scale digits from the left (aka significant digits)
1444and pads the rest with zeros. If the number is between 1 and -1, the
1445significant digits count from the first non-zero after the '.'
5dc6f178
JH
1446
1447=item fround ( -$scale ) and fround ( 0 )
1448
0716bf9b 1449These are effetively no-ops.
5d7098d5 1450
a5f75d66
AD
1451=back
1452
0716bf9b
JH
1453All rounding functions take as a second parameter a rounding mode from one of
1454the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
58cde26e
JH
1455
1456The default rounding mode is 'even'. By using
ee15d750
JH
1457C<< Math::BigFloat::round_mode($round_mode); >> you can get and set the default
1458mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
0716bf9b 1459no longer supported.
b22b3e31 1460The second parameter to the round functions then overrides the default
0716bf9b 1461temporarily.
58cde26e
JH
1462
1463The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
1464'trunc' as rounding mode to make it equivalent to:
1465
1466 $x = 2.5;
1467 $y = int($x) + 2;
1468
1469You can override this by passing the desired rounding mode as parameter to
1470C<as_number()>:
1471
1472 $x = Math::BigFloat->new(2.5);
1473 $y = $x->as_number('odd'); # $y = 3
1474
1475=head1 EXAMPLES
1476
58cde26e 1477 # not ready yet
58cde26e
JH
1478
1479=head1 Autocreating constants
1480
1481After C<use Math::BigFloat ':constant'> all the floating point constants
1482in the given scope are converted to C<Math::BigFloat>. This conversion
1483happens at compile time.
1484
1485In particular
1486
1487 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
1488
1489prints the value of C<2E-100>. Note that without conversion of
1490constants the expression 2E-100 will be calculated as normal floating point
1491number.
1492
a5f75d66
AD
1493=head1 BUGS
1494
58cde26e
JH
1495=over 2
1496
1497=item *
1498
1499The following does not work yet:
1500
1501 $m = $x->mantissa();
1502 $e = $x->exponent();
1503 $y = $m * ( 10 ** $e );
1504 print "ok\n" if $x == $y;
1505
1506=item *
1507
1508There is no fmod() function yet.
1509
1510=back
1511
1512=head1 CAVEAT
1513
1514=over 1
1515
1516=item stringify, bstr()
1517
1518Both stringify and bstr() now drop the leading '+'. The old code would return
1519'+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
1520reasoning and details.
1521
1522=item bdiv
1523
1524The following will probably not do what you expect:
1525
1526 print $c->bdiv(123.456),"\n";
1527
1528It prints both quotient and reminder since print works in list context. Also,
1529bdiv() will modify $c, so be carefull. You probably want to use
1530
1531 print $c / 123.456,"\n";
1532 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
1533
1534instead.
1535
1536=item Modifying and =
1537
1538Beware of:
1539
1540 $x = Math::BigFloat->new(5);
1541 $y = $x;
1542
1543It will not do what you think, e.g. making a copy of $x. Instead it just makes
1544a second reference to the B<same> object and stores it in $y. Thus anything
1545that modifies $x will modify $y, and vice versa.
1546
1547 $x->bmul(2);
1548 print "$x, $y\n"; # prints '10, 10'
1549
1550If you want a true copy of $x, use:
1551
1552 $y = $x->copy();
1553
1554See also the documentation in L<overload> regarding C<=>.
1555
1556=item bpow
1557
1558C<bpow()> now modifies the first argument, unlike the old code which left
1559it alone and only returned the result. This is to be consistent with
1560C<badd()> etc. The first will modify $x, the second one won't:
1561
1562 print bpow($x,$i),"\n"; # modify $x
1563 print $x->bpow($i),"\n"; # ditto
1564 print $x ** $i,"\n"; # leave $x alone
1565
1566=back
1567
1568=head1 LICENSE
a5f75d66 1569
58cde26e
JH
1570This program is free software; you may redistribute it and/or modify it under
1571the same terms as Perl itself.
5d7098d5 1572
58cde26e 1573=head1 AUTHORS
5d7098d5 1574
58cde26e
JH
1575Mark Biggar, overloaded interface by Ilya Zakharevich.
1576Completely rewritten by Tels http://bloodgate.com in 2001.
a5f75d66 1577
a5f75d66 1578=cut