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