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