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 either:
   9
  10   * the GNU Lesser General Public License as published by the Free
  11     Software Foundation; either version 3 of the License, or (at your
  12     option) any later version.
  13
  14 or
  15
  16   * the GNU General Public License as published by the Free Software
  17     Foundation; either version 2 of the License, or (at your option) any
  18     later version.
  19
  20 or both in parallel, as here.
  21
  22 The GNU MP Library is distributed in the hope that it will be useful, but
  23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  24 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  25 for more details.
  26
  27 You should have received copies of the GNU General Public License and the
  28 GNU Lesser General Public License along with the GNU MP Library.  If not,
  29 see https://www.gnu.org/licenses/.  */
  30
  31
  32 #include "mini-gmp/mini-gmp.c"
  33
  34 #define MIN(l,o) ((l) < (o) ? (l) : (o))
  35 #define PTR(x)   ((x)->_mp_d)
  36 #define SIZ(x)   ((x)->_mp_size)
  37
  38 #define xmalloc gmp_default_alloc
  39
  40 int
  41 isprime (unsigned long int t)
  42 {
  43   unsigned long int q, r, d;
  44
  45   if (t < 32)
  46     return (0xa08a28acUL >> t) & 1;
  47   if ((t & 1) == 0)
  48     return 0;
  49
  50   if (t % 3 == 0)
  51     return 0;
  52   if (t % 5 == 0)
  53     return 0;
  54   if (t % 7 == 0)
  55     return 0;
  56
  57   for (d = 11;;)
  58     {
  59       q = t / d;
  60       r = t - q * d;
  61       if (q < d)
  62         return 1;
  63       if (r == 0)
  64         break;
  65       d += 2;
  66       q = t / d;
  67       r = t - q * d;
  68       if (q < d)
  69         return 1;
  70       if (r == 0)
  71         break;
  72       d += 4;
  73     }
  74   return 0;
  75 }
  76
  77 int
  78 log2_ceil (int n)
  79 {
  80   int  e;
  81   assert (n >= 1);
  82   for (e = 0; ; e++)
  83     if ((1 << e) >= n)
  84       break;
  85   return e;
  86 }
  87
  88 /* Set inv to the inverse of d, in the style of invert_limb, ie. for
  89    udiv_qrnnd_preinv.  */
  90 void
  91 mpz_preinv_invert (mpz_t inv, mpz_t d, int numb_bits)
  92 {
  93   mpz_t  t;
  94   int    norm;
  95   assert (SIZ(d) > 0);
  96
  97   norm = numb_bits - mpz_sizeinbase (d, 2);
  98   assert (norm >= 0);
  99   mpz_init_set_ui (t, 1L);
 100   mpz_mul_2exp (t, t, 2*numb_bits - norm);
 101   mpz_tdiv_q (inv, t, d);
 102   mpz_set_ui (t, 1L);
 103   mpz_mul_2exp (t, t, numb_bits);
 104   mpz_sub (inv, inv, t);
 105
 106   mpz_clear (t);
 107 }
 108
 109 /* Calculate r satisfying r*d == 1 mod 2^n. */
 110 void
 111 mpz_invert_2exp (mpz_t r, mpz_t a, unsigned long n)
 112 {
 113   unsigned long  i;
 114   mpz_t  inv, prod;
 115
 116   assert (mpz_odd_p (a));
 117
 118   mpz_init_set_ui (inv, 1L);
 119   mpz_init (prod);
 120
 121   for (i = 1; i < n; i++)
 122     {
 123       mpz_mul (prod, inv, a);
 124       if (mpz_tstbit (prod, i) != 0)
 125         mpz_setbit (inv, i);
 126     }
 127
 128   mpz_mul (prod, inv, a);
 129   mpz_tdiv_r_2exp (prod, prod, n);
 130   assert (mpz_cmp_ui (prod, 1L) == 0);
 131
 132   mpz_set (r, inv);
 133
 134   mpz_clear (inv);
 135   mpz_clear (prod);
 136 }
 137
 138 /* Calculate inv satisfying r*a == 1 mod 2^n. */
 139 void
 140 mpz_invert_ui_2exp (mpz_t r, unsigned long a, unsigned long n)
 141 {
 142   mpz_t  az;
 143   mpz_init_set_ui (az, a);
 144   mpz_invert_2exp (r, az, n);
 145   mpz_clear (az);
 146 }