From d209a2df9d0dce6c50e91d90d1ebf448d12a75cf Mon Sep 17 00:00:00 2001
From: John Cremona <john.cremona@gmail.com>
Date: Mon, 22 Apr 2024 16:31:35 +0100
Subject: [PATCH] improve precision handling for special L-value computations

---
 libsrc/eclib/interface.h          |   3 +-
 libsrc/eclib/periods.h            |  15 +-
 libsrc/interface.cc               |   4 +-
 libsrc/periods.cc                 | 295 ++++++++++++------------------
 progs/out_no_ntl/h1bsd.out        |   6 +-
 progs/out_no_ntl/h1bsdcurisog.out |   6 +-
 progs/out_ntl/h1bsd.out           |   4 +-
 progs/out_ntl/h1bsdcurisog.out    |   4 +-
 8 files changed, 135 insertions(+), 202 deletions(-)

diff --git a/libsrc/eclib/interface.h b/libsrc/eclib/interface.h
index 148c6b01..28e8e92a 100644
--- 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
index e2e75433..17debeab 100644
--- 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
index 92fc1fc3..573c0354 100644
--- 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();}
diff --git a/libsrc/periods.cc b/libsrc/periods.cc
index b9d44f47..193acc00 100644
--- a/libsrc/periods.cc
+++ b/libsrc/periods.cc
@@ -618,6 +618,18 @@ void ldash1::compute(void)
     }
 }
 
+bigfloat ldash1::func1(long n)
+{
+#ifdef MPFP
+  long l = bit_precision();
+  set_precision(l+20);
+#endif
+  bigfloat ans = -G(factor1*to_bigfloat(n));
+#ifdef MPFP
+  set_precision(l);
+#endif
+  return ans;
+}
 
 /////////////////////////////////
 //  functions for lfchi class  //
@@ -861,102 +873,26 @@ Curve newforms::getcurve(long i, int method, bigfloat& rperiod, int verbose)
   return C;
 }
 
-// avoid underflow: log(MINDOUBLE)=-708.396
-
-bigfloat myg0(bigfloat x)
-{
-#ifndef MPFP // Multi-Precision Floating Point
-  if(x>708) return to_bigfloat(0);
-#endif
-  return exp(-x);
-}
-
-//#define TRACEG1
-
-bigfloat myg1(bigfloat x)
-{
-#ifndef MPFP // Multi-Precision Floating Point
-  if(x>708) return to_bigfloat(0);
-#endif
-  if(x<2)
-    {
-      bigfloat ans = -log(x) - Euler();
-      bigfloat p = to_bigfloat(-1);
-#ifdef TRACEG1
-      cout<<"Computing myg1 for x = "<<x<<" using series"<<endl;
-#endif
-      int ok=0;
-      bigfloat term;
-//The following does not work for large x!  (The terms get too big before
-//they get small, and overflow destroys the result.)
-      for (long n=1; (n<5000) && !ok; n++)
-	{
-	  p /=n; p*= -x;
-	  term = p/n;
-	  ans += term;
-	  ok=is_approx_zero(term/ans);
-#ifdef TRACEG1
-	  cout<<"   n="<<n<<": term = " << term << ", ans = "<<ans<<endl;
-#endif
-	}
-#ifdef TRACEG1
-      cout << "returning g1 = " << ans << endl;
-#endif
-      return ans;
-    }
-  //  else  x>2, use continued fraction form from B-G-Z p.478
-#ifdef TRACEG1
-  cout<<"Computing myg1 for x = "<<x<<" using continued fraction"<<endl;
-#endif
-  bigfloat a0, a1, a2, b0, b1=x, b2; a0=0; b0=1;
-  bigfloat ans, newans; ans=0;
-  bigfloat ca;
-  a1=exp(-x);
-  for(long k=2; k<10000; k++)
-    {
-      if (k&1) //then k is odd
-	{ ca = (k-1)/2;
-	  a2=x*a1+ca*a0; a0=a1; a1=a2;
-	  b2=x*b1+ca*b0; b0=b1; b1=b2;
-	}
-      else //k is even
-	{
-	  ca = k/2;
-	  a2=a1+ca*a0; a0=a1; a1=a2;
-	  b2=b1+ca*b0; b0=b1; b1=b2;
-	}
-      newans=a2/b2;
-#ifdef TRACEG1
-      cout<<"   k="<<k<<": approx = " << newans << endl;
-#endif
-      if (is_approx_zero(ans-newans))
-	{
-//	  cout << "returning g1 = " << newans << endl;
-	  return newans;
-	}
-      ans=newans;
-    }
-  cout << "In function g1, continued fraction method, reached end of loop!\n";
-  ans=0;
-  return ans;
-}
-
 bigfloat CG(int r, bigfloat x) // Cohen's G_r(x)
 {
+  // cout<<"CG("<<r<<","<<x<<") = ..."<<flush;
   static const bigfloat one = to_bigfloat(1);
-  bigfloat emx=exp(-x), ans=x, term=x;
-  vector<bigfloat> Av(r+1);  // indexed from 0 to r
-  int j; bigfloat n=one;
-  for(j=0;j<=r;j++) Av[j]=one;
+  bigfloat n=one, emx=exp(-x), ans=x, term=x, y;
+  vector<bigfloat> Av(r+1, one);  // indexed from 0 to r
   while(!is_approx_zero(emx*term*Av[r]))
-    //while(!is_approx_zero(term*Av[r]))
     {
       n++;
       // update A-vector, term and sum
-      for(j=1;j<=r;j++) Av[j]+=(Av[j-1]/n);
-      term*=(x/n);
-      ans+=(Av[r]*term);
+      for(int j=1;j<=r;j++)
+        Av[j]+=(Av[j-1]/n);
+      term*=x;
+      term/=n;
+      y = Av[r]*term;
+      ans+=y;
+      // cout<<"\n"<<y<<"\t"<<ans<<"\t"<<emx*ans;
+      if (is_approx_zero(y/ans)) break;
     }
+  // cout<<"\n"<<emx*ans<<" (using "<<n<<" terms)"<<endl;
   return emx*ans;
 }
 
