ulps: Fix long lines
[mathlib.git] / ulps / ulp_complex.c
CommitLineData
3ec3219c
SK
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>
8b95f70b 12#include <mpfr.h>
3ec3219c
SK
13
14#include "gen.h"
f46d092c 15#include "mytypes.h"
4dfa4f2a 16#include "subr_mpc.h"
3ec3219c
SK
17#include "subr_random.h"
18#include "ulp.h"
19
f46d092c 20static double complex
3ec3219c
SK
21calculp_double_complex(double complex computed, double complex exact)
22{
f46d092c 23 double_complex z;
3ec3219c
SK
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
f46d092c 30 z.parts[0] = ulp_real;
31 z.parts[1] = ulp_imag;
32
33 return (z.z);
3ec3219c
SK
34}
35
f46d092c 36static long double complex
4e3bb96b
SK
37calculp_long_double_complex(long double complex computedl,
38 long double complex exactl)
3ec3219c 39{
f46d092c 40 long_double_complex z;
7a97944a 41 long double ulp_reall;
42 long double ulp_imagl;
3ec3219c 43
7a97944a 44 ulp_reall = calculp_long_double(creall(computedl), creall(exactl));
45 ulp_imagl = calculp_long_double(cimagl(computedl), cimagl(exactl));
3ec3219c 46
f46d092c 47 z.parts[0] = ulp_reall;
48 z.parts[1] = ulp_imagl;
49
50 return (z.z);
3ec3219c
SK
51}
52
7a22afa6 53static void
3ec3219c
SK
54populate_complex_vars(const struct fentry *f,
55 double complex *x, double complex *y,
f91bd554 56 long double complex *xl, long double complex *yl,
3ec3219c
SK
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
da06af21 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
3ec3219c
SK
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
7a97944a 100static int
101isbogus(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
ca3e8900 120void
dedfed29 121getfunctionulp_complex(const struct fentry *f, struct ulp_complex *uc)
3ec3219c 122{
f91bd554 123 const mpc_rnd_t tonearest = MPC_RNDNN;
8b95f70b 124 const mpfr_rnd_t mpfr_tonearest = MPFR_RNDN;
3ec3219c
SK
125 mpc_t mp_exactl, mp_xl, mp_yl;
126 mpc_t mp_exact, mp_x, mp_y;
f91bd554
SK
127 long double complex xl, yl, computedl, exactl, myulpl;
128 double complex x, y, computed, exact, myulp;
3ec3219c
SK
129 size_t i;
130
3ec3219c
SK
131 /* Initialize high precision variables */
132 mpc_init2(mp_exact, 200);
7a97944a 133 mpc_init2(mp_x, 200);
134 mpc_init2(mp_y, 200);
3ec3219c
SK
135 mpc_init2(mp_exactl,200);
136 mpc_init2(mp_xl, 200);
137 mpc_init2(mp_yl, 200);
138
f91bd554 139 ULP_COMPLEX_INIT(uc);
3ec3219c
SK
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 */
51b30d33 149 DPRINTF(("populate_complex_vars: entering\n"));
f91bd554 150 populate_complex_vars(f, &x, &y, &xl, &yl, mp_x, mp_y, mp_xl, mp_yl);
51b30d33 151 DPRINTF(("populate_complex_vars: left\n"));
3ec3219c
SK
152
153 /*
da06af21 154 * Ready to call the mpc*()
3ec3219c
SK
155 * The double version ALWAYS exist!
156 */
0da540fe 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));
3ec3219c 160 if(f->f_narg == 1) {
f91bd554 161 f->f_mpc(mp_exact, mp_x, tonearest);
3ec3219c 162 } else {
f91bd554 163 f->f_mpc(mp_exact, mp_x, mp_y, tonearest);
3ec3219c 164 }
8a710f61 165 exact = mpc_get_dc(mp_exact, mpfr_tonearest);
da06af21 166 DPRINTF(("f->f_mpc: left\n"));
167 mpc_free_str(str_x);
168 mpc_free_str(str_y);
3ec3219c
SK
169
170 /* We can't tell the same for long double functions though */
da06af21 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));
24d04885 174 if (f->f_libml_complex) {
3ec3219c 175 if(f->f_narg == 1) {
f91bd554 176 f->f_mpc(mp_exactl, mp_xl, tonearest);
3ec3219c 177 } else {
f91bd554 178 f->f_mpc(mp_exactl, mp_xl, mp_yl, tonearest);
3ec3219c 179 }
8a710f61 180 exactl = mpc_get_ldc(mp_exactl, mpfr_tonearest);
3ec3219c 181 }
da06af21 182 DPRINTF(("f->f_mpcl: left\n"));
183 mpc_free_str(str_xl);
184 mpc_free_str(str_yl);
3ec3219c
SK
185
186 /* Ready to call the libm*() */
187 if (f->f_narg == 1) {
24d04885 188 computed = f->f_libm_complex(x);
3ec3219c 189 } else {
24d04885 190 computed = f->f_libm_complex(x, y);
3ec3219c
SK
191 }
192
24d04885 193 if (f->f_libml_complex) {
3ec3219c 194 if (f->f_narg == 1) {
24d04885 195 computedl = f->f_libml_complex(xl);
3ec3219c 196 } else {
24d04885 197 computedl = f->f_libml_complex(xl, yl);
3ec3219c
SK
198 }
199 }
200
201 /* Skip bogus results */
7a97944a 202 if (!isbogus(computed) && !isbogus(exact)) {
f91bd554
SK
203 myulp = calculp_double_complex(computed, exact);
204 ULP_COMPLEX_UPDATE(uc, myulp);
3ec3219c 205 } else {
4e3bb96b 206 uc->ulp_total.ulp_skipped++;
3ec3219c
SK
207 }
208
24d04885 209 if (f->f_libml_complex) {
7a97944a 210 if (!isbogus(computed) && !isbogus(exact)) {
f91bd554
SK
211 myulpl = calculp_long_double_complex(computedl, exactl);
212 ULP_COMPLEX_UPDATEL(uc, myulpl);
3ec3219c
SK
213 }
214 } else {
4e3bb96b 215 uc->ulp_total.ulp_skippedl++;
3ec3219c
SK
216 }
217 }
218
4e3bb96b 219 uc->ulp_total.ulp_avg /= (i - uc->ulp_total.ulp_skipped);
3ec3219c
SK
220
221 /* Free resources */
f91bd554
SK
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);
3ec3219c
SK
228}
229
230void
f91bd554 231printulps_double_complex(struct ulp_complex uc)
3ec3219c 232{
4e3bb96b 233 if (uc.ulp_total.ulp_max > 9.9 || uc.ulp_total.ulp_min > 9.9) {
ca3e8900 234 printf("%-10.4e %-10.4e %-10.4e ",
4e3bb96b
SK
235 uc.ulp_total.ulp_max,
236 uc.ulp_total.ulp_min,
237 uc.ulp_total.ulp_avg);
238 } else {
7a97944a 239 printf("%-10.4f %-10.4f %-10.4f ",
240 uc.ulp_total.ulp_max,
4e3bb96b
SK
241 uc.ulp_total.ulp_min,
242 uc.ulp_total.ulp_avg);
243 }
244 printf("%5u\n", uc.ulp_total.ulp_skipped);
3ec3219c
SK
245}
246
247void
f91bd554 248printulps_long_double_complex(struct ulp_complex uc)
3ec3219c 249{
4e3bb96b 250 if (uc.ulp_total.ulp_maxl > 9.9 || uc.ulp_total.ulp_minl > 9.9) {
ca3e8900 251 printf("%-10.4e %-10.4e %-10.4e ",
4e3bb96b
SK
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);
3ec3219c 258}