Skip to content

Commit

Permalink
code simplification using newer C++ idioms e.g. auto
Browse files Browse the repository at this point in the history
  • Loading branch information
@JohnCremona
JohnCremona committed Mar 4, 2024
1 parent 68f94a2 commit bc5bb34
Show file tree
Hide file tree
Showing 33 changed files with 471 additions and 932 deletions.
62 changes: 22 additions & 40 deletions libsrc/arith.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@ long primeclass::number(long n)
}
if(!ok)
{
cout<<"Not enough primes in primeclass.number("<<n<<") !"<<endl;
cerr<<"Not enough primes in primeclass.number("<<n<<") !"<<endl;
}
return p_val;
}
Expand All @@ -135,7 +135,7 @@ vector<long> primeclass::getfirst (long n) /* returns list of first n primes */
}
if(!ok)
{
cout<<"Not enough primes in primeclass.getfirst("<<n<<") !"<<endl;
cerr<<"Not enough primes in primeclass.getfirst("<<n<<") !"<<endl;
}
return ans;
}
Expand Down Expand Up @@ -194,16 +194,15 @@ vector<long> posdivs(long a, const vector<long>& plist)
// cout << plist.size() << " primes: " <<endl; cout << plist << endl;
long j,k,p,e,nd = 1;
vector<long> dlist(1,1); // cout << "Divisor 0 = 1" << endl;
vector<long>::const_iterator pr = plist.begin();
while(pr!=plist.end())
for (const auto& p : plist)
{
p=*pr++; e = val(p,a);
e = val(p,a);
dlist.resize((e+1)*dlist.size());
for (j=0; j<e; j++)
for (k=0; k<nd; k++)
dlist[nd*(j+1)+k] = p*dlist[nd*j+k];
nd*=(e+1);
}
}
return dlist;
}

