[flintqs] 01/02: Imported Upstream version 20070817
Dominique Ingoglia
godel108-guest at alioth.debian.org
Tue Sep 3 22:12:08 UTC 2013
This is an automated email from the git hooks/post-receive script.
godel108-guest pushed a commit to branch master
in repository flintqs.
commit 472f7ce708608d6a452fffe150431009b108435e
Author: Dominique Ingoglia <Dominique.Ingoglia at gmail.com>
Date: Mon Sep 2 18:31:35 2013 -0600
Imported Upstream version 20070817
---
._F2matrix.cpp | Bin 0 -> 229 bytes
._F2matrix.h | Bin 0 -> 229 bytes
._ModuloArith.cpp | Bin 0 -> 229 bytes
._ModuloArith.h | Bin 0 -> 229 bytes
._QS.cpp | Bin 0 -> 229 bytes
._QS.cpp.orig | Bin 0 -> 229 bytes
._QStodo.txt | Bin 0 -> 229 bytes
._TonelliShanks.cpp | Bin 0 -> 229 bytes
._TonelliShanks.h | Bin 0 -> 229 bytes
._gpl.txt | Bin 0 -> 229 bytes
._lanczos.c | Bin 0 -> 229 bytes
._lanczos.h | Bin 0 -> 229 bytes
._lprels.c | Bin 0 -> 229 bytes
._lprels.h | Bin 0 -> 229 bytes
._makefile | Bin 0 -> 229 bytes
._makefile.opteron | Bin 0 -> 229 bytes
._makefile.sage | Bin 0 -> 229 bytes
F2matrix.cpp | 197 ++++++
F2matrix.h | 46 ++
ModuloArith.cpp | 69 +++
ModuloArith.h | 42 ++
QS.cpp | 1665 +++++++++++++++++++++++++++++++++++++++++++++++++++
QStodo.txt | 50 ++
TonelliShanks.cpp | 147 +++++
TonelliShanks.h | 45 ++
gpl.txt | 339 +++++++++++
lanczos.c | 975 ++++++++++++++++++++++++++++++
lanczos.h | 95 +++
lprels.c | 756 +++++++++++++++++++++++
lprels.h | 51 ++
makefile | 64 ++
makefile.opteron | 60 ++
makefile.osx64 | 60 ++
makefile.sage | 60 ++
34 files changed, 4721 insertions(+)
diff --git a/._F2matrix.cpp b/._F2matrix.cpp
new file mode 100644
index 0000000..81618a0
Binary files /dev/null and b/._F2matrix.cpp differ
diff --git a/._F2matrix.h b/._F2matrix.h
new file mode 100644
index 0000000..7b51d3c
Binary files /dev/null and b/._F2matrix.h differ
diff --git a/._ModuloArith.cpp b/._ModuloArith.cpp
new file mode 100644
index 0000000..cf91239
Binary files /dev/null and b/._ModuloArith.cpp differ
diff --git a/._ModuloArith.h b/._ModuloArith.h
new file mode 100644
index 0000000..ad4f6b6
Binary files /dev/null and b/._ModuloArith.h differ
diff --git a/._QS.cpp b/._QS.cpp
new file mode 100644
index 0000000..775bcde
Binary files /dev/null and b/._QS.cpp differ
diff --git a/._QS.cpp.orig b/._QS.cpp.orig
new file mode 100644
index 0000000..4fff163
Binary files /dev/null and b/._QS.cpp.orig differ
diff --git a/._QStodo.txt b/._QStodo.txt
new file mode 100644
index 0000000..0418b86
Binary files /dev/null and b/._QStodo.txt differ
diff --git a/._TonelliShanks.cpp b/._TonelliShanks.cpp
new file mode 100644
index 0000000..208ee3e
Binary files /dev/null and b/._TonelliShanks.cpp differ
diff --git a/._TonelliShanks.h b/._TonelliShanks.h
new file mode 100644
index 0000000..61ddb95
Binary files /dev/null and b/._TonelliShanks.h differ
diff --git a/._gpl.txt b/._gpl.txt
new file mode 100644
index 0000000..0a90c9e
Binary files /dev/null and b/._gpl.txt differ
diff --git a/._lanczos.c b/._lanczos.c
new file mode 100644
index 0000000..42ca6c0
Binary files /dev/null and b/._lanczos.c differ
diff --git a/._lanczos.h b/._lanczos.h
new file mode 100644
index 0000000..0398ec8
Binary files /dev/null and b/._lanczos.h differ
diff --git a/._lprels.c b/._lprels.c
new file mode 100644
index 0000000..cb4f9e7
Binary files /dev/null and b/._lprels.c differ
diff --git a/._lprels.h b/._lprels.h
new file mode 100644
index 0000000..234ec21
Binary files /dev/null and b/._lprels.h differ
diff --git a/._makefile b/._makefile
new file mode 100644
index 0000000..b612b6b
Binary files /dev/null and b/._makefile differ
diff --git a/._makefile.opteron b/._makefile.opteron
new file mode 100644
index 0000000..1efda60
Binary files /dev/null and b/._makefile.opteron differ
diff --git a/._makefile.sage b/._makefile.sage
new file mode 100644
index 0000000..f796829
Binary files /dev/null and b/._makefile.sage differ
diff --git a/F2matrix.cpp b/F2matrix.cpp
new file mode 100644
index 0000000..be30f08
--- /dev/null
+++ b/F2matrix.cpp
@@ -0,0 +1,197 @@
+/*============================================================================
+ Copyright 2006 William Hart
+
+ This file is part of SIMPQS.
+
+ SIMPQS 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.
+
+ SIMPQS 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 SIMPQS; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+============================================================================*/
+
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+#include "F2matrix.h"
+
+static u_int32_t bitPattern[] =
+{
+ 0x80000000, 0x40000000, 0x20000000, 0x10000000,
+ 0x08000000, 0x04000000, 0x02000000, 0x01000000,
+ 0x00800000, 0x00400000, 0x00200000, 0x00100000,
+ 0x00080000, 0x00040000, 0x00020000, 0x00010000,
+ 0x00008000, 0x00004000, 0x00002000, 0x00001000,
+ 0x00000800, 0x00000400, 0x00000200, 0x00000100,
+ 0x00000080, 0x00000040, 0x00000020, 0x00000010,
+ 0x00000008, 0x00000004, 0x00000002, 0x00000001
+};
+
+void insertEntry(matrix m, u_int32_t i, u_int32_t j)
+{
+ m[i][j / 32] |= bitPattern[j % 32];
+
+ return;
+}
+
+void xorEntry(matrix m, u_int32_t i, u_int32_t j)
+{
+ m[i][j / 32] ^= bitPattern[j % 32];
+
+ return;
+}
+
+u_int32_t getEntry(matrix m, u_int32_t i, u_int32_t j)
+{
+ return m[i][j / 32] & bitPattern[j % 32];
+}
+
+void swapRows(matrix m, u_int32_t x, u_int32_t y)
+{
+ row temp;
+
+ temp = m[x];
+ m[x] = m[y];
+ m[y] = temp;
+
+ return;
+}
+
+
+void clearRow(matrix m, u_int32_t numcols, u_int32_t row)
+{
+ int32_t dwords = numcols/32;
+
+ if (numcols%32) dwords++;
+ memset(m[row],0,dwords*4);
+
+ return;
+}
+
+void displayRow(matrix m, u_int32_t row, u_int32_t numPrimes)
+{
+ int32_t length = numPrimes/32;
+ if (numPrimes%32) length++;
+ length*=64;
+
+ printf("[");
+ for (int32_t j = 0; j < length/2; j++)
+ {
+ if (getEntry(m,row,j)) printf("1");
+ else printf("0");
+ }
+ printf(" ");
+ for (int32_t j = length/2; j < length; j++)
+ {
+ if (getEntry(m,row,j)) printf("1");
+ else printf("0");
+ }
+ printf("]\n");
+ return;
+}
+
+void xorRows(matrix m, u_int32_t source, u_int32_t dest, u_int32_t length)
+{
+ u_int32_t i, q, r;
+ row x = m[dest];
+ row y = m[source];
+
+ r = length % 8; q = length - r;
+ for (i=0; i < q; i += 8)
+ {
+ x[i] ^= y[i]; x[1+i] ^= y[1+i]; x[2+i] ^= y[2+i]; x[3+i] ^= y[3+i];
+ x[4+i] ^= y[4+i]; x[5+i] ^= y[5+i]; x[6+i] ^= y[6+i]; x[7+i] ^= y[7+i];
+ }
+ switch (r)
+ {
+ case 7: x[i] ^= y[i]; i++;
+ case 6: x[i] ^= y[i]; i++;
+ case 5: x[i] ^= y[i]; i++;
+ case 4: x[i] ^= y[i]; i++;
+ case 3: x[i] ^= y[i]; i++;
+ case 2: x[i] ^= y[i]; i++;
+ case 1: x[i] ^= y[i]; i++;
+ }
+
+ return;
+}
+
+matrix constructMat(u_int32_t cols, u_int32_t rows)
+{
+ static matrix m;
+
+ u_int32_t dwords = cols/32;
+ if (cols%32) dwords++;
+ m = (row *) calloc(sizeof(row),rows);
+ if (m==NULL)
+ {
+ printf("Unable to allocate memory for matrix!\n");
+ exit(1);
+ }
+
+ for (u_int32_t i = 0; i < rows; i++)
+ {
+ m[i] = (row) calloc(2*dwords,sizeof(u_int32_t)); //two matrices, side by side
+ }
+ if (m[rows-1]==NULL)
+ {
+ printf("Unable to allocate memory for matrix!\n");
+ exit(1);
+ }
+
+ for (int32_t i = 0; i < rows; i++) //make second matrix identity, i.e. 1's along diagonal
+ {
+ insertEntry(m,i,i+32*dwords);
+ }
+
+ return m;
+}
+
+//===========================================================================
+// gaussReduce:
+//
+// Function: Apply Gaussian elimination to a matrix.
+//
+//===========================================================================
+u_int32_t gaussReduce(matrix m, u_int32_t numPrimes, u_int32_t relSought,int32_t extras)
+{
+ static u_int32_t rowUpto = 0;
+ static u_int32_t irow;
+ static u_int32_t length = (numPrimes+extras)/32;
+
+ if (numPrimes%32) length++;
+ length*=2;
+
+ for (int32_t icol = numPrimes-1; icol >= 0; icol--)
+ {
+ irow = rowUpto;
+ while ((irow < relSought)&&(getEntry(m,irow,icol)==0UL)) irow++;
+ if (irow < relSought)
+ {
+
+ swapRows(m,rowUpto,irow);
+
+ for (u_int32_t checkRow = rowUpto+1; checkRow<relSought; checkRow++)
+ {
+ if (getEntry(m,checkRow,icol)!=0UL)
+ {
+ xorRows(m,rowUpto,checkRow,length);
+ }
+ }
+
+ rowUpto++;
+ }
+ }
+
+ return rowUpto;
+}
+
diff --git a/F2matrix.h b/F2matrix.h
new file mode 100644
index 0000000..7142470
--- /dev/null
+++ b/F2matrix.h
@@ -0,0 +1,46 @@
+/*============================================================================
+ Copyright 2006 William Hart
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+============================================================================*/
+#ifndef F2MATRIX_H
+#define F2MATRIX_H
+
+#include <stdlib.h>
+#include "lanczos.h"
+
+typedef u_int32_t * row; //row of an F2 matrix
+typedef row * matrix; //matrix as a list of pointers to rows
+
+extern void insertEntry(matrix, u_int32_t, u_int32_t);
+extern void xorEntry(matrix, u_int32_t, u_int32_t);
+extern u_int32_t getEntry(matrix, u_int32_t, u_int32_t);
+extern matrix constructMat(u_int32_t, u_int32_t);
+extern void xorRows(row, row, u_int32_t);
+extern void clearRow(matrix, u_int32_t, u_int32_t);
+extern void swapRows(row, row);
+extern u_int32_t gaussReduce(matrix, u_int32_t, u_int32_t, int32_t);
+extern void displayRow(matrix, u_int32_t, u_int32_t);
+
+#endif
+
+
+
+
+
+
diff --git a/ModuloArith.cpp b/ModuloArith.cpp
new file mode 100644
index 0000000..95f17c5
--- /dev/null
+++ b/ModuloArith.cpp
@@ -0,0 +1,69 @@
+/*============================================================================
+ Copyright 2006 William Hart
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+============================================================================*/
+
+// -------------------------------------------------------
+//
+// ModuloArith.cpp
+//
+// Provides Functions for doing Modulo Arithmetic
+//
+// -------------------------------------------------------
+
+#include <gmp.h>
+#include "ModuloArith.h"
+
+mpz_t restemp; //chinese variables
+mpz_t ntemp;
+mpz_t chtemp;
+
+void modmul(mpz_t ab, mpz_t a, mpz_t b, mpz_t p)
+{
+ mpz_mul(ab,a,b);
+ mpz_fdiv_r(ab,ab,p);
+}
+
+void ChineseInit(void)
+{
+ mpz_init(restemp);
+ mpz_init(ntemp);
+ mpz_init(chtemp);
+
+ return;
+}
+
+void chinese(mpz_t res, mpz_t n, mpz_t x1, mpz_t x2, mpz_t n1, mpz_t n2)
+{
+ mpz_mul(ntemp,n1,n2);
+ mpz_invert(restemp,n2,n1);
+ modmul(restemp,restemp,n2,ntemp);
+ modmul(restemp,restemp,x1,ntemp);
+
+ mpz_invert(chtemp,n1,n2);
+ modmul(chtemp,chtemp,n1,ntemp);
+ modmul(chtemp,chtemp,x2,ntemp);
+
+ mpz_add(res,restemp,chtemp);
+ mpz_fdiv_r(res,res,ntemp);
+ mpz_set(n,ntemp);
+
+ return;
+}
+
diff --git a/ModuloArith.h b/ModuloArith.h
new file mode 100644
index 0000000..8e2e157
--- /dev/null
+++ b/ModuloArith.h
@@ -0,0 +1,42 @@
+/*============================================================================
+ Copyright 2006 William Hart
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+============================================================================*/
+
+// =====================================================================
+// INTERFACE:
+//
+// void ChineseInit(void)
+// - Initialise variables for chinese
+//
+// void modmul(mpz_t ab, mpz_t a, mpz_t b, mpz_t p)
+// - sets ab to a*b modulo p
+//
+// void chinese(mpz_t res, mpz_t n, mpz_t x1, mpz_t x2, mpz_t n1, mpz_t n2)
+// - sets n to n1*n2
+// - sets res mod n to a value congruent to x1 mod n1 and x2 mod n2
+//
+// ======================================================================
+
+extern void ChineseInit(void);
+extern void modmul(mpz_t, mpz_t, mpz_t, mpz_t);
+extern void chinese(mpz_t, mpz_t, mpz_t, mpz_t, mpz_t, mpz_t);
+
+
+
diff --git a/QS.cpp b/QS.cpp
new file mode 100644
index 0000000..4e02b8d
--- /dev/null
+++ b/QS.cpp
@@ -0,0 +1,1665 @@
+/*============================================================================
+ Copyright 2006 William Hart
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description:
+
+This is a relatively fast implementation of the self-initialising quadratic sieve.
+If you manage to improve the code, the author would like to hear about it.
+
+Contact: hart_wb {at-thingy} yahoo.com
+================================================================================*/
+
+#include <stdio.h>
+#include <stdarg.h>
+#include <stdlib.h>
+#include <math.h>
+#include <gmp.h>
+#include <string.h>
+#include <sys/times.h>
+#include <limits.h>
+
+#include "TonelliShanks.h"
+#include "ModuloArith.h"
+#include "F2matrix.h"
+#include "lanczos.h"
+#include "lprels.h"
+
+//===========================================================================
+//Uncomment these for various pieces of debugging information
+
+#define COUNT // Shows the number of relations generated and curves used during sieving
+//#define RELPRINT // Shows the actual factorizations of the relations
+//#define ERRORS // Error if relation should be divisible by a prime but isn't
+//#define POLS // Shows the polynomials being used by the sieve
+//#define ADETAILS // Prints some details about the factors of the A coefficients of the polys
+//#define LARGESTP // Prints the size of the largest factorbase prime
+//#define CURPARTS // Prints the number of curves used and number of partial relations
+//#define TIMING //displays some relative timings, if feature is available
+//#define REPORT //report sieve size, multiplier and number of primes used
+
+//===========================================================================
+//Architecture dependent fudge factors
+
+#if ULONG_MAX == 4294967295U
+#define SIEVEMASK 0xC0C0C0C0U
+#define MIDPRIME 1500
+#define SIEVEDIV 1
+#elif ULONG_MAX == 18446744073709551615U
+#define SIEVEMASK 0xC0C0C0C0C0C0C0C0U
+#define MIDPRIME 1500
+#define SIEVEDIV 1
+#endif
+
+#define CACHEBLOCKSIZE 64000 //Should be a little less than the L1/L2 cache size
+ //and a multiple of 64000
+#define MEDIUMPRIME 900
+#define SECONDPRIME 6000 //This should be lower for slower machines
+#define FUDGE 0.15 //Every program needs a mysterious fudge factor
+
+#define MINDIG 40 //Will not factor numbers with less than this number of decimal digits
+
+#define PREFETCH(addr,n) __builtin_prefetch((unsigned long*)addr+n,0,1)
+
+//===========================================================================
+//Knuth-Schroeppel multipliers and a macro to count them
+
+static const unsigned long multipliers[] = {1, 2, 3, 5, 7, 11, 13, 17, 19,
+ 23, 29, 31, 37, 41, 43};
+
+#define NUMMULTS (sizeof(multipliers)/sizeof(unsigned long))
+
+//===========================================================================
+// Large prime cutoffs
+
+unsigned long largeprimes[] =
+{
+ 250000, 300000, 370000, 440000, 510000, 580000, 650000, 720000, 790000, 8600000, //40-49
+ 930000, 1000000, 1700000, 2400000, 3100000, 3800000, 4500000, 5200000, 5900000, 6600000, //50-59
+ 7300000, 8000000, 8900000, 10000000, 11300000, 12800000, 14500000, 16300000, 18100000, 20000000, //60-69
+ 22000000, 24000000, 27000000, 32000000, 39000000, //70-74
+ 53000000, 65000000, 75000000, 87000000, 100000000, //75-79
+ 114000000, 130000000, 150000000, 172000000, 195000000, //80-84
+ 220000000, 250000000, 300000000, 350000000, 400000000, //85-89
+ 450000000, 500000000 //90-91
+};
+
+//============================================================================
+// Number of primes to use in factor base, given the number of decimal digits specified
+unsigned long primesNo[] =
+{
+ 1500, 1500, 1600, 1700, 1750, 1800, 1900, 2000, 2050, 2100, //40-49
+ 2150, 2200, 2250, 2300, 2400, 2500, 2600, 2700, 2800, 2900, //50-59
+ 3000, 3150, 5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000, //60-69
+ 9500, 10000, 11500, 13000, 15000, //70-74
+ 17000, 24000, 27000, 30000, 37000, //75-79
+ 45000, 47000, 53000, 57000, 58000, //80-84
+ 59000, 60000, 64000, 68000, 72000, //85-89
+ 76000, 80000 //90-91
+};
+
+//============================================================================
+// First prime actually sieved for
+unsigned long firstPrimes[] =
+{
+ 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, //40-49
+ 9, 8, 9, 9, 9, 9, 10, 10, 10, 10, //50-59
+ 10, 10, 11, 11, 12, 12, 13, 14, 15, 17, //60-69 //10
+ 19, 21, 22, 22, 23, //70-74
+ 24, 25, 25, 26, 26, //75-79
+ 27, 27, 27, 27, 28, //80-84
+ 28, 28, 28, 29, 29, //85-89
+ 29, 29 //90-91
+};
+
+//============================================================================
+// Logs of primes are rounded and errors accumulate; this specifies how great an error to allow
+unsigned long errorAmounts[] =
+{
+ 16, 17, 17, 18, 18, 19, 19, 19, 20, 20, //40-49
+ 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, //50-59
+ 24, 24, 25, 25, 25, 25, 26, 26, 26, 26, //60-69 //24
+ 27, 27, 28, 28, 29, //70-74
+ 29, 30, 30, 30, 31, //75-79
+ 31, 31, 31, 32, 32, //80-84
+ 32, 32, 32, 33, 33, //85-89
+ 33, 33 //90-91
+};
+
+//============================================================================
+// This is the threshold the sieve value must exceed in order to be considered for smoothness
+unsigned long thresholds[] =
+{
+ 66, 67, 67, 68, 68, 68, 69, 69, 69, 69, //40-49
+ 70, 70, 70, 71, 71, 71, 72, 72, 73, 73, //50-59
+ 74, 74, 75, 75, 76, 76, 77, 77, 78, 79, //60-69 //74
+ 80, 81, 82, 83, 84, //70-74
+ 85, 86, 87, 88, 89, //75-79
+ 91, 92, 93, 93, 94, //80-84
+ 95, 96, 97, 98, 100, //85-89
+ 101, 102 //90-91
+};
+
+//============================================================================
+// Size of sieve to use divided by 2, given the number of decimal digits specified
+//N.B: probably optimal if chosen to be a multiple of 32000, though other sizes are supported
+unsigned long sieveSize[] =
+{
+ 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, //40-49
+ 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, //50-59
+ 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, //60-69
+ 32000, 32000, 64000, 64000, 64000, //70-74
+ 96000, 96000, 96000, 128000, 128000, //75-79
+ 160000, 160000, 160000, 160000, 160000, //80-84
+ 192000, 192000, 192000, 192000, 192000, //85-89
+ 192000, 192000 //90-91
+};
+
+// Athlon tuning parameters
+/*unsigned long sieveSize[] =
+{
+ 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, //40-49
+ 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, 32000, //50-59
+ 64000, 64000, 64000, 64000, 64000, 64000, 64000, 64000, 64000, 64000, //60-69
+ 64000, 64000, 64000, 64000, 64000, //70-74
+ 128000, 128000, 128000, 128000, 128000, //75-79
+ 160000, 160000, 160000, 160000, 160000, //80-84
+ 192000, 192000, 192000, 192000, 192000, //85-89
+ 192000, 192000 //90-91
+};*/
+
+//============================================================================
+long decdigits; //number of decimal digits of n
+unsigned long secondprime; //min(numprimes, SECONDPRIME) = cutoff for using flags when sieving
+unsigned long firstprime; //first prime actually sieved with
+unsigned char errorbits; //first prime actually sieved with
+unsigned char threshold; //sieve threshold cutoff for smooth relations
+unsigned long midprime;
+unsigned long largeprime;
+
+unsigned long * factorBase; //array of factor base primes
+unsigned long numPrimes; //number of primes in factor base
+unsigned long relSought; //number of relations sought, i.e. a "few" more than numPrimes
+unsigned char * primeSizes; //array of sizes in bits, of the factor base primes
+unsigned char * sieve; //actual array where sieving takes place
+unsigned char * * offsets; //offsets for each prime to use in sieve
+unsigned char * * offsets2; //offsets for each prime to use in sieve (we switch between these)
+unsigned long relsFound =0; //number of relations found so far
+unsigned long potrels = 0; //potential relations (including duplicates)
+unsigned char * flags; //flags used for speeding up sieving for large primes
+unsigned long partials = 0; //number of partial relations
+unsigned long Mdiv2; //size of sieving interval divide 2
+unsigned long mat2off; //offset of second square block in matrix
+
+mpz_t * sqrts; //square roots of n modulo each prime in the factor base
+
+mpz_t n; //number to be factored
+mpz_t res; //smooth values which are trial factored
+
+mpz_t temp, temp2, temp3; //temporary variables
+mpz_t q,r; //quotient and remainder
+
+//Variables used for keeping time
+
+unsigned long clockstart;
+unsigned long clocktotal = 0;
+
+//Variables used by the modular inversion macro function
+long u1, u3;
+long v1, v3;
+long t1, t3, quot;
+
+//Variable used for random function
+unsigned long randval = 2994439072U;
+
+//==========================================================================================
+//Timing: provides some relative timings on X86 machines running gcc
+
+#if defined(__GNUC__) && (defined(__i386__) || defined(__x86_64__))
+
+#ifdef TIMING
+#define TIMES
+#endif
+
+double counterfirst[4];
+double countertotal[4] = {0,0,0,0};
+
+static unsigned counthi = 0;
+static unsigned countlo = 0;
+
+void counterasm(unsigned *hi, unsigned *lo)
+{
+ asm("rdtsc; movl %%edx,%0; movl %%eax,%1"
+ : "=r" (*hi), "=r" (*lo)
+ :
+ : "%edx", "%eax");
+}
+
+double getcounter()
+{
+ double total;
+
+ counterasm(&counthi, &countlo);
+
+ total = (double) counthi * (1 << 30) * 4 + countlo;
+ return total;
+}
+
+#endif
+
+/*========================================================================
+ Modular Inversion:
+
+ Function: GMP has a modular inverse function, but believe it or not,
+ this clumsy implementation is apparently quite a bit faster.
+ It inverts the value a, modulo the prime p, using the extended
+ gcd algorithm.
+
+========================================================================*/
+
+inline unsigned long modinverse(unsigned long a, unsigned long p)
+{
+ u1=1; u3=a;
+ v1=0; v3=p;
+ t1=0; t3=0;
+ while (v3)
+ {
+ quot=u3-v3;
+ if (u3 < (v3<<2))
+ {
+ if (quot < v3)
+ {
+ if (quot < 0)
+ {
+ t1 = u1; u1 = v1; v1 = t1;
+ t3 = u3; u3 = v3; v3 = t3;
+ } else
+ {
+ t1 = u1 - v1; u1 = v1; v1 = t1;
+ t3 = u3 - v3; u3 = v3; v3 = t3;
+ }
+ } else if (quot < (v3<<1))
+ {
+ t1 = u1 - (v1<<1); u1 = v1; v1 = t1;
+ t3 = u3 - (v3<<1); u3 = v3; v3 = t3;
+ } else
+ {
+ t1 = u1 - v1*3; u1 = v1; v1 = t1;
+ t3 = u3 - v3*3; u3 = v3; v3 = t3;
+ }
+ } else
+ {
+ quot=u3/v3;
+ t1 = u1 - v1*quot; u1 = v1; v1 = t1;
+ t3 = u3 - v3*quot; u3 = v3; v3 = t3;
+ }
+ }
+
+ if (u1<0) u1+=p;
+
+ return u1;
+}
+
+/*=========================================================================
+ Knuth_Schroeppel Multiplier:
+
+ Function: Find the best multiplier to use (allows 2 as a multiplier).
+ The general idea is to find a multiplier k such that kn will
+ be faster to factor. This is achieved by making kn a square
+ modulo lots of small primes. These primes will then be factor
+ base primes, and the more small factor base primes, the faster
+ relations will accumulate, since they hit the sieving interval
+ more often.
+
+==========================================================================*/
+unsigned long knuthSchroeppel(mpz_t n)
+{
+ float bestFactor = -10.0f;
+ unsigned long multiplier = 1;
+ unsigned long nmod8;
+ float factors[NUMMULTS];
+ float logpdivp;
+ mpz_t prime, r, mult;
+ long kron, multindex;
+
+ mpz_init(prime);
+ mpz_init(r);
+ mpz_init(mult);
+
+ nmod8 = mpz_fdiv_r_ui(r,n,8);
+
+ for (multindex = 0; multindex < NUMMULTS; multindex++)
+ {
+ long mod = nmod8*multipliers[multindex]%8;
+ factors[multindex] = 0.34657359; // ln2/2
+ if (mod == 1) factors[multindex] *= 4.0;
+ if (mod == 5) factors[multindex] *= 2.0;
+ factors[multindex] -= (log((float) multipliers[multindex]) / 2.0);
+ }
+
+ mpz_set_ui(prime,3);
+ while (mpz_cmp_ui(prime,10000)<0)
+ {
+ logpdivp = log((float)mpz_get_ui(prime)) / mpz_get_ui(prime);
+ kron = mpz_kronecker(n,prime);
+ for (multindex = 0; multindex < NUMMULTS; multindex++)
+ {
+ mpz_set_ui(mult,multipliers[multindex]);
+ switch (kron*mpz_kronecker(mult,prime))
+ {
+ case 0:
+ {
+ factors[multindex] += logpdivp;
+ } break;
+ case 1:
+ {
+ factors[multindex] += 2.0*logpdivp;
+ } break;
+ default: break;
+ }
+ }
+
+ mpz_nextprime(prime,prime);
+ }
+
+ for (multindex=0; multindex<NUMMULTS; multindex++)
+ {
+ if (factors[multindex] > bestFactor)
+ {
+ bestFactor = factors[multindex];
+ multiplier = multipliers[multindex];
+ }
+ }
+
+ mpz_clear(prime);
+ mpz_clear(r);
+ mpz_clear(mult);
+
+ return multiplier;
+}
+
+
+
+/*========================================================================
+ Initialize Quadratic Sieve:
+
+ Function: Initialises the global gmp variables.
+
+========================================================================*/
+void initSieve(void)
+{
+ mpz_init(n);
+ mpz_init(temp);
+ mpz_init(temp2);
+ mpz_init(temp3);
+ mpz_init(res);
+ mpz_init(q);
+ mpz_init(r);
+
+ return;
+}
+
+/*========================================================================
+ Compute Factor Base:
+
+ Function: Computes primes p up to B for which n is a square mod p,
+ allocates memory and stores them in an array pointed to by factorBase
+ Returns: number of primes actually in the factor base
+
+========================================================================*/
+void computeFactorBase(mpz_t n, unsigned long B,unsigned long multiplier)
+{
+ mpz_t currentPrime;
+ unsigned long primesinbase = 0;
+
+ factorBase = (unsigned long *) calloc(sizeof(unsigned long),B);
+
+ factorBase[primesinbase] = multiplier;
+ primesinbase++;
+ if (multiplier!=2)
+ {
+ factorBase[primesinbase] = 2;
+ primesinbase++;
+ }
+ mpz_init_set_ui(currentPrime,3);
+ while (primesinbase < B)
+ {
+ if (mpz_kronecker(n,currentPrime)==1)
+ {
+ factorBase[primesinbase] = mpz_get_ui(currentPrime);
+ primesinbase++;
+ }
+ mpz_nextprime(currentPrime,currentPrime);
+ }
+#ifdef LARGESTP
+ gmp_printf("Largest prime less than %Zd\n",currentPrime);
+#endif
+
+ mpz_clear(currentPrime);
+ return;
+}
+
+/*===========================================================================
+ Compute Prime Sizes:
+
+ Function: Computes the size in bits of each prime in the factor base
+ allocates memory for an array, primeSizes, to store the sizes
+ stores the size for each of the numPrimes primes in the array
+
+===========================================================================*/
+void computeSizes(unsigned long numPrimes)
+{
+ primeSizes = (unsigned char *) calloc(sizeof(unsigned char),numPrimes);
+ for (unsigned long i = 0; i<numPrimes; i++)
+ {
+ primeSizes[i]=(unsigned char)floor(log(factorBase[i])/log(2.0)-FUDGE+0.5);
+ }
+
+ return;
+}
+
+/*===========================================================================
+ Tonelli-Shanks:
+
+ Function: Performs Tonelli-Shanks on n mod every prime in the factor base
+ allocates memory for the results to be stored in the array sqrts
+
+===========================================================================*/
+void tonelliShanks(unsigned long numPrimes,mpz_t n)
+{
+ sqrts = (mpz_t *) calloc(sizeof(mpz_t),numPrimes);
+ mpz_array_init(sqrts[0],numPrimes,8*sizeof(unsigned long));
+
+ for (unsigned long i = 1; i<numPrimes; i++)
+ {
+ mpz_set_ui(temp,factorBase[i]);
+ sqrtmod(sqrts[i],n,temp);
+ }
+
+ return;
+}
+
+/*==========================================================================
+ evaluateSieve:
+
+ Function: searches sieve for relations and sticks them into a matrix, then
+ sticks their X and Y values into two arrays XArr and YArr
+
+===========================================================================*/
+void evaluateSieve(unsigned long ** relations, unsigned long ctimesreps, unsigned long M, unsigned char * sieve, mpz_t A, mpz_t B, mpz_t C, unsigned long * soln1, unsigned long * soln2, long polyadd, unsigned long * polycorr, mpz_t * XArr, unsigned long * aind, long min, long s,unsigned long multiplier, long * exponents, la_col_t* colarray, unsigned long * factors, char * rel_str, FILE* LPNEW,FILE* RELS)
+{
+ long i,j;
+ register unsigned long k;
+ unsigned long exponent, vv;
+ unsigned char extra;
+ register unsigned long modp;
+ unsigned long * sieve2;
+ unsigned char bits;
+ long numfactors;
+ unsigned long factnum;
+ char * last_ptr;
+ char Q_str[200];
+ char X_str[200];
+
+ i = 0;
+ j=0;
+ sieve2 = (unsigned long *) sieve;
+#ifdef POLS
+ gmp_printf("%Zdx^2%+Zdx\n%+Zd\n",A,B,C);
+#endif
+
+ while (j<M/sizeof(unsigned long))
+ {
+ do
+ {
+ while (!(sieve2[j] & SIEVEMASK)) j++;
+ i=j*sizeof(unsigned long);
+ j++;
+ while ((i<j*sizeof(unsigned long))&&(sieve[i] < threshold)) i++;
+ } while (sieve[i] < threshold);
+
+ if (i<M)
+ {
+ mpz_set_ui(temp,i+ctimesreps);
+ mpz_sub_ui(temp,temp,Mdiv2); //X
+
+ mpz_set(temp3,B); //B
+ mpz_addmul(temp3,A,temp); //AX+B
+ mpz_add(temp2,temp3,B); //AX+2B
+ mpz_mul(temp2,temp2,temp); //AX^2+2BX
+ mpz_add(res,temp2,C); //AX^2+2BX+C
+
+ bits=mpz_sizeinbase(res,2);
+ bits-=errorbits;
+
+ numfactors=0;
+
+ extra = 0;
+ if (factorBase[0]!=1)
+ {
+ mpz_set_ui(temp,factorBase[0]);
+ exponent = mpz_remove(res,res,temp);
+ exponents[0] = exponent;
+ if (exponent)
+ {
+ extra+=primeSizes[0];
+ }
+ }
+
+ mpz_set_ui(temp,factorBase[1]);
+ exponent = mpz_remove(res,res,temp);
+ exponents[1] = exponent;
+ extra+=exponent;
+
+ for (k = 2; k<firstprime; k++)
+ {
+ modp=(i+ctimesreps)%factorBase[k];
+
+ if (soln2[k]!=0xFFFFFFFFl)
+ {
+ if ((modp==soln1[k]) || (modp==soln2[k]))
+ {
+ mpz_set_ui(temp,factorBase[k]);
+ exponent = mpz_remove(res,res,temp);
+
+#ifdef ERRORS
+ if (exponent==0) printf("Error!\n");
+#endif
+ extra+=primeSizes[k];
+#ifdef RELPRINT
+ if (exponent > 0) printf(" %ld",(long)factorBase[k]);
+ if (exponent > 1) printf("^%ld",exponent);
+#endif
+ exponents[k] = exponent;
+ } else exponents[k] = 0;
+ } else
+ {
+ mpz_set_ui(temp,factorBase[k]);
+ exponent = mpz_remove(res,res,temp);
+ if (exponent) extra+=primeSizes[k];
+#ifdef RELPRINT
+ if (exponent > 0) gmp_printf(" %Zd",factorBase[k]);
+ if (exponent > 1) printf("^%ld",exponent);
+#endif
+ exponents[k] = exponent;
+ }
+ }
+ factnum = 0;
+ sieve[i]+=extra;
+ if (sieve[i] >= bits)
+ {
+ vv=((unsigned char)1<<(i&7));
+ for (k = firstprime; (k<secondprime)&&(extra<sieve[i]); k++)
+ {
+ modp=(i+ctimesreps)%factorBase[k];
+ if (soln2[k]!=0xFFFFFFFFl)
+ {
+ if ((modp==soln1[k]) || (modp==soln2[k]))
+ {
+ extra+=primeSizes[k];
+ mpz_set_ui(temp,factorBase[k]);
+ exponent = mpz_remove(res,res,temp);
+
+#ifdef ERRORS
+ if (exponent==0) printf("Error!\n");
+#endif
+
+#ifdef RELPRINT
+ if (exponent > 0) printf(" %ld",(long)factorBase[k]);
+ if (exponent > 1) printf("^%ld",exponent);
+#endif
+ factors[factnum+1] = k;
+ factors[factnum] = exponent;
+ factnum+=2;
+ }
+ } else
+ {
+ mpz_set_ui(temp,factorBase[k]);
+ exponent = mpz_remove(res,res,temp);
+ if (exponent) extra+=primeSizes[k];
+
+#ifdef RELPRINT
+ if (exponent > 0) printf(" %ld",(long)factorBase[k]);
+ if (exponent > 1) printf("^%ld",exponent);
+#endif
+ if (exponent)
+ {
+ factors[factnum+1] = k;
+ factors[factnum] = exponent;
+ factnum+=2;
+ }
+ }
+ }
+
+ for (k = secondprime; (k<numPrimes)&&(extra<sieve[i]); k++)
+ {
+ if (flags[k]&vv)
+ {
+ modp=(i+ctimesreps)%factorBase[k];
+ if ((modp==soln1[k]) || (modp==soln2[k]))
+ {
+ mpz_set_ui(temp,factorBase[k]);
+ exponent = mpz_remove(res,res,temp);
+#ifdef ERRORS
+ if (exponent==0) printf("Error!\n");
+#endif
+ extra+=primeSizes[k];
+#ifdef RELPRINT
+ if (exponent > 0) printf(" %ld",(long)factorBase[k]);
+ if (exponent > 1) printf("^%ld",exponent);
+#endif
+ factors[factnum+1] = k;
+ factors[factnum] = exponent;
+ factnum+=2;
+ }
+ }
+ }
+
+ last_ptr = rel_str;
+ if (mpz_cmp_ui(res,1000)>0)
+ {
+ if (mpz_cmp_ui(res,largeprime)<0)
+ {
+ for (unsigned long i = 0; i < firstprime; i++)
+ {
+ if (exponents[i]) add_factor(&last_ptr, (unsigned long) exponents[i], (unsigned long) i);
+ }
+ for (unsigned long i = 0; i < factnum; i+=2)
+ {
+ add_factor(&last_ptr, (unsigned long) factors[i], (unsigned long) factors[i+1]);
+ }
+ for (long i =0; i<s; i++)
+ {
+ add_factor(&last_ptr, (unsigned long) 1, (unsigned long) aind[i]+min);
+ }
+
+ add_0(&last_ptr);
+ gmp_sprintf(X_str, "%Zd\0", temp3);
+ gmp_sprintf(Q_str, "%Zd\0", res);
+ fprintf(LPNEW, "%s @ %s :%s\n", Q_str, X_str, rel_str);
+ partials++;
+ }
+#ifdef RELPRINT
+ gmp_printf(" %Zd\n",res);
+#endif
+ } else
+ {
+ mpz_neg(res,res);
+ if (mpz_cmp_ui(res,1000)>0)
+ {
+ if (mpz_cmp_ui(res,largeprime)<0)
+ {
+ for (unsigned long i = 0; i < firstprime; i++)
+ {
+ if (exponents[i]) add_factor(&last_ptr, (unsigned long) exponents[i], (unsigned long) i);
+ }
+ for (unsigned long i = 0; i < factnum; i+=2)
+ {
+ add_factor(&last_ptr, (unsigned long) factors[i], (unsigned long) factors[i+1]);
+ }
+ for (long i =0; i<s; i++)
+ {
+ add_factor(&last_ptr, (unsigned long) 1, (unsigned long) aind[i]+min);
+ }
+
+ add_0(&last_ptr);
+ gmp_sprintf(X_str, "%Zd\0", temp3);
+ gmp_sprintf(Q_str, "%Zd\0", res);
+ fprintf(LPNEW, "%s @ %s :%s\n", Q_str, X_str, rel_str);
+
+ partials++;
+ }
+#ifdef RELPRINT
+ gmp_printf(" %Zd\n",res);
+#endif
+ } else
+ {
+#ifdef RELPRINT
+ printf("....R\n");
+#endif
+ for (long i = 0; i<firstprime; i++)
+ {
+ if (exponents[i]) add_factor(&last_ptr, (unsigned long) exponents[i], (unsigned long) i);
+ }
+ for (unsigned long i = 0; i < factnum; i+=2)
+ {
+ add_factor(&last_ptr, (unsigned long) factors[i], (unsigned long) factors[i+1]);
+ }
+ for (long i =0; i<s; i++)
+ {
+ add_factor(&last_ptr, (unsigned long) 1, (unsigned long) aind[i]+min);
+ }
+
+ add_0(&last_ptr);
+ gmp_sprintf(X_str, "%Zd\0", temp3);
+ fprintf(RELS, "%s :%s\n", X_str, rel_str);
+
+ potrels++;
+ }
+ }
+ } else
+ {
+#ifdef RELPRINT
+ printf("\r \r");
+#endif
+
+ }
+ i++;
+
+ };
+ }
+
+ return;
+}
+
+
+/*=============================================================================
+ Sieve:
+
+ Function: Allocates space for a sieve of M integers and sieves the interval
+ starting at start
+
+=============================================================================*/
+void sieveInterval(unsigned long M, unsigned long numPrimes, unsigned char * sieve, long last, long first, long polyadd, unsigned long * soln1, unsigned long * soln2, unsigned long * polycorr, unsigned char * * offsets, unsigned char * * offsets2)
+{
+ register unsigned char currentprimesize;
+ register unsigned long currentprime;
+ unsigned char * position2;
+ register unsigned char * position;
+ register long diff;
+ unsigned char * end;
+ unsigned long ptimes4;
+ long correction;
+
+ end = sieve+M;
+
+ if (first)
+ {
+ for (unsigned long prime=1; prime<firstprime; prime++)
+ {
+ if (soln2[prime] == 0xFFFFFFFF) continue;
+ currentprime = factorBase[prime];
+ correction = polyadd ? -polycorr[prime]+currentprime : polycorr[prime];
+ soln1[prime]+=correction;
+ while (soln1[prime]>=currentprime) soln1[prime]-=currentprime;
+ soln2[prime]+=correction;
+ while (soln2[prime]>=currentprime) soln2[prime]-=currentprime;
+ }
+ }
+
+ for (unsigned long prime=firstprime; prime<MEDIUMPRIME; prime++)
+ {
+ if (soln2[prime] == 0xFFFFFFFF) continue;
+ currentprime = factorBase[prime];
+ currentprimesize = primeSizes[prime];
+
+ if (first)
+ {
+ correction = polyadd ? -polycorr[prime]+currentprime : polycorr[prime];
+ soln1[prime]+=correction;
+ while (soln1[prime]>=currentprime) soln1[prime]-=currentprime;
+ soln2[prime]+=correction;
+ while (soln2[prime]>=currentprime) soln2[prime]-=currentprime;
+
+ position = sieve+soln1[prime];
+ position2 = sieve+soln2[prime];
+ } else
+ {
+ position = offsets[prime];
+ position2 = offsets2[prime];
+ }
+ diff=position2-position;
+
+ ptimes4 = currentprime*4;
+ register unsigned char * bound=end-ptimes4;
+ while (bound - position > 0)
+ {
+ (* position)+=currentprimesize,(* (position+diff))+=currentprimesize, position+=currentprime;
+ (* position)+=currentprimesize,(* (position+diff))+=currentprimesize, position+=currentprime;
+ (* position)+=currentprimesize,(* (position+diff))+=currentprimesize, position+=currentprime;
+ (* position)+=currentprimesize,(* (position+diff))+=currentprimesize, position+=currentprime;
+ }
+ while ((end - position > 0)&&(end - position - diff > 0))
+ {
+ (* position)+=currentprimesize,(* (position+diff))+=currentprimesize, position+=currentprime;
+
+ }
+ position2 = position+diff;
+ if (end - position2 > 0)
+ {
+ (* position2)+=currentprimesize, position2+=currentprime;
+ }
+ if (end - position > 0)
+ {
+ (* position)+=currentprimesize, position+=currentprime;
+ }
+
+ if (!last)
+ {
+ offsets[prime] = position;
+ offsets2[prime] = position2;
+ }
+ }
+
+ for (unsigned long prime=MEDIUMPRIME; prime<midprime; prime++)
+ {
+ currentprime = factorBase[prime];
+ currentprimesize = primeSizes[prime];
+
+ if (first)
+ {
+ correction = polyadd ? -polycorr[prime]+currentprime : polycorr[prime];
+ soln1[prime]+=correction;
+ while (soln1[prime]>=currentprime) soln1[prime]-=currentprime;
+ soln2[prime]+=correction;
+ while (soln2[prime]>=currentprime) soln2[prime]-=currentprime;
+
+ position = sieve+soln1[prime];
+ position2 = sieve+soln2[prime];
+ } else
+ {
+ position = offsets[prime];
+ position2 = offsets2[prime];
+ }
+ diff=position2-position;
+
+ ptimes4 = 2*currentprime;
+ register unsigned char * bound=end-ptimes4;
+ while (bound - position > 0)
+ {
+ (* position)+=currentprimesize,(* (position+diff))+=currentprimesize, position+=currentprime;
+ (* position)+=currentprimesize,(* (position+diff))+=currentprimesize, position+=currentprime;
+ }
+ position2 = position+diff;
+ while ((end - position > 0)&&(end - position2 > 0))
+ {
+ (* position)+=currentprimesize, position+=currentprime, (* position2)+=currentprimesize, position2+=currentprime;
+
+ }
+
+ if (end - position2 > 0)
+ {
+ (* position2)+=currentprimesize, position2+=currentprime;
+ }
+ if (end - position > 0)
+ {
+ (* position)+=currentprimesize, position+=currentprime;
+ }
+ if (!last)
+ {
+ offsets[prime] = position;
+ offsets2[prime] = position2;
+ }
+ }
+
+ return;
+}
+
+/*===========================================================================
+ Sieve 2:
+
+ Function: Second sieve for larger primes
+
+=========================================================================== */
+void sieve2(unsigned long M, unsigned long numPrimes, unsigned char * sieve, long last, long first, long polyadd, unsigned long * soln1, unsigned long * soln2, unsigned long * polycorr, unsigned char * * offsets, unsigned char * * offsets2)
+{
+ register unsigned char currentprimesize;
+ register unsigned long currentprime;
+ register unsigned char * position2;
+ register unsigned char * position;
+ unsigned char * end;
+ long correction;
+
+ memset(sieve,0,M*sizeof(unsigned char));
+ memset(flags,0,numPrimes*sizeof(unsigned char));
+ end = sieve+M;
+ *end = 255; //sentinel to speed up sieve evaluators inner loop
+
+ for (unsigned long prime=midprime; prime<secondprime; prime++)
+ {
+ currentprime = factorBase[prime];
+ currentprimesize = primeSizes[prime];
+
+ correction = polyadd ? -polycorr[prime]+currentprime : polycorr[prime];
+ soln1[prime]+=correction;
+ while (soln1[prime]>=currentprime) soln1[prime]-=currentprime;
+ soln2[prime]+=correction;
+ while (soln2[prime]>=currentprime) soln2[prime]-=currentprime;
+
+ position = sieve+soln1[prime];
+ position2 = sieve+soln2[prime];
+
+ while ((end - position > 0)&&(end - position2 > 0))
+ {
+ (* position)+=currentprimesize, position+=currentprime, (* position2)+=currentprimesize, position2+=currentprime;
+ }
+
+ if (end - position2 > 0)
+ {
+ (* position2)+=currentprimesize;
+ }
+ if (end - position > 0)
+ {
+ (* position)+=currentprimesize;
+ }
+ }
+
+ for (unsigned long prime=secondprime; prime<numPrimes; prime++)
+ {
+ currentprime = factorBase[prime];
+ currentprimesize = primeSizes[prime];
+
+ correction = polyadd ? -polycorr[prime]+currentprime : polycorr[prime];
+ soln1[prime]+=correction;
+ while (soln1[prime]>=currentprime) soln1[prime]-=currentprime;
+ soln2[prime]+=correction;
+ while (soln2[prime]>=currentprime) soln2[prime]-=currentprime;
+
+ position = sieve+soln1[prime];
+ position2 = sieve+soln2[prime];
+
+ while (end - position > 0)
+ {
+ flags[prime]|=((unsigned char)1<<((position-sieve)&7)), (* position)+=currentprimesize, position+=currentprime;
+ }
+
+ while (end - position2 > 0)
+ {
+ flags[prime]|=((unsigned char)1<<((position2-sieve)&7)), (* position2)+=currentprimesize, position2+=currentprime;
+ }
+ }
+
+ return;
+}
+
+/*============================================================================
+
+ random:
+
+ Function: Generates a pseudo-random integer between 0 and n-1 inclusive
+
+============================================================================*/
+unsigned long random(unsigned long upto)
+{
+ randval = ((u_int64_t)randval*1025416097U+286824428U)%(u_int64_t)4294967291U;
+ return randval%upto;
+}
+
+
+/*============================================================================
+ mainRoutine:
+
+ Function: Generates the polynomials, initialises and calls the sieve,
+ implementing cache blocking (breaking the sieve interval into
+ small blocks for the small primes.
+
+============================================================================*/
+void mainRoutine(unsigned long Mdiv2, mpz_t n, unsigned long multiplier)
+{
+ mpz_t A; mpz_init(A);
+ mpz_t B; mpz_init(B);
+ mpz_t C; mpz_init(C);
+ mpz_t D; mpz_init(D);
+ mpz_t temp; mpz_init(temp);
+ mpz_t temp2; mpz_init(temp2);
+ mpz_t q; mpz_init(q);
+ mpz_t r; mpz_init(r);
+ mpz_t Bdivp2; mpz_init(Bdivp2);
+ mpz_t factor; mpz_init(factor);
+
+ unsigned long u1;
+
+ long s, fact, span, min;
+ unsigned long p;
+ unsigned long reps;
+
+ unsigned long curves = 0;
+
+ unsigned long ** relations;
+ long * primecount;
+
+ long * exponents = (long *) calloc(firstprime,sizeof(long));
+ if (exponents==NULL)
+ {
+ printf("Unable to allocate memory!\n");
+ abort();
+ }
+ unsigned long factors[200];
+ char rel_str[MPQS_STRING_LENGTH];
+
+ unsigned long totcomb = 0;
+
+ unsigned long next_cutoff = (relSought - 1)/40 +1;
+ unsigned long next_inc = next_cutoff;
+
+ FILE * LPNEW;
+ FILE * LPRELS;
+ FILE * COMB;
+ FILE * FNEW;
+ FILE * RELS;
+ FILE * FRELS;
+ FILE * FLPRELS;
+ LPNEW = flint_fopen("lpnew","w");
+ LPRELS = flint_fopen("lprels","w");
+ RELS = flint_fopen("rels","w");
+ FNEW = flint_fopen("fnew","w");
+ fclose(FNEW);
+ FLPRELS = flint_fopen("flprels","w");
+ fclose(FLPRELS);
+ FRELS = flint_fopen("frels","w");
+ fclose(LPRELS);
+ fclose(FRELS);
+
+
+#ifdef TIMES
+ counterfirst[2] = getcounter();
+#endif
+ s = mpz_sizeinbase(n,2)/28+1;
+
+ unsigned long * aind = (unsigned long*) calloc(sizeof(unsigned long),s);
+ unsigned long * amodp = (unsigned long*) calloc(sizeof(unsigned long),s);
+ unsigned long * Ainv = (unsigned long*) calloc(sizeof(unsigned long),numPrimes);
+ unsigned long * soln1 = (unsigned long*) calloc(sizeof(unsigned long),numPrimes);
+ unsigned long * soln2 = (unsigned long*) calloc(sizeof(unsigned long),numPrimes);
+ unsigned long ** Ainv2B = (unsigned long**) calloc(sizeof(unsigned long*),s);
+ if (Ainv2B==NULL)
+ {
+ printf("Unable to allocate memory!\n");
+ abort();
+ }
+ for (long i=0; i<s; i++)
+ {
+ Ainv2B[i] = (unsigned long *) calloc(sizeof(unsigned long),numPrimes);
+ if (Ainv2B[i]==NULL)
+ {
+ printf("Unable to allocate memory!\n");
+ abort();
+ }
+ }
+
+ mpz_t * Bterms = (mpz_t *)calloc(sizeof(mpz_t),s);
+ mpz_array_init(*Bterms,s,mpz_sizeinbase(n,2));
+
+ la_col_t* colarray = (la_col_t*)calloc(relSought,sizeof(la_col_t)); //initialise a zeroed array of column structures
+
+ relsFound = 0;
+
+ mpz_t XArr[relSought];
+ for(unsigned long i = 0; i < relSought; i++) mpz_init(XArr[i]);
+
+ sieve = (unsigned char *) calloc(sizeof(unsigned char),Mdiv2*2+4); //one dword extra for sentinel to speed up sieve evaluation loop
+ if (sieve==NULL)
+ {
+ printf("Unable to allocate memory for sieve!\n");
+ abort();
+ }
+
+
+ flags = (unsigned char*) calloc(sizeof(unsigned char),numPrimes);
+
+ offsets = (unsigned char * *)calloc(numPrimes,sizeof(unsigned char *));
+ offsets2 = (unsigned char * *)calloc(numPrimes,sizeof(unsigned char *));
+
+ relations = (unsigned long * *) calloc(relSought,sizeof(unsigned long *));
+ for (long i = 0; i < relSought; i++)
+ {
+ relations[i] = (unsigned long *) calloc(200, sizeof(unsigned long));
+ }
+
+ primecount = (long *) calloc(numPrimes, sizeof(long));
+
+//Compute min A_prime and A_span
+
+ mpz_mul_ui(temp,n,2);
+ mpz_sqrt(temp,temp);
+ mpz_div_ui(temp2,temp,Mdiv2);
+ mpz_root(temp,temp2,s);
+ for (fact = 0; mpz_cmp_ui(temp,factorBase[fact])>=0; fact++);
+ span = numPrimes/s/s/2;
+ min=fact-span/2;
+ while ((fact*fact)/min - min < span) {min--;}
+
+#ifdef ADETAILS
+ printf("s = %ld, fact = %ld, min = %ld, span = %ld\n",s,fact,min,span);
+#endif
+
+//Compute first polynomial and adjustments
+
+ while (relsFound + totcomb < relSought)
+ {
+ long i,j;
+ long ran;
+
+ mpz_set_ui(A,1);
+
+ for (i = 0; i < s-1; )
+ {
+ j=-1L;
+ ran = span/2+random(span/2);
+ while (j!=i)
+ {
+ ran++;
+ for (j=0;((j<i)&&(aind[j]!=ran));j++);
+ }
+ aind[i] = ran;
+ mpz_mul_ui(A,A,factorBase[ran+min]);
+ i++;
+ if (i < s-1)
+ {
+ j=-1L;
+ ran = ((min+span/2)*(min+span/2))/(ran+min) - random(10)-min;
+ while (j!=i)
+ {
+ ran++;
+ for (j=0;((j<i)&&(aind[j]!=ran));j++);
+ }
+ aind[i] = ran;
+ mpz_mul_ui(A,A,factorBase[ran+min]);
+ i++;
+ }
+ }
+
+ mpz_div(temp,temp2,A);
+ for (fact = 1; mpz_cmp_ui(temp,factorBase[fact])>=0; fact++);
+ fact-=min;
+ do
+ {
+ for (j=0;((j<i)&&(aind[j]!=fact));j++);
+ fact++;
+ } while (j!=i);
+ fact--;
+ aind[i] = fact;
+ mpz_mul_ui(A,A,factorBase[fact+min]);
+
+ for (long i=0; i<s; i++)
+ {
+ p = factorBase[aind[i]+min];
+ mpz_div_ui(temp,A,p);
+ amodp[i] = mpz_fdiv_r_ui(temp,temp,p);
+
+ mpz_set_ui(temp,modinverse(mpz_get_ui(temp),p));
+ mpz_mul(temp,temp,sqrts[aind[i]+min]);
+ mpz_fdiv_r_ui(temp,temp,p);
+ if (mpz_cmp_ui(temp,p/2)>0)
+ {
+ mpz_sub_ui(temp,temp,p);
+ mpz_neg(temp,temp);
+ }
+ mpz_mul(temp,temp,A);
+ mpz_div_ui(Bterms[i],temp,p);
+ }
+
+ mpz_set(B,Bterms[0]);
+ for (long i=1; i<s; i++)
+ {
+ mpz_add(B,B,Bterms[i]);
+ }
+
+ for (long i=0; i<numPrimes; i++)
+ {
+ p = factorBase[i];
+ Ainv[i] = modinverse(mpz_fdiv_r_ui(temp,A,p),p);
+
+ for (long j=0; j<s; j++)
+ {
+ mpz_fdiv_r_ui(temp,Bterms[j],p);
+ mpz_mul_ui(temp,temp,2*Ainv[i]);
+ Ainv2B[j][i] = mpz_fdiv_r_ui(temp,temp,p);
+ }
+
+ mpz_fdiv_r_ui(temp,B,p);
+ mpz_sub(temp,sqrts[i],temp);
+ mpz_add_ui(temp,temp,p);
+ mpz_mul_ui(temp,temp,Ainv[i]);
+ mpz_add_ui(temp,temp,Mdiv2);
+ soln1[i] = mpz_fdiv_r_ui(temp,temp,p);
+ mpz_sub_ui(temp,sqrts[i],p);
+ mpz_neg(temp,temp);
+ mpz_mul_ui(temp,temp,2*Ainv[i]);
+ soln2[i] = mpz_fdiv_r_ui(temp,temp,p)+soln1[i];
+ }
+
+ for (long polyindex=1; polyindex<(1<<(s-1))-1; polyindex++)
+ {
+ long j,polyadd;
+ unsigned long * polycorr;
+ for (j=0; j<s; j++)
+ {
+ if (((polyindex>>j)&1)!=0) break;
+ }
+ if ((polyadd = (((polyindex>>j)&2)!=0)))
+ {
+ mpz_add(B,B,Bterms[j]);
+ mpz_add(B,B,Bterms[j]);
+ } else
+ {
+ mpz_sub(B,B,Bterms[j]);
+ mpz_sub(B,B,Bterms[j]);
+ }
+ polycorr = Ainv2B[j];
+
+ long index;
+ for (long j=0; j<s; j++)
+ {
+ index = aind[j]+min;
+ p = factorBase[index];
+ mpz_fdiv_r_ui(D,n,p*p);
+ mpz_fdiv_r_ui(Bdivp2,B,p*p);
+ mpz_mul_ui(temp,Bdivp2,amodp[j]);
+ mpz_realloc2(temp3,64);
+ mpz_fdiv_r_ui(temp,temp,p);
+ u1 = modinverse(mpz_fdiv_r_ui(temp,temp,p),p);
+ mpz_mul(temp,Bdivp2,Bdivp2);
+ mpz_sub(temp,temp,D);
+ mpz_neg(temp,temp);
+ mpz_div_ui(temp,temp,p);
+ mpz_mul_ui(temp,temp,u1);
+ mpz_add_ui(temp,temp,Mdiv2);
+ mpz_add_ui(temp,temp,p);
+ soln1[index]=mpz_fdiv_r_ui(temp,temp,p);
+ soln2[index]=0xFFFFFFFFl;
+ }
+
+// Count the number of polynomial curves used so far and compute the C coefficient of our polynomial
+
+ curves++;
+
+ mpz_mul(C,B,B);
+ mpz_sub(C,C,n);
+ mpz_divexact(C,C,A);
+
+// Do the sieving and relation collection
+
+ mpz_set_ui(temp,Mdiv2*2);
+ mpz_fdiv_qr_ui(q,r,temp,CACHEBLOCKSIZE);
+#ifdef TIMES
+ counterfirst[3] = getcounter();
+#endif
+ sieve2(mpz_get_ui(temp),numPrimes,sieve,1,1,polyadd,soln1,soln2,polycorr,offsets,offsets2);
+#ifdef TIMES
+ countertotal[3]+=(getcounter()-counterfirst[3]);
+ counterfirst[0] = getcounter();
+#endif
+ sieveInterval(CACHEBLOCKSIZE,numPrimes,sieve,0,1,polyadd,soln1,soln2,polycorr,offsets,offsets2);
+ if (mpz_cmp_ui(q,1)>0)
+ {
+ for (reps = 1;reps < mpz_get_ui(q)-1; reps++)
+ {
+ sieveInterval(CACHEBLOCKSIZE,numPrimes,sieve+CACHEBLOCKSIZE*reps,0,0,polyadd,soln1,soln2,polycorr,offsets,offsets2);
+ }
+ if (mpz_cmp_ui(r,0)==0)
+ {
+ sieveInterval(CACHEBLOCKSIZE,numPrimes,sieve+CACHEBLOCKSIZE*reps,1,0,polyadd,soln1,soln2,polycorr,offsets,offsets2);
+ } else
+ {
+ sieveInterval(CACHEBLOCKSIZE,numPrimes,sieve+CACHEBLOCKSIZE*reps,0,0,polyadd,soln1,soln2,polycorr,offsets,offsets2);
+ reps++;
+ sieveInterval(mpz_get_ui(r),numPrimes,sieve+CACHEBLOCKSIZE*reps,1,0,polyadd,soln1,soln2,polycorr,offsets,offsets2);
+ }
+ }
+
+#ifdef TIMES
+ countertotal[0]+=(getcounter()-counterfirst[0]);
+ counterfirst[1] = getcounter();
+#endif
+ evaluateSieve(relations,0,mpz_get_ui(temp),sieve,A,B,C,soln1, soln2, polyadd, polycorr,XArr,aind,min,s,multiplier,exponents,colarray,factors,rel_str,LPNEW,RELS);
+
+ if (2*potrels >= next_cutoff)
+ {
+ fclose(LPNEW);
+ sort_lp_file("lpnew");
+ COMB = flint_fopen("comb","w");
+ mergesort_lp_file("lprels", "lpnew", "tmp", COMB);
+ fclose(COMB);
+ LPNEW = flint_fopen("lpnew","w");
+
+ fclose(RELS);
+ sort_lp_file("rels");
+ relsFound = mergesort_lp_file("frels","rels","tmp2",NULL);
+ RELS = flint_fopen("rels","w");
+
+ COMB = flint_fopen("comb", "r");
+ FNEW = flint_fopen("fnew","w");
+ combine_large_primes(numPrimes, COMB, FNEW, n, factor);
+ fclose(FNEW);
+ fclose(COMB);
+ sort_lp_file("fnew");
+ totcomb = mergesort_lp_file("flprels","fnew","tmp3",NULL);
+#ifdef COUNT
+ printf("%ld full relations, %ld combined relations\n",relsFound,totcomb);
+#endif
+ if ((next_cutoff < relSought) && (next_cutoff + next_inc/2 >= relSought))
+ next_inc = next_inc/2;
+ next_cutoff += next_inc;
+ }
+#ifdef TIMES
+ countertotal[1]+=(getcounter()-counterfirst[1]);
+#endif
+ }
+
+#ifdef COUNT
+ if (curves%20==0) printf("%ld curves.\n",(long)curves);
+#endif
+ }
+
+#ifdef CURPARTS
+ printf("%ld curves, %ld partials.\n",(long)curves,(long)partials);
+#endif
+
+#ifdef REPORT
+ printf("Done with sieving!\n");
+#endif
+
+ unsigned long ncols = relSought;
+ unsigned long nrows = numPrimes;
+
+#ifdef ERRORS
+ for (unsigned long j = relsFound; j<relSought; j++)
+ {
+ if (colarray[j].weight != 0) printf("Dirty at col %ld\n",j);
+ }
+#endif
+
+ fclose(LPNEW);
+ fclose(RELS);
+
+#ifdef COUNT
+ printf("%ld relations found in total!\n",totcomb+relsFound);
+#endif
+
+ FRELS = flint_fopen("frels","r");
+ relsFound = 0;
+ read_matrix(relations, FRELS, colarray, &relsFound, relSought, XArr, n, factorBase);
+ fclose(FRELS);
+ FLPRELS = flint_fopen("flprels","r");
+ read_matrix(relations, FLPRELS, colarray, &relsFound, relSought, XArr, n, factorBase);
+ fclose(FLPRELS);
+
+#ifdef ERRORS
+ for (unsigned long j = 0; j<relSought; j++)
+ {
+ if (colarray[j].orig != j)
+ {
+ printf("Column numbering error, %ld\n",j);
+ colarray[j].orig = j;
+ }
+ }
+
+ for (unsigned long j = 0; j<relSought; j++)
+ for (unsigned long i = 0; i<colarray[j].weight; i++)
+ if (colarray[j].data[i] > numPrimes) printf("Error prime too large: %ld\n",colarray[j].data[i]);
+
+ mpz_t test1;
+ mpz_init(test1);
+ mpz_t test2;
+ mpz_init(test2);
+ mpz_t test3;
+ mpz_init(test3);
+ unsigned long * exps = (unsigned long *) malloc(numPrimes*sizeof(unsigned long));
+ for (unsigned long j = 0; j<relSought; j++)
+ {
+ for (unsigned long i = 0; i < numPrimes; i++) exps[i] = 0;
+ mpz_set_ui(test1,1);
+ for (unsigned long i = 1; i<=relations[j][0]; i++)
+ {
+ mpz_mul_ui(test1,test1,factorBase[relations[j][i]]);
+ exps[relations[j][i]]++;
+ }
+ mpz_mod(test1,test1,n);
+ mpz_mul(test2,XArr[j],XArr[j]);
+ mpz_mod(test2,test2,n);
+ if (mpz_cmp(test1,test2)!=0)
+ {
+ mpz_add(test3,test1,test2);
+ if (mpz_cmp(test3,n)!=0)
+ {
+ gmp_printf("%Zd !=\n%Zd\nin column %ld and\n",test3,n,j);
+ gmp_printf("%Zd !=\n%Zd\n\n",test1,test2);
+ }
+ }
+ for (unsigned long i = 0; i < colarray[j].weight; i++)
+ {
+ if (exps[colarray[j].data[i]]%2 != 1) printf("Col %ld, row %ld incorrect\n",j,i);
+ exps[colarray[j].data[i]] = 0;
+ }
+ for (unsigned long i = 0; i < numPrimes; i++) if (exps[i]%2==1) printf("exps[%ld] is not even in row %ld\n",i,j);
+ }
+ free(exps);
+ mpz_clear(test1);
+ mpz_clear(test2);
+ mpz_clear(test3);
+#endif
+
+ reduce_matrix(&nrows, &ncols, colarray);
+
+#ifdef ERRORS
+ exps = (unsigned long *) malloc(numPrimes*sizeof(unsigned long));
+ for (unsigned long j = 0; j<ncols; j++)
+ {
+ for (unsigned long i = 0; i < numPrimes; i++) exps[i] = 0;
+ for (unsigned long i = 1; i<=relations[colarray[j].orig][0]; i++)
+ {
+ exps[relations[colarray[j].orig][i]]++;
+ }
+ for (unsigned long i = 0; i < colarray[j].weight; i++)
+ {
+ for (unsigned long k = 0; k < i; k++)
+ if (colarray[j].data[i] == colarray[j].data[k]) printf("Duplicate in column %ld: %ld, %ld\n",j,i,k);
+ if (exps[colarray[j].data[i]]%2 != 1) printf("Col %ld, row %ld incorrect\n",j,i);
+ exps[colarray[j].data[i]] = 0;
+ }
+ for (unsigned long i = 0; i < numPrimes; i++) if (exps[i]%2==1) printf("exps[%ld] is not even in row %ld\n",i,j);
+
+ }
+#endif
+
+ u_int64_t* nullrows;
+ do {
+ nullrows = block_lanczos(nrows, 0, ncols, colarray);
+ } while (nullrows == NULL);
+
+ long i,j, mask;
+
+ for (i = 0, mask = 0; i < ncols; i++)
+ mask |= nullrows[i];
+
+ for (i = j = 0; i < 64; i++) {
+ if (mask & ((u_int64_t)(1) << i))
+ j++;
+ }
+#ifdef REPORT
+ printf ("%ld nullspace vectors found.\n",j);
+#endif
+
+
+#ifdef ERRORS
+ exps = (unsigned long *) malloc(numPrimes*sizeof(unsigned long));
+ for (unsigned long j = 0; j<ncols; j++)
+ {
+ for (unsigned long i = 0; i < numPrimes; i++) exps[i] = 0;
+ for (unsigned long i = 1; i<=relations[colarray[j].orig][0]; i++)
+ {
+ exps[relations[colarray[j].orig][i]]++;
+ }
+ for (unsigned long i = 0; i < colarray[j].weight; i++)
+ {
+ if ((exps[colarray[j].data[i]]%2) != 1) printf("Col %ld, row %ld incorrect\n",j,i);
+ exps[colarray[j].data[i]] = 0;
+ }
+ for (unsigned long i = 0; i < numPrimes; i++) if ((exps[i]%2)==1) printf("exps[%ld] is not even in row %ld\n",i,j);
+ }
+#endif
+
+
+#ifdef TIMES
+ countertotal[2]+=(getcounter()-counterfirst[2]);
+ printf("Total time = %.0f\n",countertotal[2]);
+ printf("Polynomial generation time = %.0f\n",countertotal[2]-countertotal[0]-countertotal[1]-countertotal[3]);
+ printf("Small prime sieve time = %.0f\n",countertotal[0]);
+ printf("Large prime sieve time = %.0f\n",countertotal[3]);
+ printf("Evaluate sieve time = %.0f\n",countertotal[1]);
+#endif
+
+// We want factors of n, not kn, so divide out by the multiplier
+
+ mpz_div_ui(n,n,multiplier);
+
+// Now do the "square root" and GCD steps hopefully obtaining factors of n
+ printf("FACTORS:\n");
+ for (long l = 0;l<64;l++)
+ {
+ while (!(mask & ((u_int64_t)(1) << l))) l++;
+ mpz_set_ui(temp,1);
+ mpz_set_ui(temp2,1);
+ memset(primecount,0,numPrimes*sizeof(unsigned long));
+ for (long i = 0; i<ncols; i++)
+ {
+ if (getNullEntry(nullrows,i,l))
+ {
+ mpz_mul(temp2,temp2,XArr[colarray[i].orig]);
+ for (long j=1; j<=relations[colarray[i].orig][0]; j++)
+ {
+ primecount[relations[colarray[i].orig][j]]++;
+ }
+ }
+ if (i%30==0) mpz_mod(temp2,temp2,n);
+ }
+ for (long j = 0; j<numPrimes; j++)
+ {
+ mpz_set_ui(temp3,factorBase[j]);
+ mpz_pow_ui(temp3,temp3,primecount[j]/2);
+ mpz_mul(temp,temp,temp3);
+ if (j%30==0) mpz_mod(temp,temp,n);
+ }
+ mpz_sub(temp,temp2,temp);
+ mpz_gcd(temp,temp,n);
+ if ((mpz_cmp(temp,n)!=0)&&(mpz_cmp_ui(temp,1)!=0)) //print only nontrivial factors
+ {
+ gmp_printf("%Zd\n",temp);
+ }
+ }
+ mpz_clear(factor);
+ return;
+}
+
+/*===========================================================================
+ Main Program:
+
+ Function: Factors a user specified number using a quadratic sieve
+
+===========================================================================*/
+int main(int argc, unsigned char *argv[])
+{
+ unsigned long multiplier;
+
+ initSieve();
+
+ printf("Input number to factor [ >=40 decimal digits]: ");
+ gmp_scanf("%Zd",n);getchar();
+
+ decdigits = mpz_sizeinbase(n,10);
+ if (decdigits < 40)
+ {
+ printf("Error in input or number has too few digits.\n");
+ abort();
+ }
+
+ multiplier = knuthSchroeppel(n);
+ mpz_mul_ui(n,n,multiplier);
+
+ if (decdigits<=91)
+ {
+ numPrimes=primesNo[decdigits-MINDIG];
+
+ Mdiv2 = sieveSize[decdigits-MINDIG]/SIEVEDIV;
+ if (Mdiv2*2 < CACHEBLOCKSIZE) Mdiv2 = CACHEBLOCKSIZE/2;
+ largeprime = largeprimes[decdigits-MINDIG];
+
+#ifdef REPORT
+ printf("Using multiplier: %ld\n",(long)multiplier);
+ printf("%ld primes in factor base.\n",(long)numPrimes);
+ printf("Sieving interval M = %ld\n",(long)Mdiv2*2);
+ printf("Large prime cutoff = factorBase[%ld]\n",largeprime);
+#endif
+
+ if (numPrimes < SECONDPRIME) secondprime = numPrimes;
+ else secondprime = SECONDPRIME;
+ if (numPrimes < MIDPRIME) midprime = numPrimes;
+ else midprime = MIDPRIME;
+
+ firstprime = firstPrimes[decdigits-MINDIG];
+ errorbits = errorAmounts[decdigits-MINDIG];
+ threshold = thresholds[decdigits-MINDIG];
+
+ } else //all bets are off
+ {
+ numPrimes = 64000;
+ Mdiv2 = 192000/SIEVEDIV;
+ largeprime = numPrimes*10*decdigits;
+
+#ifdef REPORT
+ printf("Using multiplier: %ld\n",(long)multiplier);
+ printf("%ld primes in factor base.\n",(long)numPrimes);
+ printf("Sieving interval M = %ld\n",(long)Mdiv2*2);
+ printf("Large prime cutoff = factorBase[%ld]\n",largeprime);
+#endif
+
+ secondprime = SECONDPRIME;
+ midprime = MIDPRIME;
+ firstprime = 30;
+ errorbits = decdigits/4 + 2;
+ threshold = 43+(7*decdigits)/10;
+
+ }
+ relSought = numPrimes+64;
+ computeFactorBase(n, numPrimes, multiplier);
+
+ computeSizes(numPrimes);
+
+ TonelliInit();
+ tonelliShanks(numPrimes,n);
+
+ mainRoutine(Mdiv2, n,multiplier);
+
+ getchar();
+#if defined(WINCE) || defined(macintosh)
+ char * tmp_dir = NULL;
+#else
+ char * tmp_dir = getenv("TMPDIR");
+#endif
+ if (tmp_dir == NULL) tmp_dir = "./";
+ char * delfile;
+
+ delfile = get_filename(tmp_dir,unique_filename("comb"));
+ remove(delfile);
+ delfile = get_filename(tmp_dir,unique_filename("frels"));
+ remove(delfile);
+ delfile = get_filename(tmp_dir,unique_filename("flprels"));
+ remove(delfile);
+ delfile = get_filename(tmp_dir,unique_filename("lpnew"));
+ remove(delfile);
+ delfile = get_filename(tmp_dir,unique_filename("rels"));
+ remove(delfile);
+ delfile = get_filename(tmp_dir,unique_filename("fnew"));
+ remove(delfile);
+ delfile = get_filename(tmp_dir,unique_filename("lprels"));
+ remove(delfile);
+
+ return 0;
+}
+
diff --git a/QStodo.txt b/QStodo.txt
new file mode 100644
index 0000000..06929bb
--- /dev/null
+++ b/QStodo.txt
@@ -0,0 +1,50 @@
+QStodo
+=======
+
+1. ***Check for duplicate ordinary relations
+2. ***Make large prime cutoff table
+3. ***Change sievediv to 1 for all but 60 digit factorizations or separate into tuning files
+4. Make sieve look for additional relations if factorization fails
+5. ***Set all LESS values to 0
+6. ***Update old version of sieve to include new corrected polynomial choosing
+7. ***Make better polynomial choosing code
+8. ***Write ordinary relations to file
+9. ***Fix all printf %d's etc
+10. Make structures in memory smaller by using shorter data types
+11. Pass pointers to structs to functions
+12. Use structs instead of separated data types
+13. Remove second copy of large prime check for negated value
+14. ***Remove printf's from lanczos code
+15. Ensure polynomial chooser can't pick wrong A factors
+16. ***Count actual number of relations, not including duplicates
+17. ***Check more often near end of sieving if we are done
+18. ***Replace all errors with aborts
+19. Replace all file operations with safe ones
+20. Clean up after a run, i.e. free memory allocated
+21. Introduce hash tables for computers with large caches
+22. In candidate evaluation, multiply by inverses modulo a large prime instead of dividing out by each prime. (maybe not for version 0.99)
+23. Adjust 4X sieve eval code for factorisations under 40 digits
+24. Make reentrant and into a single callable function
+25. Make parallelisable
+26. Package output into factor tree
+27. Integrate into Pari
+28. Rename .cpp to .c and make compile with gcc
+29. Integrate into FLINT makefile
+30. Remove sievediv
+31. ***Use unique filenames
+32. Write README file
+
+a.
+b. Optimise for other architectures
+c. Comment code
+d.
+e.
+f.
+
+I. ***Add single prime variant
+II. Add double prime variant
+III. Optimize cache usage for larger factor bases
+IV. Implement SQUFOf, elliptic curve, pollard-brent, etc and make single factoring algorithm
+
+*** Already done.
+
diff --git a/TonelliShanks.cpp b/TonelliShanks.cpp
new file mode 100644
index 0000000..e845640
--- /dev/null
+++ b/TonelliShanks.cpp
@@ -0,0 +1,147 @@
+/*============================================================================
+ Copyright 2006 William Hart
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+============================================================================*/
+
+// -------------------------------------------------------
+//
+// TonelliShanks.cpp
+//
+// Provides Tonelli-Shanks square root mod p, and mod p^k
+//
+// -------------------------------------------------------
+
+#include <gmp.h>
+#include "TonelliShanks.h"
+#include "ModuloArith.h" //need multiplication mod p
+
+mpz_t two; //variables for sqrtmod
+mpz_t p1;
+mpz_t b;
+mpz_t g;
+mpz_t xsq;
+mpz_t mk;
+mpz_t bpow;
+mpz_t gpow;
+
+mpz_t inv; //variables for sqrtmodpow
+mpz_t tempsqpow;
+
+mpz_t pk; //variable for sqrtmodpk
+
+void TonelliInit(void)
+{
+ mpz_init(two);
+ mpz_init(p1);
+ mpz_init(b);
+ mpz_init(g);
+ mpz_init(xsq);
+ mpz_init(mk);
+ mpz_init(bpow);
+ mpz_init(gpow);
+
+ mpz_init(inv);
+ mpz_init(tempsqpow);
+
+ mpz_init(pk);
+
+ return;
+}
+
+int32_t sqrtmod(mpz_t asqrt, mpz_t a, mpz_t p)
+{
+ int32_t r,k,m;
+
+ if (mpz_kronecker(a,p)!=1)
+ {
+ mpz_set_ui(asqrt,0);
+ return 0; //return 0 if a is not a square mod p
+ }
+
+ mpz_set_ui(two,2);
+
+ mpz_sub_ui(p1,p,1);
+ r = mpz_remove(p1,p1,two);
+ mpz_powm(b,a,p1,p);
+ for (k=2; ;k++)
+ {
+ if (mpz_ui_kronecker(k,p) == -1) break;
+ }
+ mpz_set_ui(mk,k);
+ mpz_powm(g,mk,p1,p);
+ mpz_add_ui(p1,p1,1);
+ mpz_divexact_ui(p1,p1,2);
+ mpz_powm(xsq,a,p1,p);
+ if (!mpz_cmp_ui(b,1))
+ {
+ mpz_set(asqrt,xsq);
+
+ return 1;
+ }
+
+ while (mpz_cmp_ui(b,1))
+ {
+ mpz_set(bpow,b);
+ for (m=1; (m<=r-1) && (mpz_cmp_ui(bpow,1));m++)
+ {
+ mpz_powm_ui(bpow,bpow,2,p);
+ }
+ mpz_set(gpow,g);
+ for (int32_t i = 1;i<= r-m-1;i++)
+ {
+ mpz_powm_ui(gpow,gpow,2,p);
+ };
+ modmul(xsq,xsq,gpow,p);
+ mpz_powm_ui(gpow,gpow,2,p);
+ modmul(b,b,gpow,p);
+ mpz_set(gpow,g);
+ r = m;
+ }
+
+ mpz_set(asqrt,xsq);
+
+ return 1;
+}
+
+inline void sqrtmodpow(mpz_t res, mpz_t z, mpz_t a, mpz_t pk)
+{
+ mpz_mul_ui(inv,z,2);
+ mpz_invert(inv,inv,pk);
+ mpz_set(tempsqpow,a);
+ mpz_submul(tempsqpow,z,z);
+ mpz_fdiv_r(tempsqpow,tempsqpow,pk);
+ modmul(tempsqpow,tempsqpow,inv,pk);
+ mpz_add(tempsqpow,tempsqpow,z);
+ mpz_set(res,tempsqpow);
+
+ return;
+}
+
+void sqrtmodpk(mpz_t res, mpz_t z, mpz_t a, mpz_t p, int32_t k)
+{
+ mpz_set(res,z);
+ mpz_set(pk,p);
+ for (int32_t i = 2;i<=k;i++)
+ {
+ mpz_mul(pk,pk,p);
+ sqrtmodpow(res,res,a,pk);
+ }
+
+ return;
+}
diff --git a/TonelliShanks.h b/TonelliShanks.h
new file mode 100644
index 0000000..92048e2
--- /dev/null
+++ b/TonelliShanks.h
@@ -0,0 +1,45 @@
+/*============================================================================
+ Copyright 2006 William Hart
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+============================================================================*/
+
+#include <stdlib.h>
+
+// =====================================================================
+// INTERFACE:
+//
+// void TonelliInit(void)
+// - Initialises variables
+//
+// int sqrtmod(mpz_t asqrt, mpz_t a, mpz_t p)
+// - Tonelli-Shanks: sets asqrt to a square root of a modulo p
+// - Return: 0 if a is not a square mod p, 1 otherwise.
+//
+// void sqrtmodpk(mpz_t res, mpz_t z, mpz_t a, mpz_t p, int k)
+// - Given a square root z, of a mod p (from Tonelli-Shanks say)
+// - sets res to a square root of a mod p^k
+//
+// ========================================================================
+
+extern void TonelliInit(void);
+extern int32_t sqrtmod(mpz_t, mpz_t, mpz_t);
+extern inline void sqrtmodpow(mpz_t, mpz_t, mpz_t, mpz_t);
+extern void sqrtmodpk(mpz_t, mpz_t, mpz_t, mpz_t, int32_t);
+
+
diff --git a/gpl.txt b/gpl.txt
new file mode 100644
index 0000000..d511905
--- /dev/null
+++ b/gpl.txt
@@ -0,0 +1,339 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The licenses for most software are designed to take away your
+freedom to share and change it. By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users. This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it. (Some other Free Software Foundation software is covered by
+the GNU Lesser General Public License instead.) You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+ To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have. You must make sure that they, too, receive or can get the
+source code. And you must show them these terms so they know their
+rights.
+
+ We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+ Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software. If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+ Finally, any free program is threatened constantly by software
+patents. We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary. To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ GNU GENERAL PUBLIC LICENSE
+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+ 0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License. The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language. (Hereinafter, translation is included without limitation in
+the term "modification".) Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope. The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+ 1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+ 2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+ a) You must cause the modified files to carry prominent notices
+ stating that you changed the files and the date of any change.
+
+ b) You must cause any work that you distribute or publish, that in
+ whole or in part contains or is derived from the Program or any
+ part thereof, to be licensed as a whole at no charge to all third
+ parties under the terms of this License.
+
+ c) If the modified program normally reads commands interactively
+ when run, you must cause it, when started running for such
+ interactive use in the most ordinary way, to print or display an
+ announcement including an appropriate copyright notice and a
+ notice that there is no warranty (or else, saying that you provide
+ a warranty) and that users may redistribute the program under
+ these conditions, and telling the user how to view a copy of this
+ License. (Exception: if the Program itself is interactive but
+ does not normally print such an announcement, your work based on
+ the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole. If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works. But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+ 3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+ a) Accompany it with the complete corresponding machine-readable
+ source code, which must be distributed under the terms of Sections
+ 1 and 2 above on a medium customarily used for software interchange; or,
+
+ b) Accompany it with a written offer, valid for at least three
+ years, to give any third party, for a charge no more than your
+ cost of physically performing source distribution, a complete
+ machine-readable copy of the corresponding source code, to be
+ distributed under the terms of Sections 1 and 2 above on a medium
+ customarily used for software interchange; or,
+
+ c) Accompany it with the information you received as to the offer
+ to distribute corresponding source code. (This alternative is
+ allowed only for noncommercial distribution and only if you
+ received the program in object code or executable form with such
+ an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it. For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable. However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+ 4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License. Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+ 5. You are not required to accept this License, since you have not
+signed it. However, nothing else grants you permission to modify or
+distribute the Program or its derivative works. These actions are
+prohibited by law if you do not accept this License. Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+ 6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions. You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+ 7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all. For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices. Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+ 8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded. In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+ 9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number. If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation. If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+ 10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission. For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this. Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+ NO WARRANTY
+
+ 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+ 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+ <one line to give the program's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+ 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; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+ Gnomovision version 69, Copyright (C) year name of author
+ Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary. Here is a sample; alter the names:
+
+ Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+ `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+ <signature of Ty Coon>, 1 April 1989
+ Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs. If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library. If this is what you want to do, use the GNU Lesser General
+Public License instead of this License.
diff --git a/lanczos.c b/lanczos.c
new file mode 100644
index 0000000..891551a
--- /dev/null
+++ b/lanczos.c
@@ -0,0 +1,975 @@
+/*============================================================================
+ Copyright 2006 Jason Papadopoulos
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+===============================================================================
+
+Optionally, please be nice and tell me if you find this source to be
+useful. Again optionally, if you add to the functionality present here
+please consider making those additions public too, so that others may
+benefit from your work.
+ --jasonp at boo.net 9/8/06
+
+The following modifications were made by William Hart:
+ -addition of a random generator and max function
+ -added the utility function getNullEntry
+ -reformatted original code so it would operate as a standalone
+ filter and block Lanczos module
+--------------------------------------------------------------------*/
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "lanczos.h"
+
+#define NUM_EXTRA_RELATIONS 64
+
+#define BIT(x) (((u_int64_t)(1)) << (x))
+
+static const u_int64_t bitmask[64] = {
+ BIT( 0), BIT( 1), BIT( 2), BIT( 3), BIT( 4), BIT( 5), BIT( 6), BIT( 7),
+ BIT( 8), BIT( 9), BIT(10), BIT(11), BIT(12), BIT(13), BIT(14), BIT(15),
+ BIT(16), BIT(17), BIT(18), BIT(19), BIT(20), BIT(21), BIT(22), BIT(23),
+ BIT(24), BIT(25), BIT(26), BIT(27), BIT(28), BIT(29), BIT(30), BIT(31),
+ BIT(32), BIT(33), BIT(34), BIT(35), BIT(36), BIT(37), BIT(38), BIT(39),
+ BIT(40), BIT(41), BIT(42), BIT(43), BIT(44), BIT(45), BIT(46), BIT(47),
+ BIT(48), BIT(49), BIT(50), BIT(51), BIT(52), BIT(53), BIT(54), BIT(55),
+ BIT(56), BIT(57), BIT(58), BIT(59), BIT(60), BIT(61), BIT(62), BIT(63),
+};
+
+/*--------------------------------------------------------------------*/
+u_int64_t getNullEntry(u_int64_t * nullrows, long i, long l) {
+
+ /* Returns true if the entry with indices i,l is 1 in the
+ supplied 64xN matrix. This is used to read the nullspace
+ vectors which are output by the Lanczos routine */
+
+ return nullrows[i]&bitmask[l];
+}
+
+u_long random32(void) {
+
+ /* Poor man's random number generator. It satisfies no
+ particularly good randomness properties, but is good
+ enough for this application */
+
+ static unsigned long randval = 4035456057U;
+ randval = ((u_int64_t)randval*1025416097U+286824428U)%(u_int64_t)4294967291U;
+
+ return (unsigned long)randval;
+}
+
+unsigned long max(unsigned long a, unsigned long b) {
+
+ /* Returns the maximum of two unsigned long's */
+
+ return (a<b) ? b : a;
+}
+
+/*--------------------------------------------------------------------*/
+void reduce_matrix(unsigned long *nrows,
+ unsigned long *ncols, la_col_t *cols) {
+
+ /* Perform light filtering on the nrows x ncols
+ matrix specified by cols[]. The processing here is
+ limited to deleting columns that contain a singleton
+ row, then resizing the matrix to have a few more
+ columns than rows. Because deleting a column reduces
+ the counts in several different rows, the process
+ must iterate to convergence.
+
+ Note that this step is not intended to make the Lanczos
+ iteration run any faster (though it will); it's just
+ that if we don't go to this trouble then there are
+ factorizations for which the matrix step will fail
+ outright */
+
+ unsigned long r, c, i, j, k;
+ unsigned long passes;
+ unsigned long *counts;
+ unsigned long reduced_rows;
+ unsigned long reduced_cols;
+
+ /* count the number of nonzero entries in each row */
+
+ counts = (unsigned long *)calloc((size_t)*nrows, sizeof(unsigned long));
+ for (i = 0; i < *ncols; i++) {
+ for (j = 0; j < cols[i].weight; j++)
+ counts[cols[i].data[j]]++;
+ }
+
+ reduced_rows = *nrows;
+ reduced_cols = *ncols;
+ passes = 0;
+
+ do {
+ r = reduced_rows;
+
+ /* remove any columns that contain the only entry
+ in one or more rows, then update the row counts
+ to reflect the missing column. Iterate until
+ no more columns can be deleted */
+
+ do {
+ c = reduced_cols;
+ for (i = j = 0; i < reduced_cols; i++) {
+ la_col_t *col = cols + i;
+ for (k = 0; k < col->weight; k++) {
+ if (counts[col->data[k]] < 2)
+ break;
+ }
+
+ if (k < col->weight) {
+ for (k = 0; k < col->weight; k++) {
+ counts[col->data[k]]--;
+ }
+ free(col->data);
+ }
+ else {
+ cols[j++] = cols[i];
+ }
+ }
+ reduced_cols = j;
+ } while (c != reduced_cols);
+
+ /* count the number of rows that contain a
+ nonzero entry */
+
+ for (i = reduced_rows = 0; i < *nrows; i++) {
+ if (counts[i])
+ reduced_rows++;
+ }
+
+ /* Because deleting a column reduces the weight
+ of many rows, the number of nonzero rows may
+ be much less than the number of columns. Delete
+ more columns until the matrix has the correct
+ aspect ratio. Columns at the end of cols[] are
+ the heaviest, so delete those (and update the
+ row counts again) */
+
+ if (reduced_cols > reduced_rows + NUM_EXTRA_RELATIONS) {
+ for (i = reduced_rows + NUM_EXTRA_RELATIONS;
+ i < reduced_cols; i++) {
+
+ la_col_t *col = cols + i;
+ for (j = 0; j < col->weight; j++) {
+ counts[col->data[j]]--;
+ }
+ free(col->data);
+ }
+ reduced_cols = reduced_rows + NUM_EXTRA_RELATIONS;
+ }
+
+ /* if any columns were deleted in the previous step,
+ then the matrix is less dense and more columns
+ can be deleted; iterate until no further deletions
+ are possible */
+
+ passes++;
+
+ } while (r != reduced_rows);
+
+ printf("reduce to %ld x %ld in %ld passes\n",
+ reduced_rows, reduced_cols, passes);
+
+ free(counts);
+
+ /* record the final matrix size. Note that we can't touch
+ nrows because all the column data (and the sieving relations
+ that produced it) would have to be updated */
+
+ *ncols = reduced_cols;
+}
+
+/*-------------------------------------------------------------------*/
+static void mul_64x64_64x64(u_int64_t *a, u_int64_t *b, u_int64_t *c ) {
+
+ /* c[][] = x[][] * y[][], where all operands are 64 x 64
+ (i.e. contain 64 words of 64 bits each). The result
+ may overwrite a or b. */
+
+ u_int64_t ai, bj, accum;
+ u_int64_t tmp[64];
+ unsigned long i, j;
+
+ for (i = 0; i < 64; i++) {
+ j = 0;
+ accum = 0;
+ ai = a[i];
+
+ while (ai) {
+ bj = b[j];
+ if( ai & 1 )
+ accum ^= bj;
+ ai >>= 1;
+ j++;
+ }
+
+ tmp[i] = accum;
+ }
+ memcpy(c, tmp, sizeof(tmp));
+}
+
+/*-----------------------------------------------------------------------*/
+static void precompute_Nx64_64x64(u_int64_t *x, u_int64_t *c) {
+
+ /* Let x[][] be a 64 x 64 matrix in GF(2), represented
+ as 64 words of 64 bits each. Let c[][] be an 8 x 256
+ matrix of 64-bit words. This code fills c[][] with
+ a bunch of "partial matrix multiplies". For 0<=i<256,
+ the j_th row of c[][] contains the matrix product
+
+ ( i << (8*j) ) * x[][]
+
+ where the quantity in parentheses is considered a
+ 1 x 64 vector of elements in GF(2). The resulting
+ table can dramatically speed up matrix multiplies
+ by x[][]. */
+
+ u_int64_t accum, xk;
+ unsigned long i, j, k, index;
+
+ for (j = 0; j < 8; j++) {
+ for (i = 0; i < 256; i++) {
+ k = 0;
+ index = i;
+ accum = 0;
+ while (index) {
+ xk = x[k];
+ if (index & 1)
+ accum ^= xk;
+ index >>= 1;
+ k++;
+ }
+ c[i] = accum;
+ }
+
+ x += 8;
+ c += 256;
+ }
+}
+
+/*-------------------------------------------------------------------*/
+static void mul_Nx64_64x64_acc(u_int64_t *v, u_int64_t *x, u_int64_t *c,
+ u_int64_t *y, unsigned long n) {
+
+ /* let v[][] be a n x 64 matrix with elements in GF(2),
+ represented as an array of n 64-bit words. Let c[][]
+ be an 8 x 256 scratch matrix of 64-bit words.
+ This code multiplies v[][] by the 64x64 matrix
+ x[][], then XORs the n x 64 result into y[][] */
+
+ unsigned long i;
+ u_int64_t word;
+
+ precompute_Nx64_64x64(x, c);
+
+ for (i = 0; i < n; i++) {
+ word = v[i];
+ y[i] ^= c[ 0*256 + ((word>> 0) & 0xff) ]
+ ^ c[ 1*256 + ((word>> 8) & 0xff) ]
+ ^ c[ 2*256 + ((word>>16) & 0xff) ]
+ ^ c[ 3*256 + ((word>>24) & 0xff) ]
+ ^ c[ 4*256 + ((word>>32) & 0xff) ]
+ ^ c[ 5*256 + ((word>>40) & 0xff) ]
+ ^ c[ 6*256 + ((word>>48) & 0xff) ]
+ ^ c[ 7*256 + ((word>>56) ) ];
+ }
+}
+
+/*-------------------------------------------------------------------*/
+static void mul_64xN_Nx64(u_int64_t *x, u_int64_t *y,
+ u_int64_t *c, u_int64_t *xy, unsigned long n) {
+
+ /* Let x and y be n x 64 matrices. This routine computes
+ the 64 x 64 matrix xy[][] given by transpose(x) * y.
+ c[][] is a 256 x 8 scratch matrix of 64-bit words. */
+
+ unsigned long i;
+
+ memset(c, 0, 256 * 8 * sizeof(u_int64_t));
+ memset(xy, 0, 64 * sizeof(u_int64_t));
+
+ for (i = 0; i < n; i++) {
+ u_int64_t xi = x[i];
+ u_int64_t yi = y[i];
+ c[ 0*256 + ( xi & 0xff) ] ^= yi;
+ c[ 1*256 + ((xi >> 8) & 0xff) ] ^= yi;
+ c[ 2*256 + ((xi >> 16) & 0xff) ] ^= yi;
+ c[ 3*256 + ((xi >> 24) & 0xff) ] ^= yi;
+ c[ 4*256 + ((xi >> 32) & 0xff) ] ^= yi;
+ c[ 5*256 + ((xi >> 40) & 0xff) ] ^= yi;
+ c[ 6*256 + ((xi >> 48) & 0xff) ] ^= yi;
+ c[ 7*256 + ((xi >> 56) ) ] ^= yi;
+ }
+
+
+ for(i = 0; i < 8; i++) {
+
+ unsigned long j;
+ u_int64_t a0, a1, a2, a3, a4, a5, a6, a7;
+
+ a0 = a1 = a2 = a3 = 0;
+ a4 = a5 = a6 = a7 = 0;
+
+ for (j = 0; j < 256; j++) {
+ if ((j >> i) & 1) {
+ a0 ^= c[0*256 + j];
+ a1 ^= c[1*256 + j];
+ a2 ^= c[2*256 + j];
+ a3 ^= c[3*256 + j];
+ a4 ^= c[4*256 + j];
+ a5 ^= c[5*256 + j];
+ a6 ^= c[6*256 + j];
+ a7 ^= c[7*256 + j];
+ }
+ }
+
+ xy[ 0] = a0; xy[ 8] = a1; xy[16] = a2; xy[24] = a3;
+ xy[32] = a4; xy[40] = a5; xy[48] = a6; xy[56] = a7;
+ xy++;
+ }
+}
+
+/*-------------------------------------------------------------------*/
+static unsigned long find_nonsingular_sub(u_int64_t *t, unsigned long *s,
+ unsigned long *last_s, unsigned long last_dim,
+ u_int64_t *w) {
+
+ /* given a 64x64 matrix t[][] (i.e. sixty-four
+ 64-bit words) and a list of 'last_dim' column
+ indices enumerated in last_s[]:
+
+ - find a submatrix of t that is invertible
+ - invert it and copy to w[][]
+ - enumerate in s[] the columns represented in w[][] */
+
+ unsigned long i, j;
+ unsigned long dim;
+ unsigned long cols[64];
+ u_int64_t M[64][2];
+ u_int64_t mask, *row_i, *row_j;
+ u_int64_t m0, m1;
+
+ /* M = [t | I] for I the 64x64 identity matrix */
+
+ for (i = 0; i < 64; i++) {
+ M[i][0] = t[i];
+ M[i][1] = bitmask[i];
+ }
+
+ /* put the column indices from last_s[] into the
+ back of cols[], and copy to the beginning of cols[]
+ any column indices not in last_s[] */
+
+ mask = 0;
+ for (i = 0; i < last_dim; i++) {
+ cols[63 - i] = last_s[i];
+ mask |= bitmask[last_s[i]];
+ }
+ for (i = j = 0; i < 64; i++) {
+ if (!(mask & bitmask[i]))
+ cols[j++] = i;
+ }
+
+ /* compute the inverse of t[][] */
+
+ for (i = dim = 0; i < 64; i++) {
+
+ /* find the next pivot row and put in row i */
+
+ mask = bitmask[cols[i]];
+ row_i = M[cols[i]];
+
+ for (j = i; j < 64; j++) {
+ row_j = M[cols[j]];
+ if (row_j[0] & mask) {
+ m0 = row_j[0];
+ m1 = row_j[1];
+ row_j[0] = row_i[0];
+ row_j[1] = row_i[1];
+ row_i[0] = m0;
+ row_i[1] = m1;
+ break;
+ }
+ }
+
+ /* if a pivot row was found, eliminate the pivot
+ column from all other rows */
+
+ if (j < 64) {
+ for (j = 0; j < 64; j++) {
+ row_j = M[cols[j]];
+ if ((row_i != row_j) && (row_j[0] & mask)) {
+ row_j[0] ^= row_i[0];
+ row_j[1] ^= row_i[1];
+ }
+ }
+
+ /* add the pivot column to the list of
+ accepted columns */
+
+ s[dim++] = cols[i];
+ continue;
+ }
+
+ /* otherwise, use the right-hand half of M[]
+ to compensate for the absence of a pivot column */
+
+ for (j = i; j < 64; j++) {
+ row_j = M[cols[j]];
+ if (row_j[1] & mask) {
+ m0 = row_j[0];
+ m1 = row_j[1];
+ row_j[0] = row_i[0];
+ row_j[1] = row_i[1];
+ row_i[0] = m0;
+ row_i[1] = m1;
+ break;
+ }
+ }
+
+ if (j == 64) {
+#ifdef ERRORS
+ printf("lanczos error: submatrix "
+ "is not invertible\n");
+#endif
+ return 0;
+ }
+
+ /* eliminate the pivot column from the other rows
+ of the inverse */
+
+ for (j = 0; j < 64; j++) {
+ row_j = M[cols[j]];
+ if ((row_i != row_j) && (row_j[1] & mask)) {
+ row_j[0] ^= row_i[0];
+ row_j[1] ^= row_i[1];
+ }
+ }
+
+ /* wipe out the pivot row */
+
+ row_i[0] = row_i[1] = 0;
+ }
+
+ /* the right-hand half of M[] is the desired inverse */
+
+ for (i = 0; i < 64; i++)
+ w[i] = M[i][1];
+
+ /* The block Lanczos recurrence depends on all columns
+ of t[][] appearing in s[] and/or last_s[].
+ Verify that condition here */
+
+ mask = 0;
+ for (i = 0; i < dim; i++)
+ mask |= bitmask[s[i]];
+ for (i = 0; i < last_dim; i++)
+ mask |= bitmask[last_s[i]];
+
+ if (mask != (u_int64_t)(-1)) {
+#ifdef ERRORS
+ printf("lanczos error: not all columns used\n");
+#endif
+ return 0;
+ }
+
+ return dim;
+}
+
+/*-------------------------------------------------------------------*/
+void mul_MxN_Nx64(unsigned long vsize, unsigned long dense_rows,
+ unsigned long ncols, la_col_t *A,
+ u_int64_t *x, u_int64_t *b) {
+
+ /* Multiply the vector x[] by the matrix A (stored
+ columnwise) and put the result in b[]. vsize
+ refers to the number of u_int64_t's allocated for
+ x[] and b[]; vsize is probably different from ncols */
+
+ unsigned long i, j;
+
+ memset(b, 0, vsize * sizeof(u_int64_t));
+
+ for (i = 0; i < ncols; i++) {
+ la_col_t *col = A + i;
+ unsigned long *row_entries = col->data;
+ u_int64_t tmp = x[i];
+
+ for (j = 0; j < col->weight; j++) {
+ b[row_entries[j]] ^= tmp;
+ }
+ }
+
+ if (dense_rows) {
+ for (i = 0; i < ncols; i++) {
+ la_col_t *col = A + i;
+ unsigned long *row_entries = col->data + col->weight;
+ u_int64_t tmp = x[i];
+
+ for (j = 0; j < dense_rows; j++) {
+ if (row_entries[j / 32] &
+ ((unsigned long)1 << (j % 32))) {
+ b[j] ^= tmp;
+ }
+ }
+ }
+ }
+}
+
+/*-------------------------------------------------------------------*/
+void mul_trans_MxN_Nx64(unsigned long dense_rows, unsigned long ncols,
+ la_col_t *A, u_int64_t *x, u_int64_t *b) {
+
+ /* Multiply the vector x[] by the transpose of the
+ matrix A and put the result in b[]. Since A is stored
+ by columns, this is just a matrix-vector product */
+
+ unsigned long i, j;
+
+ for (i = 0; i < ncols; i++) {
+ la_col_t *col = A + i;
+ unsigned long *row_entries = col->data;
+ u_int64_t accum = 0;
+
+ for (j = 0; j < col->weight; j++) {
+ accum ^= x[row_entries[j]];
+ }
+ b[i] = accum;
+ }
+
+ if (dense_rows) {
+ for (i = 0; i < ncols; i++) {
+ la_col_t *col = A + i;
+ unsigned long *row_entries = col->data + col->weight;
+ u_int64_t accum = b[i];
+
+ for (j = 0; j < dense_rows; j++) {
+ if (row_entries[j / 32] &
+ ((unsigned long)1 << (j % 32))) {
+ accum ^= x[j];
+ }
+ }
+ b[i] = accum;
+ }
+ }
+}
+
+/*-----------------------------------------------------------------------*/
+static void transpose_vector(unsigned long ncols, u_int64_t *v, u_int64_t **trans) {
+
+ /* Hideously inefficent routine to transpose a
+ vector v[] of 64-bit words into a 2-D array
+ trans[][] of 64-bit words */
+
+ unsigned long i, j;
+ unsigned long col;
+ u_int64_t mask, word;
+
+ for (i = 0; i < ncols; i++) {
+ col = i / 64;
+ mask = bitmask[i % 64];
+ word = v[i];
+ j = 0;
+ while (word) {
+ if (word & 1)
+ trans[j][col] |= mask;
+ word = word >> 1;
+ j++;
+ }
+ }
+}
+
+/*-----------------------------------------------------------------------*/
+void combine_cols(unsigned long ncols,
+ u_int64_t *x, u_int64_t *v,
+ u_int64_t *ax, u_int64_t *av) {
+
+ /* Once the block Lanczos iteration has finished,
+ x[] and v[] will contain mostly nullspace vectors
+ between them, as well as possibly some columns
+ that are linear combinations of nullspace vectors.
+ Given vectors ax[] and av[] that are the result of
+ multiplying x[] and v[] by the matrix, this routine
+ will use Gauss elimination on the columns of [ax | av]
+ to find all of the linearly dependent columns. The
+ column operations needed to accomplish this are mir-
+ rored in [x | v] and the columns that are independent
+ are skipped. Finally, the dependent columns are copied
+ back into x[] and represent the nullspace vector output
+ of the block Lanczos code.
+
+ v[] and av[] can be NULL, in which case the elimination
+ process assumes 64 dependencies instead of 128 */
+
+ unsigned long i, j, k, bitpos, col, col_words, num_deps;
+ u_int64_t mask;
+ u_int64_t *matrix[128], *amatrix[128], *tmp;
+
+ num_deps = 128;
+ if (v == NULL || av == NULL)
+ num_deps = 64;
+
+ col_words = (ncols + 63) / 64;
+
+ for (i = 0; i < num_deps; i++) {
+ matrix[i] = (u_int64_t *)calloc((size_t)col_words,
+ sizeof(u_int64_t));
+ amatrix[i] = (u_int64_t *)calloc((size_t)col_words,
+ sizeof(u_int64_t));
+ }
+
+ /* operations on columns can more conveniently become
+ operations on rows if all the vectors are first
+ transposed */
+
+ transpose_vector(ncols, x, matrix);
+ transpose_vector(ncols, ax, amatrix);
+ if (num_deps == 128) {
+ transpose_vector(ncols, v, matrix + 64);
+ transpose_vector(ncols, av, amatrix + 64);
+ }
+
+ /* Keep eliminating rows until the unprocessed part
+ of amatrix[][] is all zero. The rows where this
+ happens correspond to linearly dependent vectors
+ in the nullspace */
+
+ for (i = bitpos = 0; i < num_deps && bitpos < ncols; bitpos++) {
+
+ /* find the next pivot row */
+
+ mask = bitmask[bitpos % 64];
+ col = bitpos / 64;
+ for (j = i; j < num_deps; j++) {
+ if (amatrix[j][col] & mask) {
+ tmp = matrix[i];
+ matrix[i] = matrix[j];
+ matrix[j] = tmp;
+ tmp = amatrix[i];
+ amatrix[i] = amatrix[j];
+ amatrix[j] = tmp;
+ break;
+ }
+ }
+ if (j == num_deps)
+ continue;
+
+ /* a pivot was found; eliminate it from the
+ remaining rows */
+
+ for (j++; j < num_deps; j++) {
+ if (amatrix[j][col] & mask) {
+
+ /* Note that the entire row, *not*
+ just the nonzero part of it, must
+ be eliminated; this is because the
+ corresponding (dense) row of matrix[][]
+ must have the same operation applied */
+
+ for (k = 0; k < col_words; k++) {
+ amatrix[j][k] ^= amatrix[i][k];
+ matrix[j][k] ^= matrix[i][k];
+ }
+ }
+ }
+ i++;
+ }
+
+ /* transpose rows i to 64 back into x[] */
+
+ for (j = 0; j < ncols; j++) {
+ u_int64_t word = 0;
+
+ col = j / 64;
+ mask = bitmask[j % 64];
+
+ for (k = i; k < 64; k++) {
+ if (matrix[k][col] & mask)
+ word |= bitmask[k];
+ }
+ x[j] = word;
+ }
+
+ for (i = 0; i < num_deps; i++) {
+ free(matrix[i]);
+ free(amatrix[i]);
+ }
+}
+
+/*-----------------------------------------------------------------------*/
+u_int64_t * block_lanczos(unsigned long nrows,
+ unsigned long dense_rows, unsigned long ncols, la_col_t *B) {
+
+ /* Solve Bx = 0 for some nonzero x; the computed
+ solution, containing up to 64 of these nullspace
+ vectors, is returned */
+
+ u_int64_t *vnext, *v[3], *x, *v0;
+ u_int64_t *winv[3];
+ u_int64_t *vt_a_v[2], *vt_a2_v[2];
+ u_int64_t *scratch;
+ u_int64_t *d, *e, *f, *f2;
+ u_int64_t *tmp;
+ unsigned long s[2][64];
+ unsigned long i, iter;
+ unsigned long n = ncols;
+ unsigned long dim0, dim1;
+ u_int64_t mask0, mask1;
+ unsigned long vsize;
+
+ /* allocate all of the size-n variables. Note that because
+ B has been preprocessed to ignore singleton rows, the
+ number of rows may really be less than nrows and may
+ be greater than ncols. vsize is the maximum of these
+ two numbers */
+
+ vsize = max(nrows, ncols);
+ v[0] = (u_int64_t *)malloc(vsize * sizeof(u_int64_t));
+ v[1] = (u_int64_t *)malloc(vsize * sizeof(u_int64_t));
+ v[2] = (u_int64_t *)malloc(vsize * sizeof(u_int64_t));
+ vnext = (u_int64_t *)malloc(vsize * sizeof(u_int64_t));
+ x = (u_int64_t *)malloc(vsize * sizeof(u_int64_t));
+ v0 = (u_int64_t *)malloc(vsize * sizeof(u_int64_t));
+ scratch = (u_int64_t *)malloc(max(vsize, 256 * 8) * sizeof(u_int64_t));
+
+ /* allocate all the 64x64 variables */
+
+ winv[0] = (u_int64_t *)malloc(64 * sizeof(u_int64_t));
+ winv[1] = (u_int64_t *)malloc(64 * sizeof(u_int64_t));
+ winv[2] = (u_int64_t *)malloc(64 * sizeof(u_int64_t));
+ vt_a_v[0] = (u_int64_t *)malloc(64 * sizeof(u_int64_t));
+ vt_a_v[1] = (u_int64_t *)malloc(64 * sizeof(u_int64_t));
+ vt_a2_v[0] = (u_int64_t *)malloc(64 * sizeof(u_int64_t));
+ vt_a2_v[1] = (u_int64_t *)malloc(64 * sizeof(u_int64_t));
+ d = (u_int64_t *)malloc(64 * sizeof(u_int64_t));
+ e = (u_int64_t *)malloc(64 * sizeof(u_int64_t));
+ f = (u_int64_t *)malloc(64 * sizeof(u_int64_t));
+ f2 = (u_int64_t *)malloc(64 * sizeof(u_int64_t));
+
+ /* The iterations computes v[0], vt_a_v[0],
+ vt_a2_v[0], s[0] and winv[0]. Subscripts larger
+ than zero represent past versions of these
+ quantities, which start off empty (except for
+ the past version of s[], which contains all
+ the column indices */
+
+ memset(v[1], 0, vsize * sizeof(u_int64_t));
+ memset(v[2], 0, vsize * sizeof(u_int64_t));
+ for (i = 0; i < 64; i++) {
+ s[1][i] = i;
+ vt_a_v[1][i] = 0;
+ vt_a2_v[1][i] = 0;
+ winv[1][i] = 0;
+ winv[2][i] = 0;
+ }
+ dim0 = 0;
+ dim1 = 64;
+ mask1 = (u_int64_t)(-1);
+ iter = 0;
+
+ /* The computed solution 'x' starts off random,
+ and v[0] starts off as B*x. This initial copy
+ of v[0] must be saved off separately */
+
+ for (i = 0; i < n; i++)
+ v[0][i] = (u_int64_t)(random32()) << 32 |
+ (u_int64_t)(random32());
+
+ memcpy(x, v[0], vsize * sizeof(u_int64_t));
+ mul_MxN_Nx64(vsize, dense_rows, ncols, B, v[0], scratch);
+ mul_trans_MxN_Nx64(dense_rows, ncols, B, scratch, v[0]);
+ memcpy(v0, v[0], vsize * sizeof(u_int64_t));
+
+ /* perform the iteration */
+
+ while (1) {
+ iter++;
+
+ /* multiply the current v[0] by a symmetrized
+ version of B, or B'B (apostrophe means
+ transpose). Use "A" to refer to B'B */
+
+ mul_MxN_Nx64(vsize, dense_rows, ncols, B, v[0], scratch);
+ mul_trans_MxN_Nx64(dense_rows, ncols, B, scratch, vnext);
+
+ /* compute v0'*A*v0 and (A*v0)'(A*v0) */
+
+ mul_64xN_Nx64(v[0], vnext, scratch, vt_a_v[0], n);
+ mul_64xN_Nx64(vnext, vnext, scratch, vt_a2_v[0], n);
+
+ /* if the former is orthogonal to itself, then
+ the iteration has finished */
+
+ for (i = 0; i < 64; i++) {
+ if (vt_a_v[0][i] != 0)
+ break;
+ }
+ if (i == 64) {
+ break;
+ }
+
+ /* Find the size-'dim0' nonsingular submatrix
+ of v0'*A*v0, invert it, and list the column
+ indices present in the submatrix */
+
+ dim0 = find_nonsingular_sub(vt_a_v[0], s[0],
+ s[1], dim1, winv[0]);
+ if (dim0 == 0)
+ break;
+
+ /* mask0 contains one set bit for every column
+ that participates in the inverted submatrix
+ computed above */
+
+ mask0 = 0;
+ for (i = 0; i < dim0; i++)
+ mask0 |= bitmask[s[0][i]];
+
+ /* compute d */
+
+ for (i = 0; i < 64; i++)
+ d[i] = (vt_a2_v[0][i] & mask0) ^ vt_a_v[0][i];
+
+ mul_64x64_64x64(winv[0], d, d);
+
+ for (i = 0; i < 64; i++)
+ d[i] = d[i] ^ bitmask[i];
+
+ /* compute e */
+
+ mul_64x64_64x64(winv[1], vt_a_v[0], e);
+
+ for (i = 0; i < 64; i++)
+ e[i] = e[i] & mask0;
+
+ /* compute f */
+
+ mul_64x64_64x64(vt_a_v[1], winv[1], f);
+
+ for (i = 0; i < 64; i++)
+ f[i] = f[i] ^ bitmask[i];
+
+ mul_64x64_64x64(winv[2], f, f);
+
+ for (i = 0; i < 64; i++)
+ f2[i] = ((vt_a2_v[1][i] & mask1) ^
+ vt_a_v[1][i]) & mask0;
+
+ mul_64x64_64x64(f, f2, f);
+
+ /* compute the next v */
+
+ for (i = 0; i < n; i++)
+ vnext[i] = vnext[i] & mask0;
+
+ mul_Nx64_64x64_acc(v[0], d, scratch, vnext, n);
+ mul_Nx64_64x64_acc(v[1], e, scratch, vnext, n);
+ mul_Nx64_64x64_acc(v[2], f, scratch, vnext, n);
+
+ /* update the computed solution 'x' */
+
+ mul_64xN_Nx64(v[0], v0, scratch, d, n);
+ mul_64x64_64x64(winv[0], d, d);
+ mul_Nx64_64x64_acc(v[0], d, scratch, x, n);
+
+ /* rotate all the variables */
+
+ tmp = v[2];
+ v[2] = v[1];
+ v[1] = v[0];
+ v[0] = vnext;
+ vnext = tmp;
+
+ tmp = winv[2];
+ winv[2] = winv[1];
+ winv[1] = winv[0];
+ winv[0] = tmp;
+
+ tmp = vt_a_v[1]; vt_a_v[1] = vt_a_v[0]; vt_a_v[0] = tmp;
+
+ tmp = vt_a2_v[1]; vt_a2_v[1] = vt_a2_v[0]; vt_a2_v[0] = tmp;
+
+ memcpy(s[1], s[0], 64 * sizeof(unsigned long));
+ mask1 = mask0;
+ dim1 = dim0;
+ }
+
+ printf("lanczos halted after %ld iterations\n", iter);
+
+ /* free unneeded storage */
+
+ free(vnext);
+ free(scratch);
+ free(v0);
+ free(vt_a_v[0]);
+ free(vt_a_v[1]);
+ free(vt_a2_v[0]);
+ free(vt_a2_v[1]);
+ free(winv[0]);
+ free(winv[1]);
+ free(winv[2]);
+ free(d);
+ free(e);
+ free(f);
+ free(f2);
+
+ /* if a recoverable failure occurred, start everything
+ over again */
+
+ if (dim0 == 0) {
+#ifdef ERRORS
+ printf("linear algebra failed; retrying...\n");
+#endif
+ free(x);
+ free(v[0]);
+ free(v[1]);
+ free(v[2]);
+ return NULL;
+ }
+
+ /* convert the output of the iteration to an actual
+ collection of nullspace vectors */
+
+ mul_MxN_Nx64(vsize, dense_rows, ncols, B, x, v[1]);
+ mul_MxN_Nx64(vsize, dense_rows, ncols, B, v[0], v[2]);
+
+ combine_cols(ncols, x, v[0], v[1], v[2]);
+
+ /* verify that these really are linear dependencies of B */
+
+ mul_MxN_Nx64(vsize, dense_rows, ncols, B, x, v[0]);
+
+ for (i = 0; i < ncols; i++) {
+ if (v[0][i] != 0)
+ break;
+ }
+ if (i < ncols) {
+ printf("lanczos error: dependencies don't work %ld\n",i);
+ abort();
+ }
+
+ free(v[0]);
+ free(v[1]);
+ free(v[2]);
+ return x;
+}
diff --git a/lanczos.h b/lanczos.h
new file mode 100644
index 0000000..c07e6d1
--- /dev/null
+++ b/lanczos.h
@@ -0,0 +1,95 @@
+/*============================================================================
+ Copyright 2006 William Hart
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+============================================================================*/
+#ifndef LANCZOS_H
+#define LANCZOS_H
+
+#include <stdlib.h>
+
+typedef struct {
+ unsigned long *data; /* The list of occupied rows in this column */
+ unsigned long weight; /* Number of nonzero entries in this column */
+ unsigned long orig; /* Original relation number */
+} la_col_t;
+
+u_int64_t getNullEntry(u_int64_t *, long, long);
+void reduce_matrix(unsigned long *, unsigned long *, la_col_t *);
+u_int64_t * block_lanczos(unsigned long, unsigned long, unsigned long, la_col_t*);
+
+/*==========================================================================
+ insertColEntry:
+
+ Function: insert an entry into a column of the matrix,
+ reallocating the space for the column if necessary
+
+===========================================================================*/
+static inline void insertColEntry(la_col_t* colarray, unsigned long colNum, unsigned long entry)
+{
+ unsigned long* temp;
+
+ if ((((colarray[colNum].weight)>>4)<<4)==colarray[colNum].weight) //need more space
+ {
+ temp = colarray[colNum].data;
+ colarray[colNum].data = (unsigned long*)malloc((colarray[colNum].weight+16)*sizeof(unsigned long));
+ for (long i = 0; i<colarray[colNum].weight; i++)
+ {
+ colarray[colNum].data[i] = temp[i];
+ }
+ if (colarray[colNum].weight!=0) free(temp);
+ }
+
+ colarray[colNum].data[colarray[colNum].weight] = entry;
+ colarray[colNum].weight++;
+ colarray[colNum].orig = colNum;
+}
+
+/*==========================================================================
+ xorColEntry:
+
+ Function: xor entry corresponding to a prime dividing A, which will be
+ either the last entry in the column, or not there at all, so we either
+ add it in or take it out
+
+===========================================================================*/
+static inline void xorColEntry(la_col_t* colarray, unsigned long colNum, unsigned long entry)
+{
+ for (long i = 0; i < colarray[colNum].weight; i++)
+ if (colarray[colNum].data[i] == entry)
+ {
+ for (unsigned long j = i; j < colarray[colNum].weight - 1; j++)
+ colarray[colNum].data[j] = colarray[colNum].data[j+1];
+ colarray[colNum].weight--;
+ return;
+ }
+ insertColEntry(colarray,colNum,entry);
+}
+
+/*==========================================================================
+ clearCol:
+
+ Function: clear a column
+
+===========================================================================*/
+static inline void clearCol(la_col_t* colarray, unsigned long colNum)
+{
+ colarray[colNum].weight =0;
+}
+
+#endif
diff --git a/lprels.c b/lprels.c
new file mode 100644
index 0000000..5b6613e
--- /dev/null
+++ b/lprels.c
@@ -0,0 +1,756 @@
+/*============================================================================
+ Copyright 2006 William Hart
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+===============================================================================*/
+
+/*
+ This code has been adapted for FLINT from mpqs.c in the Pari/GP package.
+ See http://pari.math.u-bordeaux.fr/
+*/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <gmp.h>
+#include <unistd.h>
+#include "lprels.h"
+
+#define min_bufspace 120UL /* use new buffer when < min_bufspace left */
+#define buflist_size 4096UL /* size of list-of-buffers blocks */
+#define sort_table_size 100000UL
+
+/*********************************************************************
+
+ File based large prime routines
+
+*********************************************************************/
+
+/*
+ Concatenates a filename and directory name to give a full path
+*/
+
+char * get_filename(char *dir, char *s)
+{
+ char *buf = (char *) malloc(strlen(dir) + strlen(s) + 2);
+#if defined(__EMX__) || defined(WINCE)
+ sprintf(buf, "%s\\%s", dir,s);
+#else
+ sprintf(buf, "%s/%s", dir,s);
+#endif
+ return buf;
+}
+
+char * unique_filename(char *s)
+{
+ char *buf, suf[64];
+ size_t lsuf;
+
+ sprintf(suf,".%ld.%ld", (long)getuid(), (long)getpid());
+
+ lsuf = strlen(suf);
+ /* room for s + suffix '\0' */
+ buf = (char*) malloc(8 + lsuf + 1);
+
+ sprintf(buf, "%.8s%s", s, suf);
+ return buf;
+}
+
+
+FILE * flint_fopen(char * name, char * mode)
+{
+#if defined(WINCE) || defined(macintosh)
+ char * tmp_dir = NULL;
+#else
+ char * tmp_dir = getenv("TMPDIR");
+#endif
+ if (tmp_dir == NULL) tmp_dir = "./";
+ FILE * temp_file = fopen(get_filename(tmp_dir,unique_filename(name)),mode);
+ if (!temp_file)
+ {
+ printf("Unable to open temporary file\n");
+ abort();
+ }
+ return temp_file;
+}
+
+/*
+ Compares two large prime relations according to their first
+ element (the large prime). Used by qsort.
+*/
+
+int relations_cmp(const void *a, const void *b)
+{
+ char **sa = (char**) a;
+ char **sb = (char**) b;
+ long qa = strtol(*sa, NULL, 10);
+ long qb = strtol(*sb, NULL, 10);
+
+ if (qa < qb) return -1;
+ else if (qa > qb) return 1;
+ else return strcmp(*sa, *sb);
+}
+
+
+/*
+ Writes the given string to the given file and aborts upon error
+*/
+
+void flint_fputs(char *s, FILE *file)
+{
+ if (fputs(s, file) < 0)
+ {
+ printf("Error whilst writing to large prime file!");
+ abort();
+ }
+}
+
+/*
+ Given a file "filename" containing full or large prime relations,
+ rearrange the file so that relations are sorted by their first elements.
+ Works in memory, discards duplicate lines, and overwrites the original
+ file. Returns the number of relations after sorting and discarding.
+*/
+
+long sort_lp_file(char *filename)
+{
+ FILE *TMP;
+ char *old_s, *buf, *cur_line;
+ char **s_table, **sort_table, **buflist, **buflist_head;
+ long i, j, count;
+ size_t length, bufspace;
+
+ buflist_head = (char**) malloc(buflist_size * sizeof(char*));
+ buflist = buflist_head;
+
+ *buflist++ = NULL; /* flag this as last and only buflist block */
+
+ TMP = flint_fopen(filename, "r");
+
+ /* allocate first buffer and read first line, if any, into it */
+ buf = (char*) malloc(MPQS_STRING_LENGTH * sizeof(char));
+ cur_line = buf;
+ bufspace = MPQS_STRING_LENGTH;
+
+ if (fgets(cur_line, bufspace, TMP) == NULL)
+ { /* file empty */
+ free(buf); free(buflist_head);
+ fclose(TMP);
+ return 0;
+ }
+ /* enter first buffer into buflist */
+ *buflist++ = buf; /* can't overflow the buflist block */
+ length = strlen(cur_line) + 1; /* count the \0 byte as well */
+ bufspace -= length;
+
+ s_table = (char**) malloc(sort_table_size*sizeof(char*));
+ sort_table = s_table+sort_table_size;
+
+ /* at start of loop, one line from the file is sitting in cur_line inside buf,
+ the next will go into cur_line + length, and there's room for bufspace
+ further characters in buf. The loop reads another line if one exists, and
+ if this overruns the current buffer, it allocates a fresh one --GN */
+
+ for (i = 0, sort_table--; /* until end of file */; i++, sort_table--)
+ {
+
+ *sort_table = cur_line;
+ cur_line += length;
+
+ /* if little room is left, allocate a fresh buffer before attempting to
+ * read a line, and remember to free it if no further line is forthcoming.
+ * This avoids some copying of partial lines --GN */
+ if (bufspace < min_bufspace)
+ {
+ buf = (char*) malloc(MPQS_STRING_LENGTH * sizeof(char));
+ cur_line = buf;
+ bufspace = MPQS_STRING_LENGTH;
+ if (fgets(cur_line, bufspace, TMP) == NULL) { free(buf); break; }
+
+ if (buflist - buflist_head >= buflist_size) abort();
+
+ /* remember buffer for later deallocation */
+ *buflist++ = buf;
+ length = strlen(cur_line) + 1;
+ bufspace -= length; continue;
+ }
+
+ /* normal case: try fitting another line into the current buffer */
+ if (fgets(cur_line, bufspace, TMP) == NULL) break; /* none exists */
+ length = strlen(cur_line) + 1;
+ bufspace -= length;
+
+ /* check whether we got the entire line or only part of it */
+ if (bufspace == 0 && cur_line[length-2] != '\n')
+ {
+ size_t lg1;
+ buf = (char*) malloc(MPQS_STRING_LENGTH * sizeof(char));
+ if (buflist - buflist_head >= buflist_size) abort();
+ /* remember buffer for later deallocation */
+ *buflist++ = buf;
+
+ /* copy what we've got to the new buffer */
+ (void)strcpy(buf, cur_line); /* cannot overflow */
+ cur_line = buf + length - 1; /* point at the \0 byte */
+ bufspace = MPQS_STRING_LENGTH - length + 1;
+
+ /* read remainder of line */
+ if (fgets(cur_line, bufspace, TMP) == NULL)
+ {
+ printf("MPQS: relations file truncated?!\n");
+ abort();
+ }
+ lg1 = strlen(cur_line);
+ length += lg1; /* we already counted the \0 once */
+ bufspace -= (lg1 + 1); /* but here we must take it into account */
+ cur_line = buf; /* back up to the beginning of the line */
+ }
+ } /* for */
+
+ fclose(TMP);
+
+ /* sort the whole lot in place by swapping pointers */
+ qsort(sort_table, i, sizeof(char*), relations_cmp);
+
+ /* copy results back to the original file, skipping exact duplicates */
+ TMP = flint_fopen(filename, "w");
+ old_s = sort_table[0];
+ flint_fputs(sort_table[0], TMP);
+ count = 1;
+ for(j = 1; j < i; j++)
+ {
+ if (strcmp(old_s, sort_table[j]))
+ {
+ flint_fputs(sort_table[j], TMP);
+ count++;
+ }
+ old_s = sort_table[j];
+ }
+ fflush(TMP);
+ fclose(TMP);
+ /* deallocate buffers */
+ while (*--buflist)
+ {
+ if (buflist != buflist_head)
+ free(*buflist); /* free a buffer */
+ }
+ free(buflist_head);
+ free(s_table);
+ return count;
+}
+
+/*
+ Appends contents of file fp1 to fp (auxiliary routine for merge sort) and
+ returns number of lines copied. Closes fp afterwards.
+*/
+
+long append_file(FILE *fp, FILE *fp1)
+{
+ char line[MPQS_STRING_LENGTH];
+ long c = 0;
+ while (fgets(line, MPQS_STRING_LENGTH, fp1)) { flint_fputs(line, fp); c++; }
+ if (fflush(fp))
+ {
+ printf("Error while flushing file.\n");
+ abort();
+ }
+ fclose(fp); return c;
+}
+
+/*
+ Merge-sort on the files LPREL and LPNEW; assumes that LPREL and LPNEW are
+ already sorted. Creates/truncates the TMP file, writes result to it and
+ closes it (via append_file()). Instead of LPREL, LPNEW we may also call
+ this with FREL, FNEW. In the latter case COMB should be NULL (and we
+ return the count of all full relations), in the former case it should be
+ non-NULL (and we return the count of frels we expect to be able to combine
+ out of the present lprels). If COMB is non-NULL, the combinable lprels
+ are written out to this separate file.
+
+ We retain only one occurrence of each large prime in TMP (i.e. in the
+ future LPREL file). --GN
+*/
+
+#define swap_lines() { char *line_tmp;\
+ line_tmp = line_new_old; \
+ line_new_old = line_new; \
+ line_new = line_tmp; }
+
+long mergesort_lp_file_internal(FILE *LPREL, FILE *LPNEW, FILE *COMB, FILE *TMP)
+{
+ char line1[MPQS_STRING_LENGTH], line2[MPQS_STRING_LENGTH];
+ char line[MPQS_STRING_LENGTH];
+ char *line_new = line1, *line_new_old = line2;
+ long q_new, q_new_old = -1, q, i = 0, c = 0;
+ long comb_in_progress;
+
+ if ( !fgets(line_new, MPQS_STRING_LENGTH, LPNEW) )
+ {
+ /* LPNEW is empty: copy LPREL to TMP. Could be done by a rename if we
+ didn't want to count the lines (again)... however, this case will not
+ normally happen */
+ i = append_file(TMP, LPREL);
+ return COMB ? 0 : i;
+ }
+ /* we now have a line_new from LPNEW */
+
+ if (!fgets(line, MPQS_STRING_LENGTH, LPREL))
+ {
+ /* LPREL is empty: copy LPNEW to TMP... almost. */
+ flint_fputs(line_new, TMP);
+
+ if (!COMB)
+ {
+ /* full relations mode */
+ i = append_file(TMP, LPNEW);
+ return i + 1;
+ }
+
+ /* LP mode: check for combinable relations */
+ q_new_old = atol(line_new);
+ /* we need to retain a copy of the old line just for a moment, because we
+ may yet have to write it to COMB. Do this by swapping the two buffers */
+ swap_lines();
+ comb_in_progress = 0;
+ i = 0;
+
+ while (fgets(line_new, MPQS_STRING_LENGTH, LPNEW))
+ {
+ q_new = atol(line_new);
+ if (q_new_old == q_new)
+ {
+ /* found combinables, check whether we're already busy on this
+ particular large prime */
+ if (!comb_in_progress)
+ {
+ /* if not, write first line to COMB, creating and opening the
+ file first if it isn't open yet */
+ flint_fputs(line_new_old, COMB);
+ comb_in_progress = 1;
+ }
+ /* in any case, write the current line, and count it */
+ flint_fputs(line_new, COMB);
+ i++;
+ }
+ else
+ {
+ /* not combinable */
+ q_new_old = q_new;
+ comb_in_progress = 0;
+ /* and dump it to the TMP file */
+ flint_fputs(line_new, TMP);
+ /* and stash it away for a moment */
+ swap_lines();
+ comb_in_progress = 0;
+ }
+ } /* while */
+ fclose(TMP); return i;
+ }
+
+ /* normal case: both LPNEW and LPREL are not empty */
+ q_new = atol(line_new);
+ q = atol(line);
+
+ for(;;)
+ {
+ /* main merging loop */
+ i = comb_in_progress = 0;
+
+ /* first the harder case: let LPNEW catch up with LPREL, and possibly
+ overtake it, checking for combinables coming from LPNEW alone */
+ while (q > q_new)
+ {
+ if (!COMB || !comb_in_progress) flint_fputs(line_new, TMP);
+ if (!COMB) c++; /* in FREL mode, count lines written */
+ else if (!comb_in_progress)
+ {
+ q_new_old = q_new;
+ swap_lines();
+ }
+ if (!fgets(line_new, MPQS_STRING_LENGTH, LPNEW))
+ {
+ flint_fputs(line, TMP);
+ if (!COMB) c++; else c += i;
+ i = append_file(TMP, LPREL);
+ return COMB ? c : c + i;
+ }
+ q_new = atol(line_new);
+ if (!COMB) continue;
+
+ /* LP mode only: */
+ if (q_new_old != q_new) /* not combinable */
+ comb_in_progress = 0; /* next loop will deal with it, or loop may end */
+ else
+ {
+ /* found combinables, check whether we're already busy on this
+ large prime */
+ if (!comb_in_progress)
+ {
+ flint_fputs(line_new_old, COMB);
+ comb_in_progress = 1;
+ }
+ /* in any case, write the current line, and count it */
+ flint_fputs(line_new, COMB);
+ i++;
+ }
+ } /* while q > q_new */
+
+ /* q <= q_new */
+
+ if (COMB) c += i; /* accumulate count of combinables */
+ i = 0; /* and clear it */
+ comb_in_progress = 0;/* redundant */
+
+ /* now let LPREL catch up with LPNEW, and possibly overtake it */
+ while (q < q_new)
+ {
+ flint_fputs(line, TMP);
+ if (!COMB) c++;
+ if (!fgets(line, MPQS_STRING_LENGTH, LPREL))
+ {
+ flint_fputs(line_new, TMP);
+ i = append_file(TMP, LPNEW);
+ return COMB ? c : c + i + 1;
+ }
+ else
+ q = atol(line);
+ }
+
+ /* q >= q_new */
+
+ /* Finally, it may happen that q == q_new, indicating combinables whose
+ large prime is already in LPREL, and appears now one or more times in
+ LPNEW. Thus in this sub-loop we advance LPNEW. The `line' from LPREL is
+ left alone, and will be written to TMP the next time around the main for
+ loop; we only write it to COMB here -- unless all we find is an exact
+ duplicate of the line we already have, that is. (There can be at most
+ one such, and if so it is simply discarded.) */
+ while (q == q_new)
+ {
+ if (!strcmp(line_new, line))
+ {
+ /* duplicate -- move right ahead to the next LPNEW line */
+ ;/* do nothing here */
+ }
+ else if (!COMB)
+ { /* full relations mode: write line_new out first, keep line */
+ flint_fputs(line_new, TMP);
+ c++;
+ }
+ else
+ {
+ /* LP mode, and combinable relation */
+ if (!comb_in_progress)
+ {
+ flint_fputs(line, COMB);
+ comb_in_progress = 1;
+ }
+ flint_fputs(line_new, COMB);
+ i++;
+ }
+ /* NB comb_in_progress is cleared by q_new becoming bigger than q, thus
+ the current while loop terminating, the next time through the main for
+ loop */
+
+ /* common ending: get another line_new, if any */
+ if (!fgets(line_new, MPQS_STRING_LENGTH, LPNEW))
+ {
+ flint_fputs(line, TMP);
+ if (!COMB) c++; else c += i;
+ i = append_file(TMP, LPREL);
+ return COMB ? c : c + i;
+ }
+ else
+ q_new = atol(line_new);
+ } /* while */
+
+ if (COMB) c += i; /* accumulate count of combinables */
+ }
+}
+
+/*
+ Perform mergesort of large prime files
+*/
+
+long mergesort_lp_file(char *REL_str, char *NEW_str, char *TMP_str, FILE *COMB)
+{
+ FILE *NEW = flint_fopen(NEW_str, "r");
+
+#if defined(WINCE) || defined(macintosh)
+ char * tmp_dir = NULL;
+#else
+ char * tmp_dir = getenv("TMPDIR");
+#endif
+ if (tmp_dir == NULL) tmp_dir = "./";
+ char * TMP_name = get_filename(tmp_dir,unique_filename(TMP_str));
+ char * REL_name = get_filename(tmp_dir,unique_filename(REL_str));
+ FILE * TMP = fopen(TMP_name,"w");
+ FILE * REL = fopen(REL_name,"r");
+ if ((!TMP) || (!REL))
+ {
+ printf("Unable to open temporary file\n");
+ abort();
+ }
+
+ long tp = mergesort_lp_file_internal(REL, NEW, COMB, TMP);
+ fclose(REL);
+ fclose(NEW);
+
+ if (rename(TMP_name,REL_name))
+ {
+ printf("Cannot rename file %s to %s", TMP_str, REL_str);
+ abort();
+ }
+ return tp;
+}
+
+void read_matrix(unsigned long ** relations, FILE *FREL, la_col_t* colarray, unsigned long * relsFound, unsigned long relSought, mpz_t * XArr, mpz_t n, unsigned long * factorBase)
+{
+ long e, p;
+ char buf[MPQS_STRING_LENGTH], *s;
+ //char buf2[MPQS_STRING_LENGTH];
+ unsigned long numfactors;
+ mpz_t test1, test2;
+ mpz_init(test1);
+ mpz_init(test2);
+
+
+ if (ftell(FREL) < 0)
+ {
+ printf("Error on full relations file\n");
+ abort();
+ }
+ while ((fgets(buf, MPQS_STRING_LENGTH, FREL)) && ((*relsFound) < relSought))
+ {
+ numfactors = 0;
+ gmp_sscanf(buf,"%Zd",XArr[*relsFound]);
+ s = strchr(buf, ':') + 2;
+ s = strtok(s, " \n");
+ while (s != NULL)
+ {
+ e = atol(s); if (!e) break;
+ s = strtok(NULL, " \n");
+ p = atol(s);
+ if (e & 1) xorColEntry(colarray,*relsFound,p);
+ for (long i = 0; i < e; i++) relations[*relsFound][++numfactors] = p;
+ s = strtok(NULL, " \n");
+ }
+ relations[*relsFound][0] = numfactors;
+
+ mpz_set_ui(test1,1);
+ for (unsigned long i = 1; i<=relations[*relsFound][0]; i++)
+ {
+ mpz_mul_ui(test1,test1,factorBase[relations[*relsFound][i]]);
+ if (i%30 == 0) mpz_mod(test1,test1,n);
+ }
+ mpz_mod(test1,test1,n);
+ mpz_mul(test2,XArr[*relsFound],XArr[*relsFound]);
+ mpz_mod(test2,test2,n);
+ if (mpz_cmp(test1,test2)!=0)
+ {
+ mpz_add(test1,test1,test2);
+ if (mpz_cmp(test1,n)!=0)
+ {
+ clearCol(colarray,*relsFound);
+ }
+ else (*relsFound)++;
+ } else (*relsFound)++;
+ }
+
+ mpz_clear(test1);
+ mpz_clear(test2);
+
+ return;
+}
+
+/*********************************************************************
+
+ Routines for writing relations as strings
+
+*********************************************************************/
+
+/*
+ Writes a factor pi^ei into a string as " ei pi"
+*/
+
+void add_factor(char **last, unsigned long ei, unsigned long pi)
+{
+ sprintf(*last, " %ld %ld", ei, pi);
+ *last += strlen(*last);
+}
+
+/*
+ Concatenate " 0" to string
+*/
+
+void add_0(char **last)
+{
+ char *s = *last;
+ *s++ = ' ';
+ *s++ = '0';
+ *s++ = 0; *last = s;
+}
+
+/*********************************************************************
+
+ Large prime relation combining
+
+*********************************************************************/
+
+/*
+ Add to an array of unsigned longs the exponents from a relation string
+*/
+
+void set_exponents(unsigned long *ei, char *r)
+{
+ char *s, b[MPQS_STRING_LENGTH];
+ long e;
+
+ strcpy(b, r);
+ s = strtok(b, " \n");
+ while (s != NULL)
+ {
+ e = atol(s); if (!e) break;
+ s = strtok(NULL, " \n");
+ ei[atol(s)] += e;
+ s = strtok(NULL, " \n");
+ }
+}
+
+/* Writes an lp_entry from a string */
+
+void set_lp_entry(mpqs_lp_entry *e, char *buf)
+{
+ char *s1, *s2;
+ s1 = buf; s2 = strchr(s1, ' '); *s2 = '\0';
+ e->q = atol(s1);
+ s1 = s2 + 3; s2 = strchr(s1, ' '); *s2 = '\0';
+ strcpy(e->Y, s1);
+ s1 = s2 + 3; s2 = strchr(s1, '\n'); *s2 = '\0';
+ strcpy(e->E, s1);
+}
+
+/*
+ Combines the large prime relations in COMB to full relations in FNEW.
+ FNEW is assumed to be open for writing / appending.
+*/
+
+long combine_large_primes(unsigned long numPrimes,
+ FILE *COMB, FILE *FNEW, mpz_t N, mpz_t factor)
+{
+ char new_relation[MPQS_STRING_LENGTH], buf[MPQS_STRING_LENGTH];
+ mpqs_lp_entry e[2]; /* we'll use the two alternatingly */
+ unsigned long *ei;
+ long ei_size = numPrimes;
+ long old_q;
+ mpz_t inv_q, Y1, Y2, new_Y, new_Y1;
+ mpz_init(inv_q);mpz_init(Y1);mpz_init(Y2);mpz_init(new_Y);mpz_init(new_Y1);
+ long i, l, c = 0;
+
+ if (!fgets(buf, MPQS_STRING_LENGTH, COMB)) return 0; /* should not happen */
+
+ ei = (unsigned long *) malloc(sizeof(unsigned long)*ei_size);
+
+ /* put first lp relation in row 0 of e */
+ set_lp_entry(&e[0], buf);
+
+ i = 1; /* second relation will go into row 1 */
+ old_q = e[0].q;
+ mpz_set_ui(inv_q, old_q);
+ while (!mpz_invert(inv_q, inv_q, N)) /* can happen */
+ {
+ /* We have found a factor. It could be N when N is quite small;
+ or we might just have found a divisor by sheer luck. */
+ mpz_gcd_ui(inv_q, N, old_q);
+ if (!mpz_cmp(inv_q, N)) /* pity */
+ {
+ if (!fgets(buf, MPQS_STRING_LENGTH, COMB)) { return 0; }
+ set_lp_entry(&e[0], buf);
+ old_q = e[0].q;
+ mpz_set_ui(inv_q, old_q);
+ continue;
+ }
+ mpz_set(factor, inv_q);
+ free(ei);
+ return c;
+ }
+ gmp_sscanf(e[0].Y, "%Zd", Y1);
+
+ while (fgets(buf, MPQS_STRING_LENGTH, COMB))
+ {
+ set_lp_entry(&e[i], buf);
+ if (e[i].q != old_q)
+ {
+ /* switch to combining a new bunch, swapping the rows */
+ old_q = e[i].q;
+ mpz_set_ui(inv_q, old_q);
+ while (!mpz_invert(inv_q, inv_q, N)) /* can happen */
+ {
+ mpz_gcd_ui(inv_q, N, old_q);
+ if (!mpz_cmp(inv_q, N)) /* pity */
+ {
+ old_q = -1; /* sentinel */
+ continue; /* discard this combination */
+ }
+ mpz_set(factor, inv_q);
+ free(ei);
+ return c;
+ }
+ gmp_sscanf(e[i].Y, "%Zd", Y1);
+ i = 1 - i; /* subsequent relations go to other row */
+ continue;
+ }
+ /* count and combine the two we've got, and continue in the same row */
+ memset((void *)ei, 0, ei_size * sizeof(long));
+ set_exponents(ei, e[0].E);
+ set_exponents(ei, e[1].E);
+ gmp_sscanf(e[i].Y, "%Zd", Y2);
+
+ if (mpz_cmpabs(Y1,Y2)!=0)
+ {
+ c++;
+ mpz_mul(new_Y, Y1, Y2);
+ mpz_mul(new_Y, new_Y, inv_q);
+ mpz_mod(new_Y, new_Y, N);
+
+ mpz_sub(new_Y1, N, new_Y);
+ if (mpz_cmpabs(new_Y1, new_Y) < 0) mpz_set(new_Y, new_Y1);
+
+ gmp_sprintf(buf, "%Zd\0", new_Y);
+ strcpy(new_relation, buf);
+ strcat(new_relation, " :");
+ for (l = 0; l < ei_size; l++)
+ if (ei[l])
+ {
+ sprintf(buf, " %ld %ld", ei[l], l);
+ strcat(new_relation, buf);
+ }
+ strcat(new_relation, " 0");
+ strcat(new_relation, "\n");
+
+ flint_fputs(new_relation, FNEW);
+ }
+ } /* while */
+
+ free(ei);
+ mpz_clear(inv_q);mpz_clear(Y1);mpz_clear(Y2);mpz_clear(new_Y);mpz_clear(new_Y1);
+
+ return c;
+}
+
+
diff --git a/lprels.h b/lprels.h
new file mode 100644
index 0000000..8c43b3e
--- /dev/null
+++ b/lprels.h
@@ -0,0 +1,51 @@
+/*============================================================================
+ Copyright 2006 William Hart
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+===============================================================================*/
+
+#ifndef LPRELS_H
+#define LPRELS_H
+
+#include "lanczos.h"
+
+#define MPQS_STRING_LENGTH (4 * 1024UL)
+
+typedef struct {
+ long q;
+ char Y[MPQS_STRING_LENGTH];
+ char E[MPQS_STRING_LENGTH];
+} mpqs_lp_entry;
+
+char * get_filename(char *dir, char *s);
+int mpqs_relations_cmp(const void *a, const void *b);
+void flint_fputs(char *s, FILE *file);
+long sort_lp_file(char *filename);
+long append_file(FILE *fp, FILE *fp1);
+long mpqs_mergesort_lp_file_internal(FILE *LPREL, FILE *LPNEW, FILE *COMB, FILE *TMP);
+long mergesort_lp_file(char *REL_str, char *NEW_str, char *TMP_str, FILE *COMB);
+void add_factor(char **last, unsigned long ei, unsigned long pi);
+void add_0(char **last);
+void set_exponents(unsigned long *ei, char *r);
+void set_lp_entry(mpqs_lp_entry *e, char *buf);
+long combine_large_primes(unsigned long numPrimes, FILE *COMB, FILE *FNEW, mpz_t N, mpz_t factor);
+void read_matrix(unsigned long ** relations, FILE *FREL, la_col_t* colarray, unsigned long * relsFound, unsigned long relSought, mpz_t * XArr, mpz_t n, unsigned long * factorBase);
+FILE * flint_fopen(char * name, char * mode);
+char * unique_filename(char *s);
+
+#endif
diff --git a/makefile b/makefile
new file mode 100644
index 0000000..0d92518
--- /dev/null
+++ b/makefile
@@ -0,0 +1,64 @@
+#============================================================================
+# Copyright 2006 William Hart
+#
+# This file is part of FLINT.
+#
+# FLINT 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.
+#
+# FLINT 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 FLINT; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+#
+#============================================================================
+
+CPP = g++
+OBJ = lprels.o ModuloArith.o TonelliShanks.o F2matrix.o lanczos.o QS.o $(RES)
+LINKOBJ = lprels.o ModuloArith.o TonelliShanks.o F2matrix.o lanczos.o QS.o $(RES)
+#LIBS = -L"home/dmharvey/gmp/install/lib" -lgmp
+#CXXINCS = -I"home/dmharvey/gmp/install/include"
+LIBS = -L"/usr/lib" -lgmp
+CXXINCS = -I"/usr/include"
+BIN = QuadraticSieve QuadraticSieve.exe
+CXXFLAGS = $(CXXINCS) -Wall -Wno-sign-compare -fomit-frame-pointer -O2
+CXXFLAGS2 = $(CXXINCS) -Wall -Wno-sign-compare -fomit-frame-pointer -O3
+#CXXFLAGS = $(CXXINCS) -ansi -Wall -Wno-sign-compare -march=athlon-xp -fomit-frame-pointer -O2
+#CXXFLAGS2 = $(CXXINCS) -ansi -Wall -Wno-sign-compare -march=athlon-xp -fomit-frame-pointer -O3
+#CXXFLAGS = $(CXXINCS) -Wall -Wno-sign-compare -march=opteron -fomit-frame-pointer -O2
+#CXXFLAGS2 = $(CXXINCS) -Wall -Wno-sign-compare -march=opteron -fomit-frame-pointer -O3
+RM = rm -f
+
+.PHONY: all clean clean-custom
+
+all: QuadraticSieve
+
+clean: clean-custom
+ ${RM} $(OBJ) $(BIN)
+
+$(BIN): $(OBJ)
+ $(CPP) -ansi $(LINKOBJ) -o "QuadraticSieve" $(LIBS)
+
+ModuloArith.o: ModuloArith.cpp
+ $(CPP) -ansi -c ModuloArith.cpp -o ModuloArith.o $(CXXFLAGS)
+
+TonelliShanks.o: TonelliShanks.cpp
+ $(CPP) -ansi -c TonelliShanks.cpp -o TonelliShanks.o $(CXXFLAGS)
+
+F2matrix.o: F2matrix.cpp
+ $(CPP) -ansi -c F2matrix.cpp -o F2matrix.o $(CXXFLAGS2)
+
+lanczos.o: lanczos.c
+ $(CPP) -ansi -c lanczos.c -o lanczos.o $(CXXFLAGS2)
+
+lprels.o: lprels.c
+ $(CPP) -ansi -c lprels.c -o lprels.o $(CXXFLAGS2)
+
+QS.o: QS.cpp
+ $(CPP) -ansi -c QS.cpp -o QS.o $(CXXFLAGS)
diff --git a/makefile.opteron b/makefile.opteron
new file mode 100644
index 0000000..94ef78b
--- /dev/null
+++ b/makefile.opteron
@@ -0,0 +1,60 @@
+#============================================================================
+# Copyright 2006 William Hart
+#
+# This file is part of FLINT.
+#
+# FLINT 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.
+#
+# FLINT 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 FLINT; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+#
+#============================================================================
+
+CPP = g++
+OBJ = lprels.o ModuloArith.o TonelliShanks.o F2matrix.o lanczos.o QS.o $(RES)
+LINKOBJ = lprels.o ModuloArith.o TonelliShanks.o F2matrix.o lanczos.o QS.o $(RES)
+LIBS = -L$(SAGE_LOCAL)/lib -lgmp
+CXXINCS = -I$(SAGE_LOCAL)/include
+BIN = QuadraticSieve QuadraticSieve.exe
+#CXXFLAGS = $(CXXINCS) -Wall -Wno-sign-compare -march=athlon-xp -fomit-frame-pointer -O2
+#CXXFLAGS2 = $(CXXINCS) -Wall -Wno-sign-compare -march=athlon-xp -fomit-frame-pointer -O3
+CXXFLAGS = $(CXXINCS) -Wall -Wno-sign-compare -march=opteron -fomit-frame-pointer -O2
+CXXFLAGS2 = $(CXXINCS) -Wall -Wno-sign-compare -march=opteron -fomit-frame-pointer -O3
+RM = rm -f
+
+.PHONY: all clean clean-custom
+
+all: QuadraticSieve
+
+clean: clean-custom
+ ${RM} $(OBJ) $(BIN)
+
+$(BIN): $(OBJ)
+ $(CPP) -ansi $(LINKOBJ) -o "QuadraticSieve" $(LIBS)
+
+ModuloArith.o: ModuloArith.cpp
+ $(CPP) -ansi -c ModuloArith.cpp -o ModuloArith.o $(CXXFLAGS)
+
+TonelliShanks.o: TonelliShanks.cpp
+ $(CPP) -ansi -c TonelliShanks.cpp -o TonelliShanks.o $(CXXFLAGS)
+
+F2matrix.o: F2matrix.cpp
+ $(CPP) -ansi -c F2matrix.cpp -o F2matrix.o $(CXXFLAGS2)
+
+lanczos.o: lanczos.c
+ $(CPP) -ansi -c lanczos.c -o lanczos.o $(CXXFLAGS2)
+
+lprels.o: lprels.c
+ $(CPP) -ansi -c lprels.c -o lprels.o $(CXXFLAGS2)
+
+QS.o: QS.cpp
+ $(CPP) -ansi -c QS.cpp -o QS.o $(CXXFLAGS)
diff --git a/makefile.osx64 b/makefile.osx64
new file mode 100644
index 0000000..eb22aab
--- /dev/null
+++ b/makefile.osx64
@@ -0,0 +1,60 @@
+#============================================================================
+# Copyright 2006 William Hart
+#
+# This file is part of FLINT.
+#
+# FLINT 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.
+#
+# FLINT 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 FLINT; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+#
+#============================================================================
+
+CPP = g++
+OBJ = lprels.o ModuloArith.o TonelliShanks.o F2matrix.o lanczos.o QS.o $(RES)
+LINKOBJ = lprels.o ModuloArith.o TonelliShanks.o F2matrix.o lanczos.o QS.o $(RES)
+LIBS = -L$(SAGE_LOCAL)/lib -lgmp
+CXXINCS = -I$(SAGE_LOCAL)/include
+BIN = QuadraticSieve QuadraticSieve.exe
+#CXXFLAGS = $(CXXINCS) -Wall -Wno-sign-compare -march=athlon-xp -fomit-frame-pointer -O2
+#CXXFLAGS2 = $(CXXINCS) -Wall -Wno-sign-compare -march=athlon-xp -fomit-frame-pointer -O3
+CXXFLAGS = $(CXXINCS) -Wall -Wno-sign-compare -fomit-frame-pointer -O2 -m64
+CXXFLAGS2 = $(CXXINCS) -Wall -Wno-sign-compare -fomit-frame-pointer -O3 -m64
+RM = rm -f
+
+.PHONY: all clean clean-custom
+
+all: QuadraticSieve
+
+clean: clean-custom
+ ${RM} $(OBJ) $(BIN)
+
+$(BIN): $(OBJ)
+ $(CPP) -ansi -m64 $(LINKOBJ) -o "QuadraticSieve" $(LIBS)
+
+ModuloArith.o: ModuloArith.cpp
+ $(CPP) -ansi -c ModuloArith.cpp -o ModuloArith.o $(CXXFLAGS)
+
+TonelliShanks.o: TonelliShanks.cpp
+ $(CPP) -ansi -c TonelliShanks.cpp -o TonelliShanks.o $(CXXFLAGS)
+
+F2matrix.o: F2matrix.cpp
+ $(CPP) -ansi -c F2matrix.cpp -o F2matrix.o $(CXXFLAGS2)
+
+lanczos.o: lanczos.c
+ $(CPP) -ansi -c lanczos.c -o lanczos.o $(CXXFLAGS2)
+
+lprels.o: lprels.c
+ $(CPP) -ansi -c lprels.c -o lprels.o $(CXXFLAGS2)
+
+QS.o: QS.cpp
+ $(CPP) -ansi -c QS.cpp -o QS.o $(CXXFLAGS)
diff --git a/makefile.sage b/makefile.sage
new file mode 100644
index 0000000..839f592
--- /dev/null
+++ b/makefile.sage
@@ -0,0 +1,60 @@
+#============================================================================
+# Copyright 2006 William Hart
+#
+# This file is part of FLINT.
+#
+# FLINT 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.
+#
+# FLINT 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 FLINT; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+#
+#============================================================================
+
+CPP = g++
+OBJ = lprels.o ModuloArith.o TonelliShanks.o F2matrix.o lanczos.o QS.o $(RES)
+LINKOBJ = lprels.o ModuloArith.o TonelliShanks.o F2matrix.o lanczos.o QS.o $(RES)
+LIBS = -L$(SAGE_LOCAL)/lib -lgmp
+CXXINCS = -I$(SAGE_LOCAL)/include
+BIN = QuadraticSieve QuadraticSieve.exe
+#CXXFLAGS = $(CXXINCS) -Wall -Wno-sign-compare -march=athlon-xp -fomit-frame-pointer -O2
+#CXXFLAGS2 = $(CXXINCS) -Wall -Wno-sign-compare -march=athlon-xp -fomit-frame-pointer -O3
+CXXFLAGS = $(CXXINCS) -Wall -Wno-sign-compare -fomit-frame-pointer -O2
+CXXFLAGS2 = $(CXXINCS) -Wall -Wno-sign-compare -fomit-frame-pointer -O3
+RM = rm -f
+
+.PHONY: all clean clean-custom
+
+all: QuadraticSieve
+
+clean: clean-custom
+ ${RM} $(OBJ) $(BIN)
+
+$(BIN): $(OBJ)
+ $(CPP) -ansi $(LINKOBJ) -o "QuadraticSieve" $(LIBS)
+
+ModuloArith.o: ModuloArith.cpp
+ $(CPP) -ansi -c ModuloArith.cpp -o ModuloArith.o $(CXXFLAGS)
+
+TonelliShanks.o: TonelliShanks.cpp
+ $(CPP) -ansi -c TonelliShanks.cpp -o TonelliShanks.o $(CXXFLAGS)
+
+F2matrix.o: F2matrix.cpp
+ $(CPP) -ansi -c F2matrix.cpp -o F2matrix.o $(CXXFLAGS2)
+
+lanczos.o: lanczos.c
+ $(CPP) -ansi -c lanczos.c -o lanczos.o $(CXXFLAGS2)
+
+lprels.o: lprels.c
+ $(CPP) -ansi -c lprels.c -o lprels.o $(CXXFLAGS2)
+
+QS.o: QS.cpp
+ $(CPP) -ansi -c QS.cpp -o QS.o $(CXXFLAGS)
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/flintqs.git
More information about the debian-science-commits
mailing list