Skip to content

Commit

Permalink
more of the same, WIP
Browse files Browse the repository at this point in the history
  • Loading branch information
@JohnCremona
JohnCremona committed Mar 6, 2024
1 parent 1d5eb75 commit 4c404b8
Show file tree
Hide file tree
Showing 10 changed files with 156 additions and 178 deletions.
13 changes: 2 additions & 11 deletions libsrc/bitspace.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,9 @@
//
//////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <eclib/interface.h>
#include <eclib/bitspace.h>



bitspace::bitspace(long d)
{
if(d<0)
Expand All @@ -43,18 +40,12 @@ bitspace::bitspace(long d)
d=NTL_BITS_PER_LONG;
}
maxdim=d;
pivs = new long[maxdim];
gens = new unsigned long[maxdim];
pivs.resize(maxdim);
gens.resize(maxdim);
dim=0;
bitmask=0;
}

bitspace::~bitspace()
{
delete[]pivs;
delete[]gens;
}

long bitspace::reduce(unsigned long& v, long start) const
{
long j;
Expand Down
9 changes: 5 additions & 4 deletions libsrc/eclib/bitspace.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,16 +24,17 @@
#if !defined(_ECLIB_BITSPACE_H)
#define _ECLIB_BITSPACE_H 1 //flags that this file has been included

#include <eclib/templates.h>

class bitspace {
private:
long maxdim;
long dim;
long * pivs; // holds the position of the ith pivot
unsigned long * gens; // holds the ith basis element
unsigned long bitmask; // holds the bits of the pivs
vector<long> pivs; // holds the position of the ith pivot
vector<unsigned long> gens; // holds the ith basis element
unsigned long bitmask; // holds the bits of the pivs
public:
bitspace(long d);
~bitspace();
unsigned long getbitmask() {return bitmask;}
long reduce(unsigned long& v, long start=0) const;
// reduces v mod this, returns minimal i such that the reduced v has
Expand Down
4 changes: 2 additions & 2 deletions libsrc/eclib/marith.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,8 +173,8 @@ int modsqrt(const bigint& a, const vector<bigint>& bplist, bigint& x);
// root-finding functions:

// find the number of roots of X^3 + bX^2 + cX + d = 0 (mod p)
// roots are put in r which should be allocated of size 3
int nrootscubic(long b, long c, long d, long p, long* roots);
// and assign roots to a list of these
int nrootscubic(long b, long c, long d, long p, vector<long>& roots);

void ratapprox(bigfloat x, bigint& a, bigint& b, const bigint& maxd=BIGINT(0));
void ratapprox(bigfloat x, long& a, long& b, long maxd=0);
Expand Down
11 changes: 6 additions & 5 deletions libsrc/eclib/mrank1.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,13 +50,14 @@ class rank1 : public rank12 {
//
//
// Sieving stuff:
int ipivot, pivflag;
long * auxs; long ** phimod; int * aux_flags; int * aux_types;
int ipivot, pivflag;
vector<long> auxs, amod, astepmod, hmod, hstepmod, hscalemod;
vector<vector<long>> phimod;
vector<int> aux_flags, aux_types;

int **squares;
int **squares;
int ***flags;
int **flaga; int *flagah;
long *amod, *astepmod, *ascalemod, *hmod, *hstepmod, *hscalemod;
int **flaga;
long ah_count, ah_sieve_0, ah_sieve_1, ah_sieve_2;
long ah_rfail, ah_dfail, ah_efail, ah_extra2fail, ah_pass;
void aux_init(); // define auxiliary moduli and squares
Expand Down
2 changes: 1 addition & 1 deletion libsrc/eclib/sifter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ class sifter {
int ** eps_mat;
int * pivcols;
long * auxs; long * all_p; int * nroots;
long ** thetamod; int**squares;
vector<vector<long>> thetamod; int**squares;
public:
sifter(Curvedata* EE, int na, int verb=0);
~sifter();
Expand Down
44 changes: 20 additions & 24 deletions libsrc/marith.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1115,32 +1115,28 @@ int modrat(const bigint& n, const bigint& m, const bigint& lim,
return 0;
}

// find the number of roots of X^3 + bX^2 + cX + d = 0 (mod p)
// roots are put in r which should be allocated of size 3
int nrootscubic(long b, long c, long d, long p, long* roots)
// Find the number of roots of X^3 + bX^2 + cX + d = 0 (mod p) and
// assign roots to a list of these. A stupid search is good enough
// since we only use this for very small p! Also it is tacitly assumed
// that the roots are distinct.
int nrootscubic(long b, long c, long d, long p, vector<long>& roots)
{
long r, nr=0;
int found = 0;
for (r = 0; (r<p)&&!found ; r++)
int nr=0;
long r=0;
roots.clear();
for (r = 0; r<p; r++)
{
found = (((((r+b)*r+c)*r+d)%p)==0 );
}
if (!found) return 0;

r--; // because it was incremented one extra time in the loop!
roots[nr++]=r;
long e = b + r;
long f = c + e*r;
long e0 = (-e*((p+1)/2))%p;
long dd = posmod(e0*e0-f,p);
if(legendre(dd,p)==1)
{ // stupid search is good enough since we only use this for very small p!
for (r = 1; r<p ; r++)
{
if((r*r-dd)%p==0) break;
}
roots[nr++] = (e0+r)%p;
roots[nr++] = (e0-r)%p;
if (((((r+b)*r+c)*r+d)%p)==0)
{
roots.push_back(r);
nr++;
if (nr==2) // find 3rd root the easy way
{
roots.push_back(posmod(-b-roots[0]-roots[1],p));
std::sort(roots.begin(), roots.end());
return 3;
}
}
}
return nr;
}
Expand Down
98 changes: 45 additions & 53 deletions libsrc/mrank1.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -341,7 +341,7 @@ void rank1::addquartic(const bigint& a, const bigint& b, const bigint& c,
{
if(!qlistbflag[i]) continue; // i'th already redundant
// compute (a,h) of i'th quartic:
long aa = posmod(qlist[i].geta(),auxpiv);
long aa = posmod(qlist[i].geta(),auxpiv);
long H = posmod(qlist[i].getH(),auxpiv);
if(extra2) if(i>=nfl) H=posmod(hscale*H,auxpiv);

Expand Down Expand Up @@ -913,27 +913,29 @@ void rank1::gettype(int t) // new hybrid version 13/2/96
}
#endif

long iaux;
long *amodi, *hmodi, *auxi, *hstepmodi, *hscalemodi, *astepmodi;

long iaux=num_aux;
int ***flagsi; int **flagai;
iaux=num_aux; amodi=amod; auxi=auxs; astepmodi=astepmod;
while(iaux--)
auto auxi=auxs.begin();
auto amodi=amod.begin();
auto hmodi=hmod.begin();
auto astepmodi=astepmod.begin();
auto hstepmodi=hstepmod.begin();
auto hscalemodi=hscalemod.begin();
while(iaux--)
{
*amodi++ = posmod(a, *auxi);
*amodi++ = posmod(a, *auxi);
*astepmodi++ = posmod(astep,*auxi);
auxi++;
}

while(a!=lasta)
{
a+=astep;
for(iaux=0, amodi=amod, auxi=auxs, flagai=flaga, flagsi=flags,
astepmodi=astepmod;
iaux<num_aux;
iaux++, amodi++, auxi++, flagai++, flagsi++, astepmodi++)
for(iaux=0, amodi=amod.begin(), auxi=auxs.begin(), flagai=flaga, flagsi=flags, astepmodi=astepmod.begin();
iaux<num_aux;
iaux++, amodi++, auxi++, flagai++, flagsi++, astepmodi++)
{
(*amodi)+=(*astepmodi);
(*amodi)+=(*astepmodi);
if((*amodi)>=(*auxi)) (*amodi)-=(*auxi);
*flagai = (*flagsi)[*amodi];
}
Expand Down Expand Up @@ -972,10 +974,13 @@ void rank1::gettype(int t) // new hybrid version 13/2/96
bigfloat xa=to_bigfloat(a);
bigfloat xa4 = 4*xa, xa8=8*xa;
bigfloat oneover4a = 1/xa4;

// hstep does not depend on b so can be set here:
long hstep = 8*a*cstep;
iaux=num_aux; hstepmodi=hstepmod; auxi=auxs; hscalemodi=hscalemod;
iaux=num_aux;
auxi=auxs.begin();
hstepmodi=hstepmod.begin();
hscalemodi=hscalemod.begin();
if(extra2)
while(iaux--)
{
Expand Down Expand Up @@ -1124,8 +1129,9 @@ void rank1::gettype(int t) // new hybrid version 13/2/96

c=cmin-cstep; // So first value used is cmin

iaux=num_aux; hmodi=hmod; auxi=auxs; hscalemodi=hscalemod;
while(iaux--)
iaux=num_aux;
hmodi=hmod.begin(); auxi=auxs.begin(); hscalemodi=hscalemod.begin();
while(iaux--)
{
long aux=(*auxi);
long cmod = c%aux, bmod = b%aux;
Expand All @@ -1142,9 +1148,9 @@ void rank1::gettype(int t) // new hybrid version 13/2/96
{
ah_count++;
c+=cstep;
for(iaux=0, hmodi=hmod, hstepmodi=hstepmod, auxi=auxs;
iaux<num_aux;
iaux++, hmodi++, hstepmodi++, auxi++)
for(iaux=0, hmodi=hmod.begin(), hstepmodi=hstepmod.begin(), auxi=auxs.begin();
iaux<num_aux;
iaux++, hmodi++, hstepmodi++, auxi++)
{
(*hmodi) += (*hstepmodi);
if((*hmodi)>=(*auxi)) (*hmodi)-=(*auxi);
Expand All @@ -1155,7 +1161,7 @@ void rank1::gettype(int t) // new hybrid version 13/2/96

ipivot=-1;

for(iaux=0, hmodi=hmod, flagai=flaga;
for(iaux=0, hmodi=hmod.begin(), flagai=flaga;
flagok&&(iaux<num_aux);
iaux++, hmodi++, flagai++)
{
Expand Down Expand Up @@ -1272,7 +1278,7 @@ void rank1::gettype(int t) // new hybrid version 13/2/96
if (xxe>abceps) {ah_efail++; continue;}
d = Iround(xd);
e = Iround(xe);

#endif
#ifdef DEBUG_AH
cout << ":\n [" << a<<","<<sb << "," << c
Expand Down Expand Up @@ -1753,22 +1759,22 @@ void rank1::aux_init() // define auxiliary moduli and squares
{
long i, j, a;

auxs = new long[num_aux];
phimod = new long*[num_aux];
aux_flags = new int[num_aux];
aux_types = new int[num_aux];
auxs.resize(num_aux);
aux_flags.resize(num_aux);
aux_types.resize(num_aux);
phimod.resize(num_aux, vector<long>(3));
squares = new int*[num_aux];
flags = new int**[num_aux];
flaga = new int*[num_aux];
amod = new long[num_aux];
hmod = new long[num_aux];
hstepmod = new long[num_aux];
astepmod = new long[num_aux];
hscalemod = new long[num_aux];
amod.resize(num_aux);
hmod.resize(num_aux);
hstepmod.resize(num_aux);
astepmod.resize(num_aux);
hscalemod.resize(num_aux);

auxs[0]=9; // treated specially
aux_flags[0]=1; aux_types[0]=0;
for(i=0; i<num_aux; i++) phimod[i] = new long[3];
aux_flags[0]=1;
aux_types[0]=0;
i=1;

// the rest of the auxs must be chosen carefully: if possible they should
Expand Down Expand Up @@ -1862,19 +1868,14 @@ void rank1::flag_init() // set up flag array
int thisflag;
int ***flagsi=flags;
int **squaresi=squares;
long * a4phi= new long[3];
long * eps = new long[3];

vector<long> a4phi(3);
vector<long> eps(3);
#ifdef COUNT_CODES
long * code_count = new long[5]; long icc;
vector<long> code_count(5, 0);
#endif

for(long i=0; i<num_aux; i++, squaresi++, flagsi++)
for(long i=0; i<num_aux; i++, squaresi++, flagsi++)
{
#ifdef COUNT_CODES
for(icc=0; icc<5; icc++) code_count[icc]=0;
#endif

int case1 = (aux_types[i]==1); // phi cubic has 1 root mod p
int case2 = !case1; // phi cubic has 3 roots mod p
long a, h;
Expand Down Expand Up @@ -1961,7 +1962,7 @@ void rank1::flag_init() // set up flag array

//(+,+,+) maps to flag 15=8+4+2+1 (+,-,-) maps to flag 5= 4 + 1
//(-,+,-) maps to flag 3= 2+1 (-,-,+) maps to flag 1= 1

if(eps[0]==1) thisflag=(eps[1]==1? 15: 5);
else thisflag=(eps[1]==1? 3: 1);
}
Expand All @@ -1982,17 +1983,14 @@ void rank1::flag_init() // set up flag array
{
cout << "Code count for p = " << aux << ":\n";
cout << 0 << "\t"<< 15 << "\t"<< 5 << "\t"<< 3 << "\t"<< 1 << "\n";
for(icc=0; icc<5; icc++) cout<<code_count[icc]<<"\t";
for (const auto& cc : code_count)
cout<<cc<<"\t";
cout<<endl;
double ratio = ((double)(code_count[0]))/(aux*aux);
cout<<"Percentage of (a,H) pairs failing sieve = "<<100*ratio<<endl;
}
#endif
}
delete [] a4phi; delete [] eps;
#ifdef COUNT_CODES
delete [] code_count;
#endif

if((verbose>1)&&(num_aux>0))
cout<<"finished flag_init()"<<endl;
Expand All @@ -2010,14 +2008,8 @@ void rank1::clear_sieve() // free memory related to sieve;
delete[] flags[i][a];
}
delete[] flags[i];
delete[] phimod[i];
}
num_aux=0;
delete[] auxs;
delete[] phimod;
delete[] squares;
delete[] aux_flags; delete[] aux_types;
delete[] flags; delete[] flaga;
delete[] amod; delete[] hstepmod; delete[] hscalemod;
delete[] hmod; delete[] astepmod;
}
5 changes: 1 addition & 4 deletions libsrc/sifter.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,11 +42,10 @@ sifter::sifter(Curvedata* EE, int na, int verb)

auxs = new long[num_aux];
nroots = new int[num_aux];
thetamod = new long*[num_aux];
thetamod.resize(num_aux, vector<long>(3));
squares = new int*[num_aux];
all_p = new long[2*num_aux];

for(i=0; i<num_aux; i++) thetamod[i] = new long[3];
iaux=0;
max_dim_im=0;

Expand Down Expand Up @@ -130,8 +129,6 @@ sifter::~sifter()
long i;
for(i=0; i<max_dim_im; i++) delete[]eps_mat[i];
delete[]eps_mat;
for(i=0; i<num_aux; i++) delete[]thetamod[i];
delete[]thetamod;
for(i=0; i<num_aux; i++) delete[]squares[i];
delete[]squares;
delete[]auxs;
Expand Down
Loading

0 comments on commit 4c404b8

Please sign in to comment.