Skip to content

Commit

Permalink
Merge pull request #75 from JohnCremona/bigcomplex
Browse files Browse the repository at this point in the history
Implementation of bigcomplex class
  • Loading branch information
JohnCremona authored Apr 24, 2023
2 parents 6639c6b + 7558e4a commit 0e8a341
Show file tree
Hide file tree
Showing 20 changed files with 220 additions and 170 deletions.
7 changes: 4 additions & 3 deletions configure.ac
Original file line number Diff line number Diff line change
Expand Up @@ -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], [[email protected]])
AC_INIT([eclib], [20230424], [[email protected]])

AM_INIT_AUTOMAKE([-Wall])
AC_MSG_NOTICE([Configuring eclib...])
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion libsrc/Makefile.am
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
2 changes: 1 addition & 1 deletion libsrc/cperiods.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand Down
2 changes: 1 addition & 1 deletion libsrc/cubic.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down
6 changes: 3 additions & 3 deletions libsrc/curvedata.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
155 changes: 155 additions & 0 deletions libsrc/eclib/bigcomplex.h
Original file line number Diff line number Diff line change
@@ -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<<re;
return s.str();
}
if (IsImaginary())
{
if (IsOne(im))
s << "";
else
if (IsOne(-im))
s << "-";
else s << im;
}
else
{
s<<re;
if(im>0)
s<<"+";
else
s<<"-";
RR abs_im(::NTL::abs(im));
if (abs_im>1)
s<<abs_im;
}
s<<"i";
return s.str();
}

friend ostream& operator<<(ostream& s, const bigcomplex& z)
{
s << string(z);
return s;
}

// Implementation
private:
RR re, im;
};


// inline functions
inline RR real(const bigcomplex& z) {return z.real();};
inline RR imag(const bigcomplex& z) {return z.imag();};
inline RR norm(const bigcomplex& z) {return z.norm();};
inline RR arg(const bigcomplex& z) { return z.arg();};
inline RR abs(const bigcomplex& z) {return z.abs();};
inline bigcomplex conj(const bigcomplex& z) {return z.conj();};
inline int IsReal(const bigcomplex& z) {return z.IsReal();}
inline int IsImaginary(const bigcomplex& z) {return z.IsImaginary();}
inline int IsZero(const bigcomplex& z) {return z.IsZero();}

// operators with left operand RR
inline bigcomplex operator+(const RR& r, const bigcomplex& z) {return bigcomplex(r)+z;};
inline bigcomplex operator-(const RR& r, const bigcomplex& z) {return bigcomplex(r)-z;};
inline bigcomplex operator*(const RR& r, const bigcomplex& z) {return bigcomplex(r)*z;};
inline bigcomplex operator/(const RR& r, const bigcomplex& z) {return bigcomplex(r)/z;};

// standard method functions as functions:
inline bigcomplex sqrt(const bigcomplex& z) {bigcomplex w(z); return w.sqrt();};
inline bigcomplex cos(const bigcomplex& z) {bigcomplex w(z); return w.cos();};
inline bigcomplex sin(const bigcomplex& z) {bigcomplex w(z); return w.sin();};
inline bigcomplex exp(const bigcomplex& z) {bigcomplex w(z); return w.exp();};
inline bigcomplex log(const bigcomplex& z) {bigcomplex w(z); return w.log();};


inline void swap(bigcomplex& z1, bigcomplex& z2) {bigcomplex z3(z1); z1=z2; z2=z3;}

#endif
6 changes: 6 additions & 0 deletions libsrc/eclib/compproc.h
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,13 @@ void getc4c6(const bigcomplex& w1, const bigcomplex& w2, bigcomplex& c4, bigcomp
bigcomplex discriminant(const bigcomplex& b, const bigcomplex& c, const bigcomplex& d);

vector<bigcomplex> solvecubic(const bigcomplex& c1, const bigcomplex& c2, const bigcomplex& c3);
inline vector<bigcomplex> solvecubic(const bigfloat& c1, const bigfloat& c2, const bigfloat& c3)
{return solvecubic(bigcomplex(c1), bigcomplex(c2), bigcomplex(c3));}

vector<bigcomplex> solvequartic(const bigcomplex& a, const bigcomplex& b, const bigcomplex& c, const bigcomplex& d);
inline vector<bigcomplex> solvequartic(const bigfloat& a, const bigfloat& b, const bigfloat& c, const bigfloat& d)
{return solvequartic(bigcomplex(a), bigcomplex(b), bigcomplex(c), bigcomplex(d));}

vector<bigcomplex> 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);
Expand Down
4 changes: 2 additions & 2 deletions libsrc/eclib/htconst.h
Original file line number Diff line number Diff line change
Expand Up @@ -163,9 +163,9 @@ class CurveHeightConst : public CurveRed, Cperiods {
bigfloat psi(const bigfloat& x);
vector<bigfloat> ordinates(const bigfloat& x);
vector<bigcomplex> 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(); }
Expand Down
79 changes: 6 additions & 73 deletions libsrc/eclib/interface.h
Original file line number Diff line number Diff line change
Expand Up @@ -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_
Expand Down Expand Up @@ -162,75 +162,8 @@ inline int is_approx_zero(const RR& x)
return abs(x.mantissa())<power2_ZZ(-n);
}
} // namespace NTL
#ifdef _LIBCPP_VERSION
namespace NTL {
inline bool isinf(const RR& z) {
return false;
}
inline bool isnan(const RR& z) {
return false;
}
inline RR copysign(const RR& x, const RR& y) {
if (sign(x) != sign(y)) {
return -y;
}
return y;
}
inline bool signbit(const RR& x) {
return sign(x) < 0;
}
inline RR fmax(const RR& x, const RR& y)
{
return x < y ? y : x;
}
} // namespace NTL
#endif

#include <complex>
typedef complex<RR> 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
Expand Down Expand Up @@ -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();}


//////////////////////////////////////////////////////////////////
Expand Down
Loading

0 comments on commit 0e8a341

Please sign in to comment.