| Gerhard Niklasch on Mon, 3 Aug 1998 19:13:26 +0200 (MET DST) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| MPQS patch for 2.0.11.beta |
Dear beta testers and other friends of pari-2.0.*,
Please find appended a medium-long diff for src/modules/mpqs.c which fixes
one problem of my own making in pari-2.0.11.beta (we were leaking at least
one file descriptor per MPQS factorization), and adds some other touches,
some of which make it run more efficiently (especially as far as combining
Full Relations from Large Prime Relations is concerned) whilst others merely
improve, extend or correct the diagnostics.
If you don't want to patch, you can fetch the entire gzipped mpqs.c file
from
ftp://hasse.mathematik.tu-muenchen.de/incoming/mpqs.c.gz
(beware that this _might_ get replaced at any time by yet another updated
version -- I hope not before the 2.0.x.RELEASE, but Murphy's Law has a way
of sneaking gremlins into such code... and it will be removed when 2.0.12
comes out). NB you cannot list the incoming directory, of course.
Karim: This is the definitive one; when your mail system comes up again,
you can ignore the incremental patches to this file from Friday, Saturday
and earlier today (except perhaps for the detailed list of changes).
Thanks to all testers and happy factoring,
Enjoy, Gerhard
bash$ diff -u src/modules/mpqs.c.orig src/modules/mpqs.c
--- src/modules/mpqs.c.orig Thu Jul 30 00:30:05 1998
+++ src/modules/mpqs.c Mon Aug 3 18:30:38 1998
@@ -125,6 +125,12 @@
FILE *fp = fopen(s, mode);
if (fp == NULL)
err(talker, "MPQS: could not open requested file %s", s);
+ else if (
+#ifdef MPQS_DEBUG
+ 1 ||
+#endif
+ DEBUGLEVEL >= 6)
+ fprintferr("MPQS: opening file %s (mode %s)\n", s, mode);
return fp;
}
@@ -133,10 +139,16 @@
**/
static void
-mpqs_close_file(FILE *fp)
+mpqs_close_file(FILE *fp, char *s)
{
if (fclose(fp))
- err(warner, "error trying to close temporary file");
+ err(warner, "error trying to close temporary file %s", s);
+ else if (
+#ifdef MPQS_DEBUG
+ 1 ||
+#endif
+ DEBUGLEVEL >= 6)
+ fprintferr("MPQS: closed file %s\n", s);
}
static void
@@ -144,6 +156,12 @@
{
if (unlink(f))
err(warner, "can\'t remove file %s", f);
+ else if (
+#ifdef MPQS_DEBUG
+ 1 ||
+#endif
+ DEBUGLEVEL >= 6)
+ fprintferr("MPQS: removed file %s\n", f);
}
/**
@@ -156,8 +174,10 @@
{
char **sa = (char **) a;
char **sb = (char **) b;
- long qa = atol(*sa);
- long qb = atol(*sb);
+ long qa = strtol(*sa, NULL, 10);
+ long qb = strtol(*sb, NULL, 10);
+ /* atol() isn't entirely portable for the Full Relations case where the
+ strings of digits are too long to fit into a long --GN */
if (qa < qb) return -1;
else if (qa > qb) return 1;
else return strcmp(*sa, *sb);
@@ -213,7 +233,7 @@
{ /* file empty */
avma = av;
free(buf);
- mpqs_close_file(TMP);
+ mpqs_close_file(TMP, filename);
return 0;
}
/* enter first buffer into buflist */
@@ -246,7 +266,7 @@
#ifdef MPQS_DEBUG
1 ||
#endif
- DEBUGLEVEL >= 6)
+ DEBUGLEVEL >= 7)
fprintferr("MQPS: short of space -- another buffer for sorting\n");
buf = (char *) gpmalloc(MPQS_STRING_LENGTH * sizeof(char));
cur_line = buf;
@@ -288,7 +308,7 @@
#ifdef MPQS_DEBUG
1 ||
#endif
- DEBUGLEVEL >= 6)
+ DEBUGLEVEL >= 7)
fprintferr("MQPS: line wrap -- another buffer for sorting\n");
buf = (char *) gpmalloc(MPQS_STRING_LENGTH * sizeof(char));
/* remember buffer for later deallocation */
@@ -310,7 +330,7 @@
/* read remainder of line */
if (fgets(cur_line, bufspace, TMP) == NULL)
{
- mpqs_close_file(TMP);
+ mpqs_close_file(TMP, filename);
err(talker,"MPQS: relations file truncated?!\n");
}
lg1 = strlen(cur_line);
@@ -320,7 +340,7 @@
}
} /* for */
- mpqs_close_file(TMP);
+ mpqs_close_file(TMP, filename);
#if 0 /* caught above, can no longer happen */
if (i == 0)
@@ -348,7 +368,7 @@
}
old_s = sort_table[j];
}
- mpqs_close_file(TMP);
+ mpqs_close_file(TMP, filename);
if (
#ifdef MPQS_DEBUG
1 ||
@@ -401,8 +421,7 @@
}
if (fflush(fp))
err(warner, "error whilst flushing file %s", filename);
- if (fclose(fp))
- err(warner, "error trying to close file %s", filename);
+ mpqs_close_file(fp, filename);
return c;
}
@@ -436,7 +455,7 @@
char line1[MPQS_STRING_LENGTH], line2[MPQS_STRING_LENGTH];
char line[MPQS_STRING_LENGTH];
char *line_new = line1, *line_new_old = line2, *line_tmp;
- long q_new, q_new_old, q, i = 0, c = 0;
+ long q_new, q_new_old = -1, q, i = 0, c = 0;
long comb_in_progress;
/* if LPNEW is empty, copy LPREL to TMP. Could be done by a rename
@@ -459,7 +478,7 @@
if (mode) /* full relations mode */
{
i = mpqs_append_file(TMP, LPNEW, LPTMP_str);
- return i; /* COMB cannot be open here */
+ return i + 1; /* COMB cannot be open here */
}
/* LP mode: check for combinable relations */
@@ -505,7 +524,8 @@
comb_in_progress = 0;
}
} /* while */
- if (COMB) mpqs_close_file(COMB);
+ if (COMB) mpqs_close_file(COMB, COMB_str);
+ mpqs_close_file(TMP, LPTMP_str);
return i;
}
@@ -541,7 +561,7 @@
bummer(LPTMP_str);
if (mode) c++;
i = mpqs_append_file(TMP, LPREL, LPTMP_str);
- if (COMB) mpqs_close_file(COMB);
+ if (COMB) mpqs_close_file(COMB, COMB_str);
return (mode ? c + i : c);
}
q_new = atol(line_new);
@@ -589,8 +609,8 @@
if (fputs(line_new, TMP) < 0)
bummer(LPTMP_str);
i = mpqs_append_file(TMP, LPNEW, LPTMP_str);
- if (COMB) mpqs_close_file(COMB);
- return (mode ? c + i : c);
+ if (COMB) mpqs_close_file(COMB, COMB_str);
+ return (mode ? c + i + 1 : c);
}
else
q = atol(line);
@@ -647,7 +667,7 @@
bummer(LPTMP_str);
if (mode) c++; else c += i;
i = mpqs_append_file(TMP, LPREL, LPTMP_str);
- if (COMB) mpqs_close_file(COMB);
+ if (COMB) mpqs_close_file(COMB, COMB_str);
return (mode ? c + i : c);
}
else
@@ -1030,8 +1050,8 @@
*maxB = 1l << (*fA - 1);
#endif
*si = (ulong)mpqs_parameters[i][5];
- *isoi = (ulong)mpqs_parameters[i][6];
- *fsoi = (ulong)mpqs_parameters[i][7];
+ *isoi = 10 * (ulong)mpqs_parameters[i][6]; /* percentages scaled by 10 */
+ *fsoi = 10 * (ulong)mpqs_parameters[i][7];
}
/******************************/
@@ -1995,169 +2015,169 @@
{
char buf [MPQS_STRING_LENGTH], *s1, *s2;
char new_relation [MPQS_STRING_LENGTH];
- mpqs_lp_entry e[50 + 1];
+ mpqs_lp_entry e[2]; /* we'll use the two alternatingly */
long *ei, ei_size = size_of_FB + 2;
- long old_q = -1;
- long i, j, k, c = 0;
+ long old_q;
+ GEN inv_q, Y1, Y2, new_Y, new_Y1;
+ long i, l, c = 0;
+ long av = avma, av2;
*f = NULL;
- if (fgets(buf, MPQS_STRING_LENGTH, COMB) != NULL)
- {
- ei = (long *) gpmalloc(ei_size * sizeof(long));
- s1 = buf; s2 = strchr(s1, ' '); *s2 = '\0';
- e[0].q = atol(s1);
- s1 = s2 + 3; s2 = strchr(s1, ' '); *s2 = '\0';
- strcpy(e[0].Y, s1);
- s1 = s2 + 3; s2 = strchr(s1, '\n'); *s2 = '\0';
- strcpy(e[0].ep, s1);
+ if (fgets(buf, MPQS_STRING_LENGTH, COMB) == NULL)
+ return 0; /* nothing there -- should not happen */
- i = 1;
- old_q = e[0].q;
+ /* put first lp relation in row 0 of e */
+ s1 = buf; s2 = strchr(s1, ' '); *s2 = '\0';
+ e[0].q = atol(s1);
+ s1 = s2 + 3; s2 = strchr(s1, ' '); *s2 = '\0';
+ strcpy(e[0].Y, s1);
+ s1 = s2 + 3; s2 = strchr(s1, '\n'); *s2 = '\0';
+ strcpy(e[0].ep, s1);
- while (fgets(buf, MPQS_STRING_LENGTH, COMB) != NULL)
- {
+ i = 1; /* second relation will go into row 1 */
+ old_q = e[0].q;
+ while (!inversemodulo(stoi(old_q), kN, &inv_q)) /* can happen --GN */
+ {
+ inv_q = mppgcd(inv_q, N);
+ if (is_pm1(inv_q) || egalii(inv_q, N)) /* pity */
+ {
+#ifdef MPQS_DEBUG
+ fprintferr("MPQS: skipping relation with non-invertible q\n");
+#endif
+ avma = av;
+ if (fgets(buf, MPQS_STRING_LENGTH, COMB) == NULL)
+ return 0;
s1 = buf; s2 = strchr(s1, ' '); *s2 = '\0';
- e[i].q = atol(s1);
+ e[0].q = atol(s1);
s1 = s2 + 3; s2 = strchr(s1, ' '); *s2 = '\0';
- strcpy(e[i].Y, s1);
+ strcpy(e[0].Y, s1);
s1 = s2 + 3; s2 = strchr(s1, '\n'); *s2 = '\0';
- strcpy(e[i].ep, s1);
+ strcpy(e[0].ep, s1);
+ old_q = e[0].q;
+ }
+ *f = gerepileupto(av, inv_q);
+ return c;
+ }
+ Y1 = lisexpr(e[0].Y);
+ av2 = avma; /* preserve inv_q and Y1 */
- if (e[i].q == old_q)
- {
- /* we only handle at most 50 large prime relations
- with the same q */
- if (i < 50)
- i++;
- }
- else
+ ei = (long *) gpmalloc(ei_size * sizeof(long));
+
+ while (fgets(buf, MPQS_STRING_LENGTH, COMB) != NULL)
+ {
+ s1 = buf; s2 = strchr(s1, ' '); *s2 = '\0';
+ e[i].q = atol(s1);
+ s1 = s2 + 3; s2 = strchr(s1, ' '); *s2 = '\0';
+ strcpy(e[i].Y, s1);
+ s1 = s2 + 3; s2 = strchr(s1, '\n'); *s2 = '\0';
+ strcpy(e[i].ep, s1);
+
+ if (e[i].q != old_q)
+ {
+ /* switch to combining a new bunch, swapping the rows */
+ old_q = e[i].q;
+ avma = av; /* discard old inv_q and Y1 */
+ if (!inversemodulo(stoi(old_q), kN, &inv_q)) /* can happen --GN */
{
- if (i == 1)
- {
- old_q = e[1].q;
- e[0].q = e[1].q;
- strcpy(e[0].Y, e[1].Y);
- strcpy(e[0].ep, e[1].ep);
- }
- else
+ inv_q = mppgcd(inv_q, N);
+ if (is_pm1(inv_q) || egalii(inv_q, N)) /* pity */
{
- long av = avma;
- GEN inv_q;
- if (!inversemodulo(stoi(e[0].q), kN, &inv_q)) /* can happen --GN */
- {
- inv_q = mppgcd(inv_q, N);
- if (is_pm1(inv_q) || egalii(inv_q, N)) /* pity */
- {
#ifdef MPQS_DEBUG
- fprintferr("MPQS: skipping relation with non-invertible q\n");
+ fprintferr("MPQS: skipping relation with non-invertible q\n");
#endif
- /* reset table */
- old_q = e[i].q;
- e[0].q = e[i].q;
- strcpy(e[0].Y, e[i].Y);
- strcpy(e[0].ep, e[i].ep);
- i = 1;
- avma = av;
- continue; /* discard this combination */
- }
- *f = gerepileupto(av, inv_q);
- free (ei);
- return c;
- }
- for (j = 0; j < i; j++)
- {
- for(k = j + 1; k < i; k++)
- {
- long av2 = avma, l;
- GEN Y1, Y2, new_Y, new_Y1;
- c++;
-
- /* compute relation [j, k] */
- mpqs_combine_exponents(ei, ei_size, e[j].ep, e[k].ep);
- if (DEBUGLEVEL >= 6)
- {
- fprintf(stderr, "MPQS: combining\n");
- fprintf(stderr, " [%ld] {%ld @ %s : %s}\n",
- j, e[j].q, e[j].Y, e[j].ep);
- fprintf(stderr, " [%ld] {%ld @ %s : %s}\n",
- k, e[k].q, e[k].Y, e[k].ep);
- }
- Y1 = lisexpr(e[j].Y);
- Y2 = lisexpr(e[k].Y);
- new_Y = modii(mulii(mulii(Y1, Y2), inv_q), kN);
- new_Y1 = subii(kN, new_Y);
- if (absi_cmp(new_Y1, new_Y) < 0) new_Y = new_Y1;
- s1 = gen2str(new_Y);
- strcpy(new_relation, s1);
- strcat(new_relation, " :");
- if (ei[1] % 2)
- strcat(new_relation, " 1 1");
- for (l = 2; l < ei_size; l++)
- {
- if (ei[l])
- {
- sprintf(buf, " %ld %ld", ei[l], l);
- strcat(new_relation, buf);
- }
- }
- strcat(new_relation, " 0");
- if (DEBUGLEVEL >= 6)
- fprintf(stderr, " == {%s}\n", new_relation);
- strcat(new_relation, "\n");
-
-#ifdef MPQS_DEBUG
- {
- char ejk [MPQS_STRING_LENGTH];
- GEN Qx_2, prod_pi_ei, pi_ei;
- char *s;
- long pi, ei, av;
-
- av = avma;
- Qx_2 = modii(sqri(new_Y), kN);
-
- strcpy(ejk, new_relation);
- s = strchr(ejk, ':') + 2;
- s = strtok(s, " \n");
-
- prod_pi_ei = gun;
- while (s != NULL)
- {
- ei = atol(s);
- if (ei == 0) break;
- s = strtok(NULL, " \n");
- pi = atol(s);
- pi_ei = powmodulo(stoi(FB[pi]), stoi(ei), kN);
- prod_pi_ei = modii(mulii(prod_pi_ei, pi_ei), kN);
- s = strtok(NULL, " \n");
- }
- avma = av;
-
- if (gcmp(Qx_2, prod_pi_ei) != 0)
- err(talker,
- "MPQS: combined large prime relation is false!!");
- }
-#endif
-
- if (fputs(new_relation, FNEW) < 0)
- bummer(FNEW_str);
- free(s1);
- avma = av2;
- }
- }
-
- /* reset table */
- old_q = e[i].q;
- e[0].q = e[i].q;
- strcpy(e[0].Y, e[i].Y);
- strcpy(e[0].ep, e[i].ep);
- i = 1;
avma = av;
+ old_q = -1; /* sentinel */
+ av2 = avma;
+ continue; /* discard this combination */
}
+ *f = gerepileupto(av, inv_q);
+ free (ei);
+ return c;
+ }
+ Y1 = lisexpr(e[i].Y);
+ i = 1 - i; /* subsequent relations go to other row */
+ av2 = avma; /* preserve inv_q and Y1 */
+ continue;
+ }
+ /* count and combine the two we've got, and continue in the same row */
+ c++;
+
+ mpqs_combine_exponents(ei, ei_size, e[1-i].ep, e[i].ep);
+ if (DEBUGLEVEL >= 6)
+ {
+ fprintf(stderr, "MPQS: combining\n");
+ fprintf(stderr, " {%ld @ %s : %s}\n",
+ old_q, e[1-i].Y, e[1-i].ep);
+ fprintf(stderr, " * {%ld @ %s : %s}\n",
+ e[i].q, e[i].Y, e[i].ep);
+ }
+ Y2 = lisexpr(e[i].Y);
+ new_Y = modii(mulii(mulii(Y1, Y2), inv_q), kN);
+ new_Y1 = subii(kN, new_Y);
+ if (absi_cmp(new_Y1, new_Y) < 0) new_Y = new_Y1;
+ s1 = gen2str(new_Y);
+ strcpy(new_relation, s1);
+ strcat(new_relation, " :");
+ if (ei[1] % 2)
+ strcat(new_relation, " 1 1");
+ for (l = 2; l < ei_size; l++)
+ {
+ if (ei[l])
+ {
+ sprintf(buf, " %ld %ld", ei[l], l);
+ strcat(new_relation, buf);
}
- } /* while */
- free(ei);
- }
+ }
+ strcat(new_relation, " 0");
+ if (DEBUGLEVEL >= 6)
+ fprintf(stderr, " == {%s}\n", new_relation);
+ strcat(new_relation, "\n");
+
+#ifdef MPQS_DEBUG
+ {
+ char ejk [MPQS_STRING_LENGTH];
+ GEN Qx_2, prod_pi_ei, pi_ei;
+ char *s;
+ long pi, exi, av1 = avma;
+ Qx_2 = modii(sqri(new_Y), kN);
+
+ strcpy(ejk, new_relation);
+ s = strchr(ejk, ':') + 2;
+ s = strtok(s, " \n");
+
+ prod_pi_ei = gun;
+ while (s != NULL)
+ {
+ exi = atol(s);
+ if (exi == 0) break;
+ s = strtok(NULL, " \n");
+ pi = atol(s);
+ pi_ei = powmodulo(stoi(FB[pi]), stoi(exi), kN);
+ prod_pi_ei = modii(mulii(prod_pi_ei, pi_ei), kN);
+ s = strtok(NULL, " \n");
+ }
+ avma = av1;
+
+ if (!egalii(Qx_2, prod_pi_ei))
+ {
+ free(ei);
+ err(talker,
+ "MPQS: combined large prime relation is false");
+ }
+ }
+#endif
+
+ if (fputs(new_relation, FNEW) < 0)
+ {
+ free(ei);
+ bummer(FNEW_str);
+ }
+ free(s1);
+ avma = av2;
+ } /* while */
+
+ free(ei);
if (DEBUGLEVEL >= 4)
fprintferr("MPQS: combined %ld full relation%s\n",
@@ -2322,19 +2342,24 @@
}
}
-/* read and create an mpqs_gauss_matrix from a relation file FREL
- rows = size_of_FB, cols = rel */
+/* create and read an mpqs_gauss_matrix from a relation file FREL (opened
+ by caller). Also record the position of each relation in the file for
+ later use
+ rows = size_of_FB+1, cols = rel */
static mpqs_gauss_matrix
-mpqs_gauss_read_matrix(long rows, long cols)
+mpqs_gauss_read_matrix(FILE *FREL, long rows, long cols, long *fpos)
{
mpqs_gauss_matrix m;
long i = 0, e, p;
char buf[MPQS_STRING_LENGTH], *s;
- FILE *fp;
m = mpqs_gauss_create_matrix(rows, cols);
- fp = mpqs_open_file(FREL_str, "r");
- while (fgets(buf, MPQS_STRING_LENGTH, fp) != NULL)
+ if ((fpos[0] = ftell(FREL)) < 0)
+ {
+ mpqs_close_file(FREL, FREL_str);
+ err(talker, "ftell error on full relations file");
+ }
+ while (fgets(buf, MPQS_STRING_LENGTH, FREL) != NULL)
{
s = strchr(buf, ':') + 2;
s = strtok(s, " \n");
@@ -2350,8 +2375,19 @@
s = strtok(NULL, " \n");
}
i++;
+ if (i < cols && (fpos[i] = ftell(FREL)) < 0)
+ {
+ mpqs_close_file(FREL, FREL_str);
+ err(talker, "ftell error on full relations file");
+ }
+ }
+ if (i != cols)
+ {
+ mpqs_close_file(FREL, FREL_str);
+ fprintferr("MPQS: full relations file %s than expected",
+ i > cols ? "longer" : "shorter");
+ err(talker, "MPQS panicking");
}
- mpqs_close_file(fp);
return m;
}
@@ -2460,19 +2496,19 @@
}
static char*
-mpqs_get_relation(long pos, FILE *fp)
+mpqs_get_relation(long pos, FILE *FREL)
{
static char buf [MPQS_STRING_LENGTH];
- if (fseek(fp, pos, SEEK_SET))
+ if (fseek(FREL, pos, SEEK_SET))
{
- mpqs_close_file(fp);
+ mpqs_close_file(FREL, FREL_str); /* it's always FREL */
/* may also fail, but we're giving up anyway... GN */
- err(talker, "can\'t seek relations file");
+ err(talker, "can\'t seek full relations file");
}
- if (fgets(buf, MPQS_STRING_LENGTH, fp) == NULL)
+ if (fgets(buf, MPQS_STRING_LENGTH, FREL) == NULL)
{
- mpqs_close_file(fp);
- err(talker, "relations file truncated?!");
+ mpqs_close_file(FREL, FREL_str);
+ err(talker, "full relations file truncated?!");
}
return buf;
}
@@ -2486,9 +2522,8 @@
static GEN
mpqs_solve_linear_system(GEN kN, GEN N, long rel, long *FB, long size_of_FB)
{
- FILE *fp;
+ FILE *FREL;
GEN X, Y_prod, N_or_kN, D1, base, res, new_res;
- char buf [MPQS_STRING_LENGTH];
long av=avma, av2, av3, lim, tetpil;
long *fpos, *ei;
long i, j, H_cols, H_rows;
@@ -2501,27 +2536,10 @@
#else
N_or_kN = (DEBUGLEVEL >= 4 ? kN : N);
#endif /* --GN */
- fp = mpqs_open_file(FREL_str, "r");
+ FREL = mpqs_open_file(FREL_str, "r");
fpos = (long *) gpmalloc(rel * sizeof(long));
- for (i = 0; i < rel; i++)
- {
- if ((fpos[i] = ftell(fp)) < 0)
- {
- mpqs_close_file(fp);
- err(talker, "ftell error on full relations file");
- }
- if (fgets(buf, MPQS_STRING_LENGTH, fp) == NULL)
- {
- mpqs_close_file(fp);
- err(talker, "full relations file truncated?!");
- }
- }
- mpqs_close_file(fp);
- /* somewhat inefficient... we could rewind it, and m'gauss_read_matrix()
- * could start and end with the file already open (and then we'll rewind
- * it again, or some such thing) --GN */
- m = mpqs_gauss_read_matrix(size_of_FB+1, rel);
+ m = mpqs_gauss_read_matrix(FREL, size_of_FB+1, rel, fpos);
if (DEBUGLEVEL >= 7)
{
fprintferr("\\\\ MATRIX READ BY MPQS\nFREL=");
@@ -2557,6 +2575,7 @@
mpqs_gauss_destroy_matrix(m, size_of_FB+1, rel);
mpqs_gauss_destroy_matrix(ker_m, rel, rank);
free(fpos);
+ mpqs_close_file(FREL, FREL_str);
avma = av;
return NULL; /* no factors found */
}
@@ -2564,7 +2583,7 @@
/* If the rank is r, we can expect up to 2^r pairwise coprime factors,
but it may well happen that a kernel basis vector contributes nothing
new to the decomposition. We will allocate room for max(2,r) factors
- initially (certainly adequare when one or two basis vectors work),
+ initially (certainly adequate when one or two basis vectors work),
adjusting this down at the end to what we actually found, or up if
we are very lucky and find more factors. In the upper half of our
vector, we store information about which factors we know to be composite
@@ -2582,43 +2601,25 @@
for (i=2*res_size; i; i--) res[i] = 0; /* make it safe for gcopy() */
res_next = res_last = 1;
- /* re-opens FREL file a third time for reading (see above) --GN */
- fp = mpqs_open_file(FREL_str, "r");
- /* Recomputation of N from kN removed --GN */
-
- /* use PARI idiom -- create results on the stack; only clean up when we
- absolutely have to (not very likely; numbers treated by mpqs() are
- short after all, compared to a typical PARI stack) --GN */
-#if 0
- Y_prod = cgeti(lg(N_or_kN) + 1);
- X = cgeti(lg(N_or_kN) + 1);
-#endif
ei = (long *) gpmalloc((size_of_FB + 2) * sizeof(long));
for (i = 0; i < H_cols; i++) /* loop over basis of kernel */
{
-#if 0
- affii(gun, Y_prod);
-#else
X = Y_prod = gun;
-#endif
memset(ei, 0, (size_of_FB + 2) * sizeof(long));
for (j = 0; j < H_rows; j++)
{
if (mpqs_gauss_get_bit(ker_m, j, i))
Y_prod = mpqs_add_relation(Y_prod, N_or_kN, ei,
- mpqs_get_relation(fpos[j], fp));
+ mpqs_get_relation(fpos[j], FREL));
}
-#if 0
- affii(gun, X);
-#endif
for (j = 2; j <= (size_of_FB + 1); j++)
{
if (ei[j] & 1)
{ /* this cannot happen --GN */
- mpqs_close_file(fp);
+ mpqs_close_file(FREL, FREL_str);
mpqs_gauss_destroy_matrix(m, size_of_FB+1, rel);
mpqs_gauss_destroy_matrix(ker_m, rel, rank);
free(ei); free(fpos);
@@ -2639,12 +2640,6 @@
{ /* this shouldn't happen either --GN */
if (signe(modii(subii(sqri(X), sqri(Y_prod)), N_or_kN)))
{
-#if 0
- mpqs_close_file(fp);
- mpqs_gauss_destroy_matrix(m, size_of_FB+1, rel);
- mpqs_gauss_destroy_matrix(ker_m, rel, rank);
- free(ei); free(fpos);
-#endif
fprintferr("MPQS: X^2 - Y^2 != 0 mod %s\n",
(N_or_kN == kN ? "kN" : "N"));
fprintferr("\tindex i = %ld\n", i);
@@ -2852,7 +2847,7 @@
if (res_next < 3) /* still no factors seen */
{
- mpqs_close_file(fp);
+ mpqs_close_file(FREL, FREL_str);
mpqs_gauss_destroy_matrix(m, size_of_FB+1, rel);
mpqs_gauss_destroy_matrix(ker_m, rel, rank);
free(ei); free(fpos);
@@ -2862,7 +2857,7 @@
/* normal case: convert internal format to ifac format as described in
src/basemath/ifactor1.c (triples of components: value, exponent, class
[unknown, known composite, known prime]), clean up and return result */
- mpqs_close_file(fp);
+ mpqs_close_file(FREL, FREL_str);
mpqs_gauss_destroy_matrix(m, size_of_FB+1, rel);
mpqs_gauss_destroy_matrix(ker_m, rel, rank);
free(ei); free(fpos);
@@ -2897,7 +2892,10 @@
/******************************/
-#define max_percentage 150 /* give up when nothing found after ~1.5
+/* All percentages below are actually fixed-point quantities scaled by 10
+ (value of 1 means 0.1%, 1000 means 100%). --GN */
+
+#define max_percentage 1500 /* give up when nothing found after ~1.5
times the required number of relations has
been computed (n might be a prime power
or the parameters might be exceptionally
@@ -2911,6 +2909,11 @@
** seem to make any headway.
**/
+/* TODO: this function to be renamed mpqs_main() with several extra para-
+ meters, with mpqs() as a wrapper for the standard case, so we can do
+ partial runs distributed across several machines etc. (not necessarily
+ from gp, but certainly from a small dedicated C program). --GN */
+
GEN
mpqs(GEN N)
{
@@ -2975,7 +2978,7 @@
probably factor over the FB */
long lp_bound; /* size limit for large primes */
ulong sort_interval; /* these determine when to sort and merge */
- ulong followup_sort_interval; /* the temporary files */
+ ulong followup_sort_interval; /* the temporary files (scaled percentages) */
long total_full_relations = 0;
long total_partial_relations = 0;
FILE *FREL;
@@ -2983,9 +2986,9 @@
FILE *LPREL;
FILE *LPNEW;
FILE *COMB;
- long percentage = 0;
+ long percentage = 0; /* scaled by 10, see comment above */
long total_candidates_number = 0;
- long useless_candidates = 0;
+ long useless_candidates = 0; /* another scaled percentage */
long vain_iterations = 0;
long good_iterations = 0;
long iterations = 0;
@@ -3053,9 +3056,9 @@
size_of_FB);
fprintferr("MPQS: striving for %ld relations\n",
FB_overshoot);
- fprintferr("MPQS: first sorting at %ld%%, then every %ld%% / %ld%%\n",
- sort_interval, followup_sort_interval,
- followup_sort_interval>>1);
+ fprintferr("MPQS: first sorting at %ld%%, then every %3.1f%% / %3.1f%%\n",
+ sort_interval/10, followup_sort_interval/10.,
+ followup_sort_interval/20.);
fprintferr("MPQS: initial sieving index = %ld\n",
starting_sieving_index);
}
@@ -3263,9 +3266,9 @@
FREL = mpqs_open_file(FREL_str, "w");
/* This was just for truncating */
- mpqs_close_file(FREL);
+ mpqs_close_file(FREL, FREL_str);
LPREL = mpqs_open_file(LPREL_str, "w");
- mpqs_close_file(LPREL);
+ mpqs_close_file(LPREL, LPREL_str);
LPNEW = mpqs_open_file(LPNEW_str, "w");
FNEW = mpqs_open_file(FNEW_str, "w");
@@ -3291,8 +3294,9 @@
{
/* TODO: increase some parameters. For the moment, simply give up */
CLEANUP();
- mpqs_close_file(FNEW);
- mpqs_close_file(LPNEW); /* FREL, LPREL are closed at this point */
+ mpqs_close_file(FNEW, FNEW_str);
+ mpqs_close_file(LPNEW, LPNEW_str);
+ /* FREL, LPREL are closed at this point */
mpqs_unlink(FREL_str);
mpqs_unlink(FNEW_str);
mpqs_unlink(LPREL_str);
@@ -3306,8 +3310,8 @@
if (is_pm1(fact) || egalii(fact,N))
continue; /* cannot use the current polynomial */
CLEANUP();
- mpqs_close_file(FNEW);
- mpqs_close_file(LPNEW);
+ mpqs_close_file(FNEW, FNEW_str);
+ mpqs_close_file(LPNEW, LPNEW_str);
mpqs_unlink(FREL_str);
mpqs_unlink(FNEW_str);
mpqs_unlink(LPREL_str);
@@ -3366,9 +3370,9 @@
total_full_relations += t;
percentage =
- (long)((100.0 * total_full_relations) / FB_overshoot);
+ (long)((1000.0 * total_full_relations) / FB_overshoot);
useless_candidates =
- (long)((100.0 * (total_candidates_number - total_full_relations))
+ (long)((1000.0 * (total_candidates_number - total_full_relations))
/ (total_candidates_number ? total_candidates_number : 1));
if (percentage < sort_interval)
@@ -3378,25 +3382,32 @@
if (DEBUGLEVEL >= 3)
{
if (DEBUGLEVEL >= 4)
- fprintferr("\nMPQS: passing the %ld%% checkpoint, time = %ld ms\n",
- sort_interval, timer2());
+ fprintferr("\nMPQS: passing the %3.1f%% checkpoint, time = %ld ms\n",
+ sort_interval/10., timer2());
else
- fprintferr("\nMPQS: passing the %ld%% mark\n", sort_interval);
+ fprintferr("\nMPQS: passing the %3.1f%% checkpoint\n",
+ sort_interval/10.);
flusherr();
}
/* sort LPNEW and merge it into LPREL, diverting combinables into COMB */
- mpqs_close_file(LPNEW);
+ mpqs_close_file(LPNEW, LPNEW_str);
mpqs_sort_lp_file(LPNEW_str);
LPREL = mpqs_open_file(LPREL_str, "r");
LPNEW = mpqs_open_file(LPNEW_str, "r");
tp = mpqs_mergesort_lp_file(LPREL, LPNEW, 0);
- mpqs_close_file(LPREL);
- mpqs_close_file(LPNEW);
+ mpqs_close_file(LPREL, LPREL_str);
+ mpqs_close_file(LPNEW, LPNEW_str);
mpqs_unlink(LPREL_str);
if (rename(LPTMP_str, LPREL_str))
err(talker, "can\'t rename file %s to %s", LPTMP_str, LPREL_str);
+ else if (
+#ifdef MPQS_DEBUG
+ 1 ||
+#endif
+ DEBUGLEVEL >= 6)
+ fprintferr("MPQS: renamed file %s to %s\n", LPTMP_str, LPREL_str);
LPNEW = mpqs_open_file(LPNEW_str, "w");
/* combine whatever there is to be combined */
@@ -3409,15 +3420,15 @@
, FB
#endif
);
- mpqs_close_file(COMB);
+ mpqs_close_file(COMB, COMB_str);
mpqs_unlink(COMB_str);
/* now FREL, LPREL are closed and FNEW, LPNEW are still open */
if (fact) /* factor found during combining */
{
/* clean up */
CLEANUP();
- mpqs_close_file(LPNEW);
- mpqs_close_file(FNEW);
+ mpqs_close_file(LPNEW, LPNEW_str);
+ mpqs_close_file(FNEW, FNEW_str);
mpqs_unlink(FREL_str);
mpqs_unlink(FNEW_str);
mpqs_unlink(LPREL_str);
@@ -3436,7 +3447,7 @@
}
/* sort FNEW and merge it into FREL */
- mpqs_close_file(FNEW);
+ mpqs_close_file(FNEW, FNEW_str);
mpqs_sort_lp_file(FNEW_str);
FREL = mpqs_open_file(FREL_str, "r");
@@ -3444,12 +3455,18 @@
total_full_relations = mpqs_mergesort_lp_file(FREL, FNEW, 1);
/* this being the definitive count (combinables combined, and
duplicates removed) */
- mpqs_close_file(FREL);
- mpqs_close_file(FNEW);
+ mpqs_close_file(FREL, FREL_str);
+ mpqs_close_file(FNEW, FNEW_str);
mpqs_unlink(FREL_str);
if (rename(LPTMP_str, FREL_str))
err(talker, "can\'t rename file %s to %s", LPTMP_str, FREL_str);
+ else if (
+#ifdef MPQS_DEBUG
+ 1 ||
+#endif
+ DEBUGLEVEL >= 6)
+ fprintferr("MPQS: renamed file %s to %s\n", LPTMP_str, FREL_str);
/* FNEW stays closed until we know whether we need to reopen it for
another iteration */
@@ -3460,21 +3477,22 @@
is nothing to worry about, since we _are_ making progress neverthe-
less. --GN */
percentage =
- (long)((100.0 * total_full_relations) / FB_overshoot);
+ (long)((1000.0 * total_full_relations) / FB_overshoot);
vain_iterations =
- (long)((100.0 * (iterations - good_iterations))
+ (long)((1000.0 * (iterations - good_iterations))
/ iterations);
- if (percentage >= 84)
+ /* TODO: the following could be improved (extrapolate from past
+ experience how many combined full relations we can expect in
+ addition to those we're finding directly) --GN */
+ if (percentage >= 840)
{
- if (percentage >= 98)
- sort_interval = percentage + 1;
+ if (percentage >= 980)
+ sort_interval = percentage + 10;
else
- {
sort_interval = percentage + (followup_sort_interval >> 1);
- if (sort_interval > 100 && percentage < 100)
- sort_interval = 100;
- }
+ if (sort_interval >= 1000 && percentage < 1000)
+ sort_interval = 1000;
}
else
{
@@ -3486,8 +3504,8 @@
fprintferr("MPQS: done sorting%s, time = %ld ms\n",
tp > 0 ? " and combining" : "",
timer2());
- fprintferr("MPQS: found %ld%% of the required relations\n",
- percentage);
+ fprintferr("MPQS: found %3.1f%% of the required relations\n",
+ percentage/10.);
if (DEBUGLEVEL >= 5)
{
/* total_full_relations are always plural */
@@ -3495,23 +3513,23 @@
total_full_relations);
fprintferr("MPQS: (%ld of these from partial relations)\n",
total_partial_relations);
- fprintferr("MPQS: %ld%% useless candidates\n",
- useless_candidates);
- fprintferr("MPQS: %ld%% of the iterations yielded no candidates\n",
- vain_iterations);
- fprintferr("MPQS: next checkpoint at %ld%%\n",
- sort_interval);
+ fprintferr("MPQS: %4.1f%% useless candidates\n",
+ useless_candidates/10.);
+ fprintferr("MPQS: %4.1f%% of the iterations yielded no candidates\n",
+ vain_iterations/10.);
+ fprintferr("MPQS: next checkpoint at %3.1f%%\n",
+ sort_interval/10.);
}
}
- if (percentage < 100)
+ if (percentage < 1000)
{
FNEW = mpqs_open_file(FNEW_str, "w");
/* at this point, LPNEW and FNEW are again open for writing */
continue; /* main loop */
}
- /* percentage >= 100, which implies total_full_relations > size_of_FB:
+ /* percentage >= 1000, which implies total_full_relations > size_of_FB:
try finishing it off */
/* solve the system over F_2 */
@@ -3527,7 +3545,7 @@
{
/* clean up */
CLEANUP();
- mpqs_close_file(LPNEW);
+ mpqs_close_file(LPNEW, LPNEW_str);
mpqs_unlink(FREL_str);
mpqs_unlink(FNEW_str);
mpqs_unlink(LPREL_str);
@@ -3580,7 +3598,7 @@
{
/* clean up */
CLEANUP();
- mpqs_close_file(LPNEW);
+ mpqs_close_file(LPNEW, LPNEW_str);
mpqs_unlink(FREL_str);
mpqs_unlink(FNEW_str);
mpqs_unlink(LPREL_str);