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