165 // decomposed into non-overlapping portions. Both of these
166 // uses could be eliminated. The magnitude comparison
167 // could be eliminated by extracting and comparing the
168 // exponents of a and b or just be performing a
169 // floating-point divide. Splitting a floating-point
170 // number into non-overlapping portions can be
171 // accomplished by judicious use of multiplies and
172 // additions. For details see T. J. Dekker, A Floating
173 // Point Technique for Extending the Available Precision ,
174 // Numerische Mathematik, vol. 18, 1971, pp.224-242 and
175 // subsequent work.
176
177 int ha = __HI(a); // high word of a
178 int hb = __HI(b); // high word of b
179
180 if ((ha - hb) > 0x3c00000) {
181 return a + b; // x / y > 2**60
182 }
183
184 int k = 0;
185 if (a > 0x1.0p500) { // a > 2**500
186 // scale a and b by 2**-600
187 ha -= 0x25800000;
188 hb -= 0x25800000;
189 a = a * TWO_MINUS_600;
190 b = b * TWO_MINUS_600;
191 k += 600;
192 }
193 double t1, t2;
194 if (b < 0x1.0p-500) { // b < 2**-500
195 if (b < Double.MIN_NORMAL) { // subnormal b or 0 */
196 if (b == 0.0)
197 return a;
198 t1 = 0x1.0p1022; // t1 = 2^1022
199 b *= t1;
200 a *= t1;
201 k -= 1022;
202 } else { // scale a and b by 2^600
203 ha += 0x25800000; // a *= 2^600
204 hb += 0x25800000; // b *= 2^600
205 a = a * TWO_PLUS_600;
|
165 // decomposed into non-overlapping portions. Both of these
166 // uses could be eliminated. The magnitude comparison
167 // could be eliminated by extracting and comparing the
168 // exponents of a and b or just be performing a
169 // floating-point divide. Splitting a floating-point
170 // number into non-overlapping portions can be
171 // accomplished by judicious use of multiplies and
172 // additions. For details see T. J. Dekker, A Floating
173 // Point Technique for Extending the Available Precision ,
174 // Numerische Mathematik, vol. 18, 1971, pp.224-242 and
175 // subsequent work.
176
177 int ha = __HI(a); // high word of a
178 int hb = __HI(b); // high word of b
179
180 if ((ha - hb) > 0x3c00000) {
181 return a + b; // x / y > 2**60
182 }
183
184 int k = 0;
185 if (a > 0x1.00000_ffff_ffffp500) { // a > ~2**500
186 // scale a and b by 2**-600
187 ha -= 0x25800000;
188 hb -= 0x25800000;
189 a = a * TWO_MINUS_600;
190 b = b * TWO_MINUS_600;
191 k += 600;
192 }
193 double t1, t2;
194 if (b < 0x1.0p-500) { // b < 2**-500
195 if (b < Double.MIN_NORMAL) { // subnormal b or 0 */
196 if (b == 0.0)
197 return a;
198 t1 = 0x1.0p1022; // t1 = 2^1022
199 b *= t1;
200 a *= t1;
201 k -= 1022;
202 } else { // scale a and b by 2^600
203 ha += 0x25800000; // a *= 2^600
204 hb += 0x25800000; // b *= 2^600
205 a = a * TWO_PLUS_600;
|