From c8885d42fe4799a265c69f4a6bbd4cdb3cf4daea Mon Sep 17 00:00:00 2001
From: drh <>
Date: Mon, 3 Jul 2023 09:07:23 +0000
Subject: [PATCH] Replace the dodgy error estimating logic in the previous
 check-in with full-up Dekker double-double multiplication, and this idea
 works much better. There are still a few inaccuracies, but it is much closer.

FossilOrigin-Name: 4fa6938dac2d3d813a37664053db31451a2a065f78dd212488f5f7f8d583ddc5
---
 manifest        | 17 +++++------
 manifest.uuid   |  2 +-
 src/util.c      | 77 +++++++++++++++++++++++++++++++++++++------------
 test/atof1.test |  6 ++--
 4 files changed, 69 insertions(+), 33 deletions(-)

diff --git a/manifest b/manifest
index f969dca2f8..a420e218ba 100644
--- a/manifest
+++ b/manifest
@@ -1,5 +1,5 @@
-C An\sidea\sfor\simproving\saccuracy\sof\sfp-to-decimal\sconversion\sfor\ssystems\sthat\ndo\snot\shave\s"long\sdouble".\s\sDoes\snot\squite\swork.\s\sRetained\sonly\sfor\ndocumentation\spurposes.
-D 2023-07-03T00:40:37.359
+C Replace\sthe\sdodgy\serror\sestimating\slogic\sin\sthe\sprevious\scheck-in\swith\nfull-up\sDekker\sdouble-double\smultiplication,\sand\sthis\sidea\sworks\smuch\sbetter.\nThere\sare\sstill\sa\sfew\sinaccuracies,\sbut\sit\sis\smuch\scloser.
+D 2023-07-03T09:07:23.117
 F .fossil-settings/empty-dirs dbb81e8fc0401ac46a1491ab34a7f2c7c0452f2f06b54ebb845d024ca8283ef1
 F .fossil-settings/ignore-glob 35175cdfcf539b2318cb04a9901442804be81cd677d8b889fcc9149c21f239ea
 F LICENSE.md df5091916dbb40e6e9686186587125e1b2ff51f022cc334e886c19a0e9982724
@@ -705,7 +705,7 @@ F src/trigger.c ad6ab9452715fa9a8075442e15196022275b414b9141b566af8cdb7a1605f2b0
 F src/update.c 0aa36561167a7c40d01163238c297297962f31a15a8d742216b3c37cdf25f731
 F src/upsert.c 5303dc6c518fa7d4b280ec65170f465c7a70b7ac2b22491598f6d0b4875b3145
 F src/utf.c ee39565f0843775cc2c81135751ddd93eceb91a673ea2c57f61c76f288b041a0
-F src/util.c 12741de9969af5134488b7c89a97ecc22144f93f62c6255f797abb9b3121608a
+F src/util.c 1feadae2da54a59ab7ec82d4b94ed22f94bb5d4df25e41e73c05c8f3c1328607
 F src/vacuum.c 604fcdaebe76f3497c855afcbf91b8fa5046b32de3045bab89cc008d68e40104
 F src/vdbe.c 74282a947234513872a83b0bab1b8c644ece64b3e27b053ef17677c8ff9c81e0
 F src/vdbe.h 41485521f68e9437fdb7ec4a90f9d86ab294e9bb8281e33b235915e29122cfc0
@@ -775,7 +775,7 @@ F test/async2.test c0a9bd20816d7d6a2ceca7b8c03d3d69c28ffb8b
 F test/async3.test d73a062002376d7edc1fe3edff493edbec1fc2f7
 F test/async4.test 1787e3952128aa10238bf39945126de7ca23685a
 F test/async5.test 383ab533fdb9f7ad228cc99ee66e1acb34cc0dc0
-F test/atof1.test 1f24a6f16e44892848298bd8e488159d4a60cf47be003d32747a136103c9bbac
+F test/atof1.test b4dbc37187649e23ba92e06e0d2c076304516967b5a8b0af50f8e4f2323e654c
 F test/atomic.test 065a453dde33c77ff586d91ccaa6ed419829d492dbb1a5694b8a09f3f9d7d061
 F test/atomic2.test b6863b4aa552543874f80b42fb3063f1c8c2e3d8e56b6562f00a3cc347b5c1da
 F test/atrc.c c388fac43dbba05c804432a7135ae688b32e8f25818e9994ffba4b64cf60c27c
@@ -2041,11 +2041,8 @@ F vsixtest/vsixtest.tcl 6a9a6ab600c25a91a7acc6293828957a386a8a93
 F vsixtest/vsixtest.vcxproj.data 2ed517e100c66dc455b492e1a33350c1b20fbcdc
 F vsixtest/vsixtest.vcxproj.filters 37e51ffedcdb064aad6ff33b6148725226cd608e
 F vsixtest/vsixtest_TemporaryKey.pfx e5b1b036facdb453873e7084e1cae9102ccc67a0
-P ce06982f880339cf46704e95c907249827c3e44af2b9420005200ca8abd3f371
-R ae7783ab51788007ac81a8600cb6b0a8
-T *branch * fp-to-decimal-branch1
-T *sym-fp-to-decimal-branch1 *
-T -sym-fp-to-decimal-refactor *
+P 702243333843219f8904ee0fd12283080d6c6b3fc1ffb36c534cdefb3563c00d
+R 0bc2956a2833c937de8e3462962e6784
 U drh
-Z 0f216e6dd96bccb8187eb2a64c983082
+Z 2f85ae3f7bbd29aeb71b560da07ddd48
 # Remove this line to create a well-formed Fossil manifest.