@@ -965,26 +901,27 @@ bigfloat Glarge(int r, bigfloat x) // Cohen's Gamma_r(x) for large x
   static const bigfloat zero = to_bigfloat(0);
   static const bigfloat one = to_bigfloat(1);
   static const bigfloat two = to_bigfloat(2);
-  bigfloat emx=exp(-x), ans=zero, mxinv=-one/x;
-  bigfloat term = mxinv;
-  vector<bigfloat> Cv(r+1);  // indexed from 0 to r
-  int j; bigfloat n=zero;
+  bigfloat emx2=2*exp(-x), ans=zero, mxinv=-one/x;
+  bigfloat y, term = mxinv, n=zero;
+  if (is_approx_zero(abs(emx2*term)))
+    return zero;
+  // cout<<"Glarge("<<r<<","<<x<<") = ..."<<flush;
+  vector<bigfloat> Cv(r+1, zero);  // indexed from 0 to r
   Cv[0]=one;
-  for(j=1;j<=r;j++) Cv[j]=zero;
-  //cout<<"In Glarge with r="<<r<<", x="<<x<<": ";
-  //cout<<"emx*term="<<emx*term<<endl;
-  while((n<1000) && !is_approx_zero(abs(emx*term)))
-    //while(!is_approx_zero(abs(term)))
+  while((n<1000) && !is_approx_zero(abs(emx2*term)))
     {
       n++;
       // update C-vector, term and sum
-      for(j=r;j>0;j--) Cv[j]+=(Cv[j-1]/n);
-      term*=(n*mxinv);
-      ans+=(Cv[r]*term);
-      //cout<<"term="<<term<<"; ans="<<ans<<endl;
+      for(int j=r;j>0;j--) Cv[j]+=(Cv[j-1]/n);
+      term*=n;
+      term*=mxinv;
+      y = Cv[r]*term;
+      ans+=y;
+      // cout<<"\n"<<y<<"\t"<<ans<<"\t"<<emx2*ans;
+      if (is_approx_zero(y/ans)) break;
     }
-  //cout<<"Glarge returns "<<two*emx*ans<<" after using n="<<n<<" terms"<<endl;
-  return two*emx*ans;
+  // cout<<"\n"<<emx2*ans<<" (using "<<n<<" terms)"<<endl;
+  return emx2*ans;
 }
 
 bigfloat Q(int r, bigfloat x)  // Q_r(x) polynomial from AMEC p.44
