contrib/gmp/gen-trialdivtab.c

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