gen-trialdivtab.c

   1 /* gen-trialdivtab.c
   2
   3    Contributed to the GNU project by Torbjorn Granlund.
   4
   5 Copyright 2009, 2012, 2013 Free Software Foundation, Inc.
   6
   7 This file is part of the GNU MP Library.
   8
   9 The GNU MP Library is free software; you can redistribute it and/or modify
  10 it under the terms of either:
  11
  12   * the GNU Lesser General Public License as published by the Free
  13     Software Foundation; either version 3 of the License, or (at your
  14     option) any later version.
  15
  16 or
  17
  18   * the GNU General Public License as published by the Free Software
  19     Foundation; either version 2 of the License, or (at your option) any
  20     later version.
  21
  22 or both in parallel, as here.
  23
  24 The GNU MP Library is distributed in the hope that it will be useful, but
  25 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  26 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  27 for more details.
  28
  29 You should have received copies of the GNU General Public License and the
  30 GNU Lesser General Public License along with the GNU MP Library.  If not,
  31 see https://www.gnu.org/licenses/.  */
  32
  33 /*
  34   Generate tables for fast, division-free trial division for GMP.
  35
  36   There is one main table, ptab.  It contains primes, multiplied together, and
  37   several types of pre-computed inverses.  It refers to tables of the type
  38   dtab, via the last two indices.  That table contains the individual primes in
  39   the range, except that the primes are not actually included in the table (see
  40   the P macro; it sneakingly excludes the primes themselves).  Instead, the
  41   dtab tables contains tuples for each prime (modular-inverse, limit) used for
  42   divisibility checks.
  43
  44   This interface is not intended for division of very many primes, since then
  45   other algorithms apply.
  46 */
  47
  48 #include <stdlib.h>
  49 #include <stdio.h>
  50 #include "bootstrap.c"
  51
  52 int sumspills (mpz_t, mpz_t *, int);
  53 void mpn_mod_1s_4p_cps (mpz_t [7], mpz_t);
  54
  55 int limb_bits;
  56
  57 mpz_t B;
  58
  59 int
  60 main (int argc, char *argv[])
  61 {
  62   unsigned long t, p;
  63   mpz_t ppp, acc, inv, gmp_numb_max, tmp, Bhalf;
  64   mpz_t pre[7];
  65   int i;
  66   int start_p, end_p, interval_start, interval_end, omitted_p;
  67   char *endtok;
  68   int stop;
  69   int np, start_idx;
  70
  71   if (argc < 2)
  72     {
  73       fprintf (stderr, "usage: %s bits endprime\n", argv[0]);
  74       exit (1);
  75     }
  76
  77   limb_bits = atoi (argv[1]);
  78
  79   end_p = 1290;                 /* default end prime */
  80   if (argc == 3)
  81     end_p = atoi (argv[2]);
  82
  83   printf ("#if GMP_LIMB_BITS != %d\n", limb_bits);
  84   printf ("#error This table is for GMP_LIMB_BITS = %d\n", limb_bits);
  85   printf ("#endif\n\n");
  86
  87   printf ("#if GMP_NAIL_BITS != 0\n");
  88   printf ("#error This table does not support nails\n");
  89   printf ("#endif\n\n");
  90
  91   for (i = 0; i < 7; i++)
  92     mpz_init (pre[i]);
  93
  94   mpz_init_set_ui (gmp_numb_max, 1);
  95   mpz_mul_2exp (gmp_numb_max, gmp_numb_max, limb_bits);
  96   mpz_sub_ui (gmp_numb_max, gmp_numb_max, 1);
  97
  98   mpz_init (tmp);
  99   mpz_init (inv);
 100
 101   mpz_init_set_ui (B, 1);  mpz_mul_2exp (B, B, limb_bits);
 102   mpz_init_set_ui (Bhalf, 1);  mpz_mul_2exp (Bhalf, Bhalf, limb_bits - 1);
 103
 104   start_p = 3;
 105
 106   mpz_init_set_ui (ppp, 1);
 107   mpz_init (acc);
 108   interval_start = start_p;
 109   omitted_p = 3;
 110   interval_end = 0;
 111
 112 /*  printf ("static struct gmp_primes_dtab gmp_primes_dtab[] = {\n"); */
 113
 114   printf ("#ifdef WANT_dtab\n");
 115
 116   for (t = start_p; t <= end_p; t += 2)
 117     {
 118       if (! isprime (t))
 119         continue;
 120
 121       mpz_mul_ui (acc, ppp, t);
 122       stop = mpz_cmp (acc, Bhalf) >= 0;
 123       if (!stop)
 124         {
 125           mpn_mod_1s_4p_cps (pre, acc);
 126           stop = sumspills (acc, pre + 2, 5);
 127         }
 128
 129       if (stop)
 130         {
 131           for (p = interval_start; p <= interval_end; p += 2)
 132             {
 133               if (! isprime (p))
 134                 continue;
 135
 136               printf ("  P(%d,", (int) p);
 137               mpz_invert_ui_2exp (inv, p, limb_bits);
 138               printf ("CNST_LIMB(0x");  mpz_out_str (stdout, 16, inv);  printf ("),");
 139
 140               mpz_tdiv_q_ui (tmp, gmp_numb_max, p);
 141               printf ("CNST_LIMB(0x");  mpz_out_str (stdout, 16, tmp);
 142               printf (")),\n");
 143             }
 144           mpz_set_ui (ppp, t);
 145           interval_start = t;
 146           omitted_p = t;
 147         }
 148       else
 149         {
 150           mpz_set (ppp, acc);
 151         }
 152       interval_end = t;
 153     }
 154   printf ("#define SMALLEST_OMITTED_PRIME %d\n", (int) omitted_p);
 155   printf ("#endif\n");
 156
 157   printf ("#ifdef WANT_ptab\n");
 158
 159 /*  printf ("static struct gmp_primes_ptab gmp_primes_ptab[] = {\n"); */
 160
 161   endtok = "";
 162
 163   mpz_set_ui (ppp, 1);
 164   interval_start = start_p;
 165   interval_end = 0;
 166   np = 0;
 167   start_idx = 0;
 168   for (t = start_p; t <= end_p; t += 2)
 169     {
 170       if (! isprime (t))
 171         continue;
 172
 173       mpz_mul_ui (acc, ppp, t);
 174
 175       stop = mpz_cmp (acc, Bhalf) >= 0;
 176       if (!stop)
 177         {
 178           mpn_mod_1s_4p_cps (pre, acc);
 179           stop = sumspills (acc, pre + 2, 5);
 180         }
 181
 182       if (stop)
 183         {
 184           mpn_mod_1s_4p_cps (pre, ppp);
 185           printf ("%s", endtok);
 186           printf ("  {CNST_LIMB(0x");  mpz_out_str (stdout, 16, ppp);
 187           printf ("),{CNST_LIMB(0x");  mpz_out_str (stdout, 16, pre[0]);
 188           printf ("),%d", (int) PTR(pre[1])[0]);
 189           for (i = 0; i < 5; i++)
 190             {
 191               printf (",");
 192               printf ("CNST_LIMB(0x");  mpz_out_str (stdout, 16, pre[2 + i]);
 193               printf (")");
 194             }
 195           printf ("},");
 196           printf ("%d,", start_idx);
 197           printf ("%d}", np - start_idx);
 198
 199           endtok = ",\n";
 200           mpz_set_ui (ppp, t);
 201           interval_start = t;
 202           start_idx = np;
 203         }
 204       else
 205         {
 206           mpz_set (ppp, acc);
 207         }
 208       interval_end = t;
 209       np++;
 210     }
 211
 212   printf ("\n");
 213   printf ("#endif\n");
 214
 215   return 0;
 216 }
 217
 218 unsigned long
 219 mpz_log2 (mpz_t x)
 220 {
 221   return mpz_sgn (x) ? mpz_sizeinbase (x, 2) : 0;
 222 }
 223
 224 void
 225 mpn_mod_1s_4p_cps (mpz_t cps[7], mpz_t bparm)
 226 {
 227   mpz_t b, bi;
 228   mpz_t B1modb, B2modb, B3modb, B4modb, B5modb;
 229   mpz_t t;
 230   int cnt;
 231
 232   mpz_init_set (b, bparm);
 233
 234   cnt = limb_bits - mpz_log2 (b);
 235
 236   mpz_init (bi);
 237   mpz_init (t);
 238   mpz_init (B1modb);
 239   mpz_init (B2modb);
 240   mpz_init (B3modb);
 241   mpz_init (B4modb);
 242   mpz_init (B5modb);
 243
 244   mpz_set_ui (t, 1);
 245   mpz_mul_2exp (t, t, limb_bits - cnt);
 246   mpz_sub (t, t, b);
 247   mpz_mul_2exp (t, t, limb_bits);
 248   mpz_tdiv_q (bi, t, b);                /* bi = B^2/b, except msb */
 249
 250   mpz_set_ui (t, 1);
 251   mpz_mul_2exp (t, t, limb_bits);       /* t = B */
 252   mpz_tdiv_r (B1modb, t, b);
 253
 254   mpz_mul_2exp (t, B1modb, limb_bits);
 255   mpz_tdiv_r (B2modb, t, b);
 256
 257   mpz_mul_2exp (t, B2modb, limb_bits);
 258   mpz_tdiv_r (B3modb, t, b);
 259
 260   mpz_mul_2exp (t, B3modb, limb_bits);
 261   mpz_tdiv_r (B4modb, t, b);
 262
 263   mpz_mul_2exp (t, B4modb, limb_bits);
 264   mpz_tdiv_r (B5modb, t, b);
 265
 266   mpz_set (cps[0], bi);
 267   mpz_set_ui (cps[1], cnt);
 268   mpz_tdiv_q_2exp (cps[2], B1modb, 0);
 269   mpz_tdiv_q_2exp (cps[3], B2modb, 0);
 270   mpz_tdiv_q_2exp (cps[4], B3modb, 0);
 271   mpz_tdiv_q_2exp (cps[5], B4modb, 0);
 272   mpz_tdiv_q_2exp (cps[6], B5modb, 0);
 273
 274   mpz_clear (b);
 275   mpz_clear (bi);
 276   mpz_clear (t);
 277   mpz_clear (B1modb);
 278   mpz_clear (B2modb);
 279   mpz_clear (B3modb);
 280   mpz_clear (B4modb);
 281   mpz_clear (B5modb);
 282 }
 283
 284 int
 285 sumspills (mpz_t ppp, mpz_t *a, int n)
 286 {
 287   mpz_t s;
 288   int i, ret;
 289
 290   mpz_init_set (s, a[0]);
 291
 292   for (i = 1; i < n; i++)
 293     {
 294       mpz_add (s, s, a[i]);
 295     }
 296   ret = mpz_cmp (s, B) >= 0;
 297   mpz_clear (s);
 298
 299   return ret;
 300 }