gen-psqr.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:17k
源码类别:

数学计算

开发平台:

Unix_Linux

  1. /* Generate perfect square testing data.
  2. Copyright 2002, 2003, 2004 Free Software Foundation, Inc.
  3. This file is part of the GNU MP Library.
  4. The GNU MP Library is free software; you can redistribute it and/or modify
  5. it under the terms of the GNU Lesser General Public License as published by
  6. the Free Software Foundation; either version 3 of the License, or (at your
  7. option) any later version.
  8. The GNU MP Library is distributed in the hope that it will be useful, but
  9. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  10. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  11. License for more details.
  12. You should have received a copy of the GNU Lesser General Public License
  13. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  14. #include <stdio.h>
  15. #include <stdlib.h>
  16. #include "dumbmp.c"
  17. /* The aim of this program is to choose either mpn_mod_34lsub1 or mpn_mod_1
  18.    (plus a PERFSQR_PP modulus), and generate tables indicating quadratic
  19.    residues and non-residues modulo small factors of that modulus.
  20.    For the usual 32 or 64 bit cases mpn_mod_34lsub1 gets used.  That
  21.    function exists specifically because 2^24-1 and 2^48-1 have nice sets of
  22.    prime factors.  For other limb sizes it's considered, but if it doesn't
  23.    have good factors then mpn_mod_1 will be used instead.
  24.    When mpn_mod_1 is used, the modulus PERFSQR_PP is created from a
  25.    selection of small primes, chosen to fill PERFSQR_MOD_BITS of a limb,
  26.    with that bit count chosen so (2*GMP_LIMB_BITS)*2^PERFSQR_MOD_BITS <=
  27.    GMP_LIMB_MAX, allowing PERFSQR_MOD_IDX in mpn/generic/perfsqr.c to do its
  28.    calculation within a single limb.
  29.    In either case primes can be combined to make divisors.  The table data
  30.    then effectively indicates remainders which are quadratic residues mod
  31.    all the primes.  This sort of combining reduces the number of steps
  32.    needed after mpn_mod_34lsub1 or mpn_mod_1, saving code size and time.
  33.    Nothing is gained or lost in terms of detections, the same total fraction
  34.    of non-residues will be identified.
  35.    Nothing particularly sophisticated is attempted for combining factors to
  36.    make divisors.  This is probably a kind of knapsack problem so it'd be
  37.    too hard to attempt anything completely general.  For the usual 32 and 64
  38.    bit limbs we get a good enough result just pairing the biggest and
  39.    smallest which fit together, repeatedly.
  40.    Another aim is to get powerful combinations, ie. divisors which identify
  41.    biggest fraction of non-residues, and have those run first.  Again for
  42.    the usual 32 and 64 bits it seems good enough just to pair for big
  43.    divisors then sort according to the resulting fraction of non-residues
  44.    identified.
  45.    Also in this program, a table sq_res_0x100 of residues modulo 256 is
  46.    generated.  This simply fills bits into limbs of the appropriate
  47.    build-time GMP_LIMB_BITS each.
  48. */
  49. /* Normally we aren't using const in gen*.c programs, so as not to have to
  50.    bother figuring out if it works, but using it with f_cmp_divisor and
  51.    f_cmp_fraction avoids warnings from the qsort calls. */
  52. /* Same tests as gmp.h. */
  53. #if  defined (__STDC__)                                 
  54.   || defined (__cplusplus)                              
  55.   || defined (_AIX)                                     
  56.   || defined (__DECC)                                   
  57.   || (defined (__mips) && defined (_SYSTYPE_SVR4))      
  58.   || defined (_MSC_VER)                                 
  59.   || defined (_WIN32)
  60. #define HAVE_CONST        1
  61. #endif
  62. #if ! HAVE_CONST
  63. #define const
  64. #endif
  65. mpz_t  *sq_res_0x100;          /* table of limbs */
  66. int    nsq_res_0x100;          /* elements in sq_res_0x100 array */
  67. int    sq_res_0x100_num;       /* squares in sq_res_0x100 */
  68. double sq_res_0x100_fraction;  /* sq_res_0x100_num / 256 */
  69. int     mod34_bits;        /* 3*GMP_NUMB_BITS/4 */
  70. int     mod_bits;          /* bits from PERFSQR_MOD_34 or MOD_PP */
  71. int     max_divisor;       /* all divisors <= max_divisor */
  72. int     max_divisor_bits;  /* ceil(log2(max_divisor)) */
  73. double  total_fraction;    /* of squares */
  74. mpz_t   pp;                /* product of primes, or 0 if mod_34lsub1 used */
  75. mpz_t   pp_norm;           /* pp shifted so NUMB high bit set */
  76. mpz_t   pp_inverted;       /* invert_limb style inverse */
  77. mpz_t   mod_mask;          /* 2^mod_bits-1 */
  78. char    mod34_excuse[128]; /* why mod_34lsub1 not used (if it's not) */
  79. /* raw list of divisors of 2^mod34_bits-1 or pp, just to show in a comment */
  80. struct rawfactor_t {
  81.   int     divisor;
  82.   int     multiplicity;
  83. };
  84. struct rawfactor_t  *rawfactor;
  85. int                 nrawfactor;
  86. /* factors of 2^mod34_bits-1 or pp and associated data, after combining etc */
  87. struct factor_t {
  88.   int     divisor;
  89.   mpz_t   inverse;   /* 1/divisor mod 2^mod_bits */
  90.   mpz_t   mask;      /* indicating squares mod divisor */
  91.   double  fraction;  /* squares/total */
  92. };
  93. struct factor_t  *factor;
  94. int              nfactor;       /* entries in use in factor array */
  95. int              factor_alloc;  /* entries allocated to factor array */
  96. int
  97. f_cmp_divisor (const void *parg, const void *qarg)
  98. {
  99.   const struct factor_t *p, *q;
  100.   p = parg;
  101.   q = qarg;
  102.   if (p->divisor > q->divisor)
  103.     return 1;
  104.   else if (p->divisor < q->divisor)
  105.     return -1;
  106.   else
  107.     return 0;
  108. }
  109. int
  110. f_cmp_fraction (const void *parg, const void *qarg)
  111. {
  112.   const struct factor_t *p, *q;
  113.   p = parg;
  114.   q = qarg;
  115.   if (p->fraction > q->fraction)
  116.     return 1;
  117.   else if (p->fraction < q->fraction)
  118.     return -1;
  119.   else
  120.     return 0;
  121. }
  122. /* Remove array[idx] by copying the remainder down, and adjust narray
  123.    accordingly.  */
  124. #define COLLAPSE_ELEMENT(array, idx, narray)                    
  125.   do {                                                          
  126.     mem_copyi ((char *) &(array)[idx],                          
  127.                (char *) &(array)[idx+1],                        
  128.                ((narray)-((idx)+1)) * sizeof (array[0]));       
  129.     (narray)--;                                                 
  130.   } while (0)
  131. /* return n*2^p mod m */
  132. int
  133. mul_2exp_mod (int n, int p, int m)
  134. {
  135.   int  i;
  136.   for (i = 0; i < p; i++)
  137.     n = (2 * n) % m;
  138.   return n;
  139. }
  140. /* return -n mod m */
  141. int
  142. neg_mod (int n, int m)
  143. {
  144.   ASSERT (n >= 0 && n < m);
  145.   return (n == 0 ? 0 : m-n);
  146. }
  147. /* Set "mask" to a value such that "mask & (1<<idx)" is non-zero if
  148.    "-(idx<<mod_bits)" can be a square modulo m.  */
  149. void
  150. square_mask (mpz_t mask, int m)
  151. {
  152.   int    p, i, r, idx;
  153.   p = mul_2exp_mod (1, mod_bits, m);
  154.   p = neg_mod (p, m);
  155.   mpz_set_ui (mask, 0L);
  156.   for (i = 0; i < m; i++)
  157.     {
  158.       r = (i * i) % m;
  159.       idx = (r * p) % m;
  160.       mpz_setbit (mask, (unsigned long) idx);
  161.     }
  162. }
  163. void
  164. generate_sq_res_0x100 (int limb_bits)
  165. {
  166.   int  i, res;
  167.   nsq_res_0x100 = (0x100 + limb_bits - 1) / limb_bits;
  168.   sq_res_0x100 = (mpz_t *) xmalloc (nsq_res_0x100 * sizeof (*sq_res_0x100));
  169.   for (i = 0; i < nsq_res_0x100; i++)
  170.     mpz_init_set_ui (sq_res_0x100[i], 0L);
  171.   for (i = 0; i < 0x100; i++)
  172.     {
  173.       res = (i * i) % 0x100;
  174.       mpz_setbit (sq_res_0x100[res / limb_bits],
  175.                   (unsigned long) (res % limb_bits));
  176.     }
  177.   sq_res_0x100_num = 0;
  178.   for (i = 0; i < nsq_res_0x100; i++)
  179.     sq_res_0x100_num += mpz_popcount (sq_res_0x100[i]);
  180.   sq_res_0x100_fraction = (double) sq_res_0x100_num / 256.0;
  181. }
  182. void
  183. generate_mod (int limb_bits, int nail_bits)
  184. {
  185.   int    numb_bits = limb_bits - nail_bits;
  186.   int    i, divisor;
  187.   mpz_init_set_ui (pp, 0L);
  188.   mpz_init_set_ui (pp_norm, 0L);
  189.   mpz_init_set_ui (pp_inverted, 0L);
  190.   /* no more than limb_bits many factors in a one limb modulus (and of
  191.      course in reality nothing like that many) */
  192.   factor_alloc = limb_bits;
  193.   factor = (struct factor_t *) xmalloc (factor_alloc * sizeof (*factor));
  194.   rawfactor = (struct rawfactor_t *)
  195.     xmalloc (factor_alloc * sizeof (*rawfactor));
  196.   if (numb_bits % 4 != 0)
  197.     {
  198.       strcpy (mod34_excuse, "GMP_NUMB_BITS % 4 != 0");
  199.       goto use_pp;
  200.     }
  201.   max_divisor = 2*limb_bits;
  202.   max_divisor_bits = log2_ceil (max_divisor);
  203.   if (numb_bits / 4 < max_divisor_bits)
  204.     {
  205.       /* Wind back to one limb worth of max_divisor, if that will let us use
  206.          mpn_mod_34lsub1.  */
  207.       max_divisor = limb_bits;
  208.       max_divisor_bits = log2_ceil (max_divisor);
  209.       if (numb_bits / 4 < max_divisor_bits)
  210.         {
  211.           strcpy (mod34_excuse, "GMP_NUMB_BITS / 4 too small");
  212.           goto use_pp;
  213.         }
  214.     }
  215.   {
  216.     /* Can use mpn_mod_34lsub1, find small factors of 2^mod34_bits-1. */
  217.     mpz_t  m, q, r;
  218.     int    multiplicity;
  219.     mod34_bits = (numb_bits / 4) * 3;
  220.     /* mpn_mod_34lsub1 returns a full limb value, PERFSQR_MOD_34 folds it at
  221.        the mod34_bits mark, adding the two halves for a remainder of at most
  222.        mod34_bits+1 many bits */
  223.     mod_bits = mod34_bits + 1;
  224.     mpz_init_set_ui (m, 1L);
  225.     mpz_mul_2exp (m, m, mod34_bits);
  226.     mpz_sub_ui (m, m, 1L);
  227.     mpz_init (q);
  228.     mpz_init (r);
  229.     for (i = 3; i <= max_divisor; i++)
  230.       {
  231.         if (! isprime (i))
  232.           continue;
  233.         mpz_tdiv_qr_ui (q, r, m, (unsigned long) i);
  234.         if (mpz_sgn (r) != 0)
  235.           continue;
  236.         /* if a repeated prime is found it's used as an i^n in one factor */
  237.         divisor = 1;
  238.         multiplicity = 0;
  239.         do
  240.           {
  241.             if (divisor > max_divisor / i)
  242.               break;
  243.             multiplicity++;
  244.             mpz_set (m, q);
  245.             mpz_tdiv_qr_ui (q, r, m, (unsigned long) i);
  246.           }
  247.         while (mpz_sgn (r) == 0);
  248.         ASSERT (nrawfactor < factor_alloc);
  249.         rawfactor[nrawfactor].divisor = i;
  250.         rawfactor[nrawfactor].multiplicity = multiplicity;
  251.         nrawfactor++;
  252.       }
  253.     mpz_clear (m);
  254.     mpz_clear (q);
  255.     mpz_clear (r);
  256.   }
  257.   if (nrawfactor <= 2)
  258.     {
  259.       mpz_t  new_pp;
  260.       sprintf (mod34_excuse, "only %d small factor%s",
  261.                nrawfactor, nrawfactor == 1 ? "" : "s");
  262.     use_pp:
  263.       /* reset to two limbs of max_divisor, in case the mpn_mod_34lsub1 code
  264.          tried with just one */
  265.       max_divisor = 2*limb_bits;
  266.       max_divisor_bits = log2_ceil (max_divisor);
  267.       mpz_init (new_pp);
  268.       nrawfactor = 0;
  269.       mod_bits = MIN (numb_bits, limb_bits - max_divisor_bits);
  270.       /* one copy of each small prime */
  271.       mpz_set_ui (pp, 1L);
  272.       for (i = 3; i <= max_divisor; i++)
  273.         {
  274.           if (! isprime (i))
  275.             continue;
  276.           mpz_mul_ui (new_pp, pp, (unsigned long) i);
  277.           if (mpz_sizeinbase (new_pp, 2) > mod_bits)
  278.             break;
  279.           mpz_set (pp, new_pp);
  280.           ASSERT (nrawfactor < factor_alloc);
  281.           rawfactor[nrawfactor].divisor = i;
  282.           rawfactor[nrawfactor].multiplicity = 1;
  283.           nrawfactor++;
  284.         }
  285.       /* Plus an extra copy of one or more of the primes selected, if that
  286.          still fits in max_divisor and the total in mod_bits.  Usually only
  287.          3 or 5 will be candidates */
  288.       for (i = nrawfactor-1; i >= 0; i--)
  289.         {
  290.           if (rawfactor[i].divisor > max_divisor / rawfactor[i].divisor)
  291.             continue;
  292.           mpz_mul_ui (new_pp, pp, (unsigned long) rawfactor[i].divisor);
  293.           if (mpz_sizeinbase (new_pp, 2) > mod_bits)
  294.             continue;
  295.           mpz_set (pp, new_pp);
  296.           rawfactor[i].multiplicity++;
  297.         }
  298.       mod_bits = mpz_sizeinbase (pp, 2);
  299.       mpz_set (pp_norm, pp);
  300.       while (mpz_sizeinbase (pp_norm, 2) < numb_bits)
  301.         mpz_add (pp_norm, pp_norm, pp_norm);
  302.       mpz_preinv_invert (pp_inverted, pp_norm, numb_bits);
  303.       mpz_clear (new_pp);
  304.     }
  305.   /* start the factor array */
  306.   for (i = 0; i < nrawfactor; i++)
  307.     {
  308.       int  j;
  309.       ASSERT (nfactor < factor_alloc);
  310.       factor[nfactor].divisor = 1;
  311.       for (j = 0; j < rawfactor[i].multiplicity; j++)
  312.         factor[nfactor].divisor *= rawfactor[i].divisor;
  313.       nfactor++;
  314.     }
  315.  combine:
  316.   /* Combine entries in the factor array.  Combine the smallest entry with
  317.      the biggest one that will fit with it (ie. under max_divisor), then
  318.      repeat that with the new smallest entry. */
  319.   qsort (factor, nfactor, sizeof (factor[0]), f_cmp_divisor);
  320.   for (i = nfactor-1; i >= 1; i--)
  321.     {
  322.       if (factor[i].divisor <= max_divisor / factor[0].divisor)
  323.         {
  324.           factor[0].divisor *= factor[i].divisor;
  325.           COLLAPSE_ELEMENT (factor, i, nfactor);
  326.           goto combine;
  327.         }
  328.     }
  329.   total_fraction = 1.0;
  330.   for (i = 0; i < nfactor; i++)
  331.     {
  332.       mpz_init (factor[i].inverse);
  333.       mpz_invert_ui_2exp (factor[i].inverse,
  334.                           (unsigned long) factor[i].divisor,
  335.                           (unsigned long) mod_bits);
  336.       mpz_init (factor[i].mask);
  337.       square_mask (factor[i].mask, factor[i].divisor);
  338.       /* fraction of possible squares */
  339.       factor[i].fraction = (double) mpz_popcount (factor[i].mask)
  340.         / factor[i].divisor;
  341.       /* total fraction of possible squares */
  342.       total_fraction *= factor[i].fraction;
  343.     }
  344.   /* best tests first (ie. smallest fraction) */
  345.   qsort (factor, nfactor, sizeof (factor[0]), f_cmp_fraction);
  346. }
  347. void
  348. print (int limb_bits, int nail_bits)
  349. {
  350.   int    i;
  351.   mpz_t  mhi, mlo;
  352.   printf ("/* This file generated by gen-psqr.c - DO NOT EDIT. */n");
  353.   printf ("n");
  354.   printf ("#if GMP_LIMB_BITS != %d || GMP_NAIL_BITS != %dn",
  355.           limb_bits, nail_bits);
  356.   printf ("Error, error, this data is for %d bit limb and %d bit nailn",
  357.           limb_bits, nail_bits);
  358.   printf ("#endifn");
  359.   printf ("n");
  360.   printf ("/* Non-zero bit indicates a quadratic residue mod 0x100.n");
  361.   printf ("   This test identifies %.2f%% as non-squares (%d/256). */n",
  362.           (1.0 - sq_res_0x100_fraction) * 100.0,
  363.           0x100 - sq_res_0x100_num);
  364.   printf ("static const mp_limb_tn");
  365.   printf ("sq_res_0x100[%d] = {n", nsq_res_0x100);
  366.   for (i = 0; i < nsq_res_0x100; i++)
  367.     {
  368.       printf ("  CNST_LIMB(0x");
  369.       mpz_out_str (stdout, 16, sq_res_0x100[i]);
  370.       printf ("),n");
  371.     }
  372.   printf ("};n");
  373.   printf ("n");
  374.   if (mpz_sgn (pp) != 0)
  375.     {
  376.       printf ("/* mpn_mod_34lsub1 not used due to %s */n", mod34_excuse);
  377.       printf ("/* PERFSQR_PP = ");
  378.     }
  379.   else
  380.     printf ("/* 2^%d-1 = ", mod34_bits);
  381.   for (i = 0; i < nrawfactor; i++)
  382.     {
  383.       if (i != 0)
  384.         printf (" * ");
  385.       printf ("%d", rawfactor[i].divisor);
  386.       if (rawfactor[i].multiplicity != 1)
  387.         printf ("^%d", rawfactor[i].multiplicity);
  388.     }
  389.   printf (" %s*/n", mpz_sgn (pp) == 0 ? "... " : "");
  390.   printf ("#define PERFSQR_MOD_BITS  %dn", mod_bits);
  391.   if (mpz_sgn (pp) != 0)
  392.     {
  393.       printf ("#define PERFSQR_PP            CNST_LIMB(0x");
  394.       mpz_out_str (stdout, 16, pp);
  395.       printf (")n");
  396.       printf ("#define PERFSQR_PP_NORM       CNST_LIMB(0x");
  397.       mpz_out_str (stdout, 16, pp_norm);
  398.       printf (")n");
  399.       printf ("#define PERFSQR_PP_INVERTED   CNST_LIMB(0x");
  400.       mpz_out_str (stdout, 16, pp_inverted);
  401.       printf (")n");
  402.     }
  403.   printf ("n");
  404.   mpz_init (mhi);
  405.   mpz_init (mlo);
  406.   printf ("/* This test identifies %.2f%% as non-squares. */n",
  407.           (1.0 - total_fraction) * 100.0);
  408.   printf ("#define PERFSQR_MOD_TEST(up, usize) \n");
  409.   printf ("  do {                              \n");
  410.   printf ("    mp_limb_t  r;                   \n");
  411.   if (mpz_sgn (pp) != 0)
  412.     printf ("    PERFSQR_MOD_PP (r, up, usize);  \n");
  413.   else
  414.     printf ("    PERFSQR_MOD_34 (r, up, usize);  \n");
  415.   for (i = 0; i < nfactor; i++)
  416.     {
  417.       printf ("                                    \n");
  418.       printf ("    /* %5.2f%% */                    \n",
  419.               (1.0 - factor[i].fraction) * 100.0);
  420.       printf ("    PERFSQR_MOD_%d (r, CNST_LIMB(%2d), CNST_LIMB(0x",
  421.               factor[i].divisor <= limb_bits ? 1 : 2,
  422.               factor[i].divisor);
  423.       mpz_out_str (stdout, 16, factor[i].inverse);
  424.       printf ("), \n");
  425.       printf ("                   CNST_LIMB(0x");
  426.       if ( factor[i].divisor <= limb_bits)
  427.         {
  428.           mpz_out_str (stdout, 16, factor[i].mask);
  429.         }
  430.       else
  431.         {
  432.           mpz_tdiv_r_2exp (mlo, factor[i].mask, (unsigned long) limb_bits);
  433.           mpz_tdiv_q_2exp (mhi, factor[i].mask, (unsigned long) limb_bits);
  434.           mpz_out_str (stdout, 16, mhi);
  435.           printf ("), CNST_LIMB(0x");
  436.           mpz_out_str (stdout, 16, mlo);
  437.         }
  438.       printf (")); \n");
  439.     }
  440.   printf ("  } while (0)n");
  441.   printf ("n");
  442.   printf ("/* Grand total sq_res_0x100 and PERFSQR_MOD_TEST, %.2f%% non-squares. */n",
  443.           (1.0 - (total_fraction * 44.0/256.0)) * 100.0);
  444.   printf ("n");
  445.   printf ("/* helper for tests/mpz/t-perfsqr.c */n");
  446.   printf ("#define PERFSQR_DIVISORS  { 256,");
  447.   for (i = 0; i < nfactor; i++)
  448.       printf (" %d,", factor[i].divisor);
  449.   printf (" }n");
  450.   mpz_clear (mhi);
  451.   mpz_clear (mlo);
  452. }
  453. int
  454. main (int argc, char *argv[])
  455. {
  456.   int  limb_bits, nail_bits;
  457.   if (argc != 3)
  458.     {
  459.       fprintf (stderr, "Usage: gen-psqr <limbbits> <nailbits>n");
  460.       exit (1);
  461.     }
  462.   limb_bits = atoi (argv[1]);
  463.   nail_bits = atoi (argv[2]);
  464.   if (limb_bits <= 0
  465.       || nail_bits < 0
  466.       || nail_bits >= limb_bits)
  467.     {
  468.       fprintf (stderr, "Invalid limb/nail bits: %d %dn",
  469.                limb_bits, nail_bits);
  470.       exit (1);
  471.     }
  472.   generate_sq_res_0x100 (limb_bits);
  473.   generate_mod (limb_bits, nail_bits);
  474.   print (limb_bits, nail_bits);
  475.   return 0;
  476. }