@@ -44,140 +44,145 @@ void CNAME(void *VDA, void *VDB, FLOAT *C, void *VS) {
4444
4545 FUNCTION_PROFILE_START ();
4646
47- if (db_r == ZERO && db_i == ZERO ) {
48- * C = ONE ;
49- * (S + 0 ) = ZERO ;
50- * (S + 1 ) = ZERO ;
51- return ;
52- }
47+ do {
48+ if (db_r == ZERO && db_i == ZERO ) {
49+ * C = ONE ;
50+ * (S + 0 ) = ZERO ;
51+ * (S + 1 ) = ZERO ;
52+ break ;
53+ }
5354
54- long double safmax = 1. /safmin ;
55- #if defined DOUBLE
56- long double rtmax = safmax /DBL_EPSILON ;
57- #else
58- long double rtmax = safmax /FLT_EPSILON ;
59- #endif
60- * (S1 + 0 ) = * (DB + 0 );
61- * (S1 + 1 ) = * (DB + 1 ) * -1 ;
62- if (da_r == ZERO && da_i == ZERO ) {
63- * C = ZERO ;
64- if (db_r == ZERO ) {
65- (* DA ) = fabsl (db_i );
66- * S = * S1 /(* DA );
67- * (S + 1 ) = * (S1 + 1 ) /(* DA );
68- return ;
69- } else if ( db_i == ZERO ) {
70- * DA = fabsl (db_r );
71- * S = * S1 /(* DA );
72- * (S + 1 ) = * (S1 + 1 ) /(* DA );
73- return ;
74- } else {
75- long double g1 = MAX ( fabsl (db_r ), fabsl (db_i ));
76- rtmax = sqrt (safmax /2. );
77- if (g1 > rtmin && g1 < rtmax ) { // unscaled
78- d = sqrt (adb );
79- * S = * S1 /d ;
80- * (S + 1 ) = * (S1 + 1 ) /d ;
81- * DA = d ;
82- * (DA + 1 ) = ZERO ;
83- return ;
84- } else { // scaled algorithm
85- long double u = MIN ( safmax , MAX ( safmin , g1 ));
86- FLOAT gs_r = db_r /u ;
87- FLOAT gs_i = db_i /u ;
88- d = sqrt ( gs_r * gs_r + gs_i * gs_i );
89- * S = gs_r / d ;
90- * (S + 1 ) = (gs_i * -1 ) / d ;
91- * DA = d * u ;
92- * (DA + 1 ) = ZERO ;
93- return ;
94- }
95- }
96- } else {
97- FLOAT f1 = MAX ( fabsl (da_r ), fabsl (da_i ));
98- FLOAT g1 = MAX ( fabsl (db_r ), fabsl (db_i ));
99- rtmax = sqrt (safmax / 4. );
100- if ( f1 > rtmin && f1 < rtmax && g1 > rtmin && g1 < rtmax ) { //unscaled
101- long double h = ada + adb ;
102- double adahsq = sqrt (ada * h );
103- if (ada >= h * safmin ) {
104- * C = sqrt (ada /h );
105- * R = * DA / * C ;
106- * (R + 1 ) = * (DA + 1 ) / * C ;
107- rtmax *= 2. ;
108- if ( ada > rtmin && h < rtmax ) { // no risk of intermediate overflow
109- * S = * S1 * (* DA / adahsq ) - * (S1 + 1 )* (* (DA + 1 )/adahsq );
110- * (S + 1 ) = * S1 * (* (DA + 1 ) / adahsq ) + * (S1 + 1 ) * (* DA /adahsq );
111- } else {
112- * S = * S1 * (* R /h ) - * (S1 + 1 ) * (* (R + 1 )/h );
113- * (S + 1 ) = * S1 * (* (R + 1 )/h ) + * (S1 + 1 ) * (* (R )/h );
55+ long double safmax = 1. /safmin ;
56+ #if defined DOUBLE
57+ long double rtmax = safmax /DBL_EPSILON ;
58+ #else
59+ long double rtmax = safmax /FLT_EPSILON ;
60+ #endif
61+ * (S1 + 0 ) = * (DB + 0 );
62+ * (S1 + 1 ) = * (DB + 1 ) * -1 ;
63+ if (da_r == ZERO && da_i == ZERO ) {
64+ * C = ZERO ;
65+ if (db_r == ZERO ) {
66+ (* DA ) = fabsl (db_i );
67+ * S = * S1 /(* DA );
68+ * (S + 1 ) = * (S1 + 1 ) /(* DA );
69+ break ;
70+ } else if ( db_i == ZERO ) {
71+ * DA = fabsl (db_r );
72+ * S = * S1 /(* DA );
73+ * (S + 1 ) = * (S1 + 1 ) /(* DA );
74+ break ;
75+ } else {
76+ long double g1 = MAX ( fabsl (db_r ), fabsl (db_i ));
77+ rtmax = sqrt (safmax /2. );
78+ if (g1 > rtmin && g1 < rtmax ) { // unscaled
79+ d = sqrt (adb );
80+ * S = * S1 /d ;
81+ * (S + 1 ) = * (S1 + 1 ) /d ;
82+ * DA = d ;
83+ * (DA + 1 ) = ZERO ;
84+ break ;
85+ } else { // scaled algorithm
86+ long double u = MIN ( safmax , MAX ( safmin , g1 ));
87+ FLOAT gs_r = db_r /u ;
88+ FLOAT gs_i = db_i /u ;
89+ d = sqrt ( gs_r * gs_r + gs_i * gs_i );
90+ * S = gs_r / d ;
91+ * (S + 1 ) = (gs_i * -1 ) / d ;
92+ * DA = d * u ;
93+ * (DA + 1 ) = ZERO ;
94+ break ;
95+ }
96+ }
97+ } else {
98+ FLOAT f1 = MAX ( fabsl (da_r ), fabsl (da_i ));
99+ FLOAT g1 = MAX ( fabsl (db_r ), fabsl (db_i ));
100+ rtmax = sqrt (safmax / 4. );
101+ if ( f1 > rtmin && f1 < rtmax && g1 > rtmin && g1 < rtmax ) { //unscaled
102+ long double h = ada + adb ;
103+ double adahsq = sqrt (ada * h );
104+ if (ada >= h * safmin ) {
105+ * C = sqrt (ada /h );
106+ * R = * DA / * C ;
107+ * (R + 1 ) = * (DA + 1 ) / * C ;
108+ rtmax *= 2. ;
109+ if ( ada > rtmin && h < rtmax ) { // no risk of intermediate overflow
110+ * S = * S1 * (* DA / adahsq ) - * (S1 + 1 )* (* (DA + 1 )/adahsq );
111+ * (S + 1 ) = * S1 * (* (DA + 1 ) / adahsq ) + * (S1 + 1 ) * (* DA /adahsq );
112+ } else {
113+ * S = * S1 * (* R /h ) - * (S1 + 1 ) * (* (R + 1 )/h );
114+ * (S + 1 ) = * S1 * (* (R + 1 )/h ) + * (S1 + 1 ) * (* (R )/h );
115+ }
116+ } else {
117+ * C = ada / adahsq ;
118+ if (* C >= safmin ) {
119+ * R = * DA / * C ;
120+ * (R + 1 ) = * (DA + 1 ) / * C ;
121+ } else {
122+ * R = * DA * (h / adahsq );
123+ * (R + 1 ) = * (DA + 1 ) * (h / adahsq );
124+ }
125+ * S = * S1 * ada / adahsq ;
126+ * (S + 1 ) = * (S1 + 1 ) * ada / adahsq ;
127+ }
128+ * DA = * R ;
129+ * (DA + 1 )= * (R + 1 );
130+ break ;
131+ } else { // scaled
132+ FLOAT fs_r , fs_i , gs_r , gs_i ;
133+ long double v ,w ,f2 ,g2 ,h ;
134+ long double u = MIN ( safmax , MAX ( safmin , MAX (f1 ,g1 )));
135+ gs_r = db_r /u ;
136+ gs_i = db_i /u ;
137+ g2 = sqrt ( gs_r * gs_r + gs_i * gs_i );
138+ if (f1 /u < rtmin ) {
139+ v = MIN (safmax , MAX (safmin , f1 ));
140+ w = v / u ;
141+ fs_r = * DA / v ;
142+ fs_i = * (DA + 1 ) / v ;
143+ f2 = sqrt ( fs_r * fs_r + fs_i * fs_i );
144+ h = f2 * w * w + g2 ;
145+ } else { // use same scaling for both
146+ w = 1. ;
147+ fs_r = * DA / u ;
148+ fs_i = * (DA + 1 ) / u ;
149+ f2 = sqrt ( fs_r * fs_r + fs_i * fs_i );
150+ h = f2 + g2 ;
151+ }
152+ if ( f2 >= h * safmin ) {
153+ * C = sqrt ( f2 / h );
154+ * DA = fs_r / * C ;
155+ * (DA + 1 ) = fs_i / * C ;
156+ rtmax *= 2 ;
157+ if ( f2 > rtmin && h < rtmax ) {
158+ * S = gs_r * (fs_r /sqrt (f2 * h )) - gs_i * (fs_i / sqrt (f2 * h ));
159+ * (S + 1 ) = gs_r * (fs_i /sqrt (f2 * h )) + gs_i * -1. * (fs_r / sqrt (f2 * h ));
160+ } else {
161+ * S = gs_r * (* DA /h ) - gs_i * (* (DA + 1 ) / h );
162+ * (S + 1 ) = gs_r * (* (DA + 1 ) /h ) + gs_i * -1. * (* DA / h );
163+ }
164+ } else { // intermediates might overflow
165+ d = sqrt ( f2 * h );
166+ * C = f2 /d ;
167+ if (* C >= safmin ) {
168+ * DA = fs_r / * C ;
169+ * (DA + 1 ) = fs_i / * C ;
170+ } else {
171+ * DA = fs_r * (h / d );
172+ * (DA + 1 ) = fs_i / (h / d );
173+ }
174+ * S = gs_r * (fs_r /d ) - gs_i * (fs_i / d );
175+ * (S + 1 ) = gs_r * (fs_i /d ) + gs_i * -1. * (fs_r / d );
176+ }
177+ * C *= w ;
178+ * DA *= u ;
179+ * (DA + 1 ) *= u ;
180+ break ;
114181 }
115- } else {
116- * C = ada / adahsq ;
117- if (* C >= safmin ) {
118- * R = * DA / * C ;
119- * (R + 1 ) = * (DA + 1 ) / * C ;
120- } else {
121- * R = * DA * (h / adahsq );
122- * (R + 1 ) = * (DA + 1 ) * (h / adahsq );
123- }
124- * S = * S1 * ada / adahsq ;
125- * (S + 1 ) = * (S1 + 1 ) * ada / adahsq ;
126- }
127- * DA = * R ;
128- * (DA + 1 )= * (R + 1 );
129- return ;
130- } else { // scaled
131- FLOAT fs_r , fs_i , gs_r , gs_i ;
132- long double v ,w ,f2 ,g2 ,h ;
133- long double u = MIN ( safmax , MAX ( safmin , MAX (f1 ,g1 )));
134- gs_r = db_r /u ;
135- gs_i = db_i /u ;
136- g2 = sqrt ( gs_r * gs_r + gs_i * gs_i );
137- if (f1 /u < rtmin ) {
138- v = MIN (safmax , MAX (safmin , f1 ));
139- w = v / u ;
140- fs_r = * DA / v ;
141- fs_i = * (DA + 1 ) / v ;
142- f2 = sqrt ( fs_r * fs_r + fs_i * fs_i );
143- h = f2 * w * w + g2 ;
144- } else { // use same scaling for both
145- w = 1. ;
146- fs_r = * DA / u ;
147- fs_i = * (DA + 1 ) / u ;
148- f2 = sqrt ( fs_r * fs_r + fs_i * fs_i );
149- h = f2 + g2 ;
150- }
151- if ( f2 >= h * safmin ) {
152- * C = sqrt ( f2 / h );
153- * DA = fs_r / * C ;
154- * (DA + 1 ) = fs_i / * C ;
155- rtmax *= 2 ;
156- if ( f2 > rtmin && h < rtmax ) {
157- * S = gs_r * (fs_r /sqrt (f2 * h )) - gs_i * (fs_i / sqrt (f2 * h ));
158- * (S + 1 ) = gs_r * (fs_i /sqrt (f2 * h )) + gs_i * -1. * (fs_r / sqrt (f2 * h ));
159- } else {
160- * S = gs_r * (* DA /h ) - gs_i * (* (DA + 1 ) / h );
161- * (S + 1 ) = gs_r * (* (DA + 1 ) /h ) + gs_i * -1. * (* DA / h );
162- }
163- } else { // intermediates might overflow
164- d = sqrt ( f2 * h );
165- * C = f2 /d ;
166- if (* C >= safmin ) {
167- * DA = fs_r / * C ;
168- * (DA + 1 ) = fs_i / * C ;
169- } else {
170- * DA = fs_r * (h / d );
171- * (DA + 1 ) = fs_i / (h / d );
172- }
173- * S = gs_r * (fs_r /d ) - gs_i * (fs_i / d );
174- * (S + 1 ) = gs_r * (fs_i /d ) + gs_i * -1. * (fs_r / d );
175- }
176- * C *= w ;
177- * DA *= u ;
178- * (DA + 1 ) *= u ;
179- return ;
180182 }
181- }
183+ } while (0 );
184+
185+ FUNCTION_PROFILE_END (4 , 4 , 4 );
186+ IDEBUG_END ;
182187}
183188
0 commit comments