Attempt to use the Dekker high-precision floating point algorithm in the

ext/misc/totype.c extension.

FossilOrigin-Name: dd41db3b106345c510f1f64aebcedf162701e9e21bd157441331b454ce11b353

diff --git a/ext/misc/totype.c b/ext/misc/totype.c
index 31c497a5675a63b73af256480f90baaac63f2b7f..0f6c3265b88521d5779cd4b224150706d41d6769 100644 (file)
--- a/ext/misc/totype.c
+++ b/ext/misc/totype.c
@@ -180,6 +180,62 @@ static int totypeAtoi64(const char *zNum, sqlite3_int64 *pNum, int length){
   }
 }
 
+#ifndef LONGDOUBLE_TYPE
+# define LONGDOUBLE_TYPE long double
+#endif
+typedef sqlite3_uint64 u64;
+#define LARGEST_UINT64 (0xffffffff|(((sqlite3_uint64)0xffffffff)<<32))
+#define EXP754 (((u64)0x7ff)<<52)
+#define MAN754 ((((u64)1)<<52)-1)
+#define IsNaN(X) (((X)&EXP754)==EXP754 && ((X)&MAN754)!=0)
+
+/*
+** Return true if the floating point value is Not a Number (NaN).
+*/
+int totypeIsNaN(double x){
+  int rc;   /* The value return */
+  u64 y;
+  memcpy(&y,&x,sizeof(y));
+  rc = IsNaN(y);
+  return rc;
+}
+
+/* Double-Double multiplication.  (x[0],x[1]) *= (y,yy)
+**
+** Reference:
+**   T. J. Dekker, "A Floating-Point Technique for Extending the
+**   Available Precision".  1971-07-26.
+*/
+static void dekkerMul2(volatile double *x, double y, double yy){
+  /*
+  ** The "volatile" keywords on parameter x[] and on local variables
+  ** below are needed force intermediate results to be truncated to
+  ** binary64 rather than be carried around in an extended-precision
+  ** format.  The truncation is necessary for the Dekker algorithm to
+  ** work.  Intel x86 floating point might omit the truncation without
+  ** the use of volatile. 
+  */
+  volatile double tx, ty, p, q, c, cc;
+  double hx, hy;
+  u64 m;
+  memcpy(&m, (void*)&x[0], 8);
+  m &= 0xfffffffffc000000LL;
+  memcpy(&hx, &m, 8);
+  tx = x[0] - hx;
+  memcpy(&m, &y, 8);
+  m &= 0xfffffffffc000000LL;
+  memcpy(&hy, &m, 8);
+  ty = y - hy;
+  p = hx*hy;
+  q = hx*ty + tx*hy;
+  c = p+q;
+  cc = p - c + q + tx*ty;
+  cc = x[0]*yy + x[1]*y + cc;
+  x[0] = c + cc;
+  x[1] = c - x[0];
+  x[1] += cc;
+}
+
 /*
 ** The string z[] is an text representation of a real number.
 ** Convert this string to a double and write it into *pResult.
@@ -213,6 +269,7 @@ static int totypeAtoF(const char *z, double *pResult, int length){
   double result;
   int nDigits = 0;
   int nonNum = 0;
+  int bUseLongDouble = sizeof(LONGDOUBLE_TYPE)>8;
 
   *pResult = 0.0;   /* Default return value, in case of an error */
 
@@ -282,70 +339,90 @@ static int totypeAtoF(const char *z, double *pResult, int length){
   }
 
 totype_atof_calc:
+  /* Zero is a special case */
+  if( s==0 ){
+    *pResult = sign<0 ? -0.0 : +0.0;
+    goto atof_return;
+  }
+
   /* adjust exponent by d, and update sign */
   e = (e*esign) + d;