@@ -1020,11 +957,13 @@ bigfloat Q(int r, bigfloat x)  // Q_r(x) polynomial from AMEC p.44
   case 4: return const1 +x*(-const2+x*(const3+x*x*twenty4th));
   }
 }
+
 bigfloat P(int r, bigfloat x)  // P_r(x) polynomial from AMEC p.44
 {
   static bigfloat gamma=Euler();
   return Q(r,x-gamma);
 }
+
 bigfloat Gsmall(int r, bigfloat x) // Cohen's Gamma_r(x) for small x
 {
   bigfloat a=P(r,-log(x)), b = CG(r,x);
@@ -1034,94 +973,98 @@ bigfloat Gsmall(int r, bigfloat x) // Cohen's Gamma_r(x) for small x
 bigfloat G(int r, bigfloat x)  // G_r(x)
 {
 #ifndef MPFP // Multi-Precision Floating Point
-    static const  bigfloat x0 = to_bigfloat(14);
+  static const  bigfloat x0 = to_bigfloat(14);
 #else
-    // log(2) = 0.69314718
-    static const  bigfloat x0 = to_bigfloat(0.69314718*bit_precision());
+  // log(2) = 0.693147180599..
+  // this value only determines whether to use Gsmall or Glarge:
+  static const  bigfloat x0 = to_bigfloat(0.693147*bit_precision());
 #endif
-    //  cout<<"switch point = "<<x0<<endl;
-    //  cout<<"x="<<x<<endl;
-  if(x<x0)
-    {
-      return Gsmall(r,x);
-    }
-  else
-    {
-      return Glarge(r,x);
-    }
+  //  cout<<"switch point = "<<x0<<endl;
+  //  cout<<"x="<<x<<endl;
+
+  return x<x0? Gsmall(r,x) : Glarge(r,x);
 }
 
 bigfloat ldash1::G(bigfloat x)  // G_r(x)
 {
-  switch(r) {
-  case 0: return myg0(x);
-  case 1: return myg1(x);
-  default: return ::G(r,x);
-  }
+#ifndef MPFP // if without Multi-Precision Floating Point, avoid underflow: log(MINDOUBLE)=-708.396
+  if((r<2) && (x>708))
+    return to_bigfloat(0);
+  if(r==1)
+    return myg1(x);
+#endif
+  return r==0? exp(-x) : ::G(r,x);
 }
 
-#if(0) // myg2 and myg3 were inaccurate and now replaced by general
-       // G(r,x) -- which also works for r>3!
-bigfloat myg2(bigfloat x)
+// myg1(x) is the same as G_1(x) = G(1,x)
+
+//#define TRACEG1
+
+bigfloat myg1(bigfloat x)
 {
-  static bigfloat zero=to_bigfloat(0);
-  static bigfloat one=to_bigfloat(1);
-  static bigfloat two=to_bigfloat(2);
-  static bigfloat twelve=to_bigfloat(12);
   static bigfloat gamma=Euler();
-  if (x>20) return zero;
-//#define TRACEG2
-  bigfloat ans = -log(x) - gamma;
-#ifdef TRACEG2
-  cout<<"Computing myg2 for x = "<<x<<endl;
+#ifndef MPFP // Multi-Precision Floating Point
+  if(x>708) return to_bigfloat(0);
 #endif
-  ans=ans*ans/two + PI*PI/twelve;
-  bigfloat p = one;
-  int ok=0;
-  bigfloat term;
-  for (long n=1; (n<500) && !ok; n++)
+  if(x<2)
     {
-      p /=n;  p*= -x;
-      term = (p/n)/n;
-      ans += term;
-      ok=is_approx_zero(term/ans);
-#ifdef TRACEG2
-      // cout<<"\tn="<<n<<" \tans = "<<ans<<endl;
+      bigfloat ans = -log(x) - gamma;
+      bigfloat p = to_bigfloat(-1);
+#ifdef TRACEG1
+      cout<<"Computing myg1 for x = "<<x<<" using series"<<endl;
 #endif
-    }
-#ifdef TRACEG2
-  cout<<"...returning ans = "<<ans<<endl;
+      int ok=0;
+      bigfloat term;
+//The following does not work for large x!  (The terms get too big before
+//they get small, and overflow destroys the result.)
+      for (long n=1; (n<5000) && !ok; n++)
+	{
+	  p /=n; p*= -x;
+	  term = p/n;
+	  ans += term;
+	  ok=is_approx_zero(term/ans);
+#ifdef TRACEG1
+	  cout<<"   n="<<n<<": term = " << term << ", ans = "<<ans<<endl;
 #endif
-   return ans;
-}
-
-bigfloat myg3(bigfloat x)
-{
-  static bigfloat zero=to_bigfloat(0);
-  static bigfloat one=to_bigfloat(1);
-  static bigfloat three=to_bigfloat(3);
-  static bigfloat twelve=to_bigfloat(12);
-#ifdef MPFP // Multi-Precision Floating Point
-  static const bigint nz3=atoI("2756915843737840912679655856873838262816890298077497105924627168570887325226967786076589016002130138897164");
-  bigfloat zeta3; MakeRR(zeta3,nz3,-350);
-#else
- const bigfloat zeta3 = 1.20205690315959428539973816151144999076498629;
+	}
+#ifdef TRACEG1
+      cout << "returning g1 = " << ans << endl;
+#endif
+      return ans;
+    }
+  //  else  x>2, use continued fraction form from B-G-Z p.478
+#ifdef TRACEG1
+  cout<<"Computing myg1 for x = "<<x<<" using continued fraction"<<endl;
 #endif
-  if (x>20) return zero;
-  bigfloat ans = -log(x) - Euler();
-//cout<<"Computing myg3 for x = "<<x<<endl;
-  ans = (PI*PI + 2*ans*ans) * ans/twelve - zeta3/three;
-  bigfloat p = -one;
-  bigfloat term;
-  int ok = 0;
-  for (long n=1; (n<500) && !ok; n++)
+  bigfloat a0, a1, a2, b0, b1=x, b2; a0=0; b0=1;
+  bigfloat ans, newans; ans=0;
+  bigfloat ca;
+  a1=exp(-x);
+  for(long k=2; k<10000; k++)
     {
-      p *= -x/n;
-      term = ((p/n)/n)/n;
-      ans += term;
-      ok=is_approx_zero(term/ans);
-      //cout<<"   N="<<n<<" ans = "<<ans<<endl;
+      if (k&1) //then k is odd
+	{ ca = (k-1)/2;
+	  a2=x*a1+ca*a0; a0=a1; a1=a2;
+	  b2=x*b1+ca*b0; b0=b1; b1=b2;
+	}
+      else //k is even
+	{
+	  ca = k/2;
+	  a2=a1+ca*a0; a0=a1; a1=a2;
+	  b2=b1+ca*b0; b0=b1; b1=b2;
+	}
+      newans=a2/b2;
+#ifdef TRACEG1
+      cout<<"   k="<<k<<": approx = " << newans << endl;
+#endif
+      if (is_approx_zero(ans-newans))
+	{
+//	  cout << "returning g1 = " << newans << endl;
+	  return newans;
+	}
+      ans=newans;
     }
-   return ans;
+  cout << "In function g1, continued fraction method, reached end of loop!\n";
+  ans=0;
+  return ans;
 }
