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

数学计算

开发平台:

Unix_Linux

  1. /* mpn/gcd.c: mpn_gcd for gcd of two odd integers.
  2. Copyright 1991, 1993, 1994, 1995, 1996, 1997, 1998, 2000, 2001, 2002, 2003,
  3. 2004, 2005, 2008 Free Software Foundation, Inc.
  4. This file is part of the GNU MP Library.
  5. The GNU MP Library is free software; you can redistribute it and/or modify
  6. it under the terms of the GNU Lesser General Public License as published by
  7. the Free Software Foundation; either version 3 of the License, or (at your
  8. option) any later version.
  9. The GNU MP Library is distributed in the hope that it will be useful, but
  10. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  11. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  12. License for more details.
  13. You should have received a copy of the GNU Lesser General Public License
  14. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  15. #include "gmp.h"
  16. #include "gmp-impl.h"
  17. #include "longlong.h"
  18. /* Uses the HGCD operation described in
  19.      N. M鰈ler, On Sch鰊hage's algorithm and subquadratic integer gcd
  20.      computation, Math. Comp. 77 (2008), 589-607.
  21.   to reduce inputs until they are of size below GCD_DC_THRESHOLD, and
  22.   then uses Lehmer's algorithm.
  23. */
  24. /* Some reasonable choices are n / 2 (same as in hgcd), and p = (n +
  25.  * 2)/3, which gives a balanced multiplication in
  26.  * mpn_hgcd_matrix_adjust. However, p = 2 n/3 gives slightly better
  27.  * performance. The matrix-vector multiplication is then
  28.  * 4:1-unbalanced, with matrix elements of size n/6, and vector
  29.  * elements of size p = 2n/3. */
  30. /* From analysis of the theoretical running time, it appears that when
  31.  * multiplication takes time O(n^alpha), p should be chosen so that
  32.  * the ratio of the time for the mpn_hgcd call, and the time for the
  33.  * multiplication in mpn_hgcd_matrix_adjust, is roughly 1/(alpha -
  34.  * 1). */
  35. #ifdef TUNE_GCD_P
  36. #define P_TABLE_SIZE 10000
  37. mp_size_t p_table[P_TABLE_SIZE];
  38. #define CHOOSE_P(n) ( (n) < P_TABLE_SIZE ? p_table[n] : 2*(n)/3)
  39. #else
  40. #define CHOOSE_P(n) (2*(n) / 3)
  41. #endif
  42. mp_size_t
  43. mpn_gcd (mp_ptr gp, mp_ptr up, mp_size_t usize, mp_ptr vp, mp_size_t n)
  44. {
  45.   mp_size_t talloc;
  46.   mp_size_t scratch;
  47.   mp_size_t matrix_scratch;
  48.   mp_size_t gn;
  49.   mp_ptr tp;
  50.   TMP_DECL;
  51.   /* FIXME: Check for small sizes first, before setting up temporary
  52.      storage etc. */
  53.   talloc = MPN_GCD_LEHMER_N_ITCH(n);
  54.   /* For initial division */
  55.   scratch = usize - n + 1;
  56.   if (scratch > talloc)
  57.     talloc = scratch;
  58. #if TUNE_GCD_P
  59.   if (CHOOSE_P (n) > 0)
  60. #else
  61.   if (ABOVE_THRESHOLD (n, GCD_DC_THRESHOLD))
  62. #endif
  63.     {
  64.       mp_size_t hgcd_scratch;
  65.       mp_size_t update_scratch;
  66.       mp_size_t p = CHOOSE_P (n);
  67.       mp_size_t scratch;
  68. #if TUNE_GCD_P
  69.       /* Worst case, since we don't guarantee that n - CHOOSE_P(n)
  70.  is increasing */
  71.       matrix_scratch = MPN_HGCD_MATRIX_INIT_ITCH (n);
  72.       hgcd_scratch = mpn_hgcd_itch (n);
  73.       update_scratch = 2*(n - 1);
  74. #else
  75.       matrix_scratch = MPN_HGCD_MATRIX_INIT_ITCH (n - p);
  76.       hgcd_scratch = mpn_hgcd_itch (n - p);
  77.       update_scratch = p + n - 1;
  78. #endif
  79.       scratch = matrix_scratch + MAX(hgcd_scratch, update_scratch);
  80.       if (scratch > talloc)
  81. talloc = scratch;
  82.     }
  83.   TMP_MARK;
  84.   tp = TMP_ALLOC_LIMBS(talloc);
  85.   if (usize > n)
  86.     {
  87.       mpn_tdiv_qr (tp, up, 0, up, usize, vp, n);
  88.       if (mpn_zero_p (up, n))
  89. {
  90.   MPN_COPY (gp, vp, n);
  91.   TMP_FREE;
  92.   return n;
  93. }
  94.     }
  95. #if TUNE_GCD_P
  96.   while (CHOOSE_P (n) > 0)
  97. #else
  98.   while (ABOVE_THRESHOLD (n, GCD_DC_THRESHOLD))
  99. #endif
  100.     {
  101.       struct hgcd_matrix M;
  102.       mp_size_t p = CHOOSE_P (n);
  103.       mp_size_t matrix_scratch = MPN_HGCD_MATRIX_INIT_ITCH (n - p);
  104.       mp_size_t nn;
  105.       mpn_hgcd_matrix_init (&M, n - p, tp);
  106.       nn = mpn_hgcd (up + p, vp + p, n - p, &M, tp + matrix_scratch);
  107.       if (nn > 0)
  108. {
  109.   ASSERT (M.n <= (n - p - 1)/2);
  110.   ASSERT (M.n + p <= (p + n - 1) / 2);
  111.   /* Temporary storage 2 (p + M->n) <= p + n - 1. */
  112.   n = mpn_hgcd_matrix_adjust (&M, p + nn, up, vp, p, tp + matrix_scratch);
  113. }
  114.       else
  115. {
  116.   /* Temporary storage n */
  117.   n = mpn_gcd_subdiv_step (gp, &gn, up, vp, n, tp);
  118.   if (n == 0)
  119.     {
  120.       TMP_FREE;
  121.       return gn;
  122.     }
  123. }
  124.     }
  125.   gn = mpn_gcd_lehmer_n (gp, up, vp, n, tp);
  126.   TMP_FREE;
  127.   return gn;
  128. }
  129. #ifdef TUNE_GCD_P
  130. #include <stdio.h>
  131. #include <string.h>
  132. #include <time.h>
  133. #include "speed.h"
  134. static int
  135. compare_double(const void *ap, const void *bp)
  136. {
  137.   double a = * (const double *) ap;
  138.   double b = * (const double *) bp;
  139.   if (a < b)
  140.     return -1;
  141.   else if (a > b)
  142.     return 1;
  143.   else
  144.     return 0;
  145. }
  146. static double
  147. median (double *v, size_t n)
  148. {
  149.   qsort(v, n, sizeof(*v), compare_double);
  150.   return v[n/2];
  151. }
  152. #define TIME(res, code) do {
  153.   double time_measurement[5];
  154.   unsigned time_i;
  155.   for (time_i = 0; time_i < 5; time_i++)
  156.     {
  157.       speed_starttime();
  158.       code;
  159.       time_measurement[time_i] = speed_endtime();
  160.     }
  161.   res = median(time_measurement, 5);
  162. } while (0)
  163. int
  164. main(int argc, char *argv)
  165. {
  166.   gmp_randstate_t rands;
  167.   mp_size_t n;
  168.   mp_ptr ap;
  169.   mp_ptr bp;
  170.   mp_ptr up;
  171.   mp_ptr vp;
  172.   mp_ptr gp;
  173.   mp_ptr tp;
  174.   TMP_DECL;
  175.   /* Unbuffered so if output is redirected to a file it isn't lost if the
  176.      program is killed part way through.  */
  177.   setbuf (stdout, NULL);
  178.   setbuf (stderr, NULL);
  179.   gmp_randinit_default (rands);
  180.   TMP_MARK;
  181.   ap = TMP_ALLOC_LIMBS (P_TABLE_SIZE);
  182.   bp = TMP_ALLOC_LIMBS (P_TABLE_SIZE);
  183.   up = TMP_ALLOC_LIMBS (P_TABLE_SIZE);
  184.   vp = TMP_ALLOC_LIMBS (P_TABLE_SIZE);
  185.   gp = TMP_ALLOC_LIMBS (P_TABLE_SIZE);
  186.   tp = TMP_ALLOC_LIMBS (MPN_GCD_LEHMER_N_ITCH (P_TABLE_SIZE));
  187.   mpn_random (ap, P_TABLE_SIZE);
  188.   mpn_random (bp, P_TABLE_SIZE);
  189.   memset (p_table, 0, sizeof(p_table));
  190.   for (n = 100; n++; n < P_TABLE_SIZE)
  191.     {
  192.       mp_size_t p;
  193.       mp_size_t best_p;
  194.       double best_time;
  195.       double lehmer_time;
  196.       if (ap[n-1] == 0)
  197. ap[n-1] = 1;
  198.       if (bp[n-1] == 0)
  199. bp[n-1] = 1;
  200.       p_table[n] = 0;
  201.       TIME(lehmer_time, {
  202.   MPN_COPY (up, ap, n);
  203.   MPN_COPY (vp, bp, n);
  204.   mpn_gcd_lehmer_n (gp, up, vp, n, tp);
  205. });
  206.       best_time = lehmer_time;
  207.       best_p = 0;
  208.       for (p = n * 0.48; p < n * 0.77; p++)
  209. {
  210.   double t;
  211.   p_table[n] = p;
  212.   TIME(t, {
  213.       MPN_COPY (up, ap, n);
  214.       MPN_COPY (vp, bp, n);
  215.       mpn_gcd (gp, up, n, vp, n);
  216.     });
  217.   if (t < best_time)
  218.     {
  219.       best_time = t;
  220.       best_p = p;
  221.     }
  222. }
  223.       printf("%6d %6d %5.3g", n, best_p, (double) best_p / n);
  224.       if (best_p > 0)
  225. {
  226.   double speedup = 100 * (lehmer_time - best_time) / lehmer_time;
  227.   printf(" %5.3g%%", speedup);
  228.   if (speedup < 1.0)
  229.     {
  230.       printf(" (ignored)");
  231.       best_p = 0;
  232.     }
  233. }
  234.       printf("n");
  235.       p_table[n] = best_p;
  236.     }
  237.   TMP_FREE;
  238.   gmp_randclear(rands);
  239.   return 0;
  240. }
  241. #endif /* TUNE_GCD_P */