1 #include "blas_extended.h"
2 #include "blas_extended_private.h"
3
s_r_truth2(enum blas_conj_type conj,int n,float alpha,const float * x,int incx,float beta,const float * head_y,const float * tail_y,int incy,float * r,double * head_r_true,double * tail_r_true)4 void s_r_truth2(enum blas_conj_type conj, int n, float alpha,
5 const float *x, int incx, float beta,
6 const float *head_y, const float *tail_y, int incy, float *r,
7 double *head_r_true, double *tail_r_true)
8
9 /*
10 * Purpose
11 * =======
12 *
13 * This routine computes the inner product:
14 *
15 * r_true <- beta * r + alpha * SUM_{i=0, n-1} x[i] * (head_y[i] + tail_y[i])
16 *
17 * Arguments
18 * =========
19 *
20 * conj (input) enum blas_conj_type
21 * When x and y are complex vectors, specifies whether vector
22 * components x[i] are used unconjugated or conjugated.
23 *
24 * n (input) int
25 * The length of vectors x and y.
26 *
27 * alpha (input) float
28 *
29 * x (input) const float*
30 * Array of length n.
31 *
32 * incx (input) int
33 * The stride used to access components x[i].
34 *
35 * beta (input) float
36 *
37 * head_y
38 * tail_y (input) const float*
39 * Array of length n.
40 *
41 * incy (input) int
42 * The stride used to access components y[i].
43 *
44 * r (input) float*
45 *
46 * head_r_true
47 * tail_r_true (output) float*
48 * The truth.
49 *
50 */
51 {
52
53 int i, ix = 0, iy = 0;
54 float *r_i = r;
55 const float *x_i = x;
56 const float *head_y_i = head_y;
57 const float *tail_y_i = tail_y;
58 float alpha_i = alpha;
59 float beta_i = beta;
60 float x_ii;
61 float y_ii;
62 float r_v;
63 double head_prod, tail_prod;
64 double head_sum, tail_sum;
65 double head_tmp1, tail_tmp1;
66 double head_tmp2, tail_tmp2;
67 FPU_FIX_DECL;
68
69 /* Immediate return */
70 if (n < 0) {
71 *head_r_true = *tail_r_true = 0.0;
72 return;
73 }
74
75 FPU_FIX_START;
76
77 r_v = r_i[0];
78 head_sum = tail_sum = 0.0;
79
80
81 if (incx < 0)
82 ix = (-n + 1) * incx;
83 if (incy < 0)
84 iy = (-n + 1) * incy;
85
86 for (i = 0; i < n; ++i) {
87 x_ii = x_i[ix];
88 y_ii = head_y_i[iy];
89
90 head_prod = (double) x_ii *y_ii;
91 tail_prod = 0.0; /* prod = x[i] * head_y[i] */
92 {
93 /* Compute double-double = double-double + double-double. */
94 double bv;
95 double s1, s2, t1, t2;
96
97 /* Add two hi words. */
98 s1 = head_sum + head_prod;
99 bv = s1 - head_sum;
100 s2 = ((head_prod - bv) + (head_sum - (s1 - bv)));
101
102 /* Add two lo words. */
103 t1 = tail_sum + tail_prod;
104 bv = t1 - tail_sum;
105 t2 = ((tail_prod - bv) + (tail_sum - (t1 - bv)));
106
107 s2 += t1;
108
109 /* Renormalize (s1, s2) to (t1, s2) */
110 t1 = s1 + s2;
111 s2 = s2 - (t1 - s1);
112
113 t2 += s2;
114
115 /* Renormalize (t1, t2) */
116 head_sum = t1 + t2;
117 tail_sum = t2 - (head_sum - t1);
118 } /* sum = sum+prod */
119 y_ii = tail_y_i[iy];
120 head_prod = (double) x_ii *y_ii;
121 tail_prod = 0.0; /* prod = x[i] * tail_y[i] */
122 {
123 /* Compute double-double = double-double + double-double. */
124 double bv;
125 double s1, s2, t1, t2;
126
127 /* Add two hi words. */
128 s1 = head_sum + head_prod;
129 bv = s1 - head_sum;
130 s2 = ((head_prod - bv) + (head_sum - (s1 - bv)));
131
132 /* Add two lo words. */
133 t1 = tail_sum + tail_prod;
134 bv = t1 - tail_sum;
135 t2 = ((tail_prod - bv) + (tail_sum - (t1 - bv)));
136
137 s2 += t1;
138
139 /* Renormalize (s1, s2) to (t1, s2) */
140 t1 = s1 + s2;
141 s2 = s2 - (t1 - s1);
142
143 t2 += s2;
144
145 /* Renormalize (t1, t2) */
146 head_sum = t1 + t2;
147 tail_sum = t2 - (head_sum - t1);
148 } /* sum = sum+prod */
149 ix += incx;
150 iy += incy;
151 } /* endfor */
152
153
154 {
155 double dt = (double) alpha_i;
156 {
157 /* Compute double-double = double-double * double. */
158 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
159
160 con = head_sum * split;
161 a11 = con - head_sum;
162 a11 = con - a11;
163 a21 = head_sum - a11;
164 con = dt * split;
165 b1 = con - dt;
166 b1 = con - b1;
167 b2 = dt - b1;
168
169 c11 = head_sum * dt;
170 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
171
172 c2 = tail_sum * dt;
173 t1 = c11 + c2;
174 t2 = (c2 - (t1 - c11)) + c21;
175
176 head_tmp1 = t1 + t2;
177 tail_tmp1 = t2 - (head_tmp1 - t1);
178 }
179 } /* tmp1 = sum*alpha */
180 head_tmp2 = (double) r_v *beta_i;
181 tail_tmp2 = 0.0; /* tmp2 = r*beta */
182 {
183 /* Compute double-double = double-double + double-double. */
184 double bv;
185 double s1, s2, t1, t2;
186
187 /* Add two hi words. */
188 s1 = head_tmp1 + head_tmp2;
189 bv = s1 - head_tmp1;
190 s2 = ((head_tmp2 - bv) + (head_tmp1 - (s1 - bv)));
191
192 /* Add two lo words. */
193 t1 = tail_tmp1 + tail_tmp2;
194 bv = t1 - tail_tmp1;
195 t2 = ((tail_tmp2 - bv) + (tail_tmp1 - (t1 - bv)));
196
197 s2 += t1;
198
199 /* Renormalize (s1, s2) to (t1, s2) */
200 t1 = s1 + s2;
201 s2 = s2 - (t1 - s1);
202
203 t2 += s2;
204
205 /* Renormalize (t1, t2) */
206 head_tmp1 = t1 + t2;
207 tail_tmp1 = t2 - (head_tmp1 - t1);
208 } /* tmp1 = tmp1+tmp2 */
209 *head_r_true = head_tmp1;
210 *tail_r_true = tail_tmp1; /* *r = tmp1 */
211
212 FPU_FIX_STOP;
213
214 }
d_r_truth2(enum blas_conj_type conj,int n,double alpha,const double * x,int incx,double beta,const double * head_y,const double * tail_y,int incy,double * r,double * head_r_true,double * tail_r_true)215 void d_r_truth2(enum blas_conj_type conj, int n, double alpha,
216 const double *x, int incx, double beta,
217 const double *head_y, const double *tail_y, int incy,
218 double *r, double *head_r_true, double *tail_r_true)
219
220 /*
221 * Purpose
222 * =======
223 *
224 * This routine computes the inner product:
225 *
226 * r_true <- beta * r + alpha * SUM_{i=0, n-1} x[i] * (head_y[i] + tail_y[i])
227 *
228 * Arguments
229 * =========
230 *
231 * conj (input) enum blas_conj_type
232 * When x and y are complex vectors, specifies whether vector
233 * components x[i] are used unconjugated or conjugated.
234 *
235 * n (input) int
236 * The length of vectors x and y.
237 *
238 * alpha (input) double
239 *
240 * x (input) const double*
241 * Array of length n.
242 *
243 * incx (input) int
244 * The stride used to access components x[i].
245 *
246 * beta (input) double
247 *
248 * head_y
249 * tail_y (input) const double*
250 * Array of length n.
251 *
252 * incy (input) int
253 * The stride used to access components y[i].
254 *
255 * r (input) double*
256 *
257 * head_r_true
258 * tail_r_true (output) double*
259 * The truth.
260 *
261 */
262 {
263
264 int i, ix = 0, iy = 0;
265 double *r_i = r;
266 const double *x_i = x;
267 const double *head_y_i = head_y;
268 const double *tail_y_i = tail_y;
269 double alpha_i = alpha;
270 double beta_i = beta;
271 double x_ii;
272 double y_ii;
273 double r_v;
274 double head_prod, tail_prod;
275 double head_sum, tail_sum;
276 double head_tmp1, tail_tmp1;
277 double head_tmp2, tail_tmp2;
278 FPU_FIX_DECL;
279
280 /* Immediate return */
281 if (n < 0) {
282 *head_r_true = *tail_r_true = 0.0;
283 return;
284 }
285
286 FPU_FIX_START;
287
288 r_v = r_i[0];
289 head_sum = tail_sum = 0.0;
290
291
292 if (incx < 0)
293 ix = (-n + 1) * incx;
294 if (incy < 0)
295 iy = (-n + 1) * incy;
296
297 for (i = 0; i < n; ++i) {
298 x_ii = x_i[ix];
299 y_ii = head_y_i[iy];
300
301 {
302 /* Compute double_double = double * double. */
303 double a1, a2, b1, b2, con;
304
305 con = x_ii * split;
306 a1 = con - x_ii;
307 a1 = con - a1;
308 a2 = x_ii - a1;
309 con = y_ii * split;
310 b1 = con - y_ii;
311 b1 = con - b1;
312 b2 = y_ii - b1;
313
314 head_prod = x_ii * y_ii;
315 tail_prod = (((a1 * b1 - head_prod) + a1 * b2) + a2 * b1) + a2 * b2;
316 } /* prod = x[i] * head_y[i] */
317 {
318 /* Compute double-double = double-double + double-double. */
319 double bv;
320 double s1, s2, t1, t2;
321
322 /* Add two hi words. */
323 s1 = head_sum + head_prod;
324 bv = s1 - head_sum;
325 s2 = ((head_prod - bv) + (head_sum - (s1 - bv)));
326
327 /* Add two lo words. */
328 t1 = tail_sum + tail_prod;
329 bv = t1 - tail_sum;
330 t2 = ((tail_prod - bv) + (tail_sum - (t1 - bv)));
331
332 s2 += t1;
333
334 /* Renormalize (s1, s2) to (t1, s2) */
335 t1 = s1 + s2;
336 s2 = s2 - (t1 - s1);
337
338 t2 += s2;
339
340 /* Renormalize (t1, t2) */
341 head_sum = t1 + t2;
342 tail_sum = t2 - (head_sum - t1);
343 } /* sum = sum+prod */
344 y_ii = tail_y_i[iy];
345 {
346 /* Compute double_double = double * double. */
347 double a1, a2, b1, b2, con;
348
349 con = x_ii * split;
350 a1 = con - x_ii;
351 a1 = con - a1;
352 a2 = x_ii - a1;
353 con = y_ii * split;
354 b1 = con - y_ii;
355 b1 = con - b1;
356 b2 = y_ii - b1;
357
358 head_prod = x_ii * y_ii;
359 tail_prod = (((a1 * b1 - head_prod) + a1 * b2) + a2 * b1) + a2 * b2;
360 } /* prod = x[i] * tail_y[i] */
361 {
362 /* Compute double-double = double-double + double-double. */
363 double bv;
364 double s1, s2, t1, t2;
365
366 /* Add two hi words. */
367 s1 = head_sum + head_prod;
368 bv = s1 - head_sum;
369 s2 = ((head_prod - bv) + (head_sum - (s1 - bv)));
370
371 /* Add two lo words. */
372 t1 = tail_sum + tail_prod;
373 bv = t1 - tail_sum;
374 t2 = ((tail_prod - bv) + (tail_sum - (t1 - bv)));
375
376 s2 += t1;
377
378 /* Renormalize (s1, s2) to (t1, s2) */
379 t1 = s1 + s2;
380 s2 = s2 - (t1 - s1);
381
382 t2 += s2;
383
384 /* Renormalize (t1, t2) */
385 head_sum = t1 + t2;
386 tail_sum = t2 - (head_sum - t1);
387 } /* sum = sum+prod */
388 ix += incx;
389 iy += incy;
390 } /* endfor */
391
392
393 {
394 /* Compute double-double = double-double * double. */
395 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
396
397 con = head_sum * split;
398 a11 = con - head_sum;
399 a11 = con - a11;
400 a21 = head_sum - a11;
401 con = alpha_i * split;
402 b1 = con - alpha_i;
403 b1 = con - b1;
404 b2 = alpha_i - b1;
405
406 c11 = head_sum * alpha_i;
407 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
408
409 c2 = tail_sum * alpha_i;
410 t1 = c11 + c2;
411 t2 = (c2 - (t1 - c11)) + c21;
412
413 head_tmp1 = t1 + t2;
414 tail_tmp1 = t2 - (head_tmp1 - t1);
415 } /* tmp1 = sum*alpha */
416 {
417 /* Compute double_double = double * double. */
418 double a1, a2, b1, b2, con;
419
420 con = r_v * split;
421 a1 = con - r_v;
422 a1 = con - a1;
423 a2 = r_v - a1;
424 con = beta_i * split;
425 b1 = con - beta_i;
426 b1 = con - b1;
427 b2 = beta_i - b1;
428
429 head_tmp2 = r_v * beta_i;
430 tail_tmp2 = (((a1 * b1 - head_tmp2) + a1 * b2) + a2 * b1) + a2 * b2;
431 } /* tmp2 = r*beta */
432 {
433 /* Compute double-double = double-double + double-double. */
434 double bv;
435 double s1, s2, t1, t2;
436
437 /* Add two hi words. */
438 s1 = head_tmp1 + head_tmp2;
439 bv = s1 - head_tmp1;
440 s2 = ((head_tmp2 - bv) + (head_tmp1 - (s1 - bv)));
441
442 /* Add two lo words. */
443 t1 = tail_tmp1 + tail_tmp2;
444 bv = t1 - tail_tmp1;
445 t2 = ((tail_tmp2 - bv) + (tail_tmp1 - (t1 - bv)));
446
447 s2 += t1;
448
449 /* Renormalize (s1, s2) to (t1, s2) */
450 t1 = s1 + s2;
451 s2 = s2 - (t1 - s1);
452
453 t2 += s2;
454
455 /* Renormalize (t1, t2) */
456 head_tmp1 = t1 + t2;
457 tail_tmp1 = t2 - (head_tmp1 - t1);
458 } /* tmp1 = tmp1+tmp2 */
459 *head_r_true = head_tmp1;
460 *tail_r_true = tail_tmp1; /* *r = tmp1 */
461
462 FPU_FIX_STOP;
463
464 }
c_r_truth2(enum blas_conj_type conj,int n,const void * alpha,const void * x,int incx,const void * beta,const void * head_y,const void * tail_y,int incy,void * r,double * head_r_true,double * tail_r_true)465 void c_r_truth2(enum blas_conj_type conj, int n, const void *alpha,
466 const void *x, int incx, const void *beta,
467 const void *head_y, const void *tail_y, int incy, void *r,
468 double *head_r_true, double *tail_r_true)
469
470 /*
471 * Purpose
472 * =======
473 *
474 * This routine computes the inner product:
475 *
476 * r_true <- beta * r + alpha * SUM_{i=0, n-1} x[i] * (head_y[i] + tail_y[i])
477 *
478 * Arguments
479 * =========
480 *
481 * conj (input) enum blas_conj_type
482 * When x and y are complex vectors, specifies whether vector
483 * components x[i] are used unconjugated or conjugated.
484 *
485 * n (input) int
486 * The length of vectors x and y.
487 *
488 * alpha (input) const void*
489 *
490 * x (input) const void*
491 * Array of length n.
492 *
493 * incx (input) int
494 * The stride used to access components x[i].
495 *
496 * beta (input) const void*
497 *
498 * head_y
499 * tail_y (input) const void*
500 * Array of length n.
501 *
502 * incy (input) int
503 * The stride used to access components y[i].
504 *
505 * r (input) void*
506 *
507 * head_r_true
508 * tail_r_true (output) void*
509 * The truth.
510 *
511 */
512 {
513
514 int i, ix = 0, iy = 0;
515 float *r_i = (float *) r;
516 const float *x_i = (float *) x;
517 const float *head_y_i = (float *) head_y;
518 const float *tail_y_i = (float *) tail_y;
519 float *alpha_i = (float *) alpha;
520 float *beta_i = (float *) beta;
521 float x_ii[2];
522 float y_ii[2];
523 float r_v[2];
524 double head_prod[2], tail_prod[2];
525 double head_sum[2], tail_sum[2];
526 double head_tmp1[2], tail_tmp1[2];
527 double head_tmp2[2], tail_tmp2[2];
528 FPU_FIX_DECL;
529
530 /* Immediate return */
531 if (n < 0) {
532 head_r_true[0] = tail_r_true[0] = head_r_true[1] = tail_r_true[1] = 0.0;
533 return;
534 }
535
536 FPU_FIX_START;
537
538 r_v[0] = r_i[0];
539 r_v[1] = r_i[0 + 1];
540 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
541 incx *= 2;
542 incy *= 2;
543 if (incx < 0)
544 ix = (-n + 1) * incx;
545 if (incy < 0)
546 iy = (-n + 1) * incy;
547
548 if (conj == blas_conj) {
549 for (i = 0; i < n; ++i) {
550 x_ii[0] = x_i[ix];
551 x_ii[1] = x_i[ix + 1];
552 y_ii[0] = head_y_i[iy];
553 y_ii[1] = head_y_i[iy + 1];
554 x_ii[1] = -x_ii[1];
555 {
556 double head_e1, tail_e1;
557 double d1;
558 double d2;
559 /* Real part */
560 d1 = (double) x_ii[0] * y_ii[0];
561 d2 = (double) -x_ii[1] * y_ii[1];
562 {
563 /* Compute double-double = double + double. */
564 double e, t1, t2;
565
566 /* Knuth trick. */
567 t1 = d1 + d2;
568 e = t1 - d1;
569 t2 = ((d2 - e) + (d1 - (t1 - e)));
570
571 /* The result is t1 + t2, after normalization. */
572 head_e1 = t1 + t2;
573 tail_e1 = t2 - (head_e1 - t1);
574 }
575 head_prod[0] = head_e1;
576 tail_prod[0] = tail_e1;
577 /* imaginary part */
578 d1 = (double) x_ii[0] * y_ii[1];
579 d2 = (double) x_ii[1] * y_ii[0];
580 {
581 /* Compute double-double = double + double. */
582 double e, t1, t2;
583
584 /* Knuth trick. */
585 t1 = d1 + d2;
586 e = t1 - d1;
587 t2 = ((d2 - e) + (d1 - (t1 - e)));
588
589 /* The result is t1 + t2, after normalization. */
590 head_e1 = t1 + t2;
591 tail_e1 = t2 - (head_e1 - t1);
592 }
593 head_prod[1] = head_e1;
594 tail_prod[1] = tail_e1;
595 } /* prod = x[i] * head_y[i] */
596 {
597 double head_t, tail_t;
598 double head_a, tail_a;
599 double head_b, tail_b;
600 /* Real part */
601 head_a = head_sum[0];
602 tail_a = tail_sum[0];
603 head_b = head_prod[0];
604 tail_b = tail_prod[0];
605 {
606 /* Compute double-double = double-double + double-double. */
607 double bv;
608 double s1, s2, t1, t2;
609
610 /* Add two hi words. */
611 s1 = head_a + head_b;
612 bv = s1 - head_a;
613 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
614
615 /* Add two lo words. */
616 t1 = tail_a + tail_b;
617 bv = t1 - tail_a;
618 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
619
620 s2 += t1;
621
622 /* Renormalize (s1, s2) to (t1, s2) */
623 t1 = s1 + s2;
624 s2 = s2 - (t1 - s1);
625
626 t2 += s2;
627
628 /* Renormalize (t1, t2) */
629 head_t = t1 + t2;
630 tail_t = t2 - (head_t - t1);
631 }
632 head_sum[0] = head_t;
633 tail_sum[0] = tail_t;
634 /* Imaginary part */
635 head_a = head_sum[1];
636 tail_a = tail_sum[1];
637 head_b = head_prod[1];
638 tail_b = tail_prod[1];
639 {
640 /* Compute double-double = double-double + double-double. */
641 double bv;
642 double s1, s2, t1, t2;
643
644 /* Add two hi words. */
645 s1 = head_a + head_b;
646 bv = s1 - head_a;
647 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
648
649 /* Add two lo words. */
650 t1 = tail_a + tail_b;
651 bv = t1 - tail_a;
652 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
653
654 s2 += t1;
655
656 /* Renormalize (s1, s2) to (t1, s2) */
657 t1 = s1 + s2;
658 s2 = s2 - (t1 - s1);
659
660 t2 += s2;
661
662 /* Renormalize (t1, t2) */
663 head_t = t1 + t2;
664 tail_t = t2 - (head_t - t1);
665 }
666 head_sum[1] = head_t;
667 tail_sum[1] = tail_t;
668 } /* sum = sum+prod */
669 y_ii[0] = tail_y_i[iy];
670 y_ii[1] = tail_y_i[iy + 1];
671 {
672 double head_e1, tail_e1;
673 double d1;
674 double d2;
675 /* Real part */
676 d1 = (double) x_ii[0] * y_ii[0];
677 d2 = (double) -x_ii[1] * y_ii[1];
678 {
679 /* Compute double-double = double + double. */
680 double e, t1, t2;
681
682 /* Knuth trick. */
683 t1 = d1 + d2;
684 e = t1 - d1;
685 t2 = ((d2 - e) + (d1 - (t1 - e)));
686
687 /* The result is t1 + t2, after normalization. */
688 head_e1 = t1 + t2;
689 tail_e1 = t2 - (head_e1 - t1);
690 }
691 head_prod[0] = head_e1;
692 tail_prod[0] = tail_e1;
693 /* imaginary part */
694 d1 = (double) x_ii[0] * y_ii[1];
695 d2 = (double) x_ii[1] * y_ii[0];
696 {
697 /* Compute double-double = double + double. */
698 double e, t1, t2;
699
700 /* Knuth trick. */
701 t1 = d1 + d2;
702 e = t1 - d1;
703 t2 = ((d2 - e) + (d1 - (t1 - e)));
704
705 /* The result is t1 + t2, after normalization. */
706 head_e1 = t1 + t2;
707 tail_e1 = t2 - (head_e1 - t1);
708 }
709 head_prod[1] = head_e1;
710 tail_prod[1] = tail_e1;
711 } /* prod = x[i] * tail_y[i] */
712 {
713 double head_t, tail_t;
714 double head_a, tail_a;
715 double head_b, tail_b;
716 /* Real part */
717 head_a = head_sum[0];
718 tail_a = tail_sum[0];
719 head_b = head_prod[0];
720 tail_b = tail_prod[0];
721 {
722 /* Compute double-double = double-double + double-double. */
723 double bv;
724 double s1, s2, t1, t2;
725
726 /* Add two hi words. */
727 s1 = head_a + head_b;
728 bv = s1 - head_a;
729 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
730
731 /* Add two lo words. */
732 t1 = tail_a + tail_b;
733 bv = t1 - tail_a;
734 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
735
736 s2 += t1;
737
738 /* Renormalize (s1, s2) to (t1, s2) */
739 t1 = s1 + s2;
740 s2 = s2 - (t1 - s1);
741
742 t2 += s2;
743
744 /* Renormalize (t1, t2) */
745 head_t = t1 + t2;
746 tail_t = t2 - (head_t - t1);
747 }
748 head_sum[0] = head_t;
749 tail_sum[0] = tail_t;
750 /* Imaginary part */
751 head_a = head_sum[1];
752 tail_a = tail_sum[1];
753 head_b = head_prod[1];
754 tail_b = tail_prod[1];
755 {
756 /* Compute double-double = double-double + double-double. */
757 double bv;
758 double s1, s2, t1, t2;
759
760 /* Add two hi words. */
761 s1 = head_a + head_b;
762 bv = s1 - head_a;
763 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
764
765 /* Add two lo words. */
766 t1 = tail_a + tail_b;
767 bv = t1 - tail_a;
768 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
769
770 s2 += t1;
771
772 /* Renormalize (s1, s2) to (t1, s2) */
773 t1 = s1 + s2;
774 s2 = s2 - (t1 - s1);
775
776 t2 += s2;
777
778 /* Renormalize (t1, t2) */
779 head_t = t1 + t2;
780 tail_t = t2 - (head_t - t1);
781 }
782 head_sum[1] = head_t;
783 tail_sum[1] = tail_t;
784 } /* sum = sum+prod */
785 ix += incx;
786 iy += incy;
787 } /* endfor */
788 } else {
789 /* do not conjugate */
790
791 for (i = 0; i < n; ++i) {
792 x_ii[0] = x_i[ix];
793 x_ii[1] = x_i[ix + 1];
794 y_ii[0] = head_y_i[iy];
795 y_ii[1] = head_y_i[iy + 1];
796
797 {
798 double head_e1, tail_e1;
799 double d1;
800 double d2;
801 /* Real part */
802 d1 = (double) x_ii[0] * y_ii[0];
803 d2 = (double) -x_ii[1] * y_ii[1];
804 {
805 /* Compute double-double = double + double. */
806 double e, t1, t2;
807
808 /* Knuth trick. */
809 t1 = d1 + d2;
810 e = t1 - d1;
811 t2 = ((d2 - e) + (d1 - (t1 - e)));
812
813 /* The result is t1 + t2, after normalization. */
814 head_e1 = t1 + t2;
815 tail_e1 = t2 - (head_e1 - t1);
816 }
817 head_prod[0] = head_e1;
818 tail_prod[0] = tail_e1;
819 /* imaginary part */
820 d1 = (double) x_ii[0] * y_ii[1];
821 d2 = (double) x_ii[1] * y_ii[0];
822 {
823 /* Compute double-double = double + double. */
824 double e, t1, t2;
825
826 /* Knuth trick. */
827 t1 = d1 + d2;
828 e = t1 - d1;
829 t2 = ((d2 - e) + (d1 - (t1 - e)));
830
831 /* The result is t1 + t2, after normalization. */
832 head_e1 = t1 + t2;
833 tail_e1 = t2 - (head_e1 - t1);
834 }
835 head_prod[1] = head_e1;
836 tail_prod[1] = tail_e1;
837 } /* prod = x[i] * head_y[i] */
838 {
839 double head_t, tail_t;
840 double head_a, tail_a;
841 double head_b, tail_b;
842 /* Real part */
843 head_a = head_sum[0];
844 tail_a = tail_sum[0];
845 head_b = head_prod[0];
846 tail_b = tail_prod[0];
847 {
848 /* Compute double-double = double-double + double-double. */
849 double bv;
850 double s1, s2, t1, t2;
851
852 /* Add two hi words. */
853 s1 = head_a + head_b;
854 bv = s1 - head_a;
855 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
856
857 /* Add two lo words. */
858 t1 = tail_a + tail_b;
859 bv = t1 - tail_a;
860 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
861
862 s2 += t1;
863
864 /* Renormalize (s1, s2) to (t1, s2) */
865 t1 = s1 + s2;
866 s2 = s2 - (t1 - s1);
867
868 t2 += s2;
869
870 /* Renormalize (t1, t2) */
871 head_t = t1 + t2;
872 tail_t = t2 - (head_t - t1);
873 }
874 head_sum[0] = head_t;
875 tail_sum[0] = tail_t;
876 /* Imaginary part */
877 head_a = head_sum[1];
878 tail_a = tail_sum[1];
879 head_b = head_prod[1];
880 tail_b = tail_prod[1];
881 {
882 /* Compute double-double = double-double + double-double. */
883 double bv;
884 double s1, s2, t1, t2;
885
886 /* Add two hi words. */
887 s1 = head_a + head_b;
888 bv = s1 - head_a;
889 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
890
891 /* Add two lo words. */
892 t1 = tail_a + tail_b;
893 bv = t1 - tail_a;
894 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
895
896 s2 += t1;
897
898 /* Renormalize (s1, s2) to (t1, s2) */
899 t1 = s1 + s2;
900 s2 = s2 - (t1 - s1);
901
902 t2 += s2;
903
904 /* Renormalize (t1, t2) */
905 head_t = t1 + t2;
906 tail_t = t2 - (head_t - t1);
907 }
908 head_sum[1] = head_t;
909 tail_sum[1] = tail_t;
910 } /* sum = sum+prod */
911 y_ii[0] = tail_y_i[iy];
912 y_ii[1] = tail_y_i[iy + 1];
913 {
914 double head_e1, tail_e1;
915 double d1;
916 double d2;
917 /* Real part */
918 d1 = (double) x_ii[0] * y_ii[0];
919 d2 = (double) -x_ii[1] * y_ii[1];
920 {
921 /* Compute double-double = double + double. */
922 double e, t1, t2;
923
924 /* Knuth trick. */
925 t1 = d1 + d2;
926 e = t1 - d1;
927 t2 = ((d2 - e) + (d1 - (t1 - e)));
928
929 /* The result is t1 + t2, after normalization. */
930 head_e1 = t1 + t2;
931 tail_e1 = t2 - (head_e1 - t1);
932 }
933 head_prod[0] = head_e1;
934 tail_prod[0] = tail_e1;
935 /* imaginary part */
936 d1 = (double) x_ii[0] * y_ii[1];
937 d2 = (double) x_ii[1] * y_ii[0];
938 {
939 /* Compute double-double = double + double. */
940 double e, t1, t2;
941
942 /* Knuth trick. */
943 t1 = d1 + d2;
944 e = t1 - d1;
945 t2 = ((d2 - e) + (d1 - (t1 - e)));
946
947 /* The result is t1 + t2, after normalization. */
948 head_e1 = t1 + t2;
949 tail_e1 = t2 - (head_e1 - t1);
950 }
951 head_prod[1] = head_e1;
952 tail_prod[1] = tail_e1;
953 } /* prod = x[i] * tail_y[i] */
954 {
955 double head_t, tail_t;
956 double head_a, tail_a;
957 double head_b, tail_b;
958 /* Real part */
959 head_a = head_sum[0];
960 tail_a = tail_sum[0];
961 head_b = head_prod[0];
962 tail_b = tail_prod[0];
963 {
964 /* Compute double-double = double-double + double-double. */
965 double bv;
966 double s1, s2, t1, t2;
967
968 /* Add two hi words. */
969 s1 = head_a + head_b;
970 bv = s1 - head_a;
971 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
972
973 /* Add two lo words. */
974 t1 = tail_a + tail_b;
975 bv = t1 - tail_a;
976 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
977
978 s2 += t1;
979
980 /* Renormalize (s1, s2) to (t1, s2) */
981 t1 = s1 + s2;
982 s2 = s2 - (t1 - s1);
983
984 t2 += s2;
985
986 /* Renormalize (t1, t2) */
987 head_t = t1 + t2;
988 tail_t = t2 - (head_t - t1);
989 }
990 head_sum[0] = head_t;
991 tail_sum[0] = tail_t;
992 /* Imaginary part */
993 head_a = head_sum[1];
994 tail_a = tail_sum[1];
995 head_b = head_prod[1];
996 tail_b = tail_prod[1];
997 {
998 /* Compute double-double = double-double + double-double. */
999 double bv;
1000 double s1, s2, t1, t2;
1001
1002 /* Add two hi words. */
1003 s1 = head_a + head_b;
1004 bv = s1 - head_a;
1005 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1006
1007 /* Add two lo words. */
1008 t1 = tail_a + tail_b;
1009 bv = t1 - tail_a;
1010 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1011
1012 s2 += t1;
1013
1014 /* Renormalize (s1, s2) to (t1, s2) */
1015 t1 = s1 + s2;
1016 s2 = s2 - (t1 - s1);
1017
1018 t2 += s2;
1019
1020 /* Renormalize (t1, t2) */
1021 head_t = t1 + t2;
1022 tail_t = t2 - (head_t - t1);
1023 }
1024 head_sum[1] = head_t;
1025 tail_sum[1] = tail_t;
1026 } /* sum = sum+prod */
1027 ix += incx;
1028 iy += incy;
1029 } /* endfor */
1030 }
1031
1032 {
1033 double cd[2];
1034 cd[0] = (double) alpha_i[0];
1035 cd[1] = (double) alpha_i[1];
1036 {
1037 /* Compute complex-extra = complex-extra * complex-double. */
1038 double head_a0, tail_a0;
1039 double head_a1, tail_a1;
1040 double head_t1, tail_t1;
1041 double head_t2, tail_t2;
1042 head_a0 = head_sum[0];
1043 tail_a0 = tail_sum[0];
1044 head_a1 = head_sum[1];
1045 tail_a1 = tail_sum[1];
1046 /* real part */
1047 {
1048 /* Compute double-double = double-double * double. */
1049 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1050
1051 con = head_a0 * split;
1052 a11 = con - head_a0;
1053 a11 = con - a11;
1054 a21 = head_a0 - a11;
1055 con = cd[0] * split;
1056 b1 = con - cd[0];
1057 b1 = con - b1;
1058 b2 = cd[0] - b1;
1059
1060 c11 = head_a0 * cd[0];
1061 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1062
1063 c2 = tail_a0 * cd[0];
1064 t1 = c11 + c2;
1065 t2 = (c2 - (t1 - c11)) + c21;
1066
1067 head_t1 = t1 + t2;
1068 tail_t1 = t2 - (head_t1 - t1);
1069 }
1070 {
1071 /* Compute double-double = double-double * double. */
1072 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1073
1074 con = head_a1 * split;
1075 a11 = con - head_a1;
1076 a11 = con - a11;
1077 a21 = head_a1 - a11;
1078 con = cd[1] * split;
1079 b1 = con - cd[1];
1080 b1 = con - b1;
1081 b2 = cd[1] - b1;
1082
1083 c11 = head_a1 * cd[1];
1084 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1085
1086 c2 = tail_a1 * cd[1];
1087 t1 = c11 + c2;
1088 t2 = (c2 - (t1 - c11)) + c21;
1089
1090 head_t2 = t1 + t2;
1091 tail_t2 = t2 - (head_t2 - t1);
1092 }
1093 head_t2 = -head_t2;
1094 tail_t2 = -tail_t2;
1095 {
1096 /* Compute double-double = double-double + double-double. */
1097 double bv;
1098 double s1, s2, t1, t2;
1099
1100 /* Add two hi words. */
1101 s1 = head_t1 + head_t2;
1102 bv = s1 - head_t1;
1103 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1104
1105 /* Add two lo words. */
1106 t1 = tail_t1 + tail_t2;
1107 bv = t1 - tail_t1;
1108 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1109
1110 s2 += t1;
1111
1112 /* Renormalize (s1, s2) to (t1, s2) */
1113 t1 = s1 + s2;
1114 s2 = s2 - (t1 - s1);
1115
1116 t2 += s2;
1117
1118 /* Renormalize (t1, t2) */
1119 head_t1 = t1 + t2;
1120 tail_t1 = t2 - (head_t1 - t1);
1121 }
1122 head_tmp1[0] = head_t1;
1123 tail_tmp1[0] = tail_t1;
1124 /* imaginary part */
1125 {
1126 /* Compute double-double = double-double * double. */
1127 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1128
1129 con = head_a1 * split;
1130 a11 = con - head_a1;
1131 a11 = con - a11;
1132 a21 = head_a1 - a11;
1133 con = cd[0] * split;
1134 b1 = con - cd[0];
1135 b1 = con - b1;
1136 b2 = cd[0] - b1;
1137
1138 c11 = head_a1 * cd[0];
1139 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1140
1141 c2 = tail_a1 * cd[0];
1142 t1 = c11 + c2;
1143 t2 = (c2 - (t1 - c11)) + c21;
1144
1145 head_t1 = t1 + t2;
1146 tail_t1 = t2 - (head_t1 - t1);
1147 }
1148 {
1149 /* Compute double-double = double-double * double. */
1150 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1151
1152 con = head_a0 * split;
1153 a11 = con - head_a0;
1154 a11 = con - a11;
1155 a21 = head_a0 - a11;
1156 con = cd[1] * split;
1157 b1 = con - cd[1];
1158 b1 = con - b1;
1159 b2 = cd[1] - b1;
1160
1161 c11 = head_a0 * cd[1];
1162 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1163
1164 c2 = tail_a0 * cd[1];
1165 t1 = c11 + c2;
1166 t2 = (c2 - (t1 - c11)) + c21;
1167
1168 head_t2 = t1 + t2;
1169 tail_t2 = t2 - (head_t2 - t1);
1170 }
1171 {
1172 /* Compute double-double = double-double + double-double. */
1173 double bv;
1174 double s1, s2, t1, t2;
1175
1176 /* Add two hi words. */
1177 s1 = head_t1 + head_t2;
1178 bv = s1 - head_t1;
1179 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1180
1181 /* Add two lo words. */
1182 t1 = tail_t1 + tail_t2;
1183 bv = t1 - tail_t1;
1184 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1185
1186 s2 += t1;
1187
1188 /* Renormalize (s1, s2) to (t1, s2) */
1189 t1 = s1 + s2;
1190 s2 = s2 - (t1 - s1);
1191
1192 t2 += s2;
1193
1194 /* Renormalize (t1, t2) */
1195 head_t1 = t1 + t2;
1196 tail_t1 = t2 - (head_t1 - t1);
1197 }
1198 head_tmp1[1] = head_t1;
1199 tail_tmp1[1] = tail_t1;
1200 }
1201
1202 } /* tmp1 = sum*alpha */
1203 {
1204 double head_e1, tail_e1;
1205 double d1;
1206 double d2;
1207 /* Real part */
1208 d1 = (double) r_v[0] * beta_i[0];
1209 d2 = (double) -r_v[1] * beta_i[1];
1210 {
1211 /* Compute double-double = double + double. */
1212 double e, t1, t2;
1213
1214 /* Knuth trick. */
1215 t1 = d1 + d2;
1216 e = t1 - d1;
1217 t2 = ((d2 - e) + (d1 - (t1 - e)));
1218
1219 /* The result is t1 + t2, after normalization. */
1220 head_e1 = t1 + t2;
1221 tail_e1 = t2 - (head_e1 - t1);
1222 }
1223 head_tmp2[0] = head_e1;
1224 tail_tmp2[0] = tail_e1;
1225 /* imaginary part */
1226 d1 = (double) r_v[0] * beta_i[1];
1227 d2 = (double) r_v[1] * beta_i[0];
1228 {
1229 /* Compute double-double = double + double. */
1230 double e, t1, t2;
1231
1232 /* Knuth trick. */
1233 t1 = d1 + d2;
1234 e = t1 - d1;
1235 t2 = ((d2 - e) + (d1 - (t1 - e)));
1236
1237 /* The result is t1 + t2, after normalization. */
1238 head_e1 = t1 + t2;
1239 tail_e1 = t2 - (head_e1 - t1);
1240 }
1241 head_tmp2[1] = head_e1;
1242 tail_tmp2[1] = tail_e1;
1243 } /* tmp2 = r*beta */
1244 {
1245 double head_t, tail_t;
1246 double head_a, tail_a;
1247 double head_b, tail_b;
1248 /* Real part */
1249 head_a = head_tmp1[0];
1250 tail_a = tail_tmp1[0];
1251 head_b = head_tmp2[0];
1252 tail_b = tail_tmp2[0];
1253 {
1254 /* Compute double-double = double-double + double-double. */
1255 double bv;
1256 double s1, s2, t1, t2;
1257
1258 /* Add two hi words. */
1259 s1 = head_a + head_b;
1260 bv = s1 - head_a;
1261 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1262
1263 /* Add two lo words. */
1264 t1 = tail_a + tail_b;
1265 bv = t1 - tail_a;
1266 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1267
1268 s2 += t1;
1269
1270 /* Renormalize (s1, s2) to (t1, s2) */
1271 t1 = s1 + s2;
1272 s2 = s2 - (t1 - s1);
1273
1274 t2 += s2;
1275
1276 /* Renormalize (t1, t2) */
1277 head_t = t1 + t2;
1278 tail_t = t2 - (head_t - t1);
1279 }
1280 head_tmp1[0] = head_t;
1281 tail_tmp1[0] = tail_t;
1282 /* Imaginary part */
1283 head_a = head_tmp1[1];
1284 tail_a = tail_tmp1[1];
1285 head_b = head_tmp2[1];
1286 tail_b = tail_tmp2[1];
1287 {
1288 /* Compute double-double = double-double + double-double. */
1289 double bv;
1290 double s1, s2, t1, t2;
1291
1292 /* Add two hi words. */
1293 s1 = head_a + head_b;
1294 bv = s1 - head_a;
1295 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1296
1297 /* Add two lo words. */
1298 t1 = tail_a + tail_b;
1299 bv = t1 - tail_a;
1300 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1301
1302 s2 += t1;
1303
1304 /* Renormalize (s1, s2) to (t1, s2) */
1305 t1 = s1 + s2;
1306 s2 = s2 - (t1 - s1);
1307
1308 t2 += s2;
1309
1310 /* Renormalize (t1, t2) */
1311 head_t = t1 + t2;
1312 tail_t = t2 - (head_t - t1);
1313 }
1314 head_tmp1[1] = head_t;
1315 tail_tmp1[1] = tail_t;
1316 } /* tmp1 = tmp1+tmp2 */
1317 head_r_true[0] = head_tmp1[0];
1318 tail_r_true[0] = tail_tmp1[0];
1319 head_r_true[1] = head_tmp1[1];
1320 tail_r_true[1] = tail_tmp1[1]; /* *r = tmp1 */
1321
1322 FPU_FIX_STOP;
1323
1324 }
z_r_truth2(enum blas_conj_type conj,int n,const void * alpha,const void * x,int incx,const void * beta,const void * head_y,const void * tail_y,int incy,void * r,double * head_r_true,double * tail_r_true)1325 void z_r_truth2(enum blas_conj_type conj, int n, const void *alpha,
1326 const void *x, int incx, const void *beta,
1327 const void *head_y, const void *tail_y, int incy, void *r,
1328 double *head_r_true, double *tail_r_true)
1329
1330 /*
1331 * Purpose
1332 * =======
1333 *
1334 * This routine computes the inner product:
1335 *
1336 * r_true <- beta * r + alpha * SUM_{i=0, n-1} x[i] * (head_y[i] + tail_y[i])
1337 *
1338 * Arguments
1339 * =========
1340 *
1341 * conj (input) enum blas_conj_type
1342 * When x and y are complex vectors, specifies whether vector
1343 * components x[i] are used unconjugated or conjugated.
1344 *
1345 * n (input) int
1346 * The length of vectors x and y.
1347 *
1348 * alpha (input) const void*
1349 *
1350 * x (input) const void*
1351 * Array of length n.
1352 *
1353 * incx (input) int
1354 * The stride used to access components x[i].
1355 *
1356 * beta (input) const void*
1357 *
1358 * head_y
1359 * tail_y (input) const void*
1360 * Array of length n.
1361 *
1362 * incy (input) int
1363 * The stride used to access components y[i].
1364 *
1365 * r (input) void*
1366 *
1367 * head_r_true
1368 * tail_r_true (output) void*
1369 * The truth.
1370 *
1371 */
1372 {
1373
1374 int i, ix = 0, iy = 0;
1375 double *r_i = (double *) r;
1376 const double *x_i = (double *) x;
1377 const double *head_y_i = (double *) head_y;
1378 const double *tail_y_i = (double *) tail_y;
1379 double *alpha_i = (double *) alpha;
1380 double *beta_i = (double *) beta;
1381 double x_ii[2];
1382 double y_ii[2];
1383 double r_v[2];
1384 double head_prod[2], tail_prod[2];
1385 double head_sum[2], tail_sum[2];
1386 double head_tmp1[2], tail_tmp1[2];
1387 double head_tmp2[2], tail_tmp2[2];
1388 FPU_FIX_DECL;
1389
1390 /* Immediate return */
1391 if (n < 0) {
1392 head_r_true[0] = tail_r_true[0] = head_r_true[1] = tail_r_true[1] = 0.0;
1393 return;
1394 }
1395
1396 FPU_FIX_START;
1397
1398 r_v[0] = r_i[0];
1399 r_v[1] = r_i[0 + 1];
1400 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
1401 incx *= 2;
1402 incy *= 2;
1403 if (incx < 0)
1404 ix = (-n + 1) * incx;
1405 if (incy < 0)
1406 iy = (-n + 1) * incy;
1407
1408 if (conj == blas_conj) {
1409 for (i = 0; i < n; ++i) {
1410 x_ii[0] = x_i[ix];
1411 x_ii[1] = x_i[ix + 1];
1412 y_ii[0] = head_y_i[iy];
1413 y_ii[1] = head_y_i[iy + 1];
1414 x_ii[1] = -x_ii[1];
1415 {
1416 /* Compute complex-extra = complex-double * complex-double. */
1417 double head_t1, tail_t1;
1418 double head_t2, tail_t2;
1419 /* Real part */
1420 {
1421 /* Compute double_double = double * double. */
1422 double a1, a2, b1, b2, con;
1423
1424 con = x_ii[0] * split;
1425 a1 = con - x_ii[0];
1426 a1 = con - a1;
1427 a2 = x_ii[0] - a1;
1428 con = y_ii[0] * split;
1429 b1 = con - y_ii[0];
1430 b1 = con - b1;
1431 b2 = y_ii[0] - b1;
1432
1433 head_t1 = x_ii[0] * y_ii[0];
1434 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
1435 }
1436 {
1437 /* Compute double_double = double * double. */
1438 double a1, a2, b1, b2, con;
1439
1440 con = x_ii[1] * split;
1441 a1 = con - x_ii[1];
1442 a1 = con - a1;
1443 a2 = x_ii[1] - a1;
1444 con = y_ii[1] * split;
1445 b1 = con - y_ii[1];
1446 b1 = con - b1;
1447 b2 = y_ii[1] - b1;
1448
1449 head_t2 = x_ii[1] * y_ii[1];
1450 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
1451 }
1452 head_t2 = -head_t2;
1453 tail_t2 = -tail_t2;
1454 {
1455 /* Compute double-double = double-double + double-double. */
1456 double bv;
1457 double s1, s2, t1, t2;
1458
1459 /* Add two hi words. */
1460 s1 = head_t1 + head_t2;
1461 bv = s1 - head_t1;
1462 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1463
1464 /* Add two lo words. */
1465 t1 = tail_t1 + tail_t2;
1466 bv = t1 - tail_t1;
1467 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1468
1469 s2 += t1;
1470
1471 /* Renormalize (s1, s2) to (t1, s2) */
1472 t1 = s1 + s2;
1473 s2 = s2 - (t1 - s1);
1474
1475 t2 += s2;
1476
1477 /* Renormalize (t1, t2) */
1478 head_t1 = t1 + t2;
1479 tail_t1 = t2 - (head_t1 - t1);
1480 }
1481 head_prod[0] = head_t1;
1482 tail_prod[0] = tail_t1;
1483 /* Imaginary part */
1484 {
1485 /* Compute double_double = double * double. */
1486 double a1, a2, b1, b2, con;
1487
1488 con = x_ii[1] * split;
1489 a1 = con - x_ii[1];
1490 a1 = con - a1;
1491 a2 = x_ii[1] - a1;
1492 con = y_ii[0] * split;
1493 b1 = con - y_ii[0];
1494 b1 = con - b1;
1495 b2 = y_ii[0] - b1;
1496
1497 head_t1 = x_ii[1] * y_ii[0];
1498 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
1499 }
1500 {
1501 /* Compute double_double = double * double. */
1502 double a1, a2, b1, b2, con;
1503
1504 con = x_ii[0] * split;
1505 a1 = con - x_ii[0];
1506 a1 = con - a1;
1507 a2 = x_ii[0] - a1;
1508 con = y_ii[1] * split;
1509 b1 = con - y_ii[1];
1510 b1 = con - b1;
1511 b2 = y_ii[1] - b1;
1512
1513 head_t2 = x_ii[0] * y_ii[1];
1514 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
1515 }
1516 {
1517 /* Compute double-double = double-double + double-double. */
1518 double bv;
1519 double s1, s2, t1, t2;
1520
1521 /* Add two hi words. */
1522 s1 = head_t1 + head_t2;
1523 bv = s1 - head_t1;
1524 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1525
1526 /* Add two lo words. */
1527 t1 = tail_t1 + tail_t2;
1528 bv = t1 - tail_t1;
1529 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1530
1531 s2 += t1;
1532
1533 /* Renormalize (s1, s2) to (t1, s2) */
1534 t1 = s1 + s2;
1535 s2 = s2 - (t1 - s1);
1536
1537 t2 += s2;
1538
1539 /* Renormalize (t1, t2) */
1540 head_t1 = t1 + t2;
1541 tail_t1 = t2 - (head_t1 - t1);
1542 }
1543 head_prod[1] = head_t1;
1544 tail_prod[1] = tail_t1;
1545 } /* prod = x[i] * head_y[i] */
1546 {
1547 double head_t, tail_t;
1548 double head_a, tail_a;
1549 double head_b, tail_b;
1550 /* Real part */
1551 head_a = head_sum[0];
1552 tail_a = tail_sum[0];
1553 head_b = head_prod[0];
1554 tail_b = tail_prod[0];
1555 {
1556 /* Compute double-double = double-double + double-double. */
1557 double bv;
1558 double s1, s2, t1, t2;
1559
1560 /* Add two hi words. */
1561 s1 = head_a + head_b;
1562 bv = s1 - head_a;
1563 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1564
1565 /* Add two lo words. */
1566 t1 = tail_a + tail_b;
1567 bv = t1 - tail_a;
1568 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1569
1570 s2 += t1;
1571
1572 /* Renormalize (s1, s2) to (t1, s2) */
1573 t1 = s1 + s2;
1574 s2 = s2 - (t1 - s1);
1575
1576 t2 += s2;
1577
1578 /* Renormalize (t1, t2) */
1579 head_t = t1 + t2;
1580 tail_t = t2 - (head_t - t1);
1581 }
1582 head_sum[0] = head_t;
1583 tail_sum[0] = tail_t;
1584 /* Imaginary part */
1585 head_a = head_sum[1];
1586 tail_a = tail_sum[1];
1587 head_b = head_prod[1];
1588 tail_b = tail_prod[1];
1589 {
1590 /* Compute double-double = double-double + double-double. */
1591 double bv;
1592 double s1, s2, t1, t2;
1593
1594 /* Add two hi words. */
1595 s1 = head_a + head_b;
1596 bv = s1 - head_a;
1597 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1598
1599 /* Add two lo words. */
1600 t1 = tail_a + tail_b;
1601 bv = t1 - tail_a;
1602 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1603
1604 s2 += t1;
1605
1606 /* Renormalize (s1, s2) to (t1, s2) */
1607 t1 = s1 + s2;
1608 s2 = s2 - (t1 - s1);
1609
1610 t2 += s2;
1611
1612 /* Renormalize (t1, t2) */
1613 head_t = t1 + t2;
1614 tail_t = t2 - (head_t - t1);
1615 }
1616 head_sum[1] = head_t;
1617 tail_sum[1] = tail_t;
1618 } /* sum = sum+prod */
1619 y_ii[0] = tail_y_i[iy];
1620 y_ii[1] = tail_y_i[iy + 1];
1621 {
1622 /* Compute complex-extra = complex-double * complex-double. */
1623 double head_t1, tail_t1;
1624 double head_t2, tail_t2;
1625 /* Real part */
1626 {
1627 /* Compute double_double = double * double. */
1628 double a1, a2, b1, b2, con;
1629
1630 con = x_ii[0] * split;
1631 a1 = con - x_ii[0];
1632 a1 = con - a1;
1633 a2 = x_ii[0] - a1;
1634 con = y_ii[0] * split;
1635 b1 = con - y_ii[0];
1636 b1 = con - b1;
1637 b2 = y_ii[0] - b1;
1638
1639 head_t1 = x_ii[0] * y_ii[0];
1640 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
1641 }
1642 {
1643 /* Compute double_double = double * double. */
1644 double a1, a2, b1, b2, con;
1645
1646 con = x_ii[1] * split;
1647 a1 = con - x_ii[1];
1648 a1 = con - a1;
1649 a2 = x_ii[1] - a1;
1650 con = y_ii[1] * split;
1651 b1 = con - y_ii[1];
1652 b1 = con - b1;
1653 b2 = y_ii[1] - b1;
1654
1655 head_t2 = x_ii[1] * y_ii[1];
1656 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
1657 }
1658 head_t2 = -head_t2;
1659 tail_t2 = -tail_t2;
1660 {
1661 /* Compute double-double = double-double + double-double. */
1662 double bv;
1663 double s1, s2, t1, t2;
1664
1665 /* Add two hi words. */
1666 s1 = head_t1 + head_t2;
1667 bv = s1 - head_t1;
1668 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1669
1670 /* Add two lo words. */
1671 t1 = tail_t1 + tail_t2;
1672 bv = t1 - tail_t1;
1673 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1674
1675 s2 += t1;
1676
1677 /* Renormalize (s1, s2) to (t1, s2) */
1678 t1 = s1 + s2;
1679 s2 = s2 - (t1 - s1);
1680
1681 t2 += s2;
1682
1683 /* Renormalize (t1, t2) */
1684 head_t1 = t1 + t2;
1685 tail_t1 = t2 - (head_t1 - t1);
1686 }
1687 head_prod[0] = head_t1;
1688 tail_prod[0] = tail_t1;
1689 /* Imaginary part */
1690 {
1691 /* Compute double_double = double * double. */
1692 double a1, a2, b1, b2, con;
1693
1694 con = x_ii[1] * split;
1695 a1 = con - x_ii[1];
1696 a1 = con - a1;
1697 a2 = x_ii[1] - a1;
1698 con = y_ii[0] * split;
1699 b1 = con - y_ii[0];
1700 b1 = con - b1;
1701 b2 = y_ii[0] - b1;
1702
1703 head_t1 = x_ii[1] * y_ii[0];
1704 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
1705 }
1706 {
1707 /* Compute double_double = double * double. */
1708 double a1, a2, b1, b2, con;
1709
1710 con = x_ii[0] * split;
1711 a1 = con - x_ii[0];
1712 a1 = con - a1;
1713 a2 = x_ii[0] - a1;
1714 con = y_ii[1] * split;
1715 b1 = con - y_ii[1];
1716 b1 = con - b1;
1717 b2 = y_ii[1] - b1;
1718
1719 head_t2 = x_ii[0] * y_ii[1];
1720 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
1721 }
1722 {
1723 /* Compute double-double = double-double + double-double. */
1724 double bv;
1725 double s1, s2, t1, t2;
1726
1727 /* Add two hi words. */
1728 s1 = head_t1 + head_t2;
1729 bv = s1 - head_t1;
1730 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1731
1732 /* Add two lo words. */
1733 t1 = tail_t1 + tail_t2;
1734 bv = t1 - tail_t1;
1735 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1736
1737 s2 += t1;
1738
1739 /* Renormalize (s1, s2) to (t1, s2) */
1740 t1 = s1 + s2;
1741 s2 = s2 - (t1 - s1);
1742
1743 t2 += s2;
1744
1745 /* Renormalize (t1, t2) */
1746 head_t1 = t1 + t2;
1747 tail_t1 = t2 - (head_t1 - t1);
1748 }
1749 head_prod[1] = head_t1;
1750 tail_prod[1] = tail_t1;
1751 } /* prod = x[i] * tail_y[i] */
1752 {
1753 double head_t, tail_t;
1754 double head_a, tail_a;
1755 double head_b, tail_b;
1756 /* Real part */
1757 head_a = head_sum[0];
1758 tail_a = tail_sum[0];
1759 head_b = head_prod[0];
1760 tail_b = tail_prod[0];
1761 {
1762 /* Compute double-double = double-double + double-double. */
1763 double bv;
1764 double s1, s2, t1, t2;
1765
1766 /* Add two hi words. */
1767 s1 = head_a + head_b;
1768 bv = s1 - head_a;
1769 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1770
1771 /* Add two lo words. */
1772 t1 = tail_a + tail_b;
1773 bv = t1 - tail_a;
1774 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1775
1776 s2 += t1;
1777
1778 /* Renormalize (s1, s2) to (t1, s2) */
1779 t1 = s1 + s2;
1780 s2 = s2 - (t1 - s1);
1781
1782 t2 += s2;
1783
1784 /* Renormalize (t1, t2) */
1785 head_t = t1 + t2;
1786 tail_t = t2 - (head_t - t1);
1787 }
1788 head_sum[0] = head_t;
1789 tail_sum[0] = tail_t;
1790 /* Imaginary part */
1791 head_a = head_sum[1];
1792 tail_a = tail_sum[1];
1793 head_b = head_prod[1];
1794 tail_b = tail_prod[1];
1795 {
1796 /* Compute double-double = double-double + double-double. */
1797 double bv;
1798 double s1, s2, t1, t2;
1799
1800 /* Add two hi words. */
1801 s1 = head_a + head_b;
1802 bv = s1 - head_a;
1803 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1804
1805 /* Add two lo words. */
1806 t1 = tail_a + tail_b;
1807 bv = t1 - tail_a;
1808 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1809
1810 s2 += t1;
1811
1812 /* Renormalize (s1, s2) to (t1, s2) */
1813 t1 = s1 + s2;
1814 s2 = s2 - (t1 - s1);
1815
1816 t2 += s2;
1817
1818 /* Renormalize (t1, t2) */
1819 head_t = t1 + t2;
1820 tail_t = t2 - (head_t - t1);
1821 }
1822 head_sum[1] = head_t;
1823 tail_sum[1] = tail_t;
1824 } /* sum = sum+prod */
1825 ix += incx;
1826 iy += incy;
1827 } /* endfor */
1828 } else {
1829 /* do not conjugate */
1830
1831 for (i = 0; i < n; ++i) {
1832 x_ii[0] = x_i[ix];
1833 x_ii[1] = x_i[ix + 1];
1834 y_ii[0] = head_y_i[iy];
1835 y_ii[1] = head_y_i[iy + 1];
1836
1837 {
1838 /* Compute complex-extra = complex-double * complex-double. */
1839 double head_t1, tail_t1;
1840 double head_t2, tail_t2;
1841 /* Real part */
1842 {
1843 /* Compute double_double = double * double. */
1844 double a1, a2, b1, b2, con;
1845
1846 con = x_ii[0] * split;
1847 a1 = con - x_ii[0];
1848 a1 = con - a1;
1849 a2 = x_ii[0] - a1;
1850 con = y_ii[0] * split;
1851 b1 = con - y_ii[0];
1852 b1 = con - b1;
1853 b2 = y_ii[0] - b1;
1854
1855 head_t1 = x_ii[0] * y_ii[0];
1856 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
1857 }
1858 {
1859 /* Compute double_double = double * double. */
1860 double a1, a2, b1, b2, con;
1861
1862 con = x_ii[1] * split;
1863 a1 = con - x_ii[1];
1864 a1 = con - a1;
1865 a2 = x_ii[1] - a1;
1866 con = y_ii[1] * split;
1867 b1 = con - y_ii[1];
1868 b1 = con - b1;
1869 b2 = y_ii[1] - b1;
1870
1871 head_t2 = x_ii[1] * y_ii[1];
1872 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
1873 }
1874 head_t2 = -head_t2;
1875 tail_t2 = -tail_t2;
1876 {
1877 /* Compute double-double = double-double + double-double. */
1878 double bv;
1879 double s1, s2, t1, t2;
1880
1881 /* Add two hi words. */
1882 s1 = head_t1 + head_t2;
1883 bv = s1 - head_t1;
1884 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1885
1886 /* Add two lo words. */
1887 t1 = tail_t1 + tail_t2;
1888 bv = t1 - tail_t1;
1889 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1890
1891 s2 += t1;
1892
1893 /* Renormalize (s1, s2) to (t1, s2) */
1894 t1 = s1 + s2;
1895 s2 = s2 - (t1 - s1);
1896
1897 t2 += s2;
1898
1899 /* Renormalize (t1, t2) */
1900 head_t1 = t1 + t2;
1901 tail_t1 = t2 - (head_t1 - t1);
1902 }
1903 head_prod[0] = head_t1;
1904 tail_prod[0] = tail_t1;
1905 /* Imaginary part */
1906 {
1907 /* Compute double_double = double * double. */
1908 double a1, a2, b1, b2, con;
1909
1910 con = x_ii[1] * split;
1911 a1 = con - x_ii[1];
1912 a1 = con - a1;
1913 a2 = x_ii[1] - a1;
1914 con = y_ii[0] * split;
1915 b1 = con - y_ii[0];
1916 b1 = con - b1;
1917 b2 = y_ii[0] - b1;
1918
1919 head_t1 = x_ii[1] * y_ii[0];
1920 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
1921 }
1922 {
1923 /* Compute double_double = double * double. */
1924 double a1, a2, b1, b2, con;
1925
1926 con = x_ii[0] * split;
1927 a1 = con - x_ii[0];
1928 a1 = con - a1;
1929 a2 = x_ii[0] - a1;
1930 con = y_ii[1] * split;
1931 b1 = con - y_ii[1];
1932 b1 = con - b1;
1933 b2 = y_ii[1] - b1;
1934
1935 head_t2 = x_ii[0] * y_ii[1];
1936 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
1937 }
1938 {
1939 /* Compute double-double = double-double + double-double. */
1940 double bv;
1941 double s1, s2, t1, t2;
1942
1943 /* Add two hi words. */
1944 s1 = head_t1 + head_t2;
1945 bv = s1 - head_t1;
1946 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1947
1948 /* Add two lo words. */
1949 t1 = tail_t1 + tail_t2;
1950 bv = t1 - tail_t1;
1951 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1952
1953 s2 += t1;
1954
1955 /* Renormalize (s1, s2) to (t1, s2) */
1956 t1 = s1 + s2;
1957 s2 = s2 - (t1 - s1);
1958
1959 t2 += s2;
1960
1961 /* Renormalize (t1, t2) */
1962 head_t1 = t1 + t2;
1963 tail_t1 = t2 - (head_t1 - t1);
1964 }
1965 head_prod[1] = head_t1;
1966 tail_prod[1] = tail_t1;
1967 } /* prod = x[i] * head_y[i] */
1968 {
1969 double head_t, tail_t;
1970 double head_a, tail_a;
1971 double head_b, tail_b;
1972 /* Real part */
1973 head_a = head_sum[0];
1974 tail_a = tail_sum[0];
1975 head_b = head_prod[0];
1976 tail_b = tail_prod[0];
1977 {
1978 /* Compute double-double = double-double + double-double. */
1979 double bv;
1980 double s1, s2, t1, t2;
1981
1982 /* Add two hi words. */
1983 s1 = head_a + head_b;
1984 bv = s1 - head_a;
1985 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1986
1987 /* Add two lo words. */
1988 t1 = tail_a + tail_b;
1989 bv = t1 - tail_a;
1990 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1991
1992 s2 += t1;
1993
1994 /* Renormalize (s1, s2) to (t1, s2) */
1995 t1 = s1 + s2;
1996 s2 = s2 - (t1 - s1);
1997
1998 t2 += s2;
1999
2000 /* Renormalize (t1, t2) */
2001 head_t = t1 + t2;
2002 tail_t = t2 - (head_t - t1);
2003 }
2004 head_sum[0] = head_t;
2005 tail_sum[0] = tail_t;
2006 /* Imaginary part */
2007 head_a = head_sum[1];
2008 tail_a = tail_sum[1];
2009 head_b = head_prod[1];
2010 tail_b = tail_prod[1];
2011 {
2012 /* Compute double-double = double-double + double-double. */
2013 double bv;
2014 double s1, s2, t1, t2;
2015
2016 /* Add two hi words. */
2017 s1 = head_a + head_b;
2018 bv = s1 - head_a;
2019 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2020
2021 /* Add two lo words. */
2022 t1 = tail_a + tail_b;
2023 bv = t1 - tail_a;
2024 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2025
2026 s2 += t1;
2027
2028 /* Renormalize (s1, s2) to (t1, s2) */
2029 t1 = s1 + s2;
2030 s2 = s2 - (t1 - s1);
2031
2032 t2 += s2;
2033
2034 /* Renormalize (t1, t2) */
2035 head_t = t1 + t2;
2036 tail_t = t2 - (head_t - t1);
2037 }
2038 head_sum[1] = head_t;
2039 tail_sum[1] = tail_t;
2040 } /* sum = sum+prod */
2041 y_ii[0] = tail_y_i[iy];
2042 y_ii[1] = tail_y_i[iy + 1];
2043 {
2044 /* Compute complex-extra = complex-double * complex-double. */
2045 double head_t1, tail_t1;
2046 double head_t2, tail_t2;
2047 /* Real part */
2048 {
2049 /* Compute double_double = double * double. */
2050 double a1, a2, b1, b2, con;
2051
2052 con = x_ii[0] * split;
2053 a1 = con - x_ii[0];
2054 a1 = con - a1;
2055 a2 = x_ii[0] - a1;
2056 con = y_ii[0] * split;
2057 b1 = con - y_ii[0];
2058 b1 = con - b1;
2059 b2 = y_ii[0] - b1;
2060
2061 head_t1 = x_ii[0] * y_ii[0];
2062 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
2063 }
2064 {
2065 /* Compute double_double = double * double. */
2066 double a1, a2, b1, b2, con;
2067
2068 con = x_ii[1] * split;
2069 a1 = con - x_ii[1];
2070 a1 = con - a1;
2071 a2 = x_ii[1] - a1;
2072 con = y_ii[1] * split;
2073 b1 = con - y_ii[1];
2074 b1 = con - b1;
2075 b2 = y_ii[1] - b1;
2076
2077 head_t2 = x_ii[1] * y_ii[1];
2078 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
2079 }
2080 head_t2 = -head_t2;
2081 tail_t2 = -tail_t2;
2082 {
2083 /* Compute double-double = double-double + double-double. */
2084 double bv;
2085 double s1, s2, t1, t2;
2086
2087 /* Add two hi words. */
2088 s1 = head_t1 + head_t2;
2089 bv = s1 - head_t1;
2090 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2091
2092 /* Add two lo words. */
2093 t1 = tail_t1 + tail_t2;
2094 bv = t1 - tail_t1;
2095 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2096
2097 s2 += t1;
2098
2099 /* Renormalize (s1, s2) to (t1, s2) */
2100 t1 = s1 + s2;
2101 s2 = s2 - (t1 - s1);
2102
2103 t2 += s2;
2104
2105 /* Renormalize (t1, t2) */
2106 head_t1 = t1 + t2;
2107 tail_t1 = t2 - (head_t1 - t1);
2108 }
2109 head_prod[0] = head_t1;
2110 tail_prod[0] = tail_t1;
2111 /* Imaginary part */
2112 {
2113 /* Compute double_double = double * double. */
2114 double a1, a2, b1, b2, con;
2115
2116 con = x_ii[1] * split;
2117 a1 = con - x_ii[1];
2118 a1 = con - a1;
2119 a2 = x_ii[1] - a1;
2120 con = y_ii[0] * split;
2121 b1 = con - y_ii[0];
2122 b1 = con - b1;
2123 b2 = y_ii[0] - b1;
2124
2125 head_t1 = x_ii[1] * y_ii[0];
2126 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
2127 }
2128 {
2129 /* Compute double_double = double * double. */
2130 double a1, a2, b1, b2, con;
2131
2132 con = x_ii[0] * split;
2133 a1 = con - x_ii[0];
2134 a1 = con - a1;
2135 a2 = x_ii[0] - a1;
2136 con = y_ii[1] * split;
2137 b1 = con - y_ii[1];
2138 b1 = con - b1;
2139 b2 = y_ii[1] - b1;
2140
2141 head_t2 = x_ii[0] * y_ii[1];
2142 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
2143 }
2144 {
2145 /* Compute double-double = double-double + double-double. */
2146 double bv;
2147 double s1, s2, t1, t2;
2148
2149 /* Add two hi words. */
2150 s1 = head_t1 + head_t2;
2151 bv = s1 - head_t1;
2152 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2153
2154 /* Add two lo words. */
2155 t1 = tail_t1 + tail_t2;
2156 bv = t1 - tail_t1;
2157 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2158
2159 s2 += t1;
2160
2161 /* Renormalize (s1, s2) to (t1, s2) */
2162 t1 = s1 + s2;
2163 s2 = s2 - (t1 - s1);
2164
2165 t2 += s2;
2166
2167 /* Renormalize (t1, t2) */
2168 head_t1 = t1 + t2;
2169 tail_t1 = t2 - (head_t1 - t1);
2170 }
2171 head_prod[1] = head_t1;
2172 tail_prod[1] = tail_t1;
2173 } /* prod = x[i] * tail_y[i] */
2174 {
2175 double head_t, tail_t;
2176 double head_a, tail_a;
2177 double head_b, tail_b;
2178 /* Real part */
2179 head_a = head_sum[0];
2180 tail_a = tail_sum[0];
2181 head_b = head_prod[0];
2182 tail_b = tail_prod[0];
2183 {
2184 /* Compute double-double = double-double + double-double. */
2185 double bv;
2186 double s1, s2, t1, t2;
2187
2188 /* Add two hi words. */
2189 s1 = head_a + head_b;
2190 bv = s1 - head_a;
2191 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2192
2193 /* Add two lo words. */
2194 t1 = tail_a + tail_b;
2195 bv = t1 - tail_a;
2196 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2197
2198 s2 += t1;
2199
2200 /* Renormalize (s1, s2) to (t1, s2) */
2201 t1 = s1 + s2;
2202 s2 = s2 - (t1 - s1);
2203
2204 t2 += s2;
2205
2206 /* Renormalize (t1, t2) */
2207 head_t = t1 + t2;
2208 tail_t = t2 - (head_t - t1);
2209 }
2210 head_sum[0] = head_t;
2211 tail_sum[0] = tail_t;
2212 /* Imaginary part */
2213 head_a = head_sum[1];
2214 tail_a = tail_sum[1];
2215 head_b = head_prod[1];
2216 tail_b = tail_prod[1];
2217 {
2218 /* Compute double-double = double-double + double-double. */
2219 double bv;
2220 double s1, s2, t1, t2;
2221
2222 /* Add two hi words. */
2223 s1 = head_a + head_b;
2224 bv = s1 - head_a;
2225 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2226
2227 /* Add two lo words. */
2228 t1 = tail_a + tail_b;
2229 bv = t1 - tail_a;
2230 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2231
2232 s2 += t1;
2233
2234 /* Renormalize (s1, s2) to (t1, s2) */
2235 t1 = s1 + s2;
2236 s2 = s2 - (t1 - s1);
2237
2238 t2 += s2;
2239
2240 /* Renormalize (t1, t2) */
2241 head_t = t1 + t2;
2242 tail_t = t2 - (head_t - t1);
2243 }
2244 head_sum[1] = head_t;
2245 tail_sum[1] = tail_t;
2246 } /* sum = sum+prod */
2247 ix += incx;
2248 iy += incy;
2249 } /* endfor */
2250 }
2251
2252 {
2253 /* Compute complex-extra = complex-extra * complex-double. */
2254 double head_a0, tail_a0;
2255 double head_a1, tail_a1;
2256 double head_t1, tail_t1;
2257 double head_t2, tail_t2;
2258 head_a0 = head_sum[0];
2259 tail_a0 = tail_sum[0];
2260 head_a1 = head_sum[1];
2261 tail_a1 = tail_sum[1];
2262 /* real part */
2263 {
2264 /* Compute double-double = double-double * double. */
2265 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2266
2267 con = head_a0 * split;
2268 a11 = con - head_a0;
2269 a11 = con - a11;
2270 a21 = head_a0 - a11;
2271 con = alpha_i[0] * split;
2272 b1 = con - alpha_i[0];
2273 b1 = con - b1;
2274 b2 = alpha_i[0] - b1;
2275
2276 c11 = head_a0 * alpha_i[0];
2277 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2278
2279 c2 = tail_a0 * alpha_i[0];
2280 t1 = c11 + c2;
2281 t2 = (c2 - (t1 - c11)) + c21;
2282
2283 head_t1 = t1 + t2;
2284 tail_t1 = t2 - (head_t1 - t1);
2285 }
2286 {
2287 /* Compute double-double = double-double * double. */
2288 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2289
2290 con = head_a1 * split;
2291 a11 = con - head_a1;
2292 a11 = con - a11;
2293 a21 = head_a1 - a11;
2294 con = alpha_i[1] * split;
2295 b1 = con - alpha_i[1];
2296 b1 = con - b1;
2297 b2 = alpha_i[1] - b1;
2298
2299 c11 = head_a1 * alpha_i[1];
2300 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2301
2302 c2 = tail_a1 * alpha_i[1];
2303 t1 = c11 + c2;
2304 t2 = (c2 - (t1 - c11)) + c21;
2305
2306 head_t2 = t1 + t2;
2307 tail_t2 = t2 - (head_t2 - t1);
2308 }
2309 head_t2 = -head_t2;
2310 tail_t2 = -tail_t2;
2311 {
2312 /* Compute double-double = double-double + double-double. */
2313 double bv;
2314 double s1, s2, t1, t2;
2315
2316 /* Add two hi words. */
2317 s1 = head_t1 + head_t2;
2318 bv = s1 - head_t1;
2319 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2320
2321 /* Add two lo words. */
2322 t1 = tail_t1 + tail_t2;
2323 bv = t1 - tail_t1;
2324 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2325
2326 s2 += t1;
2327
2328 /* Renormalize (s1, s2) to (t1, s2) */
2329 t1 = s1 + s2;
2330 s2 = s2 - (t1 - s1);
2331
2332 t2 += s2;
2333
2334 /* Renormalize (t1, t2) */
2335 head_t1 = t1 + t2;
2336 tail_t1 = t2 - (head_t1 - t1);
2337 }
2338 head_tmp1[0] = head_t1;
2339 tail_tmp1[0] = tail_t1;
2340 /* imaginary part */
2341 {
2342 /* Compute double-double = double-double * double. */
2343 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2344
2345 con = head_a1 * split;
2346 a11 = con - head_a1;
2347 a11 = con - a11;
2348 a21 = head_a1 - a11;
2349 con = alpha_i[0] * split;
2350 b1 = con - alpha_i[0];
2351 b1 = con - b1;
2352 b2 = alpha_i[0] - b1;
2353
2354 c11 = head_a1 * alpha_i[0];
2355 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2356
2357 c2 = tail_a1 * alpha_i[0];
2358 t1 = c11 + c2;
2359 t2 = (c2 - (t1 - c11)) + c21;
2360
2361 head_t1 = t1 + t2;
2362 tail_t1 = t2 - (head_t1 - t1);
2363 }
2364 {
2365 /* Compute double-double = double-double * double. */
2366 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2367
2368 con = head_a0 * split;
2369 a11 = con - head_a0;
2370 a11 = con - a11;
2371 a21 = head_a0 - a11;
2372 con = alpha_i[1] * split;
2373 b1 = con - alpha_i[1];
2374 b1 = con - b1;
2375 b2 = alpha_i[1] - b1;
2376
2377 c11 = head_a0 * alpha_i[1];
2378 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2379
2380 c2 = tail_a0 * alpha_i[1];
2381 t1 = c11 + c2;
2382 t2 = (c2 - (t1 - c11)) + c21;
2383
2384 head_t2 = t1 + t2;
2385 tail_t2 = t2 - (head_t2 - t1);
2386 }
2387 {
2388 /* Compute double-double = double-double + double-double. */
2389 double bv;
2390 double s1, s2, t1, t2;
2391
2392 /* Add two hi words. */
2393 s1 = head_t1 + head_t2;
2394 bv = s1 - head_t1;
2395 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2396
2397 /* Add two lo words. */
2398 t1 = tail_t1 + tail_t2;
2399 bv = t1 - tail_t1;
2400 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2401
2402 s2 += t1;
2403
2404 /* Renormalize (s1, s2) to (t1, s2) */
2405 t1 = s1 + s2;
2406 s2 = s2 - (t1 - s1);
2407
2408 t2 += s2;
2409
2410 /* Renormalize (t1, t2) */
2411 head_t1 = t1 + t2;
2412 tail_t1 = t2 - (head_t1 - t1);
2413 }
2414 head_tmp1[1] = head_t1;
2415 tail_tmp1[1] = tail_t1;
2416 }
2417 /* tmp1 = sum*alpha */
2418 {
2419 /* Compute complex-extra = complex-double * complex-double. */
2420 double head_t1, tail_t1;
2421 double head_t2, tail_t2;
2422 /* Real part */
2423 {
2424 /* Compute double_double = double * double. */
2425 double a1, a2, b1, b2, con;
2426
2427 con = r_v[0] * split;
2428 a1 = con - r_v[0];
2429 a1 = con - a1;
2430 a2 = r_v[0] - a1;
2431 con = beta_i[0] * split;
2432 b1 = con - beta_i[0];
2433 b1 = con - b1;
2434 b2 = beta_i[0] - b1;
2435
2436 head_t1 = r_v[0] * beta_i[0];
2437 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
2438 }
2439 {
2440 /* Compute double_double = double * double. */
2441 double a1, a2, b1, b2, con;
2442
2443 con = r_v[1] * split;
2444 a1 = con - r_v[1];
2445 a1 = con - a1;
2446 a2 = r_v[1] - a1;
2447 con = beta_i[1] * split;
2448 b1 = con - beta_i[1];
2449 b1 = con - b1;
2450 b2 = beta_i[1] - b1;
2451
2452 head_t2 = r_v[1] * beta_i[1];
2453 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
2454 }
2455 head_t2 = -head_t2;
2456 tail_t2 = -tail_t2;
2457 {
2458 /* Compute double-double = double-double + double-double. */
2459 double bv;
2460 double s1, s2, t1, t2;
2461
2462 /* Add two hi words. */
2463 s1 = head_t1 + head_t2;
2464 bv = s1 - head_t1;
2465 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2466
2467 /* Add two lo words. */
2468 t1 = tail_t1 + tail_t2;
2469 bv = t1 - tail_t1;
2470 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2471
2472 s2 += t1;
2473
2474 /* Renormalize (s1, s2) to (t1, s2) */
2475 t1 = s1 + s2;
2476 s2 = s2 - (t1 - s1);
2477
2478 t2 += s2;
2479
2480 /* Renormalize (t1, t2) */
2481 head_t1 = t1 + t2;
2482 tail_t1 = t2 - (head_t1 - t1);
2483 }
2484 head_tmp2[0] = head_t1;
2485 tail_tmp2[0] = tail_t1;
2486 /* Imaginary part */
2487 {
2488 /* Compute double_double = double * double. */
2489 double a1, a2, b1, b2, con;
2490
2491 con = r_v[1] * split;
2492 a1 = con - r_v[1];
2493 a1 = con - a1;
2494 a2 = r_v[1] - a1;
2495 con = beta_i[0] * split;
2496 b1 = con - beta_i[0];
2497 b1 = con - b1;
2498 b2 = beta_i[0] - b1;
2499
2500 head_t1 = r_v[1] * beta_i[0];
2501 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
2502 }
2503 {
2504 /* Compute double_double = double * double. */
2505 double a1, a2, b1, b2, con;
2506
2507 con = r_v[0] * split;
2508 a1 = con - r_v[0];
2509 a1 = con - a1;
2510 a2 = r_v[0] - a1;
2511 con = beta_i[1] * split;
2512 b1 = con - beta_i[1];
2513 b1 = con - b1;
2514 b2 = beta_i[1] - b1;
2515
2516 head_t2 = r_v[0] * beta_i[1];
2517 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
2518 }
2519 {
2520 /* Compute double-double = double-double + double-double. */
2521 double bv;
2522 double s1, s2, t1, t2;
2523
2524 /* Add two hi words. */
2525 s1 = head_t1 + head_t2;
2526 bv = s1 - head_t1;
2527 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2528
2529 /* Add two lo words. */
2530 t1 = tail_t1 + tail_t2;
2531 bv = t1 - tail_t1;
2532 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2533
2534 s2 += t1;
2535
2536 /* Renormalize (s1, s2) to (t1, s2) */
2537 t1 = s1 + s2;
2538 s2 = s2 - (t1 - s1);
2539
2540 t2 += s2;
2541
2542 /* Renormalize (t1, t2) */
2543 head_t1 = t1 + t2;
2544 tail_t1 = t2 - (head_t1 - t1);
2545 }
2546 head_tmp2[1] = head_t1;
2547 tail_tmp2[1] = tail_t1;
2548 } /* tmp2 = r*beta */
2549 {
2550 double head_t, tail_t;
2551 double head_a, tail_a;
2552 double head_b, tail_b;
2553 /* Real part */
2554 head_a = head_tmp1[0];
2555 tail_a = tail_tmp1[0];
2556 head_b = head_tmp2[0];
2557 tail_b = tail_tmp2[0];
2558 {
2559 /* Compute double-double = double-double + double-double. */
2560 double bv;
2561 double s1, s2, t1, t2;
2562
2563 /* Add two hi words. */
2564 s1 = head_a + head_b;
2565 bv = s1 - head_a;
2566 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2567
2568 /* Add two lo words. */
2569 t1 = tail_a + tail_b;
2570 bv = t1 - tail_a;
2571 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2572
2573 s2 += t1;
2574
2575 /* Renormalize (s1, s2) to (t1, s2) */
2576 t1 = s1 + s2;
2577 s2 = s2 - (t1 - s1);
2578
2579 t2 += s2;
2580
2581 /* Renormalize (t1, t2) */
2582 head_t = t1 + t2;
2583 tail_t = t2 - (head_t - t1);
2584 }
2585 head_tmp1[0] = head_t;
2586 tail_tmp1[0] = tail_t;
2587 /* Imaginary part */
2588 head_a = head_tmp1[1];
2589 tail_a = tail_tmp1[1];
2590 head_b = head_tmp2[1];
2591 tail_b = tail_tmp2[1];
2592 {
2593 /* Compute double-double = double-double + double-double. */
2594 double bv;
2595 double s1, s2, t1, t2;
2596
2597 /* Add two hi words. */
2598 s1 = head_a + head_b;
2599 bv = s1 - head_a;
2600 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2601
2602 /* Add two lo words. */
2603 t1 = tail_a + tail_b;
2604 bv = t1 - tail_a;
2605 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2606
2607 s2 += t1;
2608
2609 /* Renormalize (s1, s2) to (t1, s2) */
2610 t1 = s1 + s2;
2611 s2 = s2 - (t1 - s1);
2612
2613 t2 += s2;
2614
2615 /* Renormalize (t1, t2) */
2616 head_t = t1 + t2;
2617 tail_t = t2 - (head_t - t1);
2618 }
2619 head_tmp1[1] = head_t;
2620 tail_tmp1[1] = tail_t;
2621 } /* tmp1 = tmp1+tmp2 */
2622 head_r_true[0] = head_tmp1[0];
2623 tail_r_true[0] = tail_tmp1[0];
2624 head_r_true[1] = head_tmp1[1];
2625 tail_r_true[1] = tail_tmp1[1]; /* *r = tmp1 */
2626
2627 FPU_FIX_STOP;
2628
2629 }
2630