-#endif
diff --git a/progs/out_no_ntl/h1bsd.out b/progs/out_no_ntl/h1bsd.out
index e9fc2fec..966a1b0f 100644
--- a/progs/out_no_ntl/h1bsd.out
+++ b/progs/out_no_ntl/h1bsd.out
@@ -3,14 +3,14 @@ Enter first and last N:
 50 b	Rank = 0	L^(r)(f,1)/r! = 0.95681122658763 	L(f,1)/period = 1/5
 51 a	Rank = 0	L^(r)(f,1)/r! = 0.86005901612714 	L(f,1)/period = 1/3
 52 a	Rank = 0	L^(r)(f,1)/r! = 0.84548320864565 	L(f,1)/period = 1/2
-53 a	Rank = 1	L^(r)(f,1)/r! = 0.43586382412977 	L(f,1)/period = 0
+53 a	Rank = 1	L^(r)(f,1)/r! = 0.43586382417039 	L(f,1)/period = 0
 54 a	Rank = 0	L^(r)(f,1)/r! = 0.70157482532401 	L(f,1)/period = 1/3
 54 b	Rank = 0	L^(r)(f,1)/r! = 1.030521849701 	L(f,1)/period = 1/3
 55 a	Rank = 0	L^(r)(f,1)/r! = 1.0286698561782 	L(f,1)/period = 1/4
 56 b	Rank = 0	L^(r)(f,1)/r! = 0.94801911561387 	L(f,1)/period = 1/2
 56 a	Rank = 0	L^(r)(f,1)/r! = 0.87454831418834 	L(f,1)/period = 1/4
