1 #include "blas_extended.h"
2 #include "blas_extended_private.h"
BLAS_cgemv2_x(enum blas_order_type order,enum blas_trans_type trans,int m,int n,const void * alpha,const void * 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_cgemv2_x(enum blas_order_type order, enum blas_trans_type trans,
4 int m, int n, const void *alpha, const void *a, int lda,
5 const void *head_x, const void *tail_x, int incx,
6 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 void*
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_cgemv2_x";
62 switch (prec) {
63 case blas_prec_single:{
64
65 int i, j;
66 int iy, jx, kx, ky;
67 int lenx, leny;
68 int ai, aij;
69 int incai, incaij;
70
71 const float *a_i = (float *) a;
72 const float *head_x_i = (float *) head_x;
73 const float *tail_x_i = (float *) tail_x;
74 float *y_i = (float *) y;
75 float *alpha_i = (float *) alpha;
76 float *beta_i = (float *) beta;
77 float a_elem[2];
78 float x_elem[2];
79 float y_elem[2];
80 float prod[2];
81 float sum[2];
82 float sum2[2];
83 float tmp1[2];
84 float tmp2[2];
85
86
87 /* all error calls */
88 if (m < 0)
89 BLAS_error(routine_name, -3, m, 0);
90 else if (n <= 0)
91 BLAS_error(routine_name, -4, n, 0);
92 else if (incx == 0)
93 BLAS_error(routine_name, -10, incx, 0);
94 else if (incy == 0)
95 BLAS_error(routine_name, -13, incy, 0);
96
97 if ((order == blas_rowmajor) && (trans == blas_no_trans)) {
98 lenx = n;
99 leny = m;
100 incai = lda;
101 incaij = 1;
102 } else if ((order == blas_rowmajor) && (trans != blas_no_trans)) {
103 lenx = m;
104 leny = n;
105 incai = 1;
106 incaij = lda;
107 } else if ((order == blas_colmajor) && (trans == blas_no_trans)) {
108 lenx = n;
109 leny = m;
110 incai = 1;
111 incaij = lda;
112 } else { /* colmajor and blas_trans */
113 lenx = m;
114 leny = n;
115 incai = lda;
116 incaij = 1;
117 }
118
119 if (lda < leny)
120 BLAS_error(routine_name, -7, lda, NULL);
121
122
123
124 incx *= 2;
125 incy *= 2;
126 incai *= 2;
127 incaij *= 2;
128
129 if (incx > 0)
130 kx = 0;
131 else
132 kx = (1 - lenx) * incx;
133 if (incy > 0)
134 ky = 0;
135 else
136 ky = (1 - leny) * incy;
137
138 /* No extra-precision needed for alpha = 0 */
139 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0) {
140 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
141 iy = ky;
142 for (i = 0; i < leny; i++) {
143 y_i[iy] = 0.0;
144 y_i[iy + 1] = 0.0;
145 iy += incy;
146 }
147 } else if (!(beta_i[0] == 0.0 && beta_i[1] == 0.0)) {
148 iy = ky;
149 for (i = 0; i < leny; i++) {
150 y_elem[0] = y_i[iy];
151 y_elem[1] = y_i[iy + 1];
152 {
153 tmp1[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
154 tmp1[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
155 }
156
157 y_i[iy] = tmp1[0];
158 y_i[iy + 1] = tmp1[1];
159 iy += incy;
160 }
161 }
162 } else { /* alpha != 0 */
163 if (trans == blas_conj_trans) {
164
165 /* if beta = 0, we can save m multiplies:
166 y = alpha*A*head_x + alpha*A*tail_x */
167 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
168 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
169 /* save m more multiplies if alpha = 1 */
170 ai = 0;
171 iy = ky;
172 for (i = 0; i < leny; i++) {
173 sum[0] = sum[1] = 0.0;
174 sum2[0] = sum2[1] = 0.0;
175 aij = ai;
176 jx = kx;
177 for (j = 0; j < lenx; j++) {
178 a_elem[0] = a_i[aij];
179 a_elem[1] = a_i[aij + 1];
180 a_elem[1] = -a_elem[1];
181 x_elem[0] = head_x_i[jx];
182 x_elem[1] = head_x_i[jx + 1];
183 {
184 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
185 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
186 }
187
188 sum[0] = sum[0] + prod[0];
189 sum[1] = sum[1] + prod[1];
190 x_elem[0] = tail_x_i[jx];
191 x_elem[1] = tail_x_i[jx + 1];
192 {
193 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
194 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
195 }
196
197 sum2[0] = sum2[0] + prod[0];
198 sum2[1] = sum2[1] + prod[1];
199 aij += incaij;
200 jx += incx;
201 }
202 sum[0] = sum[0] + sum2[0];
203 sum[1] = sum[1] + sum2[1];
204 y_i[iy] = sum[0];
205 y_i[iy + 1] = sum[1];
206 ai += incai;
207 iy += incy;
208 } /* end for */
209 } else { /* alpha != 1 */
210 ai = 0;
211 iy = ky;
212 for (i = 0; i < leny; i++) {
213 sum[0] = sum[1] = 0.0;
214 sum2[0] = sum2[1] = 0.0;
215 aij = ai;
216 jx = kx;
217 for (j = 0; j < lenx; j++) {
218 a_elem[0] = a_i[aij];
219 a_elem[1] = a_i[aij + 1];
220 a_elem[1] = -a_elem[1];
221 x_elem[0] = head_x_i[jx];
222 x_elem[1] = head_x_i[jx + 1];
223 {
224 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
225 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
226 }
227
228 sum[0] = sum[0] + prod[0];
229 sum[1] = sum[1] + prod[1];
230 x_elem[0] = tail_x_i[jx];
231 x_elem[1] = tail_x_i[jx + 1];
232 {
233 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
234 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
235 }
236
237 sum2[0] = sum2[0] + prod[0];
238 sum2[1] = sum2[1] + prod[1];
239 aij += incaij;
240 jx += incx;
241 }
242 {
243 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
244 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
245 }
246
247 {
248 tmp2[0] = sum2[0] * alpha_i[0] - sum2[1] * alpha_i[1];
249 tmp2[1] = sum2[0] * alpha_i[1] + sum2[1] * alpha_i[0];
250 }
251
252 tmp1[0] = tmp1[0] + tmp2[0];
253 tmp1[1] = tmp1[1] + tmp2[1];
254 y_i[iy] = tmp1[0];
255 y_i[iy + 1] = tmp1[1];
256 ai += incai;
257 iy += incy;
258 }
259 }
260 } else { /* beta != 0 */
261 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
262 /* save m multiplies if alpha = 1 */
263 ai = 0;
264 iy = ky;
265 for (i = 0; i < leny; i++) {
266 sum[0] = sum[1] = 0.0;;
267 sum2[0] = sum2[1] = 0.0;;
268 aij = ai;
269 jx = kx;
270 for (j = 0; j < lenx; j++) {
271 a_elem[0] = a_i[aij];
272 a_elem[1] = a_i[aij + 1];
273 a_elem[1] = -a_elem[1];
274 x_elem[0] = head_x_i[jx];
275 x_elem[1] = head_x_i[jx + 1];
276 {
277 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
278 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
279 }
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] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
287 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
288 }
289
290 sum2[0] = sum2[0] + prod[0];
291 sum2[1] = sum2[1] + prod[1];
292 aij += incaij;
293 jx += incx;
294 }
295 sum[0] = sum[0] + sum2[0];
296 sum[1] = sum[1] + sum2[1];
297 y_elem[0] = y_i[iy];
298 y_elem[1] = y_i[iy + 1];
299 {
300 tmp1[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
301 tmp1[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
302 }
303
304 tmp2[0] = sum[0] + tmp1[0];
305 tmp2[1] = sum[1] + tmp1[1];
306 y_i[iy] = tmp2[0];
307 y_i[iy + 1] = tmp2[1];
308 ai += incai;
309 iy += incy;
310 }
311 } else { /* alpha != 1, the most general form:
312 y = alpha*A*head_x + alpha*A*tail_x + beta*y */
313 ai = 0;
314 iy = ky;
315 for (i = 0; i < leny; i++) {
316 sum[0] = sum[1] = 0.0;;
317 sum2[0] = sum2[1] = 0.0;;
318 aij = ai;
319 jx = kx;
320 for (j = 0; j < lenx; j++) {
321 a_elem[0] = a_i[aij];
322 a_elem[1] = a_i[aij + 1];
323 a_elem[1] = -a_elem[1];
324 x_elem[0] = head_x_i[jx];
325 x_elem[1] = head_x_i[jx + 1];
326 {
327 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
328 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
329 }
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] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
337 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
338 }
339
340 sum2[0] = sum2[0] + prod[0];
341 sum2[1] = sum2[1] + prod[1];
342 aij += incaij;
343 jx += incx;
344 }
345 {
346 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
347 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
348 }
349
350 {
351 tmp2[0] = sum2[0] * alpha_i[0] - sum2[1] * alpha_i[1];
352 tmp2[1] = sum2[0] * alpha_i[1] + sum2[1] * alpha_i[0];
353 }
354
355 tmp1[0] = tmp1[0] + tmp2[0];
356 tmp1[1] = tmp1[1] + tmp2[1];
357 y_elem[0] = y_i[iy];
358 y_elem[1] = y_i[iy + 1];
359 {
360 tmp2[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
361 tmp2[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
362 }
363
364 tmp1[0] = tmp1[0] + tmp2[0];
365 tmp1[1] = tmp1[1] + tmp2[1];
366 y_i[iy] = tmp1[0];
367 y_i[iy + 1] = tmp1[1];
368 ai += incai;
369 iy += incy;
370 }
371 }
372 }
373
374 } else {
375
376 /* if beta = 0, we can save m multiplies:
377 y = alpha*A*head_x + alpha*A*tail_x */
378 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
379 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
380 /* save m more multiplies if alpha = 1 */
381 ai = 0;
382 iy = ky;
383 for (i = 0; i < leny; i++) {
384 sum[0] = sum[1] = 0.0;
385 sum2[0] = sum2[1] = 0.0;
386 aij = ai;
387 jx = kx;
388 for (j = 0; j < lenx; j++) {
389 a_elem[0] = a_i[aij];
390 a_elem[1] = a_i[aij + 1];
391
392 x_elem[0] = head_x_i[jx];
393 x_elem[1] = head_x_i[jx + 1];
394 {
395 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
396 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
397 }
398
399 sum[0] = sum[0] + prod[0];
400 sum[1] = sum[1] + prod[1];
401 x_elem[0] = tail_x_i[jx];
402 x_elem[1] = tail_x_i[jx + 1];
403 {
404 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
405 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
406 }
407
408 sum2[0] = sum2[0] + prod[0];
409 sum2[1] = sum2[1] + prod[1];
410 aij += incaij;
411 jx += incx;
412 }
413 sum[0] = sum[0] + sum2[0];
414 sum[1] = sum[1] + sum2[1];
415 y_i[iy] = sum[0];
416 y_i[iy + 1] = sum[1];
417 ai += incai;
418 iy += incy;
419 } /* end for */
420 } else { /* alpha != 1 */
421 ai = 0;
422 iy = ky;
423 for (i = 0; i < leny; i++) {
424 sum[0] = sum[1] = 0.0;
425 sum2[0] = sum2[1] = 0.0;
426 aij = ai;
427 jx = kx;
428 for (j = 0; j < lenx; j++) {
429 a_elem[0] = a_i[aij];
430 a_elem[1] = a_i[aij + 1];
431
432 x_elem[0] = head_x_i[jx];
433 x_elem[1] = head_x_i[jx + 1];
434 {
435 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
436 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
437 }
438
439 sum[0] = sum[0] + prod[0];
440 sum[1] = sum[1] + prod[1];
441 x_elem[0] = tail_x_i[jx];
442 x_elem[1] = tail_x_i[jx + 1];
443 {
444 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
445 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
446 }
447
448 sum2[0] = sum2[0] + prod[0];
449 sum2[1] = sum2[1] + prod[1];
450 aij += incaij;
451 jx += incx;
452 }
453 {
454 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
455 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
456 }
457
458 {
459 tmp2[0] = sum2[0] * alpha_i[0] - sum2[1] * alpha_i[1];
460 tmp2[1] = sum2[0] * alpha_i[1] + sum2[1] * alpha_i[0];
461 }
462
463 tmp1[0] = tmp1[0] + tmp2[0];
464 tmp1[1] = tmp1[1] + tmp2[1];
465 y_i[iy] = tmp1[0];
466 y_i[iy + 1] = tmp1[1];
467 ai += incai;
468 iy += incy;
469 }
470 }
471 } else { /* beta != 0 */
472 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
473 /* save m multiplies if alpha = 1 */
474 ai = 0;
475 iy = ky;
476 for (i = 0; i < leny; i++) {
477 sum[0] = sum[1] = 0.0;;
478 sum2[0] = sum2[1] = 0.0;;
479 aij = ai;
480 jx = kx;
481 for (j = 0; j < lenx; j++) {
482 a_elem[0] = a_i[aij];
483 a_elem[1] = a_i[aij + 1];
484
485 x_elem[0] = head_x_i[jx];
486 x_elem[1] = head_x_i[jx + 1];
487 {
488 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
489 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
490 }
491
492 sum[0] = sum[0] + prod[0];
493 sum[1] = sum[1] + prod[1];
494 x_elem[0] = tail_x_i[jx];
495 x_elem[1] = tail_x_i[jx + 1];
496 {
497 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
498 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
499 }
500
501 sum2[0] = sum2[0] + prod[0];
502 sum2[1] = sum2[1] + prod[1];
503 aij += incaij;
504 jx += incx;
505 }
506 sum[0] = sum[0] + sum2[0];
507 sum[1] = sum[1] + sum2[1];
508 y_elem[0] = y_i[iy];
509 y_elem[1] = y_i[iy + 1];
510 {
511 tmp1[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
512 tmp1[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
513 }
514
515 tmp2[0] = sum[0] + tmp1[0];
516 tmp2[1] = sum[1] + tmp1[1];
517 y_i[iy] = tmp2[0];
518 y_i[iy + 1] = tmp2[1];
519 ai += incai;
520 iy += incy;
521 }
522 } else { /* alpha != 1, the most general form:
523 y = alpha*A*head_x + alpha*A*tail_x + beta*y */
524 ai = 0;
525 iy = ky;
526 for (i = 0; i < leny; i++) {
527 sum[0] = sum[1] = 0.0;;
528 sum2[0] = sum2[1] = 0.0;;
529 aij = ai;
530 jx = kx;
531 for (j = 0; j < lenx; j++) {
532 a_elem[0] = a_i[aij];
533 a_elem[1] = a_i[aij + 1];
534
535 x_elem[0] = head_x_i[jx];
536 x_elem[1] = head_x_i[jx + 1];
537 {
538 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
539 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
540 }
541
542 sum[0] = sum[0] + prod[0];
543 sum[1] = sum[1] + prod[1];
544 x_elem[0] = tail_x_i[jx];
545 x_elem[1] = tail_x_i[jx + 1];
546 {
547 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
548 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
549 }
550
551 sum2[0] = sum2[0] + prod[0];
552 sum2[1] = sum2[1] + prod[1];
553 aij += incaij;
554 jx += incx;
555 }
556 {
557 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
558 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
559 }
560
561 {
562 tmp2[0] = sum2[0] * alpha_i[0] - sum2[1] * alpha_i[1];
563 tmp2[1] = sum2[0] * alpha_i[1] + sum2[1] * alpha_i[0];
564 }
565
566 tmp1[0] = tmp1[0] + tmp2[0];
567 tmp1[1] = tmp1[1] + tmp2[1];
568 y_elem[0] = y_i[iy];
569 y_elem[1] = y_i[iy + 1];
570 {
571 tmp2[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
572 tmp2[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
573 }
574
575 tmp1[0] = tmp1[0] + tmp2[0];
576 tmp1[1] = tmp1[1] + tmp2[1];
577 y_i[iy] = tmp1[0];
578 y_i[iy + 1] = tmp1[1];
579 ai += incai;
580 iy += incy;
581 }
582 }
583 }
584
585 }
586 }
587
588
589
590 break;
591 }
592 case blas_prec_double:
593 case blas_prec_indigenous:{
594
595 int i, j;
596 int iy, jx, kx, ky;
597 int lenx, leny;
598 int ai, aij;
599 int incai, incaij;
600
601 const float *a_i = (float *) a;
602 const float *head_x_i = (float *) head_x;
603 const float *tail_x_i = (float *) tail_x;
604 float *y_i = (float *) y;
605 float *alpha_i = (float *) alpha;
606 float *beta_i = (float *) beta;
607 float a_elem[2];
608 float x_elem[2];
609 float y_elem[2];
610 double prod[2];
611 double sum[2];
612 double sum2[2];
613 double tmp1[2];
614 double tmp2[2];
615
616
617 /* all error calls */
618 if (m < 0)
619 BLAS_error(routine_name, -3, m, 0);
620 else if (n <= 0)
621 BLAS_error(routine_name, -4, n, 0);
622 else if (incx == 0)
623 BLAS_error(routine_name, -10, incx, 0);
624 else if (incy == 0)
625 BLAS_error(routine_name, -13, incy, 0);
626
627 if ((order == blas_rowmajor) && (trans == blas_no_trans)) {
628 lenx = n;
629 leny = m;
630 incai = lda;
631 incaij = 1;
632 } else if ((order == blas_rowmajor) && (trans != blas_no_trans)) {
633 lenx = m;
634 leny = n;
635 incai = 1;
636 incaij = lda;
637 } else if ((order == blas_colmajor) && (trans == blas_no_trans)) {
638 lenx = n;
639 leny = m;
640 incai = 1;
641 incaij = lda;
642 } else { /* colmajor and blas_trans */
643 lenx = m;
644 leny = n;
645 incai = lda;
646 incaij = 1;
647 }
648
649 if (lda < leny)
650 BLAS_error(routine_name, -7, lda, NULL);
651
652
653
654 incx *= 2;
655 incy *= 2;
656 incai *= 2;
657 incaij *= 2;
658
659 if (incx > 0)
660 kx = 0;
661 else
662 kx = (1 - lenx) * incx;
663 if (incy > 0)
664 ky = 0;
665 else
666 ky = (1 - leny) * incy;
667
668 /* No extra-precision needed for alpha = 0 */
669 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0) {
670 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
671 iy = ky;
672 for (i = 0; i < leny; i++) {
673 y_i[iy] = 0.0;
674 y_i[iy + 1] = 0.0;
675 iy += incy;
676 }
677 } else if (!(beta_i[0] == 0.0 && beta_i[1] == 0.0)) {
678 iy = ky;
679 for (i = 0; i < leny; i++) {
680 y_elem[0] = y_i[iy];
681 y_elem[1] = y_i[iy + 1];
682 {
683 tmp1[0] =
684 (double) y_elem[0] * beta_i[0] -
685 (double) y_elem[1] * beta_i[1];
686 tmp1[1] =
687 (double) y_elem[0] * beta_i[1] +
688 (double) y_elem[1] * beta_i[0];
689 }
690 y_i[iy] = tmp1[0];
691 y_i[iy + 1] = tmp1[1];
692 iy += incy;
693 }
694 }
695 } else { /* alpha != 0 */
696 if (trans == blas_conj_trans) {
697
698 /* if beta = 0, we can save m multiplies:
699 y = alpha*A*head_x + alpha*A*tail_x */
700 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
701 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
702 /* save m more multiplies if alpha = 1 */
703 ai = 0;
704 iy = ky;
705 for (i = 0; i < leny; i++) {
706 sum[0] = sum[1] = 0.0;
707 sum2[0] = sum2[1] = 0.0;
708 aij = ai;
709 jx = kx;
710 for (j = 0; j < lenx; j++) {
711 a_elem[0] = a_i[aij];
712 a_elem[1] = a_i[aij + 1];
713 a_elem[1] = -a_elem[1];
714 x_elem[0] = head_x_i[jx];
715 x_elem[1] = head_x_i[jx + 1];
716 {
717 prod[0] =
718 (double) a_elem[0] * x_elem[0] -
719 (double) a_elem[1] * x_elem[1];
720 prod[1] =
721 (double) a_elem[0] * x_elem[1] +
722 (double) a_elem[1] * x_elem[0];
723 }
724 sum[0] = sum[0] + prod[0];
725 sum[1] = sum[1] + prod[1];
726 x_elem[0] = tail_x_i[jx];
727 x_elem[1] = tail_x_i[jx + 1];
728 {
729 prod[0] =
730 (double) a_elem[0] * x_elem[0] -
731 (double) a_elem[1] * x_elem[1];
732 prod[1] =
733 (double) a_elem[0] * x_elem[1] +
734 (double) a_elem[1] * x_elem[0];
735 }
736 sum2[0] = sum2[0] + prod[0];
737 sum2[1] = sum2[1] + prod[1];
738 aij += incaij;
739 jx += incx;
740 }
741 sum[0] = sum[0] + sum2[0];
742 sum[1] = sum[1] + sum2[1];
743 y_i[iy] = sum[0];
744 y_i[iy + 1] = sum[1];
745 ai += incai;
746 iy += incy;
747 } /* end for */
748 } else { /* alpha != 1 */
749 ai = 0;
750 iy = ky;
751 for (i = 0; i < leny; i++) {
752 sum[0] = sum[1] = 0.0;
753 sum2[0] = sum2[1] = 0.0;
754 aij = ai;
755 jx = kx;
756 for (j = 0; j < lenx; j++) {
757 a_elem[0] = a_i[aij];
758 a_elem[1] = a_i[aij + 1];
759 a_elem[1] = -a_elem[1];
760 x_elem[0] = head_x_i[jx];
761 x_elem[1] = head_x_i[jx + 1];
762 {
763 prod[0] =
764 (double) a_elem[0] * x_elem[0] -
765 (double) a_elem[1] * x_elem[1];
766 prod[1] =
767 (double) a_elem[0] * x_elem[1] +
768 (double) a_elem[1] * x_elem[0];
769 }
770 sum[0] = sum[0] + prod[0];
771 sum[1] = sum[1] + prod[1];
772 x_elem[0] = tail_x_i[jx];
773 x_elem[1] = tail_x_i[jx + 1];
774 {
775 prod[0] =
776 (double) a_elem[0] * x_elem[0] -
777 (double) a_elem[1] * x_elem[1];
778 prod[1] =
779 (double) a_elem[0] * x_elem[1] +
780 (double) a_elem[1] * x_elem[0];
781 }
782 sum2[0] = sum2[0] + prod[0];
783 sum2[1] = sum2[1] + prod[1];
784 aij += incaij;
785 jx += incx;
786 }
787 {
788 tmp1[0] =
789 (double) sum[0] * alpha_i[0] -
790 (double) sum[1] * alpha_i[1];
791 tmp1[1] =
792 (double) sum[0] * alpha_i[1] +
793 (double) sum[1] * alpha_i[0];
794 }
795 {
796 tmp2[0] =
797 (double) sum2[0] * alpha_i[0] -
798 (double) sum2[1] * alpha_i[1];
799 tmp2[1] =
800 (double) sum2[0] * alpha_i[1] +
801 (double) sum2[1] * alpha_i[0];
802 }
803 tmp1[0] = tmp1[0] + tmp2[0];
804 tmp1[1] = tmp1[1] + tmp2[1];
805 y_i[iy] = tmp1[0];
806 y_i[iy + 1] = tmp1[1];
807 ai += incai;
808 iy += incy;
809 }
810 }
811 } else { /* beta != 0 */
812 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
813 /* save m multiplies if alpha = 1 */
814 ai = 0;
815 iy = ky;
816 for (i = 0; i < leny; i++) {
817 sum[0] = sum[1] = 0.0;;
818 sum2[0] = sum2[1] = 0.0;;
819 aij = ai;
820 jx = kx;
821 for (j = 0; j < lenx; j++) {
822 a_elem[0] = a_i[aij];
823 a_elem[1] = a_i[aij + 1];
824 a_elem[1] = -a_elem[1];
825 x_elem[0] = head_x_i[jx];
826 x_elem[1] = head_x_i[jx + 1];
827 {
828 prod[0] =
829 (double) a_elem[0] * x_elem[0] -
830 (double) a_elem[1] * x_elem[1];
831 prod[1] =
832 (double) a_elem[0] * x_elem[1] +
833 (double) a_elem[1] * x_elem[0];
834 }
835 sum[0] = sum[0] + prod[0];
836 sum[1] = sum[1] + prod[1];
837 x_elem[0] = tail_x_i[jx];
838 x_elem[1] = tail_x_i[jx + 1];
839 {
840 prod[0] =
841 (double) a_elem[0] * x_elem[0] -
842 (double) a_elem[1] * x_elem[1];
843 prod[1] =
844 (double) a_elem[0] * x_elem[1] +
845 (double) a_elem[1] * x_elem[0];
846 }
847 sum2[0] = sum2[0] + prod[0];
848 sum2[1] = sum2[1] + prod[1];
849 aij += incaij;
850 jx += incx;
851 }
852 sum[0] = sum[0] + sum2[0];
853 sum[1] = sum[1] + sum2[1];
854 y_elem[0] = y_i[iy];
855 y_elem[1] = y_i[iy + 1];
856 {
857 tmp1[0] =
858 (double) y_elem[0] * beta_i[0] -
859 (double) y_elem[1] * beta_i[1];
860 tmp1[1] =
861 (double) y_elem[0] * beta_i[1] +
862 (double) y_elem[1] * beta_i[0];
863 }
864 tmp2[0] = sum[0] + tmp1[0];
865 tmp2[1] = sum[1] + tmp1[1];
866 y_i[iy] = tmp2[0];
867 y_i[iy + 1] = tmp2[1];
868 ai += incai;
869 iy += incy;
870 }
871 } else { /* alpha != 1, the most general form:
872 y = alpha*A*head_x + alpha*A*tail_x + beta*y */
873 ai = 0;
874 iy = ky;
875 for (i = 0; i < leny; i++) {
876 sum[0] = sum[1] = 0.0;;
877 sum2[0] = sum2[1] = 0.0;;
878 aij = ai;
879 jx = kx;
880 for (j = 0; j < lenx; j++) {
881 a_elem[0] = a_i[aij];
882 a_elem[1] = a_i[aij + 1];
883 a_elem[1] = -a_elem[1];
884 x_elem[0] = head_x_i[jx];
885 x_elem[1] = head_x_i[jx + 1];
886 {
887 prod[0] =
888 (double) a_elem[0] * x_elem[0] -
889 (double) a_elem[1] * x_elem[1];
890 prod[1] =
891 (double) a_elem[0] * x_elem[1] +
892 (double) a_elem[1] * x_elem[0];
893 }
894 sum[0] = sum[0] + prod[0];
895 sum[1] = sum[1] + prod[1];
896 x_elem[0] = tail_x_i[jx];
897 x_elem[1] = tail_x_i[jx + 1];
898 {
899 prod[0] =
900 (double) a_elem[0] * x_elem[0] -
901 (double) a_elem[1] * x_elem[1];
902 prod[1] =
903 (double) a_elem[0] * x_elem[1] +
904 (double) a_elem[1] * x_elem[0];
905 }
906 sum2[0] = sum2[0] + prod[0];
907 sum2[1] = sum2[1] + prod[1];
908 aij += incaij;
909 jx += incx;
910 }
911 {
912 tmp1[0] =
913 (double) sum[0] * alpha_i[0] -
914 (double) sum[1] * alpha_i[1];
915 tmp1[1] =
916 (double) sum[0] * alpha_i[1] +
917 (double) sum[1] * alpha_i[0];
918 }
919 {
920 tmp2[0] =
921 (double) sum2[0] * alpha_i[0] -
922 (double) sum2[1] * alpha_i[1];
923 tmp2[1] =
924 (double) sum2[0] * alpha_i[1] +
925 (double) sum2[1] * alpha_i[0];
926 }
927 tmp1[0] = tmp1[0] + tmp2[0];
928 tmp1[1] = tmp1[1] + tmp2[1];
929 y_elem[0] = y_i[iy];
930 y_elem[1] = y_i[iy + 1];
931 {
932 tmp2[0] =
933 (double) y_elem[0] * beta_i[0] -
934 (double) y_elem[1] * beta_i[1];
935 tmp2[1] =
936 (double) y_elem[0] * beta_i[1] +
937 (double) y_elem[1] * beta_i[0];
938 }
939 tmp1[0] = tmp1[0] + tmp2[0];
940 tmp1[1] = tmp1[1] + tmp2[1];
941 y_i[iy] = tmp1[0];
942 y_i[iy + 1] = tmp1[1];
943 ai += incai;
944 iy += incy;
945 }
946 }
947 }
948
949 } else {
950
951 /* if beta = 0, we can save m multiplies:
952 y = alpha*A*head_x + alpha*A*tail_x */
953 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
954 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
955 /* save m more multiplies if alpha = 1 */
956 ai = 0;
957 iy = ky;
958 for (i = 0; i < leny; i++) {
959 sum[0] = sum[1] = 0.0;
960 sum2[0] = sum2[1] = 0.0;
961 aij = ai;
962 jx = kx;
963 for (j = 0; j < lenx; j++) {
964 a_elem[0] = a_i[aij];
965 a_elem[1] = a_i[aij + 1];
966
967 x_elem[0] = head_x_i[jx];
968 x_elem[1] = head_x_i[jx + 1];
969 {
970 prod[0] =
971 (double) a_elem[0] * x_elem[0] -
972 (double) a_elem[1] * x_elem[1];
973 prod[1] =
974 (double) a_elem[0] * x_elem[1] +
975 (double) a_elem[1] * x_elem[0];
976 }
977 sum[0] = sum[0] + prod[0];
978 sum[1] = sum[1] + prod[1];
979 x_elem[0] = tail_x_i[jx];
980 x_elem[1] = tail_x_i[jx + 1];
981 {
982 prod[0] =
983 (double) a_elem[0] * x_elem[0] -
984 (double) a_elem[1] * x_elem[1];
985 prod[1] =
986 (double) a_elem[0] * x_elem[1] +
987 (double) a_elem[1] * x_elem[0];
988 }
989 sum2[0] = sum2[0] + prod[0];
990 sum2[1] = sum2[1] + prod[1];
991 aij += incaij;
992 jx += incx;
993 }
994 sum[0] = sum[0] + sum2[0];
995 sum[1] = sum[1] + sum2[1];
996 y_i[iy] = sum[0];
997 y_i[iy + 1] = sum[1];
998 ai += incai;
999 iy += incy;
1000 } /* end for */
1001 } else { /* alpha != 1 */
1002 ai = 0;
1003 iy = ky;
1004 for (i = 0; i < leny; i++) {
1005 sum[0] = sum[1] = 0.0;
1006 sum2[0] = sum2[1] = 0.0;
1007 aij = ai;
1008 jx = kx;
1009 for (j = 0; j < lenx; j++) {
1010 a_elem[0] = a_i[aij];
1011 a_elem[1] = a_i[aij + 1];
1012
1013 x_elem[0] = head_x_i[jx];
1014 x_elem[1] = head_x_i[jx + 1];
1015 {
1016 prod[0] =
1017 (double) a_elem[0] * x_elem[0] -
1018 (double) a_elem[1] * x_elem[1];
1019 prod[1] =
1020 (double) a_elem[0] * x_elem[1] +
1021 (double) a_elem[1] * x_elem[0];
1022 }
1023 sum[0] = sum[0] + prod[0];
1024 sum[1] = sum[1] + prod[1];
1025 x_elem[0] = tail_x_i[jx];
1026 x_elem[1] = tail_x_i[jx + 1];
1027 {
1028 prod[0] =
1029 (double) a_elem[0] * x_elem[0] -
1030 (double) a_elem[1] * x_elem[1];
1031 prod[1] =
1032 (double) a_elem[0] * x_elem[1] +
1033 (double) a_elem[1] * x_elem[0];
1034 }
1035 sum2[0] = sum2[0] + prod[0];
1036 sum2[1] = sum2[1] + prod[1];
1037 aij += incaij;
1038 jx += incx;
1039 }
1040 {
1041 tmp1[0] =
1042 (double) sum[0] * alpha_i[0] -
1043 (double) sum[1] * alpha_i[1];
1044 tmp1[1] =
1045 (double) sum[0] * alpha_i[1] +
1046 (double) sum[1] * alpha_i[0];
1047 }
1048 {
1049 tmp2[0] =
1050 (double) sum2[0] * alpha_i[0] -
1051 (double) sum2[1] * alpha_i[1];
1052 tmp2[1] =
1053 (double) sum2[0] * alpha_i[1] +
1054 (double) sum2[1] * alpha_i[0];
1055 }
1056 tmp1[0] = tmp1[0] + tmp2[0];
1057 tmp1[1] = tmp1[1] + tmp2[1];
1058 y_i[iy] = tmp1[0];
1059 y_i[iy + 1] = tmp1[1];
1060 ai += incai;
1061 iy += incy;
1062 }
1063 }
1064 } else { /* beta != 0 */
1065 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
1066 /* save m multiplies if alpha = 1 */
1067 ai = 0;
1068 iy = ky;
1069 for (i = 0; i < leny; i++) {
1070 sum[0] = sum[1] = 0.0;;
1071 sum2[0] = sum2[1] = 0.0;;
1072 aij = ai;
1073 jx = kx;
1074 for (j = 0; j < lenx; j++) {
1075 a_elem[0] = a_i[aij];
1076 a_elem[1] = a_i[aij + 1];
1077
1078 x_elem[0] = head_x_i[jx];
1079 x_elem[1] = head_x_i[jx + 1];
1080 {
1081 prod[0] =
1082 (double) a_elem[0] * x_elem[0] -
1083 (double) a_elem[1] * x_elem[1];
1084 prod[1] =
1085 (double) a_elem[0] * x_elem[1] +
1086 (double) a_elem[1] * x_elem[0];
1087 }
1088 sum[0] = sum[0] + prod[0];
1089 sum[1] = sum[1] + prod[1];
1090 x_elem[0] = tail_x_i[jx];
1091 x_elem[1] = tail_x_i[jx + 1];
1092 {
1093 prod[0] =
1094 (double) a_elem[0] * x_elem[0] -
1095 (double) a_elem[1] * x_elem[1];
1096 prod[1] =
1097 (double) a_elem[0] * x_elem[1] +
1098 (double) a_elem[1] * x_elem[0];
1099 }
1100 sum2[0] = sum2[0] + prod[0];
1101 sum2[1] = sum2[1] + prod[1];
1102 aij += incaij;
1103 jx += incx;
1104 }
1105 sum[0] = sum[0] + sum2[0];
1106 sum[1] = sum[1] + sum2[1];
1107 y_elem[0] = y_i[iy];
1108 y_elem[1] = y_i[iy + 1];
1109 {
1110 tmp1[0] =
1111 (double) y_elem[0] * beta_i[0] -
1112 (double) y_elem[1] * beta_i[1];
1113 tmp1[1] =
1114 (double) y_elem[0] * beta_i[1] +
1115 (double) y_elem[1] * beta_i[0];
1116 }
1117 tmp2[0] = sum[0] + tmp1[0];
1118 tmp2[1] = sum[1] + tmp1[1];
1119 y_i[iy] = tmp2[0];
1120 y_i[iy + 1] = tmp2[1];
1121 ai += incai;
1122 iy += incy;
1123 }
1124 } else { /* alpha != 1, the most general form:
1125 y = alpha*A*head_x + alpha*A*tail_x + beta*y */
1126 ai = 0;
1127 iy = ky;
1128 for (i = 0; i < leny; i++) {
1129 sum[0] = sum[1] = 0.0;;
1130 sum2[0] = sum2[1] = 0.0;;
1131 aij = ai;
1132 jx = kx;
1133 for (j = 0; j < lenx; j++) {
1134 a_elem[0] = a_i[aij];
1135 a_elem[1] = a_i[aij + 1];
1136
1137 x_elem[0] = head_x_i[jx];
1138 x_elem[1] = head_x_i[jx + 1];
1139 {
1140 prod[0] =
1141 (double) a_elem[0] * x_elem[0] -
1142 (double) a_elem[1] * x_elem[1];
1143 prod[1] =
1144 (double) a_elem[0] * x_elem[1] +
1145 (double) a_elem[1] * x_elem[0];
1146 }
1147 sum[0] = sum[0] + prod[0];
1148 sum[1] = sum[1] + prod[1];
1149 x_elem[0] = tail_x_i[jx];
1150 x_elem[1] = tail_x_i[jx + 1];
1151 {
1152 prod[0] =
1153 (double) a_elem[0] * x_elem[0] -
1154 (double) a_elem[1] * x_elem[1];
1155 prod[1] =
1156 (double) a_elem[0] * x_elem[1] +
1157 (double) a_elem[1] * x_elem[0];
1158 }
1159 sum2[0] = sum2[0] + prod[0];
1160 sum2[1] = sum2[1] + prod[1];
1161 aij += incaij;
1162 jx += incx;
1163 }
1164 {
1165 tmp1[0] =
1166 (double) sum[0] * alpha_i[0] -
1167 (double) sum[1] * alpha_i[1];
1168 tmp1[1] =
1169 (double) sum[0] * alpha_i[1] +
1170 (double) sum[1] * alpha_i[0];
1171 }
1172 {
1173 tmp2[0] =
1174 (double) sum2[0] * alpha_i[0] -
1175 (double) sum2[1] * alpha_i[1];
1176 tmp2[1] =
1177 (double) sum2[0] * alpha_i[1] +
1178 (double) sum2[1] * alpha_i[0];
1179 }
1180 tmp1[0] = tmp1[0] + tmp2[0];
1181 tmp1[1] = tmp1[1] + tmp2[1];
1182 y_elem[0] = y_i[iy];
1183 y_elem[1] = y_i[iy + 1];
1184 {
1185 tmp2[0] =
1186 (double) y_elem[0] * beta_i[0] -
1187 (double) y_elem[1] * beta_i[1];
1188 tmp2[1] =
1189 (double) y_elem[0] * beta_i[1] +
1190 (double) y_elem[1] * beta_i[0];
1191 }
1192 tmp1[0] = tmp1[0] + tmp2[0];
1193 tmp1[1] = tmp1[1] + tmp2[1];
1194 y_i[iy] = tmp1[0];
1195 y_i[iy + 1] = tmp1[1];
1196 ai += incai;
1197 iy += incy;
1198 }
1199 }
1200 }
1201
1202 }
1203 }
1204
1205
1206
1207 break;
1208 }
1209 case blas_prec_extra:{
1210
1211 int i, j;
1212 int iy, jx, kx, ky;
1213 int lenx, leny;
1214 int ai, aij;
1215 int incai, incaij;
1216
1217 const float *a_i = (float *) a;
1218 const float *head_x_i = (float *) head_x;
1219 const float *tail_x_i = (float *) tail_x;
1220 float *y_i = (float *) y;
1221 float *alpha_i = (float *) alpha;
1222 float *beta_i = (float *) beta;
1223 float a_elem[2];
1224 float x_elem[2];
1225 float y_elem[2];
1226 double head_prod[2], tail_prod[2];
1227 double head_sum[2], tail_sum[2];
1228 double head_sum2[2], tail_sum2[2];
1229 double head_tmp1[2], tail_tmp1[2];
1230 double head_tmp2[2], tail_tmp2[2];
1231 FPU_FIX_DECL;
1232
1233 /* all error calls */
1234 if (m < 0)
1235 BLAS_error(routine_name, -3, m, 0);
1236 else if (n <= 0)
1237 BLAS_error(routine_name, -4, n, 0);
1238 else if (incx == 0)
1239 BLAS_error(routine_name, -10, incx, 0);
1240 else if (incy == 0)
1241 BLAS_error(routine_name, -13, incy, 0);
1242
1243 if ((order == blas_rowmajor) && (trans == blas_no_trans)) {
1244 lenx = n;
1245 leny = m;
1246 incai = lda;
1247 incaij = 1;
1248 } else if ((order == blas_rowmajor) && (trans != blas_no_trans)) {
1249 lenx = m;
1250 leny = n;
1251 incai = 1;
1252 incaij = lda;
1253 } else if ((order == blas_colmajor) && (trans == blas_no_trans)) {
1254 lenx = n;
1255 leny = m;
1256 incai = 1;
1257 incaij = lda;
1258 } else { /* colmajor and blas_trans */
1259 lenx = m;
1260 leny = n;
1261 incai = lda;
1262 incaij = 1;
1263 }
1264
1265 if (lda < leny)
1266 BLAS_error(routine_name, -7, lda, NULL);
1267
1268 FPU_FIX_START;
1269
1270 incx *= 2;
1271 incy *= 2;
1272 incai *= 2;
1273 incaij *= 2;
1274
1275 if (incx > 0)
1276 kx = 0;
1277 else
1278 kx = (1 - lenx) * incx;
1279 if (incy > 0)
1280 ky = 0;
1281 else
1282 ky = (1 - leny) * incy;
1283
1284 /* No extra-precision needed for alpha = 0 */
1285 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0) {
1286 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
1287 iy = ky;
1288 for (i = 0; i < leny; i++) {
1289 y_i[iy] = 0.0;
1290 y_i[iy + 1] = 0.0;
1291 iy += incy;
1292 }
1293 } else if (!(beta_i[0] == 0.0 && beta_i[1] == 0.0)) {
1294 iy = ky;
1295 for (i = 0; i < leny; i++) {
1296 y_elem[0] = y_i[iy];
1297 y_elem[1] = y_i[iy + 1];
1298 {
1299 double head_e1, tail_e1;
1300 double d1;
1301 double d2;
1302 /* Real part */
1303 d1 = (double) y_elem[0] * beta_i[0];
1304 d2 = (double) -y_elem[1] * beta_i[1];
1305 {
1306 /* Compute double-double = double + double. */
1307 double e, t1, t2;
1308
1309 /* Knuth trick. */
1310 t1 = d1 + d2;
1311 e = t1 - d1;
1312 t2 = ((d2 - e) + (d1 - (t1 - e)));
1313
1314 /* The result is t1 + t2, after normalization. */
1315 head_e1 = t1 + t2;
1316 tail_e1 = t2 - (head_e1 - t1);
1317 }
1318 head_tmp1[0] = head_e1;
1319 tail_tmp1[0] = tail_e1;
1320 /* imaginary part */
1321 d1 = (double) y_elem[0] * beta_i[1];
1322 d2 = (double) y_elem[1] * beta_i[0];
1323 {
1324 /* Compute double-double = double + double. */
1325 double e, t1, t2;
1326
1327 /* Knuth trick. */
1328 t1 = d1 + d2;
1329 e = t1 - d1;
1330 t2 = ((d2 - e) + (d1 - (t1 - e)));
1331
1332 /* The result is t1 + t2, after normalization. */
1333 head_e1 = t1 + t2;
1334 tail_e1 = t2 - (head_e1 - t1);
1335 }
1336 head_tmp1[1] = head_e1;
1337 tail_tmp1[1] = tail_e1;
1338 }
1339 y_i[iy] = head_tmp1[0];
1340 y_i[iy + 1] = head_tmp1[1];
1341 iy += incy;
1342 }
1343 }
1344 } else { /* alpha != 0 */
1345 if (trans == blas_conj_trans) {
1346
1347 /* if beta = 0, we can save m multiplies:
1348 y = alpha*A*head_x + alpha*A*tail_x */
1349 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
1350 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
1351 /* save m more multiplies if alpha = 1 */
1352 ai = 0;
1353 iy = ky;
1354 for (i = 0; i < leny; i++) {
1355 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
1356 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
1357 0.0;
1358 aij = ai;
1359 jx = kx;
1360 for (j = 0; j < lenx; j++) {
1361 a_elem[0] = a_i[aij];
1362 a_elem[1] = a_i[aij + 1];
1363 a_elem[1] = -a_elem[1];
1364 x_elem[0] = head_x_i[jx];
1365 x_elem[1] = head_x_i[jx + 1];
1366 {
1367 double head_e1, tail_e1;
1368 double d1;
1369 double d2;
1370 /* Real part */
1371 d1 = (double) a_elem[0] * x_elem[0];
1372 d2 = (double) -a_elem[1] * x_elem[1];
1373 {
1374 /* Compute double-double = double + double. */
1375 double e, t1, t2;
1376
1377 /* Knuth trick. */
1378 t1 = d1 + d2;
1379 e = t1 - d1;
1380 t2 = ((d2 - e) + (d1 - (t1 - e)));
1381
1382 /* The result is t1 + t2, after normalization. */
1383 head_e1 = t1 + t2;
1384 tail_e1 = t2 - (head_e1 - t1);
1385 }
1386 head_prod[0] = head_e1;
1387 tail_prod[0] = tail_e1;
1388 /* imaginary part */
1389 d1 = (double) a_elem[0] * x_elem[1];
1390 d2 = (double) a_elem[1] * x_elem[0];
1391 {
1392 /* Compute double-double = double + double. */
1393 double e, t1, t2;
1394
1395 /* Knuth trick. */
1396 t1 = d1 + d2;
1397 e = t1 - d1;
1398 t2 = ((d2 - e) + (d1 - (t1 - e)));
1399
1400 /* The result is t1 + t2, after normalization. */
1401 head_e1 = t1 + t2;
1402 tail_e1 = t2 - (head_e1 - t1);
1403 }
1404 head_prod[1] = head_e1;
1405 tail_prod[1] = tail_e1;
1406 }
1407 {
1408 double head_t, tail_t;
1409 double head_a, tail_a;
1410 double head_b, tail_b;
1411 /* Real part */
1412 head_a = head_sum[0];
1413 tail_a = tail_sum[0];
1414 head_b = head_prod[0];
1415 tail_b = tail_prod[0];
1416 {
1417 /* Compute double-double = double-double + double-double. */
1418 double bv;
1419 double s1, s2, t1, t2;
1420
1421 /* Add two hi words. */
1422 s1 = head_a + head_b;
1423 bv = s1 - head_a;
1424 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1425
1426 /* Add two lo words. */
1427 t1 = tail_a + tail_b;
1428 bv = t1 - tail_a;
1429 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1430
1431 s2 += t1;
1432
1433 /* Renormalize (s1, s2) to (t1, s2) */
1434 t1 = s1 + s2;
1435 s2 = s2 - (t1 - s1);
1436
1437 t2 += s2;
1438
1439 /* Renormalize (t1, t2) */
1440 head_t = t1 + t2;
1441 tail_t = t2 - (head_t - t1);
1442 }
1443 head_sum[0] = head_t;
1444 tail_sum[0] = tail_t;
1445 /* Imaginary part */
1446 head_a = head_sum[1];
1447 tail_a = tail_sum[1];
1448 head_b = head_prod[1];
1449 tail_b = tail_prod[1];
1450 {
1451 /* Compute double-double = double-double + double-double. */
1452 double bv;
1453 double s1, s2, t1, t2;
1454
1455 /* Add two hi words. */
1456 s1 = head_a + head_b;
1457 bv = s1 - head_a;
1458 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1459
1460 /* Add two lo words. */
1461 t1 = tail_a + tail_b;
1462 bv = t1 - tail_a;
1463 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1464
1465 s2 += t1;
1466
1467 /* Renormalize (s1, s2) to (t1, s2) */
1468 t1 = s1 + s2;
1469 s2 = s2 - (t1 - s1);
1470
1471 t2 += s2;
1472
1473 /* Renormalize (t1, t2) */
1474 head_t = t1 + t2;
1475 tail_t = t2 - (head_t - t1);
1476 }
1477 head_sum[1] = head_t;
1478 tail_sum[1] = tail_t;
1479 }
1480 x_elem[0] = tail_x_i[jx];
1481 x_elem[1] = tail_x_i[jx + 1];
1482 {
1483 double head_e1, tail_e1;
1484 double d1;
1485 double d2;
1486 /* Real part */
1487 d1 = (double) a_elem[0] * x_elem[0];
1488 d2 = (double) -a_elem[1] * x_elem[1];
1489 {
1490 /* Compute double-double = double + double. */
1491 double e, t1, t2;
1492
1493 /* Knuth trick. */
1494 t1 = d1 + d2;
1495 e = t1 - d1;
1496 t2 = ((d2 - e) + (d1 - (t1 - e)));
1497
1498 /* The result is t1 + t2, after normalization. */
1499 head_e1 = t1 + t2;
1500 tail_e1 = t2 - (head_e1 - t1);
1501 }
1502 head_prod[0] = head_e1;
1503 tail_prod[0] = tail_e1;
1504 /* imaginary part */
1505 d1 = (double) a_elem[0] * x_elem[1];
1506 d2 = (double) a_elem[1] * x_elem[0];
1507 {
1508 /* Compute double-double = double + double. */
1509 double e, t1, t2;
1510
1511 /* Knuth trick. */
1512 t1 = d1 + d2;
1513 e = t1 - d1;
1514 t2 = ((d2 - e) + (d1 - (t1 - e)));
1515
1516 /* The result is t1 + t2, after normalization. */
1517 head_e1 = t1 + t2;
1518 tail_e1 = t2 - (head_e1 - t1);
1519 }
1520 head_prod[1] = head_e1;
1521 tail_prod[1] = tail_e1;
1522 }
1523 {
1524 double head_t, tail_t;
1525 double head_a, tail_a;
1526 double head_b, tail_b;
1527 /* Real part */
1528 head_a = head_sum2[0];
1529 tail_a = tail_sum2[0];
1530 head_b = head_prod[0];
1531 tail_b = tail_prod[0];
1532 {
1533 /* Compute double-double = double-double + double-double. */
1534 double bv;
1535 double s1, s2, t1, t2;
1536
1537 /* Add two hi words. */
1538 s1 = head_a + head_b;
1539 bv = s1 - head_a;
1540 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1541
1542 /* Add two lo words. */
1543 t1 = tail_a + tail_b;
1544 bv = t1 - tail_a;
1545 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1546
1547 s2 += t1;
1548
1549 /* Renormalize (s1, s2) to (t1, s2) */
1550 t1 = s1 + s2;
1551 s2 = s2 - (t1 - s1);
1552
1553 t2 += s2;
1554
1555 /* Renormalize (t1, t2) */
1556 head_t = t1 + t2;
1557 tail_t = t2 - (head_t - t1);
1558 }
1559 head_sum2[0] = head_t;
1560 tail_sum2[0] = tail_t;
1561 /* Imaginary part */
1562 head_a = head_sum2[1];
1563 tail_a = tail_sum2[1];
1564 head_b = head_prod[1];
1565 tail_b = tail_prod[1];
1566 {
1567 /* Compute double-double = double-double + double-double. */
1568 double bv;
1569 double s1, s2, t1, t2;
1570
1571 /* Add two hi words. */
1572 s1 = head_a + head_b;
1573 bv = s1 - head_a;
1574 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1575
1576 /* Add two lo words. */
1577 t1 = tail_a + tail_b;
1578 bv = t1 - tail_a;
1579 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1580
1581 s2 += t1;
1582
1583 /* Renormalize (s1, s2) to (t1, s2) */
1584 t1 = s1 + s2;
1585 s2 = s2 - (t1 - s1);
1586
1587 t2 += s2;
1588
1589 /* Renormalize (t1, t2) */
1590 head_t = t1 + t2;
1591 tail_t = t2 - (head_t - t1);
1592 }
1593 head_sum2[1] = head_t;
1594 tail_sum2[1] = tail_t;
1595 }
1596 aij += incaij;
1597 jx += incx;
1598 }
1599 {
1600 double head_t, tail_t;
1601 double head_a, tail_a;
1602 double head_b, tail_b;
1603 /* Real part */
1604 head_a = head_sum[0];
1605 tail_a = tail_sum[0];
1606 head_b = head_sum2[0];
1607 tail_b = tail_sum2[0];
1608 {
1609 /* Compute double-double = double-double + double-double. */
1610 double bv;
1611 double s1, s2, t1, t2;
1612
1613 /* Add two hi words. */
1614 s1 = head_a + head_b;
1615 bv = s1 - head_a;
1616 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1617
1618 /* Add two lo words. */
1619 t1 = tail_a + tail_b;
1620 bv = t1 - tail_a;
1621 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1622
1623 s2 += t1;
1624
1625 /* Renormalize (s1, s2) to (t1, s2) */
1626 t1 = s1 + s2;
1627 s2 = s2 - (t1 - s1);
1628
1629 t2 += s2;
1630
1631 /* Renormalize (t1, t2) */
1632 head_t = t1 + t2;
1633 tail_t = t2 - (head_t - t1);
1634 }
1635 head_sum[0] = head_t;
1636 tail_sum[0] = tail_t;
1637 /* Imaginary part */
1638 head_a = head_sum[1];
1639 tail_a = tail_sum[1];
1640 head_b = head_sum2[1];
1641 tail_b = tail_sum2[1];
1642 {
1643 /* Compute double-double = double-double + double-double. */
1644 double bv;
1645 double s1, s2, t1, t2;
1646
1647 /* Add two hi words. */
1648 s1 = head_a + head_b;
1649 bv = s1 - head_a;
1650 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1651
1652 /* Add two lo words. */
1653 t1 = tail_a + tail_b;
1654 bv = t1 - tail_a;
1655 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1656
1657 s2 += t1;
1658
1659 /* Renormalize (s1, s2) to (t1, s2) */
1660 t1 = s1 + s2;
1661 s2 = s2 - (t1 - s1);
1662
1663 t2 += s2;
1664
1665 /* Renormalize (t1, t2) */
1666 head_t = t1 + t2;
1667 tail_t = t2 - (head_t - t1);
1668 }
1669 head_sum[1] = head_t;
1670 tail_sum[1] = tail_t;
1671 }
1672 y_i[iy] = head_sum[0];
1673 y_i[iy + 1] = head_sum[1];
1674 ai += incai;
1675 iy += incy;
1676 } /* end for */
1677 } else { /* alpha != 1 */
1678 ai = 0;
1679 iy = ky;
1680 for (i = 0; i < leny; i++) {
1681 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
1682 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
1683 0.0;
1684 aij = ai;
1685 jx = kx;
1686 for (j = 0; j < lenx; j++) {
1687 a_elem[0] = a_i[aij];
1688 a_elem[1] = a_i[aij + 1];
1689 a_elem[1] = -a_elem[1];
1690 x_elem[0] = head_x_i[jx];
1691 x_elem[1] = head_x_i[jx + 1];
1692 {
1693 double head_e1, tail_e1;
1694 double d1;
1695 double d2;
1696 /* Real part */
1697 d1 = (double) a_elem[0] * x_elem[0];
1698 d2 = (double) -a_elem[1] * x_elem[1];
1699 {
1700 /* Compute double-double = double + double. */
1701 double e, t1, t2;
1702
1703 /* Knuth trick. */
1704 t1 = d1 + d2;
1705 e = t1 - d1;
1706 t2 = ((d2 - e) + (d1 - (t1 - e)));
1707
1708 /* The result is t1 + t2, after normalization. */
1709 head_e1 = t1 + t2;
1710 tail_e1 = t2 - (head_e1 - t1);
1711 }
1712 head_prod[0] = head_e1;
1713 tail_prod[0] = tail_e1;
1714 /* imaginary part */
1715 d1 = (double) a_elem[0] * x_elem[1];
1716 d2 = (double) a_elem[1] * x_elem[0];
1717 {
1718 /* Compute double-double = double + double. */
1719 double e, t1, t2;
1720
1721 /* Knuth trick. */
1722 t1 = d1 + d2;
1723 e = t1 - d1;
1724 t2 = ((d2 - e) + (d1 - (t1 - e)));
1725
1726 /* The result is t1 + t2, after normalization. */
1727 head_e1 = t1 + t2;
1728 tail_e1 = t2 - (head_e1 - t1);
1729 }
1730 head_prod[1] = head_e1;
1731 tail_prod[1] = tail_e1;
1732 }
1733 {
1734 double head_t, tail_t;
1735 double head_a, tail_a;
1736 double head_b, tail_b;
1737 /* Real part */
1738 head_a = head_sum[0];
1739 tail_a = tail_sum[0];
1740 head_b = head_prod[0];
1741 tail_b = tail_prod[0];
1742 {
1743 /* Compute double-double = double-double + double-double. */
1744 double bv;
1745 double s1, s2, t1, t2;
1746
1747 /* Add two hi words. */
1748 s1 = head_a + head_b;
1749 bv = s1 - head_a;
1750 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1751
1752 /* Add two lo words. */
1753 t1 = tail_a + tail_b;
1754 bv = t1 - tail_a;
1755 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1756
1757 s2 += t1;
1758
1759 /* Renormalize (s1, s2) to (t1, s2) */
1760 t1 = s1 + s2;
1761 s2 = s2 - (t1 - s1);
1762
1763 t2 += s2;
1764
1765 /* Renormalize (t1, t2) */
1766 head_t = t1 + t2;
1767 tail_t = t2 - (head_t - t1);
1768 }
1769 head_sum[0] = head_t;
1770 tail_sum[0] = tail_t;
1771 /* Imaginary part */
1772 head_a = head_sum[1];
1773 tail_a = tail_sum[1];
1774 head_b = head_prod[1];
1775 tail_b = tail_prod[1];
1776 {
1777 /* Compute double-double = double-double + double-double. */
1778 double bv;
1779 double s1, s2, t1, t2;
1780
1781 /* Add two hi words. */
1782 s1 = head_a + head_b;
1783 bv = s1 - head_a;
1784 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1785
1786 /* Add two lo words. */
1787 t1 = tail_a + tail_b;
1788 bv = t1 - tail_a;
1789 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1790
1791 s2 += t1;
1792
1793 /* Renormalize (s1, s2) to (t1, s2) */
1794 t1 = s1 + s2;
1795 s2 = s2 - (t1 - s1);
1796
1797 t2 += s2;
1798
1799 /* Renormalize (t1, t2) */
1800 head_t = t1 + t2;
1801 tail_t = t2 - (head_t - t1);
1802 }
1803 head_sum[1] = head_t;
1804 tail_sum[1] = tail_t;
1805 }
1806 x_elem[0] = tail_x_i[jx];
1807 x_elem[1] = tail_x_i[jx + 1];
1808 {
1809 double head_e1, tail_e1;
1810 double d1;
1811 double d2;
1812 /* Real part */
1813 d1 = (double) a_elem[0] * x_elem[0];
1814 d2 = (double) -a_elem[1] * x_elem[1];
1815 {
1816 /* Compute double-double = double + double. */
1817 double e, t1, t2;
1818
1819 /* Knuth trick. */
1820 t1 = d1 + d2;
1821 e = t1 - d1;
1822 t2 = ((d2 - e) + (d1 - (t1 - e)));
1823
1824 /* The result is t1 + t2, after normalization. */
1825 head_e1 = t1 + t2;
1826 tail_e1 = t2 - (head_e1 - t1);
1827 }
1828 head_prod[0] = head_e1;
1829 tail_prod[0] = tail_e1;
1830 /* imaginary part */
1831 d1 = (double) a_elem[0] * x_elem[1];
1832 d2 = (double) a_elem[1] * x_elem[0];
1833 {
1834 /* Compute double-double = double + double. */
1835 double e, t1, t2;
1836
1837 /* Knuth trick. */
1838 t1 = d1 + d2;
1839 e = t1 - d1;
1840 t2 = ((d2 - e) + (d1 - (t1 - e)));
1841
1842 /* The result is t1 + t2, after normalization. */
1843 head_e1 = t1 + t2;
1844 tail_e1 = t2 - (head_e1 - t1);
1845 }
1846 head_prod[1] = head_e1;
1847 tail_prod[1] = tail_e1;
1848 }
1849 {
1850 double head_t, tail_t;
1851 double head_a, tail_a;
1852 double head_b, tail_b;
1853 /* Real part */
1854 head_a = head_sum2[0];
1855 tail_a = tail_sum2[0];
1856 head_b = head_prod[0];
1857 tail_b = tail_prod[0];
1858 {
1859 /* Compute double-double = double-double + double-double. */
1860 double bv;
1861 double s1, s2, t1, t2;
1862
1863 /* Add two hi words. */
1864 s1 = head_a + head_b;
1865 bv = s1 - head_a;
1866 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1867
1868 /* Add two lo words. */
1869 t1 = tail_a + tail_b;
1870 bv = t1 - tail_a;
1871 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1872
1873 s2 += t1;
1874
1875 /* Renormalize (s1, s2) to (t1, s2) */
1876 t1 = s1 + s2;
1877 s2 = s2 - (t1 - s1);
1878
1879 t2 += s2;
1880
1881 /* Renormalize (t1, t2) */
1882 head_t = t1 + t2;
1883 tail_t = t2 - (head_t - t1);
1884 }
1885 head_sum2[0] = head_t;
1886 tail_sum2[0] = tail_t;
1887 /* Imaginary part */
1888 head_a = head_sum2[1];
1889 tail_a = tail_sum2[1];
1890 head_b = head_prod[1];
1891 tail_b = tail_prod[1];
1892 {
1893 /* Compute double-double = double-double + double-double. */
1894 double bv;
1895 double s1, s2, t1, t2;
1896
1897 /* Add two hi words. */
1898 s1 = head_a + head_b;
1899 bv = s1 - head_a;
1900 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1901
1902 /* Add two lo words. */
1903 t1 = tail_a + tail_b;
1904 bv = t1 - tail_a;
1905 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1906
1907 s2 += t1;
1908
1909 /* Renormalize (s1, s2) to (t1, s2) */
1910 t1 = s1 + s2;
1911 s2 = s2 - (t1 - s1);
1912
1913 t2 += s2;
1914
1915 /* Renormalize (t1, t2) */
1916 head_t = t1 + t2;
1917 tail_t = t2 - (head_t - t1);
1918 }
1919 head_sum2[1] = head_t;
1920 tail_sum2[1] = tail_t;
1921 }
1922 aij += incaij;
1923 jx += incx;
1924 }
1925 {
1926 double cd[2];
1927 cd[0] = (double) alpha_i[0];
1928 cd[1] = (double) alpha_i[1];
1929 {
1930 /* Compute complex-extra = complex-extra * complex-double. */
1931 double head_a0, tail_a0;
1932 double head_a1, tail_a1;
1933 double head_t1, tail_t1;
1934 double head_t2, tail_t2;
1935 head_a0 = head_sum[0];
1936 tail_a0 = tail_sum[0];
1937 head_a1 = head_sum[1];
1938 tail_a1 = tail_sum[1];
1939 /* real part */
1940 {
1941 /* Compute double-double = double-double * double. */
1942 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1943
1944 con = head_a0 * split;
1945 a11 = con - head_a0;
1946 a11 = con - a11;
1947 a21 = head_a0 - a11;
1948 con = cd[0] * split;
1949 b1 = con - cd[0];
1950 b1 = con - b1;
1951 b2 = cd[0] - b1;
1952
1953 c11 = head_a0 * cd[0];
1954 c21 =
1955 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1956
1957 c2 = tail_a0 * cd[0];
1958 t1 = c11 + c2;
1959 t2 = (c2 - (t1 - c11)) + c21;
1960
1961 head_t1 = t1 + t2;
1962 tail_t1 = t2 - (head_t1 - t1);
1963 }
1964 {
1965 /* Compute double-double = double-double * double. */
1966 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1967
1968 con = head_a1 * split;
1969 a11 = con - head_a1;
1970 a11 = con - a11;
1971 a21 = head_a1 - a11;
1972 con = cd[1] * split;
1973 b1 = con - cd[1];
1974 b1 = con - b1;
1975 b2 = cd[1] - b1;
1976
1977 c11 = head_a1 * cd[1];
1978 c21 =
1979 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1980
1981 c2 = tail_a1 * cd[1];
1982 t1 = c11 + c2;
1983 t2 = (c2 - (t1 - c11)) + c21;
1984
1985 head_t2 = t1 + t2;
1986 tail_t2 = t2 - (head_t2 - t1);
1987 }
1988 head_t2 = -head_t2;
1989 tail_t2 = -tail_t2;
1990 {
1991 /* Compute double-double = double-double + double-double. */
1992 double bv;
1993 double s1, s2, t1, t2;
1994
1995 /* Add two hi words. */
1996 s1 = head_t1 + head_t2;
1997 bv = s1 - head_t1;
1998 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1999
2000 /* Add two lo words. */
2001 t1 = tail_t1 + tail_t2;
2002 bv = t1 - tail_t1;
2003 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2004
2005 s2 += t1;
2006
2007 /* Renormalize (s1, s2) to (t1, s2) */
2008 t1 = s1 + s2;
2009 s2 = s2 - (t1 - s1);
2010
2011 t2 += s2;
2012
2013 /* Renormalize (t1, t2) */
2014 head_t1 = t1 + t2;
2015 tail_t1 = t2 - (head_t1 - t1);
2016 }
2017 head_tmp1[0] = head_t1;
2018 tail_tmp1[0] = tail_t1;
2019 /* imaginary part */
2020 {
2021 /* Compute double-double = double-double * double. */
2022 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2023
2024 con = head_a1 * split;
2025 a11 = con - head_a1;
2026 a11 = con - a11;
2027 a21 = head_a1 - a11;
2028 con = cd[0] * split;
2029 b1 = con - cd[0];
2030 b1 = con - b1;
2031 b2 = cd[0] - b1;
2032
2033 c11 = head_a1 * cd[0];
2034 c21 =
2035 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2036
2037 c2 = tail_a1 * cd[0];
2038 t1 = c11 + c2;
2039 t2 = (c2 - (t1 - c11)) + c21;
2040
2041 head_t1 = t1 + t2;
2042 tail_t1 = t2 - (head_t1 - t1);
2043 }
2044 {
2045 /* Compute double-double = double-double * double. */
2046 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2047
2048 con = head_a0 * split;
2049 a11 = con - head_a0;
2050 a11 = con - a11;
2051 a21 = head_a0 - a11;
2052 con = cd[1] * split;
2053 b1 = con - cd[1];
2054 b1 = con - b1;
2055 b2 = cd[1] - b1;
2056
2057 c11 = head_a0 * cd[1];
2058 c21 =
2059 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2060
2061 c2 = tail_a0 * cd[1];
2062 t1 = c11 + c2;
2063 t2 = (c2 - (t1 - c11)) + c21;
2064
2065 head_t2 = t1 + t2;
2066 tail_t2 = t2 - (head_t2 - t1);
2067 }
2068 {
2069 /* Compute double-double = double-double + double-double. */
2070 double bv;
2071 double s1, s2, t1, t2;
2072
2073 /* Add two hi words. */
2074 s1 = head_t1 + head_t2;
2075 bv = s1 - head_t1;
2076 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2077
2078 /* Add two lo words. */
2079 t1 = tail_t1 + tail_t2;
2080 bv = t1 - tail_t1;
2081 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2082
2083 s2 += t1;
2084
2085 /* Renormalize (s1, s2) to (t1, s2) */
2086 t1 = s1 + s2;
2087 s2 = s2 - (t1 - s1);
2088
2089 t2 += s2;
2090
2091 /* Renormalize (t1, t2) */
2092 head_t1 = t1 + t2;
2093 tail_t1 = t2 - (head_t1 - t1);
2094 }
2095 head_tmp1[1] = head_t1;
2096 tail_tmp1[1] = tail_t1;
2097 }
2098
2099 }
2100 {
2101 double cd[2];
2102 cd[0] = (double) alpha_i[0];
2103 cd[1] = (double) alpha_i[1];
2104 {
2105 /* Compute complex-extra = complex-extra * complex-double. */
2106 double head_a0, tail_a0;
2107 double head_a1, tail_a1;
2108 double head_t1, tail_t1;
2109 double head_t2, tail_t2;
2110 head_a0 = head_sum2[0];
2111 tail_a0 = tail_sum2[0];
2112 head_a1 = head_sum2[1];
2113 tail_a1 = tail_sum2[1];
2114 /* real part */
2115 {
2116 /* Compute double-double = double-double * double. */
2117 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2118
2119 con = head_a0 * split;
2120 a11 = con - head_a0;
2121 a11 = con - a11;
2122 a21 = head_a0 - a11;
2123 con = cd[0] * split;
2124 b1 = con - cd[0];
2125 b1 = con - b1;
2126 b2 = cd[0] - b1;
2127
2128 c11 = head_a0 * cd[0];
2129 c21 =
2130 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2131
2132 c2 = tail_a0 * cd[0];
2133 t1 = c11 + c2;
2134 t2 = (c2 - (t1 - c11)) + c21;
2135
2136 head_t1 = t1 + t2;
2137 tail_t1 = t2 - (head_t1 - t1);
2138 }
2139 {
2140 /* Compute double-double = double-double * double. */
2141 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2142
2143 con = head_a1 * split;
2144 a11 = con - head_a1;
2145 a11 = con - a11;
2146 a21 = head_a1 - a11;
2147 con = cd[1] * split;
2148 b1 = con - cd[1];
2149 b1 = con - b1;
2150 b2 = cd[1] - b1;
2151
2152 c11 = head_a1 * cd[1];
2153 c21 =
2154 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2155
2156 c2 = tail_a1 * cd[1];
2157 t1 = c11 + c2;
2158 t2 = (c2 - (t1 - c11)) + c21;
2159
2160 head_t2 = t1 + t2;
2161 tail_t2 = t2 - (head_t2 - t1);
2162 }
2163 head_t2 = -head_t2;
2164 tail_t2 = -tail_t2;
2165 {
2166 /* Compute double-double = double-double + double-double. */
2167 double bv;
2168 double s1, s2, t1, t2;
2169
2170 /* Add two hi words. */
2171 s1 = head_t1 + head_t2;
2172 bv = s1 - head_t1;
2173 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2174
2175 /* Add two lo words. */
2176 t1 = tail_t1 + tail_t2;
2177 bv = t1 - tail_t1;
2178 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2179
2180 s2 += t1;
2181
2182 /* Renormalize (s1, s2) to (t1, s2) */
2183 t1 = s1 + s2;
2184 s2 = s2 - (t1 - s1);
2185
2186 t2 += s2;
2187
2188 /* Renormalize (t1, t2) */
2189 head_t1 = t1 + t2;
2190 tail_t1 = t2 - (head_t1 - t1);
2191 }
2192 head_tmp2[0] = head_t1;
2193 tail_tmp2[0] = tail_t1;
2194 /* imaginary part */
2195 {
2196 /* Compute double-double = double-double * double. */
2197 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2198
2199 con = head_a1 * split;
2200 a11 = con - head_a1;
2201 a11 = con - a11;
2202 a21 = head_a1 - a11;
2203 con = cd[0] * split;
2204 b1 = con - cd[0];
2205 b1 = con - b1;
2206 b2 = cd[0] - b1;
2207
2208 c11 = head_a1 * cd[0];
2209 c21 =
2210 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2211
2212 c2 = tail_a1 * cd[0];
2213 t1 = c11 + c2;
2214 t2 = (c2 - (t1 - c11)) + c21;
2215
2216 head_t1 = t1 + t2;
2217 tail_t1 = t2 - (head_t1 - t1);
2218 }
2219 {
2220 /* Compute double-double = double-double * double. */
2221 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2222
2223 con = head_a0 * split;
2224 a11 = con - head_a0;
2225 a11 = con - a11;
2226 a21 = head_a0 - a11;
2227 con = cd[1] * split;
2228 b1 = con - cd[1];
2229 b1 = con - b1;
2230 b2 = cd[1] - b1;
2231
2232 c11 = head_a0 * cd[1];
2233 c21 =
2234 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2235
2236 c2 = tail_a0 * cd[1];
2237 t1 = c11 + c2;
2238 t2 = (c2 - (t1 - c11)) + c21;
2239
2240 head_t2 = t1 + t2;
2241 tail_t2 = t2 - (head_t2 - t1);
2242 }
2243 {
2244 /* Compute double-double = double-double + double-double. */
2245 double bv;
2246 double s1, s2, t1, t2;
2247
2248 /* Add two hi words. */
2249 s1 = head_t1 + head_t2;
2250 bv = s1 - head_t1;
2251 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2252
2253 /* Add two lo words. */
2254 t1 = tail_t1 + tail_t2;
2255 bv = t1 - tail_t1;
2256 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2257
2258 s2 += t1;
2259
2260 /* Renormalize (s1, s2) to (t1, s2) */
2261 t1 = s1 + s2;
2262 s2 = s2 - (t1 - s1);
2263
2264 t2 += s2;
2265
2266 /* Renormalize (t1, t2) */
2267 head_t1 = t1 + t2;
2268 tail_t1 = t2 - (head_t1 - t1);
2269 }
2270 head_tmp2[1] = head_t1;
2271 tail_tmp2[1] = tail_t1;
2272 }
2273
2274 }
2275 {
2276 double head_t, tail_t;
2277 double head_a, tail_a;
2278 double head_b, tail_b;
2279 /* Real part */
2280 head_a = head_tmp1[0];
2281 tail_a = tail_tmp1[0];
2282 head_b = head_tmp2[0];
2283 tail_b = tail_tmp2[0];
2284 {
2285 /* Compute double-double = double-double + double-double. */
2286 double bv;
2287 double s1, s2, t1, t2;
2288
2289 /* Add two hi words. */
2290 s1 = head_a + head_b;
2291 bv = s1 - head_a;
2292 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2293
2294 /* Add two lo words. */
2295 t1 = tail_a + tail_b;
2296 bv = t1 - tail_a;
2297 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2298
2299 s2 += t1;
2300
2301 /* Renormalize (s1, s2) to (t1, s2) */
2302 t1 = s1 + s2;
2303 s2 = s2 - (t1 - s1);
2304
2305 t2 += s2;
2306
2307 /* Renormalize (t1, t2) */
2308 head_t = t1 + t2;
2309 tail_t = t2 - (head_t - t1);
2310 }
2311 head_tmp1[0] = head_t;
2312 tail_tmp1[0] = tail_t;
2313 /* Imaginary part */
2314 head_a = head_tmp1[1];
2315 tail_a = tail_tmp1[1];
2316 head_b = head_tmp2[1];
2317 tail_b = tail_tmp2[1];
2318 {
2319 /* Compute double-double = double-double + double-double. */
2320 double bv;
2321 double s1, s2, t1, t2;
2322
2323 /* Add two hi words. */
2324 s1 = head_a + head_b;
2325 bv = s1 - head_a;
2326 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2327
2328 /* Add two lo words. */
2329 t1 = tail_a + tail_b;
2330 bv = t1 - tail_a;
2331 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2332
2333 s2 += t1;
2334
2335 /* Renormalize (s1, s2) to (t1, s2) */
2336 t1 = s1 + s2;
2337 s2 = s2 - (t1 - s1);
2338
2339 t2 += s2;
2340
2341 /* Renormalize (t1, t2) */
2342 head_t = t1 + t2;
2343 tail_t = t2 - (head_t - t1);
2344 }
2345 head_tmp1[1] = head_t;
2346 tail_tmp1[1] = tail_t;
2347 }
2348 y_i[iy] = head_tmp1[0];
2349 y_i[iy + 1] = head_tmp1[1];
2350 ai += incai;
2351 iy += incy;
2352 }
2353 }
2354 } else { /* beta != 0 */
2355 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
2356 /* save m multiplies if alpha = 1 */
2357 ai = 0;
2358 iy = ky;
2359 for (i = 0; i < leny; i++) {
2360 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;;
2361 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
2362 0.0;;
2363 aij = ai;
2364 jx = kx;
2365 for (j = 0; j < lenx; j++) {
2366 a_elem[0] = a_i[aij];
2367 a_elem[1] = a_i[aij + 1];
2368 a_elem[1] = -a_elem[1];
2369 x_elem[0] = head_x_i[jx];
2370 x_elem[1] = head_x_i[jx + 1];
2371 {
2372 double head_e1, tail_e1;
2373 double d1;
2374 double d2;
2375 /* Real part */
2376 d1 = (double) a_elem[0] * x_elem[0];
2377 d2 = (double) -a_elem[1] * x_elem[1];
2378 {
2379 /* Compute double-double = double + double. */
2380 double e, t1, t2;
2381
2382 /* Knuth trick. */
2383 t1 = d1 + d2;
2384 e = t1 - d1;
2385 t2 = ((d2 - e) + (d1 - (t1 - e)));
2386
2387 /* The result is t1 + t2, after normalization. */
2388 head_e1 = t1 + t2;
2389 tail_e1 = t2 - (head_e1 - t1);
2390 }
2391 head_prod[0] = head_e1;
2392 tail_prod[0] = tail_e1;
2393 /* imaginary part */
2394 d1 = (double) a_elem[0] * x_elem[1];
2395 d2 = (double) a_elem[1] * x_elem[0];
2396 {
2397 /* Compute double-double = double + double. */
2398 double e, t1, t2;
2399
2400 /* Knuth trick. */
2401 t1 = d1 + d2;
2402 e = t1 - d1;
2403 t2 = ((d2 - e) + (d1 - (t1 - e)));
2404
2405 /* The result is t1 + t2, after normalization. */
2406 head_e1 = t1 + t2;
2407 tail_e1 = t2 - (head_e1 - t1);
2408 }
2409 head_prod[1] = head_e1;
2410 tail_prod[1] = tail_e1;
2411 }
2412 {
2413 double head_t, tail_t;
2414 double head_a, tail_a;
2415 double head_b, tail_b;
2416 /* Real part */
2417 head_a = head_sum[0];
2418 tail_a = tail_sum[0];
2419 head_b = head_prod[0];
2420 tail_b = tail_prod[0];
2421 {
2422 /* Compute double-double = double-double + double-double. */
2423 double bv;
2424 double s1, s2, t1, t2;
2425
2426 /* Add two hi words. */
2427 s1 = head_a + head_b;
2428 bv = s1 - head_a;
2429 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2430
2431 /* Add two lo words. */
2432 t1 = tail_a + tail_b;
2433 bv = t1 - tail_a;
2434 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2435
2436 s2 += t1;
2437
2438 /* Renormalize (s1, s2) to (t1, s2) */
2439 t1 = s1 + s2;
2440 s2 = s2 - (t1 - s1);
2441
2442 t2 += s2;
2443
2444 /* Renormalize (t1, t2) */
2445 head_t = t1 + t2;
2446 tail_t = t2 - (head_t - t1);
2447 }
2448 head_sum[0] = head_t;
2449 tail_sum[0] = tail_t;
2450 /* Imaginary part */
2451 head_a = head_sum[1];
2452 tail_a = tail_sum[1];
2453 head_b = head_prod[1];
2454 tail_b = tail_prod[1];
2455 {
2456 /* Compute double-double = double-double + double-double. */
2457 double bv;
2458 double s1, s2, t1, t2;
2459
2460 /* Add two hi words. */
2461 s1 = head_a + head_b;
2462 bv = s1 - head_a;
2463 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2464
2465 /* Add two lo words. */
2466 t1 = tail_a + tail_b;
2467 bv = t1 - tail_a;
2468 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2469
2470 s2 += t1;
2471
2472 /* Renormalize (s1, s2) to (t1, s2) */
2473 t1 = s1 + s2;
2474 s2 = s2 - (t1 - s1);
2475
2476 t2 += s2;
2477
2478 /* Renormalize (t1, t2) */
2479 head_t = t1 + t2;
2480 tail_t = t2 - (head_t - t1);
2481 }
2482 head_sum[1] = head_t;
2483 tail_sum[1] = tail_t;
2484 }
2485 x_elem[0] = tail_x_i[jx];
2486 x_elem[1] = tail_x_i[jx + 1];
2487 {
2488 double head_e1, tail_e1;
2489 double d1;
2490 double d2;
2491 /* Real part */
2492 d1 = (double) a_elem[0] * x_elem[0];
2493 d2 = (double) -a_elem[1] * x_elem[1];
2494 {
2495 /* Compute double-double = double + double. */
2496 double e, t1, t2;
2497
2498 /* Knuth trick. */
2499 t1 = d1 + d2;
2500 e = t1 - d1;
2501 t2 = ((d2 - e) + (d1 - (t1 - e)));
2502
2503 /* The result is t1 + t2, after normalization. */
2504 head_e1 = t1 + t2;
2505 tail_e1 = t2 - (head_e1 - t1);
2506 }
2507 head_prod[0] = head_e1;
2508 tail_prod[0] = tail_e1;
2509 /* imaginary part */
2510 d1 = (double) a_elem[0] * x_elem[1];
2511 d2 = (double) a_elem[1] * x_elem[0];
2512 {
2513 /* Compute double-double = double + double. */
2514 double e, t1, t2;
2515
2516 /* Knuth trick. */
2517 t1 = d1 + d2;
2518 e = t1 - d1;
2519 t2 = ((d2 - e) + (d1 - (t1 - e)));
2520
2521 /* The result is t1 + t2, after normalization. */
2522 head_e1 = t1 + t2;
2523 tail_e1 = t2 - (head_e1 - t1);
2524 }
2525 head_prod[1] = head_e1;
2526 tail_prod[1] = tail_e1;
2527 }
2528 {
2529 double head_t, tail_t;
2530 double head_a, tail_a;
2531 double head_b, tail_b;
2532 /* Real part */
2533 head_a = head_sum2[0];
2534 tail_a = tail_sum2[0];
2535 head_b = head_prod[0];
2536 tail_b = tail_prod[0];
2537 {
2538 /* Compute double-double = double-double + double-double. */
2539 double bv;
2540 double s1, s2, t1, t2;
2541
2542 /* Add two hi words. */
2543 s1 = head_a + head_b;
2544 bv = s1 - head_a;
2545 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2546
2547 /* Add two lo words. */
2548 t1 = tail_a + tail_b;
2549 bv = t1 - tail_a;
2550 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2551
2552 s2 += t1;
2553
2554 /* Renormalize (s1, s2) to (t1, s2) */
2555 t1 = s1 + s2;
2556 s2 = s2 - (t1 - s1);
2557
2558 t2 += s2;
2559
2560 /* Renormalize (t1, t2) */
2561 head_t = t1 + t2;
2562 tail_t = t2 - (head_t - t1);
2563 }
2564 head_sum2[0] = head_t;
2565 tail_sum2[0] = tail_t;
2566 /* Imaginary part */
2567 head_a = head_sum2[1];
2568 tail_a = tail_sum2[1];
2569 head_b = head_prod[1];
2570 tail_b = tail_prod[1];
2571 {
2572 /* Compute double-double = double-double + double-double. */
2573 double bv;
2574 double s1, s2, t1, t2;
2575
2576 /* Add two hi words. */
2577 s1 = head_a + head_b;
2578 bv = s1 - head_a;
2579 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2580
2581 /* Add two lo words. */
2582 t1 = tail_a + tail_b;
2583 bv = t1 - tail_a;
2584 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2585
2586 s2 += t1;
2587
2588 /* Renormalize (s1, s2) to (t1, s2) */
2589 t1 = s1 + s2;
2590 s2 = s2 - (t1 - s1);
2591
2592 t2 += s2;
2593
2594 /* Renormalize (t1, t2) */
2595 head_t = t1 + t2;
2596 tail_t = t2 - (head_t - t1);
2597 }
2598 head_sum2[1] = head_t;
2599 tail_sum2[1] = tail_t;
2600 }
2601 aij += incaij;
2602 jx += incx;
2603 }
2604 {
2605 double head_t, tail_t;
2606 double head_a, tail_a;
2607 double head_b, tail_b;
2608 /* Real part */
2609 head_a = head_sum[0];
2610 tail_a = tail_sum[0];
2611 head_b = head_sum2[0];
2612 tail_b = tail_sum2[0];
2613 {
2614 /* Compute double-double = double-double + double-double. */
2615 double bv;
2616 double s1, s2, t1, t2;
2617
2618 /* Add two hi words. */
2619 s1 = head_a + head_b;
2620 bv = s1 - head_a;
2621 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2622
2623 /* Add two lo words. */
2624 t1 = tail_a + tail_b;
2625 bv = t1 - tail_a;
2626 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2627
2628 s2 += t1;
2629
2630 /* Renormalize (s1, s2) to (t1, s2) */
2631 t1 = s1 + s2;
2632 s2 = s2 - (t1 - s1);
2633
2634 t2 += s2;
2635
2636 /* Renormalize (t1, t2) */
2637 head_t = t1 + t2;
2638 tail_t = t2 - (head_t - t1);
2639 }
2640 head_sum[0] = head_t;
2641 tail_sum[0] = tail_t;
2642 /* Imaginary part */
2643 head_a = head_sum[1];
2644 tail_a = tail_sum[1];
2645 head_b = head_sum2[1];
2646 tail_b = tail_sum2[1];
2647 {
2648 /* Compute double-double = double-double + double-double. */
2649 double bv;
2650 double s1, s2, t1, t2;
2651
2652 /* Add two hi words. */
2653 s1 = head_a + head_b;
2654 bv = s1 - head_a;
2655 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2656
2657 /* Add two lo words. */
2658 t1 = tail_a + tail_b;
2659 bv = t1 - tail_a;
2660 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2661
2662 s2 += t1;
2663
2664 /* Renormalize (s1, s2) to (t1, s2) */
2665 t1 = s1 + s2;
2666 s2 = s2 - (t1 - s1);
2667
2668 t2 += s2;
2669
2670 /* Renormalize (t1, t2) */
2671 head_t = t1 + t2;
2672 tail_t = t2 - (head_t - t1);
2673 }
2674 head_sum[1] = head_t;
2675 tail_sum[1] = tail_t;
2676 }
2677 y_elem[0] = y_i[iy];
2678 y_elem[1] = y_i[iy + 1];
2679 {
2680 double head_e1, tail_e1;
2681 double d1;
2682 double d2;
2683 /* Real part */
2684 d1 = (double) y_elem[0] * beta_i[0];
2685 d2 = (double) -y_elem[1] * beta_i[1];
2686 {
2687 /* Compute double-double = double + double. */
2688 double e, t1, t2;
2689
2690 /* Knuth trick. */
2691 t1 = d1 + d2;
2692 e = t1 - d1;
2693 t2 = ((d2 - e) + (d1 - (t1 - e)));
2694
2695 /* The result is t1 + t2, after normalization. */
2696 head_e1 = t1 + t2;
2697 tail_e1 = t2 - (head_e1 - t1);
2698 }
2699 head_tmp1[0] = head_e1;
2700 tail_tmp1[0] = tail_e1;
2701 /* imaginary part */
2702 d1 = (double) y_elem[0] * beta_i[1];
2703 d2 = (double) y_elem[1] * beta_i[0];
2704 {
2705 /* Compute double-double = double + double. */
2706 double e, t1, t2;
2707
2708 /* Knuth trick. */
2709 t1 = d1 + d2;
2710 e = t1 - d1;
2711 t2 = ((d2 - e) + (d1 - (t1 - e)));
2712
2713 /* The result is t1 + t2, after normalization. */
2714 head_e1 = t1 + t2;
2715 tail_e1 = t2 - (head_e1 - t1);
2716 }
2717 head_tmp1[1] = head_e1;
2718 tail_tmp1[1] = tail_e1;
2719 }
2720 {
2721 double head_t, tail_t;
2722 double head_a, tail_a;
2723 double head_b, tail_b;
2724 /* Real part */
2725 head_a = head_sum[0];
2726 tail_a = tail_sum[0];
2727 head_b = head_tmp1[0];
2728 tail_b = tail_tmp1[0];
2729 {
2730 /* Compute double-double = double-double + double-double. */
2731 double bv;
2732 double s1, s2, t1, t2;
2733
2734 /* Add two hi words. */
2735 s1 = head_a + head_b;
2736 bv = s1 - head_a;
2737 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2738
2739 /* Add two lo words. */
2740 t1 = tail_a + tail_b;
2741 bv = t1 - tail_a;
2742 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2743
2744 s2 += t1;
2745
2746 /* Renormalize (s1, s2) to (t1, s2) */
2747 t1 = s1 + s2;
2748 s2 = s2 - (t1 - s1);
2749
2750 t2 += s2;
2751
2752 /* Renormalize (t1, t2) */
2753 head_t = t1 + t2;
2754 tail_t = t2 - (head_t - t1);
2755 }
2756 head_tmp2[0] = head_t;
2757 tail_tmp2[0] = tail_t;
2758 /* Imaginary part */
2759 head_a = head_sum[1];
2760 tail_a = tail_sum[1];
2761 head_b = head_tmp1[1];
2762 tail_b = tail_tmp1[1];
2763 {
2764 /* Compute double-double = double-double + double-double. */
2765 double bv;
2766 double s1, s2, t1, t2;
2767
2768 /* Add two hi words. */
2769 s1 = head_a + head_b;
2770 bv = s1 - head_a;
2771 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2772
2773 /* Add two lo words. */
2774 t1 = tail_a + tail_b;
2775 bv = t1 - tail_a;
2776 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2777
2778 s2 += t1;
2779
2780 /* Renormalize (s1, s2) to (t1, s2) */
2781 t1 = s1 + s2;
2782 s2 = s2 - (t1 - s1);
2783
2784 t2 += s2;
2785
2786 /* Renormalize (t1, t2) */
2787 head_t = t1 + t2;
2788 tail_t = t2 - (head_t - t1);
2789 }
2790 head_tmp2[1] = head_t;
2791 tail_tmp2[1] = tail_t;
2792 }
2793 y_i[iy] = head_tmp2[0];
2794 y_i[iy + 1] = head_tmp2[1];
2795 ai += incai;
2796 iy += incy;
2797 }
2798 } else { /* alpha != 1, the most general form:
2799 y = alpha*A*head_x + alpha*A*tail_x + beta*y */
2800 ai = 0;
2801 iy = ky;
2802 for (i = 0; i < leny; i++) {
2803 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;;
2804 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
2805 0.0;;
2806 aij = ai;
2807 jx = kx;
2808 for (j = 0; j < lenx; j++) {
2809 a_elem[0] = a_i[aij];
2810 a_elem[1] = a_i[aij + 1];
2811 a_elem[1] = -a_elem[1];
2812 x_elem[0] = head_x_i[jx];
2813 x_elem[1] = head_x_i[jx + 1];
2814 {
2815 double head_e1, tail_e1;
2816 double d1;
2817 double d2;
2818 /* Real part */
2819 d1 = (double) a_elem[0] * x_elem[0];
2820 d2 = (double) -a_elem[1] * x_elem[1];
2821 {
2822 /* Compute double-double = double + double. */
2823 double e, t1, t2;
2824
2825 /* Knuth trick. */
2826 t1 = d1 + d2;
2827 e = t1 - d1;
2828 t2 = ((d2 - e) + (d1 - (t1 - e)));
2829
2830 /* The result is t1 + t2, after normalization. */
2831 head_e1 = t1 + t2;
2832 tail_e1 = t2 - (head_e1 - t1);
2833 }
2834 head_prod[0] = head_e1;
2835 tail_prod[0] = tail_e1;
2836 /* imaginary part */
2837 d1 = (double) a_elem[0] * x_elem[1];
2838 d2 = (double) a_elem[1] * x_elem[0];
2839 {
2840 /* Compute double-double = double + double. */
2841 double e, t1, t2;
2842
2843 /* Knuth trick. */
2844 t1 = d1 + d2;
2845 e = t1 - d1;
2846 t2 = ((d2 - e) + (d1 - (t1 - e)));
2847
2848 /* The result is t1 + t2, after normalization. */
2849 head_e1 = t1 + t2;
2850 tail_e1 = t2 - (head_e1 - t1);
2851 }
2852 head_prod[1] = head_e1;
2853 tail_prod[1] = tail_e1;
2854 }
2855 {
2856 double head_t, tail_t;
2857 double head_a, tail_a;
2858 double head_b, tail_b;
2859 /* Real part */
2860 head_a = head_sum[0];
2861 tail_a = tail_sum[0];
2862 head_b = head_prod[0];
2863 tail_b = tail_prod[0];
2864 {
2865 /* Compute double-double = double-double + double-double. */
2866 double bv;
2867 double s1, s2, t1, t2;
2868
2869 /* Add two hi words. */
2870 s1 = head_a + head_b;
2871 bv = s1 - head_a;
2872 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2873
2874 /* Add two lo words. */
2875 t1 = tail_a + tail_b;
2876 bv = t1 - tail_a;
2877 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2878
2879 s2 += t1;
2880
2881 /* Renormalize (s1, s2) to (t1, s2) */
2882 t1 = s1 + s2;
2883 s2 = s2 - (t1 - s1);
2884
2885 t2 += s2;
2886
2887 /* Renormalize (t1, t2) */
2888 head_t = t1 + t2;
2889 tail_t = t2 - (head_t - t1);
2890 }
2891 head_sum[0] = head_t;
2892 tail_sum[0] = tail_t;
2893 /* Imaginary part */
2894 head_a = head_sum[1];
2895 tail_a = tail_sum[1];
2896 head_b = head_prod[1];
2897 tail_b = tail_prod[1];
2898 {
2899 /* Compute double-double = double-double + double-double. */
2900 double bv;
2901 double s1, s2, t1, t2;
2902
2903 /* Add two hi words. */
2904 s1 = head_a + head_b;
2905 bv = s1 - head_a;
2906 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2907
2908 /* Add two lo words. */
2909 t1 = tail_a + tail_b;
2910 bv = t1 - tail_a;
2911 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2912
2913 s2 += t1;
2914
2915 /* Renormalize (s1, s2) to (t1, s2) */
2916 t1 = s1 + s2;
2917 s2 = s2 - (t1 - s1);
2918
2919 t2 += s2;
2920
2921 /* Renormalize (t1, t2) */
2922 head_t = t1 + t2;
2923 tail_t = t2 - (head_t - t1);
2924 }
2925 head_sum[1] = head_t;
2926 tail_sum[1] = tail_t;
2927 }
2928 x_elem[0] = tail_x_i[jx];
2929 x_elem[1] = tail_x_i[jx + 1];
2930 {
2931 double head_e1, tail_e1;
2932 double d1;
2933 double d2;
2934 /* Real part */
2935 d1 = (double) a_elem[0] * x_elem[0];
2936 d2 = (double) -a_elem[1] * x_elem[1];
2937 {
2938 /* Compute double-double = double + double. */
2939 double e, t1, t2;
2940
2941 /* Knuth trick. */
2942 t1 = d1 + d2;
2943 e = t1 - d1;
2944 t2 = ((d2 - e) + (d1 - (t1 - e)));
2945
2946 /* The result is t1 + t2, after normalization. */
2947 head_e1 = t1 + t2;
2948 tail_e1 = t2 - (head_e1 - t1);
2949 }
2950 head_prod[0] = head_e1;
2951 tail_prod[0] = tail_e1;
2952 /* imaginary part */
2953 d1 = (double) a_elem[0] * x_elem[1];
2954 d2 = (double) a_elem[1] * x_elem[0];
2955 {
2956 /* Compute double-double = double + double. */
2957 double e, t1, t2;
2958
2959 /* Knuth trick. */
2960 t1 = d1 + d2;
2961 e = t1 - d1;
2962 t2 = ((d2 - e) + (d1 - (t1 - e)));
2963
2964 /* The result is t1 + t2, after normalization. */
2965 head_e1 = t1 + t2;
2966 tail_e1 = t2 - (head_e1 - t1);
2967 }
2968 head_prod[1] = head_e1;
2969 tail_prod[1] = tail_e1;
2970 }
2971 {
2972 double head_t, tail_t;
2973 double head_a, tail_a;
2974 double head_b, tail_b;
2975 /* Real part */
2976 head_a = head_sum2[0];
2977 tail_a = tail_sum2[0];
2978 head_b = head_prod[0];
2979 tail_b = tail_prod[0];
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_sum2[0] = head_t;
3008 tail_sum2[0] = tail_t;
3009 /* Imaginary part */
3010 head_a = head_sum2[1];
3011 tail_a = tail_sum2[1];
3012 head_b = head_prod[1];
3013 tail_b = tail_prod[1];
3014 {
3015 /* Compute double-double = double-double + double-double. */
3016 double bv;
3017 double s1, s2, t1, t2;
3018
3019 /* Add two hi words. */
3020 s1 = head_a + head_b;
3021 bv = s1 - head_a;
3022 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3023
3024 /* Add two lo words. */
3025 t1 = tail_a + tail_b;
3026 bv = t1 - tail_a;
3027 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3028
3029 s2 += t1;
3030
3031 /* Renormalize (s1, s2) to (t1, s2) */
3032 t1 = s1 + s2;
3033 s2 = s2 - (t1 - s1);
3034
3035 t2 += s2;
3036
3037 /* Renormalize (t1, t2) */
3038 head_t = t1 + t2;
3039 tail_t = t2 - (head_t - t1);
3040 }
3041 head_sum2[1] = head_t;
3042 tail_sum2[1] = tail_t;
3043 }
3044 aij += incaij;
3045 jx += incx;
3046 }
3047 {
3048 double cd[2];
3049 cd[0] = (double) alpha_i[0];
3050 cd[1] = (double) alpha_i[1];
3051 {
3052 /* Compute complex-extra = complex-extra * complex-double. */
3053 double head_a0, tail_a0;
3054 double head_a1, tail_a1;
3055 double head_t1, tail_t1;
3056 double head_t2, tail_t2;
3057 head_a0 = head_sum[0];
3058 tail_a0 = tail_sum[0];
3059 head_a1 = head_sum[1];
3060 tail_a1 = tail_sum[1];
3061 /* real part */
3062 {
3063 /* Compute double-double = double-double * double. */
3064 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3065
3066 con = head_a0 * split;
3067 a11 = con - head_a0;
3068 a11 = con - a11;
3069 a21 = head_a0 - a11;
3070 con = cd[0] * split;
3071 b1 = con - cd[0];
3072 b1 = con - b1;
3073 b2 = cd[0] - b1;
3074
3075 c11 = head_a0 * cd[0];
3076 c21 =
3077 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3078
3079 c2 = tail_a0 * cd[0];
3080 t1 = c11 + c2;
3081 t2 = (c2 - (t1 - c11)) + c21;
3082
3083 head_t1 = t1 + t2;
3084 tail_t1 = t2 - (head_t1 - t1);
3085 }
3086 {
3087 /* Compute double-double = double-double * double. */
3088 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3089
3090 con = head_a1 * split;
3091 a11 = con - head_a1;
3092 a11 = con - a11;
3093 a21 = head_a1 - a11;
3094 con = cd[1] * split;
3095 b1 = con - cd[1];
3096 b1 = con - b1;
3097 b2 = cd[1] - b1;
3098
3099 c11 = head_a1 * cd[1];
3100 c21 =
3101 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3102
3103 c2 = tail_a1 * cd[1];
3104 t1 = c11 + c2;
3105 t2 = (c2 - (t1 - c11)) + c21;
3106
3107 head_t2 = t1 + t2;
3108 tail_t2 = t2 - (head_t2 - t1);
3109 }
3110 head_t2 = -head_t2;
3111 tail_t2 = -tail_t2;
3112 {
3113 /* Compute double-double = double-double + double-double. */
3114 double bv;
3115 double s1, s2, t1, t2;
3116
3117 /* Add two hi words. */
3118 s1 = head_t1 + head_t2;
3119 bv = s1 - head_t1;
3120 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
3121
3122 /* Add two lo words. */
3123 t1 = tail_t1 + tail_t2;
3124 bv = t1 - tail_t1;
3125 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
3126
3127 s2 += t1;
3128
3129 /* Renormalize (s1, s2) to (t1, s2) */
3130 t1 = s1 + s2;
3131 s2 = s2 - (t1 - s1);
3132
3133 t2 += s2;
3134
3135 /* Renormalize (t1, t2) */
3136 head_t1 = t1 + t2;
3137 tail_t1 = t2 - (head_t1 - t1);
3138 }
3139 head_tmp1[0] = head_t1;
3140 tail_tmp1[0] = tail_t1;
3141 /* imaginary part */
3142 {
3143 /* Compute double-double = double-double * double. */
3144 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3145
3146 con = head_a1 * split;
3147 a11 = con - head_a1;
3148 a11 = con - a11;
3149 a21 = head_a1 - a11;
3150 con = cd[0] * split;
3151 b1 = con - cd[0];
3152 b1 = con - b1;
3153 b2 = cd[0] - b1;
3154
3155 c11 = head_a1 * cd[0];
3156 c21 =
3157 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3158
3159 c2 = tail_a1 * cd[0];
3160 t1 = c11 + c2;
3161 t2 = (c2 - (t1 - c11)) + c21;
3162
3163 head_t1 = t1 + t2;
3164 tail_t1 = t2 - (head_t1 - t1);
3165 }
3166 {
3167 /* Compute double-double = double-double * double. */
3168 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3169
3170 con = head_a0 * split;
3171 a11 = con - head_a0;
3172 a11 = con - a11;
3173 a21 = head_a0 - a11;
3174 con = cd[1] * split;
3175 b1 = con - cd[1];
3176 b1 = con - b1;
3177 b2 = cd[1] - b1;
3178
3179 c11 = head_a0 * cd[1];
3180 c21 =
3181 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3182
3183 c2 = tail_a0 * cd[1];
3184 t1 = c11 + c2;
3185 t2 = (c2 - (t1 - c11)) + c21;
3186
3187 head_t2 = t1 + t2;
3188 tail_t2 = t2 - (head_t2 - t1);
3189 }
3190 {
3191 /* Compute double-double = double-double + double-double. */
3192 double bv;
3193 double s1, s2, t1, t2;
3194
3195 /* Add two hi words. */
3196 s1 = head_t1 + head_t2;
3197 bv = s1 - head_t1;
3198 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
3199
3200 /* Add two lo words. */
3201 t1 = tail_t1 + tail_t2;
3202 bv = t1 - tail_t1;
3203 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
3204
3205 s2 += t1;
3206
3207 /* Renormalize (s1, s2) to (t1, s2) */
3208 t1 = s1 + s2;
3209 s2 = s2 - (t1 - s1);
3210
3211 t2 += s2;
3212
3213 /* Renormalize (t1, t2) */
3214 head_t1 = t1 + t2;
3215 tail_t1 = t2 - (head_t1 - t1);
3216 }
3217 head_tmp1[1] = head_t1;
3218 tail_tmp1[1] = tail_t1;
3219 }
3220
3221 }
3222 {
3223 double cd[2];
3224 cd[0] = (double) alpha_i[0];
3225 cd[1] = (double) alpha_i[1];
3226 {
3227 /* Compute complex-extra = complex-extra * complex-double. */
3228 double head_a0, tail_a0;
3229 double head_a1, tail_a1;
3230 double head_t1, tail_t1;
3231 double head_t2, tail_t2;
3232 head_a0 = head_sum2[0];
3233 tail_a0 = tail_sum2[0];
3234 head_a1 = head_sum2[1];
3235 tail_a1 = tail_sum2[1];
3236 /* real part */
3237 {
3238 /* Compute double-double = double-double * double. */
3239 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3240
3241 con = head_a0 * split;
3242 a11 = con - head_a0;
3243 a11 = con - a11;
3244 a21 = head_a0 - a11;
3245 con = cd[0] * split;
3246 b1 = con - cd[0];
3247 b1 = con - b1;
3248 b2 = cd[0] - b1;
3249
3250 c11 = head_a0 * cd[0];
3251 c21 =
3252 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3253
3254 c2 = tail_a0 * cd[0];
3255 t1 = c11 + c2;
3256 t2 = (c2 - (t1 - c11)) + c21;
3257
3258 head_t1 = t1 + t2;
3259 tail_t1 = t2 - (head_t1 - t1);
3260 }
3261 {
3262 /* Compute double-double = double-double * double. */
3263 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3264
3265 con = head_a1 * split;
3266 a11 = con - head_a1;
3267 a11 = con - a11;
3268 a21 = head_a1 - a11;
3269 con = cd[1] * split;
3270 b1 = con - cd[1];
3271 b1 = con - b1;
3272 b2 = cd[1] - b1;
3273
3274 c11 = head_a1 * cd[1];
3275 c21 =
3276 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3277
3278 c2 = tail_a1 * cd[1];
3279 t1 = c11 + c2;
3280 t2 = (c2 - (t1 - c11)) + c21;
3281
3282 head_t2 = t1 + t2;
3283 tail_t2 = t2 - (head_t2 - t1);
3284 }
3285 head_t2 = -head_t2;
3286 tail_t2 = -tail_t2;
3287 {
3288 /* Compute double-double = double-double + double-double. */
3289 double bv;
3290 double s1, s2, t1, t2;
3291
3292 /* Add two hi words. */
3293 s1 = head_t1 + head_t2;
3294 bv = s1 - head_t1;
3295 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
3296
3297 /* Add two lo words. */
3298 t1 = tail_t1 + tail_t2;
3299 bv = t1 - tail_t1;
3300 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
3301
3302 s2 += t1;
3303
3304 /* Renormalize (s1, s2) to (t1, s2) */
3305 t1 = s1 + s2;
3306 s2 = s2 - (t1 - s1);
3307
3308 t2 += s2;
3309
3310 /* Renormalize (t1, t2) */
3311 head_t1 = t1 + t2;
3312 tail_t1 = t2 - (head_t1 - t1);
3313 }
3314 head_tmp2[0] = head_t1;
3315 tail_tmp2[0] = tail_t1;
3316 /* imaginary part */
3317 {
3318 /* Compute double-double = double-double * double. */
3319 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3320
3321 con = head_a1 * split;
3322 a11 = con - head_a1;
3323 a11 = con - a11;
3324 a21 = head_a1 - a11;
3325 con = cd[0] * split;
3326 b1 = con - cd[0];
3327 b1 = con - b1;
3328 b2 = cd[0] - b1;
3329
3330 c11 = head_a1 * cd[0];
3331 c21 =
3332 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3333
3334 c2 = tail_a1 * cd[0];
3335 t1 = c11 + c2;
3336 t2 = (c2 - (t1 - c11)) + c21;
3337
3338 head_t1 = t1 + t2;
3339 tail_t1 = t2 - (head_t1 - t1);
3340 }
3341 {
3342 /* Compute double-double = double-double * double. */
3343 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3344
3345 con = head_a0 * split;
3346 a11 = con - head_a0;
3347 a11 = con - a11;
3348 a21 = head_a0 - a11;
3349 con = cd[1] * split;
3350 b1 = con - cd[1];
3351 b1 = con - b1;
3352 b2 = cd[1] - b1;
3353
3354 c11 = head_a0 * cd[1];
3355 c21 =
3356 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3357
3358 c2 = tail_a0 * cd[1];
3359 t1 = c11 + c2;
3360 t2 = (c2 - (t1 - c11)) + c21;
3361
3362 head_t2 = t1 + t2;
3363 tail_t2 = t2 - (head_t2 - t1);
3364 }
3365 {
3366 /* Compute double-double = double-double + double-double. */
3367 double bv;
3368 double s1, s2, t1, t2;
3369
3370 /* Add two hi words. */
3371 s1 = head_t1 + head_t2;
3372 bv = s1 - head_t1;
3373 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
3374
3375 /* Add two lo words. */
3376 t1 = tail_t1 + tail_t2;
3377 bv = t1 - tail_t1;
3378 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
3379
3380 s2 += t1;
3381
3382 /* Renormalize (s1, s2) to (t1, s2) */
3383 t1 = s1 + s2;
3384 s2 = s2 - (t1 - s1);
3385
3386 t2 += s2;
3387
3388 /* Renormalize (t1, t2) */
3389 head_t1 = t1 + t2;
3390 tail_t1 = t2 - (head_t1 - t1);
3391 }
3392 head_tmp2[1] = head_t1;
3393 tail_tmp2[1] = tail_t1;
3394 }
3395
3396 }
3397 {
3398 double head_t, tail_t;
3399 double head_a, tail_a;
3400 double head_b, tail_b;
3401 /* Real part */
3402 head_a = head_tmp1[0];
3403 tail_a = tail_tmp1[0];
3404 head_b = head_tmp2[0];
3405 tail_b = tail_tmp2[0];
3406 {
3407 /* Compute double-double = double-double + double-double. */
3408 double bv;
3409 double s1, s2, t1, t2;
3410
3411 /* Add two hi words. */
3412 s1 = head_a + head_b;
3413 bv = s1 - head_a;
3414 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3415
3416 /* Add two lo words. */
3417 t1 = tail_a + tail_b;
3418 bv = t1 - tail_a;
3419 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3420
3421 s2 += t1;
3422
3423 /* Renormalize (s1, s2) to (t1, s2) */
3424 t1 = s1 + s2;
3425 s2 = s2 - (t1 - s1);
3426
3427 t2 += s2;
3428
3429 /* Renormalize (t1, t2) */
3430 head_t = t1 + t2;
3431 tail_t = t2 - (head_t - t1);
3432 }
3433 head_tmp1[0] = head_t;
3434 tail_tmp1[0] = tail_t;
3435 /* Imaginary part */
3436 head_a = head_tmp1[1];
3437 tail_a = tail_tmp1[1];
3438 head_b = head_tmp2[1];
3439 tail_b = tail_tmp2[1];
3440 {
3441 /* Compute double-double = double-double + double-double. */
3442 double bv;
3443 double s1, s2, t1, t2;
3444
3445 /* Add two hi words. */
3446 s1 = head_a + head_b;
3447 bv = s1 - head_a;
3448 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3449
3450 /* Add two lo words. */
3451 t1 = tail_a + tail_b;
3452 bv = t1 - tail_a;
3453 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3454
3455 s2 += t1;
3456
3457 /* Renormalize (s1, s2) to (t1, s2) */
3458 t1 = s1 + s2;
3459 s2 = s2 - (t1 - s1);
3460
3461 t2 += s2;
3462
3463 /* Renormalize (t1, t2) */
3464 head_t = t1 + t2;
3465 tail_t = t2 - (head_t - t1);
3466 }
3467 head_tmp1[1] = head_t;
3468 tail_tmp1[1] = tail_t;
3469 }
3470 y_elem[0] = y_i[iy];
3471 y_elem[1] = y_i[iy + 1];
3472 {
3473 double head_e1, tail_e1;
3474 double d1;
3475 double d2;
3476 /* Real part */
3477 d1 = (double) y_elem[0] * beta_i[0];
3478 d2 = (double) -y_elem[1] * beta_i[1];
3479 {
3480 /* Compute double-double = double + double. */
3481 double e, t1, t2;
3482
3483 /* Knuth trick. */
3484 t1 = d1 + d2;
3485 e = t1 - d1;
3486 t2 = ((d2 - e) + (d1 - (t1 - e)));
3487
3488 /* The result is t1 + t2, after normalization. */
3489 head_e1 = t1 + t2;
3490 tail_e1 = t2 - (head_e1 - t1);
3491 }
3492 head_tmp2[0] = head_e1;
3493 tail_tmp2[0] = tail_e1;
3494 /* imaginary part */
3495 d1 = (double) y_elem[0] * beta_i[1];
3496 d2 = (double) y_elem[1] * beta_i[0];
3497 {
3498 /* Compute double-double = double + double. */
3499 double e, t1, t2;
3500
3501 /* Knuth trick. */
3502 t1 = d1 + d2;
3503 e = t1 - d1;
3504 t2 = ((d2 - e) + (d1 - (t1 - e)));
3505
3506 /* The result is t1 + t2, after normalization. */
3507 head_e1 = t1 + t2;
3508 tail_e1 = t2 - (head_e1 - t1);
3509 }
3510 head_tmp2[1] = head_e1;
3511 tail_tmp2[1] = tail_e1;
3512 }
3513 {
3514 double head_t, tail_t;
3515 double head_a, tail_a;
3516 double head_b, tail_b;
3517 /* Real part */
3518 head_a = head_tmp1[0];
3519 tail_a = tail_tmp1[0];
3520 head_b = head_tmp2[0];
3521 tail_b = tail_tmp2[0];
3522 {
3523 /* Compute double-double = double-double + double-double. */
3524 double bv;
3525 double s1, s2, t1, t2;
3526
3527 /* Add two hi words. */
3528 s1 = head_a + head_b;
3529 bv = s1 - head_a;
3530 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3531
3532 /* Add two lo words. */
3533 t1 = tail_a + tail_b;
3534 bv = t1 - tail_a;
3535 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3536
3537 s2 += t1;
3538
3539 /* Renormalize (s1, s2) to (t1, s2) */
3540 t1 = s1 + s2;
3541 s2 = s2 - (t1 - s1);
3542
3543 t2 += s2;
3544
3545 /* Renormalize (t1, t2) */
3546 head_t = t1 + t2;
3547 tail_t = t2 - (head_t - t1);
3548 }
3549 head_tmp1[0] = head_t;
3550 tail_tmp1[0] = tail_t;
3551 /* Imaginary part */
3552 head_a = head_tmp1[1];
3553 tail_a = tail_tmp1[1];
3554 head_b = head_tmp2[1];
3555 tail_b = tail_tmp2[1];
3556 {
3557 /* Compute double-double = double-double + double-double. */
3558 double bv;
3559 double s1, s2, t1, t2;
3560
3561 /* Add two hi words. */
3562 s1 = head_a + head_b;
3563 bv = s1 - head_a;
3564 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3565
3566 /* Add two lo words. */
3567 t1 = tail_a + tail_b;
3568 bv = t1 - tail_a;
3569 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3570
3571 s2 += t1;
3572
3573 /* Renormalize (s1, s2) to (t1, s2) */
3574 t1 = s1 + s2;
3575 s2 = s2 - (t1 - s1);
3576
3577 t2 += s2;
3578
3579 /* Renormalize (t1, t2) */
3580 head_t = t1 + t2;
3581 tail_t = t2 - (head_t - t1);
3582 }
3583 head_tmp1[1] = head_t;
3584 tail_tmp1[1] = tail_t;
3585 }
3586 y_i[iy] = head_tmp1[0];
3587 y_i[iy + 1] = head_tmp1[1];
3588 ai += incai;
3589 iy += incy;
3590 }
3591 }
3592 }
3593
3594 } else {
3595
3596 /* if beta = 0, we can save m multiplies:
3597 y = alpha*A*head_x + alpha*A*tail_x */
3598 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
3599 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
3600 /* save m more multiplies if alpha = 1 */
3601 ai = 0;
3602 iy = ky;
3603 for (i = 0; i < leny; i++) {
3604 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
3605 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
3606 0.0;
3607 aij = ai;
3608 jx = kx;
3609 for (j = 0; j < lenx; j++) {
3610 a_elem[0] = a_i[aij];
3611 a_elem[1] = a_i[aij + 1];
3612
3613 x_elem[0] = head_x_i[jx];
3614 x_elem[1] = head_x_i[jx + 1];
3615 {
3616 double head_e1, tail_e1;
3617 double d1;
3618 double d2;
3619 /* Real part */
3620 d1 = (double) a_elem[0] * x_elem[0];
3621 d2 = (double) -a_elem[1] * x_elem[1];
3622 {
3623 /* Compute double-double = double + double. */
3624 double e, t1, t2;
3625
3626 /* Knuth trick. */
3627 t1 = d1 + d2;
3628 e = t1 - d1;
3629 t2 = ((d2 - e) + (d1 - (t1 - e)));
3630
3631 /* The result is t1 + t2, after normalization. */
3632 head_e1 = t1 + t2;
3633 tail_e1 = t2 - (head_e1 - t1);
3634 }
3635 head_prod[0] = head_e1;
3636 tail_prod[0] = tail_e1;
3637 /* imaginary part */
3638 d1 = (double) a_elem[0] * x_elem[1];
3639 d2 = (double) a_elem[1] * x_elem[0];
3640 {
3641 /* Compute double-double = double + double. */
3642 double e, t1, t2;
3643
3644 /* Knuth trick. */
3645 t1 = d1 + d2;
3646 e = t1 - d1;
3647 t2 = ((d2 - e) + (d1 - (t1 - e)));
3648
3649 /* The result is t1 + t2, after normalization. */
3650 head_e1 = t1 + t2;
3651 tail_e1 = t2 - (head_e1 - t1);
3652 }
3653 head_prod[1] = head_e1;
3654 tail_prod[1] = tail_e1;
3655 }
3656 {
3657 double head_t, tail_t;
3658 double head_a, tail_a;
3659 double head_b, tail_b;
3660 /* Real part */
3661 head_a = head_sum[0];
3662 tail_a = tail_sum[0];
3663 head_b = head_prod[0];
3664 tail_b = tail_prod[0];
3665 {
3666 /* Compute double-double = double-double + double-double. */
3667 double bv;
3668 double s1, s2, t1, t2;
3669
3670 /* Add two hi words. */
3671 s1 = head_a + head_b;
3672 bv = s1 - head_a;
3673 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3674
3675 /* Add two lo words. */
3676 t1 = tail_a + tail_b;
3677 bv = t1 - tail_a;
3678 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3679
3680 s2 += t1;
3681
3682 /* Renormalize (s1, s2) to (t1, s2) */
3683 t1 = s1 + s2;
3684 s2 = s2 - (t1 - s1);
3685
3686 t2 += s2;
3687
3688 /* Renormalize (t1, t2) */
3689 head_t = t1 + t2;
3690 tail_t = t2 - (head_t - t1);
3691 }
3692 head_sum[0] = head_t;
3693 tail_sum[0] = tail_t;
3694 /* Imaginary part */
3695 head_a = head_sum[1];
3696 tail_a = tail_sum[1];
3697 head_b = head_prod[1];
3698 tail_b = tail_prod[1];
3699 {
3700 /* Compute double-double = double-double + double-double. */
3701 double bv;
3702 double s1, s2, t1, t2;
3703
3704 /* Add two hi words. */
3705 s1 = head_a + head_b;
3706 bv = s1 - head_a;
3707 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3708
3709 /* Add two lo words. */
3710 t1 = tail_a + tail_b;
3711 bv = t1 - tail_a;
3712 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3713
3714 s2 += t1;
3715
3716 /* Renormalize (s1, s2) to (t1, s2) */
3717 t1 = s1 + s2;
3718 s2 = s2 - (t1 - s1);
3719
3720 t2 += s2;
3721
3722 /* Renormalize (t1, t2) */
3723 head_t = t1 + t2;
3724 tail_t = t2 - (head_t - t1);
3725 }
3726 head_sum[1] = head_t;
3727 tail_sum[1] = tail_t;
3728 }
3729 x_elem[0] = tail_x_i[jx];
3730 x_elem[1] = tail_x_i[jx + 1];
3731 {
3732 double head_e1, tail_e1;
3733 double d1;
3734 double d2;
3735 /* Real part */
3736 d1 = (double) a_elem[0] * x_elem[0];
3737 d2 = (double) -a_elem[1] * x_elem[1];
3738 {
3739 /* Compute double-double = double + double. */
3740 double e, t1, t2;
3741
3742 /* Knuth trick. */
3743 t1 = d1 + d2;
3744 e = t1 - d1;
3745 t2 = ((d2 - e) + (d1 - (t1 - e)));
3746
3747 /* The result is t1 + t2, after normalization. */
3748 head_e1 = t1 + t2;
3749 tail_e1 = t2 - (head_e1 - t1);
3750 }
3751 head_prod[0] = head_e1;
3752 tail_prod[0] = tail_e1;
3753 /* imaginary part */
3754 d1 = (double) a_elem[0] * x_elem[1];
3755 d2 = (double) a_elem[1] * x_elem[0];
3756 {
3757 /* Compute double-double = double + double. */
3758 double e, t1, t2;
3759
3760 /* Knuth trick. */
3761 t1 = d1 + d2;
3762 e = t1 - d1;
3763 t2 = ((d2 - e) + (d1 - (t1 - e)));
3764
3765 /* The result is t1 + t2, after normalization. */
3766 head_e1 = t1 + t2;
3767 tail_e1 = t2 - (head_e1 - t1);
3768 }
3769 head_prod[1] = head_e1;
3770 tail_prod[1] = tail_e1;
3771 }
3772 {
3773 double head_t, tail_t;
3774 double head_a, tail_a;
3775 double head_b, tail_b;
3776 /* Real part */
3777 head_a = head_sum2[0];
3778 tail_a = tail_sum2[0];
3779 head_b = head_prod[0];
3780 tail_b = tail_prod[0];
3781 {
3782 /* Compute double-double = double-double + double-double. */
3783 double bv;
3784 double s1, s2, t1, t2;
3785
3786 /* Add two hi words. */
3787 s1 = head_a + head_b;
3788 bv = s1 - head_a;
3789 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3790
3791 /* Add two lo words. */
3792 t1 = tail_a + tail_b;
3793 bv = t1 - tail_a;
3794 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3795
3796 s2 += t1;
3797
3798 /* Renormalize (s1, s2) to (t1, s2) */
3799 t1 = s1 + s2;
3800 s2 = s2 - (t1 - s1);
3801
3802 t2 += s2;
3803
3804 /* Renormalize (t1, t2) */
3805 head_t = t1 + t2;
3806 tail_t = t2 - (head_t - t1);
3807 }
3808 head_sum2[0] = head_t;
3809 tail_sum2[0] = tail_t;
3810 /* Imaginary part */
3811 head_a = head_sum2[1];
3812 tail_a = tail_sum2[1];
3813 head_b = head_prod[1];
3814 tail_b = tail_prod[1];
3815 {
3816 /* Compute double-double = double-double + double-double. */
3817 double bv;
3818 double s1, s2, t1, t2;
3819
3820 /* Add two hi words. */
3821 s1 = head_a + head_b;
3822 bv = s1 - head_a;
3823 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3824
3825 /* Add two lo words. */
3826 t1 = tail_a + tail_b;
3827 bv = t1 - tail_a;
3828 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3829
3830 s2 += t1;
3831
3832 /* Renormalize (s1, s2) to (t1, s2) */
3833 t1 = s1 + s2;
3834 s2 = s2 - (t1 - s1);
3835
3836 t2 += s2;
3837
3838 /* Renormalize (t1, t2) */
3839 head_t = t1 + t2;
3840 tail_t = t2 - (head_t - t1);
3841 }
3842 head_sum2[1] = head_t;
3843 tail_sum2[1] = tail_t;
3844 }
3845 aij += incaij;
3846 jx += incx;
3847 }
3848 {
3849 double head_t, tail_t;
3850 double head_a, tail_a;
3851 double head_b, tail_b;
3852 /* Real part */
3853 head_a = head_sum[0];
3854 tail_a = tail_sum[0];
3855 head_b = head_sum2[0];
3856 tail_b = tail_sum2[0];
3857 {
3858 /* Compute double-double = double-double + double-double. */
3859 double bv;
3860 double s1, s2, t1, t2;
3861
3862 /* Add two hi words. */
3863 s1 = head_a + head_b;
3864 bv = s1 - head_a;
3865 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3866
3867 /* Add two lo words. */
3868 t1 = tail_a + tail_b;
3869 bv = t1 - tail_a;
3870 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3871
3872 s2 += t1;
3873
3874 /* Renormalize (s1, s2) to (t1, s2) */
3875 t1 = s1 + s2;
3876 s2 = s2 - (t1 - s1);
3877
3878 t2 += s2;
3879
3880 /* Renormalize (t1, t2) */
3881 head_t = t1 + t2;
3882 tail_t = t2 - (head_t - t1);
3883 }
3884 head_sum[0] = head_t;
3885 tail_sum[0] = tail_t;
3886 /* Imaginary part */
3887 head_a = head_sum[1];
3888 tail_a = tail_sum[1];
3889 head_b = head_sum2[1];
3890 tail_b = tail_sum2[1];
3891 {
3892 /* Compute double-double = double-double + double-double. */
3893 double bv;
3894 double s1, s2, t1, t2;
3895
3896 /* Add two hi words. */
3897 s1 = head_a + head_b;
3898 bv = s1 - head_a;
3899 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3900
3901 /* Add two lo words. */
3902 t1 = tail_a + tail_b;
3903 bv = t1 - tail_a;
3904 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3905
3906 s2 += t1;
3907
3908 /* Renormalize (s1, s2) to (t1, s2) */
3909 t1 = s1 + s2;
3910 s2 = s2 - (t1 - s1);
3911
3912 t2 += s2;
3913
3914 /* Renormalize (t1, t2) */
3915 head_t = t1 + t2;
3916 tail_t = t2 - (head_t - t1);
3917 }
3918 head_sum[1] = head_t;
3919 tail_sum[1] = tail_t;
3920 }
3921 y_i[iy] = head_sum[0];
3922 y_i[iy + 1] = head_sum[1];
3923 ai += incai;
3924 iy += incy;
3925 } /* end for */
3926 } else { /* alpha != 1 */
3927 ai = 0;
3928 iy = ky;
3929 for (i = 0; i < leny; i++) {
3930 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
3931 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
3932 0.0;
3933 aij = ai;
3934 jx = kx;
3935 for (j = 0; j < lenx; j++) {
3936 a_elem[0] = a_i[aij];
3937 a_elem[1] = a_i[aij + 1];
3938
3939 x_elem[0] = head_x_i[jx];
3940 x_elem[1] = head_x_i[jx + 1];
3941 {
3942 double head_e1, tail_e1;
3943 double d1;
3944 double d2;
3945 /* Real part */
3946 d1 = (double) a_elem[0] * x_elem[0];
3947 d2 = (double) -a_elem[1] * x_elem[1];
3948 {
3949 /* Compute double-double = double + double. */
3950 double e, t1, t2;
3951
3952 /* Knuth trick. */
3953 t1 = d1 + d2;
3954 e = t1 - d1;
3955 t2 = ((d2 - e) + (d1 - (t1 - e)));
3956
3957 /* The result is t1 + t2, after normalization. */
3958 head_e1 = t1 + t2;
3959 tail_e1 = t2 - (head_e1 - t1);
3960 }
3961 head_prod[0] = head_e1;
3962 tail_prod[0] = tail_e1;
3963 /* imaginary part */
3964 d1 = (double) a_elem[0] * x_elem[1];
3965 d2 = (double) a_elem[1] * x_elem[0];
3966 {
3967 /* Compute double-double = double + double. */
3968 double e, t1, t2;
3969
3970 /* Knuth trick. */
3971 t1 = d1 + d2;
3972 e = t1 - d1;
3973 t2 = ((d2 - e) + (d1 - (t1 - e)));
3974
3975 /* The result is t1 + t2, after normalization. */
3976 head_e1 = t1 + t2;
3977 tail_e1 = t2 - (head_e1 - t1);
3978 }
3979 head_prod[1] = head_e1;
3980 tail_prod[1] = tail_e1;
3981 }
3982 {
3983 double head_t, tail_t;
3984 double head_a, tail_a;
3985 double head_b, tail_b;
3986 /* Real part */
3987 head_a = head_sum[0];
3988 tail_a = tail_sum[0];
3989 head_b = head_prod[0];
3990 tail_b = tail_prod[0];
3991 {
3992 /* Compute double-double = double-double + double-double. */
3993 double bv;
3994 double s1, s2, t1, t2;
3995
3996 /* Add two hi words. */
3997 s1 = head_a + head_b;
3998 bv = s1 - head_a;
3999 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4000
4001 /* Add two lo words. */
4002 t1 = tail_a + tail_b;
4003 bv = t1 - tail_a;
4004 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4005
4006 s2 += t1;
4007
4008 /* Renormalize (s1, s2) to (t1, s2) */
4009 t1 = s1 + s2;
4010 s2 = s2 - (t1 - s1);
4011
4012 t2 += s2;
4013
4014 /* Renormalize (t1, t2) */
4015 head_t = t1 + t2;
4016 tail_t = t2 - (head_t - t1);
4017 }
4018 head_sum[0] = head_t;
4019 tail_sum[0] = tail_t;
4020 /* Imaginary part */
4021 head_a = head_sum[1];
4022 tail_a = tail_sum[1];
4023 head_b = head_prod[1];
4024 tail_b = tail_prod[1];
4025 {
4026 /* Compute double-double = double-double + double-double. */
4027 double bv;
4028 double s1, s2, t1, t2;
4029
4030 /* Add two hi words. */
4031 s1 = head_a + head_b;
4032 bv = s1 - head_a;
4033 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4034
4035 /* Add two lo words. */
4036 t1 = tail_a + tail_b;
4037 bv = t1 - tail_a;
4038 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4039
4040 s2 += t1;
4041
4042 /* Renormalize (s1, s2) to (t1, s2) */
4043 t1 = s1 + s2;
4044 s2 = s2 - (t1 - s1);
4045
4046 t2 += s2;
4047
4048 /* Renormalize (t1, t2) */
4049 head_t = t1 + t2;
4050 tail_t = t2 - (head_t - t1);
4051 }
4052 head_sum[1] = head_t;
4053 tail_sum[1] = tail_t;
4054 }
4055 x_elem[0] = tail_x_i[jx];
4056 x_elem[1] = tail_x_i[jx + 1];
4057 {
4058 double head_e1, tail_e1;
4059 double d1;
4060 double d2;
4061 /* Real part */
4062 d1 = (double) a_elem[0] * x_elem[0];
4063 d2 = (double) -a_elem[1] * x_elem[1];
4064 {
4065 /* Compute double-double = double + double. */
4066 double e, t1, t2;
4067
4068 /* Knuth trick. */
4069 t1 = d1 + d2;
4070 e = t1 - d1;
4071 t2 = ((d2 - e) + (d1 - (t1 - e)));
4072
4073 /* The result is t1 + t2, after normalization. */
4074 head_e1 = t1 + t2;
4075 tail_e1 = t2 - (head_e1 - t1);
4076 }
4077 head_prod[0] = head_e1;
4078 tail_prod[0] = tail_e1;
4079 /* imaginary part */
4080 d1 = (double) a_elem[0] * x_elem[1];
4081 d2 = (double) a_elem[1] * x_elem[0];
4082 {
4083 /* Compute double-double = double + double. */
4084 double e, t1, t2;
4085
4086 /* Knuth trick. */
4087 t1 = d1 + d2;
4088 e = t1 - d1;
4089 t2 = ((d2 - e) + (d1 - (t1 - e)));
4090
4091 /* The result is t1 + t2, after normalization. */
4092 head_e1 = t1 + t2;
4093 tail_e1 = t2 - (head_e1 - t1);
4094 }
4095 head_prod[1] = head_e1;
4096 tail_prod[1] = tail_e1;
4097 }
4098 {
4099 double head_t, tail_t;
4100 double head_a, tail_a;
4101 double head_b, tail_b;
4102 /* Real part */
4103 head_a = head_sum2[0];
4104 tail_a = tail_sum2[0];
4105 head_b = head_prod[0];
4106 tail_b = tail_prod[0];
4107 {
4108 /* Compute double-double = double-double + double-double. */
4109 double bv;
4110 double s1, s2, t1, t2;
4111
4112 /* Add two hi words. */
4113 s1 = head_a + head_b;
4114 bv = s1 - head_a;
4115 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4116
4117 /* Add two lo words. */
4118 t1 = tail_a + tail_b;
4119 bv = t1 - tail_a;
4120 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4121
4122 s2 += t1;
4123
4124 /* Renormalize (s1, s2) to (t1, s2) */
4125 t1 = s1 + s2;
4126 s2 = s2 - (t1 - s1);
4127
4128 t2 += s2;
4129
4130 /* Renormalize (t1, t2) */
4131 head_t = t1 + t2;
4132 tail_t = t2 - (head_t - t1);
4133 }
4134 head_sum2[0] = head_t;
4135 tail_sum2[0] = tail_t;
4136 /* Imaginary part */
4137 head_a = head_sum2[1];
4138 tail_a = tail_sum2[1];
4139 head_b = head_prod[1];
4140 tail_b = tail_prod[1];
4141 {
4142 /* Compute double-double = double-double + double-double. */
4143 double bv;
4144 double s1, s2, t1, t2;
4145
4146 /* Add two hi words. */
4147 s1 = head_a + head_b;
4148 bv = s1 - head_a;
4149 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4150
4151 /* Add two lo words. */
4152 t1 = tail_a + tail_b;
4153 bv = t1 - tail_a;
4154 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4155
4156 s2 += t1;
4157
4158 /* Renormalize (s1, s2) to (t1, s2) */
4159 t1 = s1 + s2;
4160 s2 = s2 - (t1 - s1);
4161
4162 t2 += s2;
4163
4164 /* Renormalize (t1, t2) */
4165 head_t = t1 + t2;
4166 tail_t = t2 - (head_t - t1);
4167 }
4168 head_sum2[1] = head_t;
4169 tail_sum2[1] = tail_t;
4170 }
4171 aij += incaij;
4172 jx += incx;
4173 }
4174 {
4175 double cd[2];
4176 cd[0] = (double) alpha_i[0];
4177 cd[1] = (double) alpha_i[1];
4178 {
4179 /* Compute complex-extra = complex-extra * complex-double. */
4180 double head_a0, tail_a0;
4181 double head_a1, tail_a1;
4182 double head_t1, tail_t1;
4183 double head_t2, tail_t2;
4184 head_a0 = head_sum[0];
4185 tail_a0 = tail_sum[0];
4186 head_a1 = head_sum[1];
4187 tail_a1 = tail_sum[1];
4188 /* real part */
4189 {
4190 /* Compute double-double = double-double * double. */
4191 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4192
4193 con = head_a0 * split;
4194 a11 = con - head_a0;
4195 a11 = con - a11;
4196 a21 = head_a0 - a11;
4197 con = cd[0] * split;
4198 b1 = con - cd[0];
4199 b1 = con - b1;
4200 b2 = cd[0] - b1;
4201
4202 c11 = head_a0 * cd[0];
4203 c21 =
4204 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4205
4206 c2 = tail_a0 * cd[0];
4207 t1 = c11 + c2;
4208 t2 = (c2 - (t1 - c11)) + c21;
4209
4210 head_t1 = t1 + t2;
4211 tail_t1 = t2 - (head_t1 - t1);
4212 }
4213 {
4214 /* Compute double-double = double-double * double. */
4215 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4216
4217 con = head_a1 * split;
4218 a11 = con - head_a1;
4219 a11 = con - a11;
4220 a21 = head_a1 - a11;
4221 con = cd[1] * split;
4222 b1 = con - cd[1];
4223 b1 = con - b1;
4224 b2 = cd[1] - b1;
4225
4226 c11 = head_a1 * cd[1];
4227 c21 =
4228 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4229
4230 c2 = tail_a1 * cd[1];
4231 t1 = c11 + c2;
4232 t2 = (c2 - (t1 - c11)) + c21;
4233
4234 head_t2 = t1 + t2;
4235 tail_t2 = t2 - (head_t2 - t1);
4236 }
4237 head_t2 = -head_t2;
4238 tail_t2 = -tail_t2;
4239 {
4240 /* Compute double-double = double-double + double-double. */
4241 double bv;
4242 double s1, s2, t1, t2;
4243
4244 /* Add two hi words. */
4245 s1 = head_t1 + head_t2;
4246 bv = s1 - head_t1;
4247 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
4248
4249 /* Add two lo words. */
4250 t1 = tail_t1 + tail_t2;
4251 bv = t1 - tail_t1;
4252 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
4253
4254 s2 += t1;
4255
4256 /* Renormalize (s1, s2) to (t1, s2) */
4257 t1 = s1 + s2;
4258 s2 = s2 - (t1 - s1);
4259
4260 t2 += s2;
4261
4262 /* Renormalize (t1, t2) */
4263 head_t1 = t1 + t2;
4264 tail_t1 = t2 - (head_t1 - t1);
4265 }
4266 head_tmp1[0] = head_t1;
4267 tail_tmp1[0] = tail_t1;
4268 /* imaginary part */
4269 {
4270 /* Compute double-double = double-double * double. */
4271 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4272
4273 con = head_a1 * split;
4274 a11 = con - head_a1;
4275 a11 = con - a11;
4276 a21 = head_a1 - a11;
4277 con = cd[0] * split;
4278 b1 = con - cd[0];
4279 b1 = con - b1;
4280 b2 = cd[0] - b1;
4281
4282 c11 = head_a1 * cd[0];
4283 c21 =
4284 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4285
4286 c2 = tail_a1 * cd[0];
4287 t1 = c11 + c2;
4288 t2 = (c2 - (t1 - c11)) + c21;
4289
4290 head_t1 = t1 + t2;
4291 tail_t1 = t2 - (head_t1 - t1);
4292 }
4293 {
4294 /* Compute double-double = double-double * double. */
4295 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4296
4297 con = head_a0 * split;
4298 a11 = con - head_a0;
4299 a11 = con - a11;
4300 a21 = head_a0 - a11;
4301 con = cd[1] * split;
4302 b1 = con - cd[1];
4303 b1 = con - b1;
4304 b2 = cd[1] - b1;
4305
4306 c11 = head_a0 * cd[1];
4307 c21 =
4308 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4309
4310 c2 = tail_a0 * cd[1];
4311 t1 = c11 + c2;
4312 t2 = (c2 - (t1 - c11)) + c21;
4313
4314 head_t2 = t1 + t2;
4315 tail_t2 = t2 - (head_t2 - t1);
4316 }
4317 {
4318 /* Compute double-double = double-double + double-double. */
4319 double bv;
4320 double s1, s2, t1, t2;
4321
4322 /* Add two hi words. */
4323 s1 = head_t1 + head_t2;
4324 bv = s1 - head_t1;
4325 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
4326
4327 /* Add two lo words. */
4328 t1 = tail_t1 + tail_t2;
4329 bv = t1 - tail_t1;
4330 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
4331
4332 s2 += t1;
4333
4334 /* Renormalize (s1, s2) to (t1, s2) */
4335 t1 = s1 + s2;
4336 s2 = s2 - (t1 - s1);
4337
4338 t2 += s2;
4339
4340 /* Renormalize (t1, t2) */
4341 head_t1 = t1 + t2;
4342 tail_t1 = t2 - (head_t1 - t1);
4343 }
4344 head_tmp1[1] = head_t1;
4345 tail_tmp1[1] = tail_t1;
4346 }
4347
4348 }
4349 {
4350 double cd[2];
4351 cd[0] = (double) alpha_i[0];
4352 cd[1] = (double) alpha_i[1];
4353 {
4354 /* Compute complex-extra = complex-extra * complex-double. */
4355 double head_a0, tail_a0;
4356 double head_a1, tail_a1;
4357 double head_t1, tail_t1;
4358 double head_t2, tail_t2;
4359 head_a0 = head_sum2[0];
4360 tail_a0 = tail_sum2[0];
4361 head_a1 = head_sum2[1];
4362 tail_a1 = tail_sum2[1];
4363 /* real part */
4364 {
4365 /* Compute double-double = double-double * double. */
4366 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4367
4368 con = head_a0 * split;
4369 a11 = con - head_a0;
4370 a11 = con - a11;
4371 a21 = head_a0 - a11;
4372 con = cd[0] * split;
4373 b1 = con - cd[0];
4374 b1 = con - b1;
4375 b2 = cd[0] - b1;
4376
4377 c11 = head_a0 * cd[0];
4378 c21 =
4379 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4380
4381 c2 = tail_a0 * cd[0];
4382 t1 = c11 + c2;
4383 t2 = (c2 - (t1 - c11)) + c21;
4384
4385 head_t1 = t1 + t2;
4386 tail_t1 = t2 - (head_t1 - t1);
4387 }
4388 {
4389 /* Compute double-double = double-double * double. */
4390 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4391
4392 con = head_a1 * split;
4393 a11 = con - head_a1;
4394 a11 = con - a11;
4395 a21 = head_a1 - a11;
4396 con = cd[1] * split;
4397 b1 = con - cd[1];
4398 b1 = con - b1;
4399 b2 = cd[1] - b1;
4400
4401 c11 = head_a1 * cd[1];
4402 c21 =
4403 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4404
4405 c2 = tail_a1 * cd[1];
4406 t1 = c11 + c2;
4407 t2 = (c2 - (t1 - c11)) + c21;
4408
4409 head_t2 = t1 + t2;
4410 tail_t2 = t2 - (head_t2 - t1);
4411 }
4412 head_t2 = -head_t2;
4413 tail_t2 = -tail_t2;
4414 {
4415 /* Compute double-double = double-double + double-double. */
4416 double bv;
4417 double s1, s2, t1, t2;
4418
4419 /* Add two hi words. */
4420 s1 = head_t1 + head_t2;
4421 bv = s1 - head_t1;
4422 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
4423
4424 /* Add two lo words. */
4425 t1 = tail_t1 + tail_t2;
4426 bv = t1 - tail_t1;
4427 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
4428
4429 s2 += t1;
4430
4431 /* Renormalize (s1, s2) to (t1, s2) */
4432 t1 = s1 + s2;
4433 s2 = s2 - (t1 - s1);
4434
4435 t2 += s2;
4436
4437 /* Renormalize (t1, t2) */
4438 head_t1 = t1 + t2;
4439 tail_t1 = t2 - (head_t1 - t1);
4440 }
4441 head_tmp2[0] = head_t1;
4442 tail_tmp2[0] = tail_t1;
4443 /* imaginary part */
4444 {
4445 /* Compute double-double = double-double * double. */
4446 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4447
4448 con = head_a1 * split;
4449 a11 = con - head_a1;
4450 a11 = con - a11;
4451 a21 = head_a1 - a11;
4452 con = cd[0] * split;
4453 b1 = con - cd[0];
4454 b1 = con - b1;
4455 b2 = cd[0] - b1;
4456
4457 c11 = head_a1 * cd[0];
4458 c21 =
4459 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4460
4461 c2 = tail_a1 * cd[0];
4462 t1 = c11 + c2;
4463 t2 = (c2 - (t1 - c11)) + c21;
4464
4465 head_t1 = t1 + t2;
4466 tail_t1 = t2 - (head_t1 - t1);
4467 }
4468 {
4469 /* Compute double-double = double-double * double. */
4470 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4471
4472 con = head_a0 * split;
4473 a11 = con - head_a0;
4474 a11 = con - a11;
4475 a21 = head_a0 - a11;
4476 con = cd[1] * split;
4477 b1 = con - cd[1];
4478 b1 = con - b1;
4479 b2 = cd[1] - b1;
4480
4481 c11 = head_a0 * cd[1];
4482 c21 =
4483 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4484
4485 c2 = tail_a0 * cd[1];
4486 t1 = c11 + c2;
4487 t2 = (c2 - (t1 - c11)) + c21;
4488
4489 head_t2 = t1 + t2;
4490 tail_t2 = t2 - (head_t2 - t1);
4491 }
4492 {
4493 /* Compute double-double = double-double + double-double. */
4494 double bv;
4495 double s1, s2, t1, t2;
4496
4497 /* Add two hi words. */
4498 s1 = head_t1 + head_t2;
4499 bv = s1 - head_t1;
4500 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
4501
4502 /* Add two lo words. */
4503 t1 = tail_t1 + tail_t2;
4504 bv = t1 - tail_t1;
4505 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
4506
4507 s2 += t1;
4508
4509 /* Renormalize (s1, s2) to (t1, s2) */
4510 t1 = s1 + s2;
4511 s2 = s2 - (t1 - s1);
4512
4513 t2 += s2;
4514
4515 /* Renormalize (t1, t2) */
4516 head_t1 = t1 + t2;
4517 tail_t1 = t2 - (head_t1 - t1);
4518 }
4519 head_tmp2[1] = head_t1;
4520 tail_tmp2[1] = tail_t1;
4521 }
4522
4523 }
4524 {
4525 double head_t, tail_t;
4526 double head_a, tail_a;
4527 double head_b, tail_b;
4528 /* Real part */
4529 head_a = head_tmp1[0];
4530 tail_a = tail_tmp1[0];
4531 head_b = head_tmp2[0];
4532 tail_b = tail_tmp2[0];
4533 {
4534 /* Compute double-double = double-double + double-double. */
4535 double bv;
4536 double s1, s2, t1, t2;
4537
4538 /* Add two hi words. */
4539 s1 = head_a + head_b;
4540 bv = s1 - head_a;
4541 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4542
4543 /* Add two lo words. */
4544 t1 = tail_a + tail_b;
4545 bv = t1 - tail_a;
4546 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4547
4548 s2 += t1;
4549
4550 /* Renormalize (s1, s2) to (t1, s2) */
4551 t1 = s1 + s2;
4552 s2 = s2 - (t1 - s1);
4553
4554 t2 += s2;
4555
4556 /* Renormalize (t1, t2) */
4557 head_t = t1 + t2;
4558 tail_t = t2 - (head_t - t1);
4559 }
4560 head_tmp1[0] = head_t;
4561 tail_tmp1[0] = tail_t;
4562 /* Imaginary part */
4563 head_a = head_tmp1[1];
4564 tail_a = tail_tmp1[1];
4565 head_b = head_tmp2[1];
4566 tail_b = tail_tmp2[1];
4567 {
4568 /* Compute double-double = double-double + double-double. */
4569 double bv;
4570 double s1, s2, t1, t2;
4571
4572 /* Add two hi words. */
4573 s1 = head_a + head_b;
4574 bv = s1 - head_a;
4575 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4576
4577 /* Add two lo words. */
4578 t1 = tail_a + tail_b;
4579 bv = t1 - tail_a;
4580 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4581
4582 s2 += t1;
4583
4584 /* Renormalize (s1, s2) to (t1, s2) */
4585 t1 = s1 + s2;
4586 s2 = s2 - (t1 - s1);
4587
4588 t2 += s2;
4589
4590 /* Renormalize (t1, t2) */
4591 head_t = t1 + t2;
4592 tail_t = t2 - (head_t - t1);
4593 }
4594 head_tmp1[1] = head_t;
4595 tail_tmp1[1] = tail_t;
4596 }
4597 y_i[iy] = head_tmp1[0];
4598 y_i[iy + 1] = head_tmp1[1];
4599 ai += incai;
4600 iy += incy;
4601 }
4602 }
4603 } else { /* beta != 0 */
4604 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
4605 /* save m multiplies if alpha = 1 */
4606 ai = 0;
4607 iy = ky;
4608 for (i = 0; i < leny; i++) {
4609 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;;
4610 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
4611 0.0;;
4612 aij = ai;
4613 jx = kx;
4614 for (j = 0; j < lenx; j++) {
4615 a_elem[0] = a_i[aij];
4616 a_elem[1] = a_i[aij + 1];
4617
4618 x_elem[0] = head_x_i[jx];
4619 x_elem[1] = head_x_i[jx + 1];
4620 {
4621 double head_e1, tail_e1;
4622 double d1;
4623 double d2;
4624 /* Real part */
4625 d1 = (double) a_elem[0] * x_elem[0];
4626 d2 = (double) -a_elem[1] * x_elem[1];
4627 {
4628 /* Compute double-double = double + double. */
4629 double e, t1, t2;
4630
4631 /* Knuth trick. */
4632 t1 = d1 + d2;
4633 e = t1 - d1;
4634 t2 = ((d2 - e) + (d1 - (t1 - e)));
4635
4636 /* The result is t1 + t2, after normalization. */
4637 head_e1 = t1 + t2;
4638 tail_e1 = t2 - (head_e1 - t1);
4639 }
4640 head_prod[0] = head_e1;
4641 tail_prod[0] = tail_e1;
4642 /* imaginary part */
4643 d1 = (double) a_elem[0] * x_elem[1];
4644 d2 = (double) a_elem[1] * x_elem[0];
4645 {
4646 /* Compute double-double = double + double. */
4647 double e, t1, t2;
4648
4649 /* Knuth trick. */
4650 t1 = d1 + d2;
4651 e = t1 - d1;
4652 t2 = ((d2 - e) + (d1 - (t1 - e)));
4653
4654 /* The result is t1 + t2, after normalization. */
4655 head_e1 = t1 + t2;
4656 tail_e1 = t2 - (head_e1 - t1);
4657 }
4658 head_prod[1] = head_e1;
4659 tail_prod[1] = tail_e1;
4660 }
4661 {
4662 double head_t, tail_t;
4663 double head_a, tail_a;
4664 double head_b, tail_b;
4665 /* Real part */
4666 head_a = head_sum[0];
4667 tail_a = tail_sum[0];
4668 head_b = head_prod[0];
4669 tail_b = tail_prod[0];
4670 {
4671 /* Compute double-double = double-double + double-double. */
4672 double bv;
4673 double s1, s2, t1, t2;
4674
4675 /* Add two hi words. */
4676 s1 = head_a + head_b;
4677 bv = s1 - head_a;
4678 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4679
4680 /* Add two lo words. */
4681 t1 = tail_a + tail_b;
4682 bv = t1 - tail_a;
4683 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4684
4685 s2 += t1;
4686
4687 /* Renormalize (s1, s2) to (t1, s2) */
4688 t1 = s1 + s2;
4689 s2 = s2 - (t1 - s1);
4690
4691 t2 += s2;
4692
4693 /* Renormalize (t1, t2) */
4694 head_t = t1 + t2;
4695 tail_t = t2 - (head_t - t1);
4696 }
4697 head_sum[0] = head_t;
4698 tail_sum[0] = tail_t;
4699 /* Imaginary part */
4700 head_a = head_sum[1];
4701 tail_a = tail_sum[1];
4702 head_b = head_prod[1];
4703 tail_b = tail_prod[1];
4704 {
4705 /* Compute double-double = double-double + double-double. */
4706 double bv;
4707 double s1, s2, t1, t2;
4708
4709 /* Add two hi words. */
4710 s1 = head_a + head_b;
4711 bv = s1 - head_a;
4712 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4713
4714 /* Add two lo words. */
4715 t1 = tail_a + tail_b;
4716 bv = t1 - tail_a;
4717 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4718
4719 s2 += t1;
4720
4721 /* Renormalize (s1, s2) to (t1, s2) */
4722 t1 = s1 + s2;
4723 s2 = s2 - (t1 - s1);
4724
4725 t2 += s2;
4726
4727 /* Renormalize (t1, t2) */
4728 head_t = t1 + t2;
4729 tail_t = t2 - (head_t - t1);
4730 }
4731 head_sum[1] = head_t;
4732 tail_sum[1] = tail_t;
4733 }
4734 x_elem[0] = tail_x_i[jx];
4735 x_elem[1] = tail_x_i[jx + 1];
4736 {
4737 double head_e1, tail_e1;
4738 double d1;
4739 double d2;
4740 /* Real part */
4741 d1 = (double) a_elem[0] * x_elem[0];
4742 d2 = (double) -a_elem[1] * x_elem[1];
4743 {
4744 /* Compute double-double = double + double. */
4745 double e, t1, t2;
4746
4747 /* Knuth trick. */
4748 t1 = d1 + d2;
4749 e = t1 - d1;
4750 t2 = ((d2 - e) + (d1 - (t1 - e)));
4751
4752 /* The result is t1 + t2, after normalization. */
4753 head_e1 = t1 + t2;
4754 tail_e1 = t2 - (head_e1 - t1);
4755 }
4756 head_prod[0] = head_e1;
4757 tail_prod[0] = tail_e1;
4758 /* imaginary part */
4759 d1 = (double) a_elem[0] * x_elem[1];
4760 d2 = (double) a_elem[1] * x_elem[0];
4761 {
4762 /* Compute double-double = double + double. */
4763 double e, t1, t2;
4764
4765 /* Knuth trick. */
4766 t1 = d1 + d2;
4767 e = t1 - d1;
4768 t2 = ((d2 - e) + (d1 - (t1 - e)));
4769
4770 /* The result is t1 + t2, after normalization. */
4771 head_e1 = t1 + t2;
4772 tail_e1 = t2 - (head_e1 - t1);
4773 }
4774 head_prod[1] = head_e1;
4775 tail_prod[1] = tail_e1;
4776 }
4777 {
4778 double head_t, tail_t;
4779 double head_a, tail_a;
4780 double head_b, tail_b;
4781 /* Real part */
4782 head_a = head_sum2[0];
4783 tail_a = tail_sum2[0];
4784 head_b = head_prod[0];
4785 tail_b = tail_prod[0];
4786 {
4787 /* Compute double-double = double-double + double-double. */
4788 double bv;
4789 double s1, s2, t1, t2;
4790
4791 /* Add two hi words. */
4792 s1 = head_a + head_b;
4793 bv = s1 - head_a;
4794 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4795
4796 /* Add two lo words. */
4797 t1 = tail_a + tail_b;
4798 bv = t1 - tail_a;
4799 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4800
4801 s2 += t1;
4802
4803 /* Renormalize (s1, s2) to (t1, s2) */
4804 t1 = s1 + s2;
4805 s2 = s2 - (t1 - s1);
4806
4807 t2 += s2;
4808
4809 /* Renormalize (t1, t2) */
4810 head_t = t1 + t2;
4811 tail_t = t2 - (head_t - t1);
4812 }
4813 head_sum2[0] = head_t;
4814 tail_sum2[0] = tail_t;
4815 /* Imaginary part */
4816 head_a = head_sum2[1];
4817 tail_a = tail_sum2[1];
4818 head_b = head_prod[1];
4819 tail_b = tail_prod[1];
4820 {
4821 /* Compute double-double = double-double + double-double. */
4822 double bv;
4823 double s1, s2, t1, t2;
4824
4825 /* Add two hi words. */
4826 s1 = head_a + head_b;
4827 bv = s1 - head_a;
4828 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4829
4830 /* Add two lo words. */
4831 t1 = tail_a + tail_b;
4832 bv = t1 - tail_a;
4833 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4834
4835 s2 += t1;
4836
4837 /* Renormalize (s1, s2) to (t1, s2) */
4838 t1 = s1 + s2;
4839 s2 = s2 - (t1 - s1);
4840
4841 t2 += s2;
4842
4843 /* Renormalize (t1, t2) */
4844 head_t = t1 + t2;
4845 tail_t = t2 - (head_t - t1);
4846 }
4847 head_sum2[1] = head_t;
4848 tail_sum2[1] = tail_t;
4849 }
4850 aij += incaij;
4851 jx += incx;
4852 }
4853 {
4854 double head_t, tail_t;
4855 double head_a, tail_a;
4856 double head_b, tail_b;
4857 /* Real part */
4858 head_a = head_sum[0];
4859 tail_a = tail_sum[0];
4860 head_b = head_sum2[0];
4861 tail_b = tail_sum2[0];
4862 {
4863 /* Compute double-double = double-double + double-double. */
4864 double bv;
4865 double s1, s2, t1, t2;
4866
4867 /* Add two hi words. */
4868 s1 = head_a + head_b;
4869 bv = s1 - head_a;
4870 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4871
4872 /* Add two lo words. */
4873 t1 = tail_a + tail_b;
4874 bv = t1 - tail_a;
4875 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4876
4877 s2 += t1;
4878
4879 /* Renormalize (s1, s2) to (t1, s2) */
4880 t1 = s1 + s2;
4881 s2 = s2 - (t1 - s1);
4882
4883 t2 += s2;
4884
4885 /* Renormalize (t1, t2) */
4886 head_t = t1 + t2;
4887 tail_t = t2 - (head_t - t1);
4888 }
4889 head_sum[0] = head_t;
4890 tail_sum[0] = tail_t;
4891 /* Imaginary part */
4892 head_a = head_sum[1];
4893 tail_a = tail_sum[1];
4894 head_b = head_sum2[1];
4895 tail_b = tail_sum2[1];
4896 {
4897 /* Compute double-double = double-double + double-double. */
4898 double bv;
4899 double s1, s2, t1, t2;
4900
4901 /* Add two hi words. */
4902 s1 = head_a + head_b;
4903 bv = s1 - head_a;
4904 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4905
4906 /* Add two lo words. */
4907 t1 = tail_a + tail_b;
4908 bv = t1 - tail_a;
4909 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4910
4911 s2 += t1;
4912
4913 /* Renormalize (s1, s2) to (t1, s2) */
4914 t1 = s1 + s2;
4915 s2 = s2 - (t1 - s1);
4916
4917 t2 += s2;
4918
4919 /* Renormalize (t1, t2) */
4920 head_t = t1 + t2;
4921 tail_t = t2 - (head_t - t1);
4922 }
4923 head_sum[1] = head_t;
4924 tail_sum[1] = tail_t;
4925 }
4926 y_elem[0] = y_i[iy];
4927 y_elem[1] = y_i[iy + 1];
4928 {
4929 double head_e1, tail_e1;
4930 double d1;
4931 double d2;
4932 /* Real part */
4933 d1 = (double) y_elem[0] * beta_i[0];
4934 d2 = (double) -y_elem[1] * beta_i[1];
4935 {
4936 /* Compute double-double = double + double. */
4937 double e, t1, t2;
4938
4939 /* Knuth trick. */
4940 t1 = d1 + d2;
4941 e = t1 - d1;
4942 t2 = ((d2 - e) + (d1 - (t1 - e)));
4943
4944 /* The result is t1 + t2, after normalization. */
4945 head_e1 = t1 + t2;
4946 tail_e1 = t2 - (head_e1 - t1);
4947 }
4948 head_tmp1[0] = head_e1;
4949 tail_tmp1[0] = tail_e1;
4950 /* imaginary part */
4951 d1 = (double) y_elem[0] * beta_i[1];
4952 d2 = (double) y_elem[1] * beta_i[0];
4953 {
4954 /* Compute double-double = double + double. */
4955 double e, t1, t2;
4956
4957 /* Knuth trick. */
4958 t1 = d1 + d2;
4959 e = t1 - d1;
4960 t2 = ((d2 - e) + (d1 - (t1 - e)));
4961
4962 /* The result is t1 + t2, after normalization. */
4963 head_e1 = t1 + t2;
4964 tail_e1 = t2 - (head_e1 - t1);
4965 }
4966 head_tmp1[1] = head_e1;
4967 tail_tmp1[1] = tail_e1;
4968 }
4969 {
4970 double head_t, tail_t;
4971 double head_a, tail_a;
4972 double head_b, tail_b;
4973 /* Real part */
4974 head_a = head_sum[0];
4975 tail_a = tail_sum[0];
4976 head_b = head_tmp1[0];
4977 tail_b = tail_tmp1[0];
4978 {
4979 /* Compute double-double = double-double + double-double. */
4980 double bv;
4981 double s1, s2, t1, t2;
4982
4983 /* Add two hi words. */
4984 s1 = head_a + head_b;
4985 bv = s1 - head_a;
4986 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4987
4988 /* Add two lo words. */
4989 t1 = tail_a + tail_b;
4990 bv = t1 - tail_a;
4991 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4992
4993 s2 += t1;
4994
4995 /* Renormalize (s1, s2) to (t1, s2) */
4996 t1 = s1 + s2;
4997 s2 = s2 - (t1 - s1);
4998
4999 t2 += s2;
5000
5001 /* Renormalize (t1, t2) */
5002 head_t = t1 + t2;
5003 tail_t = t2 - (head_t - t1);
5004 }
5005 head_tmp2[0] = head_t;
5006 tail_tmp2[0] = tail_t;
5007 /* Imaginary part */
5008 head_a = head_sum[1];
5009 tail_a = tail_sum[1];
5010 head_b = head_tmp1[1];
5011 tail_b = tail_tmp1[1];
5012 {
5013 /* Compute double-double = double-double + double-double. */
5014 double bv;
5015 double s1, s2, t1, t2;
5016
5017 /* Add two hi words. */
5018 s1 = head_a + head_b;
5019 bv = s1 - head_a;
5020 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
5021
5022 /* Add two lo words. */
5023 t1 = tail_a + tail_b;
5024 bv = t1 - tail_a;
5025 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
5026
5027 s2 += t1;
5028
5029 /* Renormalize (s1, s2) to (t1, s2) */
5030 t1 = s1 + s2;
5031 s2 = s2 - (t1 - s1);
5032
5033 t2 += s2;
5034
5035 /* Renormalize (t1, t2) */
5036 head_t = t1 + t2;
5037 tail_t = t2 - (head_t - t1);
5038 }
5039 head_tmp2[1] = head_t;
5040 tail_tmp2[1] = tail_t;
5041 }
5042 y_i[iy] = head_tmp2[0];
5043 y_i[iy + 1] = head_tmp2[1];
5044 ai += incai;
5045 iy += incy;
5046 }
5047 } else { /* alpha != 1, the most general form:
5048 y = alpha*A*head_x + alpha*A*tail_x + beta*y */
5049 ai = 0;
5050 iy = ky;
5051 for (i = 0; i < leny; i++) {
5052 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;;
5053 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
5054 0.0;;
5055 aij = ai;
5056 jx = kx;
5057 for (j = 0; j < lenx; j++) {
5058 a_elem[0] = a_i[aij];
5059 a_elem[1] = a_i[aij + 1];
5060
5061 x_elem[0] = head_x_i[jx];
5062 x_elem[1] = head_x_i[jx + 1];
5063 {
5064 double head_e1, tail_e1;
5065 double d1;
5066 double d2;
5067 /* Real part */
5068 d1 = (double) a_elem[0] * x_elem[0];
5069 d2 = (double) -a_elem[1] * x_elem[1];
5070 {
5071 /* Compute double-double = double + double. */
5072 double e, t1, t2;
5073
5074 /* Knuth trick. */
5075 t1 = d1 + d2;
5076 e = t1 - d1;
5077 t2 = ((d2 - e) + (d1 - (t1 - e)));
5078
5079 /* The result is t1 + t2, after normalization. */
5080 head_e1 = t1 + t2;
5081 tail_e1 = t2 - (head_e1 - t1);
5082 }
5083 head_prod[0] = head_e1;
5084 tail_prod[0] = tail_e1;
5085 /* imaginary part */
5086 d1 = (double) a_elem[0] * x_elem[1];
5087 d2 = (double) a_elem[1] * x_elem[0];
5088 {
5089 /* Compute double-double = double + double. */
5090 double e, t1, t2;
5091
5092 /* Knuth trick. */
5093 t1 = d1 + d2;
5094 e = t1 - d1;
5095 t2 = ((d2 - e) + (d1 - (t1 - e)));
5096
5097 /* The result is t1 + t2, after normalization. */
5098 head_e1 = t1 + t2;
5099 tail_e1 = t2 - (head_e1 - t1);
5100 }
5101 head_prod[1] = head_e1;
5102 tail_prod[1] = tail_e1;
5103 }
5104 {
5105 double head_t, tail_t;
5106 double head_a, tail_a;
5107 double head_b, tail_b;
5108 /* Real part */
5109 head_a = head_sum[0];
5110 tail_a = tail_sum[0];
5111 head_b = head_prod[0];
5112 tail_b = tail_prod[0];
5113 {
5114 /* Compute double-double = double-double + double-double. */
5115 double bv;
5116 double s1, s2, t1, t2;
5117
5118 /* Add two hi words. */
5119 s1 = head_a + head_b;
5120 bv = s1 - head_a;
5121 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
5122
5123 /* Add two lo words. */
5124 t1 = tail_a + tail_b;
5125 bv = t1 - tail_a;
5126 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
5127
5128 s2 += t1;
5129
5130 /* Renormalize (s1, s2) to (t1, s2) */
5131 t1 = s1 + s2;
5132 s2 = s2 - (t1 - s1);
5133
5134 t2 += s2;
5135
5136 /* Renormalize (t1, t2) */
5137 head_t = t1 + t2;
5138 tail_t = t2 - (head_t - t1);
5139 }
5140 head_sum[0] = head_t;
5141 tail_sum[0] = tail_t;
5142 /* Imaginary part */
5143 head_a = head_sum[1];
5144 tail_a = tail_sum[1];
5145 head_b = head_prod[1];
5146 tail_b = tail_prod[1];
5147 {
5148 /* Compute double-double = double-double + double-double. */
5149 double bv;
5150 double s1, s2, t1, t2;
5151
5152 /* Add two hi words. */
5153 s1 = head_a + head_b;
5154 bv = s1 - head_a;
5155 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
5156
5157 /* Add two lo words. */
5158 t1 = tail_a + tail_b;
5159 bv = t1 - tail_a;
5160 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
5161
5162 s2 += t1;
5163
5164 /* Renormalize (s1, s2) to (t1, s2) */
5165 t1 = s1 + s2;
5166 s2 = s2 - (t1 - s1);
5167
5168 t2 += s2;
5169
5170 /* Renormalize (t1, t2) */
5171 head_t = t1 + t2;
5172 tail_t = t2 - (head_t - t1);
5173 }
5174 head_sum[1] = head_t;
5175 tail_sum[1] = tail_t;
5176 }
5177 x_elem[0] = tail_x_i[jx];
5178 x_elem[1] = tail_x_i[jx + 1];
5179 {
5180 double head_e1, tail_e1;
5181 double d1;
5182 double d2;
5183 /* Real part */
5184 d1 = (double) a_elem[0] * x_elem[0];
5185 d2 = (double) -a_elem[1] * x_elem[1];
5186 {
5187 /* Compute double-double = double + double. */
5188 double e, t1, t2;
5189
5190 /* Knuth trick. */
5191 t1 = d1 + d2;
5192 e = t1 - d1;
5193 t2 = ((d2 - e) + (d1 - (t1 - e)));
5194
5195 /* The result is t1 + t2, after normalization. */
5196 head_e1 = t1 + t2;
5197 tail_e1 = t2 - (head_e1 - t1);
5198 }
5199 head_prod[0] = head_e1;
5200 tail_prod[0] = tail_e1;
5201 /* imaginary part */
5202 d1 = (double) a_elem[0] * x_elem[1];
5203 d2 = (double) a_elem[1] * x_elem[0];
5204 {
5205 /* Compute double-double = double + double. */
5206 double e, t1, t2;
5207
5208 /* Knuth trick. */
5209 t1 = d1 + d2;
5210 e = t1 - d1;
5211 t2 = ((d2 - e) + (d1 - (t1 - e)));
5212
5213 /* The result is t1 + t2, after normalization. */
5214 head_e1 = t1 + t2;
5215 tail_e1 = t2 - (head_e1 - t1);
5216 }
5217 head_prod[1] = head_e1;
5218 tail_prod[1] = tail_e1;
5219 }
5220 {
5221 double head_t, tail_t;
5222 double head_a, tail_a;
5223 double head_b, tail_b;
5224 /* Real part */
5225 head_a = head_sum2[0];
5226 tail_a = tail_sum2[0];
5227 head_b = head_prod[0];
5228 tail_b = tail_prod[0];
5229 {
5230 /* Compute double-double = double-double + double-double. */
5231 double bv;
5232 double s1, s2, t1, t2;
5233
5234 /* Add two hi words. */
5235 s1 = head_a + head_b;
5236 bv = s1 - head_a;
5237 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
5238
5239 /* Add two lo words. */
5240 t1 = tail_a + tail_b;
5241 bv = t1 - tail_a;
5242 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
5243
5244 s2 += t1;
5245
5246 /* Renormalize (s1, s2) to (t1, s2) */
5247 t1 = s1 + s2;
5248 s2 = s2 - (t1 - s1);
5249
5250 t2 += s2;
5251
5252 /* Renormalize (t1, t2) */
5253 head_t = t1 + t2;
5254 tail_t = t2 - (head_t - t1);
5255 }
5256 head_sum2[0] = head_t;
5257 tail_sum2[0] = tail_t;
5258 /* Imaginary part */
5259 head_a = head_sum2[1];
5260 tail_a = tail_sum2[1];
5261 head_b = head_prod[1];
5262 tail_b = tail_prod[1];
5263 {
5264 /* Compute double-double = double-double + double-double. */
5265 double bv;
5266 double s1, s2, t1, t2;
5267
5268 /* Add two hi words. */
5269 s1 = head_a + head_b;
5270 bv = s1 - head_a;
5271 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
5272
5273 /* Add two lo words. */
5274 t1 = tail_a + tail_b;
5275 bv = t1 - tail_a;
5276 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
5277
5278 s2 += t1;
5279
5280 /* Renormalize (s1, s2) to (t1, s2) */
5281 t1 = s1 + s2;
5282 s2 = s2 - (t1 - s1);
5283
5284 t2 += s2;
5285
5286 /* Renormalize (t1, t2) */
5287 head_t = t1 + t2;
5288 tail_t = t2 - (head_t - t1);
5289 }
5290 head_sum2[1] = head_t;
5291 tail_sum2[1] = tail_t;
5292 }
5293 aij += incaij;
5294 jx += incx;
5295 }
5296 {
5297 double cd[2];
5298 cd[0] = (double) alpha_i[0];
5299 cd[1] = (double) alpha_i[1];
5300 {
5301 /* Compute complex-extra = complex-extra * complex-double. */
5302 double head_a0, tail_a0;
5303 double head_a1, tail_a1;
5304 double head_t1, tail_t1;
5305 double head_t2, tail_t2;
5306 head_a0 = head_sum[0];
5307 tail_a0 = tail_sum[0];
5308 head_a1 = head_sum[1];
5309 tail_a1 = tail_sum[1];
5310 /* real part */
5311 {
5312 /* Compute double-double = double-double * double. */
5313 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
5314
5315 con = head_a0 * split;
5316 a11 = con - head_a0;
5317 a11 = con - a11;
5318 a21 = head_a0 - a11;
5319 con = cd[0] * split;
5320 b1 = con - cd[0];
5321 b1 = con - b1;
5322 b2 = cd[0] - b1;
5323
5324 c11 = head_a0 * cd[0];
5325 c21 =
5326 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
5327
5328 c2 = tail_a0 * cd[0];
5329 t1 = c11 + c2;
5330 t2 = (c2 - (t1 - c11)) + c21;
5331
5332 head_t1 = t1 + t2;
5333 tail_t1 = t2 - (head_t1 - t1);
5334 }
5335 {
5336 /* Compute double-double = double-double * double. */
5337 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
5338
5339 con = head_a1 * split;
5340 a11 = con - head_a1;
5341 a11 = con - a11;
5342 a21 = head_a1 - a11;
5343 con = cd[1] * split;
5344 b1 = con - cd[1];
5345 b1 = con - b1;
5346 b2 = cd[1] - b1;
5347
5348 c11 = head_a1 * cd[1];
5349 c21 =
5350 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
5351
5352 c2 = tail_a1 * cd[1];
5353 t1 = c11 + c2;
5354 t2 = (c2 - (t1 - c11)) + c21;
5355
5356 head_t2 = t1 + t2;
5357 tail_t2 = t2 - (head_t2 - t1);
5358 }
5359 head_t2 = -head_t2;
5360 tail_t2 = -tail_t2;
5361 {
5362 /* Compute double-double = double-double + double-double. */
5363 double bv;
5364 double s1, s2, t1, t2;
5365
5366 /* Add two hi words. */
5367 s1 = head_t1 + head_t2;
5368 bv = s1 - head_t1;
5369 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
5370
5371 /* Add two lo words. */
5372 t1 = tail_t1 + tail_t2;
5373 bv = t1 - tail_t1;
5374 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
5375
5376 s2 += t1;
5377
5378 /* Renormalize (s1, s2) to (t1, s2) */
5379 t1 = s1 + s2;
5380 s2 = s2 - (t1 - s1);
5381
5382 t2 += s2;
5383
5384 /* Renormalize (t1, t2) */
5385 head_t1 = t1 + t2;
5386 tail_t1 = t2 - (head_t1 - t1);
5387 }
5388 head_tmp1[0] = head_t1;
5389 tail_tmp1[0] = tail_t1;
5390 /* imaginary part */
5391 {
5392 /* Compute double-double = double-double * double. */
5393 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
5394
5395 con = head_a1 * split;
5396 a11 = con - head_a1;
5397 a11 = con - a11;
5398 a21 = head_a1 - a11;
5399 con = cd[0] * split;
5400 b1 = con - cd[0];
5401 b1 = con - b1;
5402 b2 = cd[0] - b1;
5403
5404 c11 = head_a1 * cd[0];
5405 c21 =
5406 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
5407
5408 c2 = tail_a1 * cd[0];
5409 t1 = c11 + c2;
5410 t2 = (c2 - (t1 - c11)) + c21;
5411
5412 head_t1 = t1 + t2;
5413 tail_t1 = t2 - (head_t1 - t1);
5414 }
5415 {
5416 /* Compute double-double = double-double * double. */
5417 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
5418
5419 con = head_a0 * split;
5420 a11 = con - head_a0;
5421 a11 = con - a11;
5422 a21 = head_a0 - a11;
5423 con = cd[1] * split;
5424 b1 = con - cd[1];
5425 b1 = con - b1;
5426 b2 = cd[1] - b1;
5427
5428 c11 = head_a0 * cd[1];
5429 c21 =
5430 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
5431
5432 c2 = tail_a0 * cd[1];
5433 t1 = c11 + c2;
5434 t2 = (c2 - (t1 - c11)) + c21;
5435
5436 head_t2 = t1 + t2;
5437 tail_t2 = t2 - (head_t2 - t1);
5438 }
5439 {
5440 /* Compute double-double = double-double + double-double. */
5441 double bv;
5442 double s1, s2, t1, t2;
5443
5444 /* Add two hi words. */
5445 s1 = head_t1 + head_t2;
5446 bv = s1 - head_t1;
5447 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
5448
5449 /* Add two lo words. */
5450 t1 = tail_t1 + tail_t2;
5451 bv = t1 - tail_t1;
5452 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
5453
5454 s2 += t1;
5455
5456 /* Renormalize (s1, s2) to (t1, s2) */
5457 t1 = s1 + s2;
5458 s2 = s2 - (t1 - s1);
5459
5460 t2 += s2;
5461
5462 /* Renormalize (t1, t2) */
5463 head_t1 = t1 + t2;
5464 tail_t1 = t2 - (head_t1 - t1);
5465 }
5466 head_tmp1[1] = head_t1;
5467 tail_tmp1[1] = tail_t1;
5468 }
5469
5470 }
5471 {
5472 double cd[2];
5473 cd[0] = (double) alpha_i[0];
5474 cd[1] = (double) alpha_i[1];
5475 {
5476 /* Compute complex-extra = complex-extra * complex-double. */
5477 double head_a0, tail_a0;
5478 double head_a1, tail_a1;
5479 double head_t1, tail_t1;
5480 double head_t2, tail_t2;
5481 head_a0 = head_sum2[0];
5482 tail_a0 = tail_sum2[0];
5483 head_a1 = head_sum2[1];
5484 tail_a1 = tail_sum2[1];
5485 /* real part */
5486 {
5487 /* Compute double-double = double-double * double. */
5488 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
5489
5490 con = head_a0 * split;
5491 a11 = con - head_a0;
5492 a11 = con - a11;
5493 a21 = head_a0 - a11;
5494 con = cd[0] * split;
5495 b1 = con - cd[0];
5496 b1 = con - b1;
5497 b2 = cd[0] - b1;
5498
5499 c11 = head_a0 * cd[0];
5500 c21 =
5501 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
5502
5503 c2 = tail_a0 * cd[0];
5504 t1 = c11 + c2;
5505 t2 = (c2 - (t1 - c11)) + c21;
5506
5507 head_t1 = t1 + t2;
5508 tail_t1 = t2 - (head_t1 - t1);
5509 }
5510 {
5511 /* Compute double-double = double-double * double. */
5512 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
5513
5514 con = head_a1 * split;
5515 a11 = con - head_a1;
5516 a11 = con - a11;
5517 a21 = head_a1 - a11;
5518 con = cd[1] * split;
5519 b1 = con - cd[1];
5520 b1 = con - b1;
5521 b2 = cd[1] - b1;
5522
5523 c11 = head_a1 * cd[1];
5524 c21 =
5525 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
5526
5527 c2 = tail_a1 * cd[1];
5528 t1 = c11 + c2;
5529 t2 = (c2 - (t1 - c11)) + c21;
5530
5531 head_t2 = t1 + t2;
5532 tail_t2 = t2 - (head_t2 - t1);
5533 }
5534 head_t2 = -head_t2;
5535 tail_t2 = -tail_t2;
5536 {
5537 /* Compute double-double = double-double + double-double. */
5538 double bv;
5539 double s1, s2, t1, t2;
5540
5541 /* Add two hi words. */
5542 s1 = head_t1 + head_t2;
5543 bv = s1 - head_t1;
5544 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
5545
5546 /* Add two lo words. */
5547 t1 = tail_t1 + tail_t2;
5548 bv = t1 - tail_t1;
5549 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
5550
5551 s2 += t1;
5552
5553 /* Renormalize (s1, s2) to (t1, s2) */
5554 t1 = s1 + s2;
5555 s2 = s2 - (t1 - s1);
5556
5557 t2 += s2;
5558
5559 /* Renormalize (t1, t2) */
5560 head_t1 = t1 + t2;
5561 tail_t1 = t2 - (head_t1 - t1);
5562 }
5563 head_tmp2[0] = head_t1;
5564 tail_tmp2[0] = tail_t1;
5565 /* imaginary part */
5566 {
5567 /* Compute double-double = double-double * double. */
5568 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
5569
5570 con = head_a1 * split;
5571 a11 = con - head_a1;
5572 a11 = con - a11;
5573 a21 = head_a1 - a11;
5574 con = cd[0] * split;
5575 b1 = con - cd[0];
5576 b1 = con - b1;
5577 b2 = cd[0] - b1;
5578
5579 c11 = head_a1 * cd[0];
5580 c21 =
5581 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
5582
5583 c2 = tail_a1 * cd[0];
5584 t1 = c11 + c2;
5585 t2 = (c2 - (t1 - c11)) + c21;
5586
5587 head_t1 = t1 + t2;
5588 tail_t1 = t2 - (head_t1 - t1);
5589 }
5590 {
5591 /* Compute double-double = double-double * double. */
5592 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
5593
5594 con = head_a0 * split;
5595 a11 = con - head_a0;
5596 a11 = con - a11;
5597 a21 = head_a0 - a11;
5598 con = cd[1] * split;
5599 b1 = con - cd[1];
5600 b1 = con - b1;
5601 b2 = cd[1] - b1;
5602
5603 c11 = head_a0 * cd[1];
5604 c21 =
5605 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
5606
5607 c2 = tail_a0 * cd[1];
5608 t1 = c11 + c2;
5609 t2 = (c2 - (t1 - c11)) + c21;
5610
5611 head_t2 = t1 + t2;
5612 tail_t2 = t2 - (head_t2 - t1);
5613 }
5614 {
5615 /* Compute double-double = double-double + double-double. */
5616 double bv;
5617 double s1, s2, t1, t2;
5618
5619 /* Add two hi words. */
5620 s1 = head_t1 + head_t2;
5621 bv = s1 - head_t1;
5622 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
5623
5624 /* Add two lo words. */
5625 t1 = tail_t1 + tail_t2;
5626 bv = t1 - tail_t1;
5627 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
5628
5629 s2 += t1;
5630
5631 /* Renormalize (s1, s2) to (t1, s2) */
5632 t1 = s1 + s2;
5633 s2 = s2 - (t1 - s1);
5634
5635 t2 += s2;
5636
5637 /* Renormalize (t1, t2) */
5638 head_t1 = t1 + t2;
5639 tail_t1 = t2 - (head_t1 - t1);
5640 }
5641 head_tmp2[1] = head_t1;
5642 tail_tmp2[1] = tail_t1;
5643 }
5644
5645 }
5646 {
5647 double head_t, tail_t;
5648 double head_a, tail_a;
5649 double head_b, tail_b;
5650 /* Real part */
5651 head_a = head_tmp1[0];
5652 tail_a = tail_tmp1[0];
5653 head_b = head_tmp2[0];
5654 tail_b = tail_tmp2[0];
5655 {
5656 /* Compute double-double = double-double + double-double. */
5657 double bv;
5658 double s1, s2, t1, t2;
5659
5660 /* Add two hi words. */
5661 s1 = head_a + head_b;
5662 bv = s1 - head_a;
5663 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
5664
5665 /* Add two lo words. */
5666 t1 = tail_a + tail_b;
5667 bv = t1 - tail_a;
5668 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
5669
5670 s2 += t1;
5671
5672 /* Renormalize (s1, s2) to (t1, s2) */
5673 t1 = s1 + s2;
5674 s2 = s2 - (t1 - s1);
5675
5676 t2 += s2;
5677
5678 /* Renormalize (t1, t2) */
5679 head_t = t1 + t2;
5680 tail_t = t2 - (head_t - t1);
5681 }
5682 head_tmp1[0] = head_t;
5683 tail_tmp1[0] = tail_t;
5684 /* Imaginary part */
5685 head_a = head_tmp1[1];
5686 tail_a = tail_tmp1[1];
5687 head_b = head_tmp2[1];
5688 tail_b = tail_tmp2[1];
5689 {
5690 /* Compute double-double = double-double + double-double. */
5691 double bv;
5692 double s1, s2, t1, t2;
5693
5694 /* Add two hi words. */
5695 s1 = head_a + head_b;
5696 bv = s1 - head_a;
5697 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
5698
5699 /* Add two lo words. */
5700 t1 = tail_a + tail_b;
5701 bv = t1 - tail_a;
5702 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
5703
5704 s2 += t1;
5705
5706 /* Renormalize (s1, s2) to (t1, s2) */
5707 t1 = s1 + s2;
5708 s2 = s2 - (t1 - s1);
5709
5710 t2 += s2;
5711
5712 /* Renormalize (t1, t2) */
5713 head_t = t1 + t2;
5714 tail_t = t2 - (head_t - t1);
5715 }
5716 head_tmp1[1] = head_t;
5717 tail_tmp1[1] = tail_t;
5718 }
5719 y_elem[0] = y_i[iy];
5720 y_elem[1] = y_i[iy + 1];
5721 {
5722 double head_e1, tail_e1;
5723 double d1;
5724 double d2;
5725 /* Real part */
5726 d1 = (double) y_elem[0] * beta_i[0];
5727 d2 = (double) -y_elem[1] * beta_i[1];
5728 {
5729 /* Compute double-double = double + double. */
5730 double e, t1, t2;
5731
5732 /* Knuth trick. */
5733 t1 = d1 + d2;
5734 e = t1 - d1;
5735 t2 = ((d2 - e) + (d1 - (t1 - e)));
5736
5737 /* The result is t1 + t2, after normalization. */
5738 head_e1 = t1 + t2;
5739 tail_e1 = t2 - (head_e1 - t1);
5740 }
5741 head_tmp2[0] = head_e1;
5742 tail_tmp2[0] = tail_e1;
5743 /* imaginary part */
5744 d1 = (double) y_elem[0] * beta_i[1];
5745 d2 = (double) y_elem[1] * beta_i[0];
5746 {
5747 /* Compute double-double = double + double. */
5748 double e, t1, t2;
5749
5750 /* Knuth trick. */
5751 t1 = d1 + d2;
5752 e = t1 - d1;
5753 t2 = ((d2 - e) + (d1 - (t1 - e)));
5754
5755 /* The result is t1 + t2, after normalization. */
5756 head_e1 = t1 + t2;
5757 tail_e1 = t2 - (head_e1 - t1);
5758 }
5759 head_tmp2[1] = head_e1;
5760 tail_tmp2[1] = tail_e1;
5761 }
5762 {
5763 double head_t, tail_t;
5764 double head_a, tail_a;
5765 double head_b, tail_b;
5766 /* Real part */
5767 head_a = head_tmp1[0];
5768 tail_a = tail_tmp1[0];
5769 head_b = head_tmp2[0];
5770 tail_b = tail_tmp2[0];
5771 {
5772 /* Compute double-double = double-double + double-double. */
5773 double bv;
5774 double s1, s2, t1, t2;
5775
5776 /* Add two hi words. */
5777 s1 = head_a + head_b;
5778 bv = s1 - head_a;
5779 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
5780
5781 /* Add two lo words. */
5782 t1 = tail_a + tail_b;
5783 bv = t1 - tail_a;
5784 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
5785
5786 s2 += t1;
5787
5788 /* Renormalize (s1, s2) to (t1, s2) */
5789 t1 = s1 + s2;
5790 s2 = s2 - (t1 - s1);
5791
5792 t2 += s2;
5793
5794 /* Renormalize (t1, t2) */
5795 head_t = t1 + t2;
5796 tail_t = t2 - (head_t - t1);
5797 }
5798 head_tmp1[0] = head_t;
5799 tail_tmp1[0] = tail_t;
5800 /* Imaginary part */
5801 head_a = head_tmp1[1];
5802 tail_a = tail_tmp1[1];
5803 head_b = head_tmp2[1];
5804 tail_b = tail_tmp2[1];
5805 {
5806 /* Compute double-double = double-double + double-double. */
5807 double bv;
5808 double s1, s2, t1, t2;
5809
5810 /* Add two hi words. */
5811 s1 = head_a + head_b;
5812 bv = s1 - head_a;
5813 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
5814
5815 /* Add two lo words. */
5816 t1 = tail_a + tail_b;
5817 bv = t1 - tail_a;
5818 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
5819
5820 s2 += t1;
5821
5822 /* Renormalize (s1, s2) to (t1, s2) */
5823 t1 = s1 + s2;
5824 s2 = s2 - (t1 - s1);
5825
5826 t2 += s2;
5827
5828 /* Renormalize (t1, t2) */
5829 head_t = t1 + t2;
5830 tail_t = t2 - (head_t - t1);
5831 }
5832 head_tmp1[1] = head_t;
5833 tail_tmp1[1] = tail_t;
5834 }
5835 y_i[iy] = head_tmp1[0];
5836 y_i[iy + 1] = head_tmp1[1];
5837 ai += incai;
5838 iy += incy;
5839 }
5840 }
5841 }
5842
5843 }
5844 }
5845
5846 FPU_FIX_STOP;
5847 }
5848 break;
5849 }
5850 }
5851