1 #include "blas_extended.h"
2 #include "blas_extended_private.h"
BLAS_cdot_x(enum blas_conj_type conj,int n,const void * alpha,const void * x,int incx,const void * beta,const void * y,int incy,void * r,enum blas_prec_type prec)3 void BLAS_cdot_x(enum blas_conj_type conj, int n, const void *alpha,
4 const void *x, int incx, const void *beta,
5 const void *y, int incy, void *r, enum blas_prec_type prec)
6
7 /*
8 * Purpose
9 * =======
10 *
11 * This routine computes the inner product:
12 *
13 * r <- beta * r + alpha * SUM_{i=0, n-1} x[i] * y[i].
14 *
15 * Arguments
16 * =========
17 *
18 * conj (input) enum blas_conj_type
19 * When x and y are complex vectors, specifies whether vector
20 * components x[i] are used unconjugated or conjugated.
21 *
22 * n (input) int
23 * The length of vectors x and y.
24 *
25 * alpha (input) const void*
26 *
27 * x (input) const void*
28 * Array of length n.
29 *
30 * incx (input) int
31 * The stride used to access components x[i].
32 *
33 * beta (input) const void*
34 *
35 * y (input) const void*
36 * Array of length n.
37 *
38 * incy (input) int
39 * The stride used to access components y[i].
40 *
41 * r (input/output) void*
42 *
43 * prec (input) enum blas_prec_type
44 * Specifies the internal precision to be used.
45 * = blas_prec_single: single precision.
46 * = blas_prec_double: double precision.
47 * = blas_prec_extra : anything at least 1.5 times as accurate
48 * than double, and wider than 80-bits.
49 * We use double-double in our implementation.
50 *
51 */
52 {
53 static const char routine_name[] = "BLAS_cdot_x";
54
55 switch (prec) {
56 case blas_prec_single:{
57
58 int i, ix = 0, iy = 0;
59 float *r_i = (float *) r;
60 const float *x_i = (float *) x;
61 const float *y_i = (float *) y;
62 float *alpha_i = (float *) alpha;
63 float *beta_i = (float *) beta;
64 float x_ii[2];
65 float y_ii[2];
66 float r_v[2];
67 float prod[2];
68 float sum[2];
69 float tmp1[2];
70 float tmp2[2];
71
72
73 /* Test the input parameters. */
74 if (n < 0)
75 BLAS_error(routine_name, -2, n, NULL);
76 else if (incx == 0)
77 BLAS_error(routine_name, -5, incx, NULL);
78 else if (incy == 0)
79 BLAS_error(routine_name, -8, incy, NULL);
80
81 /* Immediate return. */
82 if (((beta_i[0] == 1.0 && beta_i[1] == 0.0))
83 && (n == 0 || (alpha_i[0] == 0.0 && alpha_i[1] == 0.0)))
84 return;
85
86
87
88 r_v[0] = r_i[0];
89 r_v[1] = r_i[0 + 1];
90 sum[0] = sum[1] = 0.0;
91 incx *= 2;
92 incy *= 2;
93 if (incx < 0)
94 ix = (-n + 1) * incx;
95 if (incy < 0)
96 iy = (-n + 1) * incy;
97
98 if (conj == blas_conj) {
99 for (i = 0; i < n; ++i) {
100 x_ii[0] = x_i[ix];
101 x_ii[1] = x_i[ix + 1];
102 y_ii[0] = y_i[iy];
103 y_ii[1] = y_i[iy + 1];
104 x_ii[1] = -x_ii[1];
105 {
106 prod[0] = x_ii[0] * y_ii[0] - x_ii[1] * y_ii[1];
107 prod[1] = x_ii[0] * y_ii[1] + x_ii[1] * y_ii[0];
108 }
109 /* prod = x[i]*y[i] */
110 sum[0] = sum[0] + prod[0];
111 sum[1] = sum[1] + prod[1]; /* sum = sum+prod */
112 ix += incx;
113 iy += incy;
114 } /* endfor */
115 } else {
116 /* do not conjugate */
117
118 for (i = 0; i < n; ++i) {
119 x_ii[0] = x_i[ix];
120 x_ii[1] = x_i[ix + 1];
121 y_ii[0] = y_i[iy];
122 y_ii[1] = y_i[iy + 1];
123
124 {
125 prod[0] = x_ii[0] * y_ii[0] - x_ii[1] * y_ii[1];
126 prod[1] = x_ii[0] * y_ii[1] + x_ii[1] * y_ii[0];
127 }
128 /* prod = x[i]*y[i] */
129 sum[0] = sum[0] + prod[0];
130 sum[1] = sum[1] + prod[1]; /* sum = sum+prod */
131 ix += incx;
132 iy += incy;
133 } /* endfor */
134 }
135
136 {
137 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
138 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
139 }
140 /* tmp1 = sum*alpha */
141 {
142 tmp2[0] = r_v[0] * beta_i[0] - r_v[1] * beta_i[1];
143 tmp2[1] = r_v[0] * beta_i[1] + r_v[1] * beta_i[0];
144 }
145 /* tmp2 = r*beta */
146 tmp1[0] = tmp1[0] + tmp2[0];
147 tmp1[1] = tmp1[1] + tmp2[1]; /* tmp1 = tmp1+tmp2 */
148 ((float *) r)[0] = tmp1[0];
149 ((float *) r)[1] = tmp1[1]; /* r = tmp1 */
150
151
152
153 break;
154 }
155 case blas_prec_double:
156 case blas_prec_indigenous:
157 {
158 int i, ix = 0, iy = 0;
159 float *r_i = (float *) r;
160 const float *x_i = (float *) x;
161 const float *y_i = (float *) y;
162 float *alpha_i = (float *) alpha;
163 float *beta_i = (float *) beta;
164 float x_ii[2];
165 float y_ii[2];
166 float r_v[2];
167 double prod[2];
168 double sum[2];
169 double tmp1[2];
170 double tmp2[2];
171
172
173 /* Test the input parameters. */
174 if (n < 0)
175 BLAS_error(routine_name, -2, n, NULL);
176 else if (incx == 0)
177 BLAS_error(routine_name, -5, incx, NULL);
178 else if (incy == 0)
179 BLAS_error(routine_name, -8, incy, NULL);
180
181 /* Immediate return. */
182 if (((beta_i[0] == 1.0 && beta_i[1] == 0.0))
183 && (n == 0 || (alpha_i[0] == 0.0 && alpha_i[1] == 0.0)))
184 return;
185
186
187
188 r_v[0] = r_i[0];
189 r_v[1] = r_i[0 + 1];
190 sum[0] = sum[1] = 0.0;
191 incx *= 2;
192 incy *= 2;
193 if (incx < 0)
194 ix = (-n + 1) * incx;
195 if (incy < 0)
196 iy = (-n + 1) * incy;
197
198 if (conj == blas_conj) {
199 for (i = 0; i < n; ++i) {
200 x_ii[0] = x_i[ix];
201 x_ii[1] = x_i[ix + 1];
202 y_ii[0] = y_i[iy];
203 y_ii[1] = y_i[iy + 1];
204 x_ii[1] = -x_ii[1];
205 {
206 prod[0] = (double) x_ii[0] * y_ii[0] - (double) x_ii[1] * y_ii[1];
207 prod[1] = (double) x_ii[0] * y_ii[1] + (double) x_ii[1] * y_ii[0];
208 } /* prod = x[i]*y[i] */
209 sum[0] = sum[0] + prod[0];
210 sum[1] = sum[1] + prod[1]; /* sum = sum+prod */
211 ix += incx;
212 iy += incy;
213 } /* endfor */
214 } else {
215 /* do not conjugate */
216
217 for (i = 0; i < n; ++i) {
218 x_ii[0] = x_i[ix];
219 x_ii[1] = x_i[ix + 1];
220 y_ii[0] = y_i[iy];
221 y_ii[1] = y_i[iy + 1];
222
223 {
224 prod[0] = (double) x_ii[0] * y_ii[0] - (double) x_ii[1] * y_ii[1];
225 prod[1] = (double) x_ii[0] * y_ii[1] + (double) x_ii[1] * y_ii[0];
226 } /* prod = x[i]*y[i] */
227 sum[0] = sum[0] + prod[0];
228 sum[1] = sum[1] + prod[1]; /* sum = sum+prod */
229 ix += incx;
230 iy += incy;
231 } /* endfor */
232 }
233
234 {
235 tmp1[0] = (double) sum[0] * alpha_i[0] - (double) sum[1] * alpha_i[1];
236 tmp1[1] = (double) sum[0] * alpha_i[1] + (double) sum[1] * alpha_i[0];
237 } /* tmp1 = sum*alpha */
238 {
239 tmp2[0] = (double) r_v[0] * beta_i[0] - (double) r_v[1] * beta_i[1];
240 tmp2[1] = (double) r_v[0] * beta_i[1] + (double) r_v[1] * beta_i[0];
241 } /* tmp2 = r*beta */
242 tmp1[0] = tmp1[0] + tmp2[0];
243 tmp1[1] = tmp1[1] + tmp2[1]; /* tmp1 = tmp1+tmp2 */
244 ((float *) r)[0] = tmp1[0];
245 ((float *) r)[1] = tmp1[1]; /* r = tmp1 */
246
247
248 }
249 break;
250 case blas_prec_extra:
251 {
252 int i, ix = 0, iy = 0;
253 float *r_i = (float *) r;
254 const float *x_i = (float *) x;
255 const float *y_i = (float *) y;
256 float *alpha_i = (float *) alpha;
257 float *beta_i = (float *) beta;
258 float x_ii[2];
259 float y_ii[2];
260 float r_v[2];
261 double head_prod[2], tail_prod[2];
262 double head_sum[2], tail_sum[2];
263 double head_tmp1[2], tail_tmp1[2];
264 double head_tmp2[2], tail_tmp2[2];
265 FPU_FIX_DECL;
266
267 /* Test the input parameters. */
268 if (n < 0)
269 BLAS_error(routine_name, -2, n, NULL);
270 else if (incx == 0)
271 BLAS_error(routine_name, -5, incx, NULL);
272 else if (incy == 0)
273 BLAS_error(routine_name, -8, incy, NULL);
274
275 /* Immediate return. */
276 if (((beta_i[0] == 1.0 && beta_i[1] == 0.0))
277 && (n == 0 || (alpha_i[0] == 0.0 && alpha_i[1] == 0.0)))
278 return;
279
280 FPU_FIX_START;
281
282 r_v[0] = r_i[0];
283 r_v[1] = r_i[0 + 1];
284 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
285 incx *= 2;
286 incy *= 2;
287 if (incx < 0)
288 ix = (-n + 1) * incx;
289 if (incy < 0)
290 iy = (-n + 1) * incy;
291
292 if (conj == blas_conj) {
293 for (i = 0; i < n; ++i) {
294 x_ii[0] = x_i[ix];
295 x_ii[1] = x_i[ix + 1];
296 y_ii[0] = y_i[iy];
297 y_ii[1] = y_i[iy + 1];
298 x_ii[1] = -x_ii[1];
299 {
300 double head_e1, tail_e1;
301 double d1;
302 double d2;
303 /* Real part */
304 d1 = (double) x_ii[0] * y_ii[0];
305 d2 = (double) -x_ii[1] * y_ii[1];
306 {
307 /* Compute double-double = double + double. */
308 double e, t1, t2;
309
310 /* Knuth trick. */
311 t1 = d1 + d2;
312 e = t1 - d1;
313 t2 = ((d2 - e) + (d1 - (t1 - e)));
314
315 /* The result is t1 + t2, after normalization. */
316 head_e1 = t1 + t2;
317 tail_e1 = t2 - (head_e1 - t1);
318 }
319 head_prod[0] = head_e1;
320 tail_prod[0] = tail_e1;
321 /* imaginary part */
322 d1 = (double) x_ii[0] * y_ii[1];
323 d2 = (double) x_ii[1] * y_ii[0];
324 {
325 /* Compute double-double = double + double. */
326 double e, t1, t2;
327
328 /* Knuth trick. */
329 t1 = d1 + d2;
330 e = t1 - d1;
331 t2 = ((d2 - e) + (d1 - (t1 - e)));
332
333 /* The result is t1 + t2, after normalization. */
334 head_e1 = t1 + t2;
335 tail_e1 = t2 - (head_e1 - t1);
336 }
337 head_prod[1] = head_e1;
338 tail_prod[1] = tail_e1;
339 } /* prod = x[i]*y[i] */
340 {
341 double head_t, tail_t;
342 double head_a, tail_a;
343 double head_b, tail_b;
344 /* Real part */
345 head_a = head_sum[0];
346 tail_a = tail_sum[0];
347 head_b = head_prod[0];
348 tail_b = tail_prod[0];
349 {
350 /* Compute double-double = double-double + double-double. */
351 double bv;
352 double s1, s2, t1, t2;
353
354 /* Add two hi words. */
355 s1 = head_a + head_b;
356 bv = s1 - head_a;
357 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
358
359 /* Add two lo words. */
360 t1 = tail_a + tail_b;
361 bv = t1 - tail_a;
362 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
363
364 s2 += t1;
365
366 /* Renormalize (s1, s2) to (t1, s2) */
367 t1 = s1 + s2;
368 s2 = s2 - (t1 - s1);
369
370 t2 += s2;
371
372 /* Renormalize (t1, t2) */
373 head_t = t1 + t2;
374 tail_t = t2 - (head_t - t1);
375 }
376 head_sum[0] = head_t;
377 tail_sum[0] = tail_t;
378 /* Imaginary part */
379 head_a = head_sum[1];
380 tail_a = tail_sum[1];
381 head_b = head_prod[1];
382 tail_b = tail_prod[1];
383 {
384 /* Compute double-double = double-double + double-double. */
385 double bv;
386 double s1, s2, t1, t2;
387
388 /* Add two hi words. */
389 s1 = head_a + head_b;
390 bv = s1 - head_a;
391 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
392
393 /* Add two lo words. */
394 t1 = tail_a + tail_b;
395 bv = t1 - tail_a;
396 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
397
398 s2 += t1;
399
400 /* Renormalize (s1, s2) to (t1, s2) */
401 t1 = s1 + s2;
402 s2 = s2 - (t1 - s1);
403
404 t2 += s2;
405
406 /* Renormalize (t1, t2) */
407 head_t = t1 + t2;
408 tail_t = t2 - (head_t - t1);
409 }
410 head_sum[1] = head_t;
411 tail_sum[1] = tail_t;
412 } /* sum = sum+prod */
413 ix += incx;
414 iy += incy;
415 } /* endfor */
416 } else {
417 /* do not conjugate */
418
419 for (i = 0; i < n; ++i) {
420 x_ii[0] = x_i[ix];
421 x_ii[1] = x_i[ix + 1];
422 y_ii[0] = y_i[iy];
423 y_ii[1] = y_i[iy + 1];
424
425 {
426 double head_e1, tail_e1;
427 double d1;
428 double d2;
429 /* Real part */
430 d1 = (double) x_ii[0] * y_ii[0];
431 d2 = (double) -x_ii[1] * y_ii[1];
432 {
433 /* Compute double-double = double + double. */
434 double e, t1, t2;
435
436 /* Knuth trick. */
437 t1 = d1 + d2;
438 e = t1 - d1;
439 t2 = ((d2 - e) + (d1 - (t1 - e)));
440
441 /* The result is t1 + t2, after normalization. */
442 head_e1 = t1 + t2;
443 tail_e1 = t2 - (head_e1 - t1);
444 }
445 head_prod[0] = head_e1;
446 tail_prod[0] = tail_e1;
447 /* imaginary part */
448 d1 = (double) x_ii[0] * y_ii[1];
449 d2 = (double) x_ii[1] * y_ii[0];
450 {
451 /* Compute double-double = double + double. */
452 double e, t1, t2;
453
454 /* Knuth trick. */
455 t1 = d1 + d2;
456 e = t1 - d1;
457 t2 = ((d2 - e) + (d1 - (t1 - e)));
458
459 /* The result is t1 + t2, after normalization. */
460 head_e1 = t1 + t2;
461 tail_e1 = t2 - (head_e1 - t1);
462 }
463 head_prod[1] = head_e1;
464 tail_prod[1] = tail_e1;
465 } /* prod = x[i]*y[i] */
466 {
467 double head_t, tail_t;
468 double head_a, tail_a;
469 double head_b, tail_b;
470 /* Real part */
471 head_a = head_sum[0];
472 tail_a = tail_sum[0];
473 head_b = head_prod[0];
474 tail_b = tail_prod[0];
475 {
476 /* Compute double-double = double-double + double-double. */
477 double bv;
478 double s1, s2, t1, t2;
479
480 /* Add two hi words. */
481 s1 = head_a + head_b;
482 bv = s1 - head_a;
483 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
484
485 /* Add two lo words. */
486 t1 = tail_a + tail_b;
487 bv = t1 - tail_a;
488 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
489
490 s2 += t1;
491
492 /* Renormalize (s1, s2) to (t1, s2) */
493 t1 = s1 + s2;
494 s2 = s2 - (t1 - s1);
495
496 t2 += s2;
497
498 /* Renormalize (t1, t2) */
499 head_t = t1 + t2;
500 tail_t = t2 - (head_t - t1);
501 }
502 head_sum[0] = head_t;
503 tail_sum[0] = tail_t;
504 /* Imaginary part */
505 head_a = head_sum[1];
506 tail_a = tail_sum[1];
507 head_b = head_prod[1];
508 tail_b = tail_prod[1];
509 {
510 /* Compute double-double = double-double + double-double. */
511 double bv;
512 double s1, s2, t1, t2;
513
514 /* Add two hi words. */
515 s1 = head_a + head_b;
516 bv = s1 - head_a;
517 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
518
519 /* Add two lo words. */
520 t1 = tail_a + tail_b;
521 bv = t1 - tail_a;
522 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
523
524 s2 += t1;
525
526 /* Renormalize (s1, s2) to (t1, s2) */
527 t1 = s1 + s2;
528 s2 = s2 - (t1 - s1);
529
530 t2 += s2;
531
532 /* Renormalize (t1, t2) */
533 head_t = t1 + t2;
534 tail_t = t2 - (head_t - t1);
535 }
536 head_sum[1] = head_t;
537 tail_sum[1] = tail_t;
538 } /* sum = sum+prod */
539 ix += incx;
540 iy += incy;
541 } /* endfor */
542 }
543
544 {
545 double cd[2];
546 cd[0] = (double) alpha_i[0];
547 cd[1] = (double) alpha_i[1];
548 {
549 /* Compute complex-extra = complex-extra * complex-double. */
550 double head_a0, tail_a0;
551 double head_a1, tail_a1;
552 double head_t1, tail_t1;
553 double head_t2, tail_t2;
554 head_a0 = head_sum[0];
555 tail_a0 = tail_sum[0];
556 head_a1 = head_sum[1];
557 tail_a1 = tail_sum[1];
558 /* real part */
559 {
560 /* Compute double-double = double-double * double. */
561 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
562
563 con = head_a0 * split;
564 a11 = con - head_a0;
565 a11 = con - a11;
566 a21 = head_a0 - a11;
567 con = cd[0] * split;
568 b1 = con - cd[0];
569 b1 = con - b1;
570 b2 = cd[0] - b1;
571
572 c11 = head_a0 * cd[0];
573 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
574
575 c2 = tail_a0 * cd[0];
576 t1 = c11 + c2;
577 t2 = (c2 - (t1 - c11)) + c21;
578
579 head_t1 = t1 + t2;
580 tail_t1 = t2 - (head_t1 - t1);
581 }
582 {
583 /* Compute double-double = double-double * double. */
584 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
585
586 con = head_a1 * split;
587 a11 = con - head_a1;
588 a11 = con - a11;
589 a21 = head_a1 - a11;
590 con = cd[1] * split;
591 b1 = con - cd[1];
592 b1 = con - b1;
593 b2 = cd[1] - b1;
594
595 c11 = head_a1 * cd[1];
596 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
597
598 c2 = tail_a1 * cd[1];
599 t1 = c11 + c2;
600 t2 = (c2 - (t1 - c11)) + c21;
601
602 head_t2 = t1 + t2;
603 tail_t2 = t2 - (head_t2 - t1);
604 }
605 head_t2 = -head_t2;
606 tail_t2 = -tail_t2;
607 {
608 /* Compute double-double = double-double + double-double. */
609 double bv;
610 double s1, s2, t1, t2;
611
612 /* Add two hi words. */
613 s1 = head_t1 + head_t2;
614 bv = s1 - head_t1;
615 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
616
617 /* Add two lo words. */
618 t1 = tail_t1 + tail_t2;
619 bv = t1 - tail_t1;
620 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
621
622 s2 += t1;
623
624 /* Renormalize (s1, s2) to (t1, s2) */
625 t1 = s1 + s2;
626 s2 = s2 - (t1 - s1);
627
628 t2 += s2;
629
630 /* Renormalize (t1, t2) */
631 head_t1 = t1 + t2;
632 tail_t1 = t2 - (head_t1 - t1);
633 }
634 head_tmp1[0] = head_t1;
635 tail_tmp1[0] = tail_t1;
636 /* imaginary part */
637 {
638 /* Compute double-double = double-double * double. */
639 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
640
641 con = head_a1 * split;
642 a11 = con - head_a1;
643 a11 = con - a11;
644 a21 = head_a1 - a11;
645 con = cd[0] * split;
646 b1 = con - cd[0];
647 b1 = con - b1;
648 b2 = cd[0] - b1;
649
650 c11 = head_a1 * cd[0];
651 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
652
653 c2 = tail_a1 * cd[0];
654 t1 = c11 + c2;
655 t2 = (c2 - (t1 - c11)) + c21;
656
657 head_t1 = t1 + t2;
658 tail_t1 = t2 - (head_t1 - t1);
659 }
660 {
661 /* Compute double-double = double-double * double. */
662 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
663
664 con = head_a0 * split;
665 a11 = con - head_a0;
666 a11 = con - a11;
667 a21 = head_a0 - a11;
668 con = cd[1] * split;
669 b1 = con - cd[1];
670 b1 = con - b1;
671 b2 = cd[1] - b1;
672
673 c11 = head_a0 * cd[1];
674 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
675
676 c2 = tail_a0 * cd[1];
677 t1 = c11 + c2;
678 t2 = (c2 - (t1 - c11)) + c21;
679
680 head_t2 = t1 + t2;
681 tail_t2 = t2 - (head_t2 - t1);
682 }
683 {
684 /* Compute double-double = double-double + double-double. */
685 double bv;
686 double s1, s2, t1, t2;
687
688 /* Add two hi words. */
689 s1 = head_t1 + head_t2;
690 bv = s1 - head_t1;
691 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
692
693 /* Add two lo words. */
694 t1 = tail_t1 + tail_t2;
695 bv = t1 - tail_t1;
696 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
697
698 s2 += t1;
699
700 /* Renormalize (s1, s2) to (t1, s2) */
701 t1 = s1 + s2;
702 s2 = s2 - (t1 - s1);
703
704 t2 += s2;
705
706 /* Renormalize (t1, t2) */
707 head_t1 = t1 + t2;
708 tail_t1 = t2 - (head_t1 - t1);
709 }
710 head_tmp1[1] = head_t1;
711 tail_tmp1[1] = tail_t1;
712 }
713
714 } /* tmp1 = sum*alpha */
715 {
716 double head_e1, tail_e1;
717 double d1;
718 double d2;
719 /* Real part */
720 d1 = (double) r_v[0] * beta_i[0];
721 d2 = (double) -r_v[1] * beta_i[1];
722 {
723 /* Compute double-double = double + double. */
724 double e, t1, t2;
725
726 /* Knuth trick. */
727 t1 = d1 + d2;
728 e = t1 - d1;
729 t2 = ((d2 - e) + (d1 - (t1 - e)));
730
731 /* The result is t1 + t2, after normalization. */
732 head_e1 = t1 + t2;
733 tail_e1 = t2 - (head_e1 - t1);
734 }
735 head_tmp2[0] = head_e1;
736 tail_tmp2[0] = tail_e1;
737 /* imaginary part */
738 d1 = (double) r_v[0] * beta_i[1];
739 d2 = (double) r_v[1] * beta_i[0];
740 {
741 /* Compute double-double = double + double. */
742 double e, t1, t2;
743
744 /* Knuth trick. */
745 t1 = d1 + d2;
746 e = t1 - d1;
747 t2 = ((d2 - e) + (d1 - (t1 - e)));
748
749 /* The result is t1 + t2, after normalization. */
750 head_e1 = t1 + t2;
751 tail_e1 = t2 - (head_e1 - t1);
752 }
753 head_tmp2[1] = head_e1;
754 tail_tmp2[1] = tail_e1;
755 } /* tmp2 = r*beta */
756 {
757 double head_t, tail_t;
758 double head_a, tail_a;
759 double head_b, tail_b;
760 /* Real part */
761 head_a = head_tmp1[0];
762 tail_a = tail_tmp1[0];
763 head_b = head_tmp2[0];
764 tail_b = tail_tmp2[0];
765 {
766 /* Compute double-double = double-double + double-double. */
767 double bv;
768 double s1, s2, t1, t2;
769
770 /* Add two hi words. */
771 s1 = head_a + head_b;
772 bv = s1 - head_a;
773 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
774
775 /* Add two lo words. */
776 t1 = tail_a + tail_b;
777 bv = t1 - tail_a;
778 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
779
780 s2 += t1;
781
782 /* Renormalize (s1, s2) to (t1, s2) */
783 t1 = s1 + s2;
784 s2 = s2 - (t1 - s1);
785
786 t2 += s2;
787
788 /* Renormalize (t1, t2) */
789 head_t = t1 + t2;
790 tail_t = t2 - (head_t - t1);
791 }
792 head_tmp1[0] = head_t;
793 tail_tmp1[0] = tail_t;
794 /* Imaginary part */
795 head_a = head_tmp1[1];
796 tail_a = tail_tmp1[1];
797 head_b = head_tmp2[1];
798 tail_b = tail_tmp2[1];
799 {
800 /* Compute double-double = double-double + double-double. */
801 double bv;
802 double s1, s2, t1, t2;
803
804 /* Add two hi words. */
805 s1 = head_a + head_b;
806 bv = s1 - head_a;
807 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
808
809 /* Add two lo words. */
810 t1 = tail_a + tail_b;
811 bv = t1 - tail_a;
812 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
813
814 s2 += t1;
815
816 /* Renormalize (s1, s2) to (t1, s2) */
817 t1 = s1 + s2;
818 s2 = s2 - (t1 - s1);
819
820 t2 += s2;
821
822 /* Renormalize (t1, t2) */
823 head_t = t1 + t2;
824 tail_t = t2 - (head_t - t1);
825 }
826 head_tmp1[1] = head_t;
827 tail_tmp1[1] = tail_t;
828 } /* tmp1 = tmp1+tmp2 */
829 ((float *) r)[0] = head_tmp1[0];
830 ((float *) r)[1] = head_tmp1[1]; /* r = tmp1 */
831
832 FPU_FIX_STOP;
833 }
834 break;
835 }
836 }
837