File Coverage

blib/lib/Math/BigInt/Lib.pm
Criterion Covered Total %
statement 219 956 22.9
branch 88 404 21.7
condition 11 42 26.1
subroutine 20 99 20.2
pod n/a
total 338 1501 22.5


line stmt bran cond sub pod time code
1             package Math::BigInt::Lib;
2              
3 51     51   1017 use 5.006001;
  51         218  
4 51     51   298 use strict;
  51         116  
  51         1006  
5 51     51   253 use warnings;
  51         100  
  51         3148  
6              
7             our $VERSION = '1.999842';
8             $VERSION =~ tr/_//d;
9              
10 51     51   362 use Carp;
  51         107  
  51         107831  
11              
12             use overload
13              
14             # overload key: with_assign
15              
16             '+' => sub {
17 0     0   0 my $class = ref $_[0];
18 0         0 my $x = $class -> _copy($_[0]);
19 0 0       0 my $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
20 0         0 return $class -> _add($x, $y);
21             },
22              
23             '-' => sub {
24 0     0   0 my $class = ref $_[0];
25 0         0 my ($x, $y);
26 0 0       0 if ($_[2]) { # if swapped
27 0         0 $y = $_[0];
28 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
29             } else {
30 0         0 $x = $class -> _copy($_[0]);
31 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
32             }
33 0         0 return $class -> _sub($x, $y);
34             },
35              
36             '*' => sub {
37 0     0   0 my $class = ref $_[0];
38 0         0 my $x = $class -> _copy($_[0]);
39 0 0       0 my $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
40 0         0 return $class -> _mul($x, $y);
41             },
42              
43             '/' => sub {
44 0     0   0 my $class = ref $_[0];
45 0         0 my ($x, $y);
46 0 0       0 if ($_[2]) { # if swapped
47 0         0 $y = $_[0];
48 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
49             } else {
50 0         0 $x = $class -> _copy($_[0]);
51 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
52             }
53 0         0 return $class -> _div($x, $y);
54             },
55              
56             '%' => sub {
57 0     0   0 my $class = ref $_[0];
58 0         0 my ($x, $y);
59 0 0       0 if ($_[2]) { # if swapped
60 0         0 $y = $_[0];
61 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
62             } else {
63 0         0 $x = $class -> _copy($_[0]);
64 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
65             }
66 0         0 return $class -> _mod($x, $y);
67             },
68              
69             '**' => sub {
70 0     0   0 my $class = ref $_[0];
71 0         0 my ($x, $y);
72 0 0       0 if ($_[2]) { # if swapped
73 0         0 $y = $_[0];
74 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
75             } else {
76 0         0 $x = $class -> _copy($_[0]);
77 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
78             }
79 0         0 return $class -> _pow($x, $y);
80             },
81              
82             '<<' => sub {
83 0     0   0 my $class = ref $_[0];
84 0         0 my ($x, $y);
85 0 0       0 if ($_[2]) { # if swapped
86 0         0 $y = $class -> _num($_[0]);
87 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
88             } else {
89 0         0 $x = $_[0];
90 0 0       0 $y = ref($_[1]) ? $class -> _num($_[1]) : $_[1];
91             }
92 0         0 return $class -> _lsft($x, $y);
93             },
94              
95             '>>' => sub {
96 0     0   0 my $class = ref $_[0];
97 0         0 my ($x, $y);
98 0 0       0 if ($_[2]) { # if swapped
99 0         0 $y = $_[0];
100 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
101             } else {
102 0         0 $x = $class -> _copy($_[0]);
103 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
104             }
105 0         0 return $class -> _rsft($x, $y);
106             },
107              
108             # overload key: num_comparison
109              
110             '<' => sub {
111 0     0   0 my $class = ref $_[0];
112 0         0 my ($x, $y);
113 0 0       0 if ($_[2]) { # if swapped
114 0         0 $y = $_[0];
115 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
116             } else {
117 0         0 $x = $class -> _copy($_[0]);
118 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
119             }
120 0         0 return $class -> _acmp($x, $y) < 0;
121             },
122              
123             '<=' => sub {
124 0     0   0 my $class = ref $_[0];
125 0         0 my ($x, $y);
126 0 0       0 if ($_[2]) { # if swapped
127 0         0 $y = $_[0];
128 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
129             } else {
130 0         0 $x = $class -> _copy($_[0]);
131 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
132             }
133 0         0 return $class -> _acmp($x, $y) <= 0;
134             },
135              
136             '>' => sub {
137 0     0   0 my $class = ref $_[0];
138 0         0 my ($x, $y);
139 0 0       0 if ($_[2]) { # if swapped
140 0         0 $y = $_[0];
141 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
142             } else {
143 0         0 $x = $class -> _copy($_[0]);
144 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
145             }
146 0         0 return $class -> _acmp($x, $y) > 0;
147             },
148              
149             '>=' => sub {
150 0     0   0 my $class = ref $_[0];
151 0         0 my ($x, $y);
152 0 0       0 if ($_[2]) { # if swapped
153 0         0 $y = $_[0];
154 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
155             } else {
156 0         0 $x = $class -> _copy($_[0]);
157 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
158             }
159 0         0 return $class -> _acmp($x, $y) >= 0;
160             },
161              
162             '==' => sub {
163 11605     11605   20842 my $class = ref $_[0];
164 11605         29896 my $x = $class -> _copy($_[0]);
165 11605 50       26408 my $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
166 11605         28389 return $class -> _acmp($x, $y) == 0;
167             },
168              
169             '!=' => sub {
170 0     0   0 my $class = ref $_[0];
171 0         0 my $x = $class -> _copy($_[0]);
172 0 0       0 my $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
173 0         0 return $class -> _acmp($x, $y) != 0;
174             },
175              
176             # overload key: 3way_comparison
177              
178             '<=>' => sub {
179 0     0   0 my $class = ref $_[0];
180 0         0 my ($x, $y);
181 0 0       0 if ($_[2]) { # if swapped
182 0         0 $y = $_[0];
183 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
184             } else {
185 0         0 $x = $class -> _copy($_[0]);
186 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
187             }
188 0         0 return $class -> _acmp($x, $y);
189             },
190              
191             # overload key: binary
192              
193             '&' => sub {
194 0     0   0 my $class = ref $_[0];
195 0         0 my ($x, $y);
196 0 0       0 if ($_[2]) { # if swapped
197 0         0 $y = $_[0];
198 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
199             } else {
200 0         0 $x = $class -> _copy($_[0]);
201 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
202             }
203 0         0 return $class -> _and($x, $y);
204             },
205              
206             '|' => sub {
207 0     0   0 my $class = ref $_[0];
208 0         0 my ($x, $y);
209 0 0       0 if ($_[2]) { # if swapped
210 0         0 $y = $_[0];
211 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
212             } else {
213 0         0 $x = $class -> _copy($_[0]);
214 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
215             }
216 0         0 return $class -> _or($x, $y);
217             },
218              
219             '^' => sub {
220 0     0   0 my $class = ref $_[0];
221 0         0 my ($x, $y);
222 0 0       0 if ($_[2]) { # if swapped
223 0         0 $y = $_[0];
224 0 0       0 $x = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
225             } else {
226 0         0 $x = $class -> _copy($_[0]);
227 0 0       0 $y = ref($_[1]) ? $_[1] : $class -> _new($_[1]);
228             }
229 0         0 return $class -> _xor($x, $y);
230             },
231              
232             # overload key: func
233              
234 0     0   0 'abs' => sub { $_[0] },
235              
236             'sqrt' => sub {
237 0     0   0 my $class = ref $_[0];
238 0         0 return $class -> _sqrt($class -> _copy($_[0]));
239             },
240              
241 0     0   0 'int' => sub { $_[0] },
242              
243             # overload key: conversion
244              
245 0 0   0   0 'bool' => sub { ref($_[0]) -> _is_zero($_[0]) ? '' : 1; },
246              
247 4     4   23 '""' => sub { ref($_[0]) -> _str($_[0]); },
248              
249 0     0   0 '0+' => sub { ref($_[0]) -> _num($_[0]); },
250              
251 0     0   0 '=' => sub { ref($_[0]) -> _copy($_[0]); },
252              
253 51     51   3076 ;
  51         2209  
  51         3015  
