← Index
NYTProf Performance Profile   « line view »
For svc/members/upsert
  Run on Tue Jan 13 11:50:22 2015
Reported on Tue Jan 13 12:09:49 2015

Filename/usr/share/perl/5.20/Math/BigInt/Calc.pm
StatementsExecuted 307 statements in 9.05ms
Subroutines
Calls P F Exclusive
Time		 Inclusive
Time		 Subroutine
111127µs203µsMath::BigInt::Calc::::BEGIN@117Math::BigInt::Calc::BEGIN@117
102150µs50µsMath::BigInt::Calc::::CORE:regcompMath::BigInt::Calc::CORE:regcomp (opcode)
11122µs22µsMath::BigInt::Calc::::BEGIN@3Math::BigInt::Calc::BEGIN@3
44212µs12µsMath::BigInt::Calc::::_newMath::BigInt::Calc::_new
11111µs55µsMath::BigInt::Calc::::BEGIN@1909Math::BigInt::Calc::BEGIN@1909
11111µs14µsMath::BigInt::Calc::::BEGIN@137Math::BigInt::Calc::BEGIN@137
11111µs13µsMath::BigInt::Calc::::BEGIN@475Math::BigInt::Calc::BEGIN@475
11111µs11µsMath::BigInt::Calc::::_base_lenMath::BigInt::Calc::_base_len
11110µs11µsMath::BigInt::Calc::::BEGIN@2186Math::BigInt::Calc::BEGIN@2186
1118µs10µsMath::BigInt::Calc::::BEGIN@2115Math::BigInt::Calc::BEGIN@2115
1118µs10µsMath::BigInt::Calc::::BEGIN@792Math::BigInt::Calc::BEGIN@792
1118µs10µsMath::BigInt::Calc::::BEGIN@165Math::BigInt::Calc::BEGIN@165
1118µs9µsMath::BigInt::Calc::::BEGIN@154Math::BigInt::Calc::BEGIN@154
1117µs8µsMath::BigInt::Calc::::BEGIN@2150Math::BigInt::Calc::BEGIN@2150
12416µs6µsMath::BigInt::Calc::::CORE:matchMath::BigInt::Calc::CORE:match (opcode)
1116µs18µsMath::BigInt::Calc::::BEGIN@4Math::BigInt::Calc::BEGIN@4
1116µs6µsMath::BigInt::Calc::::_strMath::BigInt::Calc::_str
1112µs2µsMath::BigInt::Calc::::_zeroMath::BigInt::Calc::_zero
1111µs1µsMath::BigInt::Calc::::importMath::BigInt::Calc::import
111400ns400nsMath::BigInt::Calc::::api_versionMath::BigInt::Calc::api_version (xsub)
0000s0sMath::BigInt::Calc::::_1exMath::BigInt::Calc::_1ex
0000s0sMath::BigInt::Calc::::__strip_zerosMath::BigInt::Calc::__strip_zeros
0000s0sMath::BigInt::Calc::::_acmpMath::BigInt::Calc::_acmp
0000s0sMath::BigInt::Calc::::_addMath::BigInt::Calc::_add
0000s0sMath::BigInt::Calc::::_andMath::BigInt::Calc::_and
0000s0sMath::BigInt::Calc::::_as_binMath::BigInt::Calc::_as_bin
0000s0sMath::BigInt::Calc::::_as_hexMath::BigInt::Calc::_as_hex
0000s0sMath::BigInt::Calc::::_as_octMath::BigInt::Calc::_as_oct
0000s0sMath::BigInt::Calc::::_checkMath::BigInt::Calc::_check
0000s0sMath::BigInt::Calc::::_copyMath::BigInt::Calc::_copy
0000s0sMath::BigInt::Calc::::_decMath::BigInt::Calc::_dec
0000s0sMath::BigInt::Calc::::_digitMath::BigInt::Calc::_digit
0000s0sMath::BigInt::Calc::::_div_use_divMath::BigInt::Calc::_div_use_div
0000s0sMath::BigInt::Calc::::_div_use_div_64Math::BigInt::Calc::_div_use_div_64
0000s0sMath::BigInt::Calc::::_div_use_mulMath::BigInt::Calc::_div_use_mul
0000s0sMath::BigInt::Calc::::_facMath::BigInt::Calc::_fac
0000s0sMath::BigInt::Calc::::_from_binMath::BigInt::Calc::_from_bin
0000s0sMath::BigInt::Calc::::_from_hexMath::BigInt::Calc::_from_hex
0000s0sMath::BigInt::Calc::::_from_octMath::BigInt::Calc::_from_oct
0000s0sMath::BigInt::Calc::::_gcdMath::BigInt::Calc::_gcd
0000s0sMath::BigInt::Calc::::_incMath::BigInt::Calc::_inc
0000s0sMath::BigInt::Calc::::_is_evenMath::BigInt::Calc::_is_even
0000s0sMath::BigInt::Calc::::_is_oddMath::BigInt::Calc::_is_odd
0000s0sMath::BigInt::Calc::::_is_oneMath::BigInt::Calc::_is_one
0000s0sMath::BigInt::Calc::::_is_tenMath::BigInt::Calc::_is_ten
0000s0sMath::BigInt::Calc::::_is_twoMath::BigInt::Calc::_is_two
0000s0sMath::BigInt::Calc::::_is_zeroMath::BigInt::Calc::_is_zero
0000s0sMath::BigInt::Calc::::_lenMath::BigInt::Calc::_len
0000s0sMath::BigInt::Calc::::_log_intMath::BigInt::Calc::_log_int
0000s0sMath::BigInt::Calc::::_lsftMath::BigInt::Calc::_lsft
0000s0sMath::BigInt::Calc::::_modMath::BigInt::Calc::_mod
0000s0sMath::BigInt::Calc::::_modinvMath::BigInt::Calc::_modinv
0000s0sMath::BigInt::Calc::::_modpowMath::BigInt::Calc::_modpow
0000s0sMath::BigInt::Calc::::_mul_use_divMath::BigInt::Calc::_mul_use_div
0000s0sMath::BigInt::Calc::::_mul_use_div_64Math::BigInt::Calc::_mul_use_div_64
0000s0sMath::BigInt::Calc::::_mul_use_mulMath::BigInt::Calc::_mul_use_mul
0000s0sMath::BigInt::Calc::::_nokMath::BigInt::Calc::_nok
0000s0sMath::BigInt::Calc::::_numMath::BigInt::Calc::_num
0000s0sMath::BigInt::Calc::::_oneMath::BigInt::Calc::_one
0000s0sMath::BigInt::Calc::::_orMath::BigInt::Calc::_or
0000s0sMath::BigInt::Calc::::_powMath::BigInt::Calc::_pow
0000s0sMath::BigInt::Calc::::_rootMath::BigInt::Calc::_root
0000s0sMath::BigInt::Calc::::_rsftMath::BigInt::Calc::_rsft
0000s0sMath::BigInt::Calc::::_sqrtMath::BigInt::Calc::_sqrt
0000s0sMath::BigInt::Calc::::_subMath::BigInt::Calc::_sub
0000s0sMath::BigInt::Calc::::_tenMath::BigInt::Calc::_ten
0000s0sMath::BigInt::Calc::::_twoMath::BigInt::Calc::_two
0000s0sMath::BigInt::Calc::::_xorMath::BigInt::Calc::_xor
0000s0sMath::BigInt::Calc::::_zerosMath::BigInt::Calc::_zeros
0000s0sMath::BigInt::Calc::::stepsMath::BigInt::Calc::steps
Call graph for these subroutines as a Graphviz dot language file.
Line State
ments		 Time
on line		 Calls Time
in subs		 Code
1package Math::BigInt::Calc;
2
3248µs122µs
# spent 22µs within Math::BigInt::Calc::BEGIN@3 which was called: # once (22µs+0s) by Math::BigInt::BEGIN@1 at line 3
use 5.006002;
# spent 22µs making 1 call to Math::BigInt::Calc::BEGIN@3
42393µs229µs
# spent 18µs (6+11) within Math::BigInt::Calc::BEGIN@4 which was called: # once (6µs+11µs) by Math::BigInt::BEGIN@1 at line 4
use strict;
# spent 18µs making 1 call to Math::BigInt::Calc::BEGIN@4 # spent 12µs making 1 call to strict::import
5# use warnings; # do not use warnings for older Perls
6
71500nsour $VERSION = '1.998';
8
9# Package to store unsigned big integers in decimal and do math with them
10
11# Internally the numbers are stored in an array with at least 1 element, no
12# leading zero parts (except the first) and in base 1eX where X is determined
13# automatically at loading time to be the maximum possible value
14
15# todo:
16# - fully remove funky $# stuff in div() (maybe - that code scares me...)
17
18# USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
19# instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
20# BS2000, some Crays need USE_DIV instead.
21# The BEGIN block is used to determine which of the two variants gives the
22# correct result.
23
24# Beware of things like:
25# $i = $i * $y + $car; $car = int($i / $BASE); $i = $i % $BASE;
26# This works on x86, but fails on ARM (SA1100, iPAQ) due to who knows what
27# reasons. So, use this instead (slower, but correct):
28# $i = $i * $y + $car; $car = int($i / $BASE); $i -= $BASE * $car;
29
30##############################################################################
31# global constants, flags and accessory
32
33# announce that we are compatible with MBI v1.83 and up
34sub api_version () { 2; }
35
36# constants for easier life
371100nsmy ($BASE,$BASE_LEN,$RBASE,$MAX_VAL);
38my ($AND_BITS,$XOR_BITS,$OR_BITS);
39my ($AND_MASK,$XOR_MASK,$OR_MASK);
40
41sub _base_len
42
# spent 11µs within Math::BigInt::Calc::_base_len which was called: # once (11µs+0s) by Math::BigInt::Calc::BEGIN@117 at line 152
{
43 # Set/get the BASE_LEN and assorted other, connected values.
44 # Used only by the testsuite, the set variant is used only by the BEGIN
45 # block below:
461100ns shift;
47
481400ns my ($b, $int) = @_;
491200ns if (defined $b)
50 {
51 # avoid redefinitions
5212µs undef &_mul;
531500ns undef &_div;
54
551600ns if ($] >= 5.008 && $int && $b > 7)
56 {
571300ns $BASE_LEN = $b;
5812µs *_mul = \&_mul_use_div_64;
591700ns *_div = \&_div_use_div_64;
6012µs $BASE = int("1e".$BASE_LEN);
611800ns $MAX_VAL = $BASE-1;
6214µs return $BASE_LEN unless wantarray;
63 return ($BASE_LEN, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL,);
64 }
65
66 # find whether we can use mul or div in mul()/div()
67 $BASE_LEN = $b+1;
68 my $caught = 0;
69 while (--$BASE_LEN > 5)
70 {
71 $BASE = int("1e".$BASE_LEN);
72 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
73 $caught = 0;
74 $caught += 1 if (int($BASE * $RBASE) != 1); # should be 1
75 $caught += 2 if (int($BASE / $BASE) != 1); # should be 1
76 last if $caught != 3;
77 }
78 $BASE = int("1e".$BASE_LEN);
79 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
80 $MAX_VAL = $BASE-1;
81
82 # ($caught & 1) != 0 => cannot use MUL
83 # ($caught & 2) != 0 => cannot use DIV
84 if ($caught == 2) # 2
85 {
86 # must USE_MUL since we cannot use DIV
87 *_mul = \&_mul_use_mul;
88 *_div = \&_div_use_mul;
89 }
90 else # 0 or 1
91 {
92 # can USE_DIV instead
93 *_mul = \&_mul_use_div;
94 *_div = \&_div_use_div;
95 }
96 }
97 return $BASE_LEN unless wantarray;
98 return ($BASE_LEN, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL);
99 }
100
101sub _new
102
# spent 12µs within Math::BigInt::Calc::_new which was called 4 times, avg 3µs/call: # once (5µs+0s) by Math::BigInt::Calc::BEGIN@117 at line 188 # once (3µs+0s) by Math::BigInt::new at line 642 of Math/BigInt.pm # once (2µs+0s) by Math::BigInt::Calc::BEGIN@117 at line 189 # once (2µs+0s) by Math::BigInt::Calc::BEGIN@117 at line 190
{
103 # (ref to string) return ref to num_array
104 # Convert a number from string format (without sign) to internal base
105 # 1ex format. Assumes normalized value as input.
10646µs my $il = length($_[1])-1;
107
108 # < BASE_LEN due len-1 above
109412µs return [ int($_[1]) ] if $il < $BASE_LEN; # shortcut for short numbers
110
111 # this leaves '00000' instead of int 0 and will be corrected after any op
112 [ reverse(unpack("a" . ($il % $BASE_LEN+1)
113 . ("a$BASE_LEN" x ($il / $BASE_LEN)), $_[1])) ];
114 }
115
116BEGIN
117
# spent 203µs (127+77) within Math::BigInt::Calc::BEGIN@117 which was called: # once (127µs+77µs) by Math::BigInt::BEGIN@1 at line 194
{
118 # from Daniel Pfeiffer: determine largest group of digits that is precisely
119 # multipliable with itself plus carry
120 # Test now changed to expect the proper pattern, not a result off by 1 or 2
1211500ns my ($e, $num) = 3; # lowest value we will use is 3+1-1 = 3
122 do
123171µs1440µs {
# spent 37µs making 7 calls to Math::BigInt::Calc::CORE:regcomp, avg 5µs/call # spent 4µs making 7 calls to Math::BigInt::Calc::CORE:match, avg 500ns/call
12475µs $num = ('9' x ++$e) + 0;
12573µs $num *= $num + 1.0;
126 } while ("$num" =~ /9{$e}0{$e}/); # must be a certain pattern
12714µs $e--; # last test failed, so retract one step
128 # the limits below brush the problems with the test above under the rug:
129 # the test should be able to find the proper $e automatically
13014µs12µs $e = 5 if $^O =~ /^uts/; # UTS get's some special treatment
# spent 2µs making 1 call to Math::BigInt::Calc::CORE:match
13112µs1300ns $e = 5 if $^O =~ /^unicos/; # unicos is also problematic (6 seems to work
# spent 300ns making 1 call to Math::BigInt::Calc::CORE:match
132 # there, but we play safe)
133
1341200ns my $int = 0;
1351400ns if ($e > 7)
136 {
137276µs216µs
# spent 14µs (11+2) within Math::BigInt::Calc::BEGIN@137 which was called: # once (11µs+2µs) by Math::BigInt::BEGIN@1 at line 137
use integer;
# spent 14µs making 1 call to Math::BigInt::Calc::BEGIN@137 # spent 2µs making 1 call to integer::import
1381100ns my $e1 = 7;
1391200ns $num = 7;
140 do
141124µs614µs {
# spent 13µs making 3 calls to Math::BigInt::Calc::CORE:regcomp, avg 4µs/call # spent 1µs making 3 calls to Math::BigInt::Calc::CORE:match, avg 367ns/call
14232µs $num = ('9' x ++$e1) + 0;
1433700ns $num *= $num + 1;
144 } while ("$num" =~ /9{$e1}0{$e1}/); # must be a certain pattern
1451100ns $e1--; # last test failed, so retract one step
1461300ns if ($e1 > 7)
147 {
1482300ns $int = 1; $e = $e1;
149 }
150 }
151
15213µs111µs __PACKAGE__->_base_len($e,$int); # set and store
# spent 11µs making 1 call to Math::BigInt::Calc::_base_len
153
154245µs210µs
# spent 9µs (8+1) within Math::BigInt::Calc::BEGIN@154 which was called: # once (8µs+1µs) by Math::BigInt::BEGIN@1 at line 154
use integer;
# spent 9µs making 1 call to Math::BigInt::Calc::BEGIN@154 # spent 1µs making 1 call to integer::import
155 # find out how many bits _and, _or and _xor can take (old default = 16)
156 # I don't think anybody has yet 128 bit scalars, so let's play safe.
15712µs local $^W = 0; # don't warn about 'nonportable number'
1583400ns $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15;
159
160 # find max bits, we will not go higher than numberofbits that fit into $BASE
161 # to make _and etc simpler (and faster for smaller, slower for large numbers)
1621100ns my $max = 16;
163153µs while (2 ** $max < $BASE) { $max++; }
164 {
1653139µs213µs
# spent 10µs (8+2) within Math::BigInt::Calc::BEGIN@165 which was called: # once (8µs+2µs) by Math::BigInt::BEGIN@1 at line 165
no integer;
# spent 10µs making 1 call to Math::BigInt::Calc::BEGIN@165 # spent 2µs making 1 call to integer::unimport
1661200ns $max = 16 if $] < 5.006; # older Perls might not take >16 too well
167 }
1681100ns my ($x,$y,$z);
16914µs do {
170151µs $AND_BITS++;
171308µs $x = CORE::oct('0b' . '1' x $AND_BITS); $y = $x & $x;
172153µs $z = (2 ** $AND_BITS) - 1;
173 } while ($AND_BITS < $max && $x == $z && $y == $x);
1741100ns $AND_BITS --; # retreat one step
17513µs do {
176151µs $XOR_BITS++;
177306µs $x = CORE::oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
178152µs $z = (2 ** $XOR_BITS) - 1;
179 } while ($XOR_BITS < $max && $x == $z && $y == $x);
18010s $XOR_BITS --; # retreat one step
18113µs do {
182151µs $OR_BITS++;
183306µs $x = CORE::oct('0b' . '1' x $OR_BITS); $y = $x | $x;
184152µs $z = (2 ** $OR_BITS) - 1;
185 } while ($OR_BITS < $max && $x == $z && $y == $x);
1861100ns $OR_BITS --; # retreat one step
187
18812µs15µs $AND_MASK = __PACKAGE__->_new( ( 2 ** $AND_BITS ));
# spent 5µs making 1 call to Math::BigInt::Calc::_new
18911µs12µs $XOR_MASK = __PACKAGE__->_new( ( 2 ** $XOR_BITS ));
# spent 2µs making 1 call to Math::BigInt::Calc::_new
19011µs12µs $OR_MASK = __PACKAGE__->_new( ( 2 ** $OR_BITS ));
# spent 2µs making 1 call to Math::BigInt::Calc::_new
191
192 # We can compute the approximate length no faster than the real length:
19314µs *_alen = \&_len;
1941792µs1203µs }
# spent 203µs making 1 call to Math::BigInt::Calc::BEGIN@117
195
196###############################################################################
197
198sub _zero
199
# spent 2µs within Math::BigInt::Calc::_zero which was called: # once (2µs+0s) by Math::BigInt::new at line 589 of Math/BigInt.pm
{
200 # create a zero
20114µs [ 0 ];
202 }
203
204sub _one
205 {
206 # create a one
207 [ 1 ];
208 }
209
210sub _two
211 {
212 # create a two (used internally for shifting)
213 [ 2 ];
214 }
215
216sub _ten
217 {
218 # create a 10 (used internally for shifting)
219 [ 10 ];
220 }
221
222sub _1ex
223 {
224 # create a 1Ex
225 my $rem = $_[1] % $BASE_LEN; # remainder
226 my $parts = $_[1] / $BASE_LEN; # parts
227
228 # 000000, 000000, 100
229 [ (0) x $parts, '1' . ('0' x $rem) ];
230 }
231
232sub _copy
233 {
234 # make a true copy
235 [ @{$_[1]} ];
236 }
237
238# catch and throw away
23914µs
# spent 1µs within Math::BigInt::Calc::import which was called: # once (1µs+0s) by Math::BigInt::BEGIN@1 at line 1 of (eval 48)[Math/BigInt.pm:2820]
sub import { }
240
241##############################################################################
242# convert back to string and number
243
244sub _str
245
# spent 6µs within Math::BigInt::Calc::_str which was called: # once (6µs+0s) by Math::BigInt::bstr at line 836 of Math/BigInt.pm
{
246 # (ref to BINT) return num_str
247 # Convert number from internal base 100000 format to string format.
248 # internal format is always normalized (no leading zeros, "-0" => "+0")
2491400ns my $ar = $_[1];
250
2511300ns my $l = scalar @$ar; # number of parts
2521200ns if ($l < 1) # should not happen
253 {
254 require Carp;
255 Carp::croak("$_[1] has no elements");
256 }
257
2581200ns my $ret = "";
259 # handle first one different to strip leading zeros from it (there are no
260 # leading zero parts in internal representation)
26131µs $l --; $ret .= int($ar->[$l]); $l--;
262 # Interestingly, the pre-padd method uses more time
263 # the old grep variant takes longer (14 vs. 10 sec)
26411µs my $z = '0' x ($BASE_LEN-1);
2651800ns while ($l >= 0)
266 {
267 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
268 $l--;
269 }
27012µs $ret;
271 }
272
273sub _num
274 {
275 # Make a Perl scalar number (int/float) from a BigInt object.
276 my $x = $_[1];
277
278 return 0 + $x->[0] if scalar @$x == 1; # below $BASE
279
280 # Start with the most significant element and work towards the least
281 # significant element. Avoid multiplying "inf" (which happens if the number
282 # overflows) with "0" (if there are zero elements in $x) since this gives
283 # "nan" which propagates to the output.
284
285 my $num = 0;
286 for (my $i = $#$x ; $i >= 0 ; --$i) {
287 $num *= $BASE;
288 $num += $x -> [$i];
289 }
290 return $num;
291 }
292
293##############################################################################
294# actual math code
295
296sub _add
297 {
298 # (ref to int_num_array, ref to int_num_array)
299 # routine to add two base 1eX numbers
300 # stolen from Knuth Vol 2 Algorithm A pg 231
301 # there are separate routines to add and sub as per Knuth pg 233
302 # This routine modifies array x, but not y.
303
304 my ($c,$x,$y) = @_;
305
306 return $x if (@$y == 1) && $y->[0] == 0; # $x + 0 => $x
307 if ((@$x == 1) && $x->[0] == 0) # 0 + $y => $y->copy
308 {
309 # twice as slow as $x = [ @$y ], but nec. to retain $x as ref :(
310 @$x = @$y; return $x;
311 }
312
313 # for each in Y, add Y to X and carry. If after that, something is left in
314 # X, foreach in X add carry to X and then return X, carry
315 # Trades one "$j++" for having to shift arrays
316 my $i; my $car = 0; my $j = 0;
317 for $i (@$y)
318 {
319 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
320 $j++;
321 }
322 while ($car != 0)
323 {
324 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
325 }
326 $x;
327 }
328
329sub _inc
330 {
331 # (ref to int_num_array, ref to int_num_array)
332 # Add 1 to $x, modify $x in place
333 my ($c,$x) = @_;
334
335 for my $i (@$x)
336 {
337 return $x if (($i += 1) < $BASE); # early out
338 $i = 0; # overflow, next
339 }
340 push @$x,1 if (($x->[-1] || 0) == 0); # last overflowed, so extend
341 $x;
342 }
343
344sub _dec
345 {
346 # (ref to int_num_array, ref to int_num_array)
347 # Sub 1 from $x, modify $x in place
348 my ($c,$x) = @_;
349
350 my $MAX = $BASE-1; # since MAX_VAL based on BASE
351 for my $i (@$x)
352 {
353 last if (($i -= 1) >= 0); # early out
354 $i = $MAX; # underflow, next
355 }
356 pop @$x if $x->[-1] == 0 && @$x > 1; # last underflowed (but leave 0)
357 $x;
358 }
359
360sub _sub
361 {
362 # (ref to int_num_array, ref to int_num_array, swap)
363 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
364 # subtract Y from X by modifying x in place
365 my ($c,$sx,$sy,$s) = @_;
366
367 my $car = 0; my $i; my $j = 0;
368 if (!$s)
369 {
370 for $i (@$sx)
371 {
372 last unless defined $sy->[$j] || $car;
373 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
374 }
375 # might leave leading zeros, so fix that
376 return __strip_zeros($sx);
377 }
378 for $i (@$sx)
379 {
380 # we can't do an early out if $x is < than $y, since we
381 # need to copy the high chunks from $y. Found by Bob Mathews.
382 #last unless defined $sy->[$j] || $car;
383 $sy->[$j] += $BASE
384 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
385 $j++;
386 }
387 # might leave leading zeros, so fix that
388 __strip_zeros($sy);
389 }
390
391sub _mul_use_mul
392 {
393 # (ref to int_num_array, ref to int_num_array)
394 # multiply two numbers in internal representation
395 # modifies first arg, second need not be different from first
396 my ($c,$xv,$yv) = @_;
397
398 if (@$yv == 1)
399 {
400 # shortcut for two very short numbers (improved by Nathan Zook)
401 # works also if xv and yv are the same reference, and handles also $x == 0
402 if (@$xv == 1)
403 {
404 if (($xv->[0] *= $yv->[0]) >= $BASE)
405 {
406 $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $BASE;
407 };
408 return $xv;
409 }
410 # $x * 0 => 0
411 if ($yv->[0] == 0)
412 {
413 @$xv = (0);
414 return $xv;
415 }
416 # multiply a large number a by a single element one, so speed up
417 my $y = $yv->[0]; my $car = 0;
418 foreach my $i (@$xv)
419 {
420 $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $BASE;
421 }
422 push @$xv, $car if $car != 0;
423 return $xv;
424 }
425 # shortcut for result $x == 0 => result = 0
426 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
427
428 # since multiplying $x with $x fails, make copy in this case
429 $yv = [@$xv] if $xv == $yv; # same references?
430
431 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
432
433 for $xi (@$xv)
434 {
435 $car = 0; $cty = 0;
436
437 # slow variant
438# for $yi (@$yv)
439# {
440# $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
441# $prod[$cty++] =
442# $prod - ($car = int($prod * RBASE)) * $BASE; # see USE_MUL
443# }
444# $prod[$cty] += $car if $car; # need really to check for 0?
445# $xi = shift @prod;
446
447 # faster variant
448 # looping through this if $xi == 0 is silly - so optimize it away!
449 $xi = (shift @prod || 0), next if $xi == 0;
450 for $yi (@$yv)
451 {
452 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
453## this is actually a tad slower
454## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
455 $prod[$cty++] =
456 $prod - ($car = int($prod * $RBASE)) * $BASE; # see USE_MUL
457 }
458 $prod[$cty] += $car if $car; # need really to check for 0?
459 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
460 }
461 push @$xv, @prod;
462 # can't have leading zeros
463# __strip_zeros($xv);
464 $xv;
465 }
466
467sub _mul_use_div_64
468 {
469 # (ref to int_num_array, ref to int_num_array)
470 # multiply two numbers in internal representation
471 # modifies first arg, second need not be different from first
472 # works for 64 bit integer with "use integer"
473 my ($c,$xv,$yv) = @_;
474
47521.06ms214µs
# spent 13µs (11+2) within Math::BigInt::Calc::BEGIN@475 which was called: # once (11µs+2µs) by Math::BigInt::BEGIN@1 at line 475
use integer;
# spent 13µs making 1 call to Math::BigInt::Calc::BEGIN@475 # spent 2µs making 1 call to integer::import
476 if (@$yv == 1)
477 {
478 # shortcut for two small numbers, also handles $x == 0
479 if (@$xv == 1)
480 {
481 # shortcut for two very short numbers (improved by Nathan Zook)
482 # works also if xv and yv are the same reference, and handles also $x == 0
483 if (($xv->[0] *= $yv->[0]) >= $BASE)
484 {
485 $xv->[0] =
486 $xv->[0] - ($xv->[1] = $xv->[0] / $BASE) * $BASE;
487 };
488 return $xv;
489 }
490 # $x * 0 => 0
491 if ($yv->[0] == 0)
492 {
493 @$xv = (0);
494 return $xv;
495 }
496 # multiply a large number a by a single element one, so speed up
497 my $y = $yv->[0]; my $car = 0;
498 foreach my $i (@$xv)
499 {
500 #$i = $i * $y + $car; $car = $i / $BASE; $i -= $car * $BASE;
501 $i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
502 }
503 push @$xv, $car if $car != 0;
504 return $xv;
505 }
506 # shortcut for result $x == 0 => result = 0
507 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
508
509 # since multiplying $x with $x fails, make copy in this case
510 $yv = [@$xv] if $xv == $yv; # same references?
511
512 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
513 for $xi (@$xv)
514 {
515 $car = 0; $cty = 0;
516 # looping through this if $xi == 0 is silly - so optimize it away!
517 $xi = (shift @prod || 0), next if $xi == 0;
518 for $yi (@$yv)
519 {
520 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
521 $prod[$cty++] = $prod - ($car = $prod / $BASE) * $BASE;
522 }
523 $prod[$cty] += $car if $car; # need really to check for 0?
524 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
525 }
526 push @$xv, @prod;
527 $xv;
528 }
529
530sub _mul_use_div
531 {
532 # (ref to int_num_array, ref to int_num_array)
533 # multiply two numbers in internal representation
534 # modifies first arg, second need not be different from first
535 my ($c,$xv,$yv) = @_;
536
537 if (@$yv == 1)
538 {
539 # shortcut for two small numbers, also handles $x == 0
540 if (@$xv == 1)
541 {
542 # shortcut for two very short numbers (improved by Nathan Zook)
543 # works also if xv and yv are the same reference, and handles also $x == 0
544 if (($xv->[0] *= $yv->[0]) >= $BASE)
545 {
546 $xv->[0] =
547 $xv->[0] - ($xv->[1] = int($xv->[0] / $BASE)) * $BASE;
548 };
549 return $xv;
550 }
551 # $x * 0 => 0
552 if ($yv->[0] == 0)
553 {
554 @$xv = (0);
555 return $xv;
556 }
557 # multiply a large number a by a single element one, so speed up
558 my $y = $yv->[0]; my $car = 0;
559 foreach my $i (@$xv)
560 {
561 $i = $i * $y + $car; $car = int($i / $BASE); $i -= $car * $BASE;
562 # This (together with use integer;) does not work on 32-bit Perls
563 #$i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
564 }
565 push @$xv, $car if $car != 0;
566 return $xv;
567 }
568 # shortcut for result $x == 0 => result = 0
569 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
570
571 # since multiplying $x with $x fails, make copy in this case
572 $yv = [@$xv] if $xv == $yv; # same references?
573
574 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
575 for $xi (@$xv)
576 {
577 $car = 0; $cty = 0;
578 # looping through this if $xi == 0 is silly - so optimize it away!
579 $xi = (shift @prod || 0), next if $xi == 0;
580 for $yi (@$yv)
581 {
582 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
583 $prod[$cty++] = $prod - ($car = int($prod / $BASE)) * $BASE;
584 }
585 $prod[$cty] += $car if $car; # need really to check for 0?
586 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
587 }
588 push @$xv, @prod;
589 # can't have leading zeros
590# __strip_zeros($xv);
591 $xv;
592 }
593
594sub _div_use_mul
595 {
596 # ref to array, ref to array, modify first array and return remainder if
597 # in list context
598
599 # see comments in _div_use_div() for more explanations
600
601 my ($c,$x,$yorg) = @_;
602
603 # the general div algorithm here is about O(N*N) and thus quite slow, so
604 # we first check for some special cases and use shortcuts to handle them.
605
606 # This works, because we store the numbers in a chunked format where each
607 # element contains 5..7 digits (depending on system).
608
609 # if both numbers have only one element:
610 if (@$x == 1 && @$yorg == 1)
611 {
612 # shortcut, $yorg and $x are two small numbers
613 if (wantarray)
614 {
615 my $r = [ $x->[0] % $yorg->[0] ];
616 $x->[0] = int($x->[0] / $yorg->[0]);
617 return ($x,$r);
618 }
619 else
620 {
621 $x->[0] = int($x->[0] / $yorg->[0]);
622 return $x;
623 }
624 }
625
626 # if x has more than one, but y has only one element:
627 if (@$yorg == 1)
628 {
629 my $rem;
630 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
631
632 # shortcut, $y is < $BASE
633 my $j = scalar @$x; my $r = 0;
634 my $y = $yorg->[0]; my $b;
635 while ($j-- > 0)
636 {
637 $b = $r * $BASE + $x->[$j];
638 $x->[$j] = int($b/$y);
639 $r = $b % $y;
640 }
641 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
642 return ($x,$rem) if wantarray;
643 return $x;
644 }
645
646 # now x and y have more than one element
647
648 # check whether y has more elements than x, if yet, the result will be 0
649 if (@$yorg > @$x)
650 {
651 my $rem;
652 $rem = [@$x] if wantarray; # make copy
653 splice (@$x,1); # keep ref to original array
654 $x->[0] = 0; # set to 0
655 return ($x,$rem) if wantarray; # including remainder?
656 return $x; # only x, which is [0] now
657 }
658 # check whether the numbers have the same number of elements, in that case
659 # the result will fit into one element and can be computed efficiently
660 if (@$yorg == @$x)
661 {
662 my $rem;
663 # if $yorg has more digits than $x (it's leading element is longer than
664 # the one from $x), the result will also be 0:
665 if (length(int($yorg->[-1])) > length(int($x->[-1])))
666 {
667 $rem = [@$x] if wantarray; # make copy
668 splice (@$x,1); # keep ref to org array
669 $x->[0] = 0; # set to 0
670 return ($x,$rem) if wantarray; # including remainder?
671 return $x;
672 }
673 # now calculate $x / $yorg
674 if (length(int($yorg->[-1])) == length(int($x->[-1])))
675 {
676 # same length, so make full compare
677
678 my $a = 0; my $j = scalar @$x - 1;
679 # manual way (abort if unequal, good for early ne)
680 while ($j >= 0)
681 {
682 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
683 }
684 # $a contains the result of the compare between X and Y
685 # a < 0: x < y, a == 0: x == y, a > 0: x > y
686 if ($a <= 0)
687 {
688 $rem = [ 0 ]; # a = 0 => x == y => rem 0
689 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
690 splice(@$x,1); # keep single element
691 $x->[0] = 0; # if $a < 0
692 $x->[0] = 1 if $a == 0; # $x == $y
693 return ($x,$rem) if wantarray;
694 return $x;
695 }
696 # $x >= $y, so proceed normally
697 }
698 }
699
700 # all other cases:
701
702 my $y = [ @$yorg ]; # always make copy to preserve
703
704 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
705
706 $car = $bar = $prd = 0;
707 if (($dd = int($BASE/($y->[-1]+1))) != 1)
708 {
709 for $xi (@$x)
710 {
711 $xi = $xi * $dd + $car;
712 $xi -= ($car = int($xi * $RBASE)) * $BASE; # see USE_MUL
713 }
714 push(@$x, $car); $car = 0;
715 for $yi (@$y)
716 {
717 $yi = $yi * $dd + $car;
718 $yi -= ($car = int($yi * $RBASE)) * $BASE; # see USE_MUL
719 }
720 }
721 else
722 {
723 push(@$x, 0);
724 }
725 @q = (); ($v2,$v1) = @$y[-2,-1];
726 $v2 = 0 unless $v2;
727 while ($#$x > $#$y)
728 {
729 ($u2,$u1,$u0) = @$x[-3..-1];
730 $u2 = 0 unless $u2;
731 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
732 # if $v1 == 0;
733 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
734 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
735 if ($q)
736 {
737 ($car, $bar) = (0,0);
738 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
739 {
740 $prd = $q * $y->[$yi] + $car;
741 $prd -= ($car = int($prd * $RBASE)) * $BASE; # see USE_MUL
742 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
743 }
744 if ($x->[-1] < $car + $bar)
745 {
746 $car = 0; --$q;
747 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
748 {
749 $x->[$xi] -= $BASE
750 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
751 }
752 }
753 }
754 pop(@$x);
755 unshift(@q, $q);
756 }
757 if (wantarray)
758 {
759 @d = ();
760 if ($dd != 1)
761 {
762 $car = 0;
763 for $xi (reverse @$x)
764 {
765 $prd = $car * $BASE + $xi;
766 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
767 unshift(@d, $tmp);
768 }
769 }
770 else
771 {
772 @d = @$x;
773 }
774 @$x = @q;
775 my $d = \@d;
776 __strip_zeros($x);
777 __strip_zeros($d);
778 return ($x,$d);
779 }
780 @$x = @q;
781 __strip_zeros($x);
782 $x;
783 }
784
785sub _div_use_div_64
786 {
787 # ref to array, ref to array, modify first array and return remainder if
788 # in list context
789 # This version works on 64 bit integers
790 my ($c,$x,$yorg) = @_;
791
79223.79ms211µs
# spent 10µs (8+2) within Math::BigInt::Calc::BEGIN@792 which was called: # once (8µs+2µs) by Math::BigInt::BEGIN@1 at line 792
use integer;
# spent 10µs making 1 call to Math::BigInt::Calc::BEGIN@792 # spent 2µs making 1 call to integer::import
793 # the general div algorithm here is about O(N*N) and thus quite slow, so
794 # we first check for some special cases and use shortcuts to handle them.
795
796 # This works, because we store the numbers in a chunked format where each
797 # element contains 5..7 digits (depending on system).
798
799 # if both numbers have only one element:
800 if (@$x == 1 && @$yorg == 1)
801 {
802 # shortcut, $yorg and $x are two small numbers
803 if (wantarray)
804 {
805 my $r = [ $x->[0] % $yorg->[0] ];
806 $x->[0] = int($x->[0] / $yorg->[0]);
807 return ($x,$r);
808 }
809 else
810 {
811 $x->[0] = int($x->[0] / $yorg->[0]);
812 return $x;
813 }
814 }
815 # if x has more than one, but y has only one element:
816 if (@$yorg == 1)
817 {
818 my $rem;
819 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
820
821 # shortcut, $y is < $BASE
822 my $j = scalar @$x; my $r = 0;
823 my $y = $yorg->[0]; my $b;
824 while ($j-- > 0)
825 {
826 $b = $r * $BASE + $x->[$j];
827 $x->[$j] = int($b/$y);
828 $r = $b % $y;
829 }
830 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
831 return ($x,$rem) if wantarray;
832 return $x;
833 }
834 # now x and y have more than one element
835
836 # check whether y has more elements than x, if yet, the result will be 0
837 if (@$yorg > @$x)
838 {
839 my $rem;
840 $rem = [@$x] if wantarray; # make copy
841 splice (@$x,1); # keep ref to original array
842 $x->[0] = 0; # set to 0
843 return ($x,$rem) if wantarray; # including remainder?
844 return $x; # only x, which is [0] now
845 }
846 # check whether the numbers have the same number of elements, in that case
847 # the result will fit into one element and can be computed efficiently
848 if (@$yorg == @$x)
849 {
850 my $rem;
851 # if $yorg has more digits than $x (it's leading element is longer than
852 # the one from $x), the result will also be 0:
853 if (length(int($yorg->[-1])) > length(int($x->[-1])))
854 {
855 $rem = [@$x] if wantarray; # make copy
856 splice (@$x,1); # keep ref to org array
857 $x->[0] = 0; # set to 0
858 return ($x,$rem) if wantarray; # including remainder?
859 return $x;
860 }
861 # now calculate $x / $yorg
862
863 if (length(int($yorg->[-1])) == length(int($x->[-1])))
864 {
865 # same length, so make full compare
866
867 my $a = 0; my $j = scalar @$x - 1;
868 # manual way (abort if unequal, good for early ne)
869 while ($j >= 0)
870 {
871 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
872 }
873 # $a contains the result of the compare between X and Y
874 # a < 0: x < y, a == 0: x == y, a > 0: x > y
875 if ($a <= 0)
876 {
877 $rem = [ 0 ]; # a = 0 => x == y => rem 0
878 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
879 splice(@$x,1); # keep single element
880 $x->[0] = 0; # if $a < 0
881 $x->[0] = 1 if $a == 0; # $x == $y
882 return ($x,$rem) if wantarray; # including remainder?
883 return $x;
884 }
885 # $x >= $y, so proceed normally
886
887 }
888 }
889
890 # all other cases:
891
892 my $y = [ @$yorg ]; # always make copy to preserve
893
894 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
895
896 $car = $bar = $prd = 0;
897 if (($dd = int($BASE/($y->[-1]+1))) != 1)
898 {
899 for $xi (@$x)
900 {
901 $xi = $xi * $dd + $car;
902 $xi -= ($car = int($xi / $BASE)) * $BASE;
903 }
904 push(@$x, $car); $car = 0;
905 for $yi (@$y)
906 {
907 $yi = $yi * $dd + $car;
908 $yi -= ($car = int($yi / $BASE)) * $BASE;
909 }
910 }
911 else
912 {
913 push(@$x, 0);
914 }
915
916 # @q will accumulate the final result, $q contains the current computed
917 # part of the final result
918
919 @q = (); ($v2,$v1) = @$y[-2,-1];
920 $v2 = 0 unless $v2;
921 while ($#$x > $#$y)
922 {
923 ($u2,$u1,$u0) = @$x[-3..-1];
924 $u2 = 0 unless $u2;
925 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
926 # if $v1 == 0;
927 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
928 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
929 if ($q)
930 {
931 ($car, $bar) = (0,0);
932 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
933 {
934 $prd = $q * $y->[$yi] + $car;
935 $prd -= ($car = int($prd / $BASE)) * $BASE;
936 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
937 }
938 if ($x->[-1] < $car + $bar)
939 {
940 $car = 0; --$q;
941 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
942 {
943 $x->[$xi] -= $BASE
944 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
945 }
946 }
947 }
948 pop(@$x); unshift(@q, $q);
949 }
950 if (wantarray)
951 {
952 @d = ();
953 if ($dd != 1)
954 {
955 $car = 0;
956 for $xi (reverse @$x)
957 {
958 $prd = $car * $BASE + $xi;
959 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
960 unshift(@d, $tmp);
961 }
962 }
963 else
964 {
965 @d = @$x;
966 }
967 @$x = @q;
968 my $d = \@d;
969 __strip_zeros($x);
970 __strip_zeros($d);
971 return ($x,$d);
972 }
973 @$x = @q;
974 __strip_zeros($x);
975 $x;
976 }
977
978sub _div_use_div
979 {
980 # ref to array, ref to array, modify first array and return remainder if
981 # in list context
982 my ($c,$x,$yorg) = @_;
983
984 # the general div algorithm here is about O(N*N) and thus quite slow, so
985 # we first check for some special cases and use shortcuts to handle them.
986
987 # This works, because we store the numbers in a chunked format where each
988 # element contains 5..7 digits (depending on system).
989
990 # if both numbers have only one element:
991 if (@$x == 1 && @$yorg == 1)
992 {
993 # shortcut, $yorg and $x are two small numbers
994 if (wantarray)
995 {
996 my $r = [ $x->[0] % $yorg->[0] ];
997 $x->[0] = int($x->[0] / $yorg->[0]);
998 return ($x,$r);
999 }
1000 else
1001 {
1002 $x->[0] = int($x->[0] / $yorg->[0]);
1003 return $x;
1004 }
1005 }
1006 # if x has more than one, but y has only one element:
1007 if (@$yorg == 1)
1008 {
1009 my $rem;
1010 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
1011
1012 # shortcut, $y is < $BASE
1013 my $j = scalar @$x; my $r = 0;
1014 my $y = $yorg->[0]; my $b;
1015 while ($j-- > 0)
1016 {
1017 $b = $r * $BASE + $x->[$j];
1018 $x->[$j] = int($b/$y);
1019 $r = $b % $y;
1020 }
1021 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
1022 return ($x,$rem) if wantarray;
1023 return $x;
1024 }
1025 # now x and y have more than one element
1026
1027 # check whether y has more elements than x, if yet, the result will be 0
1028 if (@$yorg > @$x)
1029 {
1030 my $rem;
1031 $rem = [@$x] if wantarray; # make copy
1032 splice (@$x,1); # keep ref to original array
1033 $x->[0] = 0; # set to 0
1034 return ($x,$rem) if wantarray; # including remainder?
1035 return $x; # only x, which is [0] now
1036 }
1037 # check whether the numbers have the same number of elements, in that case
1038 # the result will fit into one element and can be computed efficiently
1039 if (@$yorg == @$x)
1040 {
1041 my $rem;
1042 # if $yorg has more digits than $x (it's leading element is longer than
1043 # the one from $x), the result will also be 0:
1044 if (length(int($yorg->[-1])) > length(int($x->[-1])))
1045 {
1046 $rem = [@$x] if wantarray; # make copy
1047 splice (@$x,1); # keep ref to org array
1048 $x->[0] = 0; # set to 0
1049 return ($x,$rem) if wantarray; # including remainder?
1050 return $x;
1051 }
1052 # now calculate $x / $yorg
1053
1054 if (length(int($yorg->[-1])) == length(int($x->[-1])))
1055 {
1056 # same length, so make full compare
1057
1058 my $a = 0; my $j = scalar @$x - 1;
1059 # manual way (abort if unequal, good for early ne)
1060 while ($j >= 0)
1061 {
1062 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
1063 }
1064 # $a contains the result of the compare between X and Y
1065 # a < 0: x < y, a == 0: x == y, a > 0: x > y
1066 if ($a <= 0)
1067 {
1068 $rem = [ 0 ]; # a = 0 => x == y => rem 0
1069 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
1070 splice(@$x,1); # keep single element
1071 $x->[0] = 0; # if $a < 0
1072 $x->[0] = 1 if $a == 0; # $x == $y
1073 return ($x,$rem) if wantarray; # including remainder?
1074 return $x;
1075 }
1076 # $x >= $y, so proceed normally
1077
1078 }
1079 }
1080
1081 # all other cases:
1082
1083 my $y = [ @$yorg ]; # always make copy to preserve
1084
1085 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
1086
1087 $car = $bar = $prd = 0;
1088 if (($dd = int($BASE/($y->[-1]+1))) != 1)
1089 {
1090 for $xi (@$x)
1091 {
1092 $xi = $xi * $dd + $car;
1093 $xi -= ($car = int($xi / $BASE)) * $BASE;
1094 }
1095 push(@$x, $car); $car = 0;
1096 for $yi (@$y)
1097 {
1098 $yi = $yi * $dd + $car;
1099 $yi -= ($car = int($yi / $BASE)) * $BASE;
1100 }
1101 }
1102 else
1103 {
1104 push(@$x, 0);
1105 }
1106
1107 # @q will accumulate the final result, $q contains the current computed
1108 # part of the final result
1109
1110 @q = (); ($v2,$v1) = @$y[-2,-1];
1111 $v2 = 0 unless $v2;
1112 while ($#$x > $#$y)
1113 {
1114 ($u2,$u1,$u0) = @$x[-3..-1];
1115 $u2 = 0 unless $u2;
1116 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
1117 # if $v1 == 0;
1118 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
1119 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
1120 if ($q)
1121 {
1122 ($car, $bar) = (0,0);
1123 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
1124 {
1125 $prd = $q * $y->[$yi] + $car;
1126 $prd -= ($car = int($prd / $BASE)) * $BASE;
1127 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
1128 }
1129 if ($x->[-1] < $car + $bar)
1130 {
1131 $car = 0; --$q;
1132 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
1133 {
1134 $x->[$xi] -= $BASE
1135 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
1136 }
1137 }
1138 }
1139 pop(@$x); unshift(@q, $q);
1140 }
1141 if (wantarray)
1142 {
1143 @d = ();
1144 if ($dd != 1)
1145 {
1146 $car = 0;
1147 for $xi (reverse @$x)
1148 {
1149 $prd = $car * $BASE + $xi;
1150 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
1151 unshift(@d, $tmp);
1152 }
1153 }
1154 else
1155 {
1156 @d = @$x;
1157 }
1158 @$x = @q;
1159 my $d = \@d;
1160 __strip_zeros($x);
1161 __strip_zeros($d);
1162 return ($x,$d);
1163 }
1164 @$x = @q;
1165 __strip_zeros($x);
1166 $x;
1167 }
1168
1169##############################################################################
1170# testing
1171
1172sub _acmp
1173 {
1174 # internal absolute post-normalized compare (ignore signs)
1175 # ref to array, ref to array, return <0, 0, >0
1176 # arrays must have at least one entry; this is not checked for
1177 my ($c,$cx,$cy) = @_;
1178
1179 # shortcut for short numbers
1180 return (($cx->[0] <=> $cy->[0]) <=> 0)
1181 if scalar @$cx == scalar @$cy && scalar @$cx == 1;
1182
1183 # fast comp based on number of array elements (aka pseudo-length)
1184 my $lxy = (scalar @$cx - scalar @$cy)
1185 # or length of first element if same number of elements (aka difference 0)
1186 ||
1187 # need int() here because sometimes the last element is '00018' vs '18'
1188 (length(int($cx->[-1])) - length(int($cy->[-1])));
1189 return -1 if $lxy < 0; # already differs, ret
1190 return 1 if $lxy > 0; # ditto
1191
1192 # manual way (abort if unequal, good for early ne)
1193 my $a; my $j = scalar @$cx;
1194 while (--$j >= 0)
1195 {
1196 last if ($a = $cx->[$j] - $cy->[$j]);
1197 }
1198 $a <=> 0;
1199 }
1200
1201sub _len
1202 {
1203 # compute number of digits in base 10
1204
1205 # int() because add/sub sometimes leaves strings (like '00005') instead of
1206 # '5' in this place, thus causing length() to report wrong length
1207 my $cx = $_[1];
1208
1209 (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
1210 }
1211
1212sub _digit
1213 {
1214 # Return the nth digit. Zero is rightmost, so _digit(123,0) gives 3.
1215 # Negative values count from the left, so _digit(123, -1) gives 1.
1216 my ($c,$x,$n) = @_;
1217
1218 my $len = _len('',$x);
1219
1220 $n += $len if $n < 0; # -1 last, -2 second-to-last
1221 return "0" if $n < 0 || $n >= $len; # return 0 for digits out of range
1222
1223 my $elem = int($n / $BASE_LEN); # which array element
1224 my $digit = $n % $BASE_LEN; # which digit in this element
1225 substr("$x->[$elem]", -$digit-1, 1);
1226 }
1227
1228sub _zeros
1229 {
1230 # return amount of trailing zeros in decimal
1231 # check each array elem in _m for having 0 at end as long as elem == 0
1232 # Upon finding a elem != 0, stop
1233 my $x = $_[1];
1234
1235 return 0 if scalar @$x == 1 && $x->[0] == 0;
1236
1237 my $zeros = 0; my $elem;
1238 foreach my $e (@$x)
1239 {
1240 if ($e != 0)
1241 {
1242 $elem = "$e"; # preserve x
1243 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
1244 $zeros *= $BASE_LEN; # elems * 5
1245 $zeros += length($elem); # count trailing zeros
1246 last; # early out
1247 }
1248 $zeros ++; # real else branch: 50% slower!
1249 }
1250 $zeros;
1251 }
1252
1253##############################################################################
1254# _is_* routines
1255
1256sub _is_zero
1257 {
1258 # return true if arg is zero
1259 (((scalar @{$_[1]} == 1) && ($_[1]->[0] == 0))) <=> 0;
1260 }
1261
1262sub _is_even
1263 {
1264 # return true if arg is even
1265 (!($_[1]->[0] & 1)) <=> 0;
1266 }
1267
1268sub _is_odd
1269 {
1270 # return true if arg is odd
1271 (($_[1]->[0] & 1)) <=> 0;
1272 }
1273
1274sub _is_one
1275 {
1276 # return true if arg is one
1277 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 1) <=> 0;
1278 }
1279
1280sub _is_two
1281 {
1282 # return true if arg is two
1283 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 2) <=> 0;
1284 }
1285
1286sub _is_ten
1287 {
1288 # return true if arg is ten
1289 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 10) <=> 0;
1290 }
1291
1292sub __strip_zeros
1293 {
1294 # internal normalization function that strips leading zeros from the array
1295 # args: ref to array
1296 my $s = shift;
1297
1298 my $cnt = scalar @$s; # get count of parts
1299 my $i = $cnt-1;
1300 push @$s,0 if $i < 0; # div might return empty results, so fix it
1301
1302 return $s if @$s == 1; # early out
1303
1304 #print "strip: cnt $cnt i $i\n";
1305 # '0', '3', '4', '0', '0',
1306 # 0 1 2 3 4
1307 # cnt = 5, i = 4
1308 # i = 4
1309 # i = 3
1310 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1311 # >= 1: skip first part (this can be zero)
1312 while ($i > 0) { last if $s->[$i] != 0; $i--; }
1313 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1314 $s;
1315 }
1316
1317###############################################################################
1318# check routine to test internal state for corruptions
1319
1320sub _check
1321 {
1322 # used by the test suite
1323 my $x = $_[1];
1324
1325 return "$x is not a reference" if !ref($x);
1326
1327 # are all parts are valid?
1328 my $i = 0; my $j = scalar @$x; my ($e,$try);
1329 while ($i < $j)
1330 {
1331 $e = $x->[$i]; $e = 'undef' unless defined $e;
1332 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1333 last if $e !~ /^[+]?[0-9]+$/;
1334 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1335 last if "$e" !~ /^[+]?[0-9]+$/;
1336 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1337 last if '' . "$e" !~ /^[+]?[0-9]+$/;
1338 $try = ' < 0 || >= $BASE; '."($x, $e)";
1339 last if $e <0 || $e >= $BASE;
1340 # this test is disabled, since new/bnorm and certain ops (like early out
1341 # in add/sub) are allowed/expected to leave '00000' in some elements
1342 #$try = '=~ /^00+/; '."($x, $e)";
1343 #last if $e =~ /^00+/;
1344 $i++;
1345 }
1346 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1347 0;
1348 }
1349
1350
1351###############################################################################
1352
1353sub _mod
1354 {
1355 # if possible, use mod shortcut
1356 my ($c,$x,$yo) = @_;
1357
1358 # slow way since $y too big
1359 if (scalar @$yo > 1)
1360 {
1361 my ($xo,$rem) = _div($c,$x,$yo);
1362 @$x = @$rem;
1363 return $x;
1364 }
1365
1366 my $y = $yo->[0];
1367
1368 # if both are single element arrays
1369 if (scalar @$x == 1)
1370 {
1371 $x->[0] %= $y;
1372 return $x;
1373 }
1374
1375 # if @$x has more than one element, but @$y is a single element
1376 my $b = $BASE % $y;
1377 if ($b == 0)
1378 {
1379 # when BASE % Y == 0 then (B * BASE) % Y == 0
1380 # (B * BASE) % $y + A % Y => A % Y
1381 # so need to consider only last element: O(1)
1382 $x->[0] %= $y;
1383 }
1384 elsif ($b == 1)
1385 {
1386 # else need to go through all elements in @$x: O(N), but loop is a bit
1387 # simplified
1388 my $r = 0;
1389 foreach (@$x)
1390 {
1391 $r = ($r + $_) % $y; # not much faster, but heh...
1392 #$r += $_ % $y; $r %= $y;
1393 }
1394 $r = 0 if $r == $y;
1395 $x->[0] = $r;
1396 }
1397 else
1398 {
1399 # else need to go through all elements in @$x: O(N)
1400 my $r = 0;
1401 my $bm = 1;
1402 foreach (@$x)
1403 {
1404 $r = ($_ * $bm + $r) % $y;
1405 $bm = ($bm * $b) % $y;
1406
1407 #$r += ($_ % $y) * $bm;
1408 #$bm *= $b;
1409 #$bm %= $y;
1410 #$r %= $y;
1411 }
1412 $r = 0 if $r == $y;
1413 $x->[0] = $r;
1414 }
1415 @$x = $x->[0]; # keep one element of @$x
1416 return $x;
1417 }
1418
1419##############################################################################
1420# shifts
1421
1422sub _rsft
1423 {
1424 my ($c,$x,$y,$n) = @_;
1425
1426 if ($n != 10)
1427 {
1428 $n = _new($c,$n); return _div($c,$x, _pow($c,$n,$y));
1429 }
1430
1431 # shortcut (faster) for shifting by 10)
1432 # multiples of $BASE_LEN
1433 my $dst = 0; # destination
1434 my $src = _num($c,$y); # as normal int
1435 my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1])); # len of x in digits
1436 if ($src >= $xlen or ($src == $xlen and ! defined $x->[1]))
1437 {
1438 # 12345 67890 shifted right by more than 10 digits => 0
1439 splice (@$x,1); # leave only one element
1440 $x->[0] = 0; # set to zero
1441 return $x;
1442 }
1443 my $rem = $src % $BASE_LEN; # remainder to shift
1444 $src = int($src / $BASE_LEN); # source
1445 if ($rem == 0)
1446 {
1447 splice (@$x,0,$src); # even faster, 38.4 => 39.3
1448 }
1449 else
1450 {
1451 my $len = scalar @$x - $src; # elems to go
1452 my $vd; my $z = '0'x $BASE_LEN;
1453 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1454 while ($dst < $len)
1455 {
1456 $vd = $z.$x->[$src];
1457 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1458 $src++;
1459 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1460 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1461 $x->[$dst] = int($vd);
1462 $dst++;
1463 }
1464 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
1465 pop @$x if $x->[-1] == 0 && @$x > 1; # kill last element if 0
1466 } # else rem == 0
1467 $x;
1468 }
1469
1470sub _lsft
1471 {
1472 my ($c,$x,$y,$n) = @_;
1473
1474 if ($n != 10)
1475 {
1476 $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
1477 }
1478
1479 # shortcut (faster) for shifting by 10) since we are in base 10eX
1480 # multiples of $BASE_LEN:
1481 my $src = scalar @$x; # source
1482 my $len = _num($c,$y); # shift-len as normal int
1483 my $rem = $len % $BASE_LEN; # remainder to shift
1484 my $dst = $src + int($len/$BASE_LEN); # destination
1485 my $vd; # further speedup
1486 $x->[$src] = 0; # avoid first ||0 for speed
1487 my $z = '0' x $BASE_LEN;
1488 while ($src >= 0)
1489 {
1490 $vd = $x->[$src]; $vd = $z.$vd;
1491 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1492 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1493 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1494 $x->[$dst] = int($vd);
1495 $dst--; $src--;
1496 }
1497 # set lowest parts to 0
1498 while ($dst >= 0) { $x->[$dst--] = 0; }
1499 # fix spurious last zero element
1500 splice @$x,-1 if $x->[-1] == 0;
1501 $x;
1502 }
1503
1504sub _pow
1505 {
1506 # power of $x to $y
1507 # ref to array, ref to array, return ref to array
1508 my ($c,$cx,$cy) = @_;
1509
1510 if (scalar @$cy == 1 && $cy->[0] == 0)
1511 {
1512 splice (@$cx,1); $cx->[0] = 1; # y == 0 => x => 1
1513 return $cx;
1514 }
1515 if ((scalar @$cx == 1 && $cx->[0] == 1) || # x == 1
1516 (scalar @$cy == 1 && $cy->[0] == 1)) # or y == 1
1517 {
1518 return $cx;
1519 }
1520 if (scalar @$cx == 1 && $cx->[0] == 0)
1521 {
1522 splice (@$cx,1); $cx->[0] = 0; # 0 ** y => 0 (if not y <= 0)
1523 return $cx;
1524 }
1525
1526 my $pow2 = _one();
1527
1528 my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
1529 my $len = length($y_bin);
1530 while (--$len > 0)
1531 {
1532 _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1'; # is odd?
1533 _mul($c,$cx,$cx);
1534 }
1535
1536 _mul($c,$cx,$pow2);
1537 $cx;
1538 }
1539
1540sub _nok {
1541 # Return binomial coefficient (n over k).
1542 # Given refs to arrays, return ref to array.
1543 # First input argument is modified.
1544
1545 my ($c, $n, $k) = @_;
1546
1547 # If k > n/2, or, equivalently, 2*k > n, compute nok(n, k) as
1548 # nok(n, n-k), to minimize the number if iterations in the loop.
1549
1550 {
1551 my $twok = _mul($c, _two($c), _copy($c, $k)); # 2 * k
1552 if (_acmp($c, $twok, $n) > 0) { # if 2*k > n
1553 $k = _sub($c, _copy($c, $n), $k); # k = n - k
1554 }
1555 }
1556
1557 # Example:
1558 #
1559 # / 7 \ 7! 1*2*3*4 * 5*6*7 5 * 6 * 7 6 7
1560 # | | = --------- = --------------- = --------- = 5 * - * -
1561 # \ 3 / (7-3)! 3! 1*2*3*4 * 1*2*3 1 * 2 * 3 2 3
1562
1563 if (_is_zero($c, $k)) {
1564 @$n = 1;
1565 }
1566
1567 else {
1568
1569 # Make a copy of the original n, since we'll be modifying n in-place.
1570
1571 my $n_orig = _copy($c, $n);
1572
1573 # n = 5, f = 6, d = 2 (cf. example above)
1574
1575 _sub($c, $n, $k);
1576 _inc($c, $n);
1577
1578 my $f = _copy($c, $n);
1579 _inc($c, $f);
1580
1581 my $d = _two($c);
1582
1583 # while f <= n (the original n, that is) ...
1584
1585 while (_acmp($c, $f, $n_orig) <= 0) {
1586
1587 # n = (n * f / d) == 5 * 6 / 2 (cf. example above)
1588
1589 _mul($c, $n, $f);
1590 _div($c, $n, $d);
1591
1592 # f = 7, d = 3 (cf. example above)
1593
1594 _inc($c, $f);
1595 _inc($c, $d);
1596 }
1597
1598 }
1599
1600 return $n;
1601}
1602
16031900nsmy @factorials = (
1604 1,
1605 1,
1606 2,
1607 2*3,
1608 2*3*4,
1609 2*3*4*5,
1610 2*3*4*5*6,
1611 2*3*4*5*6*7,
1612);
1613
1614sub _fac
1615 {
1616 # factorial of $x
1617 # ref to array, return ref to array
1618 my ($c,$cx) = @_;
1619
1620 if ((@$cx == 1) && ($cx->[0] <= 7))
1621 {
1622 $cx->[0] = $factorials[$cx->[0]]; # 0 => 1, 1 => 1, 2 => 2 etc.
1623 return $cx;
1624 }
1625
1626 if ((@$cx == 1) && # we do this only if $x >= 12 and $x <= 7000
1627 ($cx->[0] >= 12 && $cx->[0] < 7000))
1628 {
1629
1630 # Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
1631 # See http://blogten.blogspot.com/2007/01/calculating-n.html
1632 # The above series can be expressed as factors:
1633 # k * k - (j - i) * 2
1634 # We cache k*k, and calculate (j * j) as the sum of the first j odd integers
1635
1636 # This will not work when N exceeds the storage of a Perl scalar, however,
1637 # in this case the algorithm would be way to slow to terminate, anyway.
1638
1639 # As soon as the last element of $cx is 0, we split it up and remember
1640 # how many zeors we got so far. The reason is that n! will accumulate
1641 # zeros at the end rather fast.
1642 my $zero_elements = 0;
1643
1644 # If n is even, set n = n -1
1645 my $k = _num($c,$cx); my $even = 1;
1646 if (($k & 1) == 0)
1647 {
1648 $even = $k; $k --;
1649 }
1650 # set k to the center point
1651 $k = ($k + 1) / 2;
1652# print "k $k even: $even\n";
1653 # now calculate k * k
1654 my $k2 = $k * $k;
1655 my $odd = 1; my $sum = 1;
1656 my $i = $k - 1;
1657 # keep reference to x
1658 my $new_x = _new($c, $k * $even);
1659 @$cx = @$new_x;
1660 if ($cx->[0] == 0)
1661 {
1662 $zero_elements ++; shift @$cx;
1663 }
1664# print STDERR "x = ", _str($c,$cx),"\n";
1665 my $BASE2 = int(sqrt($BASE))-1;
1666 my $j = 1;
1667 while ($j <= $i)
1668 {
1669 my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++;
1670 while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2))
1671 {
1672 $m *= ($k2 - $sum);
1673 $odd += 2; $sum += $odd; $j++;
1674# print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
1675 }
1676 if ($m < $BASE)
1677 {
1678 _mul($c,$cx,[$m]);
1679 }
1680 else
1681 {
1682 _mul($c,$cx,$c->_new($m));
1683 }
1684 if ($cx->[0] == 0)
1685 {
1686 $zero_elements ++; shift @$cx;
1687 }
1688# print STDERR "Calculate $k2 - $sum = $m (x = ", _str($c,$cx),")\n";
1689 }
1690 # multiply in the zeros again
1691 unshift @$cx, (0) x $zero_elements;
1692 return $cx;
1693 }
1694
1695 # go forward until $base is exceeded
1696 # limit is either $x steps (steps == 100 means a result always too high) or
1697 # $base.
1698 my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1699 my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1700 while ($r*$cf < $BASE && $step < $steps)
1701 {
1702 $last = $r; $r *= $cf++; $step++;
1703 }
1704 if ((@$cx == 1) && $step == $cx->[0])
1705 {
1706 # completely done, so keep reference to $x and return
1707 $cx->[0] = $r;
1708 return $cx;
1709 }
1710
1711 # now we must do the left over steps
1712 my $n; # steps still to do
1713 if (scalar @$cx == 1)
1714 {
1715 $n = $cx->[0];
1716 }
1717 else
1718 {
1719 $n = _copy($c,$cx);
1720 }
1721
1722 # Set $cx to the last result below $BASE (but keep ref to $x)
1723 $cx->[0] = $last; splice (@$cx,1);
1724 # As soon as the last element of $cx is 0, we split it up and remember
1725 # how many zeors we got so far. The reason is that n! will accumulate
1726 # zeros at the end rather fast.
1727 my $zero_elements = 0;
1728
1729 # do left-over steps fit into a scalar?
1730 if (ref $n eq 'ARRAY')
1731 {
1732 # No, so use slower inc() & cmp()
1733 # ($n is at least $BASE here)
1734 my $base_2 = int(sqrt($BASE)) - 1;
1735 #print STDERR "base_2: $base_2\n";
1736 while ($step < $base_2)
1737 {
1738 if ($cx->[0] == 0)
1739 {
1740 $zero_elements ++; shift @$cx;
1741 }
1742 my $b = $step * ($step + 1); $step += 2;
1743 _mul($c,$cx,[$b]);
1744 }
1745 $step = [$step];
1746 while (_acmp($c,$step,$n) <= 0)
1747 {
1748 if ($cx->[0] == 0)
1749 {
1750 $zero_elements ++; shift @$cx;
1751 }
1752 _mul($c,$cx,$step); _inc($c,$step);
1753 }
1754 }
1755 else
1756 {
1757 # Yes, so we can speed it up slightly
1758
1759# print "# left over steps $n\n";
1760
1761 my $base_4 = int(sqrt(sqrt($BASE))) - 2;
1762 #print STDERR "base_4: $base_4\n";
1763 my $n4 = $n - 4;
1764 while ($step < $n4 && $step < $base_4)
1765 {
1766 if ($cx->[0] == 0)
1767 {
1768 $zero_elements ++; shift @$cx;
1769 }
1770 my $b = $step * ($step + 1); $step += 2; $b *= $step * ($step + 1); $step += 2;
1771 _mul($c,$cx,[$b]);
1772 }
1773 my $base_2 = int(sqrt($BASE)) - 1;
1774 my $n2 = $n - 2;
1775 #print STDERR "base_2: $base_2\n";
1776 while ($step < $n2 && $step < $base_2)
1777 {
1778 if ($cx->[0] == 0)
1779 {
1780 $zero_elements ++; shift @$cx;
1781 }
1782 my $b = $step * ($step + 1); $step += 2;
1783 _mul($c,$cx,[$b]);
1784 }
1785 # do what's left over
1786 while ($step <= $n)
1787 {
1788 _mul($c,$cx,[$step]); $step++;
1789 if ($cx->[0] == 0)
1790 {
1791 $zero_elements ++; shift @$cx;
1792 }
1793 }
1794 }
1795 # multiply in the zeros again
1796 unshift @$cx, (0) x $zero_elements;
1797 $cx; # return result
1798 }
1799
1800#############################################################################
1801
1802sub _log_int
1803 {
1804 # calculate integer log of $x to base $base
1805 # ref to array, ref to array - return ref to array
1806 my ($c,$x,$base) = @_;
1807
1808 # X == 0 => NaN
1809 return if (scalar @$x == 1 && $x->[0] == 0);
1810 # BASE 0 or 1 => NaN
1811 return if (scalar @$base == 1 && $base->[0] < 2);
1812 my $cmp = _acmp($c,$x,$base); # X == BASE => 1
1813 if ($cmp == 0)
1814 {
1815 splice (@$x,1); $x->[0] = 1;
1816 return ($x,1)
1817 }
1818 # X < BASE
1819 if ($cmp < 0)
1820 {
1821 splice (@$x,1); $x->[0] = 0;
1822 return ($x,undef);
1823 }
1824
1825 my $x_org = _copy($c,$x); # preserve x
1826 splice(@$x,1); $x->[0] = 1; # keep ref to $x
1827
1828 # Compute a guess for the result based on:
1829 # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
1830 my $len = _len($c,$x_org);
1831 my $log = log($base->[-1]) / log(10);
1832
1833 # for each additional element in $base, we add $BASE_LEN to the result,
1834 # based on the observation that log($BASE,10) is BASE_LEN and
1835 # log(x*y) == log(x) + log(y):
1836 $log += ((scalar @$base)-1) * $BASE_LEN;
1837
1838 # calculate now a guess based on the values obtained above:
1839 my $res = int($len / $log);
1840
1841 $x->[0] = $res;
1842 my $trial = _pow ($c, _copy($c, $base), $x);
1843 my $a = _acmp($c,$trial,$x_org);
1844
1845# print STDERR "# trial ", _str($c,$x)," was: $a (0 = exact, -1 too small, +1 too big)\n";
1846
1847 # found an exact result?
1848 return ($x,1) if $a == 0;
1849
1850 if ($a > 0)
1851 {
1852 # or too big
1853 _div($c,$trial,$base); _dec($c, $x);
1854 while (($a = _acmp($c,$trial,$x_org)) > 0)
1855 {
1856# print STDERR "# big _log_int at ", _str($c,$x), "\n";
1857 _div($c,$trial,$base); _dec($c, $x);
1858 }
1859 # result is now exact (a == 0), or too small (a < 0)
1860 return ($x, $a == 0 ? 1 : 0);
1861 }
1862
1863 # else: result was to small
1864 _mul($c,$trial,$base);
1865
1866 # did we now get the right result?
1867 $a = _acmp($c,$trial,$x_org);
1868
1869 if ($a == 0) # yes, exactly
1870 {
1871 _inc($c, $x);
1872 return ($x,1);
1873 }
1874 return ($x,0) if $a > 0;
1875
1876 # Result still too small (we should come here only if the estimate above
1877 # was very off base):
1878
1879 # Now let the normal trial run obtain the real result
1880 # Simple loop that increments $x by 2 in each step, possible overstepping
1881 # the real result
1882
1883 my $base_mul = _mul($c, _copy($c,$base), $base); # $base * $base
1884
1885 while (($a = _acmp($c,$trial,$x_org)) < 0)
1886 {
1887# print STDERR "# small _log_int at ", _str($c,$x), "\n";
1888 _mul($c,$trial,$base_mul); _add($c, $x, [2]);
1889 }
1890
1891 my $exact = 1;
1892 if ($a > 0)
1893 {
1894 # overstepped the result
1895 _dec($c, $x);
1896 _div($c,$trial,$base);
1897 $a = _acmp($c,$trial,$x_org);
1898 if ($a > 0)
1899 {
1900 _dec($c, $x);
1901 }
1902 $exact = 0 if $a != 0; # a = -1 => not exact result, a = 0 => exact
1903 }
1904
1905 ($x,$exact); # return result
1906 }
1907
1908# for debugging:
19092797µs299µs
# spent 55µs (11+44) within Math::BigInt::Calc::BEGIN@1909 which was called: # once (11µs+44µs) by Math::BigInt::BEGIN@1 at line 1909
use constant DEBUG => 0;
# spent 55µs making 1 call to Math::BigInt::Calc::BEGIN@1909 # spent 44µs making 1 call to constant::import
19101200ns my $steps = 0;
1911 sub steps { $steps };
1912
1913sub _sqrt
1914 {
1915 # square-root of $x in place
1916 # Compute a guess of the result (by rule of thumb), then improve it via
1917 # Newton's method.
1918 my ($c,$x) = @_;
1919
1920 if (scalar @$x == 1)
1921 {
1922 # fits into one Perl scalar, so result can be computed directly
1923 $x->[0] = int(sqrt($x->[0]));
1924 return $x;
1925 }
1926 my $y = _copy($c,$x);
1927 # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1928 # since our guess will "grow"
1929 my $l = int((_len($c,$x)-1) / 2);
1930
1931 my $lastelem = $x->[-1]; # for guess
1932 my $elems = scalar @$x - 1;
1933 # not enough digits, but could have more?
1934 if ((length($lastelem) <= 3) && ($elems > 1))
1935 {
1936 # right-align with zero pad
1937 my $len = length($lastelem) & 1;
1938 print "$lastelem => " if DEBUG;
1939 $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1940 # former odd => make odd again, or former even to even again
1941 $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1942 print "$lastelem\n" if DEBUG;
1943 }
1944
1945 # construct $x (instead of _lsft($c,$x,$l,10)
1946 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1947 $l = int($l / $BASE_LEN);
1948 print "l = $l " if DEBUG;
1949
1950 splice @$x,$l; # keep ref($x), but modify it
1951
1952 # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1953 # that gives us:
1954 # 14400 00000 => sqrt(14400) => guess first digits to be 120
1955 # 144000 000000 => sqrt(144000) => guess 379
1956
1957 print "$lastelem (elems $elems) => " if DEBUG;
1958 $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
1959 my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
1960 $r -= 1 if $elems & 1 == 0; # 70 => 7
1961
1962 # padd with zeros if result is too short
1963 $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1964 print "now ",$x->[-1] if DEBUG;
1965 print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1966
1967 # If @$x > 1, we could compute the second elem of the guess, too, to create
1968 # an even better guess. Not implemented yet. Does it improve performance?
1969 $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
1970
1971 print "start x= ",_str($c,$x),"\n" if DEBUG;
1972 my $two = _two();
1973 my $last = _zero();
1974 my $lastlast = _zero();
1975 $steps = 0 if DEBUG;
1976 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1977 {
1978 $steps++ if DEBUG;
1979 $lastlast = _copy($c,$last);
1980 $last = _copy($c,$x);
1981 _add($c,$x, _div($c,_copy($c,$y),$x));
1982 _div($c,$x, $two );
1983 print " x= ",_str($c,$x),"\n" if DEBUG;
1984 }
1985 print "\nsteps in sqrt: $steps, " if DEBUG;
1986 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
1987 print " final ",$x->[-1],"\n" if DEBUG;
1988 $x;
1989 }
1990
1991sub _root
1992 {
1993 # take n'th root of $x in place (n >= 3)
1994 my ($c,$x,$n) = @_;
1995
1996 if (scalar @$x == 1)
1997 {
1998 if (scalar @$n > 1)
1999 {
2000 # result will always be smaller than 2 so trunc to 1 at once
2001 $x->[0] = 1;
2002 }
2003 else
2004 {
2005 # fits into one Perl scalar, so result can be computed directly
2006 # cannot use int() here, because it rounds wrongly (try
2007 # (81 ** 3) ** (1/3) to see what I mean)
2008 #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
2009 # round to 8 digits, then truncate result to integer
2010 $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
2011 }
2012 return $x;
2013 }
2014
2015 # we know now that X is more than one element long
2016
2017 # if $n is a power of two, we can repeatedly take sqrt($X) and find the
2018 # proper result, because sqrt(sqrt($x)) == root($x,4)
2019 my $b = _as_bin($c,$n);
2020 if ($b =~ /0b1(0+)$/)
2021 {
2022 my $count = CORE::length($1); # 0b100 => len('00') => 2
2023 my $cnt = $count; # counter for loop
2024 unshift (@$x, 0); # add one element, together with one
2025 # more below in the loop this makes 2
2026 while ($cnt-- > 0)
2027 {
2028 # 'inflate' $X by adding one element, basically computing
2029 # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
2030 # since len(sqrt($X)) approx == len($x) / 2.
2031 unshift (@$x, 0);
2032 # calculate sqrt($x), $x is now one element to big, again. In the next
2033 # round we make that two, again.
2034 _sqrt($c,$x);
2035 }
2036 # $x is now one element to big, so truncate result by removing it
2037 splice (@$x,0,1);
2038 }
2039 else
2040 {
2041 # trial computation by starting with 2,4,8,16 etc until we overstep
2042 my $step;
2043 my $trial = _two();
2044
2045 # while still to do more than X steps
2046 do
2047 {
2048 $step = _two();
2049 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2050 {
2051 _mul ($c, $step, [2]);
2052 _add ($c, $trial, $step);
2053 }
2054
2055 # hit exactly?
2056 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
2057 {
2058 @$x = @$trial; # make copy while preserving ref to $x
2059 return $x;
2060 }
2061 # overstepped, so go back on step
2062 _sub($c, $trial, $step);
2063 } while (scalar @$step > 1 || $step->[0] > 128);
2064
2065 # reset step to 2
2066 $step = _two();
2067 # add two, because $trial cannot be exactly the result (otherwise we would
2068 # already have found it)
2069 _add($c, $trial, $step);
2070
2071 # and now add more and more (2,4,6,8,10 etc)
2072 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2073 {
2074 _add ($c, $trial, $step);
2075 }
2076
2077 # hit not exactly? (overstepped)
2078 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2079 {
2080 _dec($c,$trial);
2081 }
2082
2083 # hit not exactly? (overstepped)
2084 # 80 too small, 81 slightly too big, 82 too big
2085 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2086 {
2087 _dec ($c, $trial);
2088 }
2089
2090 @$x = @$trial; # make copy while preserving ref to $x
2091 return $x;
2092 }
2093 $x;
2094 }
2095
2096##############################################################################
2097# binary stuff
2098
2099sub _and
2100 {
2101 my ($c,$x,$y) = @_;
2102
2103 # the shortcut makes equal, large numbers _really_ fast, and makes only a
2104 # very small performance drop for small numbers (e.g. something with less
2105 # than 32 bit) Since we optimize for large numbers, this is enabled.
2106 return $x if _acmp($c,$x,$y) == 0; # shortcut
2107
2108 my $m = _one(); my ($xr,$yr);
2109 my $mask = $AND_MASK;
2110
2111 my $x1 = $x;
2112 my $y1 = _copy($c,$y); # make copy
2113 $x = _zero();
2114 my ($b,$xrr,$yrr);
21152136µs212µs
# spent 10µs (8+2) within Math::BigInt::Calc::BEGIN@2115 which was called: # once (8µs+2µs) by Math::BigInt::BEGIN@1 at line 2115
use integer;
# spent 10µs making 1 call to Math::BigInt::Calc::BEGIN@2115 # spent 2µs making 1 call to integer::import
2116 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2117 {
2118 ($x1, $xr) = _div($c,$x1,$mask);
2119 ($y1, $yr) = _div($c,$y1,$mask);
2120
2121 # make ints() from $xr, $yr
2122 # this is when the AND_BITS are greater than $BASE and is slower for
2123 # small (<256 bits) numbers, but faster for large numbers. Disabled
2124 # due to KISS principle
2125
2126# $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2127# $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2128# _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
2129
2130 # 0+ due to '&' doesn't work in strings
2131 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
2132 _mul($c,$m,$mask);
2133 }
2134 $x;
2135 }
2136
2137sub _xor
2138 {
2139 my ($c,$x,$y) = @_;
2140
2141 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
2142
2143 my $m = _one(); my ($xr,$yr);
2144 my $mask = $XOR_MASK;
2145
2146 my $x1 = $x;
2147 my $y1 = _copy($c,$y); # make copy
2148 $x = _zero();
2149 my ($b,$xrr,$yrr);
21502145µs210µs
# spent 8µs (7+1) within Math::BigInt::Calc::BEGIN@2150 which was called: # once (7µs+1µs) by Math::BigInt::BEGIN@1 at line 2150
use integer;
# spent 8µs making 1 call to Math::BigInt::Calc::BEGIN@2150 # spent 1µs making 1 call to integer::import
2151 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2152 {
2153 ($x1, $xr) = _div($c,$x1,$mask);
2154 ($y1, $yr) = _div($c,$y1,$mask);
2155 # make ints() from $xr, $yr (see _and())
2156 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2157 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2158 #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
2159
2160 # 0+ due to '^' doesn't work in strings
2161 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
2162 _mul($c,$m,$mask);
2163 }
2164 # the loop stops when the shorter of the two numbers is exhausted
2165 # the remainder of the longer one will survive bit-by-bit, so we simple
2166 # multiply-add it in
2167 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2168 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2169
2170 $x;
2171 }
2172
2173sub _or
2174 {
2175 my ($c,$x,$y) = @_;
2176
2177 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
2178
2179 my $m = _one(); my ($xr,$yr);
2180 my $mask = $OR_MASK;
2181
2182 my $x1 = $x;
2183 my $y1 = _copy($c,$y); # make copy
2184 $x = _zero();
2185 my ($b,$xrr,$yrr);
218621.40ms213µs
# spent 11µs (10+1) within Math::BigInt::Calc::BEGIN@2186 which was called: # once (10µs+1µs) by Math::BigInt::BEGIN@1 at line 2186
use integer;
# spent 11µs making 1 call to Math::BigInt::Calc::BEGIN@2186 # spent 2µs making 1 call to integer::import
2187 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2188 {
2189 ($x1, $xr) = _div($c,$x1,$mask);
2190 ($y1, $yr) = _div($c,$y1,$mask);
2191 # make ints() from $xr, $yr (see _and())
2192# $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2193# $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2194# _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
2195
2196 # 0+ due to '|' doesn't work in strings
2197 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
2198 _mul($c,$m,$mask);
2199 }
2200 # the loop stops when the shorter of the two numbers is exhausted
2201 # the remainder of the longer one will survive bit-by-bit, so we simple
2202 # multiply-add it in
2203 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2204 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2205
2206 $x;
2207 }
2208
2209sub _as_hex
2210 {
2211 # convert a decimal number to hex (ref to array, return ref to string)
2212 my ($c,$x) = @_;
2213
2214 # fits into one element (handle also 0x0 case)
2215 return sprintf("0x%x",$x->[0]) if @$x == 1;
2216
2217 my $x1 = _copy($c,$x);
2218
2219 my $es = '';
2220 my ($xr, $h, $x10000);
2221 if ($] >= 5.006)
2222 {
2223 $x10000 = [ 0x10000 ]; $h = 'h4';
2224 }
2225 else
2226 {
2227 $x10000 = [ 0x1000 ]; $h = 'h3';
2228 }
2229 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
2230 {
2231 ($x1, $xr) = _div($c,$x1,$x10000);
2232 $es .= unpack($h,pack('V',$xr->[0]));
2233 }
2234 $es = reverse $es;
2235 $es =~ s/^[0]+//; # strip leading zeros
2236 '0x' . $es; # return result prepended with 0x
2237 }
2238
2239sub _as_bin
2240 {
2241 # convert a decimal number to bin (ref to array, return ref to string)
2242 my ($c,$x) = @_;
2243
2244 # fits into one element (and Perl recent enough), handle also 0b0 case
2245 # handle zero case for older Perls
2246 if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
2247 {
2248 my $t = '0b0'; return $t;
2249 }
2250 if (@$x == 1 && $] >= 5.006)
2251 {
2252 my $t = sprintf("0b%b",$x->[0]);
2253 return $t;
2254 }
2255 my $x1 = _copy($c,$x);
2256
2257 my $es = '';
2258 my ($xr, $b, $x10000);
2259 if ($] >= 5.006)
2260 {
2261 $x10000 = [ 0x10000 ]; $b = 'b16';
2262 }
2263 else
2264 {
2265 $x10000 = [ 0x1000 ]; $b = 'b12';
2266 }
2267 while (!(@$x1 == 1 && $x1->[0] == 0)) # _is_zero()
2268 {
2269 ($x1, $xr) = _div($c,$x1,$x10000);
2270 $es .= unpack($b,pack('v',$xr->[0]));
2271 }
2272 $es = reverse $es;
2273 $es =~ s/^[0]+//; # strip leading zeros
2274 '0b' . $es; # return result prepended with 0b
2275 }
2276
2277sub _as_oct
2278 {
2279 # convert a decimal number to octal (ref to array, return ref to string)
2280 my ($c,$x) = @_;
2281
2282 # fits into one element (handle also 0 case)
2283 return sprintf("0%o",$x->[0]) if @$x == 1;
2284
2285 my $x1 = _copy($c,$x);
2286
2287 my $es = '';
2288 my $xr;
2289 my $x1000 = [ 0100000 ];
2290 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
2291 {
2292 ($x1, $xr) = _div($c,$x1,$x1000);
2293 $es .= reverse sprintf("%05o", $xr->[0]);
2294 }
2295 $es = reverse $es;
2296 $es =~ s/^[0]+//; # strip leading zeros
2297 '0' . $es; # return result prepended with 0
2298 }
2299
2300sub _from_oct
2301 {
2302 # convert a octal number to decimal (string, return ref to array)
2303 my ($c,$os) = @_;
2304
2305 # for older Perls, play safe
2306 my $m = [ 0100000 ];
2307 my $d = 5; # 5 digits at a time
2308
2309 my $mul = _one();
2310 my $x = _zero();
2311
2312 my $len = int( (length($os)-1)/$d ); # $d digit parts, w/o the '0'
2313 my $val; my $i = -$d;
2314 while ($len >= 0)
2315 {
2316 $val = substr($os,$i,$d); # get oct digits
2317 $val = CORE::oct($val);
2318 $i -= $d; $len --;
2319 my $adder = [ $val ];
2320 _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2321 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
2322 }
2323 $x;
2324 }
2325
2326sub _from_hex
2327 {
2328 # convert a hex number to decimal (string, return ref to array)
2329 my ($c,$hs) = @_;
2330
2331 my $m = _new($c, 0x10000000); # 28 bit at a time (<32 bit!)
2332 my $d = 7; # 7 digits at a time
2333 if ($] <= 5.006)
2334 {
2335 # for older Perls, play safe
2336 $m = [ 0x10000 ]; # 16 bit at a time (<32 bit!)
2337 $d = 4; # 4 digits at a time
2338 }
2339
2340 my $mul = _one();
2341 my $x = _zero();
2342
2343 my $len = int( (length($hs)-2)/$d ); # $d digit parts, w/o the '0x'
2344 my $val; my $i = -$d;
2345 while ($len >= 0)
2346 {
2347 $val = substr($hs,$i,$d); # get hex digits
2348 $val =~ s/^0x// if $len == 0; # for last part only because
2349 $val = CORE::hex($val); # hex does not like wrong chars
2350 $i -= $d; $len --;
2351 my $adder = [ $val ];
2352 # if the resulting number was to big to fit into one element, create a
2353 # two-element version (bug found by Mark Lakata - Thanx!)
2354 if (CORE::length($val) > $BASE_LEN)
2355 {
2356 $adder = _new($c,$val);
2357 }
2358 _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2359 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
2360 }
2361 $x;
2362 }
2363
2364sub _from_bin
2365 {
2366 # convert a hex number to decimal (string, return ref to array)
2367 my ($c,$bs) = @_;
2368
2369 # instead of converting X (8) bit at a time, it is faster to "convert" the
2370 # number to hex, and then call _from_hex.
2371
2372 my $hs = $bs;
2373 $hs =~ s/^[+-]?0b//; # remove sign and 0b
2374 my $l = length($hs); # bits
2375 $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0
2376 my $h = '0x' . unpack('H*', pack ('B*', $hs)); # repack as hex
2377
2378 $c->_from_hex($h);
2379 }
2380
2381##############################################################################
2382# special modulus functions
2383
2384sub _modinv
2385 {
2386 # modular multiplicative inverse
2387 my ($c,$x,$y) = @_;
2388
2389 # modulo zero
2390 if (_is_zero($c, $y)) {
2391 return (undef, undef);
2392 }
2393
2394 # modulo one
2395 if (_is_one($c, $y)) {
2396 return (_zero($c), '+');
2397 }
2398
2399 my $u = _zero($c);
2400 my $v = _one($c);
2401 my $a = _copy($c,$y);
2402 my $b = _copy($c,$x);
2403
2404 # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the result
2405 # ($u) at the same time. See comments in BigInt for why this works.
2406 my $q;
2407 my $sign = 1;
2408 {
2409 ($a, $q, $b) = ($b, _div($c, $a, $b)); # step 1
2410 last if _is_zero($c, $b);
2411
2412 my $t = _add($c, # step 2:
2413 _mul($c, _copy($c, $v), $q) , # t = v * q
2414 $u ); # + u
2415 $u = $v; # u = v
2416 $v = $t; # v = t
2417 $sign = -$sign;
2418 redo;
2419 }
2420
2421 # if the gcd is not 1, then return NaN
2422 return (undef, undef) unless _is_one($c, $a);
2423
2424 ($v, $sign == 1 ? '+' : '-');
2425 }
2426
2427sub _modpow
2428 {
2429 # modulus of power ($x ** $y) % $z
2430 my ($c,$num,$exp,$mod) = @_;
2431
2432 # a^b (mod 1) = 0 for all a and b
2433 if (_is_one($c,$mod))
2434 {
2435 @$num = 0;
2436 return $num;
2437 }
2438
2439 # 0^a (mod m) = 0 if m != 0, a != 0
2440 # 0^0 (mod m) = 1 if m != 0
2441 if (_is_zero($c, $num)) {
2442 if (_is_zero($c, $exp)) {
2443 @$num = 1;
2444 } else {
2445 @$num = 0;
2446 }
2447 return $num;
2448 }
2449
2450# $num = _mod($c,$num,$mod); # this does not make it faster
2451
2452 my $acc = _copy($c,$num); my $t = _one();
2453
2454 my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
2455 my $len = length($expbin);
2456 while (--$len >= 0)
2457 {
2458 if ( substr($expbin,$len,1) eq '1') # is_odd
2459 {
2460 _mul($c,$t,$acc);
2461 $t = _mod($c,$t,$mod);
2462 }
2463 _mul($c,$acc,$acc);
2464 $acc = _mod($c,$acc,$mod);
2465 }
2466 @$num = @$t;
2467 $num;
2468 }
2469
2470sub _gcd {
2471 # Greatest common divisor.
2472
2473 my ($c, $x, $y) = @_;
2474
2475 # gcd(0,0) = 0
2476 # gcd(0,a) = a, if a != 0
2477
2478 if (@$x == 1 && $x->[0] == 0) {
2479 if (@$y == 1 && $y->[0] == 0) {
2480 @$x = 0;
2481 } else {
2482 @$x = @$y;
2483 }
2484 return $x;
2485 }
2486
2487 # Until $y is zero ...
2488
2489 until (@$y == 1 && $y->[0] == 0) {
2490
2491 # Compute remainder.
2492
2493 _mod($c, $x, $y);
2494
2495 # Swap $x and $y.
2496
2497 my $tmp = [ @$x ];
2498 @$x = @$y;
2499 $y = $tmp; # no deref here; that would modify input $y
2500 }
2501
2502 return $x;
2503}
2504
2505##############################################################################
2506##############################################################################
2507
250814µs1;
2509__END__
 
