gen-trialdivtab.c

   1 /* gen-trialdivtab.c
   2
   3    Contributed to the GNU project by Torbjorn Granlund.
   4
   5 Copyright 2009, 2012 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 "bootstrap.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   return mpz_sgn (x) ? mpz_sizeinbase (x, 2) : 0;
 207 }
 208
 209 void
 210 mpn_mod_1s_4p_cps (mpz_t cps[7], mpz_t bparm)
 211 {
 212   mpz_t b, bi;
 213   mpz_t B1modb, B2modb, B3modb, B4modb, B5modb;
 214   mpz_t t;
 215   int cnt;
 216
 217   mpz_init_set (b, bparm);
 218
 219   cnt = limb_bits - mpz_log2 (b);
 220
 221   mpz_init (bi);
 222   mpz_init (t);
 223   mpz_init (B1modb);
 224   mpz_init (B2modb);
 225   mpz_init (B3modb);
 226   mpz_init (B4modb);
 227   mpz_init (B5modb);
 228
 229   mpz_set_ui (t, 1);
 230   mpz_mul_2exp (t, t, limb_bits - cnt);
 231   mpz_sub (t, t, b);
 232   mpz_mul_2exp (t, t, limb_bits);
 233   mpz_tdiv_q (bi, t, b);                /* bi = B^2/b, except msb */
 234
 235   mpz_set_ui (t, 1);
 236   mpz_mul_2exp (t, t, limb_bits);       /* t = B */
 237   mpz_tdiv_r (B1modb, t, b);
 238
 239   mpz_mul_2exp (t, B1modb, limb_bits);
 240   mpz_tdiv_r (B2modb, t, b);
 241
 242   mpz_mul_2exp (t, B2modb, limb_bits);
 243   mpz_tdiv_r (B3modb, t, b);
 244
 245   mpz_mul_2exp (t, B3modb, limb_bits);
 246   mpz_tdiv_r (B4modb, t, b);
 247
 248   mpz_mul_2exp (t, B4modb, limb_bits);
 249   mpz_tdiv_r (B5modb, t, b);
 250
 251   mpz_set (cps[0], bi);
 252   mpz_set_ui (cps[1], cnt);
 253   mpz_tdiv_q_2exp (cps[2], B1modb, 0);
 254   mpz_tdiv_q_2exp (cps[3], B2modb, 0);
 255   mpz_tdiv_q_2exp (cps[4], B3modb, 0);
 256   mpz_tdiv_q_2exp (cps[5], B4modb, 0);
 257   mpz_tdiv_q_2exp (cps[6], B5modb, 0);
 258
 259   mpz_clear (b);
 260   mpz_clear (bi);
 261   mpz_clear (t);
 262   mpz_clear (B1modb);
 263   mpz_clear (B2modb);
 264   mpz_clear (B3modb);
 265   mpz_clear (B4modb);
 266   mpz_clear (B5modb);
 267 }
 268
 269 int
 270 sumspills (mpz_t ppp, mpz_t *a, int n)
 271 {
 272   mpz_t s;
 273   int i, ret;
 274
 275   mpz_init_set (s, a[0]);
 276
 277   for (i = 1; i < n; i++)
 278     {
 279       mpz_add (s, s, a[i]);
 280     }
 281   ret = mpz_cmp (s, B) >= 0;
 282   mpz_clear (s);
 283
 284   return ret;
 285 }