1package Math::BigInt::Lib;
2
3use 5.006001;
4use strict;
5use warnings;
6
7our $VERSION = '2.003002';
8$VERSION =~ tr/_//d;
9
10use Carp;
11
12use overload
13
14  # overload key: with_assign
15
16  '+'    => sub {
17                my $class = ref $_[0];
18                my $x = $class -> _copy($_[0]);
19                my $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
20                return $class -> _add($x, $y);
21            },
22
23  '-'    => sub {
24                my $class = ref $_[0];
25                my ($x, $y);
26                if ($_[2]) {            # if swapped
27                    $y = $_[0];
28                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
29                } else {
30                    $x = $class -> _copy($_[0]);
31                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
32                }
33                return $class -> _sub($x, $y);
34            },
35
36  '*'    => sub {
37                my $class = ref $_[0];
38                my $x = $class -> _copy($_[0]);
39                my $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
40                return $class -> _mul($x, $y);
41            },
42
43  '/'    => sub {
44                my $class = ref $_[0];
45                my ($x, $y);
46                if ($_[2]) {            # if swapped
47                    $y = $_[0];
48                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
49                } else {
50                    $x = $class -> _copy($_[0]);
51                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
52                }
53                return $class -> _div($x, $y);
54            },
55
56  '%'    => sub {
57                my $class = ref $_[0];
58                my ($x, $y);
59                if ($_[2]) {            # if swapped
60                    $y = $_[0];
61                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
62                } else {
63                    $x = $class -> _copy($_[0]);
64                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
65                }
66                return $class -> _mod($x, $y);
67            },
68
69  '**'   => sub {
70                my $class = ref $_[0];
71                my ($x, $y);
72                if ($_[2]) {            # if swapped
73                    $y = $_[0];
74                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
75                } else {
76                    $x = $class -> _copy($_[0]);
77                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
78                }
79                return $class -> _pow($x, $y);
80            },
81
82  '<<'   => sub {
83                my $class = ref $_[0];
84                my ($x, $y);
85                if ($_[2]) {            # if swapped
86                    $y = $class -> _num($_[0]);
87                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
88                } else {
89                    $x = $_[0];
90                    $y = ref($_[1]) ? $class -> _num($_[1]) : $_[1];
91                }
92                return $class -> _lsft($x, $y);
93            },
94
95  '>>'   => sub {
96                my $class = ref $_[0];
97                my ($x, $y);
98                if ($_[2]) {            # if swapped
99                    $y = $_[0];
100                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
101                } else {
102                    $x = $class -> _copy($_[0]);
103                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
104                }
105                return $class -> _rsft($x, $y);
106            },
107
108  # overload key: num_comparison
109
110  '<'    => sub {
111                my $class = ref $_[0];
112                my ($x, $y);
113                if ($_[2]) {            # if swapped
114                    $y = $_[0];
115                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
116                } else {
117                    $x = $class -> _copy($_[0]);
118                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
119                }
120                return $class -> _acmp($x, $y) < 0;
121            },
122
123  '<='   => sub {
124                my $class = ref $_[0];
125                my ($x, $y);
126                if ($_[2]) {            # if swapped
127                    $y = $_[0];
128                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
129                } else {
130                    $x = $class -> _copy($_[0]);
131                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
132                }
133                return $class -> _acmp($x, $y) <= 0;
134            },
135
136  '>'    => sub {
137                my $class = ref $_[0];
138                my ($x, $y);
139                if ($_[2]) {            # if swapped
140                    $y = $_[0];
141                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
142                } else {
143                    $x = $class -> _copy($_[0]);
144                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
145                }
146                return $class -> _acmp($x, $y) > 0;
147            },
148
149  '>='   => sub {
150                my $class = ref $_[0];
151                my ($x, $y);
152                if ($_[2]) {            # if swapped
153                    $y = $_[0];
154                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
155                } else {
156                    $x = $class -> _copy($_[0]);
157                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
158                }
159                return $class -> _acmp($x, $y) >= 0;
160          },
161
162  '=='   => sub {
163                my $class = ref $_[0];
164                my $x = $class -> _copy($_[0]);
165                my $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
166                return $class -> _acmp($x, $y) == 0;
167            },
168
169  '!='   => sub {
170                my $class = ref $_[0];
171                my $x = $class -> _copy($_[0]);
172                my $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
173                return $class -> _acmp($x, $y) != 0;
174            },
175
176  # overload key: 3way_comparison
177
178  '<=>'  => sub {
179                my $class = ref $_[0];
180                my ($x, $y);
181                if ($_[2]) {            # if swapped
182                    $y = $_[0];
183                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
184                } else {
185                    $x = $class -> _copy($_[0]);
186                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
187                }
188                return $class -> _acmp($x, $y);
189            },
190
191  # overload key: binary
192
193  '&'    => sub {
194                my $class = ref $_[0];
195                my ($x, $y);
196                if ($_[2]) {            # if swapped
197                    $y = $_[0];
198                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
199                } else {
200                    $x = $class -> _copy($_[0]);
201                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
202                }
203                return $class -> _and($x, $y);
204            },
205
206  '|'    => sub {
207                my $class = ref $_[0];
208                my ($x, $y);
209                if ($_[2]) {            # if swapped
210                    $y = $_[0];
211                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
212                } else {
213                    $x = $class -> _copy($_[0]);
214                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
215                }
216                return $class -> _or($x, $y);
217            },
218
219  '^'    => sub {
220                my $class = ref $_[0];
221                my ($x, $y);
222                if ($_[2]) {            # if swapped
223                    $y = $_[0];
224                    $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
225                } else {
226                    $x = $class -> _copy($_[0]);
227                    $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
228                }
229                return $class -> _xor($x, $y);
230            },
231
232  # overload key: func
233
234  'abs'  => sub { $_[0] },
235
236  'sqrt' => sub {
237                my $class = ref $_[0];
238                return $class -> _sqrt($class -> _copy($_[0]));
239            },
240
241  'int'  => sub { $_[0] },
242
243  # overload key: conversion
244
245  'bool' => sub { ref($_[0]) -> _is_zero($_[0]) ? '' : 1; },
246
247  '""'   => sub { ref($_[0]) -> _str($_[0]); },
248
249  '0+'   => sub { ref($_[0]) -> _num($_[0]); },
250
251  '='    => sub { ref($_[0]) -> _copy($_[0]); },
252
253  ;
254
255sub _new {
256    croak "@{[(caller 0)[3]]} method not implemented";
257}
258
259sub _zero {
260    my $class = shift;
261    return $class -> _new("0");
262}
263
264sub _one {
265    my $class = shift;
266    return $class -> _new("1");
267}
268
269sub _two {
270    my $class = shift;
271    return $class -> _new("2");
272
273}
274sub _ten {
275    my $class = shift;
276    return $class -> _new("10");
277}
278
279sub _1ex {
280    my ($class, $exp) = @_;
281    $exp = $class -> _num($exp) if ref($exp);
282    return $class -> _new("1" . ("0" x $exp));
283}
284
285sub _copy {
286    my ($class, $x) = @_;
287    return $class -> _new($class -> _str($x));
288}
289
290# catch and throw away
291sub import { }
292
293##############################################################################
294# convert back to string and number
295
296sub _str {
297    # Convert number from internal base 1eN format to string format. Internal
298    # format is always normalized, i.e., no leading zeros.
299    croak "@{[(caller 0)[3]]} method not implemented";
300}
301
302sub _num {
303    my ($class, $x) = @_;
304    0 + $class -> _str($x);
305}
306
307##############################################################################
308# actual math code
309
310sub _add {
311    croak "@{[(caller 0)[3]]} method not implemented";
312}
313
314sub _sub {
315    croak "@{[(caller 0)[3]]} method not implemented";
316}
317
318sub _mul {
319    my ($class, $x, $y) = @_;
320    my $sum = $class -> _zero();
321    my $i   = $class -> _zero();
322    while ($class -> _acmp($i, $y) < 0) {
323        $sum = $class -> _add($sum, $x);
324        $i   = $class -> _inc($i);
325    }
326    return $sum;
327}
328
329sub _div {
330    my ($class, $x, $y) = @_;
331
332    croak "@{[(caller 0)[3]]} requires non-zero divisor"
333      if $class -> _is_zero($y);
334
335    my $r = $class -> _copy($x);
336    my $q = $class -> _zero();
337    while ($class -> _acmp($r, $y) >= 0) {
338        $q = $class -> _inc($q);
339        $r = $class -> _sub($r, $y);
340    }
341
342    return $q, $r if wantarray;
343    return $q;
344}
345
346sub _inc {
347    my ($class, $x) = @_;
348    $class -> _add($x, $class -> _one());
349}
350
351sub _dec {
352    my ($class, $x) = @_;
353    $class -> _sub($x, $class -> _one());
354}
355
356# Signed addition. If the flag is false, $xa might be modified, but not $ya. If
357# the false is true, $ya might be modified, but not $xa.
358
359sub _sadd {
360    my $class = shift;
361    my ($xa, $xs, $ya, $ys, $flag) = @_;
362    my ($za, $zs);
363
364    # If the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
365
366    if ($xs eq $ys) {
367        if ($flag) {
368            $za = $class -> _add($ya, $xa);
369        } else {
370            $za = $class -> _add($xa, $ya);
371        }
372        $zs = $class -> _is_zero($za) ? '+' : $xs;
373        return $za, $zs;
374    }
375
376    my $acmp = $class -> _acmp($xa, $ya);       # abs(x) = abs(y)
377
378    if ($acmp == 0) {                           # x = -y or -x = y
379        $za = $class -> _zero();
380        $zs = '+';
381        return $za, $zs;
382    }
383
384    if ($acmp > 0) {                            # abs(x) > abs(y)
385        $za = $class -> _sub($xa, $ya, $flag);
386        $zs = $xs;
387    } else {                                    # abs(x) < abs(y)
388        $za = $class -> _sub($ya, $xa, !$flag);
389        $zs = $ys;
390    }
391    return $za, $zs;
392}
393
394# Signed subtraction. If the flag is false, $xa might be modified, but not $ya.
395# If the false is true, $ya might be modified, but not $xa.
396
397sub _ssub {
398    my $class = shift;
399    my ($xa, $xs, $ya, $ys, $flag) = @_;
400
401    # Swap sign of second operand and let _sadd() do the job.
402    $ys = $ys eq '+' ? '-' : '+';
403    $class -> _sadd($xa, $xs, $ya, $ys, $flag);
404}
405
406##############################################################################
407# testing
408
409sub _acmp {
410    # Compare two (absolute) values. Return -1, 0, or 1.
411    my ($class, $x, $y) = @_;
412    my $xstr = $class -> _str($x);
413    my $ystr = $class -> _str($y);
414
415    length($xstr) <=> length($ystr) || $xstr cmp $ystr;
416}
417
418sub _len {
419    my ($class, $x) = @_;
420    CORE::length($class -> _str($x));
421}
422
423sub _alen {
424    my ($class, $x) = @_;
425    $class -> _len($x);
426}
427
428sub _digit {
429    my ($class, $x, $n) = @_;
430    substr($class ->_str($x), -($n+1), 1);
431}
432
433sub _digitsum {
434    my ($class, $x) = @_;
435
436    my $len = $class -> _len($x);
437    my $sum = $class -> _zero();
438    for (my $i = 0 ; $i < $len ; ++$i) {
439        my $digit = $class -> _digit($x, $i);
440        $digit = $class -> _new($digit);
441        $sum = $class -> _add($sum, $digit);
442    }
443
444    return $sum;
445}
446
447sub _zeros {
448    my ($class, $x) = @_;
449    my $str = $class -> _str($x);
450    $str =~ /[^0](0*)\z/ ? CORE::length($1) : 0;
451}
452
453##############################################################################
454# _is_* routines
455
456sub _is_zero {
457    # return true if arg is zero
458    my ($class, $x) = @_;
459    $class -> _str($x) == 0;
460}
461
462sub _is_even {
463    # return true if arg is even
464    my ($class, $x) = @_;
465    substr($class -> _str($x), -1, 1) % 2 == 0;
466}
467
468sub _is_odd {
469    # return true if arg is odd
470    my ($class, $x) = @_;
471    substr($class -> _str($x), -1, 1) % 2 != 0;
472}
473
474sub _is_one {
475    # return true if arg is one
476    my ($class, $x) = @_;
477    $class -> _str($x) == 1;
478}
479
480sub _is_two {
481    # return true if arg is two
482    my ($class, $x) = @_;
483    $class -> _str($x) == 2;
484}
485
486sub _is_ten {
487    # return true if arg is ten
488    my ($class, $x) = @_;
489    $class -> _str($x) == 10;
490}
491
492###############################################################################
493# check routine to test internal state for corruptions
494
495sub _check {
496    # used by the test suite
497    my ($class, $x) = @_;
498    return "Input is undefined" unless defined $x;
499    return "$x is not a reference" unless ref($x);
500    return 0;
501}
502
503###############################################################################
504
505sub _mod {
506    # modulus
507    my ($class, $x, $y) = @_;
508
509    croak "@{[(caller 0)[3]]} requires non-zero second operand"
510      if $class -> _is_zero($y);
511
512    if ($class -> can('_div')) {
513        $x = $class -> _copy($x);
514        my ($q, $r) = $class -> _div($x, $y);
515        return $r;
516    } else {
517        my $r = $class -> _copy($x);
518        while ($class -> _acmp($r, $y) >= 0) {
519            $r = $class -> _sub($r, $y);
520        }
521        return $r;
522    }
523}
524
525##############################################################################
526# shifts
527
528sub _rsft {
529    my ($class, $x, $n, $b) = @_;
530    $b = $class -> _new($b) unless ref $b;
531    return scalar $class -> _div($x, $class -> _pow($class -> _copy($b), $n));
532}
533
534sub _lsft {
535    my ($class, $x, $n, $b) = @_;
536    $b = $class -> _new($b) unless ref $b;
537    return $class -> _mul($x, $class -> _pow($class -> _copy($b), $n));
538}
539
540sub _pow {
541    # power of $x to $y
542    my ($class, $x, $y) = @_;
543
544    if ($class -> _is_zero($y)) {
545        return $class -> _one();        # y == 0 => x => 1
546    }
547
548    if (($class -> _is_one($x)) ||      #    x == 1
549        ($class -> _is_one($y)))        # or y == 1
550    {
551        return $x;
552    }
553
554    if ($class -> _is_zero($x)) {
555        return $class -> _zero();       # 0 ** y => 0 (if not y <= 0)
556    }
557
558    my $pow2 = $class -> _one();
559
560    my $y_bin = $class -> _as_bin($y);
561    $y_bin =~ s/^0b//;
562    my $len = length($y_bin);
563
564    while (--$len > 0) {
565        $pow2 = $class -> _mul($pow2, $x) if substr($y_bin, $len, 1) eq '1';
566        $x = $class -> _mul($x, $x);
567    }
568
569    $x = $class -> _mul($x, $pow2);
570    return $x;
571}
572
573sub _nok {
574    # Return binomial coefficient (n over k).
575    my ($class, $n, $k) = @_;
576
577    # If k > n/2, or, equivalently, 2*k > n, compute nok(n, k) as
578    # nok(n, n-k), to minimize the number if iterations in the loop.
579
580    {
581        my $twok = $class -> _mul($class -> _two(), $class -> _copy($k));
582        if ($class -> _acmp($twok, $n) > 0) {
583            $k = $class -> _sub($class -> _copy($n), $k);
584        }
585    }
586
587    # Example:
588    #
589    # / 7 \       7!       1*2*3*4 * 5*6*7   5 * 6 * 7
590    # |   | = --------- =  --------------- = --------- = ((5 * 6) / 2 * 7) / 3
591    # \ 3 /   (7-3)! 3!    1*2*3*4 * 1*2*3   1 * 2 * 3
592    #
593    # Equivalently, _nok(11, 5) is computed as
594    #
595    # (((((((7 * 8) / 2) * 9) / 3) * 10) / 4) * 11) / 5
596
597    if ($class -> _is_zero($k)) {
598        return $class -> _one();
599    }
600
601    # Make a copy of the original n, in case the subclass modifies n in-place.
602
603    my $n_orig = $class -> _copy($n);
604
605    # n = 5, f = 6, d = 2 (cf. example above)
606
607    $n = $class -> _sub($n, $k);
608    $n = $class -> _inc($n);
609
610    my $f = $class -> _copy($n);
611    $f = $class -> _inc($f);
612
613    my $d = $class -> _two();
614
615    # while f <= n (the original n, that is) ...
616
617    while ($class -> _acmp($f, $n_orig) <= 0) {
618        $n = $class -> _mul($n, $f);
619        $n = $class -> _div($n, $d);
620        $f = $class -> _inc($f);
621        $d = $class -> _inc($d);
622    }
623
624    return $n;
625}
626
627#sub _fac {
628#    # factorial
629#    my ($class, $x) = @_;
630#
631#    my $two = $class -> _two();
632#
633#    if ($class -> _acmp($x, $two) < 0) {
634#        return $class -> _one();
635#    }
636#
637#    my $i = $class -> _copy($x);
638#    while ($class -> _acmp($i, $two) > 0) {
639#        $i = $class -> _dec($i);
640#        $x = $class -> _mul($x, $i);
641#    }
642#
643#    return $x;
644#}
645
646sub _fac {
647    # factorial
648    my ($class, $x) = @_;
649
650    # This is an implementation of the split recursive algorithm. See
651    # http://www.luschny.de/math/factorial/csharp/FactorialSplit.cs.html
652
653    my $p   = $class -> _one();
654    my $r   = $class -> _one();
655    my $two = $class -> _two();
656
657    my ($log2n) = $class -> _log_int($class -> _copy($x), $two);
658    my $h     = $class -> _zero();
659    my $shift = $class -> _zero();
660    my $k     = $class -> _one();
661
662    while ($class -> _acmp($h, $x)) {
663        $shift = $class -> _add($shift, $h);
664        $h = $class -> _rsft($class -> _copy($x), $log2n, $two);
665        $log2n = $class -> _dec($log2n) if !$class -> _is_zero($log2n);
666        my $high = $class -> _copy($h);
667        $high = $class -> _dec($high) if $class -> _is_even($h);
668        while ($class -> _acmp($k, $high)) {
669            $k = $class -> _add($k, $two);
670            $p = $class -> _mul($p, $k);
671        }
672        $r = $class -> _mul($r, $p);
673    }
674    return $class -> _lsft($r, $shift, $two);
675}
676
677sub _dfac {
678    # double factorial
679    my ($class, $x) = @_;
680
681    my $two = $class -> _two();
682
683    if ($class -> _acmp($x, $two) < 0) {
684        return $class -> _one();
685    }
686
687    my $i = $class -> _copy($x);
688    while ($class -> _acmp($i, $two) > 0) {
689        $i = $class -> _sub($i, $two);
690        $x = $class -> _mul($x, $i);
691    }
692
693    return $x;
694}
695
696sub _log_int {
697    # calculate integer log of $x to base $base
698    # calculate integer log of $x to base $base
699    # ref to array, ref to array - return ref to array
700    my ($class, $x, $base) = @_;
701
702    # X == 0 => NaN
703    return if $class -> _is_zero($x);
704
705    $base = $class -> _new(2)     unless defined($base);
706    $base = $class -> _new($base) unless ref($base);
707
708    # BASE 0 or 1 => NaN
709    return if $class -> _is_zero($base) || $class -> _is_one($base);
710
711    # X == 1 => 0 (is exact)
712    if ($class -> _is_one($x)) {
713        return $class -> _zero(), 1 if wantarray;
714        return $class -> _zero();
715    }
716
717    my $cmp = $class -> _acmp($x, $base);
718
719    # X == BASE => 1 (is exact)
720    if ($cmp == 0) {
721        return $class -> _one(), 1 if wantarray;
722        return $class -> _one();
723    }
724
725    # 1 < X < BASE => 0 (is truncated)
726    if ($cmp < 0) {
727        return $class -> _zero(), 0 if wantarray;
728        return $class -> _zero();
729    }
730
731    my $y;
732
733    # log(x) / log(b) = log(xm * 10^xe) / log(bm * 10^be)
734    #                 = (log(xm) + xe*(log(10))) / (log(bm) + be*log(10))
735
736    {
737        my $x_str = $class -> _str($x);
738        my $b_str = $class -> _str($base);
739        my $xm    = "." . $x_str;
740        my $bm    = "." . $b_str;
741        my $xe    = length($x_str);
742        my $be    = length($b_str);
743        my $log10 = log(10);
744        my $guess = int((log($xm) + $xe * $log10) / (log($bm) + $be * $log10));
745        $y = $class -> _new($guess);
746    }
747
748    my $trial = $class -> _pow($class -> _copy($base), $y);
749    my $acmp  = $class -> _acmp($trial, $x);
750
751    # Too small?
752
753    while ($acmp < 0) {
754        $trial = $class -> _mul($trial, $base);
755        $y     = $class -> _inc($y);
756        $acmp  = $class -> _acmp($trial, $x);
757    }
758
759    # Too big?
760
761    while ($acmp > 0) {
762        $trial = $class -> _div($trial, $base);
763        $y     = $class -> _dec($y);
764        $acmp  = $class -> _acmp($trial, $x);
765    }
766
767    return wantarray ? ($y, 1) : $y if $acmp == 0;      # result is exact
768    return wantarray ? ($y, 0) : $y;                    # result is too small
769}
770
771sub _ilog2 {
772    my ($class, $x) = @_;
773
774    return if $class -> _is_zero($x);
775
776    my $str = $class -> _to_hex($x);
777
778    # First do the bits in all but the most significant hex digit.
779
780    my $y = $class -> _new(length($str) - 1);
781    $y = $class -> _mul($y, $class -> _new(4));
782
783    # Now add the number of bits in the most significant hex digit.
784
785    my $n = int log(hex(substr($str, 0, 1))) / log(2);
786    $y = $class -> _add($y, $class -> _new($n));
787    return $y unless wantarray;
788
789    my $pow2 = $class -> _lsft($class -> _one(), $y, 2);
790    my $is_exact = $class -> _acmp($x, $pow2) == 0 ? 1 : 0;
791    return $y, $is_exact;
792}
793
794sub _ilog10 {
795    my ($class, $x) = @_;
796
797    return if $class -> _is_zero($x);
798
799    my $str = $class -> _str($x);
800    my $len = length($str);
801    my $y = $class -> _new($len - 1);
802    return $y unless wantarray;
803
804    #my $pow10 = $class -> _1ex($y);
805    #my $is_exact = $class -> _acmp($x, $pow10) ? 1 : 0;
806
807    my $is_exact = $str =~ /^10*$/ ? 1 : 0;
808    return $y, $is_exact;
809}
810
811sub _clog2 {
812    my ($class, $x) = @_;
813
814    return if $class -> _is_zero($x);
815
816    my $str = $class -> _to_hex($x);
817
818    # First do the bits in all but the most significant hex digit.
819
820    my $y = $class -> _new(length($str) - 1);
821    $y = $class -> _mul($y, $class -> _new(4));
822
823    # Now add the number of bits in the most significant hex digit.
824
825    my $n = int log(hex(substr($str, 0, 1))) / log(2);
826    $y = $class -> _add($y, $class -> _new($n));
827
828    # $y is now 1 too small unless $y is an exact power of 2.
829
830    my $pow2 = $class -> _lsft($class -> _one(), $y, 2);
831    my $is_exact = $class -> _acmp($x, $pow2) == 0 ? 1 : 0;
832    $y = $class -> _inc($y) if $is_exact == 0;
833    return $y, $is_exact if wantarray;
834    return $y;
835}
836
837sub _clog10 {
838    my ($class, $x) = @_;
839
840    return if $class -> _is_zero($x);
841
842    my $str = $class -> _str($x);
843    my $len = length($str);
844
845    if ($str =~ /^10*$/) {
846        my $y = $class -> _new($len - 1);
847        return $y, 1 if wantarray;
848        return $y;
849    }
850
851    my $y = $class -> _new($len);
852    return $y, 0 if wantarray;
853    return $y;
854}
855
856sub _sqrt {
857    # square-root of $y in place
858    my ($class, $y) = @_;
859
860    return $y if $class -> _is_zero($y);
861
862    my $y_str = $class -> _str($y);
863    my $y_len = length($y_str);
864
865    # Compute the guess $x.
866
867    my $xm;
868    my $xe;
869    if ($y_len % 2 == 0) {
870        $xm = sqrt("." . $y_str);
871        $xe = $y_len / 2;
872        $xm = sprintf "%.0f", int($xm * 1e15);
873        $xe -= 15;
874    } else {
875        $xm = sqrt(".0" . $y_str);
876        $xe = ($y_len + 1) / 2;
877        $xm = sprintf "%.0f", int($xm * 1e16);
878        $xe -= 16;
879    }
880
881    my $x;
882    if ($xe < 0) {
883        $x = substr $xm, 0, length($xm) + $xe;
884    } else {
885        $x = $xm . ("0" x $xe);
886    }
887
888    $x = $class -> _new($x);
889
890    # Newton's method for computing square root of y
891    #
892    # x(i+1) = x(i) - f(x(i)) / f'(x(i))
893    #        = x(i) - (x(i)^2 - y) / (2 * x(i))     # use if x(i)^2 > y
894    #        = x(i) + (y - x(i)^2) / (2 * x(i))     # use if x(i)^2 < y
895
896    # Determine if x, our guess, is too small, correct, or too large.
897
898    my $xsq = $class -> _mul($class -> _copy($x), $x);          # x(i)^2
899    my $acmp = $class -> _acmp($xsq, $y);                       # x(i)^2 <=> y
900
901    # Only assign a value to this variable if we will be using it.
902
903    my $two;
904    $two = $class -> _two() if $acmp != 0;
905
906    # If x is too small, do one iteration of Newton's method. Since the
907    # function f(x) = x^2 - y is concave and monotonically increasing, the next
908    # guess for x will either be correct or too large.
909
910    if ($acmp < 0) {
911
912        # x(i+1) = x(i) + (y - x(i)^2) / (2 * x(i))
913
914        my $numer = $class -> _sub($class -> _copy($y), $xsq);  # y - x(i)^2
915        my $denom = $class -> _mul($class -> _copy($two), $x);  # 2 * x(i)
916        my $delta = $class -> _div($numer, $denom);
917
918        unless ($class -> _is_zero($delta)) {
919            $x    = $class -> _add($x, $delta);
920            $xsq  = $class -> _mul($class -> _copy($x), $x);    # x(i)^2
921            $acmp = $class -> _acmp($xsq, $y);                  # x(i)^2 <=> y
922        }
923    }
924
925    # If our guess for x is too large, apply Newton's method repeatedly until
926    # we either have got the correct value, or the delta is zero.
927
928    while ($acmp > 0) {
929
930        # x(i+1) = x(i) - (x(i)^2 - y) / (2 * x(i))
931
932        my $numer = $class -> _sub($xsq, $y);                   # x(i)^2 - y
933        my $denom = $class -> _mul($class -> _copy($two), $x);  # 2 * x(i)
934        my $delta = $class -> _div($numer, $denom);
935        last if $class -> _is_zero($delta);
936
937        $x    = $class -> _sub($x, $delta);
938        $xsq  = $class -> _mul($class -> _copy($x), $x);        # x(i)^2
939        $acmp = $class -> _acmp($xsq, $y);                      # x(i)^2 <=> y
940    }
941
942    # When the delta is zero, our value for x might still be too large. We
943    # require that the outout is either exact or too small (i.e., rounded down
944    # to the nearest integer), so do a final check.
945
946    while ($acmp > 0) {
947        $x    = $class -> _dec($x);
948        $xsq  = $class -> _mul($class -> _copy($x), $x);        # x(i)^2
949        $acmp = $class -> _acmp($xsq, $y);                      # x(i)^2 <=> y
950    }
951
952    return $x;
953}
954
955sub _root {
956    my ($class, $y, $n) = @_;
957
958    return $y if $class -> _is_zero($y) || $class -> _is_one($y) ||
959                 $class -> _is_one($n);
960
961    # If y <= n, the result is always (truncated to) 1.
962
963    return $class -> _one() if $class -> _acmp($y, $n) <= 0;
964
965    # Compute the initial guess x of y^(1/n). When n is large, Newton's method
966    # converges slowly if the "guess" (initial value) is poor, so we need a
967    # good guess. It the guess is too small, the next guess will be too large,
968    # and from then on all guesses are too large.
969
970    my $DEBUG = 0;
971
972    # Split y into mantissa and exponent in base 10, so that
973    #
974    #   y = xm * 10^xe, where 0 < xm < 1 and xe is an integer
975
976    my $y_str  = $class -> _str($y);
977    my $ym = "." . $y_str;
978    my $ye = length($y_str);
979
980    # From this compute the approximate base 10 logarithm of y
981    #
982    #   log_10(y) = log_10(ym) + log_10(ye^10)
983    #             = log(ym)/log(10) + ye
984
985    my $log10y = log($ym) / log(10) + $ye;
986
987    # And from this compute the approximate base 10 logarithm of x, where
988    # x = y^(1/n)
989    #
990    #   log_10(x) = log_10(y)/n
991
992    my $log10x = $log10y / $class -> _num($n);
993
994    # From this compute xm and xe, the mantissa and exponent (in base 10) of x,
995    # where 1 < xm <= 10 and xe is an integer.
996
997    my $xe = int $log10x;
998    my $xm = 10 ** ($log10x - $xe);
999
1000    # Scale the mantissa and exponent to increase the integer part of ym, which
1001    # gives us better accuracy.
1002
1003    if ($DEBUG) {
1004        print "\n";
1005        print "y_str  = $y_str\n";
1006        print "ym     = $ym\n";
1007        print "ye     = $ye\n";
1008        print "log10y = $log10y\n";
1009        print "log10x = $log10x\n";
1010        print "xm     = $xm\n";
1011        print "xe     = $xe\n";
1012    }
1013
1014    my $d = $xe < 15 ? $xe : 15;
1015    $xm *= 10 ** $d;
1016    $xe -= $d;
1017
1018    if ($DEBUG) {
1019        print "\n";
1020        print "xm     = $xm\n";
1021        print "xe     = $xe\n";
1022    }
1023
1024    # If the mantissa is not an integer, round up to nearest integer, and then
1025    # convert the number to a string. It is important to always round up due to
1026    # how Newton's method behaves in this case. If the initial guess is too
1027    # small, the next guess will be too large, after which every succeeding
1028    # guess converges the correct value from above. Now, if the initial guess
1029    # is too small and n is large, the next guess will be much too large and
1030    # require a large number of iterations to get close to the solution.
1031    # Because of this, we are likely to find the solution faster if we make
1032    # sure the initial guess is not too small.
1033
1034    my $xm_int = int($xm);
1035    my $x_str = sprintf '%.0f', $xm > $xm_int ? $xm_int + 1 : $xm_int;
1036    $x_str .= "0" x $xe;
1037
1038    my $x = $class -> _new($x_str);
1039
1040    if ($DEBUG) {
1041        print "xm     = $xm\n";
1042        print "xe     = $xe\n";
1043        print "\n";
1044        print "x_str  = $x_str (initial guess)\n";
1045        print "\n";
1046    }
1047
1048    # Use Newton's method for computing n'th root of y.
1049    #
1050    # x(i+1) = x(i) - f(x(i)) / f'(x(i))
1051    #        = x(i) - (x(i)^n - y) / (n * x(i)^(n-1))   # use if x(i)^n > y
1052    #        = x(i) + (y - x(i)^n) / (n * x(i)^(n-1))   # use if x(i)^n < y
1053
1054    # Determine if x, our guess, is too small, correct, or too large. Rather
1055    # than computing x(i)^n and x(i)^(n-1) directly, compute x(i)^(n-1) and
1056    # then the same value multiplied by x.
1057
1058    my $nm1     = $class -> _dec($class -> _copy($n));           # n-1
1059    my $xpownm1 = $class -> _pow($class -> _copy($x), $nm1);     # x(i)^(n-1)
1060    my $xpown   = $class -> _mul($class -> _copy($xpownm1), $x); # x(i)^n
1061    my $acmp    = $class -> _acmp($xpown, $y);                   # x(i)^n <=> y
1062
1063    if ($DEBUG) {
1064        print "\n";
1065        print "x      = ", $class -> _str($x), "\n";
1066        print "x^n    = ", $class -> _str($xpown), "\n";
1067        print "y      = ", $class -> _str($y), "\n";
1068        print "acmp   = $acmp\n";
1069    }
1070
1071    # If x is too small, do one iteration of Newton's method. Since the
1072    # function f(x) = x^n - y is concave and monotonically increasing, the next
1073    # guess for x will either be correct or too large.
1074
1075    if ($acmp < 0) {
1076
1077        # x(i+1) = x(i) + (y - x(i)^n) / (n * x(i)^(n-1))
1078
1079        my $numer = $class -> _sub($class -> _copy($y), $xpown);    # y - x(i)^n
1080        my $denom = $class -> _mul($class -> _copy($n), $xpownm1);  # n * x(i)^(n-1)
1081        my $delta = $class -> _div($numer, $denom);
1082
1083        if ($DEBUG) {
1084            print "\n";
1085            print "numer  = ", $class -> _str($numer), "\n";
1086            print "denom  = ", $class -> _str($denom), "\n";
1087            print "delta  = ", $class -> _str($delta), "\n";
1088        }
1089
1090        unless ($class -> _is_zero($delta)) {
1091            $x       = $class -> _add($x, $delta);
1092            $xpownm1 = $class -> _pow($class -> _copy($x), $nm1);     # x(i)^(n-1)
1093            $xpown   = $class -> _mul($class -> _copy($xpownm1), $x); # x(i)^n
1094            $acmp    = $class -> _acmp($xpown, $y);                   # x(i)^n <=> y
1095
1096            if ($DEBUG) {
1097                print "\n";
1098                print "x      = ", $class -> _str($x), "\n";
1099                print "x^n    = ", $class -> _str($xpown), "\n";
1100                print "y      = ", $class -> _str($y), "\n";
1101                print "acmp   = $acmp\n";
1102            }
1103        }
1104    }
1105
1106    # If our guess for x is too large, apply Newton's method repeatedly until
1107    # we either have got the correct value, or the delta is zero.
1108
1109    while ($acmp > 0) {
1110
1111        # x(i+1) = x(i) - (x(i)^n - y) / (n * x(i)^(n-1))
1112
1113        my $numer = $class -> _sub($class -> _copy($xpown), $y);    # x(i)^n - y
1114        my $denom = $class -> _mul($class -> _copy($n), $xpownm1);  # n * x(i)^(n-1)
1115
1116        if ($DEBUG) {
1117            print "numer  = ", $class -> _str($numer), "\n";
1118            print "denom  = ", $class -> _str($denom), "\n";
1119        }
1120
1121        my $delta = $class -> _div($numer, $denom);
1122
1123        if ($DEBUG) {
1124            print "delta  = ", $class -> _str($delta), "\n";
1125        }
1126
1127        last if $class -> _is_zero($delta);
1128
1129        $x       = $class -> _sub($x, $delta);
1130        $xpownm1 = $class -> _pow($class -> _copy($x), $nm1);     # x(i)^(n-1)
1131        $xpown   = $class -> _mul($class -> _copy($xpownm1), $x); # x(i)^n
1132        $acmp    = $class -> _acmp($xpown, $y);                   # x(i)^n <=> y
1133
1134        if ($DEBUG) {
1135            print "\n";
1136            print "x      = ", $class -> _str($x), "\n";
1137            print "x^n    = ", $class -> _str($xpown), "\n";
1138            print "y      = ", $class -> _str($y), "\n";
1139            print "acmp   = $acmp\n";
1140        }
1141    }
1142
1143    # When the delta is zero, our value for x might still be too large. We
1144    # require that the outout is either exact or too small (i.e., rounded down
1145    # to the nearest integer), so do a final check.
1146
1147    while ($acmp > 0) {
1148        $x     = $class -> _dec($x);
1149        $xpown = $class -> _pow($class -> _copy($x), $n);     # x(i)^n
1150        $acmp  = $class -> _acmp($xpown, $y);                 # x(i)^n <=> y
1151    }
1152
1153    return $x;
1154}
1155
1156##############################################################################
1157# binary stuff
1158
1159sub _and {
1160    my ($class, $x, $y) = @_;
1161
1162    return $x if $class -> _acmp($x, $y) == 0;
1163
1164    my $m    = $class -> _one();
1165    my $mask = $class -> _new("32768");
1166
1167    my ($xr, $yr);                # remainders after division
1168
1169    my $xc = $class -> _copy($x);
1170    my $yc = $class -> _copy($y);
1171    my $z  = $class -> _zero();
1172
1173    until ($class -> _is_zero($xc) || $class -> _is_zero($yc)) {
1174        ($xc, $xr) = $class -> _div($xc, $mask);
1175        ($yc, $yr) = $class -> _div($yc, $mask);
1176        my $bits = $class -> _new($class -> _num($xr) & $class -> _num($yr));
1177        $z = $class -> _add($z, $class -> _mul($bits, $m));
1178        $m = $class -> _mul($m, $mask);
1179    }
1180
1181    return $z;
1182}
1183
1184sub _xor {
1185    my ($class, $x, $y) = @_;
1186
1187    return $class -> _zero() if $class -> _acmp($x, $y) == 0;
1188
1189    my $m    = $class -> _one();
1190    my $mask = $class -> _new("32768");
1191
1192    my ($xr, $yr);                # remainders after division
1193
1194    my $xc = $class -> _copy($x);
1195    my $yc = $class -> _copy($y);
1196    my $z  = $class -> _zero();
1197
1198    until ($class -> _is_zero($xc) || $class -> _is_zero($yc)) {
1199        ($xc, $xr) = $class -> _div($xc, $mask);
1200        ($yc, $yr) = $class -> _div($yc, $mask);
1201        my $bits = $class -> _new($class -> _num($xr) ^ $class -> _num($yr));
1202        $z = $class -> _add($z, $class -> _mul($bits, $m));
1203        $m = $class -> _mul($m, $mask);
1204    }
1205
1206    # The loop above stops when the smallest of the two numbers is exhausted.
1207    # The remainder of the longer one will survive bit-by-bit, so we simple
1208    # multiply-add it in.
1209
1210    $z = $class -> _add($z, $class -> _mul($xc, $m))
1211      unless $class -> _is_zero($xc);
1212    $z = $class -> _add($z, $class -> _mul($yc, $m))
1213      unless $class -> _is_zero($yc);
1214
1215    return $z;
1216}
1217
1218sub _or {
1219    my ($class, $x, $y) = @_;
1220
1221    return $x if $class -> _acmp($x, $y) == 0; # shortcut (see _and)
1222
1223    my $m    = $class -> _one();
1224    my $mask = $class -> _new("32768");
1225
1226    my ($xr, $yr);                # remainders after division
1227
1228    my $xc = $class -> _copy($x);
1229    my $yc = $class -> _copy($y);
1230    my $z  = $class -> _zero();
1231
1232    until ($class -> _is_zero($xc) || $class -> _is_zero($yc)) {
1233        ($xc, $xr) = $class -> _div($xc, $mask);
1234        ($yc, $yr) = $class -> _div($yc, $mask);
1235        my $bits = $class -> _new($class -> _num($xr) | $class -> _num($yr));
1236        $z = $class -> _add($z, $class -> _mul($bits, $m));
1237        $m = $class -> _mul($m, $mask);
1238    }
1239
1240    # The loop above stops when the smallest of the two numbers is exhausted.
1241    # The remainder of the longer one will survive bit-by-bit, so we simple
1242    # multiply-add it in.
1243
1244    $z = $class -> _add($z, $class -> _mul($xc, $m))
1245      unless $class -> _is_zero($xc);
1246    $z = $class -> _add($z, $class -> _mul($yc, $m))
1247      unless $class -> _is_zero($yc);
1248
1249    return $z;
1250}
1251
1252sub _sand {
1253    my ($class, $x, $sx, $y, $sy) = @_;
1254
1255    return ($class -> _zero(), '+')
1256      if $class -> _is_zero($x) || $class -> _is_zero($y);
1257
1258    my $sign = $sx eq '-' && $sy eq '-' ? '-' : '+';
1259
1260    my ($bx, $by);
1261
1262    if ($sx eq '-') {                   # if x is negative
1263        # two's complement: inc (dec unsigned value) and flip all "bits" in $bx
1264        $bx = $class -> _copy($x);
1265        $bx = $class -> _dec($bx);
1266        $bx = $class -> _as_hex($bx);
1267        $bx =~ s/^-?0x//;
1268        $bx =~ tr<0123456789abcdef>
1269                <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1270    } else {                            # if x is positive
1271        $bx = $class -> _as_hex($x);    # get binary representation
1272        $bx =~ s/^-?0x//;
1273        $bx =~ tr<fedcba9876543210>
1274                 <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1275    }
1276
1277    if ($sy eq '-') {                   # if y is negative
1278        # two's complement: inc (dec unsigned value) and flip all "bits" in $by
1279        $by = $class -> _copy($y);
1280        $by = $class -> _dec($by);
1281        $by = $class -> _as_hex($by);
1282        $by =~ s/^-?0x//;
1283        $by =~ tr<0123456789abcdef>
1284                <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1285    } else {
1286        $by = $class -> _as_hex($y);    # get binary representation
1287        $by =~ s/^-?0x//;
1288        $by =~ tr<fedcba9876543210>
1289                <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1290    }
1291
1292    # now we have bit-strings from X and Y, reverse them for padding
1293    $bx = reverse $bx;
1294    $by = reverse $by;
1295
1296    # padd the shorter string
1297    my $xx = "\x00"; $xx = "\x0f" if $sx eq '-';
1298    my $yy = "\x00"; $yy = "\x0f" if $sy eq '-';
1299    my $diff = CORE::length($bx) - CORE::length($by);
1300    if ($diff > 0) {
1301        # if $yy eq "\x00", we can cut $bx, otherwise we need to padd $by
1302        $by .= $yy x $diff;
1303    } elsif ($diff < 0) {
1304        # if $xx eq "\x00", we can cut $by, otherwise we need to padd $bx
1305        $bx .= $xx x abs($diff);
1306    }
1307
1308    # and the strings together
1309    my $r = $bx & $by;
1310
1311    # and reverse the result again
1312    $bx = reverse $r;
1313
1314    # One of $bx or $by was negative, so need to flip bits in the result. In both
1315    # cases (one or two of them negative, or both positive) we need to get the
1316    # characters back.
1317    if ($sign eq '-') {
1318        $bx =~ tr<\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>
1319                 <0123456789abcdef>;
1320    } else {
1321        $bx =~ tr<\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>
1322                 <fedcba9876543210>;
1323    }
1324
1325    # leading zeros will be stripped by _from_hex()
1326    $bx = '0x' . $bx;
1327    $bx = $class -> _from_hex($bx);
1328
1329    $bx = $class -> _inc($bx) if $sign eq '-';
1330
1331    # avoid negative zero
1332    $sign = '+' if $class -> _is_zero($bx);
1333
1334    return $bx, $sign;
1335}
1336
1337sub _sxor {
1338    my ($class, $x, $sx, $y, $sy) = @_;
1339
1340    return ($class -> _zero(), '+')
1341      if $class -> _is_zero($x) && $class -> _is_zero($y);
1342
1343    my $sign = $sx ne $sy ? '-' : '+';
1344
1345    my ($bx, $by);
1346
1347    if ($sx eq '-') {                   # if x is negative
1348        # two's complement: inc (dec unsigned value) and flip all "bits" in $bx
1349        $bx = $class -> _copy($x);
1350        $bx = $class -> _dec($bx);
1351        $bx = $class -> _as_hex($bx);
1352        $bx =~ s/^-?0x//;
1353        $bx =~ tr<0123456789abcdef>
1354                <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1355    } else {                            # if x is positive
1356        $bx = $class -> _as_hex($x);    # get binary representation
1357        $bx =~ s/^-?0x//;
1358        $bx =~ tr<fedcba9876543210>
1359                 <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1360    }
1361
1362    if ($sy eq '-') {                   # if y is negative
1363        # two's complement: inc (dec unsigned value) and flip all "bits" in $by
1364        $by = $class -> _copy($y);
1365        $by = $class -> _dec($by);
1366        $by = $class -> _as_hex($by);
1367        $by =~ s/^-?0x//;
1368        $by =~ tr<0123456789abcdef>
1369                <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1370    } else {
1371        $by = $class -> _as_hex($y);    # get binary representation
1372        $by =~ s/^-?0x//;
1373        $by =~ tr<fedcba9876543210>
1374                <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1375    }
1376
1377    # now we have bit-strings from X and Y, reverse them for padding
1378    $bx = reverse $bx;
1379    $by = reverse $by;
1380
1381    # padd the shorter string
1382    my $xx = "\x00"; $xx = "\x0f" if $sx eq '-';
1383    my $yy = "\x00"; $yy = "\x0f" if $sy eq '-';
1384    my $diff = CORE::length($bx) - CORE::length($by);
1385    if ($diff > 0) {
1386        # if $yy eq "\x00", we can cut $bx, otherwise we need to padd $by
1387        $by .= $yy x $diff;
1388    } elsif ($diff < 0) {
1389        # if $xx eq "\x00", we can cut $by, otherwise we need to padd $bx
1390        $bx .= $xx x abs($diff);
1391    }
1392
1393    # xor the strings together
1394    my $r = $bx ^ $by;
1395
1396    # and reverse the result again
1397    $bx = reverse $r;
1398
1399    # One of $bx or $by was negative, so need to flip bits in the result. In both
1400    # cases (one or two of them negative, or both positive) we need to get the
1401    # characters back.
1402    if ($sign eq '-') {
1403        $bx =~ tr<\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>
1404                 <0123456789abcdef>;
1405    } else {
1406        $bx =~ tr<\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>
1407                 <fedcba9876543210>;
1408    }
1409
1410    # leading zeros will be stripped by _from_hex()
1411    $bx = '0x' . $bx;
1412    $bx = $class -> _from_hex($bx);
1413
1414    $bx = $class -> _inc($bx) if $sign eq '-';
1415
1416    # avoid negative zero
1417    $sign = '+' if $class -> _is_zero($bx);
1418
1419    return $bx, $sign;
1420}
1421
1422sub _sor {
1423    my ($class, $x, $sx, $y, $sy) = @_;
1424
1425    return ($class -> _zero(), '+')
1426      if $class -> _is_zero($x) && $class -> _is_zero($y);
1427
1428    my $sign = $sx eq '-' || $sy eq '-' ? '-' : '+';
1429
1430    my ($bx, $by);
1431
1432    if ($sx eq '-') {                   # if x is negative
1433        # two's complement: inc (dec unsigned value) and flip all "bits" in $bx
1434        $bx = $class -> _copy($x);
1435        $bx = $class -> _dec($bx);
1436        $bx = $class -> _as_hex($bx);
1437        $bx =~ s/^-?0x//;
1438        $bx =~ tr<0123456789abcdef>
1439                <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1440    } else {                            # if x is positive
1441        $bx = $class -> _as_hex($x);     # get binary representation
1442        $bx =~ s/^-?0x//;
1443        $bx =~ tr<fedcba9876543210>
1444                 <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1445    }
1446
1447    if ($sy eq '-') {                   # if y is negative
1448        # two's complement: inc (dec unsigned value) and flip all "bits" in $by
1449        $by = $class -> _copy($y);
1450        $by = $class -> _dec($by);
1451        $by = $class -> _as_hex($by);
1452        $by =~ s/^-?0x//;
1453        $by =~ tr<0123456789abcdef>
1454                <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1455    } else {
1456        $by = $class -> _as_hex($y);     # get binary representation
1457        $by =~ s/^-?0x//;
1458        $by =~ tr<fedcba9876543210>
1459                <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1460    }
1461
1462    # now we have bit-strings from X and Y, reverse them for padding
1463    $bx = reverse $bx;
1464    $by = reverse $by;
1465
1466    # padd the shorter string
1467    my $xx = "\x00"; $xx = "\x0f" if $sx eq '-';
1468    my $yy = "\x00"; $yy = "\x0f" if $sy eq '-';
1469    my $diff = CORE::length($bx) - CORE::length($by);
1470    if ($diff > 0) {
1471        # if $yy eq "\x00", we can cut $bx, otherwise we need to padd $by
1472        $by .= $yy x $diff;
1473    } elsif ($diff < 0) {
1474        # if $xx eq "\x00", we can cut $by, otherwise we need to padd $bx
1475        $bx .= $xx x abs($diff);
1476    }
1477
1478    # or the strings together
1479    my $r = $bx | $by;
1480
1481    # and reverse the result again
1482    $bx = reverse $r;
1483
1484    # One of $bx or $by was negative, so need to flip bits in the result. In both
1485    # cases (one or two of them negative, or both positive) we need to get the
1486    # characters back.
1487    if ($sign eq '-') {
1488        $bx =~ tr<\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>
1489                 <0123456789abcdef>;
1490    } else {
1491        $bx =~ tr<\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>
1492                 <fedcba9876543210>;
1493    }
1494
1495    # leading zeros will be stripped by _from_hex()
1496    $bx = '0x' . $bx;
1497    $bx = $class -> _from_hex($bx);
1498
1499    $bx = $class -> _inc($bx) if $sign eq '-';
1500
1501    # avoid negative zero
1502    $sign = '+' if $class -> _is_zero($bx);
1503
1504    return $bx, $sign;
1505}
1506
1507sub _to_bin {
1508    # convert the number to a string of binary digits without prefix
1509    my ($class, $x) = @_;
1510    my $str    = '';
1511    my $tmp    = $class -> _copy($x);
1512    my $chunk = $class -> _new("16777216");     # 2^24 = 24 binary digits
1513    my $rem;
1514    until ($class -> _acmp($tmp, $chunk) < 0) {
1515        ($tmp, $rem) = $class -> _div($tmp, $chunk);
1516        $str = sprintf("%024b", $class -> _num($rem)) . $str;
1517    }
1518    unless ($class -> _is_zero($tmp)) {
1519        $str = sprintf("%b", $class -> _num($tmp)) . $str;
1520    }
1521    return length($str) ? $str : '0';
1522}
1523
1524sub _to_oct {
1525    # convert the number to a string of octal digits without prefix
1526    my ($class, $x) = @_;
1527    my $str    = '';
1528    my $tmp    = $class -> _copy($x);
1529    my $chunk = $class -> _new("16777216");     # 2^24 = 8 octal digits
1530    my $rem;
1531    until ($class -> _acmp($tmp, $chunk) < 0) {
1532        ($tmp, $rem) = $class -> _div($tmp, $chunk);
1533        $str = sprintf("%08o", $class -> _num($rem)) . $str;
1534    }
1535    unless ($class -> _is_zero($tmp)) {
1536        $str = sprintf("%o", $class -> _num($tmp)) . $str;
1537    }
1538    return length($str) ? $str : '0';
1539}
1540
1541sub _to_hex {
1542    # convert the number to a string of hexadecimal digits without prefix
1543    my ($class, $x) = @_;
1544    my $str    = '';
1545    my $tmp    = $class -> _copy($x);
1546    my $chunk = $class -> _new("16777216");     # 2^24 = 6 hexadecimal digits
1547    my $rem;
1548    until ($class -> _acmp($tmp, $chunk) < 0) {
1549        ($tmp, $rem) = $class -> _div($tmp, $chunk);
1550        $str = sprintf("%06x", $class -> _num($rem)) . $str;
1551    }
1552    unless ($class -> _is_zero($tmp)) {
1553        $str = sprintf("%x", $class -> _num($tmp)) . $str;
1554    }
1555    return length($str) ? $str : '0';
1556}
1557
1558sub _as_bin {
1559    # convert the number to a string of binary digits with prefix
1560    my ($class, $x) = @_;
1561    return '0b' . $class -> _to_bin($x);
1562}
1563
1564sub _as_oct {
1565    # convert the number to a string of octal digits with prefix
1566    my ($class, $x) = @_;
1567    return '0' . $class -> _to_oct($x);         # yes, 0 becomes "00"
1568}
1569
1570sub _as_hex {
1571    # convert the number to a string of hexadecimal digits with prefix
1572    my ($class, $x) = @_;
1573    return '0x' . $class -> _to_hex($x);
1574}
1575
1576sub _to_bytes {
1577    # convert the number to a string of bytes
1578    my ($class, $x) = @_;
1579    my $str    = '';
1580    my $tmp    = $class -> _copy($x);
1581    my $chunk = $class -> _new("65536");
1582    my $rem;
1583    until ($class -> _is_zero($tmp)) {
1584        ($tmp, $rem) = $class -> _div($tmp, $chunk);
1585        $str = pack('n', $class -> _num($rem)) . $str;
1586    }
1587    $str =~ s/^\0+//;
1588    return length($str) ? $str : "\x00";
1589}
1590
1591*_as_bytes = \&_to_bytes;
1592
1593sub _to_base {
1594    # convert the number to a string of digits in various bases
1595    my $class = shift;
1596    my $x     = shift;
1597    my $base  = shift;
1598    $base = $class -> _new($base) unless ref($base);
1599
1600    my $collseq;
1601    if (@_) {
1602        $collseq = shift;
1603        croak "The collation sequence must be a non-empty string"
1604          unless defined($collseq) && length($collseq);
1605    } else {
1606        if ($class -> _acmp($base, $class -> _new("94")) <= 0) {
1607            $collseq = '0123456789'                     #  48 ..  57
1608                     . 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'     #  65 ..  90
1609                     . 'abcdefghijklmnopqrstuvwxyz'     #  97 .. 122
1610                     . '!"#$%&\'()*+,-./'               #  33 ..  47
1611                     . ':;<=>?@'                        #  58 ..  64
1612                     . '[\\]^_`'                        #  91 ..  96
1613                     . '{|}~';                          # 123 .. 126
1614        } else {
1615            croak "When base > 94, a collation sequence must be given";
1616        }
1617    }
1618
1619    my @collseq = split '', $collseq;
1620
1621    my $str   = '';
1622    my $tmp   = $class -> _copy($x);
1623    my $rem;
1624    until ($class -> _is_zero($tmp)) {
1625        ($tmp, $rem) = $class -> _div($tmp, $base);
1626        my $num = $class -> _num($rem);
1627        croak "no character to represent '$num' in collation sequence",
1628          " (collation sequence is too short)" if $num > $#collseq;
1629        my $chr = $collseq[$num];
1630        $str = $chr . $str;
1631    }
1632    return $collseq[0] unless length $str;
1633    return $str;
1634}
1635
1636sub _to_base_num {
1637    # Convert the number to an array of integers in any base.
1638    my ($class, $x, $base) = @_;
1639
1640    # Make sure the base is an object and >= 2.
1641    $base = $class -> _new($base) unless ref($base);
1642    my $two = $class -> _two();
1643    croak "base must be >= 2" unless $class -> _acmp($base, $two) >= 0;
1644
1645    my $out   = [];
1646    my $xcopy = $class -> _copy($x);
1647    my $rem;
1648
1649    # Do all except the last (most significant) element.
1650    until ($class -> _acmp($xcopy, $base) < 0) {
1651        ($xcopy, $rem) = $class -> _div($xcopy, $base);
1652        unshift @$out, $rem;
1653    }
1654
1655    # Do the last (most significant element).
1656    unless ($class -> _is_zero($xcopy)) {
1657        unshift @$out, $xcopy;
1658    }
1659
1660    # $out is empty if $x is zero.
1661    unshift @$out, $class -> _zero() unless @$out;
1662
1663    return $out;
1664}
1665
1666sub _from_hex {
1667    # Convert a string of hexadecimal digits to a number.
1668
1669    my ($class, $hex) = @_;
1670    $hex =~ s/^0[xX]//;
1671
1672    # Find the largest number of hexadecimal digits that we can safely use with
1673    # 32 bit integers. There are 4 bits pr hexadecimal digit, and we use only
1674    # 31 bits to play safe. This gives us int(31 / 4) = 7.
1675
1676    my $len = length $hex;
1677    my $rem = 1 + ($len - 1) % 7;
1678
1679    # Do the first chunk.
1680
1681    my $ret = $class -> _new(int hex substr $hex, 0, $rem);
1682    return $ret if $rem == $len;
1683
1684    # Do the remaining chunks, if any.
1685
1686    my $shift = $class -> _new(1 << (4 * 7));
1687    for (my $offset = $rem ; $offset < $len ; $offset += 7) {
1688        my $part = int hex substr $hex, $offset, 7;
1689        $ret = $class -> _mul($ret, $shift);
1690        $ret = $class -> _add($ret, $class -> _new($part));
1691    }
1692
1693    return $ret;
1694}
1695
1696sub _from_oct {
1697    # Convert a string of octal digits to a number.
1698
1699    my ($class, $oct) = @_;
1700
1701    # Find the largest number of octal digits that we can safely use with 32
1702    # bit integers. There are 3 bits pr octal digit, and we use only 31 bits to
1703    # play safe. This gives us int(31 / 3) = 10.
1704
1705    my $len = length $oct;
1706    my $rem = 1 + ($len - 1) % 10;
1707
1708    # Do the first chunk.
1709
1710    my $ret = $class -> _new(int oct substr $oct, 0, $rem);
1711    return $ret if $rem == $len;
1712
1713    # Do the remaining chunks, if any.
1714
1715    my $shift = $class -> _new(1 << (3 * 10));
1716    for (my $offset = $rem ; $offset < $len ; $offset += 10) {
1717        my $part = int oct substr $oct, $offset, 10;
1718        $ret = $class -> _mul($ret, $shift);
1719        $ret = $class -> _add($ret, $class -> _new($part));
1720    }
1721
1722    return $ret;
1723}
1724
1725sub _from_bin {
1726    # Convert a string of binary digits to a number.
1727
1728    my ($class, $bin) = @_;
1729    $bin =~ s/^0[bB]//;
1730
1731    # The largest number of binary digits that we can safely use with 32 bit
1732    # integers is 31. We use only 31 bits to play safe.
1733
1734    my $len = length $bin;
1735    my $rem = 1 + ($len - 1) % 31;
1736
1737    # Do the first chunk.
1738
1739    my $ret = $class -> _new(int oct '0b' . substr $bin, 0, $rem);
1740    return $ret if $rem == $len;
1741
1742    # Do the remaining chunks, if any.
1743
1744    my $shift = $class -> _new(1 << 31);
1745    for (my $offset = $rem ; $offset < $len ; $offset += 31) {
1746        my $part = int oct '0b' . substr $bin, $offset, 31;
1747        $ret = $class -> _mul($ret, $shift);
1748        $ret = $class -> _add($ret, $class -> _new($part));
1749    }
1750
1751    return $ret;
1752}
1753
1754sub _from_bytes {
1755    # convert string of bytes to a number
1756    my ($class, $str) = @_;
1757    my $x    = $class -> _zero();
1758    my $base = $class -> _new("256");
1759    my $n    = length($str);
1760    for (my $i = 0 ; $i < $n ; ++$i) {
1761        $x = $class -> _mul($x, $base);
1762        my $byteval = $class -> _new(unpack 'C', substr($str, $i, 1));
1763        $x = $class -> _add($x, $byteval);
1764    }
1765    return $x;
1766}
1767
1768sub _from_base {
1769    # convert a string to a decimal number
1770    my $class = shift;
1771    my $str   = shift;
1772    my $base  = shift;
1773    $base = $class -> _new($base) unless ref($base);
1774
1775    my $n = length($str);
1776    my $x = $class -> _zero();
1777
1778    my $collseq;
1779    if (@_) {
1780        $collseq = shift();
1781    } else {
1782        if ($class -> _acmp($base, $class -> _new("36")) <= 0) {
1783            $str = uc $str;
1784            $collseq = '0123456789' . 'ABCDEFGHIJKLMNOPQRSTUVWXYZ';
1785        } elsif ($class -> _acmp($base, $class -> _new("94")) <= 0) {
1786            $collseq = '0123456789'                     #  48 ..  57
1787                     . 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'     #  65 ..  90
1788                     . 'abcdefghijklmnopqrstuvwxyz'     #  97 .. 122
1789                     . '!"#$%&\'()*+,-./'               #  33 ..  47
1790                     . ':;<=>?@'                        #  58 ..  64
1791                     . '[\\]^_`'                        #  91 ..  96
1792                     . '{|}~';                          # 123 .. 126
1793        } else {
1794            croak "When base > 94, a collation sequence must be given";
1795        }
1796        $collseq = substr $collseq, 0, $class -> _num($base);
1797    }
1798
1799    # Create a mapping from each character in the collation sequence to the
1800    # corresponding integer. Check for duplicates in the collation sequence.
1801
1802    my @collseq = split '', $collseq;
1803    my %collseq;
1804    for my $num (0 .. $#collseq) {
1805        my $chr = $collseq[$num];
1806        die "duplicate character '$chr' in collation sequence"
1807          if exists $collseq{$chr};
1808        $collseq{$chr} = $num;
1809    }
1810
1811    for (my $i = 0 ; $i < $n ; ++$i) {
1812        my $chr = substr($str, $i, 1);
1813        die "input character '$chr' does not exist in collation sequence"
1814          unless exists $collseq{$chr};
1815        $x = $class -> _mul($x, $base);
1816        my $num = $class -> _new($collseq{$chr});
1817        $x = $class -> _add($x, $num);
1818    }
1819
1820    return $x;
1821}
1822
1823sub _from_base_num {
1824    # Convert an array in the given base to a number.
1825    my ($class, $in, $base) = @_;
1826
1827    # Make sure the base is an object and >= 2.
1828    $base = $class -> _new($base) unless ref($base);
1829    my $two = $class -> _two();
1830    croak "base must be >= 2" unless $class -> _acmp($base, $two) >= 0;
1831
1832    # @$in = map { ref($_) ? $_ : $class -> _new($_) } @$in;
1833
1834    my $ele = $in -> [0];
1835
1836    $ele = $class -> _new($ele) unless ref($ele);
1837    my $x = $class -> _copy($ele);
1838
1839    for my $i (1 .. $#$in) {
1840        $x = $class -> _mul($x, $base);
1841        $ele = $in -> [$i];
1842        $ele = $class -> _new($ele) unless ref($ele);
1843        $x = $class -> _add($x, $ele);
1844    }
1845
1846    return $x;
1847}
1848
1849##############################################################################
1850# special modulus functions
1851
1852sub _modinv {
1853    # modular multiplicative inverse
1854    my ($class, $x, $y) = @_;
1855
1856    # modulo zero
1857    if ($class -> _is_zero($y)) {
1858        return;
1859    }
1860
1861    # modulo one
1862    if ($class -> _is_one($y)) {
1863        return ($class -> _zero(), '+');
1864    }
1865
1866    my $u = $class -> _zero();
1867    my $v = $class -> _one();
1868    my $a = $class -> _copy($y);
1869    my $b = $class -> _copy($x);
1870
1871    # Euclid's Algorithm for bgcd().
1872
1873    my $q;
1874    my $sign = 1;
1875    {
1876        ($a, $q, $b) = ($b, $class -> _div($a, $b));
1877        last if $class -> _is_zero($b);
1878
1879        my $vq = $class -> _mul($class -> _copy($v), $q);
1880        my $t = $class -> _add($vq, $u);
1881        $u = $v;
1882        $v = $t;
1883        $sign = -$sign;
1884        redo;
1885    }
1886
1887    # if the gcd is not 1, there exists no modular multiplicative inverse
1888    return unless $class -> _is_one($a);
1889
1890    ($v, $sign == 1 ? '+' : '-');
1891}
1892
1893sub _modpow {
1894    # modulus of power ($x ** $y) % $z
1895    my ($class, $num, $exp, $mod) = @_;
1896
1897    # a^b (mod 1) = 0 for all a and b
1898    if ($class -> _is_one($mod)) {
1899        return $class -> _zero();
1900    }
1901
1902    # 0^a (mod m) = 0 if m != 0, a != 0
1903    # 0^0 (mod m) = 1 if m != 0
1904    if ($class -> _is_zero($num)) {
1905        return $class -> _is_zero($exp) ? $class -> _one()
1906                                        : $class -> _zero();
1907    }
1908
1909    #  $num = $class -> _mod($num, $mod);   # this does not make it faster
1910
1911    my $acc = $class -> _copy($num);
1912    my $t   = $class -> _one();
1913
1914    my $expbin = $class -> _as_bin($exp);
1915    $expbin =~ s/^0b//;
1916    my $len = length($expbin);
1917
1918    while (--$len >= 0) {
1919        if (substr($expbin, $len, 1) eq '1') {
1920            $t = $class -> _mul($t, $acc);
1921            $t = $class -> _mod($t, $mod);
1922        }
1923        $acc = $class -> _mul($acc, $acc);
1924        $acc = $class -> _mod($acc, $mod);
1925    }
1926    return $t;
1927}
1928
1929sub _gcd {
1930    # Greatest common divisor.
1931
1932    my ($class, $x, $y) = @_;
1933
1934    # gcd(0, 0) = 0
1935    # gcd(0, a) = a, if a != 0
1936
1937    if ($class -> _acmp($x, $y) == 0) {
1938        return $class -> _copy($x);
1939    }
1940
1941    if ($class -> _is_zero($x)) {
1942        if ($class -> _is_zero($y)) {
1943            return $class -> _zero();
1944        } else {
1945            return $class -> _copy($y);
1946        }
1947    } else {
1948        if ($class -> _is_zero($y)) {
1949            return $class -> _copy($x);
1950        } else {
1951
1952            # Until $y is zero ...
1953
1954            $x = $class -> _copy($x);
1955            until ($class -> _is_zero($y)) {
1956
1957                # Compute remainder.
1958
1959                $x = $class -> _mod($x, $y);
1960
1961                # Swap $x and $y.
1962
1963                my $tmp = $x;
1964                $x = $class -> _copy($y);
1965                $y = $tmp;
1966            }
1967
1968            return $x;
1969        }
1970    }
1971}
1972
1973sub _lcm {
1974    # Least common multiple.
1975
1976    my ($class, $x, $y) = @_;
1977
1978    # lcm(0, x) = 0 for all x
1979
1980    return $class -> _zero()
1981      if ($class -> _is_zero($x) ||
1982          $class -> _is_zero($y));
1983
1984    my $gcd = $class -> _gcd($class -> _copy($x), $y);
1985    $x = $class -> _div($x, $gcd);
1986    $x = $class -> _mul($x, $y);
1987    return $x;
1988}
1989
1990sub _lucas {
1991    my ($class, $n) = @_;
1992
1993    $n = $class -> _num($n) if ref $n;
1994
1995    # In list context, use lucas(n) = lucas(n-1) + lucas(n-2)
1996
1997    if (wantarray) {
1998        my @y;
1999
2000        push @y, $class -> _two();
2001        return @y if $n == 0;
2002
2003        push @y, $class -> _one();
2004        return @y if $n == 1;
2005
2006        for (my $i = 2 ; $i <= $n ; ++ $i) {
2007            $y[$i] = $class -> _add($class -> _copy($y[$i - 1]), $y[$i - 2]);
2008        }
2009
2010        return @y;
2011    }
2012
2013    # In scalar context use that lucas(n) = fib(n-1) + fib(n+1).
2014    #
2015    # Remember that _fib() behaves differently in scalar context and list
2016    # context, so we must add scalar() to get the desired behaviour.
2017
2018    return $class -> _two() if $n == 0;
2019
2020    return $class -> _add(scalar($class -> _fib($n - 1)),
2021                          scalar($class -> _fib($n + 1)));
2022}
2023
2024sub _fib {
2025    my ($class, $n) = @_;
2026
2027    $n = $class -> _num($n) if ref $n;
2028
2029    # In list context, use fib(n) = fib(n-1) + fib(n-2)
2030
2031    if (wantarray) {
2032        my @y;
2033
2034        push @y, $class -> _zero();
2035        return @y if $n == 0;
2036
2037        push @y, $class -> _one();
2038        return @y if $n == 1;
2039
2040        for (my $i = 2 ; $i <= $n ; ++ $i) {
2041            $y[$i] = $class -> _add($class -> _copy($y[$i - 1]), $y[$i - 2]);
2042        }
2043
2044        return @y;
2045    }
2046
2047    # In scalar context use a fast algorithm that is much faster than the
2048    # recursive algorith used in list context.
2049
2050    my $cache = {};
2051    my $two = $class -> _two();
2052    my $fib;
2053
2054    $fib = sub {
2055        my $n = shift;
2056        return $class -> _zero() if $n <= 0;
2057        return $class -> _one()  if $n <= 2;
2058        return $cache -> {$n}    if exists $cache -> {$n};
2059
2060        my $k = int($n / 2);
2061        my $a = $fib -> ($k + 1);
2062        my $b = $fib -> ($k);
2063        my $y;
2064
2065        if ($n % 2 == 1) {
2066            # a*a + b*b
2067            $y = $class -> _add($class -> _mul($class -> _copy($a), $a),
2068                                $class -> _mul($class -> _copy($b), $b));
2069        } else {
2070            # (2*a - b)*b
2071            $y = $class -> _mul($class -> _sub($class -> _mul(
2072                   $class -> _copy($two), $a), $b), $b);
2073        }
2074
2075        $cache -> {$n} = $y;
2076        return $y;
2077    };
2078
2079    return $fib -> ($n);
2080}
2081
2082##############################################################################
2083##############################################################################
2084
20851;
2086
2087__END__
2088
2089=pod
2090
2091=head1 NAME
2092
2093Math::BigInt::Lib - virtual parent class for Math::BigInt libraries
2094
2095=head1 SYNOPSIS
2096
2097    # In the backend library for Math::BigInt et al.
2098
2099    package Math::BigInt::MyBackend;
2100
2101    use Math::BigInt::Lib;
2102    our @ISA = qw< Math::BigInt::Lib >;
2103
2104    sub _new { ... }
2105    sub _str { ... }
2106    sub _add { ... }
2107    str _sub { ... }
2108    ...
2109
2110    # In your main program.
2111
2112    use Math::BigInt lib => 'MyBackend';
2113
2114=head1 DESCRIPTION
2115
2116This module provides support for big integer calculations. It is not intended
2117to be used directly, but rather as a parent class for backend libraries used by
2118Math::BigInt, Math::BigFloat, Math::BigRat, and related modules.
2119
2120Other backend libraries include Math::BigInt::Calc, Math::BigInt::FastCalc,
2121Math::BigInt::GMP, and Math::BigInt::Pari.
2122
2123In order to allow for multiple big integer libraries, Math::BigInt was
2124rewritten to use a plug-in library for core math routines. Any module which
2125conforms to the API can be used by Math::BigInt by using this in your program:
2126
2127        use Math::BigInt lib => 'libname';
2128
2129'libname' is either the long name, like 'Math::BigInt::Pari', or only the short
2130version, like 'Pari'.
2131
2132=head2 General Notes
2133
2134A library only needs to deal with unsigned big integers. Testing of input
2135parameter validity is done by the caller, so there is no need to worry about
2136underflow (e.g., in C<_sub()> and C<_dec()>) or about division by zero (e.g.,
2137in C<_div()> and C<_mod()>)) or similar cases.
2138
2139Some libraries use methods that don't modify their argument, and some libraries
2140don't even use objects, but rather unblessed references. Because of this,
2141liberary methods are always called as class methods, not instance methods:
2142
2143    $x = Class -> method($x, $y);     # like this
2144    $x = $x -> method($y);            # not like this ...
2145    $x -> method($y);                 # ... or like this
2146
2147And with boolean methods
2148
2149    $bool = Class -> method($x, $y);  # like this
2150    $bool = $x -> method($y);         # not like this
2151
2152Return values are always objects, strings, Perl scalars, or true/false for
2153comparison routines.
2154
2155=head3 API version
2156
2157=over 4
2158
2159=item CLASS-E<gt>api_version()
2160
2161This method is no longer used and can be omitted. Methods that are not
2162implemented by a subclass will be inherited from this class.
2163
2164=back
2165
2166=head3 Constructors
2167
2168The following methods are mandatory: _new(), _str(), _add(), and _sub().
2169However, computations will be very slow without _mul() and _div().
2170
2171=over 4
2172
2173=item CLASS-E<gt>_new(STR)
2174
2175Convert a string representing an unsigned decimal number to an object
2176representing the same number. The input is normalized, i.e., it matches
2177C<^(0|[1-9]\d*)$>.
2178
2179=item CLASS-E<gt>_zero()
2180
2181Return an object representing the number zero.
2182
2183=item CLASS-E<gt>_one()
2184
2185Return an object representing the number one.
2186
2187=item CLASS-E<gt>_two()
2188
2189Return an object representing the number two.
2190
2191=item CLASS-E<gt>_ten()
2192
2193Return an object representing the number ten.
2194
2195=item CLASS-E<gt>_from_bin(STR)
2196
2197Return an object given a string representing a binary number. The input has a
2198'0b' prefix and matches the regular expression C<^0[bB](0|1[01]*)$>.
2199
2200=item CLASS-E<gt>_from_oct(STR)
2201
2202Return an object given a string representing an octal number. The input has a
2203'0' prefix and matches the regular expression C<^0[1-7]*$>.
2204
2205=item CLASS-E<gt>_from_hex(STR)
2206
2207Return an object given a string representing a hexadecimal number. The input
2208has a '0x' prefix and matches the regular expression
2209C<^0x(0|[1-9a-fA-F][\da-fA-F]*)$>.
2210
2211=item CLASS-E<gt>_from_bytes(STR)
2212
2213Returns an object given a byte string representing the number. The byte string
2214is in big endian byte order, so the two-byte input string "\x01\x00" should
2215give an output value representing the number 256.
2216
2217=item CLASS-E<gt>_from_base(STR, BASE, COLLSEQ)
2218
2219Returns an object given a string STR, a base BASE, and a collation sequence
2220COLLSEQ. Each character in STR represents a numerical value identical to the
2221character's position in COLLSEQ. All characters in STR must be present in
2222COLLSEQ.
2223
2224If BASE is less than or equal to 94, and a collation sequence is not specified,
2225the following default collation sequence is used. It contains of all the 94
2226printable ASCII characters except space/blank:
2227
2228    0123456789                  # ASCII  48 to  57
2229    ABCDEFGHIJKLMNOPQRSTUVWXYZ  # ASCII  65 to  90
2230    abcdefghijklmnopqrstuvwxyz  # ASCII  97 to 122
2231    !"#$%&'()*+,-./             # ASCII  33 to  47
2232    :;<=>?@                     # ASCII  58 to  64
2233    [\]^_`                      # ASCII  91 to  96
2234    {|}~                        # ASCII 123 to 126
2235
2236If the default collation sequence is used, and the BASE is less than or equal
2237to 36, the letter case in STR is ignored.
2238
2239For instance, with base 3 and collation sequence "-/|", the character "-"
2240represents 0, "/" represents 1, and "|" represents 2. So if STR is "/|-", the
2241output is 1 * 3**2 + 2 * 3**1 + 0 * 3**0 = 15.
2242
2243The following examples show standard binary, octal, decimal, and hexadecimal
2244conversion. All examples return 250.
2245
2246    $x = $class -> _from_base("11111010", 2)
2247    $x = $class -> _from_base("372", 8)
2248    $x = $class -> _from_base("250", 10)
2249    $x = $class -> _from_base("FA", 16)
2250
2251Some more examples, all returning 250:
2252
2253    $x = $class -> _from_base("100021", 3)
2254    $x = $class -> _from_base("3322", 4)
2255    $x = $class -> _from_base("2000", 5)
2256    $x = $class -> _from_base("caaa", 5, "abcde")
2257    $x = $class -> _from_base("42", 62)
2258    $x = $class -> _from_base("2!", 94)
2259
2260=item CLASS-E<gt>_from_base_num(ARRAY, BASE)
2261
2262Returns an object given an array of values and a base. This method is
2263equivalent to C<_from_base()>, but works on numbers in an array rather than
2264characters in a string. Unlike C<_from_base()>, all input values may be
2265arbitrarily large.
2266
2267    $x = $class -> _from_base_num([1, 1, 0, 1], 2)    # $x is 13
2268    $x = $class -> _from_base_num([3, 125, 39], 128)  # $x is 65191
2269
2270=back
2271
2272=head3 Mathematical functions
2273
2274=over 4
2275
2276=item CLASS-E<gt>_add(OBJ1, OBJ2)
2277
2278Addition. Returns the result of adding OBJ2 to OBJ1.
2279
2280=item CLASS-E<gt>_mul(OBJ1, OBJ2)
2281
2282Multiplication. Returns the result of multiplying OBJ2 and OBJ1.
2283
2284=item CLASS-E<gt>_div(OBJ1, OBJ2)
2285
2286Division. In scalar context, returns the quotient after dividing OBJ1 by OBJ2
2287and truncating the result to an integer. In list context, return the quotient
2288and the remainder.
2289
2290=item CLASS-E<gt>_sub(OBJ1, OBJ2, FLAG)
2291
2292=item CLASS-E<gt>_sub(OBJ1, OBJ2)
2293
2294Subtraction. Returns the result of subtracting OBJ2 by OBJ1. If C<flag> is false
2295or omitted, OBJ1 might be modified. If C<flag> is true, OBJ2 might be modified.
2296
2297=item CLASS-E<gt>_sadd(OBJ1, SIGN1, OBJ2, SIGN2)
2298
2299Signed addition. Returns the result of adding OBJ2 with sign SIGN2 to OBJ1 with
2300sign SIGN1.
2301
2302    ($obj3, $sign3) = $class -> _sadd($obj1, $sign1, $obj2, $sign2);
2303
2304=item CLASS-E<gt>_ssub(OBJ1, SIGN1, OBJ2, SIGN2)
2305
2306Signed subtraction. Returns the result of subtracting OBJ2 with sign SIGN2 to
2307OBJ1 with sign SIGN1.
2308
2309    ($obj3, $sign3) = $class -> _sadd($obj1, $sign1, $obj2, $sign2);
2310
2311=item CLASS-E<gt>_dec(OBJ)
2312
2313Returns the result after decrementing OBJ by one.
2314
2315=item CLASS-E<gt>_inc(OBJ)
2316
2317Returns the result after incrementing OBJ by one.
2318
2319=item CLASS-E<gt>_mod(OBJ1, OBJ2)
2320
2321Returns OBJ1 modulo OBJ2, i.e., the remainder after dividing OBJ1 by OBJ2.
2322
2323=item CLASS-E<gt>_sqrt(OBJ)
2324
2325Returns the square root of OBJ, truncated to an integer.
2326
2327=item CLASS-E<gt>_root(OBJ, N)
2328
2329Returns the Nth root of OBJ, truncated to an integer.
2330
2331=item CLASS-E<gt>_fac(OBJ)
2332
2333Returns the factorial of OBJ, i.e., the product of all positive integers up to
2334and including OBJ.
2335
2336=item CLASS-E<gt>_dfac(OBJ)
2337
2338Returns the double factorial of OBJ. If OBJ is an even integer, returns the
2339product of all positive, even integers up to and including OBJ, i.e.,
23402*4*6*...*OBJ. If OBJ is an odd integer, returns the product of all positive,
2341odd integers, i.e., 1*3*5*...*OBJ.
2342
2343=item CLASS-E<gt>_pow(OBJ1, OBJ2)
2344
2345Returns OBJ1 raised to the power of OBJ2. By convention, 0**0 = 1.
2346
2347=item CLASS-E<gt>_modinv(OBJ1, OBJ2)
2348
2349Returns the modular multiplicative inverse, i.e., return OBJ3 so that
2350
2351    (OBJ3 * OBJ1) % OBJ2 = 1 % OBJ2
2352
2353The result is returned as two arguments. If the modular multiplicative inverse
2354does not exist, both arguments are undefined. Otherwise, the arguments are a
2355number (object) and its sign ("+" or "-").
2356
2357The output value, with its sign, must either be a positive value in the range
23581,2,...,OBJ2-1 or the same value subtracted OBJ2. For instance, if the input
2359arguments are objects representing the numbers 7 and 5, the method must either
2360return an object representing the number 3 and a "+" sign, since (3*7) % 5 = 1
2361% 5, or an object representing the number 2 and a "-" sign, since (-2*7) % 5 = 1
2362% 5.
2363
2364=item CLASS-E<gt>_modpow(OBJ1, OBJ2, OBJ3)
2365
2366Returns the modular exponentiation, i.e., (OBJ1 ** OBJ2) % OBJ3.
2367
2368=item CLASS-E<gt>_rsft(OBJ, N, B)
2369
2370Returns the result after shifting OBJ N digits to thee right in base B. This is
2371equivalent to performing integer division by B**N and discarding the remainder,
2372except that it might be much faster.
2373
2374For instance, if the object $obj represents the hexadecimal number 0xabcde,
2375then C<_rsft($obj, 2, 16)> returns an object representing the number 0xabc. The
2376"remainer", 0xde, is discarded and not returned.
2377
2378=item CLASS-E<gt>_lsft(OBJ, N, B)
2379
2380Returns the result after shifting OBJ N digits to the left in base B. This is
2381equivalent to multiplying by B**N, except that it might be much faster.
2382
2383=item CLASS-E<gt>_log_int(OBJ, B)
2384
2385Returns the logarithm of OBJ to base BASE truncted to an integer. This method
2386has two output arguments, the OBJECT and a STATUS. The STATUS is Perl scalar;
2387it is 1 if OBJ is the exact result, 0 if the result was truncted to give OBJ,
2388and undef if it is unknown whether OBJ is the exact result.
2389
2390=item CLASS-E<gt>_ilog2(OBJ)
2391
2392Returns the base 2 logarithm of OBJ rounded downwards to the nearest integer,
2393i.e., C<int(log2(OBJ))>. In list context, this method returns two output
2394arguments, the OBJECT and a STATUS. The STATUS is Perl scalar; it is 1 if OBJ
2395is the exact result, 0 if the result was truncted to give OBJ, and undef if it
2396is unknown whether OBJ is the exact result.
2397
2398This method is equivalent to the more general method _log_int() when it is used
2399with base 2 argument, but _ilog2() method might be faster.
2400
2401=item CLASS-E<gt>_ilog10(OBJ)
2402
2403Returns the base 10 logarithm of OBJ rounded downwards to the nearest integer,
2404i.e., C<int(log2(OBJ))>. In list context, this method returns two output
2405arguments, the OBJECT and a STATUS. The STATUS is Perl scalar; it is 1 if OBJ
2406is the exact result, 0 if the result was truncted to give OBJ, and undef if it
2407is unknown whether OBJ is the exact result.
2408
2409This method is equivalent to the more general method _log_int() when it is used
2410with base 10 argument, but _ilog10() method might be faster.
2411
2412Also, the output from _ilog10() is always 1 smaller than the output from
2413_len().
2414
2415=item CLASS-E<gt>_clog2(OBJ)
2416
2417Returns the base 2 logarithm of OBJ rounded upwards to the nearest integer,
2418i.e., C<ceil(log2(OBJ))>. In list context, this method returns two output
2419arguments, the OBJECT and a STATUS. The STATUS is Perl scalar; it is 1 if OBJ
2420is the exact result, 0 if the result was truncted to give OBJ, and undef if it
2421is unknown whether OBJ is the exact result.
2422
2423=item CLASS-E<gt>_clog10(OBJ)
2424
2425Returns the base 10 logarithm of OBJ rounded upnwards to the nearest integer,
2426i.e., C<ceil(log2(OBJ))>. In list context, this method returns two output
2427arguments, the OBJECT and a STATUS. The STATUS is Perl scalar; it is 1 if OBJ
2428is the exact result, 0 if the result was truncted to give OBJ, and undef if it
2429is unknown whether OBJ is the exact result.
2430
2431=item CLASS-E<gt>_gcd(OBJ1, OBJ2)
2432
2433Returns the greatest common divisor of OBJ1 and OBJ2.
2434
2435=item CLASS-E<gt>_lcm(OBJ1, OBJ2)
2436
2437Return the least common multiple of OBJ1 and OBJ2.
2438
2439=item CLASS-E<gt>_fib(OBJ)
2440
2441In scalar context, returns the nth Fibonacci number: _fib(0) returns 0, _fib(1)
2442returns 1, _fib(2) returns 1, _fib(3) returns 2 etc. In list context, returns
2443the Fibonacci numbers from F(0) to F(n): 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, ...
2444
2445=item CLASS-E<gt>_lucas(OBJ)
2446
2447In scalar context, returns the nth Lucas number: _lucas(0) returns 2, _lucas(1)
2448returns 1, _lucas(2) returns 3, etc. In list context, returns the Lucas numbers
2449from L(0) to L(n): 2, 1, 3, 4, 7, 11, 18, 29,47, 76, ...
2450
2451=back
2452
2453=head3 Bitwise operators
2454
2455=over 4
2456
2457=item CLASS-E<gt>_and(OBJ1, OBJ2)
2458
2459Returns bitwise and.
2460
2461=item CLASS-E<gt>_or(OBJ1, OBJ2)
2462
2463Returns bitwise or.
2464
2465=item CLASS-E<gt>_xor(OBJ1, OBJ2)
2466
2467Returns bitwise exclusive or.
2468
2469=item CLASS-E<gt>_sand(OBJ1, OBJ2, SIGN1, SIGN2)
2470
2471Returns bitwise signed and.
2472
2473=item CLASS-E<gt>_sor(OBJ1, OBJ2, SIGN1, SIGN2)
2474
2475Returns bitwise signed or.
2476
2477=item CLASS-E<gt>_sxor(OBJ1, OBJ2, SIGN1, SIGN2)
2478
2479Returns bitwise signed exclusive or.
2480
2481=back
2482
2483=head3 Boolean operators
2484
2485=over 4
2486
2487=item CLASS-E<gt>_is_zero(OBJ)
2488
2489Returns a true value if OBJ is zero, and false value otherwise.
2490
2491=item CLASS-E<gt>_is_one(OBJ)
2492
2493Returns a true value if OBJ is one, and false value otherwise.
2494
2495=item CLASS-E<gt>_is_two(OBJ)
2496
2497Returns a true value if OBJ is two, and false value otherwise.
2498
2499=item CLASS-E<gt>_is_ten(OBJ)
2500
2501Returns a true value if OBJ is ten, and false value otherwise.
2502
2503=item CLASS-E<gt>_is_even(OBJ)
2504
2505Return a true value if OBJ is an even integer, and a false value otherwise.
2506
2507=item CLASS-E<gt>_is_odd(OBJ)
2508
2509Return a true value if OBJ is an even integer, and a false value otherwise.
2510
2511=item CLASS-E<gt>_acmp(OBJ1, OBJ2)
2512
2513Compare OBJ1 and OBJ2 and return -1, 0, or 1, if OBJ1 is numerically less than,
2514equal to, or larger than OBJ2, respectively.
2515
2516=back
2517
2518=head3 String conversion
2519
2520=over 4
2521
2522=item CLASS-E<gt>_str(OBJ)
2523
2524Returns a string representing OBJ in decimal notation. The returned string
2525should have no leading zeros, i.e., it should match C<^(0|[1-9]\d*)$>.
2526
2527=item CLASS-E<gt>_to_bin(OBJ)
2528
2529Returns the binary string representation of OBJ.
2530
2531=item CLASS-E<gt>_to_oct(OBJ)
2532
2533Returns the octal string representation of the number.
2534
2535=item CLASS-E<gt>_to_hex(OBJ)
2536
2537Returns the hexadecimal string representation of the number.
2538
2539=item CLASS-E<gt>_to_bytes(OBJ)
2540
2541Returns a byte string representation of OBJ. The byte string is in big endian
2542byte order, so if OBJ represents the number 256, the output should be the
2543two-byte string "\x01\x00".
2544
2545=item CLASS-E<gt>_to_base(OBJ, BASE, COLLSEQ)
2546
2547Returns a string representation of OBJ in base BASE with collation sequence
2548COLLSEQ.
2549
2550    $val = $class -> _new("210");
2551    $str = $class -> _to_base($val, 10, "xyz")  # $str is "zyx"
2552
2553    $val = $class -> _new("32");
2554    $str = $class -> _to_base($val, 2, "-|")  # $str is "|-----"
2555
2556See _from_base() for more information.
2557
2558=item CLASS-E<gt>_to_base_num(OBJ, BASE)
2559
2560Converts the given number to the given base. This method is equivalent to
2561C<_to_base()>, but returns numbers in an array rather than characters in a
2562string. In the output, the first element is the most significant. Unlike
2563C<_to_base()>, all input values may be arbitrarily large.
2564
2565    $x = $class -> _to_base_num(13, 2)        # $x is [1, 1, 0, 1]
2566    $x = $class -> _to_base_num(65191, 128)   # $x is [3, 125, 39]
2567
2568=item CLASS-E<gt>_as_bin(OBJ)
2569
2570Like C<_to_bin()> but with a '0b' prefix.
2571
2572=item CLASS-E<gt>_as_oct(OBJ)
2573
2574Like C<_to_oct()> but with a '0' prefix.
2575
2576=item CLASS-E<gt>_as_hex(OBJ)
2577
2578Like C<_to_hex()> but with a '0x' prefix.
2579
2580=item CLASS-E<gt>_as_bytes(OBJ)
2581
2582This is an alias to C<_to_bytes()>.
2583
2584=back
2585
2586=head3 Numeric conversion
2587
2588=over 4
2589
2590=item CLASS-E<gt>_num(OBJ)
2591
2592Returns a Perl scalar number representing the number OBJ as close as
2593possible. Since Perl scalars have limited precision, the returned value might
2594not be exactly the same as OBJ.
2595
2596=back
2597
2598=head3 Miscellaneous
2599
2600=over 4
2601
2602=item CLASS-E<gt>_copy(OBJ)
2603
2604Returns a true copy OBJ.
2605
2606=item CLASS-E<gt>_len(OBJ)
2607
2608Returns the number of the decimal digits in OBJ. The output is a Perl scalar.
2609
2610=item CLASS-E<gt>_zeros(OBJ)
2611
2612Returns the number of trailing decimal zeros. The output is a Perl scalar. The
2613number zero has no trailing decimal zeros.
2614
2615=item CLASS-E<gt>_digit(OBJ, N)
2616
2617Returns the Nth digit in OBJ as a Perl scalar. N is a Perl scalar, where zero
2618refers to the rightmost (least significant) digit, and negative values count
2619from the left (most significant digit). If $obj represents the number 123, then
2620
2621    CLASS->_digit($obj,  0)     # returns 3
2622    CLASS->_digit($obj,  1)     # returns 2
2623    CLASS->_digit($obj,  2)     # returns 1
2624    CLASS->_digit($obj, -1)     # returns 1
2625
2626=item CLASS-E<gt>_digitsum(OBJ)
2627
2628Returns the sum of the base 10 digits.
2629
2630=item CLASS-E<gt>_check(OBJ)
2631
2632Returns true if the object is invalid and false otherwise. Preferably, the true
2633value is a string describing the problem with the object. This is a check
2634routine to test the internal state of the object for corruption.
2635
2636=item CLASS-E<gt>_set(OBJ)
2637
2638xxx
2639
2640=back
2641
2642=head2 API version 2
2643
2644The following methods are required for an API version of 2 or greater.
2645
2646=head3 Constructors
2647
2648=over 4
2649
2650=item CLASS-E<gt>_1ex(N)
2651
2652Return an object representing the number 10**N where N E<gt>= 0 is a Perl
2653scalar.
2654
2655=back
2656
2657=head3 Mathematical functions
2658
2659=over 4
2660
2661=item CLASS-E<gt>_nok(OBJ1, OBJ2)
2662
2663Return the binomial coefficient OBJ1 over OBJ1.
2664
2665=back
2666
2667=head3 Miscellaneous
2668
2669=over 4
2670
2671=item CLASS-E<gt>_alen(OBJ)
2672
2673Return the approximate number of decimal digits of the object. The output is a
2674Perl scalar.
2675
2676=back
2677
2678=head1 WRAP YOUR OWN
2679
2680If you want to port your own favourite C library for big numbers to the
2681Math::BigInt interface, you can take any of the already existing modules as a
2682rough guideline. You should really wrap up the latest Math::BigInt and
2683Math::BigFloat testsuites with your module, and replace in them any of the
2684following:
2685
2686        use Math::BigInt;
2687
2688by this:
2689
2690        use Math::BigInt lib => 'yourlib';
2691
2692This way you ensure that your library really works 100% within Math::BigInt.
2693
2694=head1 BUGS
2695
2696Please report any bugs or feature requests to
2697C<bug-math-bigint at rt.cpan.org>, or through the web interface at
2698L<https://rt.cpan.org/Ticket/Create.html?Queue=Math-BigInt>
2699(requires login).
2700We will be notified, and then you'll automatically be notified of progress on
2701your bug as I make changes.
2702
2703=head1 SUPPORT
2704
2705You can find documentation for this module with the perldoc command.
2706
2707    perldoc Math::BigInt::Calc
2708
2709You can also look for information at:
2710
2711=over 4
2712
2713=item * GitHub Source Repository
2714
2715L<https://github.com/pjacklam/p5-Math-BigInt>
2716
2717=item * RT: CPAN's request tracker
2718
2719L<https://rt.cpan.org/Public/Dist/Display.html?Name=Math-BigInt>
2720
2721=item * MetaCPAN
2722
2723L<https://metacpan.org/release/Math-BigInt>
2724
2725=item * CPAN Testers Matrix
2726
2727L<http://matrix.cpantesters.org/?dist=Math-BigInt>
2728
2729=back
2730
2731=head1 LICENSE
2732
2733This program is free software; you may redistribute it and/or modify it under
2734the same terms as Perl itself.
2735
2736=head1 AUTHOR
2737
2738Peter John Acklam, E<lt>pjacklam@gmail.comE<gt>
2739
2740Code and documentation based on the Math::BigInt::Calc module by Tels
2741E<lt>nospam-abuse@bloodgate.comE<gt>
2742
2743=head1 SEE ALSO
2744
2745L<Math::BigInt>, L<Math::BigInt::Calc>, L<Math::BigInt::GMP>,
2746L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.
2747
2748=cut
2749