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

数学计算

开发平台:

Unix_Linux

  1. /* Factoring with Pollard's rho method.
  2. Copyright 1995, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2009
  3. Free Software Foundation, Inc.
  4. This program is free software; you can redistribute it and/or modify it under
  5. the terms of the GNU General Public License as published by the Free Software
  6. Foundation; either version 3 of the License, or (at your option) any later
  7. version.
  8. This program is distributed in the hope that it will be useful, but WITHOUT ANY
  9. WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
  10. PARTICULAR PURPOSE.  See the GNU General Public License for more details.
  11. You should have received a copy of the GNU General Public License along with
  12. this program.  If not, see http://www.gnu.org/licenses/.  */
  13. #include <stdlib.h>
  14. #include <stdio.h>
  15. #include <string.h>
  16. #include "gmp.h"
  17. int flag_verbose = 0;
  18. static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
  19. void
  20. factor_using_division (mpz_t t, unsigned int limit)
  21. {
  22.   mpz_t q, r;
  23.   unsigned long int f;
  24.   int ai;
  25.   unsigned *addv = add;
  26.   unsigned int failures;
  27.   if (flag_verbose > 0)
  28.     {
  29.       printf ("[trial division (%u)] ", limit);
  30.       fflush (stdout);
  31.     }
  32.   mpz_init (q);
  33.   mpz_init (r);
  34.   f = mpz_scan1 (t, 0);
  35.   mpz_div_2exp (t, t, f);
  36.   while (f)
  37.     {
  38.       printf ("2 ");
  39.       fflush (stdout);
  40.       --f;
  41.     }
  42.   for (;;)
  43.     {
  44.       mpz_tdiv_qr_ui (q, r, t, 3);
  45.       if (mpz_cmp_ui (r, 0) != 0)
  46. break;
  47.       mpz_set (t, q);
  48.       printf ("3 ");
  49.       fflush (stdout);
  50.     }
  51.   for (;;)
  52.     {
  53.       mpz_tdiv_qr_ui (q, r, t, 5);
  54.       if (mpz_cmp_ui (r, 0) != 0)
  55. break;
  56.       mpz_set (t, q);
  57.       printf ("5 ");
  58.       fflush (stdout);
  59.     }
  60.   failures = 0;
  61.   f = 7;
  62.   ai = 0;
  63.   while (mpz_cmp_ui (t, 1) != 0)
  64.     {
  65.       mpz_tdiv_qr_ui (q, r, t, f);
  66.       if (mpz_cmp_ui (r, 0) != 0)
  67. {
  68.   f += addv[ai];
  69.   if (mpz_cmp_ui (q, f) < 0)
  70.     break;
  71.   ai = (ai + 1) & 7;
  72.   failures++;
  73.   if (failures > limit)
  74.     break;
  75. }
  76.       else
  77. {
  78.   mpz_swap (t, q);
  79.   printf ("%lu ", f);
  80.   fflush (stdout);
  81.   failures = 0;
  82. }
  83.     }
  84.   mpz_clears (q, r, NULL);
  85. }
  86. void
  87. factor_using_division_2kp (mpz_t t, unsigned int limit, unsigned long p)
  88. {
  89.   mpz_t r;
  90.   mpz_t f;
  91.   unsigned int k;
  92.   if (flag_verbose > 0)
  93.     {
  94.       printf ("[trial division (%u)] ", limit);
  95.       fflush (stdout);
  96.     }
  97.   mpz_init (r);
  98.   mpz_init_set_ui (f, 2 * p);
  99.   mpz_add_ui (f, f, 1);
  100.   for (k = 1; k < limit; k++)
  101.     {
  102.       mpz_tdiv_r (r, t, f);
  103.       while (mpz_cmp_ui (r, 0) == 0)
  104. {
  105.   mpz_tdiv_q (t, t, f);
  106.   mpz_tdiv_r (r, t, f);
  107.   mpz_out_str (stdout, 10, f);
  108.   fflush (stdout);
  109.   fputc (' ', stdout);
  110. }
  111.       mpz_add_ui (f, f, 2 * p);
  112.     }
  113.   mpz_clears (f, r, NULL);
  114. }
  115. void
  116. factor_using_pollard_rho (mpz_t n, unsigned long a, unsigned long p)
  117. {
  118.   mpz_t x, x1, y, P;
  119.   mpz_t t1, t2;
  120.   unsigned long long k, l, i;
  121.   if (flag_verbose > 0)
  122.     {
  123.       printf ("[pollard-rho (%lu)] ", a);
  124.       fflush (stdout);
  125.     }
  126.   mpz_inits (t1, t2, NULL);
  127.   mpz_init_set_si (y, 2);
  128.   mpz_init_set_si (x, 2);
  129.   mpz_init_set_si (x1, 2);
  130.   mpz_init_set_ui (P, 1);
  131.   k = 1;
  132.   l = 1;
  133.   while (mpz_cmp_ui (n, 1) != 0)
  134.     {
  135.       for (;;)
  136. {
  137.   do
  138.     {
  139.       if (p != 0)
  140. {
  141.   mpz_powm_ui (x, x, p, n);
  142.   mpz_add_ui (x, x, a);
  143. }
  144.       else
  145. {
  146.   mpz_mul (t1, x, x);
  147.   mpz_mod (x, t1, n);
  148.   mpz_add_ui (x, x, a);
  149. }
  150.       mpz_sub (t1, x1, x);
  151.       mpz_mul (t2, P, t1);
  152.       mpz_mod (P, t2, n);
  153.       if (k % 32 == 1)
  154. {
  155.   mpz_gcd (t1, P, n);
  156.   if (mpz_cmp_ui (t1, 1) != 0)
  157.     goto factor_found;
  158.   mpz_set (y, x);
  159. }
  160.     }
  161.   while (--k != 0);
  162.   mpz_gcd (t1, P, n);
  163.   if (mpz_cmp_ui (t1, 1) != 0)
  164.     goto factor_found;
  165.   mpz_set (x1, x);
  166.   k = l;
  167.   l = 2 * l;
  168.   for (i = 0; i < k; i++)
  169.     {
  170.       if (p != 0)
  171. {
  172.   mpz_powm_ui (x, x, p, n);
  173.   mpz_add_ui (x, x, a);
  174. }
  175.       else
  176. {
  177.   mpz_mul (t1, x, x);
  178.   mpz_mod (x, t1, n);
  179.   mpz_add_ui (x, x, a);
  180. }
  181.     }
  182.   mpz_set (y, x);
  183. }
  184.     factor_found:
  185.       do
  186. {
  187.   if (p != 0)
  188.     {
  189.       mpz_powm_ui (y, y, p, n); mpz_add_ui (y, y, a);
  190.     }
  191.   else
  192.     {
  193.       mpz_mul (t1, y, y);
  194.       mpz_mod (y, t1, n);
  195.       mpz_add_ui (y, y, a);
  196.     }
  197.   mpz_sub (t1, x1, y);
  198.   mpz_gcd (t1, t1, n);
  199. }
  200.       while (mpz_cmp_ui (t1, 1) == 0);
  201.       mpz_divexact (n, n, t1); /* divide by t1, before t1 is overwritten */
  202.       if (!mpz_probab_prime_p (t1, 10))
  203. {
  204.   do
  205.     {
  206.       mp_limb_t a_limb;
  207.       mpn_random (&a_limb, (mp_size_t) 1);
  208.       a = a_limb;
  209.     }
  210.   while (a == 0);
  211.   if (flag_verbose > 0)
  212.     {
  213.       printf ("[composite factor--restarting pollard-rho] ");
  214.       fflush (stdout);
  215.     }
  216.   factor_using_pollard_rho (t1, a, p);
  217. }
  218.       else
  219. {
  220.   mpz_out_str (stdout, 10, t1);
  221.   fflush (stdout);
  222.   fputc (' ', stdout);
  223. }
  224.       mpz_mod (x, x, n);
  225.       mpz_mod (x1, x1, n);
  226.       mpz_mod (y, y, n);
  227.       if (mpz_probab_prime_p (n, 10))
  228. {
  229.   mpz_out_str (stdout, 10, n);
  230.   fflush (stdout);
  231.   fputc (' ', stdout);
  232.   break;
  233. }
  234.     }
  235.   mpz_clears (P, t2, t1, x1, x, y, NULL);
  236. }
  237. void
  238. factor (mpz_t t, unsigned long p)
  239. {
  240.   unsigned int division_limit;
  241.   if (mpz_sgn (t) == 0)
  242.     return;
  243.   /* Set the trial division limit according the size of t.  */
  244.   division_limit = mpz_sizeinbase (t, 2);
  245.   if (division_limit > 1000)
  246.     division_limit = 1000 * 1000;
  247.   else
  248.     division_limit = division_limit * division_limit;
  249.   if (p != 0)
  250.     factor_using_division_2kp (t, division_limit / 10, p);
  251.   else
  252.     factor_using_division (t, division_limit);
  253.   if (mpz_cmp_ui (t, 1) != 0)
  254.     {
  255.       if (flag_verbose > 0)
  256. {
  257.   printf ("[is number prime?] ");
  258.   fflush (stdout);
  259. }
  260.       if (mpz_probab_prime_p (t, 10))
  261. mpz_out_str (stdout, 10, t);
  262.       else
  263. factor_using_pollard_rho (t, 1L, p);
  264.     }
  265. }
  266. int
  267. main (int argc, char *argv[])
  268. {
  269.   mpz_t t;
  270.   unsigned long p;
  271.   int i;
  272.   if (argc > 1 && !strcmp (argv[1], "-v"))
  273.     {
  274.       flag_verbose = 1;
  275.       argv++;
  276.       argc--;
  277.     }
  278.   if (argc > 1 && !strcmp (argv[1], "-q"))
  279.     {
  280.       flag_verbose = -1;
  281.       argv++;
  282.       argc--;
  283.     }
  284.   mpz_init (t);
  285.   if (argc > 1)
  286.     {
  287.       p = 0;
  288.       for (i = 1; i < argc; i++)
  289. {
  290.   if (!strncmp (argv[i], "-Mp", 3))
  291.     {
  292.       p = atoi (argv[i] + 3);
  293.       mpz_set_ui (t, 1);
  294.       mpz_mul_2exp (t, t, p);
  295.       mpz_sub_ui (t, t, 1);
  296.     }
  297.   else if (!strncmp (argv[i], "-2kp", 4))
  298.     {
  299.       p = atoi (argv[i] + 4);
  300.       continue;
  301.     }
  302.   else
  303.     {
  304.       mpz_set_str (t, argv[i], 0);
  305.     }
  306.   if (mpz_cmp_ui (t, 0) == 0)
  307.     puts ("-");
  308.   else
  309.     {
  310.       factor (t, p);
  311.       puts ("");
  312.     }
  313. }
  314.     }
  315.   else
  316.     {
  317.       for (;;)
  318. {
  319.   mpz_inp_str (t, stdin, 0);
  320.   if (feof (stdin))
  321.     break;
  322.   if (flag_verbose >= 0)
  323.     {
  324.       mpz_out_str (stdout, 10, t); printf (" = ");
  325.     }
  326.   factor (t, 0);
  327.   puts ("");
  328. }
  329.     }
  330.   exit (0);
  331. }