1 #include <math.h>
2 #include "common.h"
3 #ifdef FUNCTION_PROFILE
4 #include "functable.h"
5 #endif
6
NAME(FLOAT * DA,FLOAT * DB,FLOAT * C,FLOAT * S)7 void NAME(FLOAT *DA, FLOAT *DB, FLOAT *C, FLOAT *S){
8
9 PRINT_DEBUG_NAME;
10
11 IDEBUG_START;
12
13 FUNCTION_PROFILE_START();
14
15 #if defined(__i386__) || defined(__x86_64__) || defined(__ia64__)
16
17 long double da_r = *(DA + 0);
18 long double da_i = *(DA + 1);
19 long double db_r = *(DB + 0);
20 long double db_i = *(DB + 1);
21 long double r;
22
23 long double ada = fabs(da_r) + fabs(da_i);
24
25 if (ada == ZERO) {
26 *C = ZERO;
27 *(S + 0) = ONE;
28 *(S + 1) = ZERO;
29 *(DA + 0) = db_r;
30 *(DA + 1) = db_i;
31 } else {
32 long double alpha_r, alpha_i;
33
34 ada = sqrt(da_r * da_r + da_i * da_i);
35
36 r = sqrt(da_r * da_r + da_i * da_i + db_r * db_r + db_i * db_i);
37
38 alpha_r = da_r / ada;
39 alpha_i = da_i / ada;
40
41 *(C + 0) = ada / r;
42 *(S + 0) = (alpha_r * db_r + alpha_i *db_i) / r;
43 *(S + 1) = (alpha_i * db_r - alpha_r *db_i) / r;
44 *(DA + 0) = alpha_r * r;
45 *(DA + 1) = alpha_i * r;
46 }
47 #else
48 FLOAT da_r = *(DA + 0);
49 FLOAT da_i = *(DA + 1);
50 FLOAT db_r = *(DB + 0);
51 FLOAT db_i = *(DB + 1);
52 FLOAT r;
53
54 FLOAT ada = fabs(da_r) + fabs(da_i);
55 FLOAT adb;
56
57 if (ada == ZERO) {
58 *C = ZERO;
59 *(S + 0) = ONE;
60 *(S + 1) = ZERO;
61 *(DA + 0) = db_r;
62 *(DA + 1) = db_i;
63 } else {
64 FLOAT scale;
65 FLOAT aa_r, aa_i, bb_r, bb_i;
66 FLOAT alpha_r, alpha_i;
67
68 aa_r = fabs(da_r);
69 aa_i = fabs(da_i);
70
71 if (aa_i > aa_r) {
72 aa_r = fabs(da_i);
73 aa_i = fabs(da_r);
74 }
75
76 scale = (aa_i / aa_r);
77 ada = aa_r * sqrt(ONE + scale * scale);
78
79 bb_r = fabs(db_r);
80 bb_i = fabs(db_i);
81
82 if (bb_i > bb_r) {
83 bb_r = fabs(bb_i);
84 bb_i = fabs(bb_r);
85 }
86
87 scale = (bb_i / bb_r);
88 adb = bb_r * sqrt(ONE + scale * scale);
89
90 scale = ada + adb;
91
92 aa_r = da_r / scale;
93 aa_i = da_i / scale;
94 bb_r = db_r / scale;
95 bb_i = db_i / scale;
96
97 r = scale * sqrt(aa_r * aa_r + aa_i * aa_i + bb_r * bb_r + bb_i * bb_i);
98
99 alpha_r = da_r / ada;
100 alpha_i = da_i / ada;
101
102 *(C + 0) = ada / r;
103 *(S + 0) = (alpha_r * db_r + alpha_i *db_i) / r;
104 *(S + 1) = (alpha_i * db_r - alpha_r *db_i) / r;
105 *(DA + 0) = alpha_r * r;
106 *(DA + 1) = alpha_i * r;
107 }
108 #endif
109
110 FUNCTION_PROFILE_END(4, 4, 4);
111
112 IDEBUG_END;
113
114 return;
115 }
116