diff --git a/manifest.uuid b/manifest.uuid
index d8fd936394..8bb5c10510 100644
--- a/manifest.uuid
+++ b/manifest.uuid
@@ -1 +1 @@
-702243333843219f8904ee0fd12283080d6c6b3fc1ffb36c534cdefb3563c00d
\ No newline at end of file
+4fa6938dac2d3d813a37664053db31451a2a065f78dd212488f5f7f8d583ddc5
\ No newline at end of file
diff --git a/src/util.c b/src/util.c
index 9bc37c4f9a..b7d6bd58e8 100644
--- a/src/util.c
+++ b/src/util.c
@@ -928,6 +928,31 @@ int sqlite3Atoi(const char *z){
   return x;
 }
 
+/* Double-Double multiplication.  *(z,zz) = (x,xx) * (y,yy)
+**
+** Reference:
+**   T. J. Dekker, "A Floating-Point Technique for Extending the
+**   Available Precision".  1971-07-26.
+*/
+static void mul2(
+  double x,  double xx,
+  double y,  double yy,
+  double *z, double *zz
+){
+  double hx, tx, hy, ty, p, q, c, cc;
+  p = 67108865.0 * x;
+  hx = x - p + p; tx = x - hx;
+  p = 67108865.0 * y;
+  hy = y - p + p; ty = y - hy;
+  p = hx*hy;
+  q = hx*ty + tx*hy;
+  c = p+q;
+  cc = p - c + q + tx*ty;
+  cc = x*yy + xx*y + cc;
+  *z = c + cc;
+  *zz = c - *z + cc;
+}
+
 /*
 ** Decode a floating-point value into an approximate decimal
 ** representation.
@@ -969,10 +994,7 @@ void sqlite3FpDecode(FpDecode *p, double r, int iRound, int mxRound){
   /* Multiply r by powers of ten until it lands somewhere in between
   ** 1.0e+19 and 1.0e+17.
   */
-  if( sizeof(long double)>8
-   || (v & 0x0000000000ffffffLL)==0
-   || (v & 0xffffffffff000000LL)==0
-  ){
+  if( sizeof(long double)>8 && 0 ){
     long double rr = r;
     if( rr>=1.0e+19 ){
       while( rr>=1.0e+119L ){ exp+=100; rr *= 1.0e-100L; }
@@ -985,23 +1007,40 @@ void sqlite3FpDecode(FpDecode *p, double r, int iRound, int mxRound){
     }
     v = (u64)rr;
   }else{
-    double r2, r3;
-    v &= 0xffffffffff000000LL;
-    memcpy(&r2, &v, 8);
-    assert( r2<r );
-    assert( r2>0.0 );
-    r3 = r - r2;
-    assert( r3>0.0 );
-    if( r2>=1.0e+19 ){
-      while( r2>=1.0e+119L ){ exp+=100; r2 *= 1.0e-100; r3 *= 1.0e-100; }
-      while( r2>=1.0e+29L  ){ exp+=10;  r2 *= 1.0e-10;  r3 *= 1.0e-10;  }
-      while( r2>=1.0e+19L  ){ exp++;    r2 *= 1.0e-1;   r3 *= 1.0e-1;   }
+    double rr = 0.0;
+    if( r>=1.0e+19 ){
+      while( r>=1.0e+119 ){
+        if( r>=1.0e+300 ){
+          exp += 300;
+          r *= 1.0e-300;
+          break;
+        }
+        exp += 100;
+        mul2(r, rr, 1.0e-100, -1.99918998026028836196e-117, &r, &rr);
+      }
+      while( r>=1.0e+29 ){
+        exp += 10;
+        mul2(r,rr, 1.0e-10, -3.6432197315497741579e-27, &r, &rr);
+      }
+      while( r>=1.0e+19 ){
+        exp += 1;
+        mul2(r,rr, 1.0e-01, -5.5511151231257827021e-18, &r, &rr);
+      }
     }else{
-      while( r2<1.0e-97L   ){ exp-=100; r2 *= 1.0e+100; r3 *= 1.0e+100; }
-      while( r2<1.0e+07L   ){ exp-=10;  r2 *= 1.0e+10;  r3 *= 1.0e+10;  }
-      while( r2<1.0e+17L   ){ exp--;    r2 *= 1.0e+1;   r3 *= 1.0e+1;   }
+      while( r<1.0e-98  ){
+        exp -= 100;
+        mul2(r, rr, 1.0e+100, -1.5902891109759918046, &r, &rr);
+      }
+      while( r<1.0e+08  ){
+        exp -= 10;
+        mul2(r, rr, 1.0e+10, 0.0, &r, &rr);
+      }
+      while( r<1.0e+18  ){
+        exp -= 1;
+        mul2(r, rr, 1.0e+01, 0.0, &r, &rr);
+      }
     }
-    v = (u64)(r2+r3);
+    v = (u64)(r)+(u64)(rr);
   }
 
 
diff --git a/test/atof1.test b/test/atof1.test
index 5959c5d5b9..ecfbdd33ee 100644
--- a/test/atof1.test
+++ b/test/atof1.test
@@ -51,10 +51,10 @@ for {set i 1} {$i<20000} {incr i} {
     set y [db eval {SELECT $x=CAST(quote($x) AS real)}]
     if {!$y} {
       db eval {SELECT real2hex($x) a, real2hex(CAST(quote($x) AS real)) b} {}
-      puts "\nIN:    $a $xf"
-      puts [format {QUOTE: %16s %s} {} [db eval {SELECT quote($x)}]]
+      puts "\n!IN:    $a $xf"
+      puts [format {!QUOTE: %16s %s} {} [db eval {SELECT quote($x)}]]
       db eval {SELECT CAST(quote($x) AS real) c} {}
-      puts "OUT:   $b [format %.32e $c]"
+      puts "!OUT:   $b [format %.32e $c]"
     }
     set y
   } {1}
-- 
2.47.2