lib/Math/BigInt/CalcEmu.pm

   1 package Math::BigInt;
   2
   3 use 5.005;
   4 use strict;
   5 # use warnings; # dont use warnings for older Perls
   6
   7 use vars qw/$VERSION/;
   8
   9 $VERSION = '0.01';
  10
  11 # See SYNOPSIS below.
  12
  13 my $CALC_EMU;
  14
  15 BEGIN
  16   {
  17   $CALC_EMU = Math::BigInt->config()->{'lib'};
  18   }
  19
  20 sub __emu_blog
  21   {
  22   my ($self,$x,$base,@r) = @_;
  23
  24   return $x->bnan() if $x->is_zero() || $base->is_zero() || $base->is_one();
  25
  26   my $acmp = $x->bacmp($base);
  27   return $x->bone('+',@r) if $acmp == 0;
  28   return $x->bzero(@r) if $acmp < 0 || $x->is_one();
  29
  30   # blog($x,$base) ** $base + $y = $x
  31
  32   # this trial multiplication is very fast, even for large counts (like for
  33   # 2 ** 1024, since this still requires only 1024 very fast steps
  34   # (multiplication of a large number by a very small number is very fast))
  35   # See Calc for an even faster algorightmn
  36   my $x_org = $x->copy();               # preserve orgx
  37   $x->bzero();                          # keep ref to $x
  38   my $trial = $base->copy();
  39   while ($trial->bacmp($x_org) <= 0)
  40     {
  41     $trial->bmul($base); $x->binc();
  42     }
  43   $x->round(@r);
  44   }
  45
  46 sub __emu_bmodinv
  47   {
  48   my ($self,$x,$y,@r) = @_;
  49
  50   my ($u, $u1) = ($self->bzero(), $self->bone());
  51   my ($a, $b) = ($y->copy(), $x->copy());
  52
  53   # first step need always be done since $num (and thus $b) is never 0
  54   # Note that the loop is aligned so that the check occurs between #2 and #1
  55   # thus saving us one step #2 at the loop end. Typical loop count is 1. Even
  56   # a case with 28 loops still gains about 3% with this layout.
  57   my $q;
  58   ($a, $q, $b) = ($b, $a->bdiv($b));                    # step #1
  59   # Euclid's Algorithm (calculate GCD of ($a,$b) in $a and also calculate
  60   # two values in $u and $u1, we use only $u1 afterwards)
  61   my $sign = 1;                                         # flip-flop
  62   while (!$b->is_zero())                                # found GCD if $b == 0
  63     {
  64     # the original algorithm had:
  65     # ($u, $u1) = ($u1, $u->bsub($u1->copy()->bmul($q))); # step #2
  66     # The following creates exact the same sequence of numbers in $u1,
  67     # except for the sign ($u1 is now always positive). Since formerly
  68     # the sign of $u1 was alternating between '-' and '+', the $sign
  69     # flip-flop will take care of that, so that at the end of the loop
  70     # we have the real sign of $u1. Keeping numbers positive gains us
  71     # speed since badd() is faster than bsub() and makes it possible
  72     # to have the algorithmn in Calc for even more speed.
  73
  74     ($u, $u1) = ($u1, $u->badd($u1->copy()->bmul($q))); # step #2
  75     $sign = - $sign;                                    # flip sign
  76
  77     ($a, $q, $b) = ($b, $a->bdiv($b));                  # step #1 again
  78     }
  79
  80   # If the gcd is not 1, then return NaN! It would be pointless to have
  81   # called bgcd to check this first, because we would then be performing
  82   # the same Euclidean Algorithm *twice* in case the gcd is 1.
  83   return $x->bnan() unless $a->is_one();
  84
  85   $u1->bneg() if $sign != 1;                            # need to flip?
  86
  87   $u1->bmod($y);                                        # calc result
  88   $x->{value} = $u1->{value};                           # and copy over to $x
  89   $x->{sign} = $u1->{sign};                             # to modify in place
  90   $x->round(@r);
  91   }
  92
  93 sub __emu_bmodpow
  94   {
  95   my ($self,$num,$exp,$mod,@r) = @_;
  96
  97   # in the trivial case,
  98   return $num->bzero(@r) if $mod->is_one();
  99   return $num->bone('+',@r) if $num->is_zero() or $num->is_one();
 100
 101   # $num->bmod($mod);           # if $x is large, make it smaller first
 102   my $acc = $num->copy();       # but this is not really faster...
 103
 104   $num->bone(); # keep ref to $num
 105
 106   my $expbin = $exp->as_bin(); $expbin =~ s/^[-]?0b//; # ignore sign and prefix
 107   my $len = CORE::length($expbin);
 108   while (--$len >= 0)
 109     {
 110     $num->bmul($acc)->bmod($mod) if substr($expbin,$len,1) eq '1';
 111     $acc->bmul($acc)->bmod($mod);
 112     }
 113
 114   $num->round(@r);
 115   }
 116
 117 sub __emu_bfac
 118   {
 119   my ($self,$x,@r) = @_;
 120
 121   return $x->bone('+',@r) if $x->is_zero() || $x->is_one();     # 0 or 1 => 1
 122
 123   my $n = $x->copy();
 124   $x->bone();
 125   # seems we need not to temp. clear A/P of $x since the result is the same
 126   my $f = $self->new(2);
 127   while ($f->bacmp($n) < 0)
 128     {
 129     $x->bmul($f); $f->binc();
 130     }
 131   $x->bmul($f,@r);                      # last step and also round result
 132   }
 133
 134 sub __emu_bpow
 135   {
 136   my ($self,$x,$y,@r) = @_;
 137
 138   return $x->bone('+',@r) if $y->is_zero();
 139   return $x->round(@r) if $x->is_one() || $y->is_one();
 140   return $x->round(@r) if $x->is_zero();  # 0**y => 0 (if not y <= 0)
 141
 142   my $pow2 = $self->bone();
 143   my $y_bin = $y->as_bin(); $y_bin =~ s/^0b//;
 144   my $len = CORE::length($y_bin);
 145   while (--$len > 0)
 146     {
 147     $pow2->bmul($x) if substr($y_bin,$len,1) eq '1';    # is odd?
 148     $x->bmul($x);
 149     }
 150   $x->bmul($pow2);
 151   $x->round(@r) if !exists $x->{_f} || $x->{_f} & MB_NEVER_ROUND == 0;
 152   $x;
 153   }
 154
 155 sub __emu_band
 156   {
 157   my ($self,$x,$y,$sx,$sy,@r) = @_;
 158
 159   return $x->bzero(@r) if $y->is_zero() || $x->is_zero();
 160
 161   my $sign = 0;                                 # sign of result
 162   $sign = 1 if $sx == -1 && $sy == -1;
 163
 164   my ($bx,$by);
 165
 166   if ($sx == -1)                                # if x is negative
 167     {
 168     # two's complement: inc and flip all "bits" in $bx
 169     $bx = $x->binc()->as_hex();                 # -1 => 0, -2 => 1, -3 => 2 etc
 170     $bx =~ s/-?0x//;
 171     $bx =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/;
 172     }
 173   else
 174     {
 175     $bx = $x->as_hex();                         # get binary representation
 176     $bx =~ s/-?0x//;
 177     $bx =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/;
 178     }
 179   if ($sy == -1)                                # if y is negative
 180     {
 181     # two's complement: inc and flip all "bits" in $by
 182     $by = $y->copy()->binc()->as_hex();         # -1 => 0, -2 => 1, -3 => 2 etc
 183     $by =~ s/-?0x//;
 184     $by =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/;
 185     }
 186   else
 187     {
 188     $by = $y->as_hex();                         # get binary representation
 189     $by =~ s/-?0x//;
 190     $by =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/;
 191     }
 192   # now we have bit-strings from X and Y, reverse them for padding
 193   $bx = reverse $bx;
 194   $by = reverse $by;
 195
 196   # cut the longer string to the length of the shorter one (the result would
 197   # be 0 due to AND anyway)
 198   my $diff = CORE::length($bx) - CORE::length($by);
 199   if ($diff > 0)
 200     {
 201     $bx = substr($bx,0,CORE::length($by));
 202     }
 203   elsif ($diff < 0)
 204     {
 205     $by = substr($by,0,CORE::length($bx));
 206     }
 207
 208   # and the strings together
 209   my $r = $bx & $by;
 210
 211   # and reverse the result again
 212   $bx = reverse $r;
 213
 214   # one of $x or $y was negative, so need to flip bits in the result
 215   # in both cases (one or two of them negative, or both positive) we need
 216   # to get the characters back.
 217   if ($sign == 1)
 218     {
 219     $bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/0123456789abcdef/;
 220     }
 221   else
 222     {
 223     $bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/fedcba9876543210/;
 224     }
 225
 226   $bx = '0x' . $bx;
 227   if ($CALC_EMU->can('_from_hex'))
 228     {
 229     $x->{value} = $CALC_EMU->_from_hex( \$bx );
 230     }
 231   else
 232     {
 233     $r = $self->new($bx);
 234     $x->{value} = $r->{value};
 235     }
 236
 237   # calculate sign of result
 238   $x->{sign} = '+';
 239   #$x->{sign} = '-' if $sx == $sy && $sx == -1 && !$x->is_zero();
 240   $x->{sign} = '-' if $sign == 1 && !$x->is_zero();
 241
 242   $x->bdec() if $sign == 1;
 243
 244   $x->round(@r);
 245   }
 246
 247 sub __emu_bior
 248   {
 249   my ($self,$x,$y,$sx,$sy,@r) = @_;
 250
 251   return $x->round(@r) if $y->is_zero();
 252
 253   my $sign = 0;                                 # sign of result
 254   $sign = 1 if ($sx == -1) || ($sy == -1);
 255
 256   my ($bx,$by);
 257
 258   if ($sx == -1)                                # if x is negative
 259     {
 260     # two's complement: inc and flip all "bits" in $bx
 261     $bx = $x->binc()->as_hex();                 # -1 => 0, -2 => 1, -3 => 2 etc
 262     $bx =~ s/-?0x//;
 263     $bx =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/;
 264     }
 265   else
 266     {
 267     $bx = $x->as_hex();                         # get binary representation
 268     $bx =~ s/-?0x//;
 269     $bx =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/;
 270     }
 271   if ($sy == -1)                                # if y is negative
 272     {
 273     # two's complement: inc and flip all "bits" in $by
 274     $by = $y->copy()->binc()->as_hex();         # -1 => 0, -2 => 1, -3 => 2 etc
 275     $by =~ s/-?0x//;
 276     $by =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/;
 277     }
 278   else
 279     {
 280     $by = $y->as_hex();                         # get binary representation
 281     $by =~ s/-?0x//;
 282     $by =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/;
 283     }
 284   # now we have bit-strings from X and Y, reverse them for padding
 285   $bx = reverse $bx;
 286   $by = reverse $by;
 287
 288   # padd the shorter string
 289   my $xx = "\x00"; $xx = "\x0f" if $sx == -1;
 290   my $yy = "\x00"; $yy = "\x0f" if $sy == -1;
 291   my $diff = CORE::length($bx) - CORE::length($by);
 292   if ($diff > 0)
 293     {
 294     $by .= $yy x $diff;
 295     }
 296   elsif ($diff < 0)
 297     {
 298     $bx .= $xx x abs($diff);
 299     }
 300
 301   # or the strings together
 302   my $r = $bx | $by;
 303
 304   # and reverse the result again
 305   $bx = reverse $r;
 306
 307   # one of $x or $y was negative, so need to flip bits in the result
 308   # in both cases (one or two of them negative, or both positive) we need
 309   # to get the characters back.
 310   if ($sign == 1)
 311     {
 312     $bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/0123456789abcdef/;
 313     }
 314   else
 315     {
 316     $bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/fedcba9876543210/;
 317     }
 318
 319   $bx = '0x' . $bx;
 320   if ($CALC_EMU->can('_from_hex'))
 321     {
 322     $x->{value} = $CALC_EMU->_from_hex( \$bx );
 323     }
 324   else
 325     {
 326     $r = $self->new($bx);
 327     $x->{value} = $r->{value};
 328     }
 329
 330   # if one of X or Y was negative, we need to decrement result
 331   $x->bdec() if $sign == 1;
 332
 333   $x->round(@r);
 334   }
 335
 336 sub __emu_bxor
 337   {
 338   my ($self,$x,$y,$sx,$sy,@r) = @_;
 339
 340   return $x->round(@r) if $y->is_zero();
 341
 342   my $sign = 0;                                 # sign of result
 343   $sign = 1 if $x->{sign} ne $y->{sign};
 344
 345   my ($bx,$by);
 346
 347   if ($sx == -1)                                # if x is negative
 348     {
 349     # two's complement: inc and flip all "bits" in $bx
 350     $bx = $x->binc()->as_hex();                 # -1 => 0, -2 => 1, -3 => 2 etc
 351     $bx =~ s/-?0x//;
 352     $bx =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/;
 353     }
 354   else
 355     {
 356     $bx = $x->as_hex();                         # get binary representation
 357     $bx =~ s/-?0x//;
 358     $bx =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/;
 359     }
 360   if ($sy == -1)                                # if y is negative
 361     {
 362     # two's complement: inc and flip all "bits" in $by
 363     $by = $y->copy()->binc()->as_hex();         # -1 => 0, -2 => 1, -3 => 2 etc
 364     $by =~ s/-?0x//;
 365     $by =~ tr/0123456789abcdef/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/;
 366     }
 367   else
 368     {
 369     $by = $y->as_hex();                         # get binary representation
 370     $by =~ s/-?0x//;
 371     $by =~ tr/fedcba9876543210/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/;
 372     }
 373   # now we have bit-strings from X and Y, reverse them for padding
 374   $bx = reverse $bx;
 375   $by = reverse $by;
 376
 377   # padd the shorter string
 378   my $xx = "\x00"; $xx = "\x0f" if $sx == -1;
 379   my $yy = "\x00"; $yy = "\x0f" if $sy == -1;
 380   my $diff = CORE::length($bx) - CORE::length($by);
 381   if ($diff > 0)
 382     {
 383     $by .= $yy x $diff;
 384     }
 385   elsif ($diff < 0)
 386     {
 387     $bx .= $xx x abs($diff);
 388     }
 389
 390   # xor the strings together
 391   my $r = $bx ^ $by;
 392
 393   # and reverse the result again
 394   $bx = reverse $r;
 395
 396   # one of $x or $y was negative, so need to flip bits in the result
 397   # in both cases (one or two of them negative, or both positive) we need
 398   # to get the characters back.
 399   if ($sign == 1)
 400     {
 401     $bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/0123456789abcdef/;
 402     }
 403   else
 404     {
 405     $bx =~ tr/\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00/fedcba9876543210/;
 406     }
 407
 408   $bx = '0x' . $bx;
 409   if ($CALC_EMU->can('_from_hex'))
 410     {
 411     $x->{value} = $CALC_EMU->_from_hex( \$bx );
 412     }
 413   else
 414     {
 415     $r = $self->new($bx);
 416     $x->{value} = $r->{value};
 417     }
 418
 419   # calculate sign of result
 420   $x->{sign} = '+';
 421   $x->{sign} = '-' if $sx != $sy && !$x->is_zero();
 422
 423   $x->bdec() if $sign == 1;
 424
 425   $x->round(@r);
 426   }
 427
 428 sub __emu_bsqrt
 429   {
 430   my ($self,$x,@r) = @_;
 431
 432   # this is slow:
 433   return $x->round(@r) if $x->is_zero();        # 0,1 => 0,1
 434
 435   return $x->bone('+',@r) if $x < 4;            # 1,2,3 => 1
 436   my $y = $x->copy();
 437   my $l = int($x->length()/2);
 438
 439   $x->bone();                                   # keep ref($x), but modify it
 440   $x->blsft($l,10) if $l != 0;                  # first guess: 1.('0' x (l/2))
 441
 442   my $last = $self->bzero();
 443   my $two = $self->new(2);
 444   my $lastlast = $self->bzero();
 445   #my $lastlast = $x+$two;
 446   while ($last != $x && $lastlast != $x)
 447     {
 448     $lastlast = $last; $last = $x->copy();
 449     $x->badd($y / $x);
 450     $x->bdiv($two);
 451     }
 452   $x->bdec() if $x * $x > $y;                   # overshot?
 453   $x->round(@r);
 454   }
 455
 456 sub __emu_broot
 457   {
 458   my ($self,$x,$y,@r) = @_;
 459
 460   return $x->bsqrt() if $y->bacmp(2) == 0;      # 2 => square root
 461
 462   # since we take at least a cubic root, and only 8 ** 1/3 >= 2 (==2):
 463   return $x->bone('+',@r) if $x < 8;            # $x=2..7 => 1
 464
 465   my $num = $x->numify();
 466
 467   if ($num <= 1000000)
 468     {
 469     $x = $self->new( int ( sprintf ("%.8f", $num ** (1 / $y->numify() ))));
 470     return $x->round(@r);
 471     }
 472
 473   # if $n is a power of two, we can repeatedly take sqrt($X) and find the
 474   # proper result, because sqrt(sqrt($x)) == root($x,4)
 475   # See Calc.pm for more details
 476   my $b = $y->as_bin();
 477   if ($b =~ /0b1(0+)/)
 478     {
 479     my $count = CORE::length($1);               # 0b100 => len('00') => 2
 480     my $cnt = $count;                           # counter for loop
 481     my $shift = $self->new(6);
 482     $x->blsft($shift);                          # add some zeros (even amount)
 483     while ($cnt-- > 0)
 484       {
 485       # 'inflate' $X by adding more zeros
 486       $x->blsft($shift);
 487       # calculate sqrt($x), $x is now a bit too big, again. In the next
 488       # round we make even bigger, again.
 489       $x->bsqrt($x);
 490       }
 491     # $x is still to big, so truncate result
 492     $x->brsft($shift);
 493     }
 494   else
 495     {
 496     # Should compute a guess of the result (by rule of thumb), then improve it
 497     # via Newton's method or something similiar.
 498     # XXX TODO
 499     warn ('broot() not fully implemented in BigInt.');
 500     }
 501   $x->round(@r);
 502   }
 503
 504 sub __emu_as_hex
 505   {
 506   my ($self,$x,$s) = @_;
 507
 508   return '0x0' if $x->is_zero();
 509
 510   my $x1 = $x->copy()->babs(); my ($xr,$x10000,$h,$es);
 511   if ($] >= 5.006)
 512     {
 513     $x10000 = $self->new (0x10000); $h = 'h4';
 514     }
 515   else
 516     {
 517     $x10000 = $self->new (0x1000); $h = 'h3';
 518     }
 519   while (!$x1->is_zero())
 520     {
 521     ($x1, $xr) = bdiv($x1,$x10000);
 522     $es .= unpack($h,pack('v',$xr->numify()));
 523     }
 524   $es = reverse $es;
 525   $es =~ s/^[0]+//;             # strip leading zeros
 526   $s . '0x' . $es;
 527   }
 528
 529 sub __emu_as_bin
 530   {
 531   my ($self,$x,$s) = @_;
 532
 533   return '0b0' if $x->is_zero();
 534
 535   my $x1 = $x->copy()->babs(); my ($xr,$x10000,$b,$es);
 536   if ($] >= 5.006)
 537     {
 538     $x10000 = $self->new (0x10000); $b = 'b16';
 539     }
 540   else
 541     {
 542     $x10000 = $self->new (0x1000); $b = 'b12';
 543     }
 544   while (!$x1->is_zero())
 545     {
 546     ($x1, $xr) = bdiv($x1,$x10000);
 547     $es .= unpack($b,pack('v',$xr->numify()));
 548     }
 549   $es = reverse $es;
 550   $es =~ s/^[0]+//;   # strip leading zeros
 551   $s . '0b' . $es;
 552   }
 553
 554 ##############################################################################
 555 ##############################################################################
 556
 557 1;
 558 __END__
 559
 560 =head1 NAME
 561
 562 Math::BigInt::CalcEmu - Emulate low-level math with BigInt code
 563
 564 =head1 SYNOPSIS
 565
 566 Contains routines that emulate low-level math functions in BigInt, e.g.
 567 optional routines the low-level math package does not provide on it's own.
 568
 569 Will be loaded on demand and automatically by BigInt.
 570
 571 Stuff here is really low-priority to optimize,
 572 since it is far better to implement the operation in the low-level math
 573 libary directly, possible even using a call to the native lib.
 574
 575 =head1 DESCRIPTION
 576
 577 =head1 METHODS
 578
 579 =head1 LICENSE
 580
 581 This program is free software; you may redistribute it and/or modify it under
 582 the same terms as Perl itself.
 583
 584 =head1 AUTHORS
 585
 586 (c) Tels http://bloodgate.com 2003 - based on BigInt code by
 587 Tels from 2001-2003.
 588
 589 =head1 SEE ALSO
 590
 591 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
 592 L<Math::BigInt::GMP> and L<Math::BigInt::Pari>.
 593
 594 =cut