# spent 6µs within Math::BigInt::Calc::CORE:match which was called 12 times, avg 533ns/call: # 7 times (4µs+0s) by Math::BigInt::Calc::BEGIN@117 at line 123, avg 500ns/call # 3 times (1µs+0s) by Math::BigInt::Calc::BEGIN@117 at line 141, avg 367ns/call # once (2µs+0s) by Math::BigInt::Calc::BEGIN@117 at line 130 # once (300ns+0s) by Math::BigInt::Calc::BEGIN@117 at line 131
sub Math::BigInt::Calc::CORE:match; # opcode
# spent 50µs within Math::BigInt::Calc::CORE:regcomp which was called 10 times, avg 5µs/call: # 7 times (37µs+0s) by Math::BigInt::Calc::BEGIN@117 at line 123, avg 5µs/call # 3 times (13µs+0s) by Math::BigInt::Calc::BEGIN@117 at line 141, avg 4µs/call
sub Math::BigInt::Calc::CORE:regcomp; # opcode
# spent 400ns within Math::BigInt::Calc::api_version which was called: # once (400ns+0s) by Math::BigInt::import at line 2826 of Math/BigInt.pm
sub Math::BigInt::Calc::api_version; # xsub
Report produced by the NYTProf 5.06 Perl profiler, developed by Tim Bunce and Adam Kaplan.