improve precision handling for special L-value computations

JohnCremona · Apr 22, 2024 · d209a2d · d209a2d
1 parent 3ea1296
commit d209a2d
Show file tree

Hide file tree

Showing 8 changed files with 135 additions and 202 deletions.
diff --git a/libsrc/eclib/interface.h b/libsrc/eclib/interface.h
@@ -152,7 +152,6 @@ inline RR sinh(const RR& x) {return (exp(x)-exp(-x))/2;}
 inline RR tan(const RR& x) {return sin(x)/cos(x);}
 RR atan2(const RR&, const RR&);
 inline int is_approx_zero(const RR& x)
-//  {return abs(x)<power2_RR(2-RR::precision());}
 {
   if (IsZero(x)) return 1;
   long n = x.exponent()+RR::precision()-1;
@@ -221,7 +220,7 @@ typedef double bigfloat;
 inline long decimal_precision() {return 14;}
 inline long bit_precision() {return 52;}
 inline int is_zero(double x) {return fabs(x)<1e-14;}
-inline int is_approx_zero(double x) {return fabs(x)<1e-10;}
+inline int is_approx_zero(double x) {return fabs(x)<1e-11;}
 inline int sign(double x) {return (double(0) < x) - (x < double(0));}
 
 // We cannot set internal bit precision in this mode, so we just set the output decimal precision

diff --git a/libsrc/eclib/periods.h b/libsrc/eclib/periods.h
@@ -26,15 +26,6 @@
 
 #include <eclib/newforms.h>
 
-#define PI Pi()
-const bigfloat eps = to_bigfloat(1.0e-16);   // ?? mindouble;
-#define EULERGAMMA Euler()
-
-bigfloat myg0(bigfloat x);
-bigfloat myg1(bigfloat x);
-bigfloat myg2(bigfloat x);
-bigfloat myg3(bigfloat x);
-
 class character {
 private:
   long modul;
@@ -123,8 +114,6 @@ class part_period :public summer {
   bigcomplex getperiod() {return bigcomplex(rp,ip);}
 };
 
-bigfloat G(int r, bigfloat x);  // G_r(x)
-
 class ldash1 : public summer {
 private:
   long r;
@@ -133,8 +122,7 @@ class ldash1 : public summer {
   bigfloat G(bigfloat x);  // G_r(x)
   void init(const level* N, const vector<long>& f_aplist, long f_sfe, const rational& f_loverp);
   void use(long n, long an) override {use1(n,an);}
-  bigfloat func1(long n) override
-  { return -G(factor1*to_bigfloat(n)); }
+  bigfloat func1(long n) override;
 public:
   ldash1 (const level* iN, const newform* f);
   ldash1 (const newforms* nf, long i);  // the i'th newform
@@ -167,5 +155,6 @@ vector<long> resort_aplist(const level* iN,
 			   const vector<long>& primelist,
 			   const vector<long>& apl);
 
+bigfloat myg1(bigfloat x);
 
 #endif
diff --git a/libsrc/interface.cc b/libsrc/interface.cc
@@ -102,6 +102,7 @@ void Compute_Euler(RR& y)
   bigfloat u, v, a, b, c;
 
   l = RR::precision();
+  RR::SetPrecision(l+20);
 
   x = 1 + static_cast<long>((0.25 * (l - 3)) * (NTL_BITS_PER_LONG * LOG2));
   n = 1 + static_cast<long>(3.591 * x);
@@ -111,7 +112,7 @@ void Compute_Euler(RR& y)
   if (sign(u) > 0)  u=-u;
   a=u;
   v=b=to_bigfloat(1);
-  
+
   for (k = 1; k <= n; k++) {
     mul(b, b, x);
     mul(b, b, x);
@@ -125,6 +126,7 @@ void Compute_Euler(RR& y)
     add(v, v, b);
   }
   div(y, u, v);
+  RR::SetPrecision(l);
 }
 
 long prec() {return RR::precision();}