-  if( e<0 ) {
-    esign = -1;
-    e *= -1;
-  } else {
-    esign = 1;
+
+  /* Try to adjust the exponent to make it smaller */
+  while( e>0 && s<(LARGEST_UINT64/10) ){
+    s *= 10;
+    e--;
+  }
+  while( e<0 && (s%10)==0 ){
+    s /= 10;
+    e++;
   }
 
-  /* if 0 significand */
-  if( !s ) {
-    /* In the IEEE 754 standard, zero is signed.
-    ** Add the sign if we've seen at least one digit */
-    result = (sign<0 && nDigits) ? -(double)0 : (double)0;
-  } else {
-    /* attempt to reduce exponent */
-    if( esign>0 ){
-      while( s<(LARGEST_INT64/10) && e>0 ) e--,s*=10;
+  if( e==0 ){
+    *pResult = s;
+  }else if( bUseLongDouble ){
+    LONGDOUBLE_TYPE r = (LONGDOUBLE_TYPE)s;
+    if( e>0 ){
+      while( e>=100  ){ e-=100; r *= 1.0e+100L; }
+      while( e>=10   ){ e-=10;  r *= 1.0e+10L;  }
+      while( e>=1    ){ e-=1;   r *= 1.0e+01L;  }
     }else{
-      while( !(s%10) && e>0 ) e--,s/=10;
+      while( e<=-100 ){ e+=100; r *= 1.0e-100L; }
+      while( e<=-10  ){ e+=10;  r *= 1.0e-10L;  }
+      while( e<=-1   ){ e+=1;   r *= 1.0e-01L;  }
     }
-
-    /* adjust the sign of significand */
-    s = sign<0 ? -s : s;
-
-    /* if exponent, scale significand as appropriate
-    ** and store in result. */
-    if( e ){
-      double scale = 1.0;
-      /* attempt to handle extremely small/large numbers better */
-      if( e>307 && e<342 ){
-        while( e%308 ) { scale *= 1.0e+1; e -= 1; }
-        if( esign<0 ){
-          result = s / scale;
-          result /= 1.0e+308;
-        }else{
-          result = s * scale;
-          result *= 1.0e+308;
-        }
-      }else if( e>=342 ){
-        if( esign<0 ){
-          result = 0.0*s;
-        }else{
-          result = 1e308*1e308*s;  /* Infinity */
-        }
-      }else{
-        /* 1.0e+22 is the largest power of 10 than can be
-        ** represented exactly. */
-        while( e%22 ) { scale *= 1.0e+1; e -= 1; }
-        while( e>0 ) { scale *= 1.0e+22; e -= 22; }
-        if( esign<0 ){
-          result = s / scale;
-        }else{
-          result = s * scale;
-        }
+    assert( r>=0.0 );
+    if( r>+1.7976931348623157081452742373e+308L ){
+#ifdef INFINITY
+      *pResult = +INFINITY;
+#else
+      *pResult = 1.0e308*10.0;
+#endif
+    }else{
+      *pResult = (double)r;
+    }
+  }else{
+    double rr[2];
+    u64 s2;
+    rr[0] = (double)s;
+    s2 = (u64)rr[0];
+#if defined(_MSC_VER) && _MSC_VER<1700
+    if( s2==0x8000000000000000LL ){ s2 = 2*(u64)(0.5*rr[0]); }
+#endif
+    rr[1] = s>=s2 ? (double)(s - s2) : -(double)(s2 - s);
+    if( e>0 ){
+      while( e>=100  ){
+        e -= 100;
+        dekkerMul2(rr, 1.0e+100, -1.5902891109759918046e+83);
+      }
+      while( e>=10   ){
+        e -= 10;
+        dekkerMul2(rr, 1.0e+10, 0.0);
+      }
+      while( e>=1    ){
+        e -= 1;
+        dekkerMul2(rr, 1.0e+01, 0.0);
+      }
+    }else{
+      while( e<=-100 ){
+        e += 100;
+        dekkerMul2(rr, 1.0e-100, -1.99918998026028836196e-117);
+      }
+      while( e<=-10  ){
+        e += 10;
+        dekkerMul2(rr, 1.0e-10, -3.6432197315497741579e-27);
+      }
+      while( e<=-1   ){
+        e += 1;
+        dekkerMul2(rr, 1.0e-01, -5.5511151231257827021e-18);
       }
-    } else {
-      result = (double)s;
     }
+    *pResult = rr[0]+rr[1];
+    if( totypeIsNaN(*pResult) ) *pResult = 1e300*1e300;
   }
+  if( sign<0 ) *pResult = -*pResult;
 
-  /* store the result */
-  *pResult = result;
-
+atof_return:
   /* return true if number and no extra non-whitespace chracters after */
   return z>=zEnd && nDigits>0 && eValid && nonNum==0;
 }

diff --git a/manifest b/manifest
index e6296bc3eec09df049ed6a3fc5a07cf30d01efd2..0e540769aad424c34e68ba2cc3a6e271bb4f4244 100644 (file)
--- a/manifest
+++ b/manifest
@@ -1,5 +1,5 @@
-C Fix\stestcase\sto\saccount\sfor\snew\sbehaviors\swith\sthis\sbranch.
-D 2024-03-14T20:39:03.615
+C Attempt\sto\suse\sthe\sDekker\shigh-precision\sfloating\spoint\salgorithm\sin\sthe\next/misc/totype.c\sextension.
+D 2024-03-15T12:59:11.249
 F .fossil-settings/empty-dirs dbb81e8fc0401ac46a1491ab34a7f2c7c0452f2f06b54ebb845d024ca8283ef1
 F .fossil-settings/ignore-glob 35175cdfcf539b2318cb04a9901442804be81cd677d8b889fcc9149c21f239ea
 F LICENSE.md df5091916dbb40e6e9686186587125e1b2ff51f022cc334e886c19a0e9982724
@@ -419,7 +419,7 @@ F ext/misc/spellfix.c c0aa7b80d6df45f7da59d912b38752bcac1af53a5766966160e6c5cdd3
 F ext/misc/sqlar.c 53e7d48f68d699a24f1a92e68e71eca8b3a9ff991fe9588c2a05bde103c6e7b7
 F ext/misc/stmt.c b090086cd6bd6281c21271d38d576eeffe662f0e6b67536352ce32bbaa438321
 F ext/misc/templatevtab.c 10f15b165b95423ddef593bc5dcb915ec4eb5e0f1066d585e5435a368b8bc22b
-F ext/misc/totype.c 75ed9827d19cc3b434fc2aeb60725d4d46e1534373615612a4d1cfdcc3d60922
+F ext/misc/totype.c e4911af239991f020017b85eaf57122be75c6416789b3946ffcc4f9e861bad22
 F ext/misc/uint.c 053fed3bce2e89583afcd4bf804d75d659879bbcedac74d0fa9ed548839a030b
 F ext/misc/unionvtab.c 716d385256d5fb4beea31b0efede640807e423e85c9784d21d22f0cce010a785
 F ext/misc/urifuncs.c f71360d14fa9e7626b563f1f781c6148109462741c5235ac63ae0f8917b9c751
@@ -2179,8 +2179,11 @@ F vsixtest/vsixtest.tcl 6a9a6ab600c25a91a7acc6293828957a386a8a93
 F vsixtest/vsixtest.vcxproj.data 2ed517e100c66dc455b492e1a33350c1b20fbcdc
 F vsixtest/vsixtest.vcxproj.filters 37e51ffedcdb064aad6ff33b6148725226cd608e
 F vsixtest/vsixtest_TemporaryKey.pfx e5b1b036facdb453873e7084e1cae9102ccc67a0
-P d543c829ef74dbd64105bd757ca660e4f02e9ce562be4f1688a701fa535351c4
-R 9a3605a04234115bfc68cd6fc0f2125f
+P 823e579362c05bb8accf6c3b158c5162a16eb23cf81d6021c9e3246e32583d1c
+R 52f253ac14fb8e14beebf78670890719
+T *branch * totype-dekker
+T *sym-totype-dekker *
+T -sym-exp-values-clause2 *
 U drh
-Z eb4b6e81254ae60e2fa77bddac123f09
+Z a6f29290405dcccc6b403f65943397c8
 # Remove this line to create a well-formed Fossil manifest.

diff --git a/manifest.uuid b/manifest.uuid
index ad17c6052ce57ecfff41202bd815b7dae5e4500a..3d4c0f90f197e537af1b5c54759c8438c0ab3837 100644 (file)
--- a/manifest.uuid
+++ b/manifest.uuid
@@ -1 +1 @@
-823e579362c05bb8accf6c3b158c5162a16eb23cf81d6021c9e3246e32583d1c
\ No newline at end of file
+dd41db3b106345c510f1f64aebcedf162701e9e21bd157441331b454ce11b353
\ No newline at end of file