bootstrap.c

   1 /* Functions needed for bootstrapping the gmp build, based on mini-gmp.
   2
   3 Copyright 2001, 2002, 2004, 2011, 2012 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   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.
  11
  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.
  16
  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/.  */
  19
  20
  21 #include "mini-gmp/mini-gmp.c"
  22
  23 #define MIN(l,o) ((l) < (o) ? (l) : (o))
  24 #define PTR(x)   ((x)->_mp_d)
  25 #define SIZ(x)   ((x)->_mp_size)
  26
  27 #define xmalloc gmp_default_alloc
  28
  29 int
  30 isprime (unsigned long int t)
  31 {
  32   unsigned long int q, r, d;
  33
  34   if (t < 32)
  35     return (0xa08a28acUL >> t) & 1;
  36   if ((t & 1) == 0)
  37     return 0;
  38
  39   if (t % 3 == 0)
  40     return 0;
  41   if (t % 5 == 0)
  42     return 0;
  43   if (t % 7 == 0)
  44     return 0;
  45
  46   for (d = 11;;)
  47     {
  48       q = t / d;
  49       r = t - q * d;
  50       if (q < d)
  51         return 1;
  52       if (r == 0)
  53         break;
  54       d += 2;
  55       q = t / d;
  56       r = t - q * d;
  57       if (q < d)
  58         return 1;
  59       if (r == 0)
  60         break;
  61       d += 4;
  62     }
  63   return 0;
  64 }
  65
  66 int
  67 log2_ceil (int n)
  68 {
  69   int  e;
  70   assert (n >= 1);
  71   for (e = 0; ; e++)
  72     if ((1 << e) >= n)
  73       break;
  74   return e;
  75 }
  76
  77 /* Set inv to the inverse of d, in the style of invert_limb, ie. for
  78    udiv_qrnnd_preinv.  */
  79 void
  80 mpz_preinv_invert (mpz_t inv, mpz_t d, int numb_bits)
  81 {
  82   mpz_t  t;
  83   int    norm;
  84   assert (SIZ(d) > 0);
  85
  86   norm = numb_bits - mpz_sizeinbase (d, 2);
  87   assert (norm >= 0);
  88   mpz_init_set_ui (t, 1L);
  89   mpz_mul_2exp (t, t, 2*numb_bits - norm);
  90   mpz_tdiv_q (inv, t, d);
  91   mpz_set_ui (t, 1L);
  92   mpz_mul_2exp (t, t, numb_bits);
  93   mpz_sub (inv, inv, t);
  94
  95   mpz_clear (t);
  96 }
  97
  98 /* Calculate r satisfying r*d == 1 mod 2^n. */
  99 void
 100 mpz_invert_2exp (mpz_t r, mpz_t a, unsigned long n)
 101 {
 102   unsigned long  i;
 103   mpz_t  inv, prod;
 104
 105   assert (mpz_odd_p (a));
 106
 107   mpz_init_set_ui (inv, 1L);
 108   mpz_init (prod);
 109
 110   for (i = 1; i < n; i++)
 111     {
 112       mpz_mul (prod, inv, a);
 113       if (mpz_tstbit (prod, i) != 0)
 114         mpz_setbit (inv, i);
 115     }
 116
 117   mpz_mul (prod, inv, a);
 118   mpz_tdiv_r_2exp (prod, prod, n);
 119   assert (mpz_cmp_ui (prod, 1L) == 0);
 120
 121   mpz_set (r, inv);
 122
 123   mpz_clear (inv);
 124   mpz_clear (prod);
 125 }
 126
 127 /* Calculate inv satisfying r*a == 1 mod 2^n. */
 128 void
 129 mpz_invert_ui_2exp (mpz_t r, unsigned long a, unsigned long n)
 130 {
 131   mpz_t  az;
 132   mpz_init_set_ui (az, a);
 133   mpz_invert_2exp (r, az, n);
 134   mpz_clear (az);
 135 }