--- a/QDouble.st	Wed Jun 21 10:34:20 2017 +0200
+++ b/QDouble.st	Wed Jun 21 20:48:59 2017 +0200
@@ -16,7 +16,7 @@
 LimitedPrecisionReal variableByteSubclass:#QDouble
 	instanceVariableNames:''
 	classVariableNames:'DefaultPrintFormat E Epsilon FMax FMin InvFact Ln10 Ln2 Pi
-		QDoubleOne QDoubleZero'
+		QDoubleOne QDoubleZero NaN'
 	poolDictionaries:''
 	category:'Magnitude-Numbers'
 !
@@ -841,6 +841,17 @@
 
 !QDouble class methodsFor:'constants'!
 
+NaN
+    "return a QDouble which represents not-a-Number (i.e. an invalid number)"
+
+    NaN isNil ifTrue:[
+        NaN := self d0:(Float NaN) d1:(Float NaN) d2:(Float NaN) d3:(Float NaN)
+    ].
+    ^ NaN
+
+    "Created: / 21-06-2017 / 20:44:57 / cg"
+!
+
 e
     "return the constant e as quad precision double.
      (returns approx. 200 bits of precision)"
@@ -1053,7 +1064,7 @@
 defaultPrintPrecision
     "return the number of decimal digits printed by default"
 
-    ^ 10
+    ^ 20
 
     "
      ShortFloat defaultPrintPrecision
@@ -1063,6 +1074,7 @@
     "
 
     "Created: / 17-06-2017 / 02:58:51 / cg"
+    "Modified: / 21-06-2017 / 13:39:08 / cg"
 !
 
 epsilon
@@ -1072,13 +1084,15 @@
     ^ 1.2154326714572500565324311366323150942261000827598106963711353e-63
 
     "
-     Float epsilon
-     ShortFloat epsilon
-     QDouble epsilon
+     Float epsilon       -> 2.22044604925031E-16
+     ShortFloat epsilon  -> 1.19209289550781E-07
+     LongFloat epsilon   -> 1.0842021724855E-19
+     QDouble epsilon     -> 1.21543267145725E-63
     "
 
     "Created: / 12-06-2017 / 18:52:44 / cg"
     "Modified: / 14-06-2017 / 19:12:23 / cg"
+    "Modified (comment): / 21-06-2017 / 13:56:34 / cg"
 !
 
 numBitsInExponent
@@ -1899,6 +1913,27 @@
 !QDouble methodsFor:'mathematical functions'!
 
 exp
+    "/ the exp_sloppy code is the algorithm from the original C++ qd package;
+    "/ however, it is inexact in the 37th digit
+    "/ Therefore, use the inherited code, which is slow, but more precise.
+
+    ^ super exp
+
+    "
+     2.0 exp                -> 7.38905609893065 
+     2.0 asQDouble exp      -> 7.38905609893065022723
+    "
+
+    "Created: / 19-06-2017 / 01:49:32 / cg"
+    "Modified (comment): / 21-06-2017 / 13:50:38 / cg"
+!
+
+exp_sloppy
+    "/ this is the algorithm from the qd C++ package;
+    "/ however, it is inexact in the 37th digit
+    "/
+    "/ the inherited exp code is much slower, but more precise.
+    
     "/ Strategy:  We first reduce the size of x by noting that
     "/
     "/  exp(kr + m * log(2)) = 2^m * exp(r)^k
@@ -1933,13 +1968,11 @@
     k := 1.0 ldexp:16.
     inv_k := 1.0 / k.
 
-    "/  double m = std::floor(a.x[0] / qd_real::_log2.x[0] + 0.5);
     m := (d0 / Float ln2 + 0.5) floor.
 
     mul_pwr2 := [:a :b | QDouble d0:(a d0 * b) d1:(a d1 * b) d2:(a d2 * b) d3:(a d3 * b) ].
     mul2 := [:a | QDouble d0:(a d0 * 2.0) d1:(a d1 * 2.0) d2:(a d2 * 2.0) d3:(a d3 * 2.0) ].
 
-    "/  qd_real r = mul_pwr2(a - qd_real::_log2 * m, inv_k);
     r := mul_pwr2 value:(self - (self class ln2 * m)) value:inv_k.
     thresh := inv_k * eps.
 
@@ -1951,7 +1984,7 @@
         t := p * (self class invFact at:i).
         i := i+1.
         s := s + t.
-    ] doWhile:[ (t asFloat abs > thresh) and:[i < 10] ]. "/ (std::abs(to_double(t)) > thresh && i < 9);
+    ] doWhile:[ (t asFloat abs > thresh) and:[i < 10] ]. 
 
     s := (mul2 value:s) + s squared.
     s := (mul2 value:s) + s squared.
@@ -1978,8 +2011,7 @@
      1.0 asQDouble exp -> 2.71828182845905
     "
 
-    "Created: / 19-06-2017 / 01:49:32 / cg"
-    "Modified: / 19-06-2017 / 19:10:15 / cg"
+    "Created: / 21-06-2017 / 13:48:27 / cg"
 !
 
 floor