--- a/QDouble.st Wed Jun 14 19:22:17 2017 +0200
+++ b/QDouble.st Thu Jun 15 09:58:16 2017 +0200
@@ -1,3 +1,5 @@
+"{ Encoding: utf8 }"
+
"
COPYRIGHT (c) 2017 by eXept Software AG
All Rights Reserved
@@ -1045,7 +1047,7 @@
- aNumber
"return the sum of the receiver and the argument, aNumber"
- ^ (aNumber negated) sumFromQDouble:self
+ ^ self + (aNumber negated)
"
(QDouble fromFloat:1e20) asDoubleArray
@@ -1053,9 +1055,13 @@
(QDouble fromFloat:1e-20) asDoubleArray
((QDouble fromFloat:1e-20) - (QDouble fromFloat:1.0)) asDoubleArray
((QDouble fromFloat:2.0) - (QDouble fromFloat:1.0)) asDoubleArray
+
+ ((QDouble fromFloat:2.0) - (QDouble fromFloat:1.0) + (QDouble fromFloat:1.0)) asDoubleArray
+ ((QDouble fromFloat:1e-20) - (QDouble fromFloat:1.0) + (QDouble fromFloat:1.0)) asDoubleArray
"
"Created: / 12-06-2017 / 23:41:39 / cg"
+ "Modified (comment): / 15-06-2017 / 00:34:41 / cg"
!
/ aNumber
@@ -1065,9 +1071,14 @@
"
((QDouble fromFloat:1e20) / (QDouble fromFloat:2.0)) asDoubleArray
+
+ ((QDouble fromFloat:1.2345) / (QDouble fromFloat:10.0)) asDoubleArray
+ ((QDouble fromFloat:1.2345) / 10.0) asDoubleArray
+
"
"Created: / 13-06-2017 / 17:59:09 / cg"
+ "Modified (comment): / 15-06-2017 / 00:14:26 / cg"
! !
!QDouble methodsFor:'coercing & converting'!
@@ -1409,19 +1420,28 @@
|q0 q1 q2 q3 r|
- q0 := self d0 / aQDouble d0.
- r := self - (aQDouble * q0).
+ q0 := aQDouble d0 / self d0.
+ "/ Stdout showCR:('q0: %1 (a[0]=%2; b[0]=%3)\n' bindWith:q0 with:self d0 with:aQDouble d0).
+ r := aQDouble - (self * q0).
+ "/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
- q1 := r d0 / aQDouble d0.
- r := r - (aQDouble * q1).
+ q1 := r d0 / self d0.
+ "/ Stdout showCR:('q1: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q1 with:r d0 with:aQDouble d0).
+ r := r - (self * q1).
+ "/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
- q2 := r d0 / aQDouble d0.
- r := r - (aQDouble * q2).
+ q2 := r d0 / self d0.
+ "/ Stdout showCR:('q2: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q2 with:r d0 with:aQDouble d0).
+ r := r - (self * q2).
+ "/ Stdout showCR:('r: %1\n' bindWith:r asDoubleArray).
- q3 := r d0 / aQDouble d0.
+ q3 := r d0 / self d0.
+ "/ Stdout showCR:('q3: %1 (r[0]=%2; b[0]=%3)\n' bindWith:q3 with:r d0 with:aQDouble d0).
r := self class d0:q0 d1:q1 d2:q2 d3:q3.
+ "/ Stdout showCR:('before renorm: %1\n' bindWith:r asDoubleArray).
r renorm.
+ "/ Stdout showCR:('after renorm: %1\n' bindWith:r asDoubleArray).
^ r
"
@@ -1437,9 +1457,33 @@
((QDouble fromFloat:1e20 / 2.0)) asFloat
((1e20 + (QDouble fromFloat:1.0) + 1e-20) / 2.0) asDoubleArray
+
+ ((QDouble fromFloat:10.0) quotientFromQDouble: (QDouble fromFloat:1.234)) asDoubleArray
+ ((QDouble fromFloat:1.234) / (QDouble fromFloat:10.0)) asDoubleArray
+
+q0: 1.234000e-01 (a[0]=1.234000e+00; b[0]=1.000000e+01)
+a: 1.234000e+00/0.000000e+00/0.000000e+00/0.000000e+00
+b: 1.000000e+01/0.000000e+00/0.000000e+00/0.000000e+00
+(b * q0): 1.234000e+00/-2.775558e-17/0.000000e+00/0.000000e+00
+r: 2.775558e-17/0.000000e+00/0.000000e+00/0.000000e+00
+
+q1: 2.775558e-18 (r[0]=2.775558e-17; b[0]=1.000000e+01)
+r: -1.540744e-33/0.000000e+00/0.000000e+00/0.000000e+00
+
+q2: -1.540744e-34 (r[0]=-1.540744e-33; b[0]=1.000000e+01)
+r: 8.552847e-50/0.000000e+00/0.000000e+00/0.000000e+00
+
+q3: 8.552847e-51 (r[0]=8.552847e-50; b[0]=1.000000e+01)
+
+before renorm: 1.234000e-01/2.775558e-18/-1.540744e-34/8.552847e-51
+after renorm: 1.234000e-01/2.775558e-18/-1.540744e-34/8.552847e-51
+1.234/10.0 is: 0.123400 / 0.000000 / -0.000000 / 0.000000
+
+
"
"Created: / 13-06-2017 / 17:50:35 / cg"
+ "Modified (comment): / 15-06-2017 / 01:02:05 / cg"
!
sumFromFloat:aFloat
@@ -1492,7 +1536,7 @@
#ifndef QD_IEEE_ADD
// sloppy_add...
- /*
+# if 0
double s0, s1, s2, s3;
double t0, t1, t2, t3;
int savedCV;
@@ -1508,10 +1552,12 @@
m_three_sum2(s3, t0, t2);
t0 = t0 + t1 + t3;
+ m_renorm5(s0, s1, s2, s3, t0);
fpu_fix_end(&savedCV);
- m_renorm5(s0, s1, s2, s3, t0);
- return qd_real(s0, s1, s2, s3, t0);
- */
+
+ __qNew_qdReal(newQD, s0, s1, s2, s3);
+ RETURN(newQD);
+# else
/* Same as above, but addition re-organized to minimize
data dependency ... unfortunately some compilers are
@@ -1564,6 +1610,7 @@
fpu_fix_end(&savedCV);
__qNew_qdReal(newQD, s0, s1, s2, s3);
RETURN(newQD);
+# endif
#else
// ieee_add...
int i, j, k;
@@ -1637,7 +1684,7 @@
"
"Created: / 12-06-2017 / 21:15:43 / cg"
- "Modified: / 14-06-2017 / 11:44:53 / cg"
+ "Modified: / 15-06-2017 / 00:49:10 / cg"
! !
!QDouble methodsFor:'inspecting'!
@@ -1784,46 +1831,460 @@
!QDouble methodsFor:'printing & storing'!
-printDigitsOn:outStream precision:precision exponentInto:expn
- |numDigits r e i d|
+digitsWithPrecision:precision
+ "generate the digits and exponent"
+
+ |numDigits r exp i d out str|
+
+ numDigits := precision+1. "/ number of digits
+ r := self abs.
+ self d0 = 0.0 ifTrue:[
+ ^ { String new:precision withAll:$0 . 0 }
+ ].
+
+ out := WriteStream on:(String new:precision+5).
+
+ "/ determine approx. exponent
+ exp := self d0 abs log10 floor.
+ exp < -300 ifTrue:[
+ r := r * (10.0 raisedToInteger:300).
+ r := r / (10.0 raisedToInteger:(exp+300)).
+ ] ifFalse:[
+ exp > 300 ifTrue:[
+self halt.
+"/ r := ldexpr(r, -53);
+"/ r := r / qd_real(10.0) ^ exp;
+"/ r := ldexp(r, 53);
+ ] ifFalse:[
+ r := r / (10.0 asQDouble raisedTo:exp).
+ ]
+ ].
+
+ "/ Fix exponent if we are off by one
+ (r >= 10.0) ifTrue:[
+ r := r / 10.0.
+ exp := exp + 1.
+ ] ifFalse:[
+ (r < 1.0) ifTrue:[
+ r := r * 10.0.
+ exp := exp - 1.
+ ]
+ ].
+
+ ((r >= 10.0) or:[ r < 1.0 ]) ifTrue:[
+ self error:'can''t compute exponent.'.
+ ].
+
+ "/ Extract the digits
+ 1 to:numDigits do:[:i |
+ d := r d0 truncated.
+ r := r - d.
+ r := r * 10.0.
+
+ out nextPut:($0 + d).
+ ].
+
+ str := out contents.
+
+ "/ Fix out-of-range digits.
+ "/ huh? what is this????
+ "/ how can an "out-of-range" digit appear???
+ numDigits to:2 by:-1 do:[:i |
+ (str at:i) < $0 ifTrue:[
+self halt.
+ str at:i-1 put:(str at:i-1) - 1.
+ str at:i put:(str at:i) + 10.
+ ] ifFalse:[
+ (str at:i) > $9 ifTrue:[
+self halt.
+ str at:i-1 put:(str at:i-1) + 1.
+ str at:i put:(str at:i) - 10.
+ ].
+ ].
+ ].
+
+ str first <= $0 ifTrue:[
+ self error:'non-positive leading digit'
+ ].
+
+ "/ Round, handle carry
+ (str at:numDigits) >= $5 ifTrue:[
+ str at:numDigits-1 put:(str at:numDigits-1) + 1.
+ i := numDigits-1.
+ [i > 1 and:[(str at:i) > $9]] whileTrue:[
+ str at:i put:(str at:i) - 10.
+ i := i - 1.
+ str at:i put:(str at:i) + 1.
+ ]
+ ].
+
+ "/ If first digit is 10, shift everything.
+ str first > $9 ifTrue:[
+ exp := exp + 1.
+ str at:1 put:$0.
+ str := '1',str
+ ].
+ ^ { (str copyTo:precision) . exp }
+
+ "
+ 1.2345 asQDouble digitsWithPrecision:5 -> #('12345' 0)
+ 12.345 asQDouble digitsWithPrecision:5 -> #('12345' 1)
+ 12345 asQDouble digitsWithPrecision:5 -> #('12345' 4)
+ 12345.1 asQDouble digitsWithPrecision:5 -> #('12345' 4)
+ 12345.9 asQDouble digitsWithPrecision:5 -> #('12346' 4)
+
+ 1.2345 asQDouble / 10.0 asQDouble
+ 1.2345 asQDouble / 10.0
+ "
+
+ "Created: / 15-06-2017 / 09:10:01 / cg"
+!
+
+digitsWithPrecision:precision exponentInto:expnHolder
+ "generate the digits and return exponent info"
+
+ |numDigits r exp i d out str|
numDigits := precision+1. "/ number of digits
r := self abs.
self d0 = 0.0 ifTrue:[
- expn value:0.
- precision timesRepeat:[ outStream nextPut:$0 ].
+ expnHolder value:0.
+ ^ String new:precision withAll:$0
+ ].
+
+ out := WriteStream on:(String new:precision+5).
+
+ "/ determine approx. exponent
+ exp := self d0 abs log10 floor.
+ exp < -300 ifTrue:[
+ r := r * (10.0 raisedToInteger:300).
+ r := r / (10.0 raisedToInteger:(exp+300)).
+ ] ifFalse:[
+ exp > 300 ifTrue:[
+self halt.
+"/ r := ldexpr(r, -53);
+"/ r := r / qd_real(10.0) ^ exp;
+"/ r := ldexp(r, 53);
+ ] ifFalse:[
+ r := r / (10.0 asQDouble raisedTo:exp).
+ ]
+ ].
+
+ "/ Fix exponent if we are off by one
+ (r >= 10.0) ifTrue:[
+ r := r / 10.0.
+ exp := exp + 1.
+ ] ifFalse:[
+ (r < 1.0) ifTrue:[
+ r := r * 10.0.
+ exp := exp - 1.
+ ]
+ ].
+
+ ((r >= 10.0) or:[ r < 1.0 ]) ifTrue:[
+ self error:'can''t compute exponent.'.
+ ].
+
+ "/ Extract the digits
+ 1 to:numDigits do:[:i |
+ d := r d0 truncated.
+ r := r - d.
+ r := r * 10.0.
+
+ out nextPut:($0 + d).
+ ].
+
+ str := out contents.
+
+ "/ Fix out-of-range digits.
+ "/ huh? what is this????
+ "/ how can an "out-of-range" digit appear???
+ i := numDigits.
+ [i > 1] whileTrue:[
+ (str at:i) < $0 ifTrue:[
+ self halt.
+ str at:i-1 put:(str at:i-1) - 1.
+ str at:i put:(str at:i) + 10.
+ ] ifFalse:[
+ (str at:i) > $9 ifTrue:[
+ self halt.
+ str at:i-1 put:(str at:i-1) + 1.
+ str at:i put:(str at:i) - 10.
+ ].
+ ].
+ i := i - 1.
+ ].
+
+ (str at:1) <= $0 ifTrue:[
+ self error:'non-positive leading digit'
+ ].
+
+ "/ Round, handle carry
+ (str at:numDigits) >= $5 ifTrue:[
+ str at:numDigits-1 put:(str at:numDigits-1) + 1.
+ i := numDigits-1.
+ [i > 1 and:[(str at:i) > $9]] whileTrue:[
+ str at:i put:(str at:i) - 10.
+ i := i - 1.
+ str at:i put:(str at:i) + 1.
+ ]
+ ].
+
+ "/ If first digit is 10, shift everything.
+ (str at:1) > $9 ifTrue:[
+ exp := exp + 1.
+ str at:1 put:$0.
+ str := '1',str
+ ].
+ expnHolder value:exp.
+ ^ (str copyTo:precision).
+
+ "
+ 1.2345 asQDouble digitsWithPrecision:5 exponentInto:[:e | Transcript showCR:e] -> '12345' / 0
+ 12.345 asQDouble digitsWithPrecision:5 exponentInto:[:e | Transcript showCR:e] -> '12345' / 1
+ 12345 asQDouble digitsWithPrecision:5 exponentInto:[:e | Transcript showCR:e]
+ 12345.1 asQDouble digitsWithPrecision:5 exponentInto:[:e | Transcript showCR:e]
+ 12345.9 asQDouble digitsWithPrecision:5 exponentInto:[:e | Transcript showCR:e]
+
+ 1.2345 asQDouble / 10.0 asQDouble
+ 1.2345 asQDouble / 10.0
+ "
+
+ "Created: / 15-06-2017 / 09:05:42 / cg"
+!
+
+printOn:aStream
+ "return a printed representation of the receiver"
+
+ self d1 = 0.0 ifTrue:[
+ self d0 printOn:aStream.
+ ^ self.
+ ].
+ self
+ printOn:aStream precision:0 width:0
+ fixed:false showPositive:false uppercase:false fillChar:(Character space)
+
+ "
+ (QDouble fromFloat:1.2345) printOn:Transcript
+ ((QDouble fromFloat:1.2345) / 10.0) printOn:Transcript
+ "
+
+ "Created: / 15-06-2017 / 01:51:36 / cg"
+ "Modified (comment): / 15-06-2017 / 08:50:25 / cg"
+!
+
+printOn:aStream precision:precision width:width
+ fixed:fixed showPositive:showPositive uppercase:uppercase fillChar:fillChar
+
+ "return a printed representation of the receiver.
+ This is a parametrized entry, which can be used by printf-like functions"
+
+ |sgn count delta exp|
+
+ self d1 = 0.0 ifTrue:[
+ self d0 printOn:aStream.
^ self.
].
- "/ determine approx. exponent
- e := self d0 abs log10 floor.
- self halt.
- e < -300 ifTrue:[
- r := r * (10.0 raisedToInteger:300).
- r := r / (10.0 raisedToInteger:(e+300)).
+ count := 0.
+ sgn := true.
+ exp := 0.
+
+ self isInfinite ifTrue:[
+ self < 0 ifTrue:[
+ aStream nextPut:$-.
+ count := 1.
+ ] ifFalse:[
+ showPositive ifTrue:[
+ aStream nextPut:$+.
+ count := 1.
+ ] ifFalse:[
+ sgn := false.
+ ].
+ ].
+ uppercase ifTrue:[
+ aStream nextPutAll:'INF'
+ ] ifFalse:[
+ aStream nextPutAll:'inf'
+ ].
+ count := count + 3.
] ifFalse:[
- e > 300 ifTrue:[
- ] ifFalse:[
+ self isNaN ifTrue:[
+ uppercase ifTrue:[
+ aStream nextPutAll:'NAN'
+ ] ifFalse:[
+ aStream nextPutAll:'nan'
+ ].
+ count := count + 3.
+ ] ifFalse:[
+ self < 0 ifTrue:[
+ aStream nextPut:$-.
+ count := count + 1.
+ ] ifFalse:[
+ showPositive ifTrue:[
+ aStream nextPut:$+.
+ count := count + 1.
+ ] ifFalse:[
+ sgn := false.
+ ].
+ ].
+ self = 0.0 ifTrue:[
+ aStream nextPut:$0.
+ count := count + 1.
+ precision > 0 ifTrue:[
+ aStream nextPut:$..
+ count := count + 1.
+ precision timesRepeat:[ aStream nextPut:$0 ].
+ count := count + precision.
+ ].
+ self halt.
+ ] ifFalse:[
+ |off d d_width_extra|
+
+ "/ non-zero case
+ off := fixed
+ ifTrue:[ 1 + self abs log10 floor ]
+ ifFalse:[1].
+ d := precision + off.
+ d_width_extra := d.
+ fixed ifTrue:[
+ d_width_extra := 120 max:d.
+ ].
+ "/ highly special case - fixed mode, precision is zero, abs(*this) < 1.0
+ "/ without this trap a number like 0.9 printed fixed with 0 precision prints as 0
+ "/ should be rounded to 1.
+ (fixed and:[ (precision == 0) and:[ (self abs < 1.0) ]]) ifTrue:[
+ (self abs >= 0.5) ifTrue:[
+ aStream nextPut:$1
+ ] ifFalse:[
+ aStream nextPut:$0
+ ].
+ ^ self
+ ].
+
+ "/ handle near zero to working precision (but not exactly zero)
+ (fixed and:[ d <= 0 ]) ifTrue:[
+ aStream nextPut:$0.
+ (precision > 0) ifTrue:[
+ aStream nextPut:$. .
+ aStream next:precision put:$0.
+ ]
+ ] ifFalse:[
+ "/ default
+
+ |t j|
+
+ fixed ifTrue:[
+ t := String streamContents:[:s |
+ self
+ printDigitsOn:s
+ precision:d_width_extra
+ exponentInto:[:e | exp := e]
+ ].
+ ] ifFalse:[
+ t := String streamContents:[:s |
+ self
+ printDigitsOn:s
+ precision:d
+ exponentInto:[:e | exp := e]
+ ].
+ ].
+
+ fixed ifTrue:[
+ "/ fix the string if it's been computed incorrectly
+ "/ round here in the decimal string if required
+ t := self round_string_qd:t at:(d + 1) offsetInto:[:o | off := o].
+
+ (off > 0) ifTrue:[
+ aStream next:off putAll:t startingAt:0.
+ (precision > 0) ifTrue:[
+ aStream nextPut:$. .
+ aStream next:precision putAll:t startingAt:off+1.
+ ]
+ ] ifFalse:[
+ aStream nextPutAll:'0.'.
+ (off < 0) ifTrue:[
+ aStream next:off negated put:$0.
+ ].
+ aStream next:d putAll:t startingAt:0.
+ ]
+ ] ifFalse:[
+ aStream nextPut:(t at:1).
+ (precision > 0) ifTrue:[
+ aStream nextPut:$. .
+ ].
+ aStream next:precision putAll:t startingAt:2.
+ ]
+ ].
+ ].
]
].
-
-
- "
- 1.2345 asQDouble printDigitsOn:Transcript precision:5 exponentInto:[:e |]
- "
+
+ "/ trap for improper offset with large values
+ "/ without this trap, output of values of the for 10^j - 1 fail for j > 28
+ "/ and are output with the point in the wrong place, leading to a dramatically off value
- "Created: / 13-06-2017 / 17:28:08 / cg"
+"/ (fixed and:[ (precision > 0) ]) ifTrue:[
+"/ "/ make sure that the value isn't dramatically larger
+"/ from_string = atof(s.c_str());
+"/
+"/ // if this ratio is large, then we've got problems
+"/ if( fabs( from_string / this->x[0] ) > 3.0 ){
+"/
+"/ int point_position;
+"/ char temp;
+"/
+"/ // loop on the string, find the point, move it up one
+"/ // don't act on the first character
+"/ for(i=1; i < s.length(); i++){
+"/ if(s[i] == '.'){
+"/ s[i] = s[i-1] ;
+"/ s[i-1] = '.' ;
+"/ break;
+"/ }
+"/ }
+"/
+"/ from_string = atof(s.c_str());
+"/ // if this ratio is large, then the string has not been fixed
+"/ if( fabs( from_string / this->x[0] ) > 3.0 ){
+"/ dd_real::error("Re-rounding unsuccessful in large number fixed point trap.") ;
+"/ }
+"/ }
+"/ }
+"/
+"/ if (!!fixed) {
+"/ /* Fill in exponent part */
+"/ s += uppercase ? 'E' : 'e';
+"/ append_expn(s, e);
+"/ }
+"/ }
+
+ "/ fill in the blanks
+ (delta := width-count) > 0 ifTrue:[
+ self halt.
+"/ if (fmt & ios_base::internal) {
+"/ if (sgn)
+"/ s.insert(static_cast<string::size_type>(1), delta, fill);
+"/ else
+"/ s.insert(static_cast<string::size_type>(0), delta, fill);
+"/ } else if (fmt & ios_base::left) {
+"/ s.append(delta, fill);
+"/ } else {
+"/ s.insert(static_cast<string::size_type>(0), delta, fill);
+"/ }
+"/ }
+ ].
+
+ "Created: / 15-06-2017 / 02:37:31 / cg"
!
printString
"return a printed representation of the receiver"
- self d1 = 0.0 ifTrue:[
- ^ self d0 printString
- ].
- ^ self d0 printString , '...'
+ ^ String streamContents:[:s | self printOn:s]
"Created: / 12-06-2017 / 23:41:04 / cg"
+ "Modified: / 15-06-2017 / 01:53:46 / cg"
! !
!QDouble methodsFor:'private accessing'!
@@ -1867,10 +2328,10 @@
!
renorm
+ "destructive renormalization"
%{
double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
double c0, c1, c2, c3;
- OBJ newQD;
c0 = a[0];
c1 = a[1];
@@ -1890,6 +2351,21 @@
"
"Created: / 13-06-2017 / 18:05:33 / cg"
+ "Modified: / 15-06-2017 / 00:12:59 / cg"
+! !
+
+!QDouble methodsFor:'testing'!
+
+isInfinite
+ ^ self d0 isInfinite
+
+ "Created: / 15-06-2017 / 01:57:57 / cg"
+!
+
+isNaN
+ ^ self d0 isNaN
+
+ "Created: / 15-06-2017 / 01:57:35 / cg"
! !
!QDouble class methodsFor:'documentation'!