254              
255             sub _new {
256 0     0   0 croak "@{[(caller 0)[3]]} method not implemented";
  0         0  
257             }
258              
259             sub _zero {
260 0     0   0 my $class = shift;
261 0         0 return $class -> _new("0");
262             }
263              
264             sub _one {
265 0     0   0 my $class = shift;
266 0         0 return $class -> _new("1");
267             }
268              
269             sub _two {
270 0     0   0 my $class = shift;
271 0         0 return $class -> _new("2");
272              
273             }
274             sub _ten {
275 0     0   0 my $class = shift;
276 0         0 return $class -> _new("10");
277             }
278              
279             sub _1ex {
280 0     0   0 my ($class, $exp) = @_;
281 0 0       0 $exp = $class -> _num($exp) if ref($exp);
282 0         0 return $class -> _new("1" . ("0" x $exp));
283             }
284              
285             sub _copy {
286 0     0   0 my ($class, $x) = @_;
287 0         0 return $class -> _new($class -> _str($x));
288             }
289              
290             # catch and throw away
291       51     sub import { }
292              
293             ##############################################################################
294             # convert back to string and number
295              
296             sub _str {
297             # Convert number from internal base 1eN format to string format. Internal
298             # format is always normalized, i.e., no leading zeros.
299 0     0   0 croak "@{[(caller 0)[3]]} method not implemented";
  0         0  
300             }
301              
302             sub _num {
303 0     0   0 my ($class, $x) = @_;
304 0         0 0 + $class -> _str($x);
305             }
306              
307             ##############################################################################
308             # actual math code
309              
310             sub _add {
311 0     0   0 croak "@{[(caller 0)[3]]} method not implemented";
  0         0  
312             }
313              
314             sub _sub {
315 0     0   0 croak "@{[(caller 0)[3]]} method not implemented";
  0         0  
316             }
317              
318             sub _mul {
319 0     0   0 my ($class, $x, $y) = @_;
320 0         0 my $sum = $class -> _zero();
321 0         0 my $i = $class -> _zero();
322 0         0 while ($class -> _acmp($i, $y) < 0) {
323 0         0 $sum = $class -> _add($sum, $x);
324 0         0 $i = $class -> _inc($i);
325             }
326 0         0 return $sum;
327             }
328              
329             sub _div {
330 0     0   0 my ($class, $x, $y) = @_;
331              
332 0 0       0 croak "@{[(caller 0)[3]]} requires non-zero divisor"
  0         0  
333             if $class -> _is_zero($y);
334              
335 0         0 my $r = $class -> _copy($x);
336 0         0 my $q = $class -> _zero();
337 0         0 while ($class -> _acmp($r, $y) >= 0) {
338 0         0 $q = $class -> _inc($q);
339 0         0 $r = $class -> _sub($r, $y);
340             }
341              
342 0 0       0 return $q, $r if wantarray;
343 0         0 return $q;
344             }
345              
346             sub _inc {
347 0     0   0 my ($class, $x) = @_;
348 0         0 $class -> _add($x, $class -> _one());
349             }
350              
351             sub _dec {
352 0     0   0 my ($class, $x) = @_;
353 0         0 $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              
359             sub _sadd {
360 65124     65124   101432 my $class = shift;
361 65124         125456 my ($xa, $xs, $ya, $ys, $flag) = @_;
362 65124         96084 my ($za, $zs);
363              
364             # If the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
365              
366 65124 100       130393 if ($xs eq $ys) {
367 36480 50       63488 if ($flag) {
368 0         0 $za = $class -> _add($ya, $xa);
369             } else {
370 36480         93429 $za = $class -> _add($xa, $ya);
371             }
372 36480 100       91273 $zs = $class -> _is_zero($za) ? '+' : $xs;
373 36480         137082 return $za, $zs;
374             }
375              
376 28644         73827 my $acmp = $class -> _acmp($xa, $ya); # abs(x) = abs(y)
377              
378 28644 100       63023 if ($acmp == 0) { # x = -y or -x = y
379 6904         16911 $za = $class -> _zero();
380 6904         11338 $zs = '+';
381 6904         25102 return $za, $zs;
382             }
383              
384 21740 100       40839 if ($acmp > 0) { # abs(x) > abs(y)
385 15783         40082 $za = $class -> _sub($xa, $ya, $flag);
386 15783         27707 $zs = $xs;
387             } else { # abs(x) < abs(y)
388 5957         16225 $za = $class -> _sub($ya, $xa, !$flag);
389 5957         10277 $zs = $ys;
390             }
391 21740         74390 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              
397             sub _ssub {
398 32947     32947   52868 my $class = shift;
399 32947         68359 my ($xa, $xs, $ya, $ys, $flag) = @_;
400              
401             # Swap sign of second operand and let _sadd() do the job.
402 32947 100       68215 $ys = $ys eq '+' ? '-' : '+';
403 32947         74031 $class -> _sadd($xa, $xs, $ya, $ys, $flag);
404             }
405              
406             ##############################################################################
407             # testing
408              
409             sub _acmp {
410             # Compare two (absolute) values. Return -1, 0, or 1.
411 0     0   0 my ($class, $x, $y) = @_;
412 0         0 my $xstr = $class -> _str($x);
413 0         0 my $ystr = $class -> _str($y);
414              
415 0 0       0 length($xstr) <=> length($ystr) || $xstr cmp $ystr;
416             }
417              
418             sub _len {
419 0     0   0 my ($class, $x) = @_;
420 0         0 CORE::length($class -> _str($x));
421             }
422              
423             sub _alen {
424 27     27   223 my ($class, $x) = @_;
425 27         72 $class -> _len($x);
426             }
427              
428             sub _digit {
429 0     0   0 my ($class, $x, $n) = @_;
430 0         0 substr($class ->_str($x), -($n+1), 1);
431             }
432              
433             sub _digitsum {
434 0     0   0 my ($class, $x) = @_;
435              
436 0         0 my $len = $class -> _len($x);
437 0         0 my $sum = $class -> _zero();
438 0         0 for (my $i = 0 ; $i < $len ; ++$i) {
439 0         0 my $digit = $class -> _digit($x, $i);
440 0         0 $digit = $class -> _new($digit);
441 0         0 $sum = $class -> _add($sum, $digit);
442             }
443              
444 0         0 return $sum;
445             }
446              
447             sub _zeros {
448 0     0   0 my ($class, $x) = @_;
449 0         0 my $str = $class -> _str($x);
450 0 0       0 $str =~ /[^0](0*)\z/ ? CORE::length($1) : 0;
451             }
452              
453             ##############################################################################
454             # _is_* routines
455              
456             sub _is_zero {
457             # return true if arg is zero
458 0     0   0 my ($class, $x) = @_;
459 0         0 $class -> _str($x) == 0;
460             }
461              
462             sub _is_even {
463             # return true if arg is even
464 0     0   0 my ($class, $x) = @_;
465 0         0 substr($class -> _str($x), -1, 1) % 2 == 0;
466             }
467              
468             sub _is_odd {
469             # return true if arg is odd
470 0     0   0 my ($class, $x) = @_;
471 0         0 substr($class -> _str($x), -1, 1) % 2 != 0;
472             }
473              
474             sub _is_one {
475             # return true if arg is one
476 0     0   0 my ($class, $x) = @_;
477 0         0 $class -> _str($x) == 1;
478             }
479              
480             sub _is_two {
481             # return true if arg is two
482 0     0   0 my ($class, $x) = @_;
483 0         0 $class -> _str($x) == 2;
484             }
485              
486             sub _is_ten {
487             # return true if arg is ten
488 0     0   0 my ($class, $x) = @_;
489 0         0 $class -> _str($x) == 10;
490             }
491              
492             ###############################################################################
493             # check routine to test internal state for corruptions
494              
495             sub _check {
496             # used by the test suite
497 5498     5498   9873 my ($class, $x) = @_;
498 5498 50       13739 return "Input is undefined" unless defined $x;
499 5498 100       12168 return "$x is not a reference" unless ref($x);
500 5497         11810 return 0;
501             }
502              
503             ###############################################################################
504              
505             sub _mod {
506             # modulus
507 0     0   0 my ($class, $x, $y) = @_;
508              
509 0 0       0 croak "@{[(caller 0)[3]]} requires non-zero second operand"
  0         0  
510             if $class -> _is_zero($y);
511              
512 0 0       0 if ($class -> can('_div')) {
513 0         0 $x = $class -> _copy($x);
514 0         0 my ($q, $r) = $class -> _div($x, $y);
515 0         0 return $r;
516             } else {
517 0         0 my $r = $class -> _copy($x);
518 0         0 while ($class -> _acmp($r, $y) >= 0) {
519 0         0 $r = $class -> _sub($r, $y);
520             }
521 0         0 return $r;
522             }
523             }
524              
525             ##############################################################################
526             # shifts
527              
528             sub _rsft {
529 0     0   0 my ($class, $x, $n, $b) = @_;
530 0 0       0 $b = $class -> _new($b) unless ref $b;
531 0         0 return scalar $class -> _div($x, $class -> _pow($class -> _copy($b), $n));
532             }
533              
534             sub _lsft {
535 0     0   0 my ($class, $x, $n, $b) = @_;
536 0 0       0 $b = $class -> _new($b) unless ref $b;
537 0         0 return $class -> _mul($x, $class -> _pow($class -> _copy($b), $n));
538             }
539              
540             sub _pow {
541             # power of $x to $y
542 0     0   0 my ($class, $x, $y) = @_;
543              
544 0 0       0 if ($class -> _is_zero($y)) {
545 0         0 return $class -> _one(); # y == 0 => x => 1
546             }
547              
548 0 0 0     0 if (($class -> _is_one($x)) || # x == 1
549             ($class -> _is_one($y))) # or y == 1
550             {
551 0         0 return $x;
552             }
553              
554 0 0       0 if ($class -> _is_zero($x)) {
555 0         0 return $class -> _zero(); # 0 ** y => 0 (if not y <= 0)
556             }
557              
558 0         0 my $pow2 = $class -> _one();
559              
560 0         0 my $y_bin = $class -> _as_bin($y);
561 0         0 $y_bin =~ s/^0b//;
562 0         0 my $len = length($y_bin);
563              
564 0         0 while (--$len > 0) {
565 0 0       0 $pow2 = $class -> _mul($pow2, $x) if substr($y_bin, $len, 1) eq '1';
566 0         0 $x = $class -> _mul($x, $x);
567             }
568              
569 0         0 $x = $class -> _mul($x, $pow2);
570 0         0 return $x;
571             }
572              
573             sub _nok {
574             # Return binomial coefficient (n over k).
575 0     0   0 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 0         0 my $twok = $class -> _mul($class -> _two(), $class -> _copy($k));
  0         0  
582 0 0       0 if ($class -> _acmp($twok, $n) > 0) {
583 0         0 $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 0 0       0 if ($class -> _is_zero($k)) {
598 0         0 return $class -> _one();
599             }
600              
601             # Make a copy of the original n, in case the subclass modifies n in-place.
602              
603 0         0 my $n_orig = $class -> _copy($n);
604              
605             # n = 5, f = 6, d = 2 (cf. example above)
606              
607 0         0 $n = $class -> _sub($n, $k);
608 0         0 $n = $class -> _inc($n);
609              
610 0         0 my $f = $class -> _copy($n);
611 0         0 $f = $class -> _inc($f);
612              
613 0         0 my $d = $class -> _two();
614              
615             # while f <= n (the original n, that is) ...
616              
617 0         0 while ($class -> _acmp($f, $n_orig) <= 0) {
618 0         0 $n = $class -> _mul($n, $f);
619 0         0 $n = $class -> _div($n, $d);
620 0         0 $f = $class -> _inc($f);
621 0         0 $d = $class -> _inc($d);
622             }
623              
624 0         0 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              
646             sub _fac {
647             # factorial
648 0     0   0 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 0         0 my $p = $class -> _one();
654 0         0 my $r = $class -> _one();
655 0         0 my $two = $class -> _two();
656              
657 0         0 my ($log2n) = $class -> _log_int($class -> _copy($x), $two);
658 0         0 my $h = $class -> _zero();
659 0         0 my $shift = $class -> _zero();
660 0         0 my $k = $class -> _one();
661              
662 0         0 while ($class -> _acmp($h, $x)) {
663 0         0 $shift = $class -> _add($shift, $h);
664 0         0 $h = $class -> _rsft($class -> _copy($x), $log2n, $two);
665 0 0       0 $log2n = $class -> _dec($log2n) if !$class -> _is_zero($log2n);
666 0         0 my $high = $class -> _copy($h);
667 0 0       0 $high = $class -> _dec($high) if $class -> _is_even($h);
668 0         0 while ($class -> _acmp($k, $high)) {
669 0         0 $k = $class -> _add($k, $two);
670 0         0 $p = $class -> _mul($p, $k);
671             }
672 0         0 $r = $class -> _mul($r, $p);
673             }
674 0         0 return $class -> _lsft($r, $shift, $two);
675             }
676              
677             sub _dfac {
678             # double factorial
679 77     77   164 my ($class, $x) = @_;
680              
681 77         214 my $two = $class -> _two();
682              
683 77 50       202 if ($class -> _acmp($x, $two) < 0) {
684 0         0 return $class -> _one();
685             }
686              
687 77         227 my $i = $class -> _copy($x);
688 77         205 while ($class -> _acmp($i, $two) > 0) {
689 210         710 $i = $class -> _sub($i, $two);
690 210         531 $x = $class -> _mul($x, $i);
691             }
692              
693 77         256 return $x;
694             }
695              
696             sub _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 0     0   0 my ($class, $x, $base) = @_;
701              
702             # X == 0 => NaN
703 0 0       0 return if $class -> _is_zero($x);
704              
705 0 0       0 $base = $class -> _new(2) unless defined($base);
706 0 0       0 $base = $class -> _new($base) unless ref($base);
707              
708             # BASE 0 or 1 => NaN
709 0 0 0     0 return if $class -> _is_zero($base) || $class -> _is_one($base);
710              
711             # X == 1 => 0 (is exact)
712 0 0       0 if ($class -> _is_one($x)) {
713 0         0 return $class -> _zero(), 1;
714             }
715              
716 0         0 my $cmp = $class -> _acmp($x, $base);
717              
718             # X == BASE => 1 (is exact)
719 0 0       0 if ($cmp == 0) {
720 0         0 return $class -> _one(), 1;
721             }
722              
723             # 1 < X < BASE => 0 (is truncated)
724 0 0       0 if ($cmp < 0) {
725 0         0 return $class -> _zero(), 0;
726             }
727              
728 0         0 my $y;
729              
730             # log(x) / log(b) = log(xm * 10^xe) / log(bm * 10^be)
731             # = (log(xm) + xe*(log(10))) / (log(bm) + be*log(10))
732              
733             {
734 0         0 my $x_str = $class -> _str($x);
  0         0  
735 0         0 my $b_str = $class -> _str($base);
736 0         0 my $xm = "." . $x_str;
737 0         0 my $bm = "." . $b_str;
738 0         0 my $xe = length($x_str);
739 0         0 my $be = length($b_str);
740 0         0 my $log10 = log(10);
741 0         0 my $guess = int((log($xm) + $xe * $log10) / (log($bm) + $be * $log10));
742 0         0 $y = $class -> _new($guess);
743             }
744              
745 0         0 my $trial = $class -> _pow($class -> _copy($base), $y);
746 0         0 my $acmp = $class -> _acmp($trial, $x);
747              
748             # Did we get the exact result?
749              
750 0 0       0 return $y, 1 if $acmp == 0;
751              
752             # Too small?
753              
754 0         0 while ($acmp < 0) {
755 0         0 $trial = $class -> _mul($trial, $base);
756 0         0 $y = $class -> _inc($y);
757 0         0 $acmp = $class -> _acmp($trial, $x);
758             }
759              
760             # Too big?
761              
762 0         0 while ($acmp > 0) {
763 0         0 $trial = $class -> _div($trial, $base);
764 0         0 $y = $class -> _dec($y);
765 0         0 $acmp = $class -> _acmp($trial, $x);
766             }
767              
768 0 0       0 return $y, 1 if $acmp == 0; # result is exact
769 0         0 return $y, 0; # result is too small
770             }
771              
772             sub _sqrt {
773             # square-root of $y in place
774 0     0   0 my ($class, $y) = @_;
775              
776 0 0       0 return $y if $class -> _is_zero($y);
777              
778 0         0 my $y_str = $class -> _str($y);
779 0         0 my $y_len = length($y_str);
780              
781             # Compute the guess $x.
782              
783 0         0 my $xm;
784             my $xe;
785 0 0       0 if ($y_len % 2 == 0) {
786 0         0 $xm = sqrt("." . $y_str);
787 0         0 $xe = $y_len / 2;
788 0         0 $xm = sprintf "%.0f", int($xm * 1e15);
789 0         0 $xe -= 15;
790             } else {
791 0         0 $xm = sqrt(".0" . $y_str);
792 0         0 $xe = ($y_len + 1) / 2;
793 0         0 $xm = sprintf "%.0f", int($xm * 1e16);
794 0         0 $xe -= 16;
795             }
796              
797 0         0 my $x;
798 0 0       0 if ($xe < 0) {
799 0         0 $x = substr $xm, 0, length($xm) + $xe;
800             } else {
801 0         0 $x = $xm . ("0" x $xe);
802             }
803              
804 0         0 $x = $class -> _new($x);
805              
806             # Newton's method for computing square root of y
807             #
808             # x(i+1) = x(i) - f(x(i)) / f'(x(i))
809             # = x(i) - (x(i)^2 - y) / (2 * x(i)) # use if x(i)^2 > y
810             # = x(i) + (y - x(i)^2) / (2 * x(i)) # use if x(i)^2 < y
811              
812             # Determine if x, our guess, is too small, correct, or too large.
813              
814 0         0 my $xsq = $class -> _mul($class -> _copy($x), $x); # x(i)^2
815 0         0 my $acmp = $class -> _acmp($xsq, $y); # x(i)^2 <=> y
816              
817             # Only assign a value to this variable if we will be using it.
818              
819 0         0 my $two;
820 0 0       0 $two = $class -> _two() if $acmp != 0;
821              
822             # If x is too small, do one iteration of Newton's method. Since the
823             # function f(x) = x^2 - y is concave and monotonically increasing, the next
824             # guess for x will either be correct or too large.
825              
826 0 0       0 if ($acmp < 0) {
827              
828             # x(i+1) = x(i) + (y - x(i)^2) / (2 * x(i))
829              
830 0         0 my $numer = $class -> _sub($class -> _copy($y), $xsq); # y - x(i)^2
831 0         0 my $denom = $class -> _mul($class -> _copy($two), $x); # 2 * x(i)
832 0         0 my $delta = $class -> _div($numer, $denom);
833              
834 0 0       0 unless ($class -> _is_zero($delta)) {
835 0         0 $x = $class -> _add($x, $delta);
836 0         0 $xsq = $class -> _mul($class -> _copy($x), $x); # x(i)^2
837 0         0 $acmp = $class -> _acmp($xsq, $y); # x(i)^2 <=> y
838             }
839             }
840              
841             # If our guess for x is too large, apply Newton's method repeatedly until
842             # we either have got the correct value, or the delta is zero.
843              
844 0         0 while ($acmp > 0) {
845              
846             # x(i+1) = x(i) - (x(i)^2 - y) / (2 * x(i))
847              
848 0         0 my $numer = $class -> _sub($xsq, $y); # x(i)^2 - y
849 0         0 my $denom = $class -> _mul($class -> _copy($two), $x); # 2 * x(i)
850 0         0 my $delta = $class -> _div($numer, $denom);
851 0 0       0 last if $class -> _is_zero($delta);
852              
853 0         0 $x = $class -> _sub($x, $delta);
854 0         0 $xsq = $class -> _mul($class -> _copy($x), $x); # x(i)^2
855 0         0 $acmp = $class -> _acmp($xsq, $y); # x(i)^2 <=> y
856             }
857              
858             # When the delta is zero, our value for x might still be too large. We
859             # require that the outout is either exact or too small (i.e., rounded down
860             # to the nearest integer), so do a final check.
861              
862 0         0 while ($acmp > 0) {
863 0         0 $x = $class -> _dec($x);
864 0         0 $xsq = $class -> _mul($class -> _copy($x), $x); # x(i)^2
865 0         0 $acmp = $class -> _acmp($xsq, $y); # x(i)^2 <=> y
866             }
867              
868 0         0 return $x;
869             }
870              
871             sub _root {
872 0     0   0 my ($class, $y, $n) = @_;
873              
874 0 0 0     0 return $y if $class -> _is_zero($y) || $class -> _is_one($y) ||
      0        
875             $class -> _is_one($n);
876              
877             # If y <= n, the result is always (truncated to) 1.
878              
879 0 0       0 return $class -> _one() if $class -> _acmp($y, $n) <= 0;
880              
881             # Compute the initial guess x of y^(1/n). When n is large, Newton's method
882             # converges slowly if the "guess" (initial value) is poor, so we need a
883             # good guess. It the guess is too small, the next guess will be too large,
884             # and from then on all guesses are too large.
885              
886 0         0 my $DEBUG = 0;
887              
888             # Split y into mantissa and exponent in base 10, so that
889             #
890             # y = xm * 10^xe, where 0 < xm < 1 and xe is an integer
891              
892 0         0 my $y_str = $class -> _str($y);
893 0         0 my $ym = "." . $y_str;
894 0         0 my $ye = length($y_str);
895              
896             # From this compute the approximate base 10 logarithm of y
897             #
898             # log_10(y) = log_10(ym) + log_10(ye^10)
899             # = log(ym)/log(10) + ye
900              
901 0         0 my $log10y = log($ym) / log(10) + $ye;
902              
903             # And from this compute the approximate base 10 logarithm of x, where
904             # x = y^(1/n)
905             #
906             # log_10(x) = log_10(y)/n
907              
908 0         0 my $log10x = $log10y / $class -> _num($n);
909              
910             # From this compute xm and xe, the mantissa and exponent (in base 10) of x,
911             # where 1 < xm <= 10 and xe is an integer.
912              
913 0         0 my $xe = int $log10x;
914 0         0 my $xm = 10 ** ($log10x - $xe);
915              
916             # Scale the mantissa and exponent to increase the integer part of ym, which
917             # gives us better accuracy.
918              
919 0 0       0 if ($DEBUG) {
920 0         0 print "\n";
921 0         0 print "y_str = $y_str\n";
922 0         0 print "ym = $ym\n";
923 0         0 print "ye = $ye\n";
924 0         0 print "log10y = $log10y\n";
925 0         0 print "log10x = $log10x\n";
926 0         0 print "xm = $xm\n";
927 0         0 print "xe = $xe\n";
928             }
929              
930 0 0       0 my $d = $xe < 15 ? $xe : 15;
931 0         0 $xm *= 10 ** $d;
932 0         0 $xe -= $d;
933              
934 0 0       0 if ($DEBUG) {
935 0         0 print "\n";
936 0         0 print "xm = $xm\n";
937 0         0 print "xe = $xe\n";
938             }
939              
940             # If the mantissa is not an integer, round up to nearest integer, and then
941             # convert the number to a string. It is important to always round up due to
942             # how Newton's method behaves in this case. If the initial guess is too
943             # small, the next guess will be too large, after which every succeeding
944             # guess converges the correct value from above. Now, if the initial guess
945             # is too small and n is large, the next guess will be much too large and
946             # require a large number of iterations to get close to the solution.
947             # Because of this, we are likely to find the solution faster if we make
948             # sure the initial guess is not too small.
949              
950 0         0 my $xm_int = int($xm);
951 0 0       0 my $x_str = sprintf '%.0f', $xm > $xm_int ? $xm_int + 1 : $xm_int;
952 0         0 $x_str .= "0" x $xe;
953              
954 0         0 my $x = $class -> _new($x_str);
955              
956 0 0       0 if ($DEBUG) {
957 0         0 print "xm = $xm\n";
958 0         0 print "xe = $xe\n";
959 0         0 print "\n";
960 0         0 print "x_str = $x_str (initial guess)\n";
961 0         0 print "\n";
962             }
963              
964             # Use Newton's method for computing n'th root of y.
965             #
966             # x(i+1) = x(i) - f(x(i)) / f'(x(i))
967             # = x(i) - (x(i)^n - y) / (n * x(i)^(n-1)) # use if x(i)^n > y
968             # = x(i) + (y - x(i)^n) / (n * x(i)^(n-1)) # use if x(i)^n < y
969              
970             # Determine if x, our guess, is too small, correct, or too large. Rather
971             # than computing x(i)^n and x(i)^(n-1) directly, compute x(i)^(n-1) and
972             # then the same value multiplied by x.
973              
974 0         0 my $nm1 = $class -> _dec($class -> _copy($n)); # n-1
975 0         0 my $xpownm1 = $class -> _pow($class -> _copy($x), $nm1); # x(i)^(n-1)
976 0         0 my $xpown = $class -> _mul($class -> _copy($xpownm1), $x); # x(i)^n
977 0         0 my $acmp = $class -> _acmp($xpown, $y); # x(i)^n <=> y
978              
979 0 0       0 if ($DEBUG) {
980 0         0 print "\n";
981 0         0 print "x = ", $class -> _str($x), "\n";
982 0         0 print "x^n = ", $class -> _str($xpown), "\n";
983 0         0 print "y = ", $class -> _str($y), "\n";
984 0         0 print "acmp = $acmp\n";
985             }
986              
987             # If x is too small, do one iteration of Newton's method. Since the
988             # function f(x) = x^n - y is concave and monotonically increasing, the next
989             # guess for x will either be correct or too large.
990              
991 0 0       0 if ($acmp < 0) {
992              
993             # x(i+1) = x(i) + (y - x(i)^n) / (n * x(i)^(n-1))
994              
995 0         0 my $numer = $class -> _sub($class -> _copy($y), $xpown); # y - x(i)^n
996 0         0 my $denom = $class -> _mul($class -> _copy($n), $xpownm1); # n * x(i)^(n-1)
997 0         0 my $delta = $class -> _div($numer, $denom);
998              
999 0 0       0 if ($DEBUG) {
1000 0         0 print "\n";
1001 0         0 print "numer = ", $class -> _str($numer), "\n";
1002 0         0 print "denom = ", $class -> _str($denom), "\n";
1003 0         0 print "delta = ", $class -> _str($delta), "\n";
1004             }
1005              
1006 0 0       0 unless ($class -> _is_zero($delta)) {
1007 0         0 $x = $class -> _add($x, $delta);
1008 0         0 $xpownm1 = $class -> _pow($class -> _copy($x), $nm1); # x(i)^(n-1)
1009 0         0 $xpown = $class -> _mul($class -> _copy($xpownm1), $x); # x(i)^n
1010 0         0 $acmp = $class -> _acmp($xpown, $y); # x(i)^n <=> y
1011              
1012 0 0       0 if ($DEBUG) {
1013 0         0 print "\n";
1014 0         0 print "x = ", $class -> _str($x), "\n";
1015 0         0 print "x^n = ", $class -> _str($xpown), "\n";
1016 0         0 print "y = ", $class -> _str($y), "\n";
1017 0         0 print "acmp = $acmp\n";
1018             }
1019             }
1020             }
1021              
1022             # If our guess for x is too large, apply Newton's method repeatedly until
1023             # we either have got the correct value, or the delta is zero.
1024              
1025 0         0 while ($acmp > 0) {
1026              
1027             # x(i+1) = x(i) - (x(i)^n - y) / (n * x(i)^(n-1))
1028              
1029 0         0 my $numer = $class -> _sub($class -> _copy($xpown), $y); # x(i)^n - y
1030 0         0 my $denom = $class -> _mul($class -> _copy($n), $xpownm1); # n * x(i)^(n-1)
1031              
1032 0 0       0 if ($DEBUG) {
1033 0         0 print "numer = ", $class -> _str($numer), "\n";
1034 0         0 print "denom = ", $class -> _str($denom), "\n";
1035             }
1036              
1037 0         0 my $delta = $class -> _div($numer, $denom);
1038              
1039 0 0       0 if ($DEBUG) {
1040 0         0 print "delta = ", $class -> _str($delta), "\n";
1041             }
1042              
1043 0 0       0 last if $class -> _is_zero($delta);
1044              
1045 0         0 $x = $class -> _sub($x, $delta);
1046 0         0 $xpownm1 = $class -> _pow($class -> _copy($x), $nm1); # x(i)^(n-1)
1047 0         0 $xpown = $class -> _mul($class -> _copy($xpownm1), $x); # x(i)^n
1048 0         0 $acmp = $class -> _acmp($xpown, $y); # x(i)^n <=> y
1049              
1050 0 0       0 if ($DEBUG) {
1051 0         0 print "\n";
1052 0         0 print "x = ", $class -> _str($x), "\n";
1053 0         0 print "x^n = ", $class -> _str($xpown), "\n";
1054 0         0 print "y = ", $class -> _str($y), "\n";
1055 0         0 print "acmp = $acmp\n";
1056             }
1057             }
1058              
1059             # When the delta is zero, our value for x might still be too large. We
1060             # require that the outout is either exact or too small (i.e., rounded down
1061             # to the nearest integer), so do a final check.
1062              
1063 0         0 while ($acmp > 0) {
1064 0         0 $x = $class -> _dec($x);
1065 0         0 $xpown = $class -> _pow($class -> _copy($x), $n); # x(i)^n
1066 0         0 $acmp = $class -> _acmp($xpown, $y); # x(i)^n <=> y
1067             }
1068              
1069 0         0 return $x;
1070             }
1071              
1072             ##############################################################################
1073             # binary stuff
1074              
1075             sub _and {
1076 0     0   0 my ($class, $x, $y) = @_;
1077              
1078 0 0       0 return $x if $class -> _acmp($x, $y) == 0;
1079              
1080 0         0 my $m = $class -> _one();
1081 0         0 my $mask = $class -> _new("32768");
1082              
1083 0         0 my ($xr, $yr); # remainders after division
1084              
1085 0         0 my $xc = $class -> _copy($x);
1086 0         0 my $yc = $class -> _copy($y);
1087 0         0 my $z = $class -> _zero();
1088              
1089 0   0     0 until ($class -> _is_zero($xc) || $class -> _is_zero($yc)) {
1090 0         0 ($xc, $xr) = $class -> _div($xc, $mask);
1091 0         0 ($yc, $yr) = $class -> _div($yc, $mask);
1092 0         0 my $bits = $class -> _new($class -> _num($xr) & $class -> _num($yr));
1093 0         0 $z = $class -> _add($z, $class -> _mul($bits, $m));
1094 0         0 $m = $class -> _mul($m, $mask);
1095             }
1096              
1097 0         0 return $z;
1098             }
1099              
1100             sub _xor {
1101 0     0   0 my ($class, $x, $y) = @_;
1102              
1103 0 0       0 return $class -> _zero() if $class -> _acmp($x, $y) == 0;
1104              
1105 0         0 my $m = $class -> _one();
1106 0         0 my $mask = $class -> _new("32768");
1107              
1108 0         0 my ($xr, $yr); # remainders after division
1109              
1110 0         0 my $xc = $class -> _copy($x);
1111 0         0 my $yc = $class -> _copy($y);
1112 0         0 my $z = $class -> _zero();
1113              
1114 0   0     0 until ($class -> _is_zero($xc) || $class -> _is_zero($yc)) {
1115 0         0 ($xc, $xr) = $class -> _div($xc, $mask);
1116 0         0 ($yc, $yr) = $class -> _div($yc, $mask);
1117 0         0 my $bits = $class -> _new($class -> _num($xr) ^ $class -> _num($yr));
1118 0         0 $z = $class -> _add($z, $class -> _mul($bits, $m));
1119 0         0 $m = $class -> _mul($m, $mask);
1120             }
1121              
1122             # The loop above stops when the smallest of the two numbers is exhausted.
1123             # The remainder of the longer one will survive bit-by-bit, so we simple
1124             # multiply-add it in.
1125              
1126 0 0       0 $z = $class -> _add($z, $class -> _mul($xc, $m))
1127             unless $class -> _is_zero($xc);
1128 0 0       0 $z = $class -> _add($z, $class -> _mul($yc, $m))
1129             unless $class -> _is_zero($yc);
1130              
1131 0         0 return $z;
1132             }
1133              
1134             sub _or {
1135 0     0   0 my ($class, $x, $y) = @_;
1136              
1137 0 0       0 return $x if $class -> _acmp($x, $y) == 0; # shortcut (see _and)
1138              
1139 0         0 my $m = $class -> _one();
1140 0         0 my $mask = $class -> _new("32768");
1141              
1142 0         0 my ($xr, $yr); # remainders after division
1143              
1144 0         0 my $xc = $class -> _copy($x);
1145 0         0 my $yc = $class -> _copy($y);
1146 0         0 my $z = $class -> _zero();
1147              
1148 0   0     0 until ($class -> _is_zero($xc) || $class -> _is_zero($yc)) {
1149 0         0 ($xc, $xr) = $class -> _div($xc, $mask);
1150 0         0 ($yc, $yr) = $class -> _div($yc, $mask);
1151 0         0 my $bits = $class -> _new($class -> _num($xr) | $class -> _num($yr));
1152 0         0 $z = $class -> _add($z, $class -> _mul($bits, $m));
1153 0         0 $m = $class -> _mul($m, $mask);
1154             }
1155              
1156             # The loop above stops when the smallest of the two numbers is exhausted.
1157             # The remainder of the longer one will survive bit-by-bit, so we simple
1158             # multiply-add it in.
1159              
1160 0 0       0 $z = $class -> _add($z, $class -> _mul($xc, $m))
1161             unless $class -> _is_zero($xc);
1162 0 0       0 $z = $class -> _add($z, $class -> _mul($yc, $m))
1163             unless $class -> _is_zero($yc);
1164              
1165 0         0 return $z;
1166             }
1167              
1168             sub _sand {
1169 33     33   87 my ($class, $x, $sx, $y, $sy) = @_;
1170              
1171 33 50 33     99 return ($class -> _zero(), '+')
1172             if $class -> _is_zero($x) || $class -> _is_zero($y);
1173              
1174 33 100 100     149 my $sign = $sx eq '-' && $sy eq '-' ? '-' : '+';
1175              
1176 33         53 my ($bx, $by);
1177              
1178 33 100       87 if ($sx eq '-') { # if x is negative
1179             # two's complement: inc (dec unsigned value) and flip all "bits" in $bx
1180 24         78 $bx = $class -> _copy($x);
1181 24         78 $bx = $class -> _dec($bx);
1182 24         86 $bx = $class -> _as_hex($bx);
1183 24         129 $bx =~ s/^-?0x//;
1184 24         58 $bx =~ tr<0123456789abcdef>
1185             <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1186             } else { # if x is positive
1187 9         37 $bx = $class -> _as_hex($x); # get binary representation
1188 9         84 $bx =~ s/^-?0x//;
1189 9         26 $bx =~ tr
1190             <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1191             }
1192              
1193 33 100       112 if ($sy eq '-') { # if y is negative
1194             # two's complement: inc (dec unsigned value) and flip all "bits" in $by
1195 25         74 $by = $class -> _copy($y);
1196 25         76 $by = $class -> _dec($by);
1197 25         72 $by = $class -> _as_hex($by);
1198 25         103 $by =~ s/^-?0x//;
1199 25         58 $by =~ tr<0123456789abcdef>
1200             <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1201             } else {
1202 8         31 $by = $class -> _as_hex($y); # get binary representation
1203 8         47 $by =~ s/^-?0x//;
1204 8         17 $by =~ tr
1205             <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1206             }
1207              
1208             # now we have bit-strings from X and Y, reverse them for padding
1209 33         73 $bx = reverse $bx;
1210 33         62 $by = reverse $by;
1211              
1212             # padd the shorter string
1213 33 100       46 my $xx = "\x00"; $xx = "\x0f" if $sx eq '-';
  33         95  
1214 33 100       53 my $yy = "\x00"; $yy = "\x0f" if $sy eq '-';
  33         79  
1215 33         56 my $diff = CORE::length($bx) - CORE::length($by);
1216 33 100       135 if ($diff > 0) {
    50          
1217             # if $yy eq "\x00", we can cut $bx, otherwise we need to padd $by
1218 9         23 $by .= $yy x $diff;
1219             } elsif ($diff < 0) {
1220             # if $xx eq "\x00", we can cut $by, otherwise we need to padd $bx
1221 0         0 $bx .= $xx x abs($diff);
1222             }
1223              
1224             # and the strings together
1225 33         81 my $r = $bx & $by;
1226              
1227             # and reverse the result again
1228 33         58 $bx = reverse $r;
1229              
1230             # One of $bx or $by was negative, so need to flip bits in the result. In both
1231             # cases (one or two of them negative, or both positive) we need to get the
1232             # characters back.
1233 33 100       80 if ($sign eq '-') {
1234 16         34 $bx =~ tr<\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>
1235             <0123456789abcdef>;
1236             } else {
1237 17         33 $bx =~ tr<\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>
1238             ;
1239             }
1240              
1241             # leading zeros will be stripped by _from_hex()
1242 33         75 $bx = '0x' . $bx;
1243 33         118 $bx = $class -> _from_hex($bx);
1244              
1245 33 100       122 $bx = $class -> _inc($bx) if $sign eq '-';
1246              
1247             # avoid negative zero
1248 33 100       91 $sign = '+' if $class -> _is_zero($bx);
1249              
1250 33         150 return $bx, $sign;
1251             }
1252              
1253             sub _sxor {
1254 40     40   118 my ($class, $x, $sx, $y, $sy) = @_;
1255              
1256 40 50 33     124 return ($class -> _zero(), '+')
1257             if $class -> _is_zero($x) && $class -> _is_zero($y);
1258              
1259 40 100       102 my $sign = $sx ne $sy ? '-' : '+';
1260              
1261 40         62 my ($bx, $by);
1262              
1263 40 100       92 if ($sx eq '-') { # if x is negative
1264             # two's complement: inc (dec unsigned value) and flip all "bits" in $bx
1265 27         69 $bx = $class -> _copy($x);
1266 27         83 $bx = $class -> _dec($bx);
1267 27         83 $bx = $class -> _as_hex($bx);
1268 27         136 $bx =~ s/^-?0x//;
1269 27         60 $bx =~ tr<0123456789abcdef>
1270             <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1271             } else { # if x is positive
1272 13         70 $bx = $class -> _as_hex($x); # get binary representation
1273 13         60 $bx =~ s/^-?0x//;
1274 13         30 $bx =~ tr
1275             <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1276             }
1277              
1278 40 100       103 if ($sy eq '-') { # if y is negative
1279             # two's complement: inc (dec unsigned value) and flip all "bits" in $by
1280 29         88 $by = $class -> _copy($y);
1281 29         81 $by = $class -> _dec($by);
1282 29         83 $by = $class -> _as_hex($by);
1283 29         120 $by =~ s/^-?0x//;
1284 29         69 $by =~ tr<0123456789abcdef>
1285             <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1286             } else {
1287 11         43 $by = $class -> _as_hex($y); # get binary representation
1288 11         70 $by =~ s/^-?0x//;
1289 11         25 $by =~ tr
1290             <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1291             }
1292              
1293             # now we have bit-strings from X and Y, reverse them for padding
1294 40         80 $bx = reverse $bx;
1295 40         70 $by = reverse $by;
1296              
1297             # padd the shorter string
1298 40 100       54 my $xx = "\x00"; $xx = "\x0f" if $sx eq '-';
  40         100  
1299 40 100       65 my $yy = "\x00"; $yy = "\x0f" if $sy eq '-';
  40         86  
1300 40         75 my $diff = CORE::length($bx) - CORE::length($by);
1301 40 100       115 if ($diff > 0) {
    100          
1302             # if $yy eq "\x00", we can cut $bx, otherwise we need to padd $by
1303 9         32 $by .= $yy x $diff;
1304             } elsif ($diff < 0) {
1305             # if $xx eq "\x00", we can cut $by, otherwise we need to padd $bx
1306 3         13 $bx .= $xx x abs($diff);
1307             }
1308              
1309             # xor the strings together
1310 40         96 my $r = $bx ^ $by;
1311              
1312             # and reverse the result again
1313 40         74 $bx = reverse $r;
1314              
1315             # One of $bx or $by was negative, so need to flip bits in the result. In both
1316             # cases (one or two of them negative, or both positive) we need to get the
1317             # characters back.
1318 40 100       86 if ($sign eq '-') {
1319 24         55 $bx =~ tr<\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>
1320             <0123456789abcdef>;
1321             } else {
1322 16         37 $bx =~ tr<\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>
1323             ;
1324             }
1325              
1326             # leading zeros will be stripped by _from_hex()
1327 40         80 $bx = '0x' . $bx;
1328 40         127 $bx = $class -> _from_hex($bx);
1329              
1330 40 100       142 $bx = $class -> _inc($bx) if $sign eq '-';
1331              
1332             # avoid negative zero
1333 40 100       105 $sign = '+' if $class -> _is_zero($bx);
1334              
1335 40         187 return $bx, $sign;
1336             }
1337              
1338             sub _sor {
1339 35     35   99 my ($class, $x, $sx, $y, $sy) = @_;
1340              
1341 35 50 33     106 return ($class -> _zero(), '+')
1342             if $class -> _is_zero($x) && $class -> _is_zero($y);
1343              
1344 35 50 66     147 my $sign = $sx eq '-' || $sy eq '-' ? '-' : '+';
1345              
1346 35         58 my ($bx, $by);
1347              
1348 35 100       79 if ($sx eq '-') { # if x is negative
1349             # two's complement: inc (dec unsigned value) and flip all "bits" in $bx
1350 23         57 $bx = $class -> _copy($x);
1351 23         75 $bx = $class -> _dec($bx);
1352 23         76 $bx = $class -> _as_hex($bx);
1353 23         111 $bx =~ s/^-?0x//;
1354 23         55 $bx =~ tr<0123456789abcdef>
1355             <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1356             } else { # if x is positive
1357 12         56 $bx = $class -> _as_hex($x); # get binary representation
1358 12         57 $bx =~ s/^-?0x//;
1359 12         27 $bx =~ tr
1360             <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1361             }
1362              
1363 35 100       91 if ($sy eq '-') { # if y is negative
1364             # two's complement: inc (dec unsigned value) and flip all "bits" in $by
1365 24         68 $by = $class -> _copy($y);
1366 24         64 $by = $class -> _dec($by);
1367 24         66 $by = $class -> _as_hex($by);
1368 24         103 $by =~ s/^-?0x//;
1369 24         48 $by =~ tr<0123456789abcdef>
1370             <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1371             } else {
1372 11         37 $by = $class -> _as_hex($y); # get binary representation
1373 11         48 $by =~ s/^-?0x//;
1374 11         26 $by =~ tr
1375             <\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>;
1376             }
1377              
1378             # now we have bit-strings from X and Y, reverse them for padding
1379 35         80 $bx = reverse $bx;
1380 35         57 $by = reverse $by;
1381              
1382             # padd the shorter string
1383 35 100       54 my $xx = "\x00"; $xx = "\x0f" if $sx eq '-';
  35         96  
1384 35 100       55 my $yy = "\x00"; $yy = "\x0f" if $sy eq '-';
  35         75  
1385 35         62 my $diff = CORE::length($bx) - CORE::length($by);
1386 35 100       116 if ($diff > 0) {
    100          
1387             # if $yy eq "\x00", we can cut $bx, otherwise we need to padd $by
1388 12         33 $by .= $yy x $diff;
1389             } elsif ($diff < 0) {
1390             # if $xx eq "\x00", we can cut $by, otherwise we need to padd $bx
1391 3         22 $bx .= $xx x abs($diff);
1392             }
1393              
1394             # or the strings together
1395 35         90 my $r = $bx | $by;
1396              
1397             # and reverse the result again
1398 35         57 $bx = reverse $r;
1399              
1400             # One of $bx or $by was negative, so need to flip bits in the result. In both
1401             # cases (one or two of them negative, or both positive) we need to get the
1402             # characters back.
1403 35 50       96 if ($sign eq '-') {
1404 35         69 $bx =~ tr<\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>
1405             <0123456789abcdef>;
1406             } else {
1407 0         0 $bx =~ tr<\x0f\x0e\x0d\x0c\x0b\x0a\x09\x08\x07\x06\x05\x04\x03\x02\x01\x00>
1408             ;
1409             }
1410              
1411             # leading zeros will be stripped by _from_hex()
1412 35         74 $bx = '0x' . $bx;
1413 35         120 $bx = $class -> _from_hex($bx);
1414              
1415 35 50       187 $bx = $class -> _inc($bx) if $sign eq '-';
1416              
1417             # avoid negative zero
1418 35 50       145 $sign = '+' if $class -> _is_zero($bx);
1419              
1420 35         162 return $bx, $sign;
1421             }
1422              
1423             sub _to_bin {
1424             # convert the number to a string of binary digits without prefix
1425 101     101   216 my ($class, $x) = @_;
1426 101         1165 my $str = '';
1427 101         308 my $tmp = $class -> _copy($x);
1428 101         295 my $chunk = $class -> _new("16777216"); # 2^24 = 24 binary digits
1429 101         190 my $rem;
1430 101         288 until ($class -> _acmp($tmp, $chunk) < 0) {
1431 34         127 ($tmp, $rem) = $class -> _div($tmp, $chunk);
1432 34         139 $str = sprintf("%024b", $class -> _num($rem)) . $str;
1433             }
1434 101 100       293 unless ($class -> _is_zero($tmp)) {
1435 93         271 $str = sprintf("%b", $class -> _num($tmp)) . $str;
1436             }
1437 101 100       465 return length($str) ? $str : '0';
1438             }
1439              
1440             sub _to_oct {
1441             # convert the number to a string of octal digits without prefix
1442 48     48   113 my ($class, $x) = @_;
1443 48         97 my $str = '';
1444 48         135 my $tmp = $class -> _copy($x);
1445 48         150 my $chunk = $class -> _new("16777216"); # 2^24 = 8 octal digits
1446 48         99 my $rem;
1447 48         156 until ($class -> _acmp($tmp, $chunk) < 0) {
1448 24         120 ($tmp, $rem) = $class -> _div($tmp, $chunk);
1449 24         125 $str = sprintf("%08o", $class -> _num($rem)) . $str;
1450             }
1451 48 100       163 unless ($class -> _is_zero($tmp)) {
1452 40         157 $str = sprintf("%o", $class -> _num($tmp)) . $str;
1453             }
1454 48 100       230 return length($str) ? $str : '0';
1455             }
1456              
1457             sub _to_hex {
1458             # convert the number to a string of hexadecimal digits without prefix
1459 40     40   142 my ($class, $x) = @_;
1460 40         80 my $str = '';
1461 40         118 my $tmp = $class -> _copy($x);
1462 40         123 my $chunk = $class -> _new("16777216"); # 2^24 = 6 hexadecimal digits
1463 40         70 my $rem;
1464 40         132 until ($class -> _acmp($tmp, $chunk) < 0) {
1465 16         75 ($tmp, $rem) = $class -> _div($tmp, $chunk);
1466 16         66 $str = sprintf("%06x", $class -> _num($rem)) . $str;
1467             }
1468 40 100       125 unless ($class -> _is_zero($tmp)) {
1469 32         110 $str = sprintf("%x", $class -> _num($tmp)) . $str;
1470             }
1471 40 100       238 return length($str) ? $str : '0';
1472             }
1473              
1474             sub _as_bin {
1475             # convert the number to a string of binary digits with prefix
1476 0     0   0 my ($class, $x) = @_;
1477 0         0 return '0b' . $class -> _to_bin($x);
1478             }
1479              
1480             sub _as_oct {
1481             # convert the number to a string of octal digits with prefix
1482 0     0   0 my ($class, $x) = @_;
1483 0         0 return '0' . $class -> _to_oct($x); # yes, 0 becomes "00"
1484             }
1485              
1486             sub _as_hex {
1487             # convert the number to a string of hexadecimal digits with prefix
1488 0     0   0 my ($class, $x) = @_;
1489 0         0 return '0x' . $class -> _to_hex($x);
1490             }
1491              
1492             sub _to_bytes {
1493             # convert the number to a string of bytes
1494 0     0   0 my ($class, $x) = @_;
1495 0         0 my $str = '';
1496 0         0 my $tmp = $class -> _copy($x);
1497 0         0 my $chunk = $class -> _new("65536");
1498 0         0 my $rem;
1499 0         0 until ($class -> _is_zero($tmp)) {
1500 0         0 ($tmp, $rem) = $class -> _div($tmp, $chunk);
1501 0         0 $str = pack('n', $class -> _num($rem)) . $str;
1502             }
1503 0         0 $str =~ s/^\0+//;
1504 0 0       0 return length($str) ? $str : "\x00";
1505             }
1506              
1507             *_as_bytes = \&_to_bytes;
1508              
1509             sub _to_base {
1510             # convert the number to a string of digits in various bases
1511 0     0   0 my $class = shift;
1512 0         0 my $x = shift;
1513 0         0 my $base = shift;
1514 0 0       0 $base = $class -> _new($base) unless ref($base);
1515              
1516 0         0 my $collseq;
1517 0 0       0 if (@_) {
1518 0         0 $collseq = shift;
1519 0 0 0     0 croak "The collation sequence must be a non-empty string"
1520             unless defined($collseq) && length($collseq);
1521             } else {
1522 0 0       0 if ($class -> _acmp($base, $class -> _new("94")) <= 0) {
1523 0         0 $collseq = '0123456789' # 48 .. 57
1524             . 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' # 65 .. 90
1525             . 'abcdefghijklmnopqrstuvwxyz' # 97 .. 122
1526             . '!"#$%&\'()*+,-./' # 33 .. 47
1527             . ':;<=>?@' # 58 .. 64
1528             . '[\\]^_`' # 91 .. 96
1529             . '{|}~'; # 123 .. 126
1530             } else {
1531 0         0 croak "When base > 94, a collation sequence must be given";
1532             }
1533             }
1534              
1535 0         0 my @collseq = split '', $collseq;
1536              
1537 0         0 my $str = '';
1538 0         0 my $tmp = $class -> _copy($x);
1539 0         0 my $rem;
1540 0         0 until ($class -> _is_zero($tmp)) {
1541 0         0 ($tmp, $rem) = $class -> _div($tmp, $base);
1542 0         0 my $num = $class -> _num($rem);
1543 0 0       0 croak "no character to represent '$num' in collation sequence",
1544             " (collation sequence is too short)" if $num > $#collseq;
1545 0         0 my $chr = $collseq[$num];
1546 0         0 $str = $chr . $str;
1547             }
1548 0 0       0 return $collseq[0] unless length $str;
1549 0         0 return $str;
1550             }
1551              
1552             sub _to_base_num {
1553             # Convert the number to an array of integers in any base.
1554 0     0   0 my ($class, $x, $base) = @_;
1555              
1556             # Make sure the base is an object and >= 2.
1557 0 0       0 $base = $class -> _new($base) unless ref($base);
1558 0         0 my $two = $class -> _two();
1559 0 0       0 croak "base must be >= 2" unless $class -> _acmp($base, $two) >= 0;
1560              
1561 0         0 my $out = [];
1562 0         0 my $xcopy = $class -> _copy($x);
1563 0         0 my $rem;
1564              
1565             # Do all except the last (most significant) element.
1566 0         0 until ($class -> _acmp($xcopy, $base) < 0) {
1567 0         0 ($xcopy, $rem) = $class -> _div($xcopy, $base);
1568 0         0 unshift @$out, $rem;
1569             }
1570              
1571             # Do the last (most significant element).
1572 0 0       0 unless ($class -> _is_zero($xcopy)) {
1573 0         0 unshift @$out, $xcopy;
1574             }
1575              
1576             # $out is empty if $x is zero.
1577 0 0       0 unshift @$out, $class -> _zero() unless @$out;
1578              
1579 0         0 return $out;
1580             }
1581              
1582             sub _from_hex {
1583             # Convert a string of hexadecimal digits to a number.
1584              
1585 0     0   0 my ($class, $hex) = @_;
1586 0         0 $hex =~ s/^0[xX]//;
1587              
1588             # Find the largest number of hexadecimal digits that we can safely use with
1589             # 32 bit integers. There are 4 bits pr hexadecimal digit, and we use only
1590             # 31 bits to play safe. This gives us int(31 / 4) = 7.
1591              
1592 0         0 my $len = length $hex;
1593 0         0 my $rem = 1 + ($len - 1) % 7;
1594              
1595             # Do the first chunk.
1596              
1597 0         0 my $ret = $class -> _new(int hex substr $hex, 0, $rem);
1598 0 0       0 return $ret if $rem == $len;
1599              
1600             # Do the remaining chunks, if any.
1601              
1602 0         0 my $shift = $class -> _new(1 << (4 * 7));
1603 0         0 for (my $offset = $rem ; $offset < $len ; $offset += 7) {
1604 0         0 my $part = int hex substr $hex, $offset, 7;
1605 0         0 $ret = $class -> _mul($ret, $shift);
1606 0         0 $ret = $class -> _add($ret, $class -> _new($part));
1607             }
1608              
1609 0         0 return $ret;
1610             }
1611              
1612             sub _from_oct {
1613             # Convert a string of octal digits to a number.
1614              
1615 0     0   0 my ($class, $oct) = @_;
1616              
1617             # Find the largest number of octal digits that we can safely use with 32
1618             # bit integers. There are 3 bits pr octal digit, and we use only 31 bits to
1619             # play safe. This gives us int(31 / 3) = 10.
1620              
1621 0         0 my $len = length $oct;
1622 0         0 my $rem = 1 + ($len - 1) % 10;
1623              
1624             # Do the first chunk.
1625              
1626 0         0 my $ret = $class -> _new(int oct substr $oct, 0, $rem);
1627 0 0       0 return $ret if $rem == $len;
1628              
1629             # Do the remaining chunks, if any.
1630              
1631 0         0 my $shift = $class -> _new(1 << (3 * 10));
1632 0         0 for (my $offset = $rem ; $offset < $len ; $offset += 10) {
1633 0         0 my $part = int oct substr $oct, $offset, 10;
1634 0         0 $ret = $class -> _mul($ret, $shift);
1635 0         0 $ret = $class -> _add($ret, $class -> _new($part));
1636             }
1637              
1638 0         0 return $ret;
1639             }
1640              
1641             sub _from_bin {
1642             # Convert a string of binary digits to a number.
1643              
1644 0     0   0 my ($class, $bin) = @_;
1645 0         0 $bin =~ s/^0[bB]//;
1646              
1647             # The largest number of binary digits that we can safely use with 32 bit
1648             # integers is 31. We use only 31 bits to play safe.
1649              
1650 0         0 my $len = length $bin;
1651 0         0 my $rem = 1 + ($len - 1) % 31;
1652              
1653             # Do the first chunk.
1654              
1655 0         0 my $ret = $class -> _new(int oct '0b' . substr $bin, 0, $rem);
1656 0 0       0 return $ret if $rem == $len;
1657              
1658             # Do the remaining chunks, if any.
1659              
1660 0         0 my $shift = $class -> _new(1 << 31);
1661 0         0 for (my $offset = $rem ; $offset < $len ; $offset += 31) {
1662 0         0 my $part = int oct '0b' . substr $bin, $offset, 31;
1663 0         0 $ret = $class -> _mul($ret, $shift);
1664 0         0 $ret = $class -> _add($ret, $class -> _new($part));
1665             }
1666              
1667 0         0 return $ret;
1668             }
1669              
1670             sub _from_bytes {
1671             # convert string of bytes to a number
1672 0     0   0 my ($class, $str) = @_;
1673 0         0 my $x = $class -> _zero();
1674 0         0 my $base = $class -> _new("256");
1675 0         0 my $n = length($str);
1676 0         0 for (my $i = 0 ; $i < $n ; ++$i) {
1677 0         0 $x = $class -> _mul($x, $base);
1678 0         0 my $byteval = $class -> _new(unpack 'C', substr($str, $i, 1));
1679 0         0 $x = $class -> _add($x, $byteval);
1680             }
1681 0         0 return $x;
1682             }
1683              
1684             sub _from_base {
1685             # convert a string to a decimal number
1686 0     0   0 my $class = shift;
1687 0         0 my $str = shift;
1688 0         0 my $base = shift;
1689 0 0       0 $base = $class -> _new($base) unless ref($base);
1690              
1691 0         0 my $n = length($str);
1692 0         0 my $x = $class -> _zero();
1693              
1694 0         0 my $collseq;
1695 0 0       0 if (@_) {
1696 0         0 $collseq = shift();
1697             } else {
1698 0 0       0 if ($class -> _acmp($base, $class -> _new("36")) <= 0) {
    0          
1699 0         0 $str = uc $str;
1700 0         0 $collseq = '0123456789' . 'ABCDEFGHIJKLMNOPQRSTUVWXYZ';
1701             } elsif ($class -> _acmp($base, $class -> _new("94")) <= 0) {
1702 0         0 $collseq = '0123456789' # 48 .. 57
1703             . 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' # 65 .. 90
1704             . 'abcdefghijklmnopqrstuvwxyz' # 97 .. 122
1705             . '!"#$%&\'()*+,-./' # 33 .. 47
1706             . ':;<=>?@' # 58 .. 64
1707             . '[\\]^_`' # 91 .. 96
1708             . '{|}~'; # 123 .. 126
1709             } else {
1710 0         0 croak "When base > 94, a collation sequence must be given";
1711             }
1712 0         0 $collseq = substr $collseq, 0, $class -> _num($base);
1713             }
1714              
1715             # Create a mapping from each character in the collation sequence to the
1716             # corresponding integer. Check for duplicates in the collation sequence.
1717              
1718 0         0 my @collseq = split '', $collseq;
1719 0         0 my %collseq;
1720 0         0 for my $num (0 .. $#collseq) {
1721 0         0 my $chr = $collseq[$num];
1722             die "duplicate character '$chr' in collation sequence"
1723 0 0       0 if exists $collseq{$chr};
1724 0         0 $collseq{$chr} = $num;
1725             }
1726              
1727 0         0 for (my $i = 0 ; $i < $n ; ++$i) {
1728 0         0 my $chr = substr($str, $i, 1);
1729             die "input character '$chr' does not exist in collation sequence"
1730 0 0       0 unless exists $collseq{$chr};
1731 0         0 $x = $class -> _mul($x, $base);
1732 0         0 my $num = $class -> _new($collseq{$chr});
1733 0         0 $x = $class -> _add($x, $num);
1734             }
1735              
1736 0         0 return $x;
1737             }
1738              
1739             sub _from_base_num {
1740             # Convert an array in the given base to a number.
1741 0     0   0 my ($class, $in, $base) = @_;
1742              
1743             # Make sure the base is an object and >= 2.
1744 0 0       0 $base = $class -> _new($base) unless ref($base);
1745 0         0 my $two = $class -> _two();
1746 0 0       0 croak "base must be >= 2" unless $class -> _acmp($base, $two) >= 0;
1747              
1748             # @$in = map { ref($_) ? $_ : $class -> _new($_) } @$in;
1749              
1750 0         0 my $ele = $in -> [0];
1751              
1752 0 0       0 $ele = $class -> _new($ele) unless ref($ele);
1753 0         0 my $x = $class -> _copy($ele);
1754              
1755 0         0 for my $i (1 .. $#$in) {
1756 0         0 $x = $class -> _mul($x, $base);
1757 0         0 $ele = $in -> [$i];
1758 0 0       0 $ele = $class -> _new($ele) unless ref($ele);
1759 0         0 $x = $class -> _add($x, $ele);
1760             }
1761              
1762 0         0 return $x;
1763             }
1764              
1765             ##############################################################################
1766             # special modulus functions
1767              
1768             sub _modinv {
1769             # modular multiplicative inverse
1770 0     0   0 my ($class, $x, $y) = @_;
1771              
1772             # modulo zero
1773 0 0       0 if ($class -> _is_zero($y)) {
1774 0         0 return;
1775             }
1776              
1777             # modulo one
1778 0 0       0 if ($class -> _is_one($y)) {
1779 0         0 return ($class -> _zero(), '+');
1780             }
1781              
1782 0         0 my $u = $class -> _zero();
1783 0         0 my $v = $class -> _one();
1784 0         0 my $a = $class -> _copy($y);
1785 0         0 my $b = $class -> _copy($x);
1786              
1787             # Euclid's Algorithm for bgcd().
1788              
1789 0         0 my $q;
1790 0         0 my $sign = 1;
1791             {
1792 0         0 ($a, $q, $b) = ($b, $class -> _div($a, $b));
  0         0  
1793 0 0       0 last if $class -> _is_zero($b);
1794              
1795 0         0 my $vq = $class -> _mul($class -> _copy($v), $q);
1796 0         0 my $t = $class -> _add($vq, $u);
1797 0         0 $u = $v;
1798 0         0 $v = $t;
1799 0         0 $sign = -$sign;
1800 0         0 redo;
1801             }
1802              
1803             # if the gcd is not 1, there exists no modular multiplicative inverse
1804 0 0       0 return unless $class -> _is_one($a);
1805              
1806 0 0       0 ($v, $sign == 1 ? '+' : '-');
1807             }
1808              
1809             sub _modpow {
1810             # modulus of power ($x ** $y) % $z
1811 0     0   0 my ($class, $num, $exp, $mod) = @_;
1812              
1813             # a^b (mod 1) = 0 for all a and b
1814 0 0       0 if ($class -> _is_one($mod)) {
1815 0         0 return $class -> _zero();
1816             }
1817              
1818             # 0^a (mod m) = 0 if m != 0, a != 0
1819             # 0^0 (mod m) = 1 if m != 0
1820 0 0       0 if ($class -> _is_zero($num)) {
1821 0 0       0 return $class -> _is_zero($exp) ? $class -> _one()
1822             : $class -> _zero();
1823             }
1824              
1825             # $num = $class -> _mod($num, $mod); # this does not make it faster
1826              
1827 0         0 my $acc = $class -> _copy($num);
1828 0         0 my $t = $class -> _one();
1829              
1830 0         0 my $expbin = $class -> _as_bin($exp);
1831 0         0 $expbin =~ s/^0b//;
1832 0         0 my $len = length($expbin);
1833              
1834 0         0 while (--$len >= 0) {
1835 0 0       0 if (substr($expbin, $len, 1) eq '1') {
1836 0         0 $t = $class -> _mul($t, $acc);
1837 0         0 $t = $class -> _mod($t, $mod);
1838             }
1839 0         0 $acc = $class -> _mul($acc, $acc);
1840 0         0 $acc = $class -> _mod($acc, $mod);
1841             }
1842 0         0 return $t;
1843             }
1844              
1845             sub _gcd {
1846             # Greatest common divisor.
1847              
1848 0     0   0 my ($class, $x, $y) = @_;
1849              
1850             # gcd(0, 0) = 0
1851             # gcd(0, a) = a, if a != 0
1852              
1853 0 0       0 if ($class -> _acmp($x, $y) == 0) {
1854 0         0 return $class -> _copy($x);
1855             }
1856              
1857 0 0       0 if ($class -> _is_zero($x)) {
1858 0 0       0 if ($class -> _is_zero($y)) {
1859 0         0 return $class -> _zero();
1860             } else {
1861 0         0 return $class -> _copy($y);
1862             }
1863             } else {
1864 0 0       0 if ($class -> _is_zero($y)) {
1865 0         0 return $class -> _copy($x);
1866             } else {
1867              
1868             # Until $y is zero ...
1869              
1870 0         0 $x = $class -> _copy($x);
1871 0         0 until ($class -> _is_zero($y)) {
1872              
1873             # Compute remainder.
1874              
1875 0         0 $x = $class -> _mod($x, $y);
1876              
1877             # Swap $x and $y.
1878              
1879 0         0 my $tmp = $x;
1880 0         0 $x = $class -> _copy($y);
1881 0         0 $y = $tmp;
1882             }
1883              
1884 0         0 return $x;
1885             }
1886             }
1887             }
1888              
1889             sub _lcm {
1890             # Least common multiple.
1891              
1892 26     26   68 my ($class, $x, $y) = @_;
1893              
1894             # lcm(0, x) = 0 for all x
1895              
1896 26 100 100     83 return $class -> _zero()
1897             if ($class -> _is_zero($x) ||
1898             $class -> _is_zero($y));
1899              
1900 14         61 my $gcd = $class -> _gcd($class -> _copy($x), $y);
1901 14         55 $x = $class -> _div($x, $gcd);
1902 14         51 $x = $class -> _mul($x, $y);
1903 14         65 return $x;
1904             }
1905              
1906             sub _lucas {
1907 0     0     my ($class, $n) = @_;
1908              
1909 0 0         $n = $class -> _num($n) if ref $n;
1910              
1911             # In list context, use lucas(n) = lucas(n-1) + lucas(n-2)
1912              
1913 0 0         if (wantarray) {
1914 0           my @y;
1915              
1916 0           push @y, $class -> _two();
1917 0 0         return @y if $n == 0;
1918              
1919 0           push @y, $class -> _one();
1920 0 0         return @y if $n == 1;
1921              
1922 0           for (my $i = 2 ; $i <= $n ; ++ $i) {
1923 0           $y[$i] = $class -> _add($class -> _copy($y[$i - 1]), $y[$i - 2]);
1924             }
1925              
1926 0           return @y;
1927             }
1928              
1929             # In scalar context use that lucas(n) = fib(n-1) + fib(n+1).
1930             #
1931             # Remember that _fib() behaves differently in scalar context and list
1932             # context, so we must add scalar() to get the desired behaviour.
1933              
1934 0 0         return $class -> _two() if $n == 0;
1935              
1936 0           return $class -> _add(scalar($class -> _fib($n - 1)),
1937             scalar($class -> _fib($n + 1)));
1938             }
1939              
1940             sub _fib {
1941 0     0     my ($class, $n) = @_;
1942              
1943 0 0         $n = $class -> _num($n) if ref $n;
1944              
1945             # In list context, use fib(n) = fib(n-1) + fib(n-2)
1946              
1947 0 0         if (wantarray) {
1948 0           my @y;
1949              
1950 0           push @y, $class -> _zero();
1951 0 0         return @y if $n == 0;
1952              
1953 0           push @y, $class -> _one();
1954 0 0         return @y if $n == 1;
1955              
1956 0           for (my $i = 2 ; $i <= $n ; ++ $i) {
1957 0           $y[$i] = $class -> _add($class -> _copy($y[$i - 1]), $y[$i - 2]);
1958             }
1959              
1960 0           return @y;
1961             }
1962              
1963             # In scalar context use a fast algorithm that is much faster than the
1964             # recursive algorith used in list context.
1965              
1966 0           my $cache = {};
1967 0           my $two = $class -> _two();
1968 0           my $fib;
1969              
1970             $fib = sub {
1971 0     0     my $n = shift;
1972 0 0         return $class -> _zero() if $n <= 0;
1973 0 0         return $class -> _one() if $n <= 2;
1974 0 0         return $cache -> {$n} if exists $cache -> {$n};
1975              
1976 0           my $k = int($n / 2);
1977 0           my $a = $fib -> ($k + 1);
1978 0           my $b = $fib -> ($k);
1979 0           my $y;
1980              
1981 0 0         if ($n % 2 == 1) {
1982             # a*a + b*b
1983 0           $y = $class -> _add($class -> _mul($class -> _copy($a), $a),
1984             $class -> _mul($class -> _copy($b), $b));
1985             } else {
1986             # (2*a - b)*b
1987 0           $y = $class -> _mul($class -> _sub($class -> _mul(
1988             $class -> _copy($two), $a), $b), $b);
1989             }
1990              
1991 0           $cache -> {$n} = $y;
1992 0           return $y;
1993 0           };
1994              
1995 0           return $fib -> ($n);
1996             }
1997              
1998             ##############################################################################
1999             ##############################################################################
2000              
2001             1;
2002              
2003             __END__