Revision e4f0108028701f575c3bf69e567316e3150f7dd1 authored by Paul Zimmermann on 13 October 2010, 08:19:28 UTC, committed by Paul Zimmermann on 13 October 2010, 08:19:28 UTC
[bwc] added "krylov" or "mksol" in ETA log lines git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/cado-nfs/trunk@352 3eaf19af-ecc0-40f6-b43f-672aa0c71c71
1 parent ebbb4f9
phi.c
/* Make SNFS polynomials for cyclotomic numbers.
Copyright 2005, 2006, 2007, 2008 Alexander Kruppa.
This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
more details.
You should have received a copy of the GNU General Public License along
with this program; see the file COPYING. If not, write to the Free
Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
02111-1307, USA.
*/
/* Version: 0.1.2.2
History: 0.1 initial release
0.1.1 Added GGNFS format, estimate of cost and time,
changed x to mpz_t type.
0.1.2 Makes better polynomials for numbers with composite
base and no useful algebraic factors
0.1.2.1 Allow reading N from stdin
0.1.2.2 Allow printing polynomial in msieve format. Make better
polynomials if exponent is 3 (mod 6). Fixed sign in
printed number (i.e. was 2^128-1 instead of 2^128+1).
Todo: Detect Aurifeullian factors and make polynomials for those.
Compute root properties that depend on x.
Test some more to see if any of this actually works.
*/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gmp.h>
#include <assert.h>
#define VERSION "0.1.2.1"
#define uint unsigned int
/* Output formats */
#define FRANKE 0
#define CWI 1
#define GGNFS 2
#define MSIEVE 3
static void
mpz_add_si (mpz_t R, mpz_t S, long i)
{
if (i < 0)
mpz_sub_ui (R, S, (unsigned long) -i);
else
mpz_add_ui (R, S, (unsigned long) i);
}
/* Compute p[0]^v[0] * ... * p[k-1]^v[k-1] */
static void
mpz_mul_list (mpz_t e, int *p, int *v, uint k)
{
int i;
uint j;
mpz_set_ui (e, 1);
for (j = 0; j < k; j++)
for (i = 0; i < v[j]; i++)
mpz_mul_ui (e, e, p[j]);
}
/* Find prime factors < maxp of R. Cofactor is put in R,
the i-th prime factor in p_i with exponent v_i.
Returns number of distinct prime factors found. */
static uint
trialdiv (mpz_t R, int *p, int *v, uint maxp, uint flen)
{
uint i, k = 0;
#ifdef DEBUG
printf ("trialdiv: ");
#endif
for (i = 2; i < maxp && k < flen; i = (i + 1) | 1) /* 2 3 5 7 9... */
if (mpz_divisible_ui_p (R, i))
{
p[k] = i;
v[k] = 0;
do
{
v[k]++;
mpz_divexact_ui (R, R, i);
} while (mpz_divisible_ui_p (R, i));
#ifdef DEBUG
printf ("%d^%d ", i, v[k]);
#endif
k++;
}
#ifdef DEBUG
printf ("\n");
#endif
return k;
}
static void
usage ()
{
printf ("phi \n" VERSION);
printf ("Usage: phi [options] [n x] N\n");
printf ("Makes SNFS polynomial for factoring the primitive part Phi_n(x) of x^n-1\n");
printf ("Cofactor N must divide x^n-1 if n is odd, or x^(n/2)+1 if n is even\n");
printf ("If you do not specify x and n, the program will try to auto-detect them.\n");
printf ("Options:\n");
printf ("-deg4 Make a degree 4 polynomial\n");
printf ("-deg5 Make a degree 5 polynomial\n");
printf ("-deg6 Make a degree 6 polynomial\n");
printf (" By default chooses the degree that minimises the SNFS difficulty\n");
printf ("-franke Use Franke file format for polynomial\n");
printf ("-ggnfs Use GGNFS file format for polynomial\n");
printf ("-cwi Use CWI file format for polynomial\n");
printf ("-msieve Use msieve file format for polynomial\n");
printf ("-short Omit unneeded lines for Franke file format\n");
printf ("- Read N from stdin\n");
}
static int
autodet (int *x, uint *n, mpz_t N)
{
const int flen = 32;
const int maxp = 2000; /* Limit for prime divisors */
const int maxbase = 1000; /* Limit for bases to try */
int p[flen], v[flen]; /* p are primes, v exponents */
int b, i, j, k;
mpz_t R; /* Remaining cofactor */
mpz_t e; /* exponent */
mpz_init (R);
mpz_sub_ui (R, N, 1);
mpz_init (e);
k = trialdiv (R, p, v, maxp, flen);
mpz_mul_list (e, p, v, k);
#ifdef DEBUG
gmp_printf ("\ne = %Zd\n", e);
#endif
for (b = 2; b <= maxbase; b++)
{
/* Skip perfect powers */
mpz_set_ui (R, b);
if (mpz_perfect_power_p (R))
continue;
#ifdef DEBUG
/* printf ("Trying base %d\n", b); */
#endif
mpz_powm (R, R, e, N);
if (mpz_cmp_ui (R, 1) == 0)
{
#ifdef DEBUG
printf ("Discovered! base = %d\n", b);
#endif
break;
}
}
if (b > maxbase)
{
/* No suitable base found */
mpz_clear (R);
mpz_clear (e);
return 0;
}
/* Find exponent */
for (j = 0; j < k; j++)
{
while (v[j] > 0)
{
/* See if lowering this exponent by one still satisfies ord|e */
v[j]--;
mpz_mul_list (e, p, v, k);
#ifdef DEBUG
gmp_printf ("Trying if ord | e = %Zd\n", e);
#endif
mpz_set_ui (R, b);
mpz_powm (R, R, e, N);
if (mpz_cmp_ui (R, 1) != 0)
{
/* No, cannot decrease this exponent any further. */
v[j]++; /* Undo decrease */
#ifdef DEBUG
printf ("No, leaving exponent of %d at %d\n", p[j], v[j]);
#endif
break; /* Try next prime */
}
#ifdef DEBUG
else
printf ("Yes, decreasing exponent of %d to %d\n", p[j], v[j]);
#endif
}
}
mpz_mul_list (e, p, v, k);
if (mpz_fits_uint_p (e))
{
if (x != NULL)
*x = b;
if (n != NULL)
*n = mpz_get_ui (e);
i = 1;
}
else
i = 0;
mpz_clear (R);
mpz_clear (e);
return i;
}
/* Difficulty must be in base 10 */
static double
snfscost (const double difficulty)
{
const double c = 1.5262856567; /* (32/9)^(1/3) */
const double logd = difficulty * 2.3025850930; /* *log(10) */
return exp (c * pow (logd, 1./3) * pow (log (logd), 2./3));
}
int
main (int argc, char **argv)
{
int format = FRANKE; /* Output format */
mpz_t x;
int n = 0; /* We factor phi_n(x) */
int k; /* k|n */
int degree = 0;
int sign, i, j;
int halved = 0; /* Did we use the degree halving trick? */
int shortform = 0; /* Short output for Franke format? */
int read_stdin = 0; /* Read number from stdin? */
mpz_t N; /* The cofactor of the number to factor */
mpz_t f[7]; /* Coefficients on algebraic side, max degree = 6 */
mpz_t g[2]; /* Coefficients on rational side, max degree = 1 */
mpz_t M; /* Common root */
mpz_t t; /* A temporary */
double dx; /* x as a double */
double difficulty;
double skewness;
double alpha;
const double timefactor = 1. / 2.1e15;
/* Parse command line parameters */
while (argc > 1 && argv[1][0] == '-')
{
if (strcmp (argv[1], "-cwi") == 0)
{
format = CWI;
argc--;
argv++;
}
else if (strcmp (argv[1], "-franke") == 0)
{
format = FRANKE;
argc--;
argv++;
}
else if (strcmp (argv[1], "-ggnfs") == 0)
{
format = GGNFS;
argc--;
argv++;
}
else if (strcmp (argv[1], "-msieve") == 0)
{
format = MSIEVE;
argc--;
argv++;
}
else if (strcmp (argv[1], "-short") == 0)
{
shortform = 1;
argc--;
argv++;
}
else if (strcmp (argv[1], "-deg4") == 0)
{
degree = 4;
argc--;
argv++;
}
else if (strcmp (argv[1], "-deg5") == 0)
{
degree = 5;
argc--;
argv++;
}
else if (strcmp (argv[1], "-deg6") == 0)
{
degree = 6;
argc--;
argv++;
}
else if (strcmp (argv[1], "-") == 0)
{
read_stdin = 1;
argc--;
argv++;
}
else
{
fprintf (stderr, "Error, unknown command line option: %s\n",
argv[1]);
usage();
exit (EXIT_FAILURE);
}
}
mpz_init (N);
mpz_init (x);
if (argc + read_stdin == 2)
{
int xi = 0;
if (read_stdin)
mpz_inp_str (N, stdin, 10);
else
mpz_set_str (N, argv[1], 10);
n = 0;
if (autodet (&xi, (uint *) &n, N))
{
mpz_set_ui (x, xi);
#ifdef DEBUG
printf ("Auto-detected Phi_%d(%d)\n", n, xi);
}
else
{
printf ("Auto-detectedion failed\n");
#endif
}
}
else if (argc + read_stdin > 3)
{
n = atoi (argv[1]);
mpz_set_str (x, argv[2], 10);
if (read_stdin)
mpz_inp_str (N, stdin, 10);
else
mpz_set_str (N, argv[3], 10);
}
if (n <= 0 || mpz_sgn (x) <= 0 || mpz_sgn (N) <= 0)
{
gmp_printf ("Error, invalid parameters: n = %d, x = %Zd\n", n, x);
usage();
exit (EXIT_FAILURE);
}
mpz_init (t);
mpz_init (g[0]);
mpz_init (g[1]);
for (i = 0; i < 7; i++)
mpz_init (f[i]);
mpz_init (M);
sign = -1;
/* Simplify even n, that is x^(n/2)+1 numbers */
if (n % 2 == 0)
{
sign = +1;
n /= 2;
}
/* Thus the number we factor is x^n+sign */
/* Check that N | x^n+sign */
mpz_powm_ui (t, x, n, N);
mpz_add_si (t, t, sign);
mpz_mod (t, t, N);
if (mpz_sgn (t) != 0)
{
gmp_fprintf (stderr, "Error, N does not divide %Zd^%d%c1\n",
x, n, sign < 0 ? '-' : '+');
exit (EXIT_FAILURE);
}
/* Find out what kind of number we have (algebraic factors etc.) and
choose polynomial */
dx = mpz_get_d (x);
if ((degree == 0 || degree == 4) && n % 15 == 0)
{
/* x^(15*k)+-1, remaining size 8/15 = 0.5333... */
degree = 4;
k = n / 15;
mpz_set_si (f[4], 1);
mpz_set_si (f[3], sign);
mpz_set_si (f[2], -4);
mpz_set_si (f[1], -sign * 4);
mpz_set_si (f[0], 1);
mpz_pow_ui (M, x, k);
halved = 1;
skewness = 1.;
difficulty = log10 (dx) * 8. * k;
alpha = 1.684;
}
else if ((degree == 0 || degree == 6) && n % 21 == 0)
{
/* x^(21*k)+-1, remaining size 12/21 = 0.571428... */
degree = 6;
k = n / 21;
mpz_set_si (f[6], 1);
mpz_set_si (f[5], sign);
mpz_set_si (f[4], -6);
mpz_set_si (f[3], -sign * 6);
mpz_set_si (f[2], 8);
mpz_set_si (f[1], sign * 8);
mpz_set_si (f[0], 1);
mpz_pow_ui (M, x, k);
halved = 1;
skewness = 1.;
difficulty = log10 (dx) * 12. * k;
alpha = 2.116;
}
else if ((degree == 0 || degree == 6) && n % 9 == 0)
{
/* x^(9*k)+-1 / x^(3*k)+-1 = x^(6k) -+ x^(3k) + 1
remaining size 2/3 = 0.666... */
degree = 6;
k = n / 9;
mpz_set_ui (f[6], 1);
mpz_set_si (f[3], -sign);
mpz_set_ui (f[0], 1);
mpz_pow_ui (M, x, k);
skewness = 1.;
difficulty = log10 (dx) * 6. * k;
alpha = 1.621;
}
else if ((degree == 0 || degree == 6) && n % 9 == 3)
{
/* x^(9*k+3)+-1, remaining size 2/3 = 0.666... */
degree = 6;
k = (n - 3) / 9;
mpz_mul (f[6], x, x);
mpz_mul_si (f[3], x, -sign);
mpz_set_ui (f[0], 1);
mpz_pow_ui (M, x, k);
skewness = pow (dx, -1./3.);
difficulty = log10 (dx) * (n / 3 * 2);
alpha = 0.;
}
else if ((degree == 0 || degree == 6) && n % 9 == 6)
{
/* x^(9*k+3)+-1, remaining size 2/3 = 0.666... */
degree = 6;
k = (n + 3) / 9;
mpz_set_ui (f[6], 1);
mpz_mul_si (f[3], x, -sign);
mpz_mul (f[0], x, x);
mpz_pow_ui (M, x, k);
skewness = pow (dx, 1./3.);
difficulty = log10 (dx) * (n / 3 * 2 + 2);
alpha = 0.;
}
else if ((degree == 0 || degree == 4) && n % 6 == 0)
{
/* x^(6*k)+1, remaining size 2/3 = 0.666... */
degree = 4;
k = n / 6;
mpz_set_si (f[4], 1);
mpz_set_si (f[2], -1);
mpz_set_si (f[0], 1);
mpz_pow_ui (M, x, k);
skewness = 1.;
difficulty = log10 (dx) * 4. * k;
alpha = 1.292;
}
else if ((degree == 0 || degree == 4) && n % 6 == 3)
{
/* x^(6*k+3)+-1, remaining size 2/3 = 0.666... */
degree = 4;
k = (n - 3) / 6;
mpz_mul (f[4], x, x);
mpz_mul_si (f[2], x, -sign);
mpz_set_si (f[0], 1);
mpz_pow_ui (M, x, k);
skewness = pow (dx, -0.5);
difficulty = log10 (dx) * 4. * k;
alpha = 0.;
}
else if ((degree == 0 || degree == 4) && n % 5 == 0)
{
/* x^(5*k)+-1, remaining size 4/5 = 0.8 */
degree = 4;
k = n / 5;
mpz_set_si (f[4], 1);
mpz_set_si (f[3], -sign);
mpz_set_si (f[2], 1);
mpz_set_si (f[1], -sign);
mpz_set_si (f[0], 1);
mpz_pow_ui (M, x, k);
skewness = 1.;
difficulty = log10 (dx) * 4. * k;
alpha = 1.453;
}
else if ((degree == 0 || degree == 6) && n % 7 == 0)
{
/* x^(7*k)+-1, remaining size 6/7 = 0.857142... */
degree = 6;
k = n / 7;
mpz_set_si (f[6], 1);
mpz_set_si (f[5], -sign);
mpz_set_si (f[4], 1);
mpz_set_si (f[3], -sign);
mpz_set_si (f[2], 1);
mpz_set_si (f[1], -sign);
mpz_set_si (f[0], 1);
mpz_pow_ui (M, x, k);
skewness = 1.;
difficulty = log10 (dx) * 6. * k;
alpha = 2.244;
}
else if ((degree == 0 || degree == 5) && n % 11 == 0)
{
/* x^(11*k)+-1, remaining size 10/11 = 0.90... */
degree = 5;
k = n / 11;
mpz_set_si (f[5], 1);
mpz_set_si (f[4], -sign*1);
mpz_set_si (f[3], -4);
mpz_set_si (f[2], sign*3);
mpz_set_si (f[1], 3);
mpz_set_si (f[0], -sign);
mpz_pow_ui (M, x, k);
halved = 1;
skewness = 1.;
difficulty = log10 (dx) * 10. * k;
alpha = 2.224;
}
else if ((degree == 0 || degree == 6) && n % 13 == 0)
{
/* x^(13*k)+-1, remaining size 12/13 = 0.923076... */
degree = 6;
k = n / 13;
mpz_set_si (f[6], 1);
mpz_set_si (f[5], -sign);
mpz_set_si (f[4], -5);
mpz_set_si (f[3], sign*4);
mpz_set_si (f[2], 6);
mpz_set_si (f[1], -sign*3);
mpz_set_si (f[0], -1);
mpz_pow_ui (M, x, k);
halved = 1;
skewness = 1.;
difficulty = log10 (dx) * 12. * k;
alpha = 3.095;
}
else
{
/* Make polynomial without using any algebraic factor */
const int d = (degree == 0) ? 5 : degree;
const int maxp = 10000;
const int flen = 20; /* Only examine 21 different factors
(including cofactor). Should easily suffice */
int p[flen], vp[flen], best = 0;
mpz_t R;
double rating, bestrating = 0.;
mpz_init_set (R, x);
k = trialdiv (R, p, vp, maxp, flen);
/* Now x = p[0] ^ v[0] * ... * p[k - 1] ^ v[k - 1] * R */
/* Here, k is the number of distinct prime factors */
for (i = 0; i < k; i++)
vp[i] *= n;
/* Now N = p[0] ^ v[0] * ... * p[k - 1] ^ v[k - 1] * R^n - 1 */
/* Make all exponents multiples of $d$.
When decreasing the exponent of $p$ by $u$, we need to multiply
$f_d$ by $p^u$.
When increasing the exponent of $p$ by $u$, we need to multiply
$f_0$ by $p^u$ and $M$ by $p$.
For each exponent $v$, we can choose between increasing by
$(-v) \bmod 5$ or decreasing by $v \bmod 5$.
So we have a search space of $2^{k+1}$ possibilities.
For now, we simply use $f_d + f_0 + c$ as the function to minimise,
where c is what we multiplied M by.
This rating function could and should be improved, but it seems to
usually produce decent polynomials.
*/
for (j = 0; j < (1 << (k + ((mpz_cmp_ui (R, 1) > 1) ? 1 : 0))); j++)
{
double fd = 1., f0 = 1., c = 1.;
for (i = 0; i < k; i++)
{
if (j & (1 << i))
{
/* Reduce exponent by v[i] % d */
fd *= pow ((double)p[i], (double)(vp[i] % d));
}
else if (vp[i] % d != 0)
{
/* Multiply by -v[i] % d */
f0 *= pow ((double)p[i], (double)(d - vp[i] % d));
c *= p[i];
}
}
if (j & (1 << k))
fd *= pow (mpz_get_d (R), (double)(n % d));
else
{
f0 *= pow (mpz_get_d (R), (double)((d - n % d) % d));
c *= mpz_get_d (R);
}
/* rating = c * c * pow (fd, 1. + 1./d) * pow (f0, 1. - 1./d); */
rating = fd + f0 + c;
#ifdef DEBUG
printf ("j = %d, %.0f * x^%d - %.0f, c = %.0f, rating = %f\n",
j, fd, d, f0, c, rating);
#endif
if (j == 0 || rating < bestrating)
{
best = j;
bestrating = rating;
#ifdef DEBUG
printf ("New best rating!\n");
#endif
}
}
degree = d;
difficulty = log10 (dx) * n;
mpz_set_ui (M, 1);
/* Now make the actual coefficients f_d and f_0 */
mpz_set_ui (f[d], 1);
mpz_set_ui (f[0], 1);
for (i = 0; i < k; i++)
{
if (best & (1 << i))
{
uint u = vp[i] % d;
mpz_ui_pow_ui (t, p[i], u);
mpz_mul (f[d], f[d], t);
mpz_ui_pow_ui (t, p[i], (vp[i] - u) / d);
mpz_mul (M, M, t);
}
else
{
uint u = (d - vp[i] % d) % d; /* u = (-vp[i]) % d */
mpz_ui_pow_ui (t, p[i], u);
mpz_mul (f[0], f[0], t);
difficulty += log10 ((double) p[i]) * u;
mpz_ui_pow_ui (t, p[i], (vp[i] + u) / d);
mpz_mul (M, M, t);
}
}
if (best & (1 << k))
{
uint u = n % d;
mpz_pow_ui (t, R, u);
mpz_mul (f[d], f[d], t);
mpz_pow_ui (t, R, (n - u) / d);
mpz_mul (M, M, t);
}
else
{
uint u = (d - n % d) % d;
mpz_pow_ui (t, R, u);
mpz_mul (f[0], f[0], t);
mpz_pow_ui (t, R, (n + u) / d);
mpz_mul (M, M, t);
difficulty += log10 (mpz_get_d (R)) * u;
}
skewness = pow (mpz_get_d (f[0]) / mpz_get_d (f[d]), 1. / d);
if (sign == -1)
mpz_neg (f[0], f[0]);
alpha = 0.; /* Alpha may depend on x, can't calculate yet */
}
if (halved)
{
mpz_neg (g[1], M); /* Y1 = -x^k */
mpz_mul (g[0], M, M);
mpz_add_ui (g[0], g[0], 1); /* Y0 = x^(2k)+1 */
mpz_invert (t, M, N);
mpz_add (M, M, t); /* M = x^k + x^(-k) */
}
else
{
mpz_set (g[0], M);
mpz_set_si (g[1], -1);
}
/* Test polynomials (check that values vanish (mod N) at root M) */
mpz_set_ui (t, 0);
for (i = degree; i >= 0; i--)
{
mpz_mul (t, t, M);
mpz_add (t, t, f[i]);
mpz_mod (t, t, N);
}
if (mpz_sgn (t) != 0)
{
gmp_fprintf (stderr, "Error: M=%Zd is not a root of f(x) % N\n", M);
fprintf (stderr, "f(x) = ");
for (i = degree; i >= 0; i--)
gmp_fprintf (stderr, "%Zd*x^%d +", f[i], i);
gmp_fprintf (stderr, "\n""Remainder is %Zd\n", t);
gmp_fprintf (stderr, "%Zd^%d%c1\n", x, n, sign < 0 ? '-' : '+');
fprintf (stderr, "Please report this bug.\n");
exit (EXIT_FAILURE);
}
mpz_mul (t, g[1], M);
mpz_add (t, t, g[0]);
mpz_mod (t, t, N);
if (mpz_sgn (t) != 0)
{
gmp_fprintf (stderr, "Error: M=%Zd is not a root of g(x) % N\n", M);
gmp_fprintf (stderr, "Remainder is %Zd\n", t);
fprintf (stderr, "Please report this bug.\n");
exit (EXIT_FAILURE);
}
/* Output polynomials */
if (format == FRANKE)
{
const double cost = snfscost (difficulty + alpha / log (10.) * degree);
gmp_printf ("%Zd\n", N);
gmp_printf ("# %Zd^%d%c1, difficulty: %.2f, skewness: %.2f, alpha: %.2f\n"
"# Cost: %g, est. time: %.2f GHz days (not accurate yet!)\n",
x, n, sign == 1 ? '+' : '-', difficulty, skewness, alpha,
cost, cost * timefactor);
for (i = degree; i >= 0; i--)
if (!shortform || mpz_sgn (f[i]) != 0)
gmp_printf ("X%d %Zd\n", i, f[i]);
if (!shortform || halved) /* Franke's siever uses g(x)=-x+M by default */
{
gmp_printf ("Y1 %Zd\n", g[1]);
gmp_printf ("Y0 %Zd\n", g[0]);
}
gmp_printf ("M %Zd\n", M);
}
else if (format == GGNFS)
{
const double cost = snfscost (difficulty);
gmp_printf ("n: %Zd\n", N);
gmp_printf ("# %Zd^%d%c1, difficulty: %.2f, skewness: %.2f, alpha: %.2f\n"
"# cost: %g, est. time: %.2f GHz days (not accurate yet!)\n",
x, n, sign == 1 ? '+' : '-', difficulty, skewness, alpha,
cost, cost * timefactor);
printf ("skew: %.3f\n", skewness);
for (i = degree; i >= 0; i--)
if (mpz_sgn (f[i]) != 0)
gmp_printf ("c%d: %Zd\n", i, f[i]);
gmp_printf ("Y1: %Zd\n", g[1]);
gmp_printf ("Y0: %Zd\n", g[0]);
gmp_printf ("m: %Zd\n", M);
printf ("type: snfs\n");
}
else if (format == MSIEVE)
{
gmp_printf ("N %Zd\n", N);
gmp_printf ("R0 %Zd\n", g[0]);
gmp_printf ("R1 %Zd\n", g[1]);
for (i = 0; i <= degree; i++)
gmp_printf ("A%d %Zd\n", i, f[i]);
}
else
{
gmp_printf ("%Zd\n", N);
gmp_printf ("%Zd\n\n2\n\n1\n", M);
gmp_printf ("%Zd %Zd\n", g[0], g[1]);
printf ("\n\n%d\n", degree);
for (i = 0; i <= degree; i++)
gmp_printf ("%Zd ", f[i]);
printf ("\n");
}
return 0;
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...