File Coverage

blib/lib/Math/BigInt/Pari.pm
Criterion Covered Total %
statement 157 162 96.9
branch 60 66 90.9
condition 2 2 100.0
subroutine 57 60 95.0
pod n/a
total 276 290 95.1


line stmt bran cond sub pod time code
1             package Math::BigInt::Pari;
2              
3 8     8   845004 use 5.006002;
  8         115  
4 8     8   56 use strict;
  8         15  
  8         239  
5 8     8   43 use warnings;
  8         12  
  8         337  
6              
7 8     8   6708 use Math::BigInt::Lib 1.999801;
  8         104689  
  8         43  
8              
9             our @ISA = qw< Math::BigInt::Lib >;
10              
11             our $VERSION = '1.3012';
12              
13 8         44 use Math::Pari qw(PARI pari2pv gdivent bittest
14             gcmp gcmp0 gcmp1 gcd ifact gpui gmul
15             binomial lcm
16 8     8   6535 );
  8         95004  
17              
18             # MBI will call this, so catch it and throw it away
19       2     sub import { }
20             sub api_version() { 2; } # we are compatible with MBI v1.83 and up
21              
22             my $zero = PARI(0); # for _copy
23             my $one = PARI(1); # for _inc and _dec
24             my $two = PARI(2); # for _is_two
25             my $ten = PARI(10); # for _digit
26              
27             sub _new {
28             # the . '' is because new($2) will give a magical scalar to us, and PARI
29             # does not like this at all :/
30             # use Devel::Peek; print Dump($_[1]);
31 46940     46940   3417506 PARI($_[1] . '')
32             }
33              
34             sub _from_hex {
35 321     321   4378 my $hex = $_[1];
36 321 50       925 $hex = '0x' . $hex unless substr($hex, 0, 2) eq '0x';
37 321         881 Math::Pari::_hex_cvt($hex);
38             }
39              
40             sub _from_oct {
41 6     6   40 Math::Pari::_hex_cvt('0' . $_[1]);
42             }
43              
44             sub _from_bin {
45 68     68   900 my $bin = $_[1];
46 68 50       207 $bin = '0b' . $bin unless substr($bin, 0, 2) eq '0b';
47 68         197 Math::Pari::_hex_cvt($bin);
48             }
49              
50             sub _as_bin {
51 17     17   661 my $v = unpack('B*', _mp2os($_[1]));
52 17 100       82 return "0b0" if $v eq '';
53 14         65 $v =~ s/^0*/0b/;
54 14         48 $v;
55             }
56              
57             sub _as_oct {
58 15     15   526 my $v = _mp2oct($_[1]);
59 15 100       49 return "00" if $v eq '';
60 12         54 $v =~ s/^0*/0/;
61 12         36 $v;
62             }
63              
64             sub _as_hex {
65 75     75   988 my $v = unpack('H*', _mp2os($_[1]));
66 75 100       213 return "0x0" if $v eq '';
67 62         248 $v =~ s/^0*/0x/;
68 62         187 $v;
69             }
70              
71             sub _to_bin {
72 83     83   1376 my $v = unpack('B*', _mp2os($_[1]));
73 83         416 $v =~ s/^0+//;
74 83 100       203 return "0" if $v eq '';
75 73         220 $v;
76             }
77              
78             sub _to_oct {
79 12     12   511 my $v = _mp2oct($_[1]);
80 12         29 $v =~ s/^0+//;
81 12 100       35 return "0" if $v eq '';
82 10         25 $v;
83             }
84              
85             sub _to_hex {
86 10     10   449 my $v = unpack('H*', _mp2os($_[1]));
87 10         32 $v =~ s/^0+//;
88 10 100       29 return "0" if $v eq '';
89 8         21 $v;
90             }
91              
92             sub _to_bytes {
93 0     0   0 my $bytes = _mp2os($_[1]);
94 0 0       0 return $bytes eq '' ? "\x00" : $bytes;
95             }
96              
97             sub _mp2os {
98 185     185   349 my($p) = @_;
99 185         436 $p = PARI($p);
100 185         1062 my $base = PARI(1) << PARI(4*8);
101 185         742 my $res = '';
102 185         789 while ($p != 0) {
103 184         664 my $r = $p % $base;
104 184         838 $p = ($p - $r) / $base;
105 184         854 my $buf = pack 'V', $r;
106 184 100       597 if ($p == 0) {
107 157 100       803 $buf = $r >= 16777216 ? $buf
    100          
    100          
108             : $r >= 65536 ? substr($buf, 0, 3)
109             : $r >= 256 ? substr($buf, 0, 2)
110             : substr($buf, 0, 1);
111             }
112 184         814 $res .= $buf;
113             }
114 185         1004 scalar reverse $res;
115             }
116              
117             sub _mp2oct {
118 27     27   63 my($p) = @_;
119 27         75 $p = PARI($p);
120 27         65 my $base = PARI(8);
121 27         48 my $res = '';
122 27         124 while ($p != 0) {
123 181         507 my $r = $p % $base;
124 181         712 $p = ($p - $r) / $base;
125 181         891 $res .= $r;
126             }
127 27         120 scalar reverse $res;
128             }
129              
130 12253     12253   1065523 sub _zero { PARI(0) }
131 732     732   52071 sub _one { PARI(1) }
132 103     103   768 sub _two { PARI(2) }
133 2     2   9 sub _ten { PARI(10) }
134              
135 5     5   49 sub _1ex { gpui(PARI(10), $_[1]) }
136              
137 33430     33430   926349 sub _copy { $_[1] + $zero; }
138              
139 52692     52692   1327572 sub _str { pari2pv($_[1]) }
140              
141 10246     10246   93852 sub _num { 0 + pari2pv($_[1]) }
142              
143 23928     23928   328166 sub _add { $_[1] += $_[2] }
144              
145             sub _sub {
146 24176 100   24176   120565 if ($_[3]) {
147 2468         8152 $_[2] = $_[1] - $_[2];
148 2468         5789 return $_[2];
149             }
150 21708         95197 $_[1] -= $_[2];
151             }
152              
153 11434     11434   316918 sub _mul { $_[1] = gmul($_[1], $_[2]) }
154              
155             sub _div {
156 4951 100   4951   53543 if (wantarray) {
157 874         3691 my $r = $_[1] % $_[2];
158 874         4576 $_[1] = gdivent($_[1], $_[2]);
159 874         3552 return ($_[1], $r);
160             }
161 4077         20958 $_[1] = gdivent($_[1], $_[2]);
162             }
163              
164 691     691   7055 sub _mod { $_[1] %= $_[2]; }
165              
166             sub _nok {
167 16     16   127 my ($class, $n, $k) = @_;
168              
169             # Math::Pari doesn't seem to be able to handle the case when n is large and
170             # k is almost as big as n. For instance, the following returns zero (at
171             # least for certain versions and configurations):
172             #
173             # $n = PARI("10000000000000000000");
174             # $k = PARI("9999999999999999999");
175             # print Math::Pari::binomial($n, $k);'
176              
177             # If k > n/2, or, equivalently, 2*k > n, compute nok(n, k) as nok(n, n-k).
178              
179             {
180 16         25 my $twok = $class -> _mul($class -> _two(), $class -> _copy($k));
  16         38  
181 16 100       52 if ($class -> _acmp($twok, $n) > 0) {
182 8         24 $k = $class -> _sub($class -> _copy($n), $k);
183             }
184             }
185              
186 16         221 binomial($n, $k);
187             }
188              
189 2175     2175   22197 sub _inc { ++$_[1]; }
190 85     85   1427 sub _dec { --$_[1]; }
191              
192 36     36   1492 sub _and { $_[1] &= $_[2] }
193              
194 53     53   2049 sub _xor { $_[1] ^= $_[2] }
195              
196 51     51   1794 sub _or { $_[1] |= $_[2] }
197              
198 281     281   7649 sub _pow { gpui($_[1], $_[2]) }
199              
200 25     25   719 sub _gcd { gcd($_[1], $_[2]) }
201              
202 11     11   307 sub _lcm { lcm($_[1], $_[2]) }
203              
204 48205     48205   867615 sub _len { length(pari2pv($_[1])) } # costly!
205              
206             # XXX TODO: calc len in base 2 then appr. in base 10
207 0     0   0 sub _alen { length(pari2pv($_[1])) }
208              
209             sub _zeros {
210 35995     35995   1181334 my ($class, $x) = @_;
211 35995 100       95432 return 0 if gcmp0($x); # 0 has no trailing zeros
212              
213 35380         69809 my $u = $class -> _str($x);
214 35380         431009 $u =~ /(0*)\z/;
215 35380         105051 return length($1);
216              
217             #my $s = pari2pv($_[1]);
218             #my $i = length($s);
219             #my $zeros = 0;
220             #while (--$i >= 0) {
221             # substr($s, $i, 1) eq '0' ? $zeros ++ : last;
222             #}
223             #$zeros;
224             }
225              
226             sub _digit {
227             # if $n < 0, we need to count from left and thus can't use the other method:
228 27 100   27   281 if ($_[2] < 0) {
229 10         165 return substr(pari2pv($_[1]), -1 - $_[2], 1);
230             }
231             # else this is faster (except for very short numbers)
232             # shift the number right by $n digits, then extract last digit via % 10
233 17         370 pari2pv(gdivent($_[1], $ten ** $_[2]) % $ten);
234             }
235              
236 125363     125363   3015748 sub _is_zero { gcmp0($_[1]) }
237              
238 2049     2049   20322 sub _is_one { gcmp1($_[1]) }
239              
240 58 100   58   2174 sub _is_two { gcmp($_[1], $two) ? 0 : 1 }
241              
242 2 100   2   13 sub _is_ten { gcmp($_[1], $ten) ? 0 : 1 }
243              
244 23 100   23   1589 sub _is_even { bittest($_[1], 0) ? 0 : 1 }
245              
246 196 100   196   4993 sub _is_odd { bittest($_[1], 0) ? 1 : 0 }
247              
248             sub _acmp {
249 30013   100 30013   402424 my $i = gcmp($_[1], $_[2]) || 0;
250             # work around bug in Pari (on 64bit systems?)
251 30013 100       60132 $i = -1 if $i == 4294967295;
252 30013         56262 $i;
253             }
254              
255             sub _check {
256 1830     1830   747520 my ($class, $x) = @_;
257 1830 50       5001 return "Undefined" unless defined $x;
258 1830 100       4428 return "$x is not a reference to Math::Pari"
259             unless ref($x) eq 'Math::Pari';
260 1829         4092 return 0;
261             }
262              
263             sub _sqrt {
264             # square root of $x
265 242     242   4517 my ($class, $x) = @_;
266 242         1268 my $y = Math::Pari::sqrtint($x);
267 242 50       2460 return $y * $y > $x ? $y - $one : $y; # bug in sqrtint()?
268             }
269              
270             sub _root {
271             # n'th root
272              
273             # Native version:
274 47     47   1743 return $_[1] = int(Math::Pari::sqrtn($_[1] + 0.5, $_[2]));
275             }
276              
277             sub _modpow {
278             # modulus of power ($x ** $y) % $z
279 144     144   1383 my ($c, $num, $exp, $mod) = @_;
280              
281             # in the trivial case,
282 144 100       425 if (gcmp1($mod)) {
283 50         132 $num = PARI(0);
284 50         119 return $num;
285             }
286              
287 94 100       241 if (gcmp1($num)) {
288 29         68 $num = PARI(1);
289 29         72 return $num;
290             }
291              
292 65 100       152 if (gcmp0($num)) {
293 12 100       27 if (gcmp0($exp)) {
294 3         15 return PARI(1);
295             } else {
296 9         37 return PARI(0);
297             }
298             }
299              
300 53         119 my $acc = $c -> _copy($num);
301 53         122 my $t = $c -> _one();
302              
303 53         119 my $expbin = $c -> _to_bin($exp);
304 53         96 my $len = length($expbin);
305 53         114 while (--$len >= 0) {
306 268 100       563 if (substr($expbin, $len, 1) eq '1') { # is_odd
307 141         346 $c -> _mul($t, $acc);
308 141         262 $t = $c -> _mod($t, $mod);
309             }
310 268         580 $c -> _mul($acc, $acc);
311 268         500 $acc = $c -> _mod($acc, $mod);
312             }
313 53         83 $num = $t;
314 53         172 $num;
315             }
316              
317             sub _rsft {
318             # (X, Y, N) = @_; means X >> Y in base N
319              
320 14272 100   14272   66826 if ($_[3] != 2) {
321 14257         136792 return $_[1] = gdivent($_[1], PARI($_[3]) ** $_[2]);
322             }
323 15         158 $_[1] >>= $_[2];
324             }
325              
326             sub _lsft {
327             # (X, Y, N) = @_; means X >> Y in base N
328              
329 9944 100   9944   48161 if ($_[3] != 2) {
330 9936         97733 return $_[1] *= PARI($_[3]) ** $_[2];
331             }
332 8         95 $_[1] <<= $_[2];
333             }
334              
335             sub _fac {
336             # factorial of argument
337 48     48   1805 $_[1] = ifact($_[1]);
338             }
339              
340             # _set() - set an already existing object to the given scalar value
341              
342             sub _set {
343 0     0     my ($c, $x, $y) = @_;
344 0           *x = \PARI($y . '');
345             }
346              
347             1;
348              
349             __END__