ad070096c701dce59c4e616c349634dd825b13b7
[perl.git] / dist / Math-BigRat / lib / Math / BigRat.pm
1
2 #
3 # "Tax the rat farms." - Lord Vetinari
4 #
5
6 # The following hash values are used:
7 #   sign : +,-,NaN,+inf,-inf
8 #   _d   : denominator
9 #   _n   : numerator (value = _n/_d)
10 #   _a   : accuracy
11 #   _p   : precision
12 # You should not look at the innards of a BigRat - use the methods for this.
13
14 package Math::BigRat;
15
16 # anything older is untested, and unlikely to work
17 use 5.006;
18 use strict;
19 use Carp ();
20
21 use Math::BigFloat;
22 use vars qw($VERSION @ISA $upgrade $downgrade
23             $accuracy $precision $round_mode $div_scale $_trap_nan $_trap_inf);
24
25 @ISA = qw(Math::BigFloat);
26
27 $VERSION = '0.26_01';
28 $VERSION = eval $VERSION;
29
30 # inherit overload from Math::BigFloat, but disable the bitwise ops that don't
31 # make much sense for rationals unless they're truncated or something first
32
33 use overload
34     map {
35         my $op = $_;
36         ($op => sub {
37             Carp::croak("bitwise operation $op not supported in Math::BigRat");
38         });
39     } qw(& | ^ ~ << >> &= |= ^= <<= >>=);
40
41 BEGIN
42   {
43   *objectify = \&Math::BigInt::objectify;       # inherit this from BigInt
44   *AUTOLOAD = \&Math::BigFloat::AUTOLOAD;       # can't inherit AUTOLOAD
45   # we inherit these from BigFloat because currently it is not possible
46   # that MBF has a different $MBI variable than we, because MBF also uses
47   # Math::BigInt::config->('lib'); (there is always only one library loaded)
48   *_e_add = \&Math::BigFloat::_e_add;
49   *_e_sub = \&Math::BigFloat::_e_sub;
50   *as_int = \&as_number;
51   *is_pos = \&is_positive;
52   *is_neg = \&is_negative;
53   }
54
55 ##############################################################################
56 # Global constants and flags. Access these only via the accessor methods!
57
58 $accuracy = $precision = undef;
59 $round_mode = 'even';
60 $div_scale = 40;
61 $upgrade = undef;
62 $downgrade = undef;
63
64 # These are internally, and not to be used from the outside at all!
65
66 $_trap_nan = 0;                         # are NaNs ok? set w/ config()
67 $_trap_inf = 0;                         # are infs ok? set w/ config()
68
69 # the package we are using for our private parts, defaults to:
70 # Math::BigInt->config()->{lib}
71 my $MBI = 'Math::BigInt::Calc';
72
73 my $nan = 'NaN';
74 my $class = 'Math::BigRat';
75
76 sub isa
77   {
78   return 0 if $_[1] =~ /^Math::Big(Int|Float)/;         # we aren't
79   UNIVERSAL::isa(@_);
80   }
81
82 ##############################################################################
83
84 sub _new_from_float
85   {
86   # turn a single float input into a rational number (like '0.1')
87   my ($self,$f) = @_;
88
89   return $self->bnan() if $f->is_nan();
90   return $self->binf($f->{sign}) if $f->{sign} =~ /^[+-]inf$/;
91
92   $self->{_n} = $MBI->_copy( $f->{_m} );        # mantissa
93   $self->{_d} = $MBI->_one();
94   $self->{sign} = $f->{sign} || '+';
95   if ($f->{_es} eq '-')
96     {
97     # something like Math::BigRat->new('0.1');
98     # 1 / 1 => 1/10
99     $MBI->_lsft ( $self->{_d}, $f->{_e} ,10);
100     }
101   else
102     {
103     # something like Math::BigRat->new('10');
104     # 1 / 1 => 10/1
105     $MBI->_lsft ( $self->{_n}, $f->{_e} ,10) unless
106       $MBI->_is_zero($f->{_e});
107     }
108   $self;
109   }
110
111 sub new
112   {
113   # create a Math::BigRat
114   my $class = shift;
115
116   my ($n,$d) = @_;
117
118   my $self = { }; bless $self,$class;
119
120   # input like (BigInt) or (BigFloat):
121   if ((!defined $d) && (ref $n) && (!$n->isa('Math::BigRat')))
122     {
123     if ($n->isa('Math::BigFloat'))
124       {
125       $self->_new_from_float($n);
126       }
127     if ($n->isa('Math::BigInt'))
128       {
129       # TODO: trap NaN, inf
130       $self->{_n} = $MBI->_copy($n->{value});           # "mantissa" = N
131       $self->{_d} = $MBI->_one();                       # d => 1
132       $self->{sign} = $n->{sign};
133       }
134     if ($n->isa('Math::BigInt::Lite'))
135       {
136       # TODO: trap NaN, inf
137       $self->{sign} = '+'; $self->{sign} = '-' if $$n < 0;
138       $self->{_n} = $MBI->_new(abs($$n));               # "mantissa" = N
139       $self->{_d} = $MBI->_one();                       # d => 1
140       }
141     return $self->bnorm();                              # normalize (120/1 => 12/10)
142     }
143
144   # input like (BigInt,BigInt) or (BigLite,BigLite):
145   if (ref($d) && ref($n))
146     {
147     # do N first (for $self->{sign}):
148     if ($n->isa('Math::BigInt'))
149       {
150       # TODO: trap NaN, inf
151       $self->{_n} = $MBI->_copy($n->{value});           # "mantissa" = N
152       $self->{sign} = $n->{sign};
153       }
154     elsif ($n->isa('Math::BigInt::Lite'))
155       {
156       # TODO: trap NaN, inf
157       $self->{sign} = '+'; $self->{sign} = '-' if $$n < 0;
158       $self->{_n} = $MBI->_new(abs($$n));               # "mantissa" = $n
159       }
160     else
161       {
162       require Carp;
163       Carp::croak(ref($n) . " is not a recognized object format for Math::BigRat->new");
164       }
165     # now D:
166     if ($d->isa('Math::BigInt'))
167       {
168       # TODO: trap NaN, inf
169       $self->{_d} = $MBI->_copy($d->{value});           # "mantissa" = D
170       # +/+ or -/- => +, +/- or -/+ => -
171       $self->{sign} = $d->{sign} ne $self->{sign} ? '-' : '+';
172       }
173     elsif ($d->isa('Math::BigInt::Lite'))
174       {
175       # TODO: trap NaN, inf
176       $self->{_d} = $MBI->_new(abs($$d));               # "mantissa" = D
177       my $ds = '+'; $ds = '-' if $$d < 0;
178       # +/+ or -/- => +, +/- or -/+ => -
179       $self->{sign} = $ds ne $self->{sign} ? '-' : '+';
180       }
181     else
182       {
183       require Carp;
184       Carp::croak(ref($d) . " is not a recognized object format for Math::BigRat->new");
185       }
186     return $self->bnorm();                              # normalize (120/1 => 12/10)
187     }
188   return $n->copy() if ref $n;                          # already a BigRat
189
190   if (!defined $n)
191     {
192     $self->{_n} = $MBI->_zero();                        # undef => 0
193     $self->{_d} = $MBI->_one();
194     $self->{sign} = '+';
195     return $self;
196     }
197
198   # string input with / delimiter
199   if ($n =~ /\s*\/\s*/)
200     {
201     return $class->bnan() if $n =~ /\/.*\//;    # 1/2/3 isn't valid
202     return $class->bnan() if $n =~ /\/\s*$/;    # 1/ isn't valid
203     ($n,$d) = split (/\//,$n);
204     # try as BigFloats first
205     if (($n =~ /[\.eE]/) || ($d =~ /[\.eE]/))
206       {
207       local $Math::BigFloat::accuracy = undef;
208       local $Math::BigFloat::precision = undef;
209
210       # one of them looks like a float
211       my $nf = Math::BigFloat->new($n,undef,undef);
212       $self->{sign} = '+';
213       return $self->bnan() if $nf->is_nan();
214
215       $self->{_n} = $MBI->_copy( $nf->{_m} );   # get mantissa
216
217       # now correct $self->{_n} due to $n
218       my $f = Math::BigFloat->new($d,undef,undef);
219       return $self->bnan() if $f->is_nan();
220       $self->{_d} = $MBI->_copy( $f->{_m} );
221
222       # calculate the difference between nE and dE
223       my $diff_e = $nf->exponent()->bsub( $f->exponent);
224       if ($diff_e->is_negative())
225         {
226         # < 0: mul d with it
227         $MBI->_lsft( $self->{_d}, $MBI->_new( $diff_e->babs()), 10);
228         }
229       elsif (!$diff_e->is_zero())
230         {
231         # > 0: mul n with it
232         $MBI->_lsft( $self->{_n}, $MBI->_new( $diff_e), 10);
233         }
234       }
235     else
236       {
237       # both d and n look like (big)ints
238
239       $self->{sign} = '+';                                      # no sign => '+'
240       $self->{_n} = undef;
241       $self->{_d} = undef;
242       if ($n =~ /^([+-]?)0*([0-9]+)\z/)                         # first part ok?
243         {
244         $self->{sign} = $1 || '+';                              # no sign => '+'
245         $self->{_n} = $MBI->_new($2 || 0);
246         }
247
248       if ($d =~ /^([+-]?)0*([0-9]+)\z/)                         # second part ok?
249         {
250         $self->{sign} =~ tr/+-/-+/ if ($1 || '') eq '-';        # negate if second part neg.
251         $self->{_d} = $MBI->_new($2 || 0);
252         }
253
254       if (!defined $self->{_n} || !defined $self->{_d})
255         {
256         $d = Math::BigInt->new($d,undef,undef) unless ref $d;
257         $n = Math::BigInt->new($n,undef,undef) unless ref $n;
258
259         if ($n->{sign} =~ /^[+-]$/ && $d->{sign} =~ /^[+-]$/)
260           {
261           # both parts are ok as integers (wierd things like ' 1e0'
262           $self->{_n} = $MBI->_copy($n->{value});
263           $self->{_d} = $MBI->_copy($d->{value});
264           $self->{sign} = $n->{sign};
265           $self->{sign} =~ tr/+-/-+/ if $d->{sign} eq '-';      # -1/-2 => 1/2
266           return $self->bnorm();
267           }
268
269         $self->{sign} = '+';                                    # a default sign
270         return $self->bnan() if $n->is_nan() || $d->is_nan();
271
272         # handle inf cases:
273         if ($n->is_inf() || $d->is_inf())
274           {
275           if ($n->is_inf())
276             {
277             return $self->bnan() if $d->is_inf();               # both are inf => NaN
278             my $s = '+';                # '+inf/+123' or '-inf/-123'
279             $s = '-' if substr($n->{sign},0,1) ne $d->{sign};
280             # +-inf/123 => +-inf
281             return $self->binf($s);
282             }
283           # 123/inf => 0
284           return $self->bzero();
285           }
286         }
287       }
288
289     return $self->bnorm();
290     }
291
292   # simple string input
293   if (($n =~ /[\.eE]/) && $n !~ /^0x/)
294     {
295     # looks like a float, quacks like a float, so probably is a float
296     $self->{sign} = 'NaN';
297     local $Math::BigFloat::accuracy = undef;
298     local $Math::BigFloat::precision = undef;
299     $self->_new_from_float(Math::BigFloat->new($n,undef,undef));
300     }
301   else
302     {
303     # for simple forms, use $MBI directly
304     if ($n =~ /^([+-]?)0*([0-9]+)\z/)
305       {
306       $self->{sign} = $1 || '+';
307       $self->{_n} = $MBI->_new($2 || 0);
308       $self->{_d} = $MBI->_one();
309       }
310     else
311       {
312       my $n = Math::BigInt->new($n,undef,undef);
313       $self->{_n} = $MBI->_copy($n->{value});
314       $self->{_d} = $MBI->_one();
315       $self->{sign} = $n->{sign};
316       return $self->bnan() if $self->{sign} eq 'NaN';
317       return $self->binf($self->{sign}) if $self->{sign} =~ /^[+-]inf$/;
318       }
319     }
320   $self->bnorm();
321   }
322
323 sub copy
324   {
325   # if two arguments, the first one is the class to "swallow" subclasses
326   my ($c,$x) = @_;
327
328   if (scalar @_ == 1)
329     {
330     $x = $_[0];
331     $c = ref($x);
332     }
333   return unless ref($x); # only for objects
334
335   my $self = bless {}, $c;
336
337   $self->{sign} = $x->{sign};
338   $self->{_d} = $MBI->_copy($x->{_d});
339   $self->{_n} = $MBI->_copy($x->{_n});
340   $self->{_a} = $x->{_a} if defined $x->{_a};
341   $self->{_p} = $x->{_p} if defined $x->{_p};
342   $self;
343   }
344
345 ##############################################################################
346
347 sub config
348   {
349   # return (later set?) configuration data as hash ref
350   my $class = shift || 'Math::BigRat';
351
352   if (@_ == 1 && ref($_[0]) ne 'HASH')
353     {
354     my $cfg = $class->SUPER::config();
355     return $cfg->{$_[0]};
356     }
357
358   my $cfg = $class->SUPER::config(@_);
359
360   # now we need only to override the ones that are different from our parent
361   $cfg->{class} = $class;
362   $cfg->{with} = $MBI;
363   $cfg;
364   }
365
366 ##############################################################################
367
368 sub bstr
369   {
370   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
371
372   if ($x->{sign} !~ /^[+-]$/)           # inf, NaN etc
373     {
374     my $s = $x->{sign}; $s =~ s/^\+//;  # +inf => inf
375     return $s;
376     }
377
378   my $s = ''; $s = $x->{sign} if $x->{sign} ne '+';     # '+3/2' => '3/2'
379
380   return $s . $MBI->_str($x->{_n}) if $MBI->_is_one($x->{_d});
381   $s . $MBI->_str($x->{_n}) . '/' . $MBI->_str($x->{_d});
382   }
383
384 sub bsstr
385   {
386   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
387
388   if ($x->{sign} !~ /^[+-]$/)           # inf, NaN etc
389     {
390     my $s = $x->{sign}; $s =~ s/^\+//;  # +inf => inf
391     return $s;
392     }
393
394   my $s = ''; $s = $x->{sign} if $x->{sign} ne '+';     # +3 vs 3
395   $s . $MBI->_str($x->{_n}) . '/' . $MBI->_str($x->{_d});
396   }
397
398 sub bnorm
399   {
400   # reduce the number to the shortest form
401   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
402
403   # Both parts must be objects of whatever we are using today.
404   if ( my $c = $MBI->_check($x->{_n}) )
405     {
406     require Carp; Carp::croak ("n did not pass the self-check ($c) in bnorm()");
407     }
408   if ( my $c = $MBI->_check($x->{_d}) )
409     {
410     require Carp; Carp::croak ("d did not pass the self-check ($c) in bnorm()");
411     }
412
413   # no normalize for NaN, inf etc.
414   return $x if $x->{sign} !~ /^[+-]$/;
415
416   # normalize zeros to 0/1
417   if ($MBI->_is_zero($x->{_n}))
418     {
419     $x->{sign} = '+';                                   # never leave a -0
420     $x->{_d} = $MBI->_one() unless $MBI->_is_one($x->{_d});
421     return $x;
422     }
423
424   return $x if $MBI->_is_one($x->{_d});                 # no need to reduce
425
426   # reduce other numbers
427   my $gcd = $MBI->_copy($x->{_n});
428   $gcd = $MBI->_gcd($gcd,$x->{_d});
429
430   if (!$MBI->_is_one($gcd))
431     {
432     $x->{_n} = $MBI->_div($x->{_n},$gcd);
433     $x->{_d} = $MBI->_div($x->{_d},$gcd);
434     }
435   $x;
436   }
437
438 ##############################################################################
439 # sign manipulation
440
441 sub bneg
442   {
443   # (BRAT or num_str) return BRAT
444   # negate number or make a negated number from string
445   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
446
447   return $x if $x->modify('bneg');
448
449   # for +0 dont negate (to have always normalized +0). Does nothing for 'NaN'
450   $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $MBI->_is_zero($x->{_n}));
451   $x;
452   }
453
454 ##############################################################################
455 # special values
456
457 sub _bnan
458   {
459   # used by parent class bnan() to initialize number to NaN
460   my $self = shift;
461
462   if ($_trap_nan)
463     {
464     require Carp;
465     my $class = ref($self);
466     # "$self" below will stringify the object, this blows up if $self is a
467     # partial object (happens under trap_nan), so fix it beforehand
468     $self->{_d} = $MBI->_zero() unless defined $self->{_d};
469     $self->{_n} = $MBI->_zero() unless defined $self->{_n};
470     Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
471     }
472   $self->{_n} = $MBI->_zero();
473   $self->{_d} = $MBI->_zero();
474   }
475
476 sub _binf
477   {
478   # used by parent class bone() to initialize number to +inf/-inf
479   my $self = shift;
480
481   if ($_trap_inf)
482     {
483     require Carp;
484     my $class = ref($self);
485     # "$self" below will stringify the object, this blows up if $self is a
486     # partial object (happens under trap_nan), so fix it beforehand
487     $self->{_d} = $MBI->_zero() unless defined $self->{_d};
488     $self->{_n} = $MBI->_zero() unless defined $self->{_n};
489     Carp::croak ("Tried to set $self to inf in $class\::_binf()");
490     }
491   $self->{_n} = $MBI->_zero();
492   $self->{_d} = $MBI->_zero();
493   }
494
495 sub _bone
496   {
497   # used by parent class bone() to initialize number to +1/-1
498   my $self = shift;
499   $self->{_n} = $MBI->_one();
500   $self->{_d} = $MBI->_one();
501   }
502
503 sub _bzero
504   {
505   # used by parent class bzero() to initialize number to 0
506   my $self = shift;
507   $self->{_n} = $MBI->_zero();
508   $self->{_d} = $MBI->_one();
509   }
510
511 ##############################################################################
512 # mul/add/div etc
513
514 sub badd
515   {
516   # add two rational numbers
517
518   # set up parameters
519   my ($self,$x,$y,@r) = (ref($_[0]),@_);
520   # objectify is costly, so avoid it
521   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
522     {
523     ($self,$x,$y,@r) = objectify(2,@_);
524     }
525
526   # +inf + +inf => +inf,  -inf + -inf => -inf
527   return $x->binf(substr($x->{sign},0,1))
528     if $x->{sign} eq $y->{sign} && $x->{sign} =~ /^[+-]inf$/;
529
530   # +inf + -inf or -inf + +inf => NaN
531   return $x->bnan() if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/);
532
533   #  1   1    gcd(3,4) = 1    1*3 + 1*4    7
534   #  - + -                  = --------- = --
535   #  4   3                      4*3       12
536
537   # we do not compute the gcd() here, but simple do:
538   #  5   7    5*3 + 7*4   43
539   #  - + -  = --------- = --
540   #  4   3       4*3      12
541
542   # and bnorm() will then take care of the rest
543
544   # 5 * 3
545   $x->{_n} = $MBI->_mul( $x->{_n}, $y->{_d});
546
547   # 7 * 4
548   my $m = $MBI->_mul( $MBI->_copy( $y->{_n} ), $x->{_d} );
549
550   # 5 * 3 + 7 * 4
551   ($x->{_n}, $x->{sign}) = _e_add( $x->{_n}, $m, $x->{sign}, $y->{sign});
552
553   # 4 * 3
554   $x->{_d} = $MBI->_mul( $x->{_d}, $y->{_d});
555
556   # normalize result, and possible round
557   $x->bnorm()->round(@r);
558   }
559
560 sub bsub
561   {
562   # subtract two rational numbers
563
564   # set up parameters
565   my ($self,$x,$y,@r) = (ref($_[0]),@_);
566   # objectify is costly, so avoid it
567   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
568     {
569     ($self,$x,$y,@r) = objectify(2,@_);
570     }
571
572   # flip sign of $x, call badd(), then flip sign of result
573   $x->{sign} =~ tr/+-/-+/
574     unless $x->{sign} eq '+' && $MBI->_is_zero($x->{_n});       # not -0
575   $x->badd($y,@r);                              # does norm and round
576   $x->{sign} =~ tr/+-/-+/
577     unless $x->{sign} eq '+' && $MBI->_is_zero($x->{_n});       # not -0
578   $x;
579   }
580
581 sub bmul
582   {
583   # multiply two rational numbers
584
585   # set up parameters
586   my ($self,$x,$y,@r) = (ref($_[0]),@_);
587   # objectify is costly, so avoid it
588   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
589     {
590     ($self,$x,$y,@r) = objectify(2,@_);
591     }
592
593   return $x->bnan() if ($x->{sign} eq 'NaN' || $y->{sign} eq 'NaN');
594
595   # inf handling
596   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
597     {
598     return $x->bnan() if $x->is_zero() || $y->is_zero();
599     # result will always be +-inf:
600     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
601     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
602     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
603     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
604     return $x->binf('-');
605     }
606
607   # x== 0 # also: or y == 1 or y == -1
608   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
609
610   # XXX TODO:
611   # According to Knuth, this can be optimized by doing gcd twice (for d and n)
612   # and reducing in one step. This would save us the bnorm() at the end.
613
614   #  1   2    1 * 2    2    1
615   #  - * - =  -----  = -  = -
616   #  4   3    4 * 3    12   6
617
618   $x->{_n} = $MBI->_mul( $x->{_n}, $y->{_n});
619   $x->{_d} = $MBI->_mul( $x->{_d}, $y->{_d});
620
621   # compute new sign
622   $x->{sign} = $x->{sign} eq $y->{sign} ? '+' : '-';
623
624   $x->bnorm()->round(@r);
625   }
626
627 sub bdiv
628   {
629   # (dividend: BRAT or num_str, divisor: BRAT or num_str) return
630   # (BRAT,BRAT) (quo,rem) or BRAT (only rem)
631
632   # set up parameters
633   my ($self,$x,$y,@r) = (ref($_[0]),@_);
634   # objectify is costly, so avoid it
635   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
636     {
637     ($self,$x,$y,@r) = objectify(2,@_);
638     }
639
640   return $self->_div_inf($x,$y)
641    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
642
643   # x== 0 # also: or y == 1 or y == -1
644   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
645
646   # XXX TODO: list context, upgrade
647   # According to Knuth, this can be optimized by doing gcd twice (for d and n)
648   # and reducing in one step. This would save us the bnorm() at the end.
649
650   # 1     1    1   3
651   # -  /  - == - * -
652   # 4     3    4   1
653
654   $x->{_n} = $MBI->_mul( $x->{_n}, $y->{_d});
655   $x->{_d} = $MBI->_mul( $x->{_d}, $y->{_n});
656
657   # compute new sign
658   $x->{sign} = $x->{sign} eq $y->{sign} ? '+' : '-';
659
660   $x->bnorm()->round(@r);
661   $x;
662   }
663
664 sub bmod
665   {
666   # compute "remainder" (in Perl way) of $x / $y
667
668   # set up parameters
669   my ($self,$x,$y,@r) = (ref($_[0]),@_);
670   # objectify is costly, so avoid it
671   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
672     {
673     ($self,$x,$y,@r) = objectify(2,@_);
674     }
675
676   return $self->_div_inf($x,$y)
677    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
678
679   return $x if $x->is_zero();           # 0 / 7 = 0, mod 0
680
681   # compute $x - $y * floor($x/$y), keeping the sign of $x
682
683   # copy x to u, make it positive and then do a normal division ($u/$y)
684   my $u = bless { sign => '+' }, $self;
685   $u->{_n} = $MBI->_mul( $MBI->_copy($x->{_n}), $y->{_d} );
686   $u->{_d} = $MBI->_mul( $MBI->_copy($x->{_d}), $y->{_n} );
687
688   # compute floor(u)
689   if (! $MBI->_is_one($u->{_d}))
690     {
691     $u->{_n} = $MBI->_div($u->{_n},$u->{_d});   # 22/7 => 3/1 w/ truncate
692     # no need to set $u->{_d} to 1, since below we set it to $y->{_d} anyway
693     }
694
695   # now compute $y * $u
696   $u->{_d} = $MBI->_copy($y->{_d});             # 1 * $y->{_d}, see floor above
697   $u->{_n} = $MBI->_mul($u->{_n},$y->{_n});
698
699   my $xsign = $x->{sign}; $x->{sign} = '+';     # remember sign and make x positive
700   # compute $x - $u
701   $x->bsub($u);
702   $x->{sign} = $xsign;                          # put sign back
703
704   $x->bnorm()->round(@r);
705   }
706
707 ##############################################################################
708 # bdec/binc
709
710 sub bdec
711   {
712   # decrement value (subtract 1)
713   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
714
715   return $x if $x->{sign} !~ /^[+-]$/;  # NaN, inf, -inf
716
717   if ($x->{sign} eq '-')
718     {
719     $x->{_n} = $MBI->_add( $x->{_n}, $x->{_d});         # -5/2 => -7/2
720     }
721   else
722     {
723     if ($MBI->_acmp($x->{_n},$x->{_d}) < 0)             # n < d?
724       {
725       # 1/3 -- => -2/3
726       $x->{_n} = $MBI->_sub( $MBI->_copy($x->{_d}), $x->{_n});
727       $x->{sign} = '-';
728       }
729     else
730       {
731       $x->{_n} = $MBI->_sub($x->{_n}, $x->{_d});        # 5/2 => 3/2
732       }
733     }
734   $x->bnorm()->round(@r);
735   }
736
737 sub binc
738   {
739   # increment value (add 1)
740   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
741
742   return $x if $x->{sign} !~ /^[+-]$/;  # NaN, inf, -inf
743
744   if ($x->{sign} eq '-')
745     {
746     if ($MBI->_acmp($x->{_n},$x->{_d}) < 0)
747       {
748       # -1/3 ++ => 2/3 (overflow at 0)
749       $x->{_n} = $MBI->_sub( $MBI->_copy($x->{_d}), $x->{_n});
750       $x->{sign} = '+';
751       }
752     else
753       {
754       $x->{_n} = $MBI->_sub($x->{_n}, $x->{_d});        # -5/2 => -3/2
755       }
756     }
757   else
758     {
759     $x->{_n} = $MBI->_add($x->{_n},$x->{_d});           # 5/2 => 7/2
760     }
761   $x->bnorm()->round(@r);
762   }
763
764 ##############################################################################
765 # is_foo methods (the rest is inherited)
766
767 sub is_int
768   {
769   # return true if arg (BRAT or num_str) is an integer
770   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
771
772   return 1 if ($x->{sign} =~ /^[+-]$/) &&       # NaN and +-inf aren't
773     $MBI->_is_one($x->{_d});                    # x/y && y != 1 => no integer
774   0;
775   }
776
777 sub is_zero
778   {
779   # return true if arg (BRAT or num_str) is zero
780   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
781
782   return 1 if $x->{sign} eq '+' && $MBI->_is_zero($x->{_n});
783   0;
784   }
785
786 sub is_one
787   {
788   # return true if arg (BRAT or num_str) is +1 or -1 if signis given
789   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
790
791   my $sign = $_[2] || ''; $sign = '+' if $sign ne '-';
792   return 1
793    if ($x->{sign} eq $sign && $MBI->_is_one($x->{_n}) && $MBI->_is_one($x->{_d}));
794   0;
795   }
796
797 sub is_odd
798   {
799   # return true if arg (BFLOAT or num_str) is odd or false if even
800   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
801
802   return 1 if ($x->{sign} =~ /^[+-]$/) &&               # NaN & +-inf aren't
803     ($MBI->_is_one($x->{_d}) && $MBI->_is_odd($x->{_n})); # x/2 is not, but 3/1
804   0;
805   }
806
807 sub is_even
808   {
809   # return true if arg (BINT or num_str) is even or false if odd
810   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
811
812   return 0 if $x->{sign} !~ /^[+-]$/;                   # NaN & +-inf aren't
813   return 1 if ($MBI->_is_one($x->{_d})                  # x/3 is never
814      && $MBI->_is_even($x->{_n}));                      # but 4/1 is
815   0;
816   }
817
818 ##############################################################################
819 # parts() and friends
820
821 sub numerator
822   {
823   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
824
825   # NaN, inf, -inf
826   return Math::BigInt->new($x->{sign}) if ($x->{sign} !~ /^[+-]$/);
827
828   my $n = Math::BigInt->new($MBI->_str($x->{_n})); $n->{sign} = $x->{sign};
829   $n;
830   }
831
832 sub denominator
833   {
834   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
835
836   # NaN
837   return Math::BigInt->new($x->{sign}) if $x->{sign} eq 'NaN';
838   # inf, -inf
839   return Math::BigInt->bone() if $x->{sign} !~ /^[+-]$/;
840
841   Math::BigInt->new($MBI->_str($x->{_d}));
842   }
843
844 sub parts
845   {
846   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
847
848   my $c = 'Math::BigInt';
849
850   return ($c->bnan(),$c->bnan()) if $x->{sign} eq 'NaN';
851   return ($c->binf(),$c->binf()) if $x->{sign} eq '+inf';
852   return ($c->binf('-'),$c->binf()) if $x->{sign} eq '-inf';
853
854   my $n = $c->new( $MBI->_str($x->{_n}));
855   $n->{sign} = $x->{sign};
856   my $d = $c->new( $MBI->_str($x->{_d}));
857   ($n,$d);
858   }
859
860 sub length
861   {
862   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
863
864   return $nan unless $x->is_int();
865   $MBI->_len($x->{_n});                         # length(-123/1) => length(123)
866   }
867
868 sub digit
869   {
870   my ($self,$x,$n) = ref($_[0]) ? (undef,$_[0],$_[1]) : objectify(1,@_);
871
872   return $nan unless $x->is_int();
873   $MBI->_digit($x->{_n},$n || 0);               # digit(-123/1,2) => digit(123,2)
874   }
875
876 ##############################################################################
877 # special calc routines
878
879 sub bceil
880   {
881   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
882
883   return $x if $x->{sign} !~ /^[+-]$/ ||        # not for NaN, inf
884             $MBI->_is_one($x->{_d});            # 22/1 => 22, 0/1 => 0
885
886   $x->{_n} = $MBI->_div($x->{_n},$x->{_d});     # 22/7 => 3/1 w/ truncate
887   $x->{_d} = $MBI->_one();                      # d => 1
888   $x->{_n} = $MBI->_inc($x->{_n})
889     if $x->{sign} eq '+';                       # +22/7 => 4/1
890   $x->{sign} = '+' if $MBI->_is_zero($x->{_n}); # -0 => 0
891   $x;
892   }
893
894 sub bfloor
895   {
896   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
897
898   return $x if $x->{sign} !~ /^[+-]$/ ||        # not for NaN, inf
899             $MBI->_is_one($x->{_d});            # 22/1 => 22, 0/1 => 0
900
901   $x->{_n} = $MBI->_div($x->{_n},$x->{_d});     # 22/7 => 3/1 w/ truncate
902   $x->{_d} = $MBI->_one();                      # d => 1
903   $x->{_n} = $MBI->_inc($x->{_n})
904     if $x->{sign} eq '-';                       # -22/7 => -4/1
905   $x;
906   }
907
908 sub bfac
909   {
910   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
911
912   # if $x is not an integer
913   if (($x->{sign} ne '+') || (!$MBI->_is_one($x->{_d})))
914     {
915     return $x->bnan();
916     }
917
918   $x->{_n} = $MBI->_fac($x->{_n});
919   # since _d is 1, we don't need to reduce/norm the result
920   $x->round(@r);
921   }
922
923 sub bpow
924   {
925   # power ($x ** $y)
926
927   # set up parameters
928   my ($self,$x,$y,@r) = (ref($_[0]),@_);
929   # objectify is costly, so avoid it
930   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
931     {
932     ($self,$x,$y,@r) = objectify(2,@_);
933     }
934
935   return $x if $x->{sign} =~ /^[+-]inf$/;       # -inf/+inf ** x
936   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
937   return $x->bone(@r) if $y->is_zero();
938   return $x->round(@r) if $x->is_one() || $y->is_one();
939
940   if ($x->{sign} eq '-' && $MBI->_is_one($x->{_n}) && $MBI->_is_one($x->{_d}))
941     {
942     # if $x == -1 and odd/even y => +1/-1
943     return $y->is_odd() ? $x->round(@r) : $x->babs()->round(@r);
944     # my Casio FX-5500L has a bug here: -1 ** 2 is -1, but -1 * -1 is 1;
945     }
946   # 1 ** -y => 1 / (1 ** |y|)
947   # so do test for negative $y after above's clause
948
949   return $x->round(@r) if $x->is_zero();  # 0**y => 0 (if not y <= 0)
950
951   # shortcut if y == 1/N (is then sqrt() respective broot())
952   if ($MBI->_is_one($y->{_n}))
953     {
954     return $x->bsqrt(@r) if $MBI->_is_two($y->{_d});    # 1/2 => sqrt
955     return $x->broot($MBI->_str($y->{_d}),@r);          # 1/N => root(N)
956     }
957
958   # shortcut y/1 (and/or x/1)
959   if ($MBI->_is_one($y->{_d}))
960     {
961     # shortcut for x/1 and y/1
962     if ($MBI->_is_one($x->{_d}))
963       {
964       $x->{_n} = $MBI->_pow($x->{_n},$y->{_n});         # x/1 ** y/1 => (x ** y)/1
965       if ($y->{sign} eq '-')
966         {
967         # 0.2 ** -3 => 1/(0.2 ** 3)
968         ($x->{_n},$x->{_d}) = ($x->{_d},$x->{_n});      # swap
969         }
970       # correct sign; + ** + => +
971       if ($x->{sign} eq '-')
972         {
973         # - * - => +, - * - * - => -
974         $x->{sign} = '+' if $MBI->_is_even($y->{_n});
975         }
976       return $x->round(@r);
977       }
978     # x/z ** y/1
979     $x->{_n} = $MBI->_pow($x->{_n},$y->{_n});           # 5/2 ** y/1 => 5 ** y / 2 ** y
980     $x->{_d} = $MBI->_pow($x->{_d},$y->{_n});
981     if ($y->{sign} eq '-')
982       {
983       # 0.2 ** -3 => 1/(0.2 ** 3)
984       ($x->{_n},$x->{_d}) = ($x->{_d},$x->{_n});        # swap
985       }
986     # correct sign; + ** + => +
987     if ($x->{sign} eq '-')
988       {
989       # - * - => +, - * - * - => -
990       $x->{sign} = '+' if $MBI->_is_even($y->{_n});
991       }
992     return $x->round(@r);
993     }
994
995 #  print STDERR "# $x $y\n";
996
997   # otherwise:
998
999   #      n/d     n  ______________
1000   # a/b       =  -\/  (a/b) ** d
1001
1002   # (a/b) ** n == (a ** n) / (b ** n)
1003   $MBI->_pow($x->{_n}, $y->{_n} );
1004   $MBI->_pow($x->{_d}, $y->{_n} );
1005
1006   return $x->broot($MBI->_str($y->{_d}),@r);            # n/d => root(n)
1007   }
1008
1009 sub blog
1010   {
1011   # set up parameters
1012   my ($self,$x,$y,@r) = (ref($_[0]),@_);
1013
1014   # objectify is costly, so avoid it
1015   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1016     {
1017     ($self,$x,$y,@r) = objectify(2,$class,@_);
1018     }
1019
1020   # blog(1,Y) => 0
1021   return $x->bzero() if $x->is_one() && $y->{sign} eq '+';
1022
1023   # $x <= 0 => NaN
1024   return $x->bnan() if $x->is_zero() || $x->{sign} ne '+' || $y->{sign} ne '+';
1025
1026   if ($x->is_int() && $y->is_int())
1027     {
1028     return $self->new($x->as_number()->blog($y->as_number(),@r));
1029     }
1030
1031   # do it with floats
1032   $x->_new_from_float( $x->_as_float()->blog(Math::BigFloat->new("$y"),@r) );
1033   }
1034
1035 sub bexp
1036   {
1037   # set up parameters
1038   my ($self,$x,$y,@r) = (ref($_[0]),@_);
1039
1040   # objectify is costly, so avoid it
1041   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1042     {
1043     ($self,$x,$y,@r) = objectify(2,$class,@_);
1044     }
1045
1046   return $x->binf(@r) if $x->{sign} eq '+inf';
1047   return $x->bzero(@r) if $x->{sign} eq '-inf';
1048
1049   # we need to limit the accuracy to protect against overflow
1050   my $fallback = 0;
1051   my ($scale,@params);
1052   ($x,@params) = $x->_find_round_parameters(@r);
1053
1054   # also takes care of the "error in _find_round_parameters?" case
1055   return $x if $x->{sign} eq 'NaN';
1056
1057   # no rounding at all, so must use fallback
1058   if (scalar @params == 0)
1059     {
1060     # simulate old behaviour
1061     $params[0] = $self->div_scale();    # and round to it as accuracy
1062     $params[1] = undef;                 # P = undef
1063     $scale = $params[0]+4;              # at least four more for proper round
1064     $params[2] = $r[2];                 # round mode by caller or undef
1065     $fallback = 1;                      # to clear a/p afterwards
1066     }
1067   else
1068     {
1069     # the 4 below is empirical, and there might be cases where it's not enough...
1070     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1071     }
1072
1073   return $x->bone(@params) if $x->is_zero();
1074
1075   # See the comments in Math::BigFloat on how this algorithm works.
1076   # Basically we calculate A and B (where B is faculty(N)) so that A/B = e
1077
1078   my $x_org = $x->copy();
1079   if ($scale <= 75)
1080     {
1081     # set $x directly from a cached string form
1082     $x->{_n} = $MBI->_new("90933395208605785401971970164779391644753259799242");
1083     $x->{_d} = $MBI->_new("33452526613163807108170062053440751665152000000000");
1084     $x->{sign} = '+';
1085     }
1086   else
1087     {
1088     # compute A and B so that e = A / B.
1089
1090     # After some terms we end up with this, so we use it as a starting point:
1091     my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
1092     my $F = $MBI->_new(42); my $step = 42;
1093
1094     # Compute how many steps we need to take to get $A and $B sufficiently big
1095     my $steps = Math::BigFloat::_len_to_steps($scale - 4);
1096 #    print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
1097     while ($step++ <= $steps)
1098       {
1099       # calculate $a * $f + 1
1100       $A = $MBI->_mul($A, $F);
1101       $A = $MBI->_inc($A);
1102       # increment f
1103       $F = $MBI->_inc($F);
1104       }
1105     # compute $B as factorial of $steps (this is faster than doing it manually)
1106     my $B = $MBI->_fac($MBI->_new($steps));
1107
1108 #  print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";
1109
1110     $x->{_n} = $A;
1111     $x->{_d} = $B;
1112     $x->{sign} = '+';
1113     }
1114
1115   # $x contains now an estimate of e, with some surplus digits, so we can round
1116   if (!$x_org->is_one())
1117     {
1118     # raise $x to the wanted power and round it in one step:
1119     $x->bpow($x_org, @params);
1120     }
1121   else
1122     {
1123     # else just round the already computed result
1124     delete $x->{_a}; delete $x->{_p};
1125     # shortcut to not run through _find_round_parameters again
1126     if (defined $params[0])
1127       {
1128       $x->bround($params[0],$params[2]);                # then round accordingly
1129       }
1130     else
1131       {
1132       $x->bfround($params[1],$params[2]);               # then round accordingly
1133       }
1134     }
1135   if ($fallback)
1136     {
1137     # clear a/p after round, since user did not request it
1138     delete $x->{_a}; delete $x->{_p};
1139     }
1140
1141   $x;
1142   }
1143
1144 sub bnok
1145   {
1146   # set up parameters
1147   my ($self,$x,$y,@r) = (ref($_[0]),@_);
1148
1149   # objectify is costly, so avoid it
1150   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1151     {
1152     ($self,$x,$y,@r) = objectify(2,$class,@_);
1153     }
1154
1155   # do it with floats
1156   $x->_new_from_float( $x->_as_float()->bnok(Math::BigFloat->new("$y"),@r) );
1157   }
1158
1159 sub _float_from_part
1160   {
1161   my $x = shift;
1162
1163   my $f = Math::BigFloat->bzero();
1164   $f->{_m} = $MBI->_copy($x);
1165   $f->{_e} = $MBI->_zero();
1166
1167   $f;
1168   }
1169
1170 sub _as_float
1171   {
1172   my $x = shift;
1173
1174   local $Math::BigFloat::upgrade = undef;
1175   local $Math::BigFloat::accuracy = undef;
1176   local $Math::BigFloat::precision = undef;
1177   # 22/7 => 3.142857143..
1178
1179   my $a = $x->accuracy() || 0;
1180   if ($a != 0 || !$MBI->_is_one($x->{_d}))
1181     {
1182     # n/d
1183     return scalar Math::BigFloat->new($x->{sign} . $MBI->_str($x->{_n}))->bdiv( $MBI->_str($x->{_d}), $x->accuracy());
1184     }
1185   # just n
1186   Math::BigFloat->new($x->{sign} . $MBI->_str($x->{_n}));
1187   }
1188
1189 sub broot
1190   {
1191   # set up parameters
1192   my ($self,$x,$y,@r) = (ref($_[0]),@_);
1193   # objectify is costly, so avoid it
1194   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1195     {
1196     ($self,$x,$y,@r) = objectify(2,@_);
1197     }
1198
1199   if ($x->is_int() && $y->is_int())
1200     {
1201     return $self->new($x->as_number()->broot($y->as_number(),@r));
1202     }
1203
1204   # do it with floats
1205   $x->_new_from_float( $x->_as_float()->broot($y->_as_float(),@r) )->bnorm()->bround(@r);
1206   }
1207
1208 sub bmodpow
1209   {
1210   # set up parameters
1211   my ($self,$x,$y,$m,@r) = (ref($_[0]),@_);
1212   # objectify is costly, so avoid it
1213   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1214     {
1215     ($self,$x,$y,$m,@r) = objectify(3,@_);
1216     }
1217
1218   # $x or $y or $m are NaN or +-inf => NaN
1219   return $x->bnan()
1220    if $x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/ ||
1221    $m->{sign} !~ /^[+-]$/;
1222
1223   if ($x->is_int() && $y->is_int() && $m->is_int())
1224     {
1225     return $self->new($x->as_number()->bmodpow($y->as_number(),$m,@r));
1226     }
1227
1228   warn ("bmodpow() not fully implemented");
1229   $x->bnan();
1230   }
1231
1232 sub bmodinv
1233   {
1234   # set up parameters
1235   my ($self,$x,$y,@r) = (ref($_[0]),@_);
1236   # objectify is costly, so avoid it
1237   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1238     {
1239     ($self,$x,$y,@r) = objectify(2,@_);
1240     }
1241
1242   # $x or $y are NaN or +-inf => NaN
1243   return $x->bnan()
1244    if $x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/;
1245
1246   if ($x->is_int() && $y->is_int())
1247     {
1248     return $self->new($x->as_number()->bmodinv($y->as_number(),@r));
1249     }
1250
1251   warn ("bmodinv() not fully implemented");
1252   $x->bnan();
1253   }
1254
1255 sub bsqrt
1256   {
1257   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1258
1259   return $x->bnan() if $x->{sign} !~ /^[+]/;    # NaN, -inf or < 0
1260   return $x if $x->{sign} eq '+inf';            # sqrt(inf) == inf
1261   return $x->round(@r) if $x->is_zero() || $x->is_one();
1262
1263   local $Math::BigFloat::upgrade = undef;
1264   local $Math::BigFloat::downgrade = undef;
1265   local $Math::BigFloat::precision = undef;
1266   local $Math::BigFloat::accuracy = undef;
1267   local $Math::BigInt::upgrade = undef;
1268   local $Math::BigInt::precision = undef;
1269   local $Math::BigInt::accuracy = undef;
1270
1271   $x->{_n} = _float_from_part( $x->{_n} )->bsqrt();
1272   $x->{_d} = _float_from_part( $x->{_d} )->bsqrt();
1273
1274   # XXX TODO: we probably can optimze this:
1275
1276   # if sqrt(D) was not integer
1277   if ($x->{_d}->{_es} ne '+')
1278     {
1279     $x->{_n}->blsft($x->{_d}->exponent()->babs(),10);   # 7.1/4.51 => 7.1/45.1
1280     $x->{_d} = $MBI->_copy( $x->{_d}->{_m} );           # 7.1/45.1 => 71/45.1
1281     }
1282   # if sqrt(N) was not integer
1283   if ($x->{_n}->{_es} ne '+')
1284     {
1285     $x->{_d}->blsft($x->{_n}->exponent()->babs(),10);   # 71/45.1 => 710/45.1
1286     $x->{_n} = $MBI->_copy( $x->{_n}->{_m} );           # 710/45.1 => 710/451
1287     }
1288
1289   # convert parts to $MBI again
1290   $x->{_n} = $MBI->_lsft( $MBI->_copy( $x->{_n}->{_m} ), $x->{_n}->{_e}, 10)
1291     if ref($x->{_n}) ne $MBI && ref($x->{_n}) ne 'ARRAY';
1292   $x->{_d} = $MBI->_lsft( $MBI->_copy( $x->{_d}->{_m} ), $x->{_d}->{_e}, 10)
1293     if ref($x->{_d}) ne $MBI && ref($x->{_d}) ne 'ARRAY';
1294
1295   $x->bnorm()->round(@r);
1296   }
1297
1298 sub blsft
1299   {
1300   my ($self,$x,$y,$b,@r) = objectify(3,@_);
1301
1302   $b = 2 unless defined $b;
1303   $b = $self->new($b) unless ref ($b);
1304   $x->bmul( $b->copy()->bpow($y), @r);
1305   $x;
1306   }
1307
1308 sub brsft
1309   {
1310   my ($self,$x,$y,$b,@r) = objectify(3,@_);
1311
1312   $b = 2 unless defined $b;
1313   $b = $self->new($b) unless ref ($b);
1314   $x->bdiv( $b->copy()->bpow($y), @r);
1315   $x;
1316   }
1317
1318 ##############################################################################
1319 # round
1320
1321 sub round
1322   {
1323   $_[0];
1324   }
1325
1326 sub bround
1327   {
1328   $_[0];
1329   }
1330
1331 sub bfround
1332   {
1333   $_[0];
1334   }
1335
1336 ##############################################################################
1337 # comparing
1338
1339 sub bcmp
1340   {
1341   # compare two signed numbers
1342
1343   # set up parameters
1344   my ($self,$x,$y) = (ref($_[0]),@_);
1345   # objectify is costly, so avoid it
1346   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1347     {
1348     ($self,$x,$y) = objectify(2,@_);
1349     }
1350
1351   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1352     {
1353     # handle +-inf and NaN
1354     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1355     return 0 if $x->{sign} eq $y->{sign} && $x->{sign} =~ /^[+-]inf$/;
1356     return +1 if $x->{sign} eq '+inf';
1357     return -1 if $x->{sign} eq '-inf';
1358     return -1 if $y->{sign} eq '+inf';
1359     return +1;
1360     }
1361   # check sign for speed first
1362   return 1 if $x->{sign} eq '+' && $y->{sign} eq '-';   # does also 0 <=> -y
1363   return -1 if $x->{sign} eq '-' && $y->{sign} eq '+';  # does also -x <=> 0
1364
1365   # shortcut
1366   my $xz = $MBI->_is_zero($x->{_n});
1367   my $yz = $MBI->_is_zero($y->{_n});
1368   return 0 if $xz && $yz;                               # 0 <=> 0
1369   return -1 if $xz && $y->{sign} eq '+';                # 0 <=> +y
1370   return 1 if $yz && $x->{sign} eq '+';                 # +x <=> 0
1371
1372   my $t = $MBI->_mul( $MBI->_copy($x->{_n}), $y->{_d});
1373   my $u = $MBI->_mul( $MBI->_copy($y->{_n}), $x->{_d});
1374
1375   my $cmp = $MBI->_acmp($t,$u);                         # signs are equal
1376   $cmp = -$cmp if $x->{sign} eq '-';                    # both are '-' => reverse
1377   $cmp;
1378   }
1379
1380 sub bacmp
1381   {
1382   # compare two numbers (as unsigned)
1383
1384   # set up parameters
1385   my ($self,$x,$y) = (ref($_[0]),@_);
1386   # objectify is costly, so avoid it
1387   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1388     {
1389     ($self,$x,$y) = objectify(2,$class,@_);
1390     }
1391
1392   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1393     {
1394     # handle +-inf and NaN
1395     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1396     return 0 if $x->{sign} =~ /^[+-]inf$/ && $y->{sign} =~ /^[+-]inf$/;
1397     return 1 if $x->{sign} =~ /^[+-]inf$/ && $y->{sign} !~ /^[+-]inf$/;
1398     return -1;
1399     }
1400
1401   my $t = $MBI->_mul( $MBI->_copy($x->{_n}), $y->{_d});
1402   my $u = $MBI->_mul( $MBI->_copy($y->{_n}), $x->{_d});
1403   $MBI->_acmp($t,$u);                                   # ignore signs
1404   }
1405
1406 ##############################################################################
1407 # output conversation
1408
1409 sub numify
1410   {
1411   # convert 17/8 => float (aka 2.125)
1412   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1413
1414   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, NaN, etc
1415
1416   # N/1 => N
1417   my $neg = ''; $neg = '-' if $x->{sign} eq '-';
1418   return $neg . $MBI->_num($x->{_n}) if $MBI->_is_one($x->{_d});
1419
1420   $x->_as_float()->numify() + 0.0;
1421   }
1422
1423 sub as_number
1424   {
1425   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1426
1427   # NaN, inf etc
1428   return Math::BigInt->new($x->{sign}) if $x->{sign} !~ /^[+-]$/;
1429
1430   my $u = Math::BigInt->bzero();
1431   $u->{sign} = $x->{sign};
1432   $u->{value} = $MBI->_div( $MBI->_copy($x->{_n}), $x->{_d});   # 22/7 => 3
1433   $u;
1434   }
1435
1436 sub as_float
1437   {
1438   # return N/D as Math::BigFloat
1439
1440   # set up parameters
1441   my ($self,$x,@r) = (ref($_[0]),@_);
1442   # objectify is costly, so avoid it
1443   ($self,$x,@r) = objectify(1,$class,@_) unless ref $_[0];
1444
1445   # NaN, inf etc
1446   return Math::BigFloat->new($x->{sign}) if $x->{sign} !~ /^[+-]$/;
1447
1448   my $u = Math::BigFloat->bzero();
1449   $u->{sign} = $x->{sign};
1450   # n
1451   $u->{_m} = $MBI->_copy($x->{_n});
1452   $u->{_e} = $MBI->_zero();
1453   $u->bdiv( $MBI->_str($x->{_d}), @r);
1454   # return $u
1455   $u;
1456   }
1457
1458 sub as_bin
1459   {
1460   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1461
1462   return $x unless $x->is_int();
1463
1464   my $s = $x->{sign}; $s = '' if $s eq '+';
1465   $s . $MBI->_as_bin($x->{_n});
1466   }
1467
1468 sub as_hex
1469   {
1470   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1471
1472   return $x unless $x->is_int();
1473
1474   my $s = $x->{sign}; $s = '' if $s eq '+';
1475   $s . $MBI->_as_hex($x->{_n});
1476   }
1477
1478 sub as_oct
1479   {
1480   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1481
1482   return $x unless $x->is_int();
1483
1484   my $s = $x->{sign}; $s = '' if $s eq '+';
1485   $s . $MBI->_as_oct($x->{_n});
1486   }
1487
1488 ##############################################################################
1489
1490 sub from_hex
1491   {
1492   my $class = shift;
1493
1494   $class->new(@_);
1495   }
1496
1497 sub from_bin
1498   {
1499   my $class = shift;
1500
1501   $class->new(@_);
1502   }
1503
1504 sub from_oct
1505   {
1506   my $class = shift;
1507
1508   my @parts;
1509   for my $c (@_)
1510     {
1511     push @parts, Math::BigInt->from_oct($c);
1512     }
1513   $class->new ( @parts );
1514   }
1515
1516 ##############################################################################
1517 # import
1518
1519 sub import
1520   {
1521   my $self = shift;
1522   my $l = scalar @_;
1523   my $lib = ''; my @a;
1524   my $try = 'try';
1525
1526   for ( my $i = 0; $i < $l ; $i++)
1527     {
1528     if ( $_[$i] eq ':constant' )
1529       {
1530       # this rest causes overlord er load to step in
1531       overload::constant float => sub { $self->new(shift); };
1532       }
1533 #    elsif ($_[$i] eq 'upgrade')
1534 #      {
1535 #     # this causes upgrading
1536 #      $upgrade = $_[$i+1];             # or undef to disable
1537 #      $i++;
1538 #      }
1539     elsif ($_[$i] eq 'downgrade')
1540       {
1541       # this causes downgrading
1542       $downgrade = $_[$i+1];            # or undef to disable
1543       $i++;
1544       }
1545     elsif ($_[$i] =~ /^(lib|try|only)\z/)
1546       {
1547       $lib = $_[$i+1] || '';            # default Calc
1548       $try = $1;                        # lib, try or only
1549       $i++;
1550       }
1551     elsif ($_[$i] eq 'with')
1552       {
1553       # this argument is no longer used
1554       #$MBI = $_[$i+1] || 'Math::BigInt::Calc'; # default Math::BigInt::Calc
1555       $i++;
1556       }
1557     else
1558       {
1559       push @a, $_[$i];
1560       }
1561     }
1562   require Math::BigInt;
1563
1564   # let use Math::BigInt lib => 'GMP'; use Math::BigRat; still have GMP
1565   if ($lib ne '')
1566     {
1567     my @c = split /\s*,\s*/, $lib;
1568     foreach (@c)
1569       {
1570       $_ =~ tr/a-zA-Z0-9://cd;                    # limit to sane characters
1571       }
1572     $lib = join(",", @c);
1573     }
1574   my @import = ('objectify');
1575   push @import, $try => $lib if $lib ne '';
1576
1577   # MBI already loaded, so feed it our lib arguments
1578   Math::BigInt->import( @import );
1579
1580   $MBI = Math::BigFloat->config()->{lib};
1581
1582   # register us with MBI to get notified of future lib changes
1583   Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
1584
1585   # any non :constant stuff is handled by our parent, Exporter (loaded
1586   # by Math::BigFloat, even if @_ is empty, to give it a chance
1587   $self->SUPER::import(@a);             # for subclasses
1588   $self->export_to_level(1,$self,@a);   # need this, too
1589   }
1590
1591 1;
1592
1593 __END__
1594
1595 =head1 NAME
1596
1597 Math::BigRat - Arbitrary big rational numbers
1598
1599 =head1 SYNOPSIS
1600
1601         use Math::BigRat;
1602
1603         my $x = Math::BigRat->new('3/7'); $x += '5/9';
1604
1605         print $x->bstr(),"\n";
1606         print $x ** 2,"\n";
1607
1608         my $y = Math::BigRat->new('inf');
1609         print "$y ", ($y->is_inf ? 'is' : 'is not') , " infinity\n";
1610
1611         my $z = Math::BigRat->new(144); $z->bsqrt();
1612
1613 =head1 DESCRIPTION
1614
1615 Math::BigRat complements Math::BigInt and Math::BigFloat by providing support
1616 for arbitrary big rational numbers.
1617
1618 =head2 MATH LIBRARY
1619
1620 You can change the underlying module that does the low-level
1621 math operations by using:
1622
1623         use Math::BigRat try => 'GMP';
1624
1625 Note: This needs Math::BigInt::GMP installed.
1626
1627 The following would first try to find Math::BigInt::Foo, then
1628 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
1629
1630         use Math::BigRat try => 'Foo,Math::BigInt::Bar';
1631
1632 If you want to get warned when the fallback occurs, replace "try" with
1633 "lib":
1634
1635         use Math::BigRat lib => 'Foo,Math::BigInt::Bar';
1636
1637 If you want the code to die instead, replace "try" with
1638 "only":
1639
1640         use Math::BigRat only => 'Foo,Math::BigInt::Bar';
1641
1642 =head1 METHODS
1643
1644 Any methods not listed here are derived from Math::BigFloat (or
1645 Math::BigInt), so make sure you check these two modules for further
1646 information.
1647
1648 =head2 new()
1649
1650         $x = Math::BigRat->new('1/3');
1651
1652 Create a new Math::BigRat object. Input can come in various forms:
1653
1654         $x = Math::BigRat->new(123);                            # scalars
1655         $x = Math::BigRat->new('inf');                          # infinity
1656         $x = Math::BigRat->new('123.3');                        # float
1657         $x = Math::BigRat->new('1/3');                          # simple string
1658         $x = Math::BigRat->new('1 / 3');                        # spaced
1659         $x = Math::BigRat->new('1 / 0.1');                      # w/ floats
1660         $x = Math::BigRat->new(Math::BigInt->new(3));           # BigInt
1661         $x = Math::BigRat->new(Math::BigFloat->new('3.1'));     # BigFloat
1662         $x = Math::BigRat->new(Math::BigInt::Lite->new('2'));   # BigLite
1663
1664         # You can also give D and N as different objects:
1665         $x = Math::BigRat->new(
1666                 Math::BigInt->new(-123),
1667                 Math::BigInt->new(7),
1668                 );                      # => -123/7
1669
1670 =head2 numerator()
1671
1672         $n = $x->numerator();
1673
1674 Returns a copy of the numerator (the part above the line) as signed BigInt.
1675
1676 =head2 denominator()
1677
1678         $d = $x->denominator();
1679
1680 Returns a copy of the denominator (the part under the line) as positive BigInt.
1681
1682 =head2 parts()
1683
1684         ($n,$d) = $x->parts();
1685
1686 Return a list consisting of (signed) numerator and (unsigned) denominator as
1687 BigInts.
1688
1689 =head2 numify()
1690
1691         my $y = $x->numify();
1692
1693 Returns the object as a scalar. This will lose some data if the object
1694 cannot be represented by a normal Perl scalar (integer or float), so
1695 use L<as_int()> or L<as_float()> instead.
1696
1697 This routine is automatically used whenever a scalar is required:
1698
1699         my $x = Math::BigRat->new('3/1');
1700         @array = (0,1,2,3);
1701         $y = $array[$x];                # set $y to 3
1702
1703 =head2 as_int()/as_number()
1704
1705         $x = Math::BigRat->new('13/7');
1706         print $x->as_int(),"\n";                # '1'
1707
1708 Returns a copy of the object as BigInt, truncated to an integer.
1709
1710 C<as_number()> is an alias for C<as_int()>.
1711
1712 =head2 as_float()
1713
1714         $x = Math::BigRat->new('13/7');
1715         print $x->as_float(),"\n";              # '1'
1716
1717         $x = Math::BigRat->new('2/3');
1718         print $x->as_float(5),"\n";             # '0.66667'
1719
1720 Returns a copy of the object as BigFloat, preserving the
1721 accuracy as wanted, or the default of 40 digits.
1722
1723 This method was added in v0.22 of Math::BigRat (April 2008).
1724
1725 =head2 as_hex()
1726
1727         $x = Math::BigRat->new('13');
1728         print $x->as_hex(),"\n";                # '0xd'
1729
1730 Returns the BigRat as hexadecimal string. Works only for integers.
1731
1732 =head2 as_bin()
1733
1734         $x = Math::BigRat->new('13');
1735         print $x->as_bin(),"\n";                # '0x1101'
1736
1737 Returns the BigRat as binary string. Works only for integers.
1738
1739 =head2 as_oct()
1740
1741         $x = Math::BigRat->new('13');
1742         print $x->as_oct(),"\n";                # '015'
1743
1744 Returns the BigRat as octal string. Works only for integers.
1745
1746 =head2 from_hex()/from_bin()/from_oct()
1747
1748         my $h = Math::BigRat->from_hex('0x10');
1749         my $b = Math::BigRat->from_bin('0b10000000');
1750         my $o = Math::BigRat->from_oct('020');
1751
1752 Create a BigRat from an hexadecimal, binary or octal number
1753 in string form.
1754
1755 =head2 length()
1756
1757         $len = $x->length();
1758
1759 Return the length of $x in digits for integer values.
1760
1761 =head2 digit()
1762
1763         print Math::BigRat->new('123/1')->digit(1);     # 1
1764         print Math::BigRat->new('123/1')->digit(-1);    # 3
1765
1766 Return the N'ths digit from X when X is an integer value.
1767
1768 =head2 bnorm()
1769
1770         $x->bnorm();
1771
1772 Reduce the number to the shortest form. This routine is called
1773 automatically whenever it is needed.
1774
1775 =head2 bfac()
1776
1777         $x->bfac();
1778
1779 Calculates the factorial of $x. For instance:
1780
1781         print Math::BigRat->new('3/1')->bfac(),"\n";    # 1*2*3
1782         print Math::BigRat->new('5/1')->bfac(),"\n";    # 1*2*3*4*5
1783
1784 Works currently only for integers.
1785
1786 =head2 bround()/round()/bfround()
1787
1788 Are not yet implemented.
1789
1790 =head2 bmod()
1791
1792         use Math::BigRat;
1793         my $x = Math::BigRat->new('7/4');
1794         my $y = Math::BigRat->new('4/3');
1795         print $x->bmod($y);
1796
1797 Set $x to the remainder of the division of $x by $y.
1798
1799 =head2 bneg()
1800
1801         $x->bneg();
1802
1803 Used to negate the object in-place.
1804
1805 =head2 is_one()
1806
1807         print "$x is 1\n" if $x->is_one();
1808
1809 Return true if $x is exactly one, otherwise false.
1810
1811 =head2 is_zero()
1812
1813         print "$x is 0\n" if $x->is_zero();
1814
1815 Return true if $x is exactly zero, otherwise false.
1816
1817 =head2 is_pos()/is_positive()
1818
1819         print "$x is >= 0\n" if $x->is_positive();
1820
1821 Return true if $x is positive (greater than or equal to zero), otherwise
1822 false. Please note that '+inf' is also positive, while 'NaN' and '-inf' aren't.
1823
1824 C<is_positive()> is an alias for C<is_pos()>.
1825
1826 =head2 is_neg()/is_negative()
1827
1828         print "$x is < 0\n" if $x->is_negative();
1829
1830 Return true if $x is negative (smaller than zero), otherwise false. Please
1831 note that '-inf' is also negative, while 'NaN' and '+inf' aren't.
1832
1833 C<is_negative()> is an alias for C<is_neg()>.
1834
1835 =head2 is_int()
1836
1837         print "$x is an integer\n" if $x->is_int();
1838
1839 Return true if $x has a denominator of 1 (e.g. no fraction parts), otherwise
1840 false. Please note that '-inf', 'inf' and 'NaN' aren't integer.
1841
1842 =head2 is_odd()
1843
1844         print "$x is odd\n" if $x->is_odd();
1845
1846 Return true if $x is odd, otherwise false.
1847
1848 =head2 is_even()
1849
1850         print "$x is even\n" if $x->is_even();
1851
1852 Return true if $x is even, otherwise false.
1853
1854 =head2 bceil()
1855
1856         $x->bceil();
1857
1858 Set $x to the next bigger integer value (e.g. truncate the number to integer
1859 and then increment it by one).
1860
1861 =head2 bfloor()
1862
1863         $x->bfloor();
1864
1865 Truncate $x to an integer value.
1866
1867 =head2 bsqrt()
1868
1869         $x->bsqrt();
1870
1871 Calculate the square root of $x.
1872
1873 =head2 broot()
1874
1875         $x->broot($n);
1876
1877 Calculate the N'th root of $x.
1878
1879 =head2 badd()/bmul()/bsub()/bdiv()/bdec()/binc()
1880
1881 Please see the documentation in L<Math::BigInt>.
1882
1883 =head2 copy()
1884
1885         my $z = $x->copy();
1886
1887 Makes a deep copy of the object.
1888
1889 Please see the documentation in L<Math::BigInt> for further details.
1890
1891 =head2 bstr()/bsstr()
1892
1893         my $x = Math::BigInt->new('8/4');
1894         print $x->bstr(),"\n";                  # prints 1/2
1895         print $x->bsstr(),"\n";                 # prints 1/2
1896
1897 Return a string representing this object.
1898
1899 =head2 bacmp()/bcmp()
1900
1901 Used to compare numbers.
1902
1903 Please see the documentation in L<Math::BigInt> for further details.
1904
1905 =head2 blsft()/brsft()
1906
1907 Used to shift numbers left/right.
1908
1909 Please see the documentation in L<Math::BigInt> for further details.
1910
1911 =head2 bpow()
1912
1913         $x->bpow($y);
1914
1915 Compute $x ** $y.
1916
1917 Please see the documentation in L<Math::BigInt> for further details.
1918
1919 =head2 bexp()
1920
1921         $x->bexp($accuracy);            # calculate e ** X
1922
1923 Calculates two integers A and B so that A/B is equal to C<e ** $x>, where C<e> is
1924 Euler's number.
1925
1926 This method was added in v0.20 of Math::BigRat (May 2007).
1927
1928 See also L<blog()>.
1929
1930 =head2 bnok()
1931
1932         $x->bnok($y);              # x over y (binomial coefficient n over k)
1933
1934 Calculates the binomial coefficient n over k, also called the "choose"
1935 function. The result is equivalent to:
1936
1937         ( n )      n!
1938         | - |  = -------
1939         ( k )    k!(n-k)!
1940
1941 This method was added in v0.20 of Math::BigRat (May 2007).
1942
1943 =head2 config()
1944
1945         use Data::Dumper;
1946
1947         print Dumper ( Math::BigRat->config() );
1948         print Math::BigRat->config()->{lib},"\n";
1949
1950 Returns a hash containing the configuration, e.g. the version number, lib
1951 loaded etc. The following hash keys are currently filled in with the
1952 appropriate information.
1953
1954         key             RO/RW   Description
1955                                 Example
1956         ============================================================
1957         lib             RO      Name of the Math library
1958                                 Math::BigInt::Calc
1959         lib_version     RO      Version of 'lib'
1960                                 0.30
1961         class           RO      The class of config you just called
1962                                 Math::BigRat
1963         version         RO      version number of the class you used
1964                                 0.10
1965         upgrade         RW      To which class numbers are upgraded
1966                                 undef
1967         downgrade       RW      To which class numbers are downgraded
1968                                 undef
1969         precision       RW      Global precision
1970                                 undef
1971         accuracy        RW      Global accuracy
1972                                 undef
1973         round_mode      RW      Global round mode
1974                                 even
1975         div_scale       RW      Fallback accuracy for div
1976                                 40
1977         trap_nan        RW      Trap creation of NaN (undef = no)
1978                                 undef
1979         trap_inf        RW      Trap creation of +inf/-inf (undef = no)
1980                                 undef
1981
1982 By passing a reference to a hash you may set the configuration values. This
1983 works only for values that a marked with a C<RW> above, anything else is
1984 read-only.
1985
1986 =head2 objectify()
1987
1988 This is an internal routine that turns scalars into objects.
1989
1990 =head1 BUGS
1991
1992 Some things are not yet implemented, or only implemented half-way:
1993
1994 =over 2
1995
1996 =item inf handling (partial)
1997
1998 =item NaN handling (partial)
1999
2000 =item rounding (not implemented except for bceil/bfloor)
2001
2002 =item $x ** $y where $y is not an integer
2003
2004 =item bmod(), blog(), bmodinv() and bmodpow() (partial)
2005
2006 =back
2007
2008 =head1 LICENSE
2009
2010 This program is free software; you may redistribute it and/or modify it under
2011 the same terms as Perl itself.
2012
2013 =head1 SEE ALSO
2014
2015 L<Math::BigFloat> and L<Math::Big> as well as L<Math::BigInt::BitVect>,
2016 L<Math::BigInt::Pari> and  L<Math::BigInt::GMP>.
2017
2018 See L<http://search.cpan.org/search?dist=bignum> for a way to use
2019 Math::BigRat.
2020
2021 The package at L<http://search.cpan.org/search?dist=Math%3A%3ABigRat>
2022 may contain more documentation and examples as well as testcases.
2023
2024 =head1 AUTHORS
2025
2026 (C) by Tels L<http://bloodgate.com/> 2001 - 2009.
2027
2028 Currently maintained by Jonathan "Duke" Leto <jonathan@leto.net> L<http://leto.net>
2029
2030 =cut