← Index
NYTProf Performance Profile   « block view • line view • sub view »
For /usr/share/koha/opac/cgi-bin/opac/opac-search.pl
  Run on Tue Oct 15 17:10:45 2013
Reported on Tue Oct 15 17:11:26 2013

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