-57 a	Rank = 1	L^(r)(f,1)/r! = 0.41749163961642 	L(f,1)/period = 0
+57 a	Rank = 1	L^(r)(f,1)/r! = 0.41749163967662 	L(f,1)/period = 0
 57 b	Rank = 0	L^(r)(f,1)/r! = 0.58650959165657 	L(f,1)/period = 2/5
 57 c	Rank = 0	L^(r)(f,1)/r! = 1.0853022722677 	L(f,1)/period = 1/4
-58 a	Rank = 1	L^(r)(f,1)/r! = 0.46370416461717 	L(f,1)/period = 0
+58 a	Rank = 1	L^(r)(f,1)/r! = 0.46370416478316 	L(f,1)/period = 0
 58 b	Rank = 0	L^(r)(f,1)/r! = 1.0332075093651 	L(f,1)/period = 2/5
diff --git a/progs/out_no_ntl/h1bsdcurisog.out b/progs/out_no_ntl/h1bsdcurisog.out
index 8e437ea3..09105e80 100644
--- a/progs/out_no_ntl/h1bsdcurisog.out
+++ b/progs/out_no_ntl/h1bsdcurisog.out
@@ -1,11 +1,11 @@
     11 a   1 [0,-1,1,-10,-20]	0	5	5	1.2692093042796	0.25384186085591	1	1
     11 a   2 [0,-1,1,-7820,-263580]	0	1	1	0.25384186085591	0.25384186085591	1	1
     11 a   3 [0,-1,1,0,0]	0	5	1	6.3460465213978	0.25384186085591	1	1
-    37 a   1 [0,0,1,-1,0]	1	1	1	5.9869172924639	0.30599977380009	0.051111408239969	0.999999999889
+    37 a   1 [0,0,1,-1,0]	1	1	1	5.9869172924639	0.30599977382942	0.051111408239969	0.99999999998486
     37 b   1 [0,1,1,-23,-50]	0	3	3	2.1770431858085	0.72568106193615	1	1
     37 b   2 [0,1,1,-3,1]	0	3	1	6.5311295574254	0.72568106193615	1	1
     37 b   3 [0,1,1,-1873,-31833]	0	1	1	0.72568106193659	0.72568106193615	1	1
-    57 a   1 [0,-1,1,-2,2]	1	1	2	5.5555045213725	0.41749163961642	0.037574592736823	0.99999999985658
+    57 a   1 [0,-1,1,-2,2]	1	1	2	5.5555045213725	0.41749163967662	0.037574592736823	1.0000000000008
     57 b   1 [0,1,1,20,-32]	0	5	10	1.4662739791414	0.58650959165657	1	1
     57 b   2 [0,1,1,-4390,-113432]	0	1	2	0.29325479582828	0.58650959165657	1	1
     57 c   1 [1,0,1,-7,5]	0	4	4	4.3412090890707	1.0853022722677	1	1
@@ -28,4 +28,4 @@
     90 c   6 [1,-1,1,-3002,63929]	0	12	144	1.3375959945886	1.3375959945886	1	1
     90 c   7 [1,-1,1,-48002,4059929]	0	6	36	1.3375959945892	1.3375959945886	1	1
     90 c   8 [1,-1,1,-4082,14681]	0	6	72	0.66879799729432	1.3375959945886	1	1
-   389 a   1 [0,1,1,-2,0]	2	1	1	4.9804251217101	0.75931650084309	0.15246017794314	1.0000000007305
+   389 a   1 [0,1,1,-2,0]	2	1	1	4.9804251217101	0.75931650083449	0.15246017794314	1.0000000007191
diff --git a/progs/out_ntl/h1bsd.out b/progs/out_ntl/h1bsd.out
index 27d6c763..6f62a1e1 100644
--- a/progs/out_ntl/h1bsd.out
+++ b/progs/out_ntl/h1bsd.out
@@ -3,7 +3,7 @@ Enter first and last N:
 50 b	Rank = 0	L^(r)(f,1)/r! = 0.95681122658762796 	L(f,1)/period = 1/5
 51 a	Rank = 0	L^(r)(f,1)/r! = 0.86005901612713837 	L(f,1)/period = 1/3
 52 a	Rank = 0	L^(r)(f,1)/r! = 0.84548320864564628 	L(f,1)/period = 1/2
