147 __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[1] = d1; \ |
147 __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[1] = d1; \ |
148 __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[2] = d2; \ |
148 __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[2] = d2; \ |
149 __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[3] = d3; \ |
149 __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[3] = d3; \ |
150 } |
150 } |
151 |
151 |
|
152 // sigh: not all compilers (borland) support inline functions; |
|
153 // therefore we have to use macros... |
|
154 // sigh2: c-macros are unhygienic - to avoid catching/hiding variable bindings, |
|
155 // use different names in each macros (i.e. a_xxx) |
|
156 |
152 // qd_real(c0, c1, c2, c3); |
157 // qd_real(c0, c1, c2, c3); |
153 |
158 |
154 #define _QD_SPLITTER 134217729.0 // = 2^27 + 1 |
159 #define _QD_SPLITTER 134217729.0 // = 2^27 + 1 |
155 #define _QD_SPLIT_THRESH 6.69692879491417e+299 // = 2^996 |
160 #define _QD_SPLIT_THRESH 6.69692879491417e+299 // = 2^996 |
156 |
161 |
157 #define m_quick_two_sum(rslt, a, b, err) {\ |
162 #define m_quick_two_sum(rslt_1, a_1, b_1, err_1)\ |
158 double s = a + b;\ |
163 {\ |
159 err = b - (s - a);\ |
164 double s_1 = (a_1) + (b_1);\ |
160 rslt = s; \ |
165 (err_1) = (b_1) - (s_1 - (a_1));\ |
161 } |
166 (rslt_1) = s_1; \ |
162 |
167 } |
163 #define m_quick_two_diff(rslt, a, b, err) {\ |
168 |
164 double s = a - b;\ |
169 #define m_quick_two_diff(rslt_2, a_2, b_2, err_2)\ |
165 err = (a - s) - b;\ |
170 {\ |
166 rslt = s;\ |
171 double s_2 = (a_2) - (b_2);\ |
167 } |
172 (err_2) = ((a_2) - s_2) - (b_2);\ |
168 |
173 (rslt_2) = s_2;\ |
169 #define m_two_sum(rslt, a, b, err) {\ |
174 } |
170 double s = a + b;\ |
175 |
171 double bb = s - a;\ |
176 #define m_two_sum(rslt_3, a_3, b_3, err_3)\ |
172 err = (a - (s - bb)) + (b - bb);\ |
177 {\ |
173 rslt = s;\ |
178 double s_3 = a_3 + b_3;\ |
|
179 double bb_3 = s_3 - a_3;\ |
|
180 err_3 = (a_3 - (s_3 - bb_3)) + (b_3 - bb_3);\ |
|
181 rslt_3 = s_3;\ |
174 } |
182 } |
175 |
183 |
176 /* Computes fl(a-b) and err(a-b). */ |
184 /* Computes fl(a-b) and err(a-b). */ |
177 #define m_two_diff(rslt, a, b, err) {\ |
185 #define m_two_diff(rslt_4, a_4, b_4, err_4)\ |
178 double s = a - b;\ |
186 {\ |
179 double bb = s - a;\ |
187 double s_4 = a_4 - b_4;\ |
180 err = (a - (s - bb)) - (b + bb);\ |
188 double bb_4 = s_4 - a_4;\ |
181 rslt = s;\ |
189 err_4 = (a_4 - (s_4 - bb_4)) - (b_4 + bb_4);\ |
182 } |
190 rslt_4 = s_4;\ |
183 |
191 } |
184 #define m_three_sum(a, b, c) { \ |
192 |
185 double t1, t2, t3; \ |
193 #define m_three_sum(a_5, b_5, c_5)\ |
186 m_two_sum(t1, a, b, t2); \ |
194 { \ |
187 m_two_sum(a, c, t1, t3); \ |
195 double t1_5, t2_5, t3_5; \ |
188 m_two_sum(b, t2, t3, c); \ |
196 m_two_sum(t1_5, a_5, b_5, t2_5); \ |
189 } |
197 m_two_sum(a_5, c_5, t1_5, t3_5); \ |
190 |
198 m_two_sum(b_5, t2_5, t3_5, c_5); \ |
191 #define m_three_sum2(a, b, c) {\ |
199 } |
192 double t1, t2, t3;\ |
200 |
193 m_two_sum(t1, a, b, t2);\ |
201 #define m_three_sum2(a_6, b_6, c_6)\ |
194 m_two_sum(a, c, t1, t3);\ |
202 {\ |
195 b = t2 + t3;\ |
203 double t1_6, t2_6, t3_6;\ |
|
204 m_two_sum(t1_6, a_6, b_6, t2_6);\ |
|
205 m_two_sum(a_6, c_6, t1_6, t3_6);\ |
|
206 b_6 = t2_6 + t3_6;\ |
196 } |
207 } |
197 |
208 |
198 #ifndef QD_FMS |
209 #ifndef QD_FMS |
199 |
210 |
200 /* Computes high word and lo word of a */ |
211 /* Computes high word and lo word of a */ |
201 #define m_split(a, hi, lo) {\ |
212 #define m_split(a_7, hi_7, lo_7)\ |
202 double temp;\ |
213 {\ |
203 if (a > _QD_SPLIT_THRESH || a < -_QD_SPLIT_THRESH) {\ |
214 double temp_7;\ |
204 a *= 3.7252902984619140625e-09;\ |
215 if (a_7 > _QD_SPLIT_THRESH || a_7 < -_QD_SPLIT_THRESH) {\ |
205 temp = _QD_SPLITTER * a;\ |
216 a_7 *= 3.7252902984619140625e-09;\ |
206 hi = temp - (temp - a);\ |
217 temp_7 = _QD_SPLITTER * a_7;\ |
207 lo = a - hi;\ |
218 hi_7 = temp_7 - (temp_7 - a_7);\ |
208 hi *= 268435456.0;\ |
219 lo_7 = a_7 - hi_7;\ |
209 lo *= 268435456.0;\ |
220 hi_7 *= 268435456.0;\ |
210 } else {\ |
221 lo_7 *= 268435456.0;\ |
211 temp = _QD_SPLITTER * a;\ |
222 } else {\ |
212 hi = temp - (temp - a);\ |
223 temp_7 = _QD_SPLITTER * a_7;\ |
213 lo = a - hi;\ |
224 hi_7 = temp_7 - (temp_7 - a_7);\ |
214 }\ |
225 lo_7 = a_7 - hi_7;\ |
|
226 }\ |
215 } |
227 } |
216 |
228 |
217 #endif |
229 #endif |
218 |
230 |
219 |
231 |
220 #ifdef QD_FMS |
232 #ifdef QD_FMS |
221 |
233 |
222 /* Computes fl(a*b) and err(a*b). */ |
234 /* Computes fl(a*b) and err(a*b). */ |
223 |
235 |
224 #define m_two_prod(rslt, a, b, err) {\ |
236 #define m_two_prod(rslt_8, a_8, b_8, err_8)\ |
225 double p = a * b;\ |
237 {\ |
226 err = QD_FMS(a, b, p);\ |
238 double p_8 = a_8 * b_8;\ |
227 rslt = p; \ |
239 err_8 = QD_FMS(a_8, b_8, p_8);\ |
|
240 rslt_8 = p_8; \ |
228 } |
241 } |
229 |
242 |
230 #else |
243 #else |
231 |
244 |
232 #define m_two_prod(rslt, a, b, err) {\ |
245 #define m_two_prod(rslt_8, a_8, b_8, err_8)\ |
233 double a_hi, a_lo, b_hi, b_lo;\ |
246 {\ |
234 double p = a * b;\ |
247 double a_hi_8, a_lo_8, b_hi_8, b_lo_8;\ |
235 m_split(a, a_hi, a_lo);\ |
248 double p_8 = a_8 * b_8;\ |
236 m_split(b, b_hi, b_lo);\ |
249 m_split(a_8, a_hi_8, a_lo_8);\ |
237 err = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo;\ |
250 m_split(b_8, b_hi_8, b_lo_8);\ |
238 rslt = p; \ |
251 err_8 = ((a_hi_8 * b_hi_8 - p_8) + a_hi_8 * b_lo_8 + a_lo_8 * b_hi_8) + a_lo_8 * b_lo_8;\ |
|
252 rslt_8 = p_8; \ |
239 } |
253 } |
240 |
254 |
241 #endif |
255 #endif |
242 |
256 |
243 #ifdef QD_FMS |
257 #ifdef QD_FMS |
244 |
258 |
245 #define m_two_sqr(rslt, a, err) {\ |
259 #define m_two_sqr(rslt_9, a_9, err_9)\ |
246 double p = a * a;\ |
260 {\ |
247 err = QD_FMS(a, a, p);\ |
261 double p_9 = a_9 * a_9;\ |
248 rslt = p;\ |
262 err_9 = QD_FMS(a_9, a_9, p_9);\ |
|
263 rslt_9 = p_9;\ |
249 } |
264 } |
250 |
265 |
251 #else |
266 #else |
252 |
267 |
253 #define m_two_sqr(rslt, a, err) {\ |
268 #define m_two_sqr(rslt_9, a_9, err_9)\ |
254 double hi, lo;\ |
269 {\ |
255 double q = a * a;\ |
270 double hi_9, lo_9;\ |
256 m_split(a, hi, lo);\ |
271 double q_9 = a_9 * a_9;\ |
257 err = ((hi * hi - q) + 2.0 * hi * lo) + lo * lo;\ |
272 m_split(a_9, hi_9, lo_9);\ |
258 rslt = q;\ |
273 err_9 = ((hi_9 * hi_9 - q_9) + 2.0 * hi_9 * lo_9) + lo_9 * lo_9;\ |
|
274 rslt_9 = q_9;\ |
259 } |
275 } |
260 |
276 |
261 #endif |
277 #endif |
262 |
278 |
263 #define m_renorm4(c0, c1, c2, c3) {\ |
279 #define m_renorm4(c0_10, c1_10, c2_10, c3_10)\ |
264 double s0, s1, s2 = 0.0, s3 = 0.0;\ |
280 {\ |
|
281 double s0_10, s1_10, s2_10 = 0.0, s3_10 = 0.0;\ |
265 \ |
282 \ |
266 if (! isinf(c0)) { \ |
283 if (! isinf(c0_10)) { \ |
267 \ |
284 \ |
268 m_quick_two_sum(s0, c2, c3, c3);\ |
285 m_quick_two_sum(s0_10, c2_10, c3_10, c3_10);\ |
269 m_quick_two_sum(s0, c1, s0, c2);\ |
286 m_quick_two_sum(s0_10, c1_10, s0_10, c2_10);\ |
270 m_quick_two_sum(c0, c0, s0, c1);\ |
287 m_quick_two_sum(c0_10, c0_10, s0_10, c1_10);\ |
271 \ |
288 \ |
272 s0 = c0;\ |
289 s0_10 = c0_10;\ |
273 s1 = c1;\ |
290 s1_10 = c1_10;\ |
274 if (s1 != 0.0) {\ |
291 if (s1_10 != 0.0) {\ |
275 m_quick_two_sum(s1, s1, c2, s2);\ |
292 m_quick_two_sum(s1_10, s1_10, c2_10, s2_10);\ |
276 if (s2 != 0.0) {\ |
293 if (s2_10 != 0.0) {\ |
277 m_quick_two_sum(s2, s2, c3, s3);\ |
294 m_quick_two_sum(s2_10, s2_10, c3_10, s3_10);\ |
278 } else {\ |
295 } else {\ |
279 m_quick_two_sum(s1, s1, c3, s2);\ |
296 m_quick_two_sum(s1_10, s1_10, c3_10, s2_10);\ |
280 }\ |
297 }\ |
281 } else {\ |
298 } else {\ |
282 m_quick_two_sum(s0, s0, c2, s1);\ |
299 m_quick_two_sum(s0_10, s0_10, c2_10, s1_10);\ |
283 if (s1 != 0.0) {\ |
300 if (s1_10 != 0.0) {\ |
284 m_quick_two_sum(s1, s1, c3, s2);\ |
301 m_quick_two_sum(s1_10, s1_10, c3_10, s2_10);\ |
285 } else {\ |
302 } else {\ |
286 m_quick_two_sum(s0, s0, c3, s1);\ |
303 m_quick_two_sum(s0_10, s0_10, c3_10, s1_10);\ |
287 }\ |
304 }\ |
288 }\ |
305 }\ |
289 \ |
306 \ |
290 c0 = s0;\ |
307 c0_10 = s0_10;\ |
291 c1 = s1;\ |
308 c1_10 = s1_10;\ |
292 c2 = s2;\ |
309 c2_10 = s2_10;\ |
293 c3 = s3;\ |
310 c3_10 = s3_10;\ |
294 }\ |
311 }\ |
295 } |
312 } |
296 |
313 |
297 #define m_renorm5(c0, c1, c2, c3, c4) { \ |
314 #define m_renorm5(c0_11, c1_11, c2_11, c3_11, c4_11)\ |
298 double s0, s1, s2 = 0.0, s3 = 0.0; \ |
315 {\ |
|
316 double s0_11, s1_11, s2_11 = 0.0, s3_11 = 0.0; \ |
299 \ |
317 \ |
300 if (! isinf(c0)) { \ |
318 if (! isinf(c0_11)) { \ |
301 m_two_sum(s0, c3, c4, c4); \ |
319 m_quick_two_sum(s0_11, c3_11, c4_11, c4_11); \ |
302 m_quick_two_sum(s0, c2, s0, c3); \ |
320 m_quick_two_sum(s0_11, c2_11, s0_11, c3_11); \ |
303 m_quick_two_sum(s0, c1, s0, c2); \ |
321 m_quick_two_sum(s0_11, c1_11, s0_11, c2_11); \ |
304 m_quick_two_sum(c0, c0, s0, c1); \ |
322 m_quick_two_sum(c0_11, c0_11, s0_11, c1_11); \ |
305 \ |
323 \ |
306 s0 = c0; \ |
324 s0_11 = c0_11; \ |
307 s1 = c1; \ |
325 s1_11 = c1_11; \ |
308 \ |
326 \ |
309 m_two_sum(s0, c0, c1, s1); \ |
327 m_quick_two_sum(s0_11, c0_11, c1_11, s1_11); \ |
310 if (s1 != 0.0) { \ |
328 if (s1_11 != 0.0) { \ |
311 m_quick_two_sum(s1, s1, c2, s2); \ |
329 m_quick_two_sum(s1_11, s1_11, c2_11, s2_11); \ |
312 if (s2 != 0.0) { \ |
330 if (s2_11 != 0.0) { \ |
313 m_quick_two_sum(s2 ,s2, c3, s3); \ |
331 m_quick_two_sum(s2_11 ,s2_11, c3_11, s3_11); \ |
314 if (s3 != 0.0) {\ |
332 if (s3_11 != 0.0) {\ |
315 s3 += c4; \ |
333 s3_11 += c4_11; \ |
316 } else {\ |
334 } else {\ |
317 s2 += c4;\ |
335 s2_11 += c4_11;\ |
318 }\ |
336 }\ |
319 } else { \ |
337 } else { \ |
320 m_quick_two_sum(s1, s1, c3, s2); \ |
338 m_quick_two_sum(s1_11, s1_11, c3_11, s2_11); \ |
321 if (s2 != 0.0) {\ |
339 if (s2_11 != 0.0) {\ |
322 m_quick_two_sum(s2, s2, c4, s3); \ |
340 m_quick_two_sum(s2_11, s2_11, c4_11, s3_11); \ |
323 } else { \ |
341 } else { \ |
324 m_quick_two_sum(s1, s1, c4, s2); \ |
342 m_quick_two_sum(s1_11, s1_11, c4_11, s2_11); \ |
325 } \ |
343 } \ |
326 } \ |
344 } \ |
327 } else { \ |
345 } else { \ |
328 m_quick_two_sum(s0,s0, c2, s1); \ |
346 m_quick_two_sum(s0_11,s0_11, c2_11, s1_11); \ |
329 if (s1 != 0.0) { \ |
347 if (s1_11 != 0.0) { \ |
330 m_quick_two_sum(s1,s1, c3, s2); \ |
348 m_quick_two_sum(s1_11,s1_11, c3_11, s2_11); \ |
331 if (s2 != 0.0) {\ |
349 if (s2_11 != 0.0) {\ |
332 m_quick_two_sum(s2,s2, c4, s3); \ |
350 m_quick_two_sum(s2_11,s2_11, c4_11, s3_11); \ |
333 } else { \ |
351 } else { \ |
334 m_quick_two_sum(s1 ,s1, c4, s2); \ |
352 m_quick_two_sum(s1_11 ,s1_11, c4_11, s2_11); \ |
335 } \ |
353 } \ |
336 } else { \ |
354 } else { \ |
337 m_quick_two_sum(s0,s0, c3, s1); \ |
355 m_quick_two_sum(s0_11,s0_11, c3_11, s1_11); \ |
338 if (s1 != 0.0) { \ |
356 if (s1_11 != 0.0) { \ |
339 m_quick_two_sum(s1,s1, c4, s2); \ |
357 m_quick_two_sum(s1_11,s1_11, c4_11, s2_11); \ |
340 } else { \ |
358 } else { \ |
341 m_quick_two_sum(s0,s0, c4, s1); \ |
359 m_quick_two_sum(s0_11,s0_11, c4_11, s1_11); \ |
342 } \ |
360 } \ |
343 } \ |
361 } \ |
344 } \ |
362 } \ |
345 \ |
363 \ |
346 c0 = s0; \ |
364 c0_11 = s0_11; \ |
347 c1 = s1; \ |
365 c1_11 = s1_11; \ |
348 c2 = s2; \ |
366 c2_11 = s2_11; \ |
349 c3 = s3; \ |
367 c3_11 = s3_11; \ |
350 } \ |
368 } \ |
351 } |
369 } |
352 |
370 |
353 %} |
371 %} |
354 ! ! |
372 ! ! |