[flintqs] 03/04: Imported Upstream version 1.0
Dominique Ingoglia
godel108-guest at alioth.debian.org
Tue Aug 20 23:25:23 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 795e9be0618fbb4431c69d172b844d57c156a51a
Author: root <root at epsilon>
Date: Tue Aug 20 16:41:10 2013 -0600
Imported Upstream version 1.0
---
._F2matrix.cpp | Bin 229 -> 0 bytes
._F2matrix.h | Bin 229 -> 0 bytes
._ModuloArith.cpp | Bin 229 -> 0 bytes
._ModuloArith.h | Bin 229 -> 0 bytes
._QS.cpp | Bin 229 -> 0 bytes
._QS.cpp.orig | Bin 229 -> 0 bytes
._QStodo.txt | Bin 229 -> 0 bytes
._TonelliShanks.cpp | Bin 229 -> 0 bytes
._TonelliShanks.h | Bin 229 -> 0 bytes
._gpl.txt | Bin 229 -> 0 bytes
._lanczos.c | Bin 229 -> 0 bytes
._lanczos.h | Bin 229 -> 0 bytes
._lprels.c | Bin 229 -> 0 bytes
._lprels.h | Bin 229 -> 0 bytes
._makefile | Bin 229 -> 0 bytes
._makefile.opteron | Bin 229 -> 0 bytes
._makefile.sage | Bin 229 -> 0 bytes
.pc/.quilt_patches | 1 +
.pc/.quilt_series | 1 +
.pc/.version | 1 +
F2matrix.cpp | 197 ----
F2matrix.h | 46 -
ModuloArith.cpp | 69 --
ModuloArith.h | 42 -
QS.cpp | 1665 --------------------------------
QStodo.txt | 50 -
TonelliShanks.cpp | 147 ---
TonelliShanks.h | 45 -
changelog | 5 +
compat | 1 +
control | 15 +
copyright | 16 +
flintqs.dirs | 2 +
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 --
patches/flintqs-gcc-4.3-fix.patch.diff | 11 +
patches/lanczos.h.patch.diff | 17 +
patches/log.patch.diff | 12 +
patches/series | 4 +
patches/stdint.patch.diff | 11 +
rules | 6 +
source/format | 1 +
ssh | 27 +
ssh.pub | 1 +
51 files changed, 132 insertions(+), 4721 deletions(-)
diff --git a/._F2matrix.cpp b/._F2matrix.cpp
deleted file mode 100644
index 81618a0..0000000
Binary files a/._F2matrix.cpp and /dev/null differ
diff --git a/._F2matrix.h b/._F2matrix.h
deleted file mode 100644
index 7b51d3c..0000000
Binary files a/._F2matrix.h and /dev/null differ
diff --git a/._ModuloArith.cpp b/._ModuloArith.cpp
deleted file mode 100644
index cf91239..0000000
Binary files a/._ModuloArith.cpp and /dev/null differ
diff --git a/._ModuloArith.h b/._ModuloArith.h
deleted file mode 100644
index ad4f6b6..0000000
Binary files a/._ModuloArith.h and /dev/null differ
diff --git a/._QS.cpp b/._QS.cpp
deleted file mode 100644
index 775bcde..0000000
Binary files a/._QS.cpp and /dev/null differ
diff --git a/._QS.cpp.orig b/._QS.cpp.orig
deleted file mode 100644
index 4fff163..0000000
Binary files a/._QS.cpp.orig and /dev/null differ
diff --git a/._QStodo.txt b/._QStodo.txt
deleted file mode 100644
index 0418b86..0000000
Binary files a/._QStodo.txt and /dev/null differ
diff --git a/._TonelliShanks.cpp b/._TonelliShanks.cpp
deleted file mode 100644
index 208ee3e..0000000
Binary files a/._TonelliShanks.cpp and /dev/null differ
diff --git a/._TonelliShanks.h b/._TonelliShanks.h
deleted file mode 100644
index 61ddb95..0000000
Binary files a/._TonelliShanks.h and /dev/null differ
diff --git a/._gpl.txt b/._gpl.txt
deleted file mode 100644
index 0a90c9e..0000000
Binary files a/._gpl.txt and /dev/null differ
diff --git a/._lanczos.c b/._lanczos.c
deleted file mode 100644
index 42ca6c0..0000000
Binary files a/._lanczos.c and /dev/null differ
diff --git a/._lanczos.h b/._lanczos.h
deleted file mode 100644
index 0398ec8..0000000
Binary files a/._lanczos.h and /dev/null differ
diff --git a/._lprels.c b/._lprels.c
deleted file mode 100644
index cb4f9e7..0000000
Binary files a/._lprels.c and /dev/null differ
diff --git a/._lprels.h b/._lprels.h
deleted file mode 100644
index 234ec21..0000000
Binary files a/._lprels.h and /dev/null differ
diff --git a/._makefile b/._makefile
deleted file mode 100644
index b612b6b..0000000
Binary files a/._makefile and /dev/null differ
diff --git a/._makefile.opteron b/._makefile.opteron
deleted file mode 100644
index 1efda60..0000000
Binary files a/._makefile.opteron and /dev/null differ
diff --git a/._makefile.sage b/._makefile.sage
deleted file mode 100644
index f796829..0000000
Binary files a/._makefile.sage and /dev/null differ
diff --git a/.pc/.quilt_patches b/.pc/.quilt_patches
new file mode 100644
index 0000000..4baccb8
--- /dev/null
+++ b/.pc/.quilt_patches
@@ -0,0 +1 @@
+patches
diff --git a/.pc/.quilt_series b/.pc/.quilt_series
new file mode 100644
index 0000000..c206706
--- /dev/null
+++ b/.pc/.quilt_series
@@ -0,0 +1 @@
+series
diff --git a/.pc/.version b/.pc/.version
new file mode 100644
index 0000000..0cfbf08
--- /dev/null
+++ b/.pc/.version
@@ -0,0 +1 @@
+2
diff --git a/F2matrix.cpp b/F2matrix.cpp
deleted file mode 100644
index be30f08..0000000
--- a/F2matrix.cpp
+++ /dev/null
@@ -1,197 +0,0 @@
-/*============================================================================
- 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
deleted file mode 100644
index 7142470..0000000
--- a/F2matrix.h
+++ /dev/null
@@ -1,46 +0,0 @@
-/*============================================================================
- 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
deleted file mode 100644
index 95f17c5..0000000
--- a/ModuloArith.cpp
+++ /dev/null
@@ -1,69 +0,0 @@
-/*============================================================================
- 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
deleted file mode 100644
index 8e2e157..0000000
--- a/ModuloArith.h
+++ /dev/null
@@ -1,42 +0,0 @@
-/*============================================================================
- 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
deleted file mode 100644
index 4e02b8d..0000000
--- a/QS.cpp
+++ /dev/null
@@ -1,1665 +0,0 @@
-/*============================================================================
- 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
deleted file mode 100644
index 06929bb..0000000
--- a/QStodo.txt
+++ /dev/null
@@ -1,50 +0,0 @@
-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
deleted file mode 100644
index e845640..0000000
--- a/TonelliShanks.cpp
+++ /dev/null
@@ -1,147 +0,0 @@
-/*============================================================================
- 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
deleted file mode 100644
index 92048e2..0000000
--- a/TonelliShanks.h
+++ /dev/null
@@ -1,45 +0,0 @@
-/*============================================================================
- 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/changelog b/changelog
new file mode 100644
index 0000000..a62b3a2
--- /dev/null
+++ b/changelog
@@ -0,0 +1,5 @@
+flintqs (20070817-1.0) UNRELEASED; urgency=low
+
+ * Initial release. (Closes: #XXXXXX)
+
+ -- Dominique Ingoglia <Dominique.ingoglia at gmail.com> Fri, 16 Aug 2013 10:25:35 -0600
diff --git a/compat b/compat
new file mode 100644
index 0000000..301160a
--- /dev/null
+++ b/compat
@@ -0,0 +1 @@
+8
\ No newline at end of file
diff --git a/control b/control
new file mode 100644
index 0000000..215d22a
--- /dev/null
+++ b/control
@@ -0,0 +1,15 @@
+Source:flintqs
+Maintainer:Debian Science Maintainers <debian-science-maintainers at lists.alioth.debian.org>
+Uploaders:Dominique Ingoglia <dominique.ingoglia at gmail.com>
+Section:math
+Priority:optional
+Standards-Version: 3.9.3
+Build-depends:debhelper (>= 8)
+Vcs-Git: git://git.debian.org/git/debian-science/packages/flintqs.git
+Vcs-Browser: http://git.debian.org/?p=debian-science/packages/flintqs.git
+
+Package:flintqs
+Architecture:any
+Depends:${shlibs:Depends}, ${misc:Depends}
+Description:for flint
+ A library for number theory
diff --git a/copyright b/copyright
new file mode 100644
index 0000000..2fc4922
--- /dev/null
+++ b/copyright
@@ -0,0 +1,16 @@
+: Copyright © 2007-2008
+License: GPL-3+
+ 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 3 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 the Debian GNU/Linux distribution in /etc/share/common-licenses/GPL.
+ If not, see <http://www.gnu.org/licenses/>.
+
diff --git a/flintqs.dirs b/flintqs.dirs
new file mode 100644
index 0000000..f57433a
--- /dev/null
+++ b/flintqs.dirs
@@ -0,0 +1,2 @@
+usr/bin
+usr/share/man/man1
\ No newline at end of file
diff --git a/gpl.txt b/gpl.txt
deleted file mode 100644
index d511905..0000000
--- a/gpl.txt
+++ /dev/null
@@ -1,339 +0,0 @@
- 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
deleted file mode 100644
index 891551a..0000000
--- a/lanczos.c
+++ /dev/null
@@ -1,975 +0,0 @@
-/*============================================================================
- 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
deleted file mode 100644
index c07e6d1..0000000
--- a/lanczos.h
+++ /dev/null
@@ -1,95 +0,0 @@
-/*============================================================================
- 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
deleted file mode 100644
index 5b6613e..0000000
--- a/lprels.c
+++ /dev/null
@@ -1,756 +0,0 @@
-/*============================================================================
- 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
deleted file mode 100644
index 8c43b3e..0000000
--- a/lprels.h
+++ /dev/null
@@ -1,51 +0,0 @@
-/*============================================================================
- 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
deleted file mode 100644
index 0d92518..0000000
--- a/makefile
+++ /dev/null
@@ -1,64 +0,0 @@
-#============================================================================
-# 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
deleted file mode 100644
index 94ef78b..0000000
--- a/makefile.opteron
+++ /dev/null
@@ -1,60 +0,0 @@
-#============================================================================
-# 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
deleted file mode 100644
index eb22aab..0000000
--- a/makefile.osx64
+++ /dev/null
@@ -1,60 +0,0 @@
-#============================================================================
-# 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
deleted file mode 100644
index 839f592..0000000
--- a/makefile.sage
+++ /dev/null
@@ -1,60 +0,0 @@
-#============================================================================
-# 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)
diff --git a/patches/flintqs-gcc-4.3-fix.patch.diff b/patches/flintqs-gcc-4.3-fix.patch.diff
new file mode 100644
index 0000000..0a28fd1
--- /dev/null
+++ b/patches/flintqs-gcc-4.3-fix.patch.diff
@@ -0,0 +1,11 @@
+--- src/QS.cpp.orig 2008-04-14 14:16:38.000000000 -0700
++++ src/QS.cpp 2008-04-14 14:16:54.000000000 -0700
+@@ -1563,7 +1563,7 @@
+ Function: Factors a user specified number using a quadratic sieve
+
+ ===========================================================================*/
+-int main(int argc, unsigned char *argv[])
++int main(int argc, char *argv[])
+ {
+ unsigned long multiplier;
+
diff --git a/patches/lanczos.h.patch.diff b/patches/lanczos.h.patch.diff
new file mode 100644
index 0000000..1e40fec
--- /dev/null
+++ b/patches/lanczos.h.patch.diff
@@ -0,0 +1,17 @@
+diff -ur src/lanczos.h b/lanczos.h
+--- src/lanczos.h 2007-05-06 00:52:39.000000000 +0200
++++ b/lanczos.h 2012-05-26 06:17:05.784874818 +0200
+@@ -21,6 +21,13 @@
+ #ifndef LANCZOS_H
+ #define LANCZOS_H
+
++#ifdef __sun
++#define u_int32_t unsigned int
++#define u_int64_t unsigned long long
++#endif
++
++#include <sys/types.h> // needed on MacOS X 10.5 for uint type
++
+ #include <stdlib.h>
+
+ typedef struct {
diff --git a/patches/log.patch.diff b/patches/log.patch.diff
new file mode 100644
index 0000000..0dcf4c4
--- /dev/null
+++ b/patches/log.patch.diff
@@ -0,0 +1,12 @@
+diff -ru src/QS.cpp b/QS.cpp
+--- src/QS.cpp 2008-04-14 23:16:38.000000000 +0200
++++ b/QS.cpp 2012-05-26 06:30:15.823360044 +0200
+@@ -467,7 +467,7 @@
+ 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);
++ primeSizes[i]=(unsigned char)floor(log((double)factorBase[i])/log(2.0)-FUDGE+0.5);
+ }
+
+ return;
diff --git a/patches/series b/patches/series
new file mode 100644
index 0000000..3f70538
--- /dev/null
+++ b/patches/series
@@ -0,0 +1,4 @@
+flintqs-gcc-4.3-fix.patch.diff
+lanczos.h.patch.diff
+log.patch.diff
+stdint.patch.diff
diff --git a/patches/stdint.patch.diff b/patches/stdint.patch.diff
new file mode 100644
index 0000000..ec0ac0a
--- /dev/null
+++ b/patches/stdint.patch.diff
@@ -0,0 +1,11 @@
+diff -ru src/TonelliShanks.h b/TonelliShanks.h
+--- src/TonelliShanks.h 2007-05-06 00:52:39.000000000 +0200
++++ b/TonelliShanks.h 2012-06-21 12:18:28.353358700 +0200
+@@ -20,6 +20,7 @@
+ ============================================================================*/
+
+ #include <stdlib.h>
++#include <stdint.h>
+
+ // =====================================================================
+ // INTERFACE:
diff --git a/rules b/rules
new file mode 100755
index 0000000..c227ae6
--- /dev/null
+++ b/rules
@@ -0,0 +1,6 @@
+#!/usr/bin/make -f
+%:
+ dh $@
+
+override_dh_auto_install:
+ $(MAKE) DESTDIR=$$(pwd)/debian/flintqs prefix=/usr install
diff --git a/source/format b/source/format
new file mode 100644
index 0000000..46ebe02
--- /dev/null
+++ b/source/format
@@ -0,0 +1 @@
+3.0 (quilt)
\ No newline at end of file
diff --git a/ssh b/ssh
new file mode 100644
index 0000000..c3b526c
--- /dev/null
+++ b/ssh
@@ -0,0 +1,27 @@
+-----BEGIN RSA PRIVATE KEY-----
+MIIEowIBAAKCAQEArAivL4DCaHLKUrdt4u2vj/RIkeiB5ADUQc9t962tNLuqqLo3
+BsWOgmwluA0AzEyyHYIpdzPcOWl9OeA55oAAGX/2qJsoUwb1vZyzNtTm3FWsxPNh
+fvKfvvPB+nNS5xzas2srhmMRDt1/Gn0XNUUih6MB7CKPcI6r3QUcEgNfMa9TKvuu
+V0eCkqFdoA6fMLUhcYdw0UIT2CwAqHj22FRzdxQP5g0H9uplD6RaooUm2Sb77Guc
+ybnQPz9An/JZmt0o8a5s1i9DwL86HlCSA6vj0t0CG2Z1Q3wuL45O5mVRQfrq4cvs
+cYgPNSlQ9GC2UVTPzhb4yy6xZDbilSDZVeTr9QIDAQABAoIBACoQFYVv3hjbuExx
+PRT3OK3h9Lx4NQoiicNtjF26wVbba+bFYR7uvuF0v+Q4ibFqL0K3yJu0umvvNwcn
+pACP23Zgq1aeWUWztfIellMZyzikWhHt0DDR8e0mfI9YEzUfAPpNgd7h6hHQZnt7
+imkj9kVjvdyWtqu2tp7b2PkuieADrzs4nQfaS1Ojei8n+M88lul7MRvwflb5hv09
+jj93pnetXXH5whbEWeEktUqTuuQecFMtfBT183vTnEpu0hxDF8FltkX3IgzjBH3P
+V7xEO2d5KoG6YtUCajeIW/mx3H0t+6W7wV10OE8gOLUGONYiA5xuTg6D9domCZZm
+sIjZPCECgYEA5IAifFuZkEzzupztW1RvsvYtKpfBHkemPMbQJ2scllJgaWZg2Unt
+yGvrmqwqnwoNQDgewVskFrw24FJH++xpEs1FYqB0hOpbBk0qSk4afNq0GhKPOQIq
+ixiHvLtA9cea1xcLct2Fb7gL6dZ91n1lUmS0H5wGeTmTywiQtU03zwsCgYEAwLzf
+Mr4WL+inJ0Bxh3vRSI25McoHgUR6jPGQFmhlJBZ0LPaswaWZS2S8GQ5dgrL2Iatj
+fbWdIWQ+0hWCCKPDIrk4OO3QmbV2tmyAzr608C0gkTTbQJw8DMGby4WJe5K/SAGX
+VcXEKRW1oAkSEdvxH6AOtfW9LJDzClBcZTk3EP8CgYEAyqeJ7lkfHQfisgMzz+hX
+GJWVAU2ODVjmasi5G/y3Yeq1b0VJZ+1VYoe0cX14X4z+q5IaVMqMe016LgFLrnbB
+ydccTpiYPrnK+Q+/Dh+vBkTBrs3/EESHjs22tQAuYM0i2tipYrps+eR1THLbMDwO
+fMCrr80lQKZ8GXoDPYi6knkCgYBdycrG827yg0ELvbVBG4RczPJIgyohwkPsYAQg
+k05cQDzqQGMSnFW7NVq+ypnAZvuUqMTyQDUlMZXMP0EWmTH0rLLqKPdwRLhuzt/j
+OzPrB9qoLlNe3mfuQSxh3ipnoqJIFNYim+j3oSPPq3pKjH+KRyXBb8JNdH+ADljX
+vP7J2wKBgGJRP4YBBfE3gjluydMsVYSw5+tH8uY+NPII8o2y+5gNjUkrfUT6CO47
+yU3Un2K8bb8pz/LtPTVu4zDqKqhI10plohMLsVW2zRJKNd1Y/uWVwU4ELYoVPyee
+sCaLnYfJlnbgSMN6TExctJ3ezuZsLTwuc1o5haKc3KjLa4+9daGx
+-----END RSA PRIVATE KEY-----
diff --git a/ssh.pub b/ssh.pub
new file mode 100644
index 0000000..4235187
--- /dev/null
+++ b/ssh.pub
@@ -0,0 +1 @@
+ssh-rsa AAAAB3NzaC1yc2EAAAADAQABAAABAQCsCK8vgMJocspSt23i7a+P9EiR6IHkANRBz233ra00u6qoujcGxY6CbCW4DQDMTLIdgil3M9w5aX054DnmgAAZf/aomyhTBvW9nLM21ObcVazE82F+8p++88H6c1LnHNqzayuGYxEO3X8afRc1RSKHowHsIo9wjqvdBRwSA18xr1Mq+65XR4KSoV2gDp8wtSFxh3DRQhPYLACoePbYVHN3FA/mDQf26mUPpFqihSbZJvvsa5zJudA/P0Cf8lma3SjxrmzWL0PAvzoeUJIDq+PS3QIbZnVDfC4vjk7mZVFB+urhy+xxiA81KVD0YLZRVM/OFvjLLrFkNuKVINlV5Ov1 godel at epsilon
--
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