--- a/QDouble.st Tue Jun 13 16:53:49 2017 +0200
+++ b/QDouble.st Tue Jun 13 17:56:31 2017 +0200
@@ -21,6 +21,7 @@
!
!QDouble primitiveDefinitions!
+
%{
#include <stdio.h>
#include <errno.h>
@@ -201,7 +202,7 @@
m_quick_two_sum(s2, s2, c3, s3);\
} else {\
m_quick_two_sum(s1, s1, c3, s2);\
- }\
+ }\
} else {\
m_quick_two_sum(s0, s0, c2, s1);\
if (s1 != 0.0) {\
@@ -916,25 +917,6 @@
"
"Created: / 12-06-2017 / 23:41:39 / cg"
-!
-
-negated
- ^ self class
- d0:(self d0) negated
- d1:(self d1) negated
- d2:(self d2) negated
- d3:(self d3) negated
-
- "
- (QDouble fromFloat:1.0) negated
- ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) negated asDoubleArray
-
- (((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0))
- + ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0))) asDoubleArray
- "
-
- "Created: / 12-06-2017 / 20:14:55 / cg"
- "Modified (comment): / 12-06-2017 / 23:46:57 / cg"
! !
!QDouble methodsFor:'coercing & converting'!
@@ -1016,6 +998,24 @@
"Created: / 13-06-2017 / 01:56:53 / cg"
! !
+!QDouble methodsFor:'comparing'!
+
+< aNumber
+ "return true, if the argument, aNumber is greater than the receiver"
+
+ ^ aNumber lessFromQDouble:self
+
+ "Created: / 13-06-2017 / 16:58:53 / cg"
+!
+
+= aNumber
+ "return true, if the argument, aNumber has the same value as than the receiver"
+
+ ^ aNumber equalFromQDouble:self
+
+ "Created: / 13-06-2017 / 17:12:09 / cg"
+! !
+
!QDouble methodsFor:'double dispatching'!
differenceFromFloat:aFloat
@@ -1080,44 +1080,37 @@
"Modified (comment): / 13-06-2017 / 08:43:45 / cg"
!
-floor
-%{
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- OBJ newQD;
+lessFromQDouble:aQDouble
+ "sent when aQDouble does not know how to compare to the receiver..
+ Return true if aQDouble < self"
- double x0, x1, x2, x3;
- x1 = x2 = x3 = 0.0;
- x0 =floor(a[0]);
-
- if (x0 == a[0]) {
- x1 = floor(a[1]);
+%{
+ if (__Class(aQDouble) == QDouble) {
+ double *a = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
+ double *b = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- if (x1 == a[1]) {
- x2 = floor(a[2]);
-
- if (x2 == a[2]) {
- x3 = floor(a[3]);
- }
- }
-
- m_renorm4(x0, x1, x2, x3);
- __qNew_qdReal(newQD, x0, x1, x2, x3);
- RETURN( newQD );
+ // now compare if a < b!
+ RETURN
+ ((a[0] < b[0] ||
+ (a[0] == b[0] && (a[1] < b[1] ||
+ (a[1] == b[1] && (a[2] < b[2] ||
+ (a[2] == b[2] && a[3] < b[3])))))) ? true : false);
}
-
- __qNew_qdReal(newQD, x0, x1, x2, x3);
- RETURN( newQD );
%}.
+ ^ super lessFromQDouble:aQDouble
"
- (QDouble fromFloat:4.0) floor
- (QDouble fromFloat:4.1) floor
- (QDouble fromFloat:0.1) floor
- (0.1 + (QDouble fromFloat:1.0)) floor
- (1e20 + (QDouble fromFloat:1.0)) floor
- "
+ (1.0 + 1e-40) > 1.0
+ ((QDouble fromFloat:1.0) + (QDouble fromFloat:1e-40)) > (QDouble fromFloat:1.0)
- "Created: / 13-06-2017 / 01:52:44 / cg"
+ (QDouble fromFloat:1.0) > (QDouble fromFloat:1.0)
+ (QDouble fromFloat:1.1) > (QDouble fromFloat:1.0)
+ (QDouble fromFloat:1.0) > 1.0
+ (QDouble fromFloat:1.1) > 1.0
+ 1.0 > (QDouble fromFloat:1.0)
+ "
+
+ "Created: / 13-06-2017 / 17:07:47 / cg"
!
productFromFloat:aFloat
@@ -1221,62 +1214,49 @@
"Created: / 13-06-2017 / 01:06:22 / cg"
!
-squared
-%{
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double p0, p1, p2, p3, p4, p5;
- double q0, q1, q2, q3;
- double s0, s1;
- double t0, t1;
- double t;
- OBJ newQD;
+quotientFromQDouble:aQDouble
+ |q0 q1 q2 q3 r|
+
+ q0 := self d0 / aQDouble d0.
+ r := self - (aQDouble * q0).
+
+ q1 := r d0 / aQDouble d0.
+ r := r - (aQDouble * q1).
- m_two_sqr(p0, a[0], q0);
- t = 2.0 * a[0];
-
- m_two_prod(p1, t, a[1], q1);
- m_two_prod(p2, t, a[2], q2);
- m_two_sqr(p3, a[1], q3);
-
- m_two_sum(p1, q0, p1, q0);
-
- m_two_sum(q0, q0, q1, q1);
- m_two_sum(p2, p2, p3, p3);
-
- m_two_sum(s0, q0, p2, t0);
- m_two_sum(s1, q1, p3, t1);
+ q2 := r d0 / aQDouble d0.
+ r := r - (aQDouble * q2).
- m_two_sum(s1, s1, t0, t0);
- t0 += t1;
-
- m_quick_two_sum(s1, s1, t0, t0);
- m_quick_two_sum(p2, s0, s1, t1);
- m_quick_two_sum(p3, t1, t0, q0);
-
- p4 = 2.0 * a[0] * a[3];
- p5 = 2.0 * a[1] * a[2];
-
- m_two_sum(p4, p4, p5, p5);
- m_two_sum(q2, q2, q3, q3);
-
- m_two_sum(t0, p4, q2, t1);
- t1 = t1 + p5 + q3;
-
- m_two_sum(p3, p3, t0, p4);
- p4 = p4 + q0 + t1;
-
- m_renorm5(p0, p1, p2, p3, p4);
-
- __qNew_qdReal(newQD, p0, p1, s0, s1);
- RETURN( newQD );
+ q3 := r d0 / aQDouble d0.
+%{
+ {
+ double _q0 = __floatVal(q0);
+ double _q1 = __floatVal(q1);
+ double _q2 = __floatVal(q2);
+ double _q3 = __floatVal(q3);;
+ OBJ newQD;
+
+ m_renorm4(_q0, _q1, _q2, _q3);
+ __qNew_qdReal(newQD, _q0, _q1, _q2, _q3);
+ RETURN( newQD );
+ }
%}.
+ ^ super quotientFromQDouble:aQDouble.
"
- (QDouble fromFloat:4.0) squared
- (1e20 + (QDouble fromFloat:1.0)) squared
+ 2.0 / (QDouble fromFloat:2.0)
+ 1e20 / (QDouble fromFloat:1.0)
+ (2.0 / (QDouble fromFloat:1.0)) asFloat
+ (1e20 / (QDouble fromFloat:1.0)) asFloat
+
+ (QDouble fromFloat:2.0) / 2.0
+ (QDouble fromFloat:1e20) / 2.0
+ ((QDouble fromFloat:1.0) / 2.0) asFloat
+ ((QDouble fromFloat:1e20 / 2.0)) asFloat
+
+ ((1e20 + (QDouble fromFloat:1.0) + 1e-20) / 2.0) asDoubleArray
"
- "Created: / 13-06-2017 / 01:27:58 / cg"
+ "Created: / 13-06-2017 / 17:50:35 / cg"
!
sumFromFloat:aFloat
@@ -1480,8 +1460,166 @@
"Created: / 12-06-2017 / 23:43:05 / cg"
! !
+!QDouble methodsFor:'mathematical functions'!
+
+floor
+ "return the receiver truncated towards negative infinity"
+
+%{
+ double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
+ OBJ newQD;
+
+ double x0, x1, x2, x3;
+ x1 = x2 = x3 = 0.0;
+ x0 =floor(a[0]);
+
+ if (x0 == a[0]) {
+ x1 = floor(a[1]);
+
+ if (x1 == a[1]) {
+ x2 = floor(a[2]);
+
+ if (x2 == a[2]) {
+ x3 = floor(a[3]);
+ }
+ }
+
+ m_renorm4(x0, x1, x2, x3);
+ __qNew_qdReal(newQD, x0, x1, x2, x3);
+ RETURN( newQD );
+ }
+
+ __qNew_qdReal(newQD, x0, x1, x2, x3);
+ RETURN( newQD );
+%}.
+
+ "
+ (QDouble fromFloat:4.0) floor
+ (QDouble fromFloat:4.1) floor
+ (QDouble fromFloat:0.1) floor
+ (0.1 + (QDouble fromFloat:1.0)) floor
+ (1e20 + (QDouble fromFloat:1.0)) floor
+
+ (QDouble fromFloat:1.5) floor
+ (QDouble fromFloat:0.5) floor
+ (QDouble fromFloat:-0.5) floor
+ (QDouble fromFloat:-1.5) floor
+ "
+
+ "Created: / 13-06-2017 / 01:52:44 / cg"
+ "Modified (comment): / 13-06-2017 / 17:33:19 / cg"
+!
+
+negated
+ ^ self class
+ d0:(self d0) negated
+ d1:(self d1) negated
+ d2:(self d2) negated
+ d3:(self d3) negated
+
+ "
+ (QDouble fromFloat:1.0) negated
+ ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0)) negated asDoubleArray
+
+ (((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0))
+ + ((QDouble fromFloat:1e20) + (QDouble fromFloat:1.0))) asDoubleArray
+ "
+
+ "Created: / 12-06-2017 / 20:14:55 / cg"
+ "Modified (comment): / 12-06-2017 / 23:46:57 / cg"
+!
+
+squared
+%{
+ double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
+ double p0, p1, p2, p3, p4, p5;
+ double q0, q1, q2, q3;
+ double s0, s1;
+ double t0, t1;
+ double t;
+ OBJ newQD;
+
+ m_two_sqr(p0, a[0], q0);
+ t = 2.0 * a[0];
+
+ m_two_prod(p1, t, a[1], q1);
+ m_two_prod(p2, t, a[2], q2);
+ m_two_sqr(p3, a[1], q3);
+
+ m_two_sum(p1, q0, p1, q0);
+
+ m_two_sum(q0, q0, q1, q1);
+ m_two_sum(p2, p2, p3, p3);
+
+ m_two_sum(s0, q0, p2, t0);
+ m_two_sum(s1, q1, p3, t1);
+
+ m_two_sum(s1, s1, t0, t0);
+ t0 += t1;
+
+ m_quick_two_sum(s1, s1, t0, t0);
+ m_quick_two_sum(p2, s0, s1, t1);
+ m_quick_two_sum(p3, t1, t0, q0);
+
+ p4 = 2.0 * a[0] * a[3];
+ p5 = 2.0 * a[1] * a[2];
+
+ m_two_sum(p4, p4, p5, p5);
+ m_two_sum(q2, q2, q3, q3);
+
+ m_two_sum(t0, p4, q2, t1);
+ t1 = t1 + p5 + q3;
+
+ m_two_sum(p3, p3, t0, p4);
+ p4 = p4 + q0 + t1;
+
+ m_renorm5(p0, p1, p2, p3, p4);
+
+ __qNew_qdReal(newQD, p0, p1, s0, s1);
+ RETURN( newQD );
+%}.
+
+ "
+ (QDouble fromFloat:4.0) squared
+ (1e20 + (QDouble fromFloat:1.0)) squared
+ "
+
+ "Created: / 13-06-2017 / 01:27:58 / cg"
+! !
+
!QDouble methodsFor:'printing & storing'!
+printDigitsOn:outStream precision:precision exponentInto:expn
+ |numDigits r e i d|
+
+ numDigits := precision+1. "/ number of digits
+ r := self abs.
+ self d0 = 0.0 ifTrue:[
+ expn value:0.
+ precision timesRepeat:[ outStream nextPut:$0 ].
+ ^ 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)).
+ ] ifFalse:[
+ e > 300 ifTrue:[
+ ] ifFalse:[
+ ]
+ ].
+
+
+ "
+ 1.2345 asQDouble printDigitsOn:Transcript precision:5 exponentInto:[:e |]
+ "
+
+ "Created: / 13-06-2017 / 17:28:08 / cg"
+!
+
printString
"return a printed representation of the receiver"