1 /* dumbmp mini GMP compatible library.
3 Copyright 2001, 2002, 2004 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of the GNU Lesser General Public License as published by
9 the Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
12 The GNU MP Library is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 License for more details.
17 You should have received a copy of the GNU Lesser General Public License
18 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
21 /* The code here implements a subset (a very limited subset) of the main GMP
22 functions. It's designed for use in a few build-time calculations and
23 will be slow, but highly portable.
25 None of the normal GMP configure things are used, nor any of the normal
26 gmp.h or gmp-impl.h. To use this file in a program just #include
29 ANSI function definitions can be used here, since ansi2knr is run if
30 necessary. But other ANSI-isms like "const" should be avoided.
32 mp_limb_t here is an unsigned long, since that's a sensible type
33 everywhere we know of, with 8*sizeof(unsigned long) giving the number of
34 bits in the type (that not being true for instance with int or short on
37 Only the low half of each mp_limb_t is used, so as to make carry handling
38 and limb multiplies easy. GMP_LIMB_BITS is the number of bits used. */
45 typedef unsigned long mp_limb_t
;
53 #define GMP_LIMB_BITS (sizeof (mp_limb_t) * 8 / 2)
55 #define ABS(x) ((x) >= 0 ? (x) : -(x))
56 #define MIN(l,o) ((l) < (o) ? (l) : (o))
57 #define MAX(h,i) ((h) > (i) ? (h) : (i))
59 #define ALLOC(x) ((x)->_mp_alloc)
60 #define PTR(x) ((x)->_mp_d)
61 #define SIZ(x) ((x)->_mp_size)
62 #define ABSIZ(x) ABS (SIZ (x))
63 #define LOMASK ((1L << GMP_LIMB_BITS) - 1)
64 #define LO(x) ((x) & LOMASK)
65 #define HI(x) ((x) >> GMP_LIMB_BITS)
67 #define ASSERT(cond) \
71 fprintf (stderr, "Assertion failure\n"); \
84 fprintf (stderr
, "Out of memory (alloc %d bytes)\n", n
);
93 return (mp_limb_t
*) xmalloc (n
* sizeof (mp_limb_t
));
97 mem_copyi (char *dst
, char *src
, int size
)
100 for (i
= 0; i
< size
; i
++)
110 for (i
= 2; i
< n
; i
++)
128 mpz_realloc (mpz_t r
, int n
)
134 PTR(r
) = (mp_limb_t
*) realloc (PTR(r
), n
* sizeof (mp_limb_t
));
137 fprintf (stderr
, "Out of memory (realloc to %d)\n", n
);
143 mpn_normalize (mp_limb_t
*rp
, int *rnp
)
146 while (rn
> 0 && rp
[rn
-1] == 0)
152 mpn_copyi (mp_limb_t
*dst
, mp_limb_t
*src
, int n
)
155 for (i
= 0; i
< n
; i
++)
160 mpn_zero (mp_limb_t
*rp
, int rn
)
163 for (i
= 0; i
< rn
; i
++)
171 PTR(r
) = xmalloc_limbs (ALLOC(r
));
181 SIZ (r
) = 0xbadcafeL
;
182 PTR (r
) = (mp_limb_t
*) 0xdeadbeefL
;
188 return (SIZ(a
) > 0 ? 1 : SIZ(a
) == 0 ? 0 : -1);
197 return (PTR(a
)[0] & 1) != 0;
206 return (PTR(a
)[0] & 1) == 0;
210 mpz_sizeinbase (mpz_t a
, int base
)
213 mp_limb_t
*ap
= PTR (a
);
224 for (hi
= ap
[an
- 1]; hi
!= 0; hi
>>= 1)
226 return (an
- 1) * GMP_LIMB_BITS
+ cnt
;
230 mpz_set (mpz_t r
, mpz_t a
)
232 mpz_realloc (r
, ABSIZ (a
));
234 mpn_copyi (PTR(r
), PTR(a
), ABSIZ (a
));
238 mpz_init_set (mpz_t r
, mpz_t a
)
245 mpz_set_ui (mpz_t r
, unsigned long ui
)
252 mpn_normalize (PTR(r
), &rn
);
257 mpz_init_set_ui (mpz_t r
, unsigned long ui
)
264 mpz_setbit (mpz_t r
, unsigned long bit
)
266 int limb
, rn
, extend
;
271 abort (); /* only r>=0 */
273 limb
= bit
/ GMP_LIMB_BITS
;
274 bit
%= GMP_LIMB_BITS
;
276 mpz_realloc (r
, limb
+1);
278 extend
= (limb
+1) - rn
;
280 mpn_zero (rp
+ rn
, extend
);
282 rp
[limb
] |= (mp_limb_t
) 1 << bit
;
283 SIZ(r
) = MAX (rn
, limb
+1);
287 mpz_tstbit (mpz_t r
, unsigned long bit
)
292 abort (); /* only r>=0 */
294 limb
= bit
/ GMP_LIMB_BITS
;
298 bit
%= GMP_LIMB_BITS
;
299 return (PTR(r
)[limb
] >> bit
) & 1;
303 popc_limb (mp_limb_t a
)
315 mpz_popcount (mpz_t a
)
324 for (i
= 0; i
< SIZ(a
); i
++)
325 ret
+= popc_limb (PTR(a
)[i
]);
330 mpz_add (mpz_t r
, mpz_t a
, mpz_t b
)
332 int an
= ABSIZ (a
), bn
= ABSIZ (b
), rn
;
333 mp_limb_t
*rp
, *ap
, *bp
;
337 if ((SIZ (a
) ^ SIZ (b
)) < 0)
338 abort (); /* really subtraction */
342 mpz_realloc (r
, MAX (an
, bn
) + 1);
343 ap
= PTR (a
); bp
= PTR (b
); rp
= PTR (r
);
346 mp_limb_t
*tp
; int tn
;
347 tn
= an
; an
= bn
; bn
= tn
;
348 tp
= ap
; ap
= bp
; bp
= tp
;
352 for (i
= 0; i
< bn
; i
++)
354 t
= ap
[i
] + bp
[i
] + cy
;
358 for (i
= bn
; i
< an
; i
++)
367 mpn_normalize (rp
, &rn
);
372 mpz_add_ui (mpz_t r
, mpz_t a
, unsigned long int ui
)
383 mpz_sub (mpz_t r
, mpz_t a
, mpz_t b
)
385 int an
= ABSIZ (a
), bn
= ABSIZ (b
), rn
;
386 mp_limb_t
*rp
, *ap
, *bp
;
390 if ((SIZ (a
) ^ SIZ (b
)) < 0)
391 abort (); /* really addition */
395 mpz_realloc (r
, MAX (an
, bn
) + 1);
396 ap
= PTR (a
); bp
= PTR (b
); rp
= PTR (r
);
399 mp_limb_t
*tp
; int tn
;
400 tn
= an
; an
= bn
; bn
= tn
;
401 tp
= ap
; ap
= bp
; bp
= tp
;
405 for (i
= 0; i
< bn
; i
++)
407 t
= ap
[i
] - bp
[i
] - cy
;
411 for (i
= bn
; i
< an
; i
++)
423 for (i
= 0; i
< rn
; i
++)
433 mpn_normalize (rp
, &rn
);
438 mpz_sub_ui (mpz_t r
, mpz_t a
, unsigned long int ui
)
449 mpz_mul (mpz_t r
, mpz_t a
, mpz_t b
)
451 int an
= ABSIZ (a
), bn
= ABSIZ (b
), rn
;
452 mp_limb_t
*scratch
, *tmp
, *ap
= PTR (a
), *bp
= PTR (b
);
456 scratch
= xmalloc_limbs (an
+ bn
);
459 for (bi
= 0; bi
< bn
; bi
++)
462 for (ai
= 0; ai
< an
; ai
++)
466 for (bi
= 0; bi
< bn
; bi
++)
468 t
= ap
[ai
] * bp
[bi
] + tmp
[bi
] + cy
;
476 mpn_normalize (scratch
, &rn
);
479 SIZ (r
) = (SIZ (a
) ^ SIZ (b
)) >= 0 ? rn
: -rn
;
484 mpz_mul_ui (mpz_t r
, mpz_t a
, unsigned long int ui
)
495 mpz_mul_2exp (mpz_t r
, mpz_t a
, unsigned long int bcnt
)
506 mpz_ui_pow_ui (mpz_t r
, unsigned long b
, unsigned long e
)
515 for (i
= 0; i
< e
; i
++)
522 mpz_tdiv_q_2exp (mpz_t r
, mpz_t a
, unsigned long int bcnt
)
527 mp_limb_t high_limb
, low_limb
;
531 lcnt
= bcnt
/ GMP_LIMB_BITS
;
532 rn
= ABS (as
) - lcnt
;
544 cnt
= bcnt
% GMP_LIMB_BITS
;
548 tnc
= GMP_LIMB_BITS
- cnt
;
550 low_limb
= high_limb
>> cnt
;
552 for (i
= rn
- 1; i
!= 0; i
--)
555 *rp
++ = low_limb
| LO (high_limb
<< tnc
);
556 low_limb
= high_limb
>> cnt
;
564 mpn_copyi (rp
, ap
, rn
);
567 SIZ (r
) = as
>= 0 ? rn
: -rn
;
572 mpz_tdiv_r_2exp (mpz_t r
, mpz_t a
, unsigned long int bcnt
)
579 bwhole
= bcnt
/ GMP_LIMB_BITS
;
580 bcnt
%= GMP_LIMB_BITS
;
584 PTR(r
)[rn
-1] &= ((mp_limb_t
) 1 << bcnt
) - 1;
585 mpn_normalize (PTR(r
), &rn
);
586 SIZ(r
) = (SIZ(r
) >= 0 ? rn
: -rn
);
591 mpz_cmp (mpz_t a
, mpz_t b
)
593 mp_limb_t
*ap
, *bp
, al
, bl
;
594 int as
= SIZ (a
), bs
= SIZ (b
);
599 return as
> bs
? 1 : -1;
601 sign
= as
> 0 ? 1 : -1;
605 for (i
= ABS (as
) - 1; i
>= 0; i
--)
610 return al
> bl
? sign
: -sign
;
616 mpz_cmp_ui (mpz_t a
, unsigned long b
)
620 mpz_init_set_ui (bz
, b
);
621 ret
= mpz_cmp (a
, bz
);
627 mpz_tdiv_qr (mpz_t q
, mpz_t r
, mpz_t a
, mpz_t b
)
632 ASSERT (mpz_sgn(a
) >= 0);
633 ASSERT (mpz_sgn(b
) > 0);
635 mpz_init_set (tmpr
, a
);
636 mpz_init_set (tmpb
, b
);
639 if (mpz_cmp (tmpr
, tmpb
) > 0)
641 cnt
= mpz_sizeinbase (tmpr
, 2) - mpz_sizeinbase (tmpb
, 2) + 1;
642 mpz_mul_2exp (tmpb
, tmpb
, cnt
);
644 for ( ; cnt
> 0; cnt
--)
646 mpz_mul_2exp (q
, q
, 1);
647 mpz_tdiv_q_2exp (tmpb
, tmpb
, 1L);
648 if (mpz_cmp (tmpr
, tmpb
) >= 0)
650 mpz_sub (tmpr
, tmpr
, tmpb
);
651 mpz_add_ui (q
, q
, 1L);
652 ASSERT (mpz_cmp (tmpr
, tmpb
) < 0);
663 mpz_tdiv_qr_ui (mpz_t q
, mpz_t r
, mpz_t a
, unsigned long b
)
666 mpz_init_set_ui (bz
, b
);
667 mpz_tdiv_qr (q
, r
, a
, bz
);
672 mpz_tdiv_q (mpz_t q
, mpz_t a
, mpz_t b
)
677 mpz_tdiv_qr (q
, r
, a
, b
);
682 mpz_tdiv_q_ui (mpz_t q
, mpz_t n
, unsigned long d
)
685 mpz_init_set_ui (dz
, d
);
686 mpz_tdiv_q (q
, n
, dz
);
690 /* Set inv to the inverse of d, in the style of invert_limb, ie. for
691 udiv_qrnnd_preinv. */
693 mpz_preinv_invert (mpz_t inv
, mpz_t d
, int numb_bits
)
699 norm
= numb_bits
- mpz_sizeinbase (d
, 2);
701 mpz_init_set_ui (t
, 1L);
702 mpz_mul_2exp (t
, t
, 2*numb_bits
- norm
);
703 mpz_tdiv_q (inv
, t
, d
);
705 mpz_mul_2exp (t
, t
, numb_bits
);
706 mpz_sub (inv
, inv
, t
);
711 /* Remove leading '0' characters from the start of a string, by copying the
714 strstrip_leading_zeros (char *s
)
731 mpz_get_str (char *buf
, int base
, mpz_t a
)
733 static char tohex
[] = "0123456789abcdef";
735 mp_limb_t alimb
, *ap
;
745 buf
= xmalloc (ABSIZ (a
) * (GMP_LIMB_BITS
/ 4) + 3);
756 bn
= an
* (GMP_LIMB_BITS
/ 4);
759 for (i
= 0; i
< an
; i
++)
762 for (j
= 0; j
< GMP_LIMB_BITS
/ 4; j
++)
765 *bp
= tohex
[alimb
& 0xF];
774 strstrip_leading_zeros (buf
);
779 mpz_out_str (FILE *file
, int base
, mpz_t a
)
786 str
= mpz_get_str (0, 16, a
);
791 /* Calculate r satisfying r*d == 1 mod 2^n. */
793 mpz_invert_2exp (mpz_t r
, mpz_t a
, unsigned long n
)
798 ASSERT (mpz_odd_p (a
));
800 mpz_init_set_ui (inv
, 1L);
803 for (i
= 1; i
< n
; i
++)
805 mpz_mul (prod
, inv
, a
);
806 if (mpz_tstbit (prod
, i
) != 0)
810 mpz_mul (prod
, inv
, a
);
811 mpz_tdiv_r_2exp (prod
, prod
, n
);
812 ASSERT (mpz_cmp_ui (prod
, 1L) == 0);
820 /* Calculate inv satisfying r*a == 1 mod 2^n. */
822 mpz_invert_ui_2exp (mpz_t r
, unsigned long a
, unsigned long n
)
825 mpz_init_set_ui (az
, a
);
826 mpz_invert_2exp (r
, az
, n
);
832 mpz_pow_ui (mpz_t x
, mpz_t y
, unsigned long z
)
836 mpz_init_set_ui (t
, 1);
845 mpz_addmul_ui (mpz_t x
, mpz_t y
, unsigned long z
)
850 mpz_mul_ui (t
, y
, z
);
855 /* x=floor(y^(1/z)) */
857 mpz_root (mpz_t x
, mpz_t y
, unsigned long z
)
863 fprintf (stderr
, "mpz_root does not accept negative values\n");
866 if (mpz_cmp_ui (y
, 1) <= 0)
875 mpz_pow_ui (t
, u
, z
- 1);
876 mpz_tdiv_q (t
, y
, t
);
877 mpz_addmul_ui (t
, u
, z
- 1);
878 mpz_tdiv_q_ui (t
, t
, z
);
879 if (mpz_cmp (t
, u
) >= 0)