--- a/QDouble.st Tue Jun 20 08:56:36 2017 +0200
+++ b/QDouble.st Tue Jun 20 09:14:05 2017 +0200
@@ -1656,144 +1656,142 @@
sumFromQDouble:aQDouble
%{
if (__Class(aQDouble) == QDouble) {
- double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
- double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
- OBJ newQD;
+ double *a = __QuadDoubleInstPtr(self)->d_quadDoubleValue;
+ double *b = __QuadDoubleInstPtr(aQDouble)->d_quadDoubleValue;
+ OBJ newQD;
#ifndef QD_IEEE_ADD
- // sloppy_add...
+ // sloppy_add...
# if 0
- double s0, s1, s2, s3;
- double t0, t1, t2, t3;
- int savedCV;
-
- fpu_fix_start(&savedCV);
- m_two_sum(s0, a[0], b[0], t0);
- m_two_sum(s1, a[1], b[1], t1);
- m_two_sum(s2, a[2], b[2], t2);
- m_two_sum(s3, a[3], b[3], t3);
-
- m_two_sum(s1, s1, t0, t0);
- m_three_sum(s2, t0, t1);
- m_three_sum2(s3, t0, t2);
- t0 = t0 + t1 + t3;
-
- m_renorm5(s0, s1, s2, s3, t0);
- fpu_fix_end(&savedCV);
-
- __qNew_qdReal(newQD, s0, s1, s2, s3);
- RETURN(newQD);
+ double s0, s1, s2, s3;
+ double t0, t1, t2, t3;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ m_two_sum(s0, a[0], b[0], t0);
+ m_two_sum(s1, a[1], b[1], t1);
+ m_two_sum(s2, a[2], b[2], t2);
+ m_two_sum(s3, a[3], b[3], t3);
+
+ m_two_sum(s1, s1, t0, t0);
+ m_three_sum(s2, t0, t1);
+ m_three_sum2(s3, t0, t2);
+ t0 = t0 + t1 + t3;
+
+ m_renorm5(s0, s1, s2, s3, t0);
+ fpu_fix_end(&savedCV);
+
+ __qNew_qdReal(newQD, s0, s1, s2, s3);
# else
- /* Same as above, but addition re-organized to minimize
- data dependency ... unfortunately some compilers are
- not very smart to do this automatically */
- double s0, s1, s2, s3;
- double t0, t1, t2, t3;
- double v0, v1, v2, v3;
- double u0, u1, u2, u3;
- double w0, w1, w2, w3;
- int savedCV;
-
- fpu_fix_start(&savedCV);
- s0 = a[0] + b[0];
- s1 = a[1] + b[1];
- s2 = a[2] + b[2];
- s3 = a[3] + b[3];
-
- v0 = s0 - a[0];
- v1 = s1 - a[1];
- v2 = s2 - a[2];
- v3 = s3 - a[3];
-
- u0 = s0 - v0;
- u1 = s1 - v1;
- u2 = s2 - v2;
- u3 = s3 - v3;
-
- w0 = a[0] - u0;
- w1 = a[1] - u1;
- w2 = a[2] - u2;
- w3 = a[3] - u3;
-
- u0 = b[0] - v0;
- u1 = b[1] - v1;
- u2 = b[2] - v2;
- u3 = b[3] - v3;
-
- t0 = w0 + u0;
- t1 = w1 + u1;
- t2 = w2 + u2;
- t3 = w3 + u3;
-
- m_two_sum(s1, s1, t0, t0);
- m_three_sum(s2, t0, t1);
- m_three_sum2(s3, t0, t2);
- t0 = t0 + t1 + t3;
-
- /* renormalize */
- m_renorm5(s0, s1, s2, s3, t0);
- fpu_fix_end(&savedCV);
- __qNew_qdReal(newQD, s0, s1, s2, s3);
- RETURN(newQD);
+ /* Same as above, but addition re-organized to minimize
+ data dependency ... unfortunately some compilers are
+ not very smart to do this automatically */
+ double s0, s1, s2, s3;
+ double t0, t1, t2, t3;
+ double v0, v1, v2, v3;
+ double u0, u1, u2, u3;
+ double w0, w1, w2, w3;
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ s0 = a[0] + b[0];
+ s1 = a[1] + b[1];
+ s2 = a[2] + b[2];
+ s3 = a[3] + b[3];
+
+ v0 = s0 - a[0];
+ v1 = s1 - a[1];
+ v2 = s2 - a[2];
+ v3 = s3 - a[3];
+
+ u0 = s0 - v0;
+ u1 = s1 - v1;
+ u2 = s2 - v2;
+ u3 = s3 - v3;
+
+ w0 = a[0] - u0;
+ w1 = a[1] - u1;
+ w2 = a[2] - u2;
+ w3 = a[3] - u3;
+
+ u0 = b[0] - v0;
+ u1 = b[1] - v1;
+ u2 = b[2] - v2;
+ u3 = b[3] - v3;
+
+ t0 = w0 + u0;
+ t1 = w1 + u1;
+ t2 = w2 + u2;
+ t3 = w3 + u3;
+
+ m_two_sum(s1, s1, t0, t0);
+ m_three_sum(s2, t0, t1);
+ m_three_sum2(s3, t0, t2);
+ t0 = t0 + t1 + t3;
+
+ /* renormalize */
+ m_renorm5(s0, s1, s2, s3, t0);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, s0, s1, s2, s3);
# endif
#else
- // ieee_add...
- int i, j, k;
- double s, t;
- double u, v; /* double-length accumulator */
- double x[4] = {0.0, 0.0, 0.0, 0.0};
- int savedCV;
-
- fpu_fix_start(&savedCV);
- i = j = k = 0;
- if (abs(a[i]) > abs(b[j]))
- u = a[i++];
- else
- u = b[j++];
- if (abs(a[i]) > abs(b[j]))
- v = a[i++];
- else
- v = b[j++];
-
- u = quick_two_sum(u, v, v);
-
- while (k < 4) {
- if (i >= 4 && j >= 4) {
- x[k] = u;
- if (k < 3)
- x[++k] = v;
- break;
- }
-
- if (i >= 4)
- t = b[j++];
- else if (j >= 4)
- t = a[i++];
- else if (abs(a[i]) > abs(b[j])) {
- t = a[i++];
- } else
- t = b[j++];
-
- s = quick_three_accum(u, v, t);
-
- if (s != 0.0) {
- x[k++] = s;
- }
- }
-
- /* add the rest. */
- for (k = i; k < 4; k++)
- x[3] += a[k];
- for (k = j; k < 4; k++)
- x[3] += b[k];
-
- m_renorm4(x[0], x[1], x[2], x[3]);
- fpu_fix_end(&savedCV);
- __qNew_qdReal(newQD, x[0], x[1], x[2], x[3]);
- RETURN(newQD);
+ // ieee_add...
+ int i, j, k;
+ double s, t;
+ double u, v; /* double-length accumulator */
+ double x[4] = {0.0, 0.0, 0.0, 0.0};
+ int savedCV;
+
+ fpu_fix_start(&savedCV);
+ i = j = k = 0;
+ if (abs(a[i]) > abs(b[j]))
+ u = a[i++];
+ else
+ u = b[j++];
+ if (abs(a[i]) > abs(b[j]))
+ v = a[i++];
+ else
+ v = b[j++];
+
+ u = quick_two_sum(u, v, v);
+
+ while (k < 4) {
+ if (i >= 4 && j >= 4) {
+ x[k] = u;
+ if (k < 3)
+ x[++k] = v;
+ break;
+ }
+
+ if (i >= 4)
+ t = b[j++];
+ else if (j >= 4)
+ t = a[i++];
+ else if (abs(a[i]) > abs(b[j])) {
+ t = a[i++];
+ } else
+ t = b[j++];
+
+ s = quick_three_accum(u, v, t);
+
+ if (s != 0.0) {
+ x[k++] = s;
+ }
+ }
+
+ /* add the rest. */
+ for (k = i; k < 4; k++)
+ x[3] += a[k];
+ for (k = j; k < 4; k++)
+ x[3] += b[k];
+
+ m_renorm4(x[0], x[1], x[2], x[3]);
+ fpu_fix_end(&savedCV);
+ __qNew_qdReal(newQD, x[0], x[1], x[2], x[3]);
#endif
+ RETURN(newQD);
}
%}.
^ super sumFromQDouble:aQDouble
@@ -1811,7 +1809,7 @@
"
"Created: / 12-06-2017 / 21:15:43 / cg"
- "Modified: / 15-06-2017 / 00:49:10 / cg"
+ "Modified: / 20-06-2017 / 00:04:06 / cg"
! !
!QDouble methodsFor:'inspecting'!
@@ -1841,34 +1839,36 @@
"/ evaluated using the familiar Taylor series. Reducing the
"/ argument substantially speeds up the convergence. */
- |k inv_k d0 m mul_pwr2 r
+ |k inv_k d0 m mul_pwr2 mul2 r
s p t thresh eps i|
eps := 1.21543267145725e-63. "/ = 2^-209.
- k := 1.0 ldexp:16.
- inv_k := 1.0 / k.
-
d0 := self d0.
(d0 <= -709.0) ifTrue:[
- ^ 0.0 asQDouble.
+ ^ 0.0 asQDouble.
].
(d0 >= 709.0) ifTrue:[
- ^ Infinity positive
+ ^ Infinity positive
+ ].
+
+ (d0 = 0.0) ifTrue:[
+ ^ 1.0 asQDouble.
].
-
- self isZero ifTrue:[
- ^ 1.0 asQDouble.
+ (d0 = 1.0) ifTrue:[
+ self isOne ifTrue:[
+ ^ self class e
+ ].
].
- self isOne ifTrue:[
- ^ self class e
- ].
+ 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.
@@ -1878,28 +1878,29 @@
s := r + (mul_pwr2 value:p value:0.5).
i := 1.
[
- p := p * r.
- t := p * (self class invFact at:i).
- i := i+1.
- s := s + t.
+ p := p * r.
+ 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);
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
- s := (mul_pwr2 value:s value:2.0) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+ s := (mul2 value:s) + s squared.
+
s := s + 1.0.
^ s ldexp:m asInteger.
@@ -1909,6 +1910,7 @@
"
"Created: / 19-06-2017 / 01:49:32 / cg"
+ "Modified: / 19-06-2017 / 19:10:15 / cg"
!
floor
@@ -2006,26 +2008,27 @@
d0 := self d0.
d0 = 1.0 ifTrue:[
- self isOne ifTrue:[ ^ self class zero ].
+ self isOne ifTrue:[ ^ self class zero ].
].
d0 < 0.0 ifTrue:[
- ^ self class
- raise:#domainErrorSignal
- receiver:self
- selector:#ln
- arguments:#()
- errorString:'bad receiver in ln'
+ ^ self class
+ raise:#domainErrorSignal
+ receiver:self
+ selector:#ln
+ arguments:#()
+ errorString:'bad receiver in ln'
].
d0 = 0.0 ifTrue:[
- ^ Infinity negative
+ ^ Infinity negative
].
"/ initial approx.
x := d0 ln asQDouble.
- "/ two iterations of taylor
+ "/ three more iterations of taylor
x := x + (self * (x negated exp)) - 1.0.
x := x + (self * (x negated exp)) - 1.0.
x := x + (self * (x negated exp)) - 1.0.
+
^ x
"
@@ -2034,7 +2037,7 @@
"
"Created: / 18-06-2017 / 23:32:54 / cg"
- "Modified: / 19-06-2017 / 01:17:02 / cg"
+ "Modified (comment): / 19-06-2017 / 19:01:25 / cg"
!
negated