diff -Nru eclib-20220621/configure.ac eclib-20221012/configure.ac --- eclib-20220621/configure.ac 2022-06-21 14:25:26.000000000 +0000 +++ eclib-20221012/configure.ac 2022-10-12 11:04:09.000000000 +0000 @@ -7,7 +7,7 @@ # The version is the concatenation of 'v' and the year, month, date: vyyyymmdd # To access this as a string or integer triple, see libsrc/eclib/version.h -AC_INIT([eclib], [20220621], [john.cremona@gmail.com]) +AC_INIT([eclib], [20221012], [john.cremona@gmail.com]) AM_INIT_AUTOMAKE([-Wall]) AC_MSG_NOTICE([Configuring eclib...]) @@ -40,9 +40,9 @@ # # NB The suffix of the library name (libec.so here) is (c-a).a.r -LT_CURRENT=10 +LT_CURRENT=11 LT_REVISION=0 -LT_AGE=0 +LT_AGE=1 AC_SUBST(LT_CURRENT) AC_SUBST(LT_REVISION) AC_SUBST(LT_AGE) diff -Nru eclib-20220621/debian/changelog eclib-20221012/debian/changelog --- eclib-20220621/debian/changelog 2022-06-22 07:03:59.000000000 +0000 +++ eclib-20221012/debian/changelog 2022-10-18 05:06:17.000000000 +0000 @@ -1,3 +1,15 @@ +eclib (20221012-1) unstable; urgency=medium + + * New upstream release. + + -- Julien Puydt Tue, 18 Oct 2022 07:06:17 +0200 + +eclib (20220621-2) unstable; urgency=medium + + * Fix d/watch. + + -- Julien Puydt Sun, 02 Oct 2022 21:48:28 +0200 + eclib (20220621-1) unstable; urgency=medium * New upstream release. diff -Nru eclib-20220621/debian/watch eclib-20221012/debian/watch --- eclib-20220621/debian/watch 2022-06-22 07:03:59.000000000 +0000 +++ eclib-20221012/debian/watch 2022-10-18 05:06:17.000000000 +0000 @@ -1,2 +1,2 @@ version=4 -https://github.com/JohnCremona/eclib/releases .*/v?(\d\S*)\.tar\.gz +https://github.com/JohnCremona/eclib/tags .*/v?(\d\S*)\.tar\.gz diff -Nru eclib-20220621/doc/sage_packaging.txt eclib-20221012/doc/sage_packaging.txt --- eclib-20220621/doc/sage_packaging.txt 1970-01-01 00:00:00.000000000 +0000 +++ eclib-20221012/doc/sage_packaging.txt 2022-10-12 11:04:09.000000000 +0000 @@ -0,0 +1,35 @@ +Notes on how to package a new version of eclib for Sage + +After making a new eclib release VERSION: + +1. make dist-bzip2 + This creates the file eclib-.tar.bz2 + Copy it to ~/sage/upstream + +2. In ~/sage, checkout the development branch, update it from trac and + run "make". + +3. Make a new sage branch, say + git checkout -b eclib- + +4. In ~/sage/build/pkgs/eclib, update the version number in + package-version.txt and spkg-configure.m4 + +5. Run (using this sage) + sage --package fix-checksum eclib + +6. Build this new eclib in sage, with checks: + sage -i -c eclib + + [Doing 'make eclib' is also OK, bu does not run eclib's checks.] + +7. Test sage to see if any doctests need changing. Start with: + sage -t --long src/sage/libs/eclib + sage -t --long src/sage/schemes/elliptic_curves + +8. Make a github release and upload eclib-.tar.bz2 to it. + +9. Commit the changes to sage, upload the branch to the trac ticket + and call for review. + + diff -Nru eclib-20220621/libsrc/cubic.cc eclib-20221012/libsrc/cubic.cc --- eclib-20220621/libsrc/cubic.cc 2022-06-21 14:25:26.000000000 +0000 +++ eclib-20221012/libsrc/cubic.cc 2022-10-12 11:04:09.000000000 +0000 @@ -28,31 +28,182 @@ #include #include -void cubic::init() -{ - coeffs = new bigint[4]; -} +bigint zero(0), one(1); -cubic::~cubic() +// comparison operator for sorting +int operator<(const cubic& F1, const cubic& F2) { - delete [] coeffs; + return std::lexicographical_compare(F1.coeffs.begin(), F1.coeffs.end(), F2.coeffs.begin(), F2.coeffs.end()); } -void cubic::transform(const unimod& m) +vector transform_helper(const bigint& a, const bigint& b, const bigint& c, const bigint& d, + const unimod& m) { bigint m112=sqr(m(1,1)); bigint m113=m112*m(1,1); bigint m212=sqr(m(2,1)); bigint m213=m212*m(2,1); bigint m222=sqr(m(2,2)); bigint m223=m222*m(2,2); bigint m122=sqr(m(1,2)); bigint m123=m122*m(1,2); - bigint A = m113*a() + m(2,1)*m112*b() + m212*m(1,1)*c() + m213*d(); - bigint B = 3*m(1,2)*m112*a() + (m(2,2)*m112 + 2*m(2,1)*m(1,2)*m(1,1))*b() - + (2*m(2,2)*m(2,1)*m(1,1) + m212*m(1,2))*c() + 3*m(2,2)*m212*d(); - bigint C = 3*m122*m(1,1)*a() + (2*m(2,2)*m(1,2)*m(1,1) + m(2,1)*m122)*b() - + (m222*m(1,1) + 2*m(2,2)*m(2,1)*m(1,2))*c() + 3*m222*m(2,1)*d(); - bigint D = m123*a() + m(2,2)*m122*b() + m222*m(1,2)*c() + m223*d(); - set(A,B,C,D); - } + bigint A = m113*a + m(2,1)*m112*b + m212*m(1,1)*c + m213*d; + bigint B = 3*m(1,2)*m112*a + (m(2,2)*m112 + 2*m(2,1)*m(1,2)*m(1,1))*b + + (2*m(2,2)*m(2,1)*m(1,1) + m212*m(1,2))*c + 3*m(2,2)*m212*d; + bigint C = 3*m122*m(1,1)*a + (2*m(2,2)*m(1,2)*m(1,1) + m(2,1)*m122)*b + + (m222*m(1,1) + 2*m(2,2)*m(2,1)*m(1,2))*c + 3*m222*m(2,1)*d; + bigint D = m123*a + m(2,2)*m122*b + m222*m(1,2)*c + m223*d; + + return {A,B,C,D}; +} + +vector transform_helper(const vector& abcd, const unimod& m) +{ + return transform_helper(abcd[0],abcd[1],abcd[2],abcd[3],m); +} + +void cubic::transform(const unimod& m) +{ + coeffs = transform_helper(coeffs, m); +} + +cubic transform(const cubic& F, const unimod& m) +{ + return cubic(transform_helper(F.coeffs, m)); +} + +void cubic::sl2_reduce(unimod& m) +{ + if (disc()<0) + jc_reduce(m); + else + hess_reduce(m); +} + +//#define DEBUG_NORMALISE + +// - for an sl2-reduced cubic, normalise w.r.t. <-I> (default) or or +// - updates m by multiplying by normalising transformation +void cubic::normalise(unimod& m) +{ + int nautos = 2; + +#ifdef DEBUG_NORMALISE + cout<<"Normalising "<<(*this)<coeffs; + autpower = *(transforms.begin() + (biggest-Flist.begin())); +#ifdef DEBUG_NORMALISE + cout<<"Largest one after sorting is "<<(*this)<<", the transform by "<& Glist) const +{ + for (auto Gi=Glist.begin(); Gi!=Glist.end(); ++Gi) + if (sl2_equivalent(*Gi)) + return 1; + return 0; +} + +int cubic::gl2_equivalent_in_list(const vector& Glist) const +{ + for (auto Gi=Glist.begin(); Gi!=Glist.end(); ++Gi) + if (gl2_equivalent(*Gi)) + return 1; + return 0; +} + +// affine roots of F mod q. +// NB rootsmod requires a non-constant polynomial + +// affine roots of F mod q, assuming leading coefficient a() is nonzero: +vector cubic::roots_mod(const bigint& q) const +{ + bigint aq(a()%q), bq(b()%q), cq(c()%q), dq(d()%q); + if (is_zero(aq) && is_zero(bq) && is_zero(cq)) + return {}; + return rootsmod({dq,cq,bq,aq}, q); +} + +// Return 1 iff F has a projective root mod q: +int cubic::has_roots_mod(const bigint& q) const +{ + return div(q,a()) || (roots_mod(q).size() > 0); +} void cubic::x_shift(const bigint& e, unimod& m) { @@ -96,61 +247,63 @@ bigint cubic::j_c1() const { - bigint b2=sqr(b()); - bigint b3=b()*b2; - bigint b4=b()*b3; - bigint b5=b()*b4; - bigint b6=b()*b5; - bigint a2=sqr(a()); - bigint a3=a()*a2; - bigint a4=a()*a3; - bigint c2=sqr(c()); - bigint c3=c()*c2; - bigint c4=c()*c3; - bigint c5=c()*c4; - bigint c6=c()*c5; - bigint d2=sqr(d()); - bigint d3=d()*d2; - bigint d4=d()*d3; - bigint ac=a()*c(), bd=b()*d(); - return - 108*b3*a2*d() - 3*b4*c2 + 54*a2*c4 + 18*b5*d() + 243*a2*d2*b2 - - 54*b3*ac*d() - 162*bd*c2*a2 - 54*a3*c3 + 486*a3*bd*c() + 3*c4*b2 - - 18*c5*a() + 54*c3*a()*bd - 243*d2*a2*c2 + 162*d2*ac*b2 + 2*c6 - + bigint a = coeffs[0], b=coeffs[1], c=coeffs[2], d=coeffs[3]; + bigint b2=sqr(b); + bigint b3=b*b2; + bigint b4=b*b3; + bigint b5=b*b4; + bigint b6=b*b5; + bigint a2=sqr(a); + bigint a3=a*a2; + bigint a4=a*a3; + bigint c2=sqr(c); + bigint c3=c*c2; + bigint c4=c*c3; + bigint c5=c*c4; + bigint c6=c*c5; + bigint d2=sqr(d); + bigint d3=d*d2; + bigint d4=d*d3; + bigint ac=a*c, bd=b*d; + return - 108*b3*a2*d - 3*b4*c2 + 54*a2*c4 + 18*b5*d + 243*a2*d2*b2 - + 54*b3*ac*d - 162*bd*c2*a2 - 54*a3*c3 + 486*a3*bd*c + 3*c4*b2 - + 18*c5*a + 54*c3*a*bd - 243*d2*a2*c2 + 162*d2*ac*b2 + 2*c6 - 729*a4*d2 - 2*b6 + 18*b4*ac - 27*a2*b2*c2 + 729*d4*a2 + 54*b3*d3 + - 108*c3*d2*a() - 18*c4*bd + 27*d2*c2*b2 - 486*d3*ac*b() - 54*d2*b4; + 108*c3*d2*a - 18*c4*bd + 27*d2*c2*b2 - 486*d3*ac*b - 54*d2*b4; } -// MINUS the quantity called C_2 in the paper, = -Norm(h0-h1) and should be -// NON-POSITIVE for a reduced form: +// The quantity called C_2 in the paper, = Norm(h0-h1) and should be +// NON-NEGATIVE for a reduced form: bigint cubic::j_c2() const { - bigint b2=sqr(b()); - bigint b3=b()*b2; - bigint b4=b()*b3; - bigint b5=b()*b4; - bigint b6=b()*b5; - bigint a2=sqr(a()); - bigint a3=a()*a2; - bigint a4=a()*a3; - bigint c2=sqr(c()); - bigint c3=c()*c2; - bigint c4=c()*c3; - bigint c5=c()*c4; - bigint c6=c()*c5; - bigint d2=sqr(d()); - bigint d3=d()*d2; - bigint d4=d()*d3; - bigint ac=a()*c(), bd=b()*d(); - - return - 108*b3*a2*d() + 12*b4*c2 - 216*a2*c4 - 72*b5*d() - 486*a3*c2*d() - + 270*a2*c3*b() - 90*b3*c2*a() - 972*a2*d2*b2 + 216*b3*ac*d() + - 648*bd*c2*a2 - 54*a3*c3 + 486*a3*bd*c() - 16*c3*b3 + 216*d2*b3*a() + - 72*d()*b4*c() + 72*c4*b()*a() + 216*d()*c3*a2 - 432*d()*b2*a()*c2 - - 729*a4*d2 - 2*b6 - + 18*b4*ac - 27*a2*b2*c2 + 6*b5*c() - 648*b2*c()*a2*d() - + 162*a()*d()*b4 + 1458*a3*d2*b(); + bigint a = coeffs[0], b=coeffs[1], c=coeffs[2], d=coeffs[3]; + bigint b2=sqr(b); + bigint b3=b*b2; + bigint b4=b*b3; + bigint b5=b*b4; + bigint b6=b*b5; + bigint a2=sqr(a); + bigint a3=a*a2; + bigint a4=a*a3; + bigint c2=sqr(c); + bigint c3=c*c2; + bigint c4=c*c3; + bigint c5=c*c4; + bigint c6=c*c5; + bigint d2=sqr(d); + bigint d3=d*d2; + bigint d4=d*d3; + bigint ac=a*c, bd=b*d; + + return 108*b3*a2*d - 12*b4*c2 + 216*a2*c4 + 72*b5*d + 486*a3*c2*d + - 270*a2*c3*b + 90*b3*c2*a + 972*a2*d2*b2 - 216*b3*ac*d - + 648*bd*c2*a2 + 54*a3*c3 - 486*a3*bd*c + 16*c3*b3 - 216*d2*b3*a - + 72*d*b4*c - 72*c4*b*a - 216*d*c3*a2 + 432*d*b2*a*c2 + + 729*a4*d2 + 2*b6 + - 18*b4*ac + 27*a2*b2*c2 - 6*b5*c + 648*b2*c*a2*d + - 162*a*d*b4 - 1458*a3*d2*b; } // The quantity called C_3 in the paper, = Norm(h0+h1) and should be @@ -201,7 +354,7 @@ return 27*d*c3*a2 + (27*d2*b3 - 54*d*c2*b2 + 9*c4*b)*a + 9*d*c*b4 - 2*c3*b3; } -//#define DEBUG +//#define DEBUG_REDUCE bigcomplex cubic::hess_root() const { @@ -218,13 +371,17 @@ return gamma; } -int cubic::is_hessian_reduced() // for positive discriminant only +int cubic::is_hessian_reduced() +// for positive discriminant only +// The condition is -P < Q <= P < R or 0 <= Q <= P=R. { bigint P = p_semi(); bigint R = r_semi(); if (P>R) return 0; + // now P<=R bigint Q = q_semi(); if (Q>P) return 0; + // now Q<=P<=R if (P==R) return (Q>=0); return (Q>-P); } @@ -233,7 +390,7 @@ { int s=1; bigint k; m.reset(); -#ifdef DEBUG +#ifdef DEBUG_REDUCE cout<<"Using hess_reduce() on "<<(*this)<r_semi()) { s=1; invert(m); -#ifdef DEBUG +#ifdef DEBUG_REDUCE cout << "invert: " << (*this) << endl; #endif } } // Now we have -P <= Q < P <= R and test for boundary condition - if((p_semi()==r_semi()) && (q_semi()<0)) + if ((p_semi()==r_semi()) && (q_semi()<0)) { invert(m); -#ifdef DEBUG +#ifdef DEBUG_REDUCE cout << "Final inversion: " << (*this) << endl; #endif } - if(a()<0) negate(m); + normalise(m); } void cubic::mathews_reduce(unimod& m) @@ -278,7 +435,7 @@ if(mat_c1()<0) { s=1; invert(m); -#ifdef DEBUG +#ifdef DEBUG_REDUCE cout << "invert: " << (*this) << endl; #endif } @@ -288,7 +445,7 @@ { s=1; x_shift(k,m); -#ifdef DEBUG +#ifdef DEBUG_REDUCE cout << "Shift by "<0) { s=1; x_shift(plus1,m); -#ifdef DEBUG +#ifdef DEBUG_REDUCE cout << "Shift by +1: "<<(*this)<=0, i.e. h0<=h2 + bigint C2 = j_c2(); if (C2<0) return 0; - bigint C4 = j_c4(); // = N(h1)/8, not in JCM paper - if (C1==0) return (C4>=0); + // now C1, C2 >=0, i.e. h1<=h0<=h2 + if (is_zero(C1)) // i.e. h0=h2 + { + bigint C4 = j_c4(); // = N(h1)/8, not in JCM paper + return (C4>=0); // i.e. h1 >= 0 + } bigint C3 = j_c3(); - return (C3>0); + return (C3>0); // i.e. h1 > -h0 } void cubic::jc_reduce(unimod& m) { int s=1; bigint k, jc2, jc3; - bigint plus1, minus1; plus1=1; minus1=-1; + bigint plus1(one), minus1(-one); bigfloat alpha, ra, rb, rc, rd, h0, h1, h2; m.reset(); -#ifdef DEBUG +#ifdef DEBUG_REDUCE cout << "\nJC-reducing " << (*this) << "...\n"; - cout<<"C1="<=bb) + q+=1; +#ifdef DEBUG_REDUCE + cout << "[a=0] shift "<< (*this); +#endif + x_shift(q,m); +#ifdef DEBUG_REDUCE + cout << " by "< " << (*this) << endl; - cout<<"C1="<0) || (j_c3()<0)) + if ((j_c2()<0) || (j_c3()<=0)) // then -h0 < h1 <= h0 fails, so shift: { s=1; alpha = real_root(); -#ifdef DEBUG +#ifdef DEBUG_REDUCE cout<<"alpha = "< "<<(*this)< "<<(*this)< "<<(*this)< " << (*this) << endl; - cout<<"C1="< cubic::rational_roots() const { - vector co(coeffs, coeffs+4); - return roots(co); + return roots(coeffs); } @@ -546,10 +724,7 @@ int sU; bigfloat i3a, a23, ra2, rb2, D, D2, D3, D32, D4, Pmax, Pmin; - int sl2_equiv, gl2_equiv; - bigint ZERO(0), ONE(1); unimod m; - const unimod m1(ONE,ZERO,ZERO,-ONE); bigfloat third = 1/to_bigfloat(3); bigfloat const1 = 2 / sqrt(to_bigfloat(27)); @@ -591,6 +766,7 @@ if ((verbose>1) && glist.size()==0) cout<1) { @@ -641,6 +817,7 @@ if (verbose>1 && glist.size()==0) cout<1) cout<<"found "<1) @@ -666,9 +843,9 @@ cout<::const_iterator gi=glist.begin(); gi!=glist.end(); gi++) + for (auto gi=glist.begin(); gi!=glist.end(); gi++) { - cubic g = *gi, gneg; + cubic g = *gi; if(verbose) cout< coeffs; // will always have length 4 public: - void set(long a, long b, long c, long d) - {coeffs[0]=a; coeffs[1]=b; coeffs[2]=c; coeffs[3]=d;} - void set(const bigint& a, const bigint& b, const bigint& c, const bigint& d) - {coeffs[0]=a; coeffs[1]=b; coeffs[2]=c; coeffs[3]=d;} - void set(const cubic& q) - {coeffs[0]=q.coeffs[0]; coeffs[1]=q.coeffs[1]; - coeffs[2]=q.coeffs[2]; coeffs[3]=q.coeffs[3];} - cubic() - {init(); set(0,0,0,0);} - ~cubic(); + // void set(long a, long b, long c, long d) + // {coeffs = {BIGINT(a),BIGINT(b),BIGINT(c),BIGINT(d)};} + // void set(const bigint& a, const bigint& b, const bigint& c, const bigint& d) + // {coeffs = {a,b,c,d};} + // void set(const vector& abcd) + // {std::copy(abcd.begin(), abcd.end(), coeffs.begin());} + // void set(const cubic& q) + // {std::copy(q.coeffs.begin(), q.coeffs.end(), coeffs.begin());} + cubic() + {coeffs.resize(4, BIGINT(0));} cubic(const bigint& a, const bigint& b, const bigint& c, const bigint& d) - {init(); set(a,b,c,d);} + {coeffs = {a,b,c,d};} cubic(long a, long b, long c, long d) - {init(); set(a,b,c,d);} - cubic(const cubic& q) - {init(); set(q);} - void operator=(const cubic& g) {set(g);} - int operator==(const cubic& g) - {return ((coeffs[0]==g.coeffs[0]) && (coeffs[1]==g.coeffs[1]) && - (coeffs[2]==g.coeffs[2]) && (coeffs[3]==g.coeffs[3]));} - inline bigint coeff(int i) - {if((i>=0)&&(i<=3)) return coeffs[i]; else return coeffs[0];} + {coeffs = {BIGINT(a),BIGINT(b),BIGINT(c),BIGINT(d)};} + cubic(const vector& abcd) + {coeffs = abcd;} + cubic(const cubic& q) + {coeffs = q.coeffs;} + int operator==(const cubic& g) const + {return (coeffs==g.coeffs);} + inline bigint coeff(int i) + {if((i>=0)&&(i<=3)) return coeffs[i]; else return coeffs[0];} inline bigint operator[](int i) const - {if((i>=0)&&(i<=3)) return coeffs[i]; else return coeffs[0];} + {if((i>=0)&&(i<=3)) return coeffs[i]; else return coeffs[0];} inline bigint a(void) const {return coeffs[0];} inline bigint b(void) const {return coeffs[1];} inline bigint c(void) const {return coeffs[2];} @@ -87,22 +85,55 @@ os<<"["<1, or 0 < Re(z) <= +1/2 and |z|=1. Note that -I always fixes + // z, so does [0,-1;1,0] (of order 4) if z=i (covariant a multiple + // of (X^2+Y^2)), and [1,1;-1,0] (of order 6) if z=(1+sqrt(-3))/2 + // (covariant a multiple of (X^2+XY+Y^2). + void sl2_reduce(unimod& m); + + // for an sl2-reduced cubic, normalise w.r.t. <-I> (default) or + // or . Two reduced and normalised forms are SL(2,Z)-equivalent + // iff they are equal + void normalise(unimod& m); + + // Test for sl2/gl2-equivalence: + int sl2_equivalent(const cubic& G) const; + int gl2_equivalent(const cubic& G) const; + // Test for sl2/gl2-equivalence to one in a list: + int sl2_equivalent_in_list(const vector& Glist) const; + int gl2_equivalent_in_list(const vector& Glist) const; + + // affine roots of F mod q, assuming leading coefficient a() is nonzero: + vector roots_mod(const bigint& q) const; + + // Return 1 iff F has a projective root mod q: + int has_roots_mod(const bigint& q) const; + + // Mathews quantities for use when disc<0: bigint mat_c1() const { return d()*(d()-b())+a()*(c()-a());} bigint mat_c2() const { return a()*d() - (a()+b())*(a()+b()+c());} bigint mat_c3() const { return a()*d() + (a()-b())*(a()-b()+c());} + + // P, Q, R: coefficients of the Hessian, used for reduction when disc>0 bigint p_semi() const { return sqr(b())-3*a()*c(); } bigint q_semi() const @@ -111,9 +142,22 @@ { return sqr(c())-3*b()*d(); } bigint u_semi() const { return 2*b()*sqr(b()) + 27*sqr(a())*d() - 9*a()*b()*c();} + + // jc_reduce uses the algebraic real quadratic covariant with coeffs + // [h0,h1,h2]. A form is reduced if -h0=0, C2>=0, C3>0 and C4>=0 if C1=0. + + // jc_c1() is the quantity denoted C1=N(h2-h0) in the paper, its + // sign is sign(h2-h0) so is >=0 for a reduced form: bigint j_c1() const; + // jc_c2() is the quantity denoted C2=N(h0-h1) in the paper, its + // sign is sign(h0-h1) so is >=0 for a reduced form: bigint j_c2() const; + // jc_c3() is the quantity denoted C3=N(h0+h1) in the paper, its + // sign is sign(h1+h0) so is >0 for a reduced form: bigint j_c3() const; + // jc_c4() is not in the paper, its sign is sign(h1) so is >=0 for a + // reduced form with C1=0, i.e. with h0=h1: bigint j_c4() const; bigcomplex hess_root() const; @@ -130,6 +174,21 @@ {return ((a()==0) || (rational_roots().size()>0));} int is_irreducible() const {return ((a()!=0) && (rational_roots().size()==0));} + + bigint content() const + { + return gcd(gcd(gcd(coeffs[0], coeffs[1]), coeffs[2]), coeffs[3]); + } + int is_primitive() const + { + return content()==1; + } + // divide by a constant factor (which should divide all the coefficients) + cubic operator/(const bigint& g) const + { + return cubic(coeffs[0]/g, coeffs[1]/g, coeffs[2]/g, coeffs[3]/g); + } + }; @@ -138,11 +197,12 @@ return os<<"["< reduced_cubics(const bigint& disc, int include_reducibles=1, int gl2=0, int verbose=0); diff -Nru eclib-20220621/progs/in/list_cubics.in eclib-20221012/progs/in/list_cubics.in --- eclib-20220621/progs/in/list_cubics.in 2022-06-21 14:25:26.000000000 +0000 +++ eclib-20221012/progs/in/list_cubics.in 2022-10-12 11:04:09.000000000 +0000 @@ -1,4 +1,4 @@ -0 +0 0 200 0 0 200 0 1 200 1 0 diff -Nru eclib-20220621/progs/list_cubics.cc eclib-20221012/progs/list_cubics.cc --- eclib-20220621/progs/list_cubics.cc 2022-06-21 14:25:26.000000000 +0000 +++ eclib-20221012/progs/list_cubics.cc 2022-10-12 11:04:09.000000000 +0000 @@ -1,7 +1,7 @@ -// list_cubics.cc: Program for listing integer cubics with given discriminant bound +// list_cubics.cc: Program for listing integer cubics with given or bounded discriminant ////////////////////////////////////////////////////////////////////////// // -// Copyright 1990-2012 John Cremona +// Copyright 1990-2022 John Cremona // // This file is part of the eclib package. // @@ -26,61 +26,85 @@ #include #include +void output_cubics(const bigint& disc, int include_reducibles=1, int gl2=0, int verbose=0) +{ + vector glist = reduced_cubics(disc, include_reducibles, gl2, verbose); + if (glist.size()==0) + { + if(verbose>1) + { + cout<< "No "; + if (!include_reducibles) cout << "irreducible"; + cout << " cubics with discriminant " << disc << endl; + } + } + else + { + if (verbose) + cout << glist.size() << " with discriminant "; + cout << disc << " " << glist << endl; + } +} + int main() { initprimes("PRIMES"); bigint disc, absdisc, maxdisc; - int neg; - int verbose=0, include_reducibles=1, gl2=0; + int neg, verbose=0, include_reducibles=1, gl2=0, single=0; + cerr << "Verbosity level (0, 1 or 2): "; cin >> verbose; - while(cerr << "Enter discriminant bound (positive or negative, 0 to stop): ", cin >> maxdisc, !is_zero(maxdisc)) + cerr << "Single discriminants (1) or a range (0)? "; + cin >> single; + + if (single) { cerr << "Include reducible cubics? (0 or 1): "; cin >> include_reducibles; + cerr << "Use GL(2,Z)-equivalence instead of SL(2,Z)? (0 or 1): "; cin >> gl2; - neg=(maxdisc<0); - if (include_reducibles) - cout << "Cubics with "; - else - cout << "Irreducible cubics with "; - if(neg) - { - ::negate(maxdisc); - cout << "negative discriminant down to -"; - } - else - { - cout << "positive discriminant up to "; - } - cout << maxdisc; - cout << " up to " << (gl2?"GL":"SL") << "(2,Z)-equivalence"; - cout << endl; - for(absdisc=1; absdisc<=maxdisc; absdisc++) + while(cerr << "Enter discriminant (positive or negative, 0 to stop): ", cin >> disc, !is_zero(disc)) + { + if (verbose) + { + cout << (include_reducibles? "Cubics": "Irreducible cubics") << " with discriminant "; + cout << disc; + cout << " up to " << (gl2?"GL":"SL") << "(2,Z)-equivalence"; + cout << endl; + } + output_cubics(disc, include_reducibles, gl2, verbose); + } + } + else + { + while(cerr << "Enter discriminant bound (positive or negative, 0 to stop): ", cin >> maxdisc, !is_zero(maxdisc)) + { + cerr << "Include reducible cubics? (0 or 1): "; + cin >> include_reducibles; + + cerr << "Use GL(2,Z)-equivalence instead of SL(2,Z)? (0 or 1): "; + cin >> gl2; + + neg=(maxdisc<0); + + cout << (include_reducibles? "Cubics with ": "Irreducible cubics with "); + cout << (neg? "negative discriminant down to ": "positive discriminant up to "); + cout << maxdisc; + cout << ", up to " << (gl2? "GL": "SL") << "(2,Z)-equivalence"; + cout << endl; + + if (neg) + ::negate(maxdisc); + for(absdisc=1; absdisc<=maxdisc; absdisc++) { disc=absdisc; if(neg) ::negate(disc); - vector glist = reduced_cubics(disc, include_reducibles, gl2, verbose); - if (glist.size()==0) - { - if(verbose>1) - { - cout<< "No "; - if (!include_reducibles) cout << "irreducible"; - cout << " cubics with discriminant " << disc << endl; - } - } - else - { - cout << glist.size(); - cout << " with discriminant " << disc; - if (glist.size()>0) cout<< " : " << glist; - cout << endl; - } + output_cubics(disc, include_reducibles, gl2, verbose); } - cout< #include #include +#include int main() { initprimes("PRIMES"); bigint a, b, c, d, disc; unimod m; + while(cout << "Enter cubic coeffs a, b, c, d: ", cin >> a >> b >> c >> d, !(is_zero(a)&&is_zero(b)&&is_zero(c)&&is_zero(d))) @@ -44,7 +46,9 @@ if(disc>0) { cout << "Using Hessian to reduce...\n"; + m.reset(); g.hess_reduce(m); + assert (transform(g0,m)==g); cout << "Hessian reduced cubic = "<