diff --git a/configure.ac b/configure.ac index a41c8f75..0405f117 100644 --- a/configure.ac +++ b/configure.ac @@ -7,7 +7,7 @@ AC_PREREQ([2.65]) # 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], [20221012], [john.cremona@gmail.com]) +AC_INIT([eclib], [20230424], [john.cremona@gmail.com]) AM_INIT_AUTOMAKE([-Wall]) AC_MSG_NOTICE([Configuring eclib...]) @@ -40,9 +40,10 @@ LT_INIT # # NB The suffix of the library name (libec.so here) is (c-a).a.r -LT_CURRENT=11 +LT_CURRENT=12 LT_REVISION=0 -LT_AGE=1 +LT_AGE=2 + AC_SUBST(LT_CURRENT) AC_SUBST(LT_REVISION) AC_SUBST(LT_AGE) diff --git a/libsrc/Makefile.am b/libsrc/Makefile.am index 3d1a2226..f0397404 100644 --- a/libsrc/Makefile.am +++ b/libsrc/Makefile.am @@ -15,7 +15,7 @@ LIBS = $(FLINT_LIBS) $(NTL_LIBS) $(PARI_LIBS) $(BOOST_LIBS) $(PTHREAD_LIB lib_LTLIBRARIES = libec.la -PROCS_DOTHS = eclib/interface.h eclib/templates.h eclib/arith.h eclib/xmod.h eclib/marith.h eclib/compproc.h eclib/vec.h eclib/vector.h eclib/mat.h eclib/matrix.h eclib/sub.h eclib/subspace.h eclib/rat.h eclib/bigrat.h eclib/kbessel.h eclib/svec.h eclib/svector.h eclib/smat.h eclib/smatrix.h eclib/smat_elim.h eclib/smatrix_elim.h eclib/mvector.h eclib/mmatrix.h eclib/msubspace.h eclib/method.h eclib/splitbase.h eclib/xsplit.h eclib/conic.h eclib/legendre.h eclib/quadratic.h eclib/unimod.h eclib/illl.h eclib/hilbert.h eclib/timer.h eclib/cubic.h eclib/gf.h eclib/polys.h eclib/realroots.h eclib/parifact.h eclib/p2points.h eclib/xsplit_data.h eclib/threadpool.h eclib/logger.h eclib/types.h +PROCS_DOTHS = eclib/interface.h eclib/templates.h eclib/arith.h eclib/xmod.h eclib/marith.h eclib/bigcomplex.h eclib/compproc.h eclib/vec.h eclib/vector.h eclib/mat.h eclib/matrix.h eclib/sub.h eclib/subspace.h eclib/rat.h eclib/bigrat.h eclib/kbessel.h eclib/svec.h eclib/svector.h eclib/smat.h eclib/smatrix.h eclib/smat_elim.h eclib/smatrix_elim.h eclib/mvector.h eclib/mmatrix.h eclib/msubspace.h eclib/method.h eclib/splitbase.h eclib/xsplit.h eclib/conic.h eclib/legendre.h eclib/quadratic.h eclib/unimod.h eclib/illl.h eclib/hilbert.h eclib/timer.h eclib/cubic.h eclib/gf.h eclib/polys.h eclib/realroots.h eclib/parifact.h eclib/p2points.h eclib/xsplit_data.h eclib/threadpool.h eclib/logger.h eclib/types.h QCURVES_DOTHS = eclib/curve.h eclib/points.h eclib/cperiods.h eclib/isogs.h eclib/reader.h eclib/mwprocs.h eclib/lambda.h eclib/sifter.h eclib/sieve_search.h eclib/htconst.h eclib/egr.h eclib/saturate.h eclib/divpol.h eclib/pointsmod.h eclib/curvemod.h eclib/ffmod.h eclib/tlss.h eclib/elog.h eclib/getcurve.h diff --git a/libsrc/cperiods.cc b/libsrc/cperiods.cc index 0f254798..4c19ce51 100644 --- a/libsrc/cperiods.cc +++ b/libsrc/cperiods.cc @@ -329,7 +329,7 @@ void Cperiods::store_sums() } w1squared = w1*w1; w1cubed = w1*w1squared; - bigcomplex term = one, qtm = qtau; + bigcomplex term(one), qtm(qtau); sum3=to_bigfloat(0); for (bigfloat m=to_bigfloat(1); ! SMALL(term); m+=1) { diff --git a/libsrc/cubic.cc b/libsrc/cubic.cc index 757fe9fc..5722c6cc 100644 --- a/libsrc/cubic.cc +++ b/libsrc/cubic.cc @@ -362,7 +362,7 @@ bigcomplex cubic::hess_root() const if(!is_positive(disc())) { cout<<"Error: hess_root called with negative dicriminant!\n"; - return to_bigfloat(0); + return bigcomplex(); // 0 } bigfloat P = I2bigfloat(p_semi()); bigfloat Q = I2bigfloat(q_semi()); diff --git a/libsrc/curvedata.cc b/libsrc/curvedata.cc index bebdb288..7894f28c 100644 --- a/libsrc/curvedata.cc +++ b/libsrc/curvedata.cc @@ -321,9 +321,9 @@ Curvedata opt_x_shift(const Curvedata& C, bigint& k) C.getbi(b2,b4,b6,b8); cubic b_cubic(four,b2,2*b4,b6); k = b_cubic.shift_reduce(); - Curvedata CC(C); - CC.transform(k,zero,zero); - return CC; + Curvedata CD(C); + CD.transform(k,zero,zero); + return CD; } #if(0) diff --git a/libsrc/eclib/bigcomplex.h b/libsrc/eclib/bigcomplex.h new file mode 100644 index 00000000..84b14429 --- /dev/null +++ b/libsrc/eclib/bigcomplex.h @@ -0,0 +1,155 @@ +// bigcomplex.h: complex class built on NTL's RR +////////////////////////////////////////////////////////////////////////// +// +// Copyright 1990-2023 John Cremona +// +// This file is part of the eclib package. +// +// eclib 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. +// +// eclib 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 eclib; if not, write to the Free Software Foundation, +// Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA +// +////////////////////////////////////////////////////////////////////////// + +#ifndef _ECLIB_BIGCOMPLEX_H_ +#define _ECLIB_BIGCOMPLEX_H_ + +class bigcomplex{ +public: + explicit bigcomplex() : re(RR()), im(RR()) {;} + explicit bigcomplex(const RR& r) : re(r), im(RR()) {;} + explicit bigcomplex(const RR& r, const RR& i) : re(r), im(i) {;} + bigcomplex(const bigcomplex&) = default; + + // standard class methods + RR real() const {return re;}; + RR imag() const {return im;}; + RR norm() const {return sqr(re)+sqr(im);}; + RR arg() const { return atan2(im, re);}; + RR abs() const {return ::NTL::sqrt(norm());}; + bigcomplex conj() const {return bigcomplex(re,-im);}; + bigcomplex timesI() const {return bigcomplex(-im,re);}; + int IsReal() const {return ::NTL::IsZero(im);} + int IsImaginary() const {return ::NTL::IsZero(re);} + int IsZero() const {return ::NTL::IsZero(re) && ::NTL::IsZero(im);} + + + bigcomplex operator= (const RR& r) {re=r; im=RR(); return *this;} + bigcomplex operator+=(const RR& r) {re+=r; return *this;}; + bigcomplex operator+(const RR& r) const {bigcomplex z(*this); z.re+=r; return z;}; + bigcomplex operator-=(const RR& r) {re -= r; return *this;}; + bigcomplex operator-(const RR& r) const {bigcomplex z(*this); z.re-=r; return z;}; + bigcomplex operator*=(const RR& r) {re*=r; im*=r; return *this;}; + bigcomplex operator*(const RR& r) const {bigcomplex z(*this); z.re*=r; z.im*=r; return z;}; + bigcomplex operator/=(const RR& r) {re/=r; im/=r; return *this;}; + bigcomplex operator/(const RR& r) const {bigcomplex z(*this); z.re/=r; z.im/=r; return z;}; + + bigcomplex operator=(const bigcomplex& z) {re=z.re; im=z.im; return *this;}; + bigcomplex operator+=(const bigcomplex& z) {re+=z.re; im+=z.im; return *this;}; + bigcomplex operator+(const bigcomplex& z) const {return bigcomplex(re+z.re, im+z.im);}; + bigcomplex operator+() {return *this;}; + bigcomplex operator-=(const bigcomplex& z) {re-=z.re; im-=z.im; return *this;}; + bigcomplex operator-(const bigcomplex& z) const {return bigcomplex(re-z.re, im-z.im);}; + bigcomplex operator-() const {return bigcomplex(-re,-im);}; + bigcomplex operator*=(const bigcomplex& z) {RR x = re; re=x*z.re-im*z.im; im = x*z.im+im*z.re; return *this;}; + bigcomplex operator*(const bigcomplex& z) const {bigcomplex w(*this); w*=z; return w;}; + bigcomplex operator/=(const bigcomplex& z) {RR r=z.norm(), x = re; re=(x*z.re+im*z.im)/r; im = (-x*z.im+im*z.re)/r; return *this;}; + bigcomplex operator/(const bigcomplex& z) const {bigcomplex w(*this); w/=z; return w;}; + + // standard functions as class methods: + bigcomplex sqrt() const {RR r = abs(); RR u = ::NTL::sqrt((r+re)/2), v=::NTL::sqrt((r-re)/2); if (im<0) v=-v; return bigcomplex(u,v);}; + bigcomplex cos() const {bigcomplex z(this->timesI()); return (z+z.conj())/to_RR(2);}; + bigcomplex sin() const {bigcomplex z(this->timesI()); return (z-z.conj())/bigcomplex(RR(),to_RR(2));}; + bigcomplex exp() const {RR e = ::NTL::exp(re); return bigcomplex(e * ::NTL::cos(im), e * ::NTL::sin(im));}; + bigcomplex log() const {return bigcomplex(::NTL::log(abs()), arg());} + + + operator string() const + { + ostringstream s; + s << "(" << re << "," << im << ")"; + return s.str(); + } + + string pretty_string() const + { + ostringstream s; + if (IsReal()) + { + s<0) + s<<"+"; + else + s<<"-"; + RR abs_im(::NTL::abs(im)); + if (abs_im>1) + s< solvecubic(const bigcomplex& c1, const bigcomplex& c2, const bigcomplex& c3); +inline vector solvecubic(const bigfloat& c1, const bigfloat& c2, const bigfloat& c3) +{return solvecubic(bigcomplex(c1), bigcomplex(c2), bigcomplex(c3));} + vector solvequartic(const bigcomplex& a, const bigcomplex& b, const bigcomplex& c, const bigcomplex& d); +inline vector solvequartic(const bigfloat& a, const bigfloat& b, const bigfloat& c, const bigfloat& d) +{return solvequartic(bigcomplex(a), bigcomplex(b), bigcomplex(c), bigcomplex(d));} + vector solverealquartic(const bigfloat& a, const bigfloat& b, const bigfloat& c, const bigfloat& d, const bigfloat& e); void quadsolve(const bigfloat& p, const bigfloat& q, bigcomplex& root1,bigcomplex& root2); diff --git a/libsrc/eclib/htconst.h b/libsrc/eclib/htconst.h index bfb40f38..95292d43 100644 --- a/libsrc/eclib/htconst.h +++ b/libsrc/eclib/htconst.h @@ -163,9 +163,9 @@ class CurveHeightConst : public CurveRed, Cperiods { bigfloat psi(const bigfloat& x); vector ordinates(const bigfloat& x); vector ztopoint(const bigcomplex& z) - {return ellztopoint(z,I2bigfloat(a1),I2bigfloat(a2),I2bigfloat(a3));} + {return ellztopoint(z, bigcomplex(I2bigfloat(a1)), bigcomplex(I2bigfloat(a2)), bigcomplex(I2bigfloat(a3)));} bigcomplex pointtoz(const bigfloat& x, const bigfloat& y) - {return ellpointtoz(*this,*this,x,y);} + {return ellpointtoz(*this, *this, x, y);} public: CurveHeightConst(CurveRed& CR); void compute() {compute_phase1(); compute_phase2(); } diff --git a/libsrc/eclib/interface.h b/libsrc/eclib/interface.h index ba5a5aa4..b8500dc0 100644 --- a/libsrc/eclib/interface.h +++ b/libsrc/eclib/interface.h @@ -24,7 +24,7 @@ // The macro NO_MPFP can optionally be set. It controls whether most // real and complex floating point functions are implemented using // doubles and complex doubles (if set) or using NTL bigfloats -// (RR) and bigcomplexes (CC) (if not set, the default) +// (RR) and bigcomplexes (if not set, the default) #ifndef _ECLIB_INTERFACE_H_ #define _ECLIB_INTERFACE_H_ @@ -162,75 +162,8 @@ inline int is_approx_zero(const RR& x) return abs(x.mantissa()) -typedef complex CC; -typedef CC bigcomplex; - -#ifdef _LIBCPP_VERSION -template <> inline RR std::abs(const CC &z) -{ - RR re = z.real(); - RR im = z.imag(); - return sqrt(re*re + im*im); -} -template <> inline CC std::exp(const CC &z) -{ - RR im = z.imag(); - RR e = exp(z.real()); - return CC(e * cos(im), e * sin(im)); -} -inline CC operator/(const CC &a, const CC &b) -{ - RR are = a.real(); - RR aim = a.imag(); - RR bre = b.real(); - RR bim = b.imag(); - if (abs(bre) <= abs(bim)) { - RR r = bre / bim; - RR den = bim + r*bre; - return CC((are*r + aim)/den, (aim*r - are)/den); - } else { - RR r = bim / bre; - RR den = bre + r*bim; - return CC((are + aim*r)/den, (aim - are*r)/den); - } -} -inline CC operator /(const RR &a, const CC &b) -{ - CC r(a); - return r/b; -} -inline CC &operator /=(CC &a, const CC &b) -{ - a = a/b; - return a; -} - -#endif +#include "bigcomplex.h" // NB Internal precision is always bit-precision. We used to use // decimal precision in the user interface with a conversion factor of @@ -271,10 +204,10 @@ inline int is_zero(bigfloat x) {return IsZero(x);} inline int is_zero(bigcomplex z) {return IsZero(z.real()) && IsZero(z.imag());} inline void Iasb(bigint& a, bigfloat x) {RoundToZZ(a,x);} inline void Iasb(long& a, bigfloat x) {ZZ n; RoundToZZ(n,x); a=I2long(n);} -istream& operator>>(istream& is, CC& z); -inline CC pow(const CC& a, int e) {return exp(to_RR(e)*log(a));} -inline CC pow(const CC& a, long e) {return exp(to_RR(e)*log(a));} -inline CC pow(const CC& a, const RR& e) {return exp(e*log(a));} +istream& operator>>(istream& is, bigcomplex& z); +inline bigcomplex pow(const bigcomplex& a, int e) {return (to_RR(e)*a.log()).exp();} +inline bigcomplex pow(const bigcomplex& a, long e) {return (to_RR(e)*a.log()).exp();} +inline bigcomplex pow(const bigcomplex& a, const RR& e) {return (e*a.log()).exp();} ////////////////////////////////////////////////////////////////// diff --git a/libsrc/interface.cc b/libsrc/interface.cc index 14b86e1c..c05fdfe2 100644 --- a/libsrc/interface.cc +++ b/libsrc/interface.cc @@ -275,7 +275,7 @@ RR atan2 (const RR & y, const RR & x) // The template version requires an automatic conversion from 0 to an // RR, so cannot be used as is. We have manually instantiated it // here. -istream& operator>>(istream& is, CC& z) +istream& operator>>(istream& is, bigcomplex& z) { RR r, i; char c; @@ -287,12 +287,12 @@ istream& operator>>(istream& is, CC& z) { is >> i >> c; if (c == ')') - z = CC(r, i); + z = bigcomplex(r, i); else is.setstate(ios_base::failbit); } else if (c == ')') - z = CC(r, to_RR(0)); + z = bigcomplex(r, to_RR(0)); else is.setstate(ios_base::failbit); } @@ -300,7 +300,7 @@ istream& operator>>(istream& is, CC& z) { is.putback(c); is >> r; - z = CC(r, to_RR(0)); + z = bigcomplex(r, to_RR(0)); } return is; } @@ -308,24 +308,14 @@ istream& operator>>(istream& is, CC& z) // return value is 1 for success, else 0 int longify(const bigfloat& x, long& a, int rounding) { - if ((x<=MAXLONG)&&(x>=MINLONG)) - { - switch(rounding) - { - case 0: // round to nearest - default: - a = I2long(RoundToZZ(x)); return 1; - case 1: // round up - a = I2long(CeilToZZ(x)); return 1; - case -1: // round down - a = I2long(FloorToZZ(x)); return 1; - } - } - else + bigint z = (rounding==0? RoundToZZ(x): (rounding==1? CeilToZZ(x): FloorToZZ(x))); + if(!is_long(z)) { cerr<<"Attempt to convert "<>nro>>nco; - long size=nro*nco; - delete[] entries; - entries = new scalar[size]; - scalar* m=entries; - while(size--) fin>>*m++; - fin.close(); -} - -#endif - istream& operator>>(istream& s, mat& m) { long n=m.nro*m.nco; diff --git a/libsrc/mwprocs.cc b/libsrc/mwprocs.cc index b5c21f69..22f370e2 100644 --- a/libsrc/mwprocs.cc +++ b/libsrc/mwprocs.cc @@ -59,9 +59,9 @@ vector roots_of_cubic(const Curve& E) ra4=I2bigfloat(a4), ra6=I2bigfloat(a6); - bigcomplex c1 = ra2 + ra1*(ra1/4) ; - bigcomplex c2 = ra4 + ra1*(ra3/2) ; - bigcomplex c3 = ra6 + ra3*(ra3/4) ; + bigfloat c1 = ra2 + ra1*(ra1/4); + bigfloat c2 = ra4 + ra1*(ra3/2); + bigfloat c3 = ra6 + ra3*(ra3/4); return solvecubic(c1,c2,c3); } diff --git a/libsrc/saturate.cc b/libsrc/saturate.cc index 26883ddc..90d38891 100644 --- a/libsrc/saturate.cc +++ b/libsrc/saturate.cc @@ -303,10 +303,12 @@ int saturator::enlarge() bigint old_index_bound = the_index_bound; set_index_bound(); if(verbose) - if (the_index_bound < old_index_bound) - cout << "Reducing index bound from " << old_index_bound <<" to " << the_index_bound << endl; - else - cout << "The index bound " << the_index_bound << " has not changed"<> i; -// The following should work, doesn't in TURBO C++, but DOES with GCC!! - { cout << "Creating an array of 3 matrices\n"; mat* matlist = new mat[3]; diff --git a/tests/mmattest.cc b/tests/mmattest.cc index 715b8664..6177c03a 100644 --- a/tests/mmattest.cc +++ b/tests/mmattest.cc @@ -42,8 +42,6 @@ int main(void) cout << "A = " << a; } -// The following should work, doesn't in TURBO C++, but DOES with GCC!! - { cout << "Creating an array of 3 matrices\n"; mat_m* matlist = new mat_m[3]; diff --git a/tests/out_no_ntl/tversion.out b/tests/out_no_ntl/tversion.out index 5d997a7f..a07953da 100644 --- a/tests/out_no_ntl/tversion.out +++ b/tests/out_no_ntl/tversion.out @@ -1,12 +1,12 @@ -This eclib version is 20221012 -[ year month day ] = [ 2022 10 12 ] +This eclib version is 20230424 +[ year month day ] = [ 2023 4 24 ] Configure options and compilation date: -eclib version 20221012, using NTL bigints but no multiprecision floating point +eclib version 20230424, using NTL bigints but no multiprecision floating point Testing comparison of version date and various (year, month, day) date triples: -Date 20200101 is before the version date 20221012 -Date 20300611 is after the version date 20221012 -Date 20210317 is before the version date 20221012 +Date 20200101 is before the version date 20230424 +Date 20300611 is after the version date 20230424 +Date 20210317 is before the version date 20230424 diff --git a/tests/out_ntl/telog.out b/tests/out_ntl/telog.out index a81a8ec1..568d8fc1 100644 --- a/tests/out_ntl/telog.out +++ b/tests/out_ntl/telog.out @@ -54,10 +54,10 @@ Checking... Curve [1,1,0,-202,1025] The point Q is = [-8:51:1] -Elliptic log of P is (2.15251566696160241628533178186449561190793579506,0) +Elliptic log of P is (2.15251566696160241628533178186449561190793579504,0) Reconstructed P = [-8:51:1] The point P3 is = [8:-3:1] -Elliptic log of P is (3.40172041025849963266962742312340879857800215119,0) +Elliptic log of P is (3.40172041025849963266962742312340879857800215114,0) Reconstructed P = [8:-3:1] (m*[8:-3:1])/m = [ [8:-3:1] ] Checking... diff --git a/tests/out_ntl/tversion.out b/tests/out_ntl/tversion.out index d40faa8e..040ff448 100644 --- a/tests/out_ntl/tversion.out +++ b/tests/out_ntl/tversion.out @@ -1,12 +1,12 @@ -This eclib version is 20221012 -[ year month day ] = [ 2022 10 12 ] +This eclib version is 20230424 +[ year month day ] = [ 2023 4 24 ] Configure options and compilation date: -eclib version 20221012, using NTL bigints and NTL real and complex multiprecision floating point +eclib version 20230424, using NTL bigints and NTL real and complex multiprecision floating point Testing comparison of version date and various (year, month, day) date triples: -Date 20200101 is before the version date 20221012 -Date 20300611 is after the version date 20221012 -Date 20210317 is before the version date 20221012 +Date 20200101 is before the version date 20230424 +Date 20300611 is after the version date 20230424 +Date 20210317 is before the version date 20230424