ulps: Fix long lines
[mathlib.git] / ulps / ulp_complex.c
1 #define _XOPEN_SOURCE 600
2
3 #include <assert.h>
4 #include <complex.h>
5 #include <float.h>
6 #include <math.h>
7 #include <stdio.h>
8 #include <string.h>
9
10 #include <gmp.h>
11 #include <mpc.h>
12 #include <mpfr.h>
13
14 #include "gen.h"
15 #include "mytypes.h"
16 #include "subr_mpc.h"
17 #include "subr_random.h"
18 #include "ulp.h"
19
20 static double complex
21 calculp_double_complex(double complex computed, double complex exact)
22 {
23         double_complex z;
24         double ulp_real;
25         double ulp_imag;
26
27         ulp_real = calculp_double(creal(computed), creal(exact));
28         ulp_imag = calculp_double(cimag(computed), cimag(exact));
29
30         z.parts[0] = ulp_real;
31         z.parts[1] = ulp_imag;
32
33         return (z.z);
34 }
35
36 static long double complex
37 calculp_long_double_complex(long double complex computedl,
38                             long double complex exactl)
39 {
40         long_double_complex z;
41         long double ulp_reall;
42         long double ulp_imagl;
43
44         ulp_reall = calculp_long_double(creall(computedl), creall(exactl));
45         ulp_imagl = calculp_long_double(cimagl(computedl), cimagl(exactl));
46
47         z.parts[0] = ulp_reall;
48         z.parts[1] = ulp_imagl;
49
50         return (z.z);
51 }
52
53 static void
54 populate_complex_vars(const struct fentry *f,
55     double complex *x, double complex *y,
56     long double complex *xl, long double complex *yl,
57     mpc_t mp_x, mpc_t mp_y, mpc_t mp_xl, mpc_t mp_yl)
58 {
59         const mpc_rnd_t tonearest = MPC_RNDNN;
60         long double complex txl, tyl;
61         double complex tx, ty;
62
63         assert(f && x && y && xl && yl);
64         txl = tx = 0.0;
65         tyl = ty = 0.0;
66
67         /* Generate random arguments */
68         if (f->f_narg == 1) {
69                 do {
70                         tx  = random_double_complex(FP_NORMAL);
71                         txl = random_long_double_complex(FP_NORMAL);
72                 } while (!f->f_uc.fp1(tx) || !f->f_uc.fp1(txl));
73         }
74         if (f->f_narg == 2) {
75                 do {
76                         tx  = random_double_complex(FP_NORMAL);
77                         ty  = random_double_complex(FP_NORMAL);
78                         txl = random_long_double_complex(FP_NORMAL);
79                         tyl = random_long_double_complex(FP_NORMAL);
80                 } while (!f->f_uc.fp2(tx, ty) || !f->f_uc.fp2(txl, tyl));
81         }
82
83         *x  = tx;
84         *y  = ty;
85         *xl = txl;
86         *yl = tyl;
87
88         DPRINTF(( "x  = % .16e + I *% .16e\n",  creal (*x),  cimag (*x)));
89         DPRINTF(( "y  = % .16e + I *% .16e\n",  creal (*y),  cimag (*y)));
90         DPRINTF(("xl = % .16Le + I *% .16Le\n", creall(*xl), cimagl(*xl)));
91         DPRINTF(("yl = % .16Le + I *% .16Le\n", creall(*yl), cimagl(*yl)));
92
93         /* Copy arguments to mpc variables */
94         mpc_set_d_d  (mp_x,  creal (tx),  cimag (tx),  tonearest);
95         mpc_set_d_d  (mp_y,  creal (ty),  cimag (ty),  tonearest);
96         mpc_set_ld_ld(mp_xl, creall(txl), cimagl(txl), tonearest);
97         mpc_set_ld_ld(mp_yl, creall(tyl), cimagl(tyl), tonearest);
98 }
99
100 static int
101 isbogus(long double complex z)
102 {
103         long double real = creall(z);
104         long double imag = cimagl(z);
105         int bogus = 0;
106
107         /* Real part */
108         if ((fpclassify(real) != FP_NORMAL) &&
109             (fpclassify(real) != FP_SUBNORMAL))
110                 bogus = 1;
111
112         /* Imaginary part */
113         if ((fpclassify(imag) != FP_NORMAL) &&
114             (fpclassify(imag) != FP_SUBNORMAL))
115                 bogus = 1;
116
117         return (bogus);
118 }
119
120 void
121 getfunctionulp_complex(const struct fentry *f, struct ulp_complex *uc)
122 {
123         const mpc_rnd_t tonearest = MPC_RNDNN;
124         const mpfr_rnd_t mpfr_tonearest = MPFR_RNDN;
125         mpc_t mp_exactl, mp_xl, mp_yl;
126         mpc_t mp_exact, mp_x, mp_y;
127         long double complex xl, yl, computedl, exactl, myulpl;
128         double complex x, y, computed, exact, myulp;
129         size_t i;
130
131         /* Initialize high precision variables */
132         mpc_init2(mp_exact, 200);
133         mpc_init2(mp_x,     200);
134         mpc_init2(mp_y,     200);
135         mpc_init2(mp_exactl,200);
136         mpc_init2(mp_xl,    200);
137         mpc_init2(mp_yl,    200);
138
139         ULP_COMPLEX_INIT(uc);
140
141         for (i = 0; i < NITERATIONS; i++) {
142                 if (i % 100 == 0) {
143                         printf("Percentage complete: %2.2f\r",
144                             (100.0 * i)/NITERATIONS);
145                         fflush(stdout);
146                 }
147
148                 /* Generate random arguments */
149                 DPRINTF(("populate_complex_vars: entering\n"));
150                 populate_complex_vars(f, &x, &y, &xl, &yl, mp_x, mp_y, mp_xl, mp_yl);
151                 DPRINTF(("populate_complex_vars: left\n"));
152
153                 /*
154                  * Ready to call the mpc*()
155                  * The double version ALWAYS exist!
156                  */
157                 char *str_x = mpc_get_str(10, 16, mp_x, tonearest);
158                 char *str_y = mpc_get_str(10, 16, mp_y, tonearest);
159                 DPRINTF(("f->f_mpc: enter x=%s, y=%s\n", str_x, str_y));
160                 if(f->f_narg == 1) {
161                         f->f_mpc(mp_exact, mp_x, tonearest);
162                 } else {
163                         f->f_mpc(mp_exact, mp_x, mp_y, tonearest);
164                 }
165                 exact = mpc_get_dc(mp_exact, mpfr_tonearest);
166                 DPRINTF(("f->f_mpc: left\n"));
167                 mpc_free_str(str_x);
168                 mpc_free_str(str_y);
169
170                 /* We can't tell the same for long double functions though */
171                 char *str_xl = mpc_get_str(10, 35, mp_xl, tonearest);
172                 char *str_yl = mpc_get_str(10, 35, mp_yl, tonearest);
173                 DPRINTF(("f->f_mpcl: enter x=%s, y=%s\n", str_xl, str_yl));
174                 if (f->f_libml_complex) {
175                         if(f->f_narg == 1) {
176                                 f->f_mpc(mp_exactl, mp_xl, tonearest);
177                         } else {
178                                 f->f_mpc(mp_exactl, mp_xl, mp_yl, tonearest);
179                         }
180                         exactl = mpc_get_ldc(mp_exactl, mpfr_tonearest);
181                 }
182                 DPRINTF(("f->f_mpcl: left\n"));
183                 mpc_free_str(str_xl);
184                 mpc_free_str(str_yl);
185
186                 /* Ready to call the libm*() */
187                 if (f->f_narg == 1) {
188                         computed = f->f_libm_complex(x);
189                 } else {
190                         computed = f->f_libm_complex(x, y);
191                 }
192
193                 if (f->f_libml_complex) {
194                         if (f->f_narg == 1) {
195                                 computedl = f->f_libml_complex(xl);
196                         } else {
197                                 computedl = f->f_libml_complex(xl, yl);
198                         }
199                 }
200
201                 /* Skip bogus results */
202                 if (!isbogus(computed) && !isbogus(exact)) {
203                         myulp = calculp_double_complex(computed, exact);
204                         ULP_COMPLEX_UPDATE(uc, myulp);
205                 } else {
206                         uc->ulp_total.ulp_skipped++;
207                 }
208
209                 if (f->f_libml_complex) {
210                         if (!isbogus(computed) && !isbogus(exact)) {
211                                 myulpl = calculp_long_double_complex(computedl, exactl);
212                                 ULP_COMPLEX_UPDATEL(uc, myulpl);
213                         }
214                 } else {
215                         uc->ulp_total.ulp_skippedl++;
216                 }
217         }
218
219         uc->ulp_total.ulp_avg  /= (i - uc->ulp_total.ulp_skipped);
220
221         /* Free resources */
222         mpc_clear(mp_exact);
223         mpc_clear(mp_x);
224         mpc_clear(mp_y);
225         mpc_clear(mp_exactl);
226         mpc_clear(mp_xl);
227         mpc_clear(mp_yl);
228 }
229
230 void
231 printulps_double_complex(struct ulp_complex uc)
232 {
233         if (uc.ulp_total.ulp_max > 9.9 || uc.ulp_total.ulp_min > 9.9) {
234                 printf("%-10.4e %-10.4e %-10.4e   ",
235                     uc.ulp_total.ulp_max,
236                     uc.ulp_total.ulp_min,
237                     uc.ulp_total.ulp_avg);
238         } else {
239                 printf("%-10.4f %-10.4f %-10.4f   ",
240                     uc.ulp_total.ulp_max,
241                     uc.ulp_total.ulp_min,
242                     uc.ulp_total.ulp_avg);
243         }
244         printf("%5u\n", uc.ulp_total.ulp_skipped);
245 }
246
247 void
248 printulps_long_double_complex(struct ulp_complex uc)
249 {
250         if (uc.ulp_total.ulp_maxl > 9.9 || uc.ulp_total.ulp_minl > 9.9) {
251                 printf("%-10.4e %-10.4e %-10.4e   ",
252                     (double)uc.ulp_real.ulp_maxl,
253                     (double)uc.ulp_real.ulp_minl,
254                     (double)uc.ulp_real.ulp_avgl);
255         } else {
256         }
257         printf("%5u\n", uc.ulp_total.ulp_skippedl);
258 }