Expand All @@ -213,10 +212,9 @@ vector<long> alldivs(long a, const vector<long>& plist)
long j,k,p,e,nd = 2;
vector<long> dlist(1,1); // cout << "Divisor 0 = 1" << endl;
dlist.push_back(-1); // cout << "Divisor 1 = -1" << endl;
vector<long>::const_iterator pr = plist.begin();
while(pr!=plist.end())
for (const auto& p : plist)
{
p = *pr++; e = val(p,a);
e = val(p,a);
dlist.resize((e+1)*dlist.size());
for (j=0; j<e; j++)
for (k=0; k<nd; k++)
Expand All @@ -230,10 +228,9 @@ vector<long> sqdivs(long a, const vector<long>& plist)
{
long j,k,p,e,nd = 1;
vector<long> dlist(1,1);
vector<long>::const_iterator pr = plist.begin();
while(pr!=plist.end())
for (const auto& p : plist)
{
p = *pr++; e = val(p,a)/2;
e = val(p,a)/2;
dlist.resize((e+1)*dlist.size());
for (j=0; j<e; j++)
for (k=0; k<nd; k++)
Expand All @@ -247,10 +244,9 @@ vector<long> sqfreedivs(long a, const vector<long>& plist)
{
long j,k,p,e,nd = 1;
vector<long> dlist(1,1);
vector<long>::const_iterator pr = plist.begin();
while(pr!=plist.end())
for (const auto& p : plist)
{
p = *pr++; e = 1;
e = 1;
dlist.resize((e+1)*dlist.size());
for (j=0; j<e; j++)
for (k=0; k<nd; k++)
Expand Down Expand Up @@ -507,40 +503,26 @@ int is_squarefree(long n)
if(n%9==0) return 0;
if(n%25==0) return 0;
if(n%49==0) return 0;
vector<long>plist = pdivs(n);
for(unsigned int i=0; i<plist.size(); i++)
if(val(plist[i],n)>1) return 0;
auto plist = pdivs(n);
for( const auto p : plist)
if (val(p,n)>1)
return 0;
return 1;
}

int is_valid_conductor(long n)
{
long m=n, p, e;
if (n<11) return 0;
long m=n, e;
e=0; while(!(m&1)) {e++; m>>=1;} if(e>8) return 0;
e=0; while(!(m%3)) {e++; m/=3;} if(e>5) return 0;
vector<long>plist = pdivs(m);
for(unsigned int i=0; i<plist.size(); i++)
auto plist = pdivs(m);
for( const auto& p : plist)
{
p = plist[i];
e=0; while(!(m%p)) {e++; m/=p;} if(e>2) return 0;
if (val(p,m)>2)
return 0;
}
return 1;
}




// The following is no longer used

long old_kronecker(long d, long n)
{ long ans=0; long m=n, d4=d%4; if(d4<0)d4+=4;
if ((gcd(d,n)==1) && ((d4==0)||(d4==1)))
{ ans=1;
while (!(m%4)) m/=4;
if (!(m%2)) {m/=2; ans*=((d-1)%8 ? -1 : 1);}
ans *= legendre(d,m);
}
return ans;
return 1;
}

/* END OF FILE */
116 changes: 29 additions & 87 deletions libsrc/conic.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,19 +60,8 @@ int solve_conic(const bigint& a, const bigint& b, const bigint& c, const bigint&
bigint& x, bigint& y, bigint& z, int method)
{
vector<bigint> factorbase = pdivs(2*d);
// cout<<"factorbase(1) = "<<factorbase<<endl;
if(is_zero(b))
{
factorbase=vector_union(factorbase,pdivs(a));
factorbase=vector_union(factorbase,pdivs(c));
}
else
{
bigint disc = b*b-4*a*c;
factorbase=vector_union(factorbase,pdivs(a));
factorbase=vector_union(factorbase,pdivs(disc));
}
// cout<<"factorbase(2) = "<<factorbase<<endl;
factorbase=vector_union(factorbase,pdivs(a));
factorbase=vector_union(factorbase,pdivs((is_zero(b)? c : b*b-4*a*c)));
return solve_conic(a,b,c,d,factorbase,x,y,z,method);
}

Expand Down Expand Up @@ -112,7 +101,7 @@ int solve_conic(const bigint& a, const bigint& b, const bigint& c, const bigint&
vector<bigint> aplist, bplist, cplist, dplist;

int nondiag=!is_zero(b);
bb=a*d;
bb=a*d;
aa=-a*c; if(nondiag) aa=sqr(b)-4*aa;
aplist=factorbase;
bplist=factorbase;
Expand All @@ -123,7 +112,7 @@ int solve_conic(const bigint& a, const bigint& b, const bigint& c, const bigint&
if(solve_conic_diag(a1,aplist,b1,bplist,x,y,z,method))
{
conic_diag_reduce(a1,b1,x,y,z,verb);
x*=(a2*b2); y*=a2; z*=b2;
x*=(a2*b2); y*=a2; z*=b2;
if(nondiag) x-=b*z;
y *= a;
z *= a; if(nondiag) zz<<=1;
Expand All @@ -138,7 +127,7 @@ int solve_conic(const bigint& a, const bigint& b, const bigint& c, const bigint&
}
}
}

int solve_conic_diag_nontriv(const bigint& a, const vector<bigint>& aplist,
const bigint& b, const vector<bigint>& bplist,
bigint& x, bigint& y, bigint& z,
Expand Down Expand Up @@ -529,56 +518,22 @@ int testlocsol(const bigint& a, const vector<bigint>& alist,
// coprime and square-free, their prime factors being in alist etc.
{
int as=sign(a), bs=sign(b), cs=sign(c);
if((as==bs)&&(bs==cs))
{
// cout<<"testlocsol("<<a<<","<<b<<","<<c<<") returning 0 because of signs\n";
return 0;
}
bigint p, two; two=2;
bigint mab=-a*b;
vector<bigint>::const_iterator pr;
pr=clist.begin();
while(pr!=clist.end())
{
p=*pr++;
if(p==two) continue;
if(legendre(mab,p)!=1)
{
// cout<<"testlocsol fails legendre(mab,p) with "
// <<"(a,b,p)=("<<a<<","<<b<<","<<p<<")\n";
return 0;
}
}
bigint mbc=-b*c;
pr=alist.begin();
while(pr!=alist.end())
{
p=*pr++;
if(p==two) continue;
if(legendre(mbc,p)!=1)
{
// cout<<"testlocsol fails legendre(mbc,p) with "
// <<"(b,c,p)=("<<b<<","<<c<<","<<p<<")\n";
return 0;
}
}
bigint mca=-c*a;
pr=blist.begin();
while(pr!=blist.end())
{
p=*pr++;
if(p==two) continue;
if(legendre(mca,p)!=1)
{
// cout<<"testlocsol fails legendre(mca,p) with "
// <<"(c,a,p)=("<<c<<","<<a<<","<<p<<")\n";
return 0;
}
}
if((as==bs)&&(bs==cs))
return 0;
auto test = [](const vector<bigint>& plist, const bigint& m)
{
return !std::any_of(plist.begin(), plist.end(),
[m](const auto& p)
{
return (p>2) && (legendre(m,p)!=1);
}
);
};
return test(clist, -a*b) && test(alist, -b*c) && test(blist, -a*c);
return 1;
}
}

int testlocsol(const bigint& a, const vector<bigint>& alist,
int testlocsol(const bigint& a, const vector<bigint>& alist,
const bigint& b, const vector<bigint>& blist)
// tests if ax^2+by^2=z^2 is soluble, where a, b are
// square-free, their prime factors being in alist and blist.
Expand All @@ -589,20 +544,17 @@ int testlocsol(const bigint& a, const vector<bigint>& alist,
vector<bigint> a0list, b0list, clist;

long sa=sign(a), sb=sign(b);
if((sa<0)&&(sb<0))
if((sa<0)&&(sb<0))
{
// cout<<"testlocsol("<<a<<","<<b<<") returning 0 because of signs\n";
return 0; // nothing more to do as no real solution
}
if(sa<0) ::negate(a0);
if(sb<0) ::negate(b0);

vector<bigint>::const_iterator pr;
pr=alist.begin();
while(pr!=alist.end())

for (const auto& p : alist)
{
p=*pr++;
if(div(p,b))
if(div(p,b))
{
c*=p; clist.push_back(p);
}
Expand All @@ -611,26 +563,16 @@ int testlocsol(const bigint& a, const vector<bigint>& alist,
a0*=p; a0list.push_back(p);
}
}
pr=blist.begin();
while(pr!=blist.end())
for (const auto& p : blist)
{
p=*pr++;
if(!div(p,c)) {b0*=p; b0list.push_back(p);}
}

#if(0)
if((a!=-a0*c)||(b!=-b0*c)||(abs(c)!=gcd(a,b)))
{
cout<<"Error: (a,b)=("<<a<<","<<b<<") gives a0="<<a0<<", b0="<<b0<<", c="<<c<<endl;
if(!div(p,c))
{
b0*=p; b0list.push_back(p);
}
}
// cout<<"Calling testlocsol(a,b,c) with a,b,c="<<a0<<","<<b0<<","<<c<<"\n";
// cout<<"alist="<<a0list<<endl;
// cout<<"blist="<<b0list<<endl;
// cout<<"clist="<<clist<<endl;
#endif

// Now a=a0*c, b=b0*c, c=gcd(a,b), and a0list, b0list, clist hold their primes

return testlocsol(a0,a0list,b0,b0list,c,clist);

}
Expand Down
43 changes: 2 additions & 41 deletions libsrc/curve.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@

int valid_invariants(const bigint& c4, const bigint& c6)
{
bigint disc = c4*c4*c4; disc -= c6*c6;
bigint disc = c4*c4*c4 - c6*c6;
if (sign(disc)==0) return 0;
if (ndiv(1728,disc)) return 0; // need c4^3-c6^2=1728D, with D|=0
long x6= mod(c6,27);
Expand Down Expand Up @@ -83,10 +83,8 @@ void minimise_c4c6(const bigint& c4, const bigint& c6, const bigint& discr,
g = gcd( c62, newdiscr ); if(is_one(g)) return;
const vector<bigint>& p_list = pdivs(g);
// cout<<"g = "<<g<<endl;
vector<bigint>::const_iterator pr = p_list.begin();
while ( pr!=p_list.end() )
for (const auto& p : p_list)
{
p = *pr++;
d = (long)floor(val(p,g)/12.0);
// cout<<"With p="<<p<<", initial d="<<d<<endl;
if (p==2)
Expand Down Expand Up @@ -268,40 +266,3 @@ void Curve::tex_print(ostream &os) const


// end of file: curve.cc

//
// The following functions from OM's code are never used
//
// /* number of solutions to y^2+ay=b mod p */
// long quadroots(const bigint& a, const bigint& b, long p )
// {
// if (p == 2) {
// if (!odd(a)) return(1) ;
// if (!odd(b)) return(2) ;
// return(0) ;
// }
// else {
// long d = I2long((a*a+4*b)%p);
// if (d==0) return 1;
// if (legendre(d,p)==1) return 2;
// return 0;
// }
// }

// // for finding number of points mod 2 and 3
// long pointsmod(long p, const Curve& E)
// {
// bigint a1,a2,a3,a4,a6;
// E.getai(a1,a2,a3,a4,a6);
// if (p == 2)
// return( 1 + quadroots(a3, a6, 2) +
// quadroots(a1 + a3, 1 + a2 + a4 + a6, 2)) ;
// if (p ==3)
// return( 1 + quadroots(a3, a6, 3) +
// quadroots(a1 + a3, 1 + a2 + a4 + a6, 3)
// + quadroots(-a1 + a3, -1 + a2 - a4 + a6, 3)) ;
// // Now p>3
// long count=0;
// for (long x=0;x<p;x++) count+=quadroots(x*a1+a3,x*(x*(x+a2)+a4)+a6, p);
// return count;
// }
Loading

0 comments on commit bc5bb34

Please sign in to comment.