1 #include "EXTERN.h"
2 #include "perl.h"
3 #include "XSUB.h"
4 #include "gmp.h"
5
6 /*
7 Math::GMP, a Perl module for high-speed arbitrary size integer
8 calculations
9 Copyright (C) 2000 James H. Turner
10
11 This library is free software; you can redistribute it and/or
12 modify it under the terms of the GNU Library General Public
13 License as published by the Free Software Foundation; either
14 version 2 of the License, or (at your option) any later version.
15
16 This library is distributed in the hope that it will be useful,
17 but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 Library General Public License for more details.
20
21 You should have received a copy of the GNU Library General Public
22 License along with this library; if not, write to the Free
23 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24
25 You can contact the author at chip@redhat.com, chipt@cpan.org, or by mail:
26
27 Chip Turner
28 Red Hat Inc.
29 2600 Meridian Park Blvd
30 Durham, NC 27713
31 */
32
33 #define SWAP_GMP if (swap) { mpz_t* t = m; m = n; n = t; }
34
35 /*
36 * mpz_rootrem() has bug with negative first argument before 5.1.0
37 */
need_rootrem_workaround(mpz_t * m,unsigned long n)38 static int need_rootrem_workaround(mpz_t* m, unsigned long n) {
39 /* workaround only valid for n odd (n even should be an error) */
40 if ((n & 1) == 0)
41 return 0;
42
43 /* workaround only relevant for m negative */
44 if (mpz_sgn(*m) >= 0)
45 return 0;
46
47 /* workaround only needed for gmp_version < 5.1.0 */
48 if ((gmp_version[0] && gmp_version[1] != '.') /* >= 10.0.0 */
49 || (gmp_version[0] > '5') /* >= 6.0.0 */
50 || (gmp_version[0] == '5' && gmp_version[2] != '0') /* >= 5.1.0 */
51 )
52 return 0;
53
54 return 1;
55 }
56
57 static int
not_here(char * s)58 not_here(char *s)
59 {
60 croak("%s not implemented on this architecture", s);
61 return -1;
62 }
63
64 #if 0
65 static double
66 constant(char *name, int arg)
67 {
68 errno = 0;
69 switch (*name) {
70 }
71 errno = EINVAL;
72 return 0;
73
74 not_there:
75 errno = ENOENT;
76 return 0;
77 }
78 #endif
79
80 mpz_t *
pv2gmp(char * pv)81 pv2gmp(char* pv)
82 {
83 mpz_t* z;
84 SV* sv;
85
86 z = malloc (sizeof(mpz_t));
87 mpz_init_set_str(*z, pv, 0);
88 sv = sv_newmortal();
89 sv_setref_pv(sv, "Math::GMP", (void*)z);
90 return z;
91 }
92
93 mpz_t *
sv2gmp(SV * sv)94 sv2gmp(SV* sv)
95 {
96 char* pv;
97
98 /* MAYCHANGE in perlguts.pod - bug in perl */
99 if (SvGMAGICAL(sv)) mg_get(sv);
100
101 if (SvROK(sv) && sv_derived_from(sv, "Math::GMP")) {
102 IV tmp = SvIV((SV*)SvRV(sv));
103 return (mpz_t *)tmp;
104 }
105
106 pv = SvPV_nolen(sv);
107 return pv2gmp(pv);
108 }
109
110
111 MODULE = Math::GMP PACKAGE = Math::GMP
112 PROTOTYPES: ENABLE
113
114 mpz_t *
115 new_from_scalar(s)
116 char * s
117
118 CODE:
119 RETVAL = malloc (sizeof(mpz_t));
120 mpz_init_set_str(*RETVAL, s, 0);
121 OUTPUT:
122 RETVAL
123
124 mpz_t *
125 new_from_scalar_with_base(s, b)
126 char * s
127 int b
128
129 CODE:
130 RETVAL = malloc (sizeof(mpz_t));
131 mpz_init_set_str(*RETVAL, s, b);
132 OUTPUT:
133 RETVAL
134
135 void
136 destroy(n)
137 mpz_t *n
138
139 CODE:
140 mpz_clear(*n);
141 free(n);
142
143 SV *
144 _gmp_build_version()
145 CODE:
146 char buf[] = STRINGIFY(__GNU_MP_VERSION)
147 "." STRINGIFY(__GNU_MP_VERSION_MINOR)
148 "." STRINGIFY(__GNU_MP_VERSION_PATCHLEVEL);
149 RETVAL = newSV(6);
150 scan_vstring(buf, buf + sizeof(buf), RETVAL);
151 OUTPUT:
152 RETVAL
153
154 SV *
155 _gmp_lib_version()
156 CODE:
157 const char* v = gmp_version;
158 RETVAL = newSV(6);
159 scan_vstring(v, v + strlen(v), RETVAL);
160 OUTPUT:
161 RETVAL
162
163 SV *
164 stringify(n)
165 mpz_t * n
166
167 PREINIT:
168 int len;
169
170 CODE:
171 len = mpz_sizeinbase(*n, 10);
172 {
173 char *buf;
174 buf = malloc(len + 2);
175 mpz_get_str(buf, 10, *n);
176 RETVAL = newSVpv(buf, strlen(buf));
177 free(buf);
178 }
179 OUTPUT:
180 RETVAL
181
182
183 SV *
get_str_gmp(n,b)184 get_str_gmp(n, b)
185 mpz_t * n
186 int b
187
188 PREINIT:
189 int len;
190
191 CODE:
192 len = mpz_sizeinbase(*n, b);
193 {
194 char *buf;
195 buf = malloc(len + 2);
196 mpz_get_str(buf, b, *n);
197 RETVAL = newSVpv(buf, strlen(buf));
198 free(buf);
199 }
200 OUTPUT:
201 RETVAL
202
203 int
204 sizeinbase_gmp(n, b)
205 mpz_t * n
206 int b
207
208 CODE:
209 RETVAL = mpz_sizeinbase(*n, b);
210 OUTPUT:
211 RETVAL
212
213 unsigned long
214 uintify(n)
215 mpz_t * n
216
217 CODE:
218 RETVAL = mpz_get_ui(*n);
219 OUTPUT:
220 RETVAL
221
222 void
223 add_ui_gmp(n, v)
224 mpz_t * n
225 unsigned long v
226
227 CODE:
228 mpz_add_ui(*n, *n, v);
229
230
231 long
232 intify(n)
233 mpz_t * n
234
235 CODE:
236 RETVAL = mpz_get_si(*n);
237 OUTPUT:
238 RETVAL
239
240 mpz_t *
241 mul_2exp_gmp(n, e)
242 mpz_t * n
243 unsigned long e
244
245 CODE:
246 RETVAL = malloc (sizeof(mpz_t));
247 mpz_init(*RETVAL);
248 mpz_mul_2exp(*RETVAL, *n, e);
249 OUTPUT:
250 RETVAL
251
252 mpz_t *
253 div_2exp_gmp(n, e)
254 mpz_t * n
255 unsigned long e
256
257 CODE:
258 RETVAL = malloc (sizeof(mpz_t));
259 mpz_init(*RETVAL);
260 mpz_div_2exp(*RETVAL, *n, e);
261 OUTPUT:
262 RETVAL
263
264
265 mpz_t *
266 powm_gmp(n, exp, mod)
267 mpz_t * n
268 mpz_t * exp
269 mpz_t * mod
270
271 CODE:
272 RETVAL = malloc (sizeof(mpz_t));
273 mpz_init(*RETVAL);
274 mpz_powm(*RETVAL, *n, *exp, *mod);
275 OUTPUT:
276 RETVAL
277
278
279 mpz_t *
280 mmod_gmp(a, b)
281 mpz_t * a
282 mpz_t * b
283
284 CODE:
285 RETVAL = malloc (sizeof(mpz_t));
286 mpz_init(*RETVAL);
287 mpz_mmod(*RETVAL, *a, *b);
288 OUTPUT:
289 RETVAL
290
291
292 mpz_t *
293 mod_2exp_gmp(in, cnt)
294 mpz_t * in
295 unsigned long cnt
296
297 CODE:
298 RETVAL = malloc (sizeof(mpz_t));
299 mpz_init(*RETVAL);
300 mpz_mod_2exp(*RETVAL, *in, cnt);
301 OUTPUT:
302 RETVAL
303
304
305 mpz_t *
306 op_add(m,n,swap)
307 mpz_t * m
308 mpz_t * n
309 bool swap
310
311 CODE:
312 RETVAL = malloc (sizeof(mpz_t));
313 mpz_init(*RETVAL);
314 mpz_add(*RETVAL, *m, *n);
315 OUTPUT:
316 RETVAL
317
318
319 mpz_t *
320 op_sub(m,n,swap)
321 mpz_t * m
322 mpz_t * n
323 bool swap
324
325 CODE:
326 SWAP_GMP
327 RETVAL = malloc (sizeof(mpz_t));
328 mpz_init(*RETVAL);
329 mpz_sub(*RETVAL, *m, *n);
330 OUTPUT:
331 RETVAL
332
333
334 mpz_t *
335 op_mul(m,n,swap)
336 mpz_t * m
337 mpz_t * n
338 bool swap
339
340 CODE:
341 RETVAL = malloc (sizeof(mpz_t));
342 mpz_init(*RETVAL);
343 mpz_mul(*RETVAL, *m, *n);
344 OUTPUT:
345 RETVAL
346
347
348 mpz_t *
349 op_div(m,n,swap)
350 mpz_t * m
351 mpz_t * n
352 bool swap
353
354 CODE:
355 SWAP_GMP
356 RETVAL = malloc (sizeof(mpz_t));
357 mpz_init(*RETVAL);
358 mpz_div(*RETVAL, *m, *n);
359 OUTPUT:
360 RETVAL
361
362
363 void
364 bdiv(m,n)
365 mpz_t * m
366 mpz_t * n
367
368 PREINIT:
369 mpz_t * quo;
370 mpz_t * rem;
371 PPCODE:
372 quo = malloc (sizeof(mpz_t));
373 rem = malloc (sizeof(mpz_t));
374 mpz_init(*quo);
375 mpz_init(*rem);
376 mpz_tdiv_qr(*quo, *rem, *m, *n);
377 EXTEND(SP, 2);
378 PUSHs(sv_setref_pv(sv_newmortal(), "Math::GMP", (void*)quo));
379 PUSHs(sv_setref_pv(sv_newmortal(), "Math::GMP", (void*)rem));
380
381
382
383 mpz_t *
384 op_mod(m,n,swap)
385 mpz_t * m
386 mpz_t * n
387 bool swap
388
389 CODE:
390 SWAP_GMP
391 RETVAL = malloc (sizeof(mpz_t));
392 mpz_init(*RETVAL);
393 mpz_mod(*RETVAL, *m, *n);
394 OUTPUT:
395 RETVAL
396
397 mpz_t *
398 bmodinv(m,n)
399 mpz_t * m
400 mpz_t * n
401
402 CODE:
403 RETVAL = malloc (sizeof(mpz_t));
404 mpz_init(*RETVAL);
405 mpz_invert(*RETVAL, *m, *n);
406 OUTPUT:
407 RETVAL
408
409
410 int
op_spaceship(m,n,swap)411 op_spaceship(m,n,swap)
412 mpz_t * m
413 mpz_t * n
414 bool swap
415
416 PREINIT:
417 int i;
418 CODE:
419 i = mpz_cmp(*m, *n);
420 if (swap) {
421 i = -i;
422 }
423 RETVAL = (i < 0) ? -1 : (i > 0) ? 1 : 0;
424 OUTPUT:
425 RETVAL
426
427 int
428 op_eq(m,n,swap)
429 mpz_t* m
430 mpz_t* n
431 bool swap
432
433 PREINIT:
434 int i;
435 CODE:
436 i = mpz_cmp(*m, *n);
437 RETVAL = (i == 0) ? 1 : 0;
438 OUTPUT:
439 RETVAL
440
441 int
442 legendre(m, n)
443 mpz_t * m
444 mpz_t * n
445
446 CODE:
447 RETVAL = mpz_legendre(*m, *n);
448 OUTPUT:
449 RETVAL
450
451 int
452 jacobi(m, n)
453 mpz_t * m
454 mpz_t * n
455
456 CODE:
457 RETVAL = mpz_jacobi(*m, *n);
458 OUTPUT:
459 RETVAL
460
461 int
462 probab_prime(m, reps)
463 mpz_t * m
464 int reps
465
466 CODE:
467 RETVAL = mpz_probab_prime_p(*m, reps);
468 OUTPUT:
469 RETVAL
470
471 mpz_t *
472 op_pow(m,n)
473 mpz_t * m
474 long n
475
476 CODE:
477 RETVAL = malloc (sizeof(mpz_t));
478 mpz_init(*RETVAL);
479 /* fprintf(stderr, "n is %ld\n", n);*/
480 mpz_pow_ui(*RETVAL, *m, n);
481 OUTPUT:
482 RETVAL
483
484
485 mpz_t *
486 bgcd(m,n)
487 mpz_t * m
488 mpz_t * n
489
490 CODE:
491 RETVAL = malloc (sizeof(mpz_t));
492 mpz_init(*RETVAL);
493 mpz_gcd(*RETVAL, *m, *n);
494 OUTPUT:
495 RETVAL
496
497
498 mpz_t *
499 blcm(m,n)
500 mpz_t * m
501 mpz_t * n
502
503 CODE:
504 RETVAL = malloc (sizeof(mpz_t));
505 mpz_init(*RETVAL);
506 mpz_lcm(*RETVAL, *m, *n);
507 OUTPUT:
508 RETVAL
509
510
511 mpz_t *
512 fibonacci(n)
513 long n
514
515 CODE:
516 RETVAL = malloc (sizeof(mpz_t));
517 mpz_init(*RETVAL);
518 mpz_fib_ui(*RETVAL, n);
519 OUTPUT:
520 RETVAL
521
522
523 mpz_t *
524 band(m,n, ...)
525 mpz_t * m
526 mpz_t * n
527
528 CODE:
529 RETVAL = malloc (sizeof(mpz_t));
530 mpz_init(*RETVAL);
531 mpz_and(*RETVAL, *m, *n);
532 OUTPUT:
533 RETVAL
534
535 mpz_t *
536 bxor(m,n, ...)
537 mpz_t * m
538 mpz_t * n
539
540 CODE:
541 RETVAL = malloc (sizeof(mpz_t));
542 mpz_init(*RETVAL);
543 mpz_xor(*RETVAL, *m, *n);
544 OUTPUT:
545 RETVAL
546
547 mpz_t *
548 bior(m,n, ...)
549 mpz_t * m
550 mpz_t * n
551
552 CODE:
553 RETVAL = malloc (sizeof(mpz_t));
554 mpz_init(*RETVAL);
555 mpz_ior(*RETVAL, *m, *n);
556 OUTPUT:
557 RETVAL
558
559 mpz_t *
560 blshift(m,n,swap)
561 mpz_t * m
562 mpz_t * n
563 bool swap
564
565 CODE:
566 SWAP_GMP
567 RETVAL = malloc (sizeof(mpz_t));
568 mpz_init(*RETVAL);
569 mpz_mul_2exp(*RETVAL, *m, mpz_get_ui(*n));
570 OUTPUT:
571 RETVAL
572
573 mpz_t *
574 brshift(m,n,swap)
575 mpz_t * m
576 mpz_t * n
577 bool swap
578
579 CODE:
580 SWAP_GMP
581 RETVAL = malloc (sizeof(mpz_t));
582 mpz_init(*RETVAL);
583 mpz_div_2exp(*RETVAL, *m, mpz_get_ui(*n));
584 OUTPUT:
585 RETVAL
586
587 mpz_t *
588 bfac(n)
589 long n
590
591 CODE:
592 RETVAL = malloc (sizeof(mpz_t));
593 mpz_init(*RETVAL);
594 mpz_fac_ui(*RETVAL, n);
595 OUTPUT:
596 RETVAL
597
598
599 mpz_t *
600 gmp_copy(m)
601 mpz_t * m
602
603 CODE:
604 RETVAL = malloc (sizeof(mpz_t));
605 mpz_init_set(*RETVAL, *m);
606 OUTPUT:
607 RETVAL
608
609 int
610 gmp_tstbit(m,n)
611 mpz_t * m
612 long n
613
614 CODE:
615 RETVAL = mpz_tstbit(*m,n);
616 OUTPUT:
617 RETVAL
618
619 mpz_t *
620 broot(m,n)
621 mpz_t * m
622 unsigned long n
623
624 CODE:
625 RETVAL = malloc (sizeof(mpz_t));
626 mpz_init(*RETVAL);
627 mpz_root(*RETVAL, *m, n);
628 OUTPUT:
629 RETVAL
630
631 void
632 brootrem(m,n)
633 mpz_t * m
634 unsigned long n
635
636 PREINIT:
637 mpz_t * root;
638 mpz_t * remainder;
639 PPCODE:
640 root = malloc (sizeof(mpz_t));
641 remainder = malloc (sizeof(mpz_t));
642 mpz_init(*root);
643 mpz_init(*remainder);
644 if (need_rootrem_workaround(m, n)) {
645 /* Older libgmp have bugs for negative m, but if we need to we can
646 * work on abs(m) then negate the results.
647 */
648 mpz_neg(*root, *m);
649 mpz_rootrem(*root, *remainder, *root, n);
650 mpz_neg(*root, *root);
651 mpz_neg(*remainder, *remainder);
652 } else {
653 mpz_rootrem(*root, *remainder, *m, n);
654 }
655 EXTEND(SP, 2);
656 PUSHs(sv_setref_pv(sv_newmortal(), "Math::GMP", (void*)root));
657 PUSHs(sv_setref_pv(sv_newmortal(), "Math::GMP", (void*)remainder));
658
659 mpz_t *
660 bsqrt(m)
661 mpz_t * m
662
663 CODE:
664 RETVAL = malloc (sizeof(mpz_t));
665 mpz_init(*RETVAL);
666 mpz_sqrt(*RETVAL, *m);
667 OUTPUT:
668 RETVAL
669
670 void
671 bsqrtrem(m)
672 mpz_t * m
673
674 PREINIT:
675 mpz_t * sqrt;
676 mpz_t * remainder;
677 PPCODE:
678 sqrt = malloc (sizeof(mpz_t));
679 remainder = malloc (sizeof(mpz_t));
680 mpz_init(*sqrt);
681 mpz_init(*remainder);
682 mpz_sqrtrem(*sqrt, *remainder, *m);
683 EXTEND(SP, 2);
684 PUSHs(sv_setref_pv(sv_newmortal(), "Math::GMP", (void*)sqrt));
685 PUSHs(sv_setref_pv(sv_newmortal(), "Math::GMP", (void*)remainder));
686
687 int
688 is_perfect_power(m)
689 mpz_t * m
690
691 CODE:
692 RETVAL = mpz_perfect_power_p(*m) ? 1 : 0;
693 OUTPUT:
694 RETVAL
695
696 int
697 is_perfect_square(m)
698 mpz_t * m
699
700 CODE:
701 RETVAL = mpz_perfect_square_p(*m) ? 1 : 0;
702 OUTPUT:
703 RETVAL
704
705