Merge pull request #35 from JohnCremona/prec

reworked precision to use bit precision everywhere except for output

Author: JohnCremona

.travis.yml

@@ -6,10 +6,9 @@ os:
 osx_image: xcode6.4
 dist: trusty
 
-matrix:
-  include:
-    - env: CONFIGURE_OPTS="--disable-mpfp"
-      os: linux
+env:
+    - CONFIGURE_OPTS=""
+    - CONFIGURE_OPTS="--disable-mpfp"
 
 install:
     - |

doc/mwrank/mwrank.options

@@ -41,10 +41,10 @@ List of command line options, their meaning and default values (set in options.h
 			options (controlled by -v); for minimal useful
 			output use "-q -v 0 -o"
 
--p n	precision (#dp)
+-p n	precision (#bits)
 			Only relevant for the multiprecision
 			floating point version.  Range 1..large,
-			default 15.  This is decimal precision.
+			default 50.  This is bit precision.
 			Stringly recommended to increase for curves
 			with 2-torsion where a second descent is used
 			since then the reduction applied to second

libsrc/compproc.cc

@@ -74,6 +74,15 @@ bigcomplex normalize(bigcomplex& w1, bigcomplex& w2)
     w1=w1-w2*round(tau.real());
     tau=w1/w2;
    }
+   if (tau.real()==-0.5)
+     {
+       w1+=w2; tau=w1/w2;
+     }
+   else
+     if (abs(tau)==1 && tau.real()<0)
+       {
+         w3=-w1; w1=w2; w2=w3; tau=w1/w2;
+       }
    return tau;
 }
 

libsrc/eclib/interface.h

@@ -230,14 +230,33 @@ inline CC &operator /=(CC &a, const CC &b)
 
 #endif
 
+// NB Internal precision is always bit-precision.  We used to use
+// decimal precision in the user interface with a conversion factor of
+// 3.33 (approx log(10)/log(2)).  NTL has a default internal bit
+// precision of 150 which can be changed using RR::SetPrecision(), and
+// RR:SetOutputPrecision(d) sets the output to d decimal places
+// (default 10).  See www.shoup.net/ntl/doc/tour-ex6.html and
+// www.shoup.net/ntl/doc/RR.cpp.html.
+
+// Set internal precision to n bits and output precision to (0.3*n)-1 decimal places
 inline void set_precision(long n)
-  {RR::SetPrecision(long(n*3.33));RR::SetOutputPrecision(n);}
+{RR::SetPrecision(n); RR::SetOutputPrecision(long(0.3*n)-1);}
+
+// Mostly for backward compatibility (used in saturate.cc) or for
+// temporarily changing internal precision when no output is relevant:
 inline void set_bit_precision(long n)
   {RR::SetPrecision(n);}
+
+// Prompt user for precision
 inline void set_precision(const string prompt)
   {long n; cerr<<prompt<<": "; cin>>n; set_precision(n);}
+
+// read current precision converted to decimal (approximately)
 inline long decimal_precision() {return long(RR::precision()*0.3);}
+
+// read current bit precision
 inline long bit_precision() {return RR::precision();}
+
 inline int is_approx_zero(const bigcomplex& z)
   {return is_approx_zero(z.real())&&is_approx_zero(z.imag());}
 inline RR to_bigfloat(const int& n) {return to_RR(n);}
@@ -265,9 +284,12 @@ inline CC pow(const CC& a, const RR& e)  {return exp(e*log(a));}
 typedef double bigfloat;
 
 inline long decimal_precision() {return 15;}
+inline long bit_precision() {return 53;}
 inline int is_zero(double x) {return fabs(x)<1e-15;}
 inline int is_approx_zero(double x) {return fabs(x)<1e-10;}
-inline void set_precision(long n) {cout.precision(n);}
+
+// We cannot set internal bit precision in this mode, so we just set the output decimal precision
+inline void set_precision(long n) {cout.precision(min(15,long(0.3*n)));}
 inline void set_precision(const string prompt)  {cout.precision(15);}
 #define Pi()    3.1415926535897932384626433832795028841
 #define Euler() (0.57721566490153286060651209008240243104)

libsrc/eclib/options.h

@@ -31,7 +31,7 @@ using namespace std;
 
 #define DEFAULT_QUIET 0
 #define DEFAULT_VERBOSE 1
-#define DEFAULT_PRECISION 15
+#define DEFAULT_PRECISION 50
 #define DEFAULT_HLIMQ 10
 #define DEFAULT_NAUX 15
 #define DEFAULT_HLIMC 0
@@ -46,7 +46,7 @@ class mrank_options {
 private:
   int quiet;            // 0/1, controls header output
   int verbose;          // 0-3, controls output verbosity
-  long precision;       // 1-\infty, controls decimal precision
+  long precision;       // 1-\infty, controls floating point precision (in bits)
   long hlimq;           // 1-20, height limit for quartic search
   long naux;            // number of primes used in syzygy sieve
   long hlimc;           // 1-15, height limit for curve search
@@ -125,7 +125,7 @@ class mrank_options {
       cerr << "-q\t""quiet""\t\tturns OFF banner display\n";
       cerr << "-v n\t""verbosity""\tsets verbosity to n (default="<<DEFAULT_VERBOSE<<")\n";
       cerr << "-o\t""PARI/GP output""\tturns ON extra PARI/GP short output (default is OFF)\n";
-      cerr << "-p n\t""precision""\tsets precision to n decimals (default="<<DEFAULT_PRECISION<<")\n";
+      cerr << "-p n\t""precision""\tsets real precision to n bits (default="<<DEFAULT_PRECISION<<")\n";
       cerr << "-b n\t""quartic bound""\tbound on quartic point search (default="<<DEFAULT_HLIMQ<<")\n";
       cerr << "-x n\t""n aux""\t\tnumber of aux primes used for sieving (default="<<DEFAULT_NAUX<<")\n";
       cerr << "-l\t""list""\t\tturns ON listing of points (default ON unless v=0)\n";
@@ -163,7 +163,7 @@ class mrank_options {
       if(quiet)cerr<<"ON"; else cerr<<"OFF"; cerr<<"\n";
       cerr << "PARI/GP output "; 
       if(output_pari)cerr<<"ON"; else cerr<<"OFF"; cerr<<"\n";
-      cerr << "Precision = " << precision << " decimal places (only relevant for multiprecision version)\n";
+      cerr << "Precision = " << precision << " bits (only relevant for multiprecision version)\n";
       cerr << "Verbosity level = " << verbose << "\n";
       cerr << "Limit on height for point search on quartics: "<<hlimq<<"\n";
       cerr << "Number of auxiliary primes for syzygy sieving: "<<naux<<"\n";

libsrc/heights.cc

@@ -156,6 +156,8 @@ bigfloat realheight(const bigfloat& x, const Curvedata* E)
   t=2*abs(b4); if(t>H) H=t;
   t=2*abs(b6); if(t>H) H=t;
   t=abs(b8);   if(t>H) H=t;
+  // NB We use decimal precision here since the formula for nlim (from
+  // Silverman) is given in terms of decimal places required.
   long precision = decimal_precision();
 #ifdef DEBUG_HEIGHT
   cout<<"decimal precision = "<<precision<<endl;

libsrc/periods.cc

@@ -371,14 +371,14 @@ periods_via_lfchi::periods_via_lfchi (const level* iN, const newform* f)
   rootmod =  sqrt(to_bigfloat(N));  
   factor1 =  exp(-(TWOPI)/ (to_bigfloat(chi1.modulus()) * rootmod));
   factor2 =  exp(-(TWOPI)/ (to_bigfloat(chi2.modulus()) * rootmod));
-  long dp = decimal_precision();
+  long bp = bit_precision();
  // MUST keep brackets around TWOPI!:
-  bigfloat rootlimit1=(dp*LOG10-log(1-factor1))*rootmod/(TWOPI) ;
-  bigfloat rootlimit2=(dp*LOG10-log(1-factor1))*rootmod/(TWOPI) ; 
+  bigfloat rootlimit1=(bp-log(1-factor1))*rootmod/(TWOPI) ;
+  bigfloat rootlimit2=(bp-log(1-factor1))*rootmod/(TWOPI) ; 
   Iasb(limit1,rootlimit1);
   Iasb(limit2,rootlimit2);
 #ifdef TRACE_USE 
-  cout << "Decimal precision = "<<dp<<endl;
+  cout << "Bit precision = "<<bp<<endl;
   cout << "Basic limits on n in sums = " << limit << endl;
 #endif
   limit1=chi1.modulus()*limit1;  
@@ -390,7 +390,7 @@ periods_via_lfchi::periods_via_lfchi (const level* iN, const newform* f)
 #endif
 /*
   if(primelist[primelist.size()-1]<limit)
-    cout<<"\nLARGEST PRIME "<<primelist[primelist.size()-1]<<" IS BELOW LIMIT "<< limit <<" required for decimal precision "<<dp<<endl;
+    cout<<"\nLARGEST PRIME "<<primelist[primelist.size()-1]<<" IS BELOW LIMIT "<< limit <<" required for bit precision "<<bp<<endl;
 */
   rootlimit=sqrt(to_bigfloat(limit));
   an_cache.resize(I2long(Ifloor(rootlimit+1)),0);
@@ -465,9 +465,9 @@ void periods_direct::compute(void)
   b = posmod(b,d);
   c = posmod(c,d);
   factor2 = factor1 * drecip;
-  long dp = decimal_precision();
-  Iasb(limit1,(-dp*LOG10-log((1-exp(factor1))/3))/(factor1)) ;
-  Iasb(limit2,(-dp*LOG10-log((1-exp(factor2))/3))/(factor2)) ;
+  long bp = bit_precision();
+  Iasb(limit1,(-bp-log((1-exp(factor1))/3))/(factor1)) ;
+  Iasb(limit2,(-bp-log((1-exp(factor2))/3))/(factor2)) ;
 
   limit = limit2;
   rootlimit=sqrt(to_bigfloat(limit));
@@ -479,7 +479,7 @@ void periods_direct::compute(void)
 #endif
 /*
   if(primelist[primelist.size()-1]<limit)
-    cout<<"\nLARGEST PRIME "<<primelist[primelist.size()-1]<<" IS BELOW LIMIT "<< limit <<" required for decimal precision "<<dp<<endl;
+    cout<<"\nLARGEST PRIME "<<primelist[primelist.size()-1]<<" IS BELOW LIMIT "<< limit <<" required for bit precision "<<bp<<endl;
 */
 #ifdef TRACE_CACHE
   cout << "Initial an_cache = "<<an_cache<<endl;
@@ -555,8 +555,8 @@ void part_period::compute(const bigcomplex& z0)
 
 void part_period::compute()
 {
-  long dp = decimal_precision();
-  Iasb(limit,dp*LOG10/y0);
+  long bp = bit_precision();
+  Iasb(limit,bp/y0);
   limit1=limit2=limit;
   rootlimit=sqrt(to_bigfloat(limit));
   an_cache.resize(I2long(Ifloor(rootlimit+1)),0);
@@ -595,7 +595,7 @@ void ldash1::init(const level* iN, const vector<long>& f_aplist, long f_sfe, con
   rootmod=sqrt(to_bigfloat(N));
   factor1 = (TWOPI)/rootmod;
   long maxp = prime_number(nap);
-  limit  =  I2long(Ifloor((15+decimal_precision())*log(10)/factor1));
+  limit  =  I2long(Ifloor((30+bit_precision())/factor1));
   if(limit>maxp) limit=maxp;
   limit1 = limit;
   rootlimit=sqrt(to_bigfloat(limit));
@@ -636,8 +636,8 @@ lfchi::lfchi (const level* iN, const newform* f)
 {
   initaplist(iN, f->aplist);
   rootmod =  sqrt(to_bigfloat(N));
-  long dp = decimal_precision();
-  Iasb(limit0,dp*LOG10*rootmod/(TWOPI)) ; // brackets essential
+  long bp = bit_precision();
+  Iasb(limit0,bp*rootmod/(TWOPI)) ; // brackets essential
 }
 
 void lfchi::compute(long ell)
@@ -843,7 +843,7 @@ Curve newforms::getcurve(long i, int method, bigfloat& rperiod, int verbose)
   if((abs((wR-wRC)/wRC)>0.0001))
     {
       cout<<"Real period of constructed curve does not match that"
-	  <<" of the newform (using decimal precision "<<decimal_precision()<<")"<<endl;
+	  <<" of the newform (using bit precision "<<bit_precision()<<")"<<endl;
       cout<<"Real period of C: "<<real(wRC)<<endl;
       cout<<"Real period of f: "<<real(wR)<<endl;
       cout<<"Ratio = "<<real(wR)/real(wRC)<<endl;
@@ -852,7 +852,7 @@ Curve newforms::getcurve(long i, int method, bigfloat& rperiod, int verbose)
   if((abs((wRI-wRIC)/wRIC)>0.0001))
     {
       cout<<"Second period of constructed curve does not match that"
-	  <<" of the newform (using decimal precision "<<decimal_precision()<<")"<<endl;
+	  <<" of the newform (using bit precision "<< bit_precision()<<")"<<endl;
       cout<<"Imag part of second period of C: "<<real(wRIC)<<endl;
       cout<<"Imag part of second period of f: "<<real(wRI)<<endl;
       cout<<"Ratio of imaginary parts = "<<real(wRI)/real(wRIC)<<endl;
@@ -1034,7 +1034,8 @@ bigfloat G(int r, bigfloat x)  // G_r(x)
 #ifndef MPFP // Multi-Precision Floating Point
     static const  bigfloat x0 = to_bigfloat(14);
 #else
-    static const  bigfloat x0 = to_bigfloat(log(10)*decimal_precision());
+    // log(2) = 0.69314718
+    static const  bigfloat x0 = to_bigfloat(0.69314718*bit_precision());
 #endif
     //  cout<<"switch point = "<<x0<<endl;
     //  cout<<"x="<<x<<endl;
@@ -1058,7 +1059,7 @@ bigfloat ldash1::G(bigfloat x)  // G_r(x)
 }
 
 #if(0) // myg2 and myg3 were inaccurate and now replaced by general
-       // G(r,x) -- whcih also works for r>3!
+       // G(r,x) -- which also works for r>3!
 bigfloat myg2(bigfloat x)
 {
   static bigfloat zero=to_bigfloat(0);

libsrc/saturate.cc

@@ -220,7 +220,7 @@ int saturator::enlarge()
         {
           cout<<"...but it isn't! (this may be due to insufficient precision:";
 #ifdef NTL_ALL
-          cout << " decimal precision " <<prec<<" was used)"<<endl;
+          cout << " bit precision " <<prec<<" was used)"<<endl;
 #else
           cout << " standard double precision was used)"<<endl;
 #endif

progs/h1bsd.cc

@@ -35,7 +35,7 @@
 
 int main(void)
 {
-  set_precision(15);
+  set_precision(50);
  int limit,n=1; 
 #ifdef AUTOLOOP
  cout<<"Enter first and last N: ";cin>>n>>limit; 

progs/h1bsdcurisog.cc

@@ -60,7 +60,7 @@
 
 int main(void)
 {
-  set_precision(10);
+  set_precision(35);
 
  long limit, n=1, hlim1=10, hlim2=15; 
  bigint nn;