-53 a	Rank = 1	L^(r)(f,1)/r! = 0.43586382417785719 	L(f,1)/period = 0
+53 a	Rank = 1	L^(r)(f,1)/r! = 0.43586382417785716 	L(f,1)/period = 0
 54 a	Rank = 0	L^(r)(f,1)/r! = 0.70157482532401025 	L(f,1)/period = 1/3
 54 b	Rank = 0	L^(r)(f,1)/r! = 1.0305218497010025 	L(f,1)/period = 1/3
 55 a	Rank = 0	L^(r)(f,1)/r! = 1.0286698561781891 	L(f,1)/period = 1/4
@@ -12,5 +12,5 @@ Enter first and last N:
 57 a	Rank = 1	L^(r)(f,1)/r! = 0.41749163967630964 	L(f,1)/period = 0
 57 b	Rank = 0	L^(r)(f,1)/r! = 0.58650959165656799 	L(f,1)/period = 2/5
 57 c	Rank = 0	L^(r)(f,1)/r! = 1.0853022722676732 	L(f,1)/period = 1/4
-58 a	Rank = 1	L^(r)(f,1)/r! = 0.46370416478825061 	L(f,1)/period = 0
+58 a	Rank = 1	L^(r)(f,1)/r! = 0.46370416478825058 	L(f,1)/period = 0
 58 b	Rank = 0	L^(r)(f,1)/r! = 1.033207509365107 	L(f,1)/period = 2/5
diff --git a/progs/out_ntl/h1bsdcurisog.out b/progs/out_ntl/h1bsdcurisog.out
index a71325f8..ac310a7f 100644
--- a/progs/out_ntl/h1bsdcurisog.out
+++ b/progs/out_ntl/h1bsdcurisog.out
@@ -1,7 +1,7 @@
     11 a   1 [0,-1,1,-10,-20]	0	5	5	1.2692093042795534216887946168	0.25384186085591068433775892335	1	1
     11 a   2 [0,-1,1,-7820,-263580]	0	1	1	0.25384186085591068433775892335	0.25384186085591068433775892335	1	1
     11 a   3 [0,-1,1,0,0]	0	5	1	6.3460465213977671084439730838	0.25384186085591068433775892335	1	1
-    37 a   1 [0,0,1,-1,0]	1	1	1	5.9869172924639192596640199589	0.30599977383405230182048368337	0.051111408239968840235886099757	1.0000000000000000000000000001
+    37 a   1 [0,0,1,-1,0]	1	1	1	5.9869172924639192596640199589	0.30599977383405230182048368332	0.051111408239968840235886099757	1
     37 b   1 [0,1,1,-23,-50]	0	3	3	2.1770431858084583470086166231	0.72568106193615278233620554103	1	1
     37 b   2 [0,1,1,-3,1]	0	3	1	6.5311295574253750410258498692	0.72568106193615278233620554103	1	1
     37 b   3 [0,1,1,-1873,-31833]	0	1	1	0.72568106193615278233620553794	0.72568106193615278233620554103	1	1
@@ -28,4 +28,4 @@
     90 c   6 [1,-1,1,-3002,63929]	0	12	144	1.3375959945886485057653424429	1.3375959945886485057653424429	1	1
     90 c   7 [1,-1,1,-48002,4059929]	0	6	36	1.3375959945886485057653474741	1.3375959945886485057653424429	1	1
     90 c   8 [1,-1,1,-4082,14681]	0	6	72	0.66879799729432425288267122147	1.3375959945886485057653424429	1	1
-   389 a   1 [0,1,1,-2,0]	2	1	1	4.9804251217101101506427155839	0.75931650028842677023019260774	0.15246017794314375162432475705	0.9999999999999999999999999998
+   389 a   1 [0,1,1,-2,0]	2	1	1	4.9804251217101101506427155839	0.7593165002884267702301926079	0.15246017794314375162432475705	1