ulps: Fix long lines
[mathlib.git] / ulps / ulp_real.c
1 #define _XOPEN_SOURCE 600
2
3 #include <assert.h>
4 #include <float.h>
5 #include <math.h>
6 #include <stdio.h>
7 #include <string.h>
8
9 #include <gmp.h>
10 #include <mpfr.h>
11
12 #include "gen.h"
13 #include "subr_random.h"
14 #include "ulp.h"
15
16 double
17 calculp_double(double computed, double exact)
18 {
19         double xbefore, xafter;
20         double ulps;
21
22         xbefore = nextafter(exact, -INFINITY);
23         xafter = nextafter(exact, +INFINITY);
24         assert(xafter != xbefore);
25
26         ulps = fabs( (exact - computed) / (xafter - xbefore) );
27
28         return (ulps);
29 }
30
31 long double
32 calculp_long_double(long double computedl, long double exactl)
33 {
34         long double xbeforel, xafterl;
35         long double ulpsl;
36
37         xbeforel = nextafterl(exactl, -INFINITY);
38         xafterl = nextafterl(exactl, +INFINITY);
39         assert(xafterl != xbeforel);
40
41         ulpsl = fabsl( (exactl - computedl) / (xafterl - xbeforel) );
42
43         return (ulpsl);
44 }
45
46
47
48 static void
49 populate_vars(const struct fentry *f,
50     double *x, double *y, long double *xl, long double *yl,
51     mpfr_ptr mp_x, mpfr_ptr mp_y, mpfr_ptr mp_xl, mpfr_ptr mp_yl)
52 {
53         const mpfr_rnd_t tonearest = GMP_RNDN;
54         long double txl, tyl;
55         double tx, ty;
56
57         assert(f && x && y && xl && yl);
58         txl = tx = 0.0;
59         tyl = ty = 0.0;
60
61         /* Generate random arguments */
62         if (f->f_narg == 1) {
63                 do {
64                         tx  = random_double(FP_NORMAL);
65                         txl = random_long_double(FP_NORMAL);
66                 } while (!f->f_u.fp1(tx) || !f->f_u.fp1(txl));
67         }
68         if (f->f_narg == 2) {
69                 do {
70                         tx  = random_double(FP_NORMAL);
71                         ty  = random_double(FP_NORMAL);
72                         txl = random_long_double(FP_NORMAL);
73                         tyl = random_long_double(FP_NORMAL);
74                 } while (!f->f_u.fp2(tx, ty) || !f->f_u.fp2(txl, tyl));
75         }
76
77         /* Hack, yikes */
78         if (!strcmp(f->f_name, "yn")) {
79                 txl = tx = 1.0;
80         }
81
82         *x  = tx;
83         *y  = ty;
84         *xl = txl;
85         *yl = tyl;
86
87         /* Copy arguments to mpfr variables */
88         mpfr_set_d (mp_x,  tx,  tonearest);
89         mpfr_set_d (mp_y,  ty,  tonearest);
90         mpfr_set_ld(mp_xl, txl, tonearest);
91         mpfr_set_ld(mp_yl, tyl, tonearest);
92 }
93
94 void
95 getfunctionulp_real(const struct fentry *f, struct ulp *u)
96 {
97         const mpfr_rnd_t tonearest = GMP_RNDN;
98         mpfr_t mp_exactl, mp_xl, mp_yl;
99         mpfr_t mp_exact, mp_x, mp_y;
100         long double xl, yl, computedl, exactl, myulpl;
101         double x, y, computed, exact, myulp;
102         size_t i;
103
104         /* Initialize high precision variables */
105         mpfr_inits2(200, mp_exact,  mp_x,  mp_y,  (mpfr_ptr)NULL);
106         mpfr_inits2(200, mp_exactl, mp_xl, mp_yl, (mpfr_ptr)NULL);
107
108         ULP_INIT(u);
109
110         for (i = 0; i < NITERATIONS; i++) {
111                 if (i % 100 == 0) {
112                         printf("Percentage complete: %2.2f\r",
113                             (100.0 * i)/NITERATIONS);
114                         fflush(stdout);
115                 }
116
117                 /* Generate random arguments */
118                 populate_vars(f, &x, &y, &xl, &yl, mp_x, mp_y, mp_xl, mp_yl);
119
120                 /*
121                  * Ready to call the mpfr*()
122                  * The double version ALWAYS exist!
123                  */
124                 if(f->f_narg == 1) {
125                         f->f_mpfr(mp_exact, mp_x, tonearest);
126                 } else {
127                         /* Hack, yikes */
128                         if (!strcmp(f->f_name, "yn"))
129                                 f->f_mpfr(mp_exact, (long)x, mp_y, tonearest);
130                         else
131                                 f->f_mpfr(mp_exact, mp_x, mp_y, tonearest);
132                 }
133                 exact = mpfr_get_d(mp_exact,  tonearest);
134
135                 /* We can't tell the same for long double functions though */
136                 if (f->f_libml_real) {
137                         if(f->f_narg == 1) {
138                                 f->f_mpfr(mp_exactl, mp_xl, tonearest);
139                         } else {
140                                 f->f_mpfr(mp_exactl, mp_xl, mp_yl, tonearest);
141                         }
142                         exactl = mpfr_get_ld(mp_exactl, tonearest);
143                 }
144
145                 /* Ready to call the libm*() */
146                 if (f->f_narg == 1) {
147                         computed = f->f_libm_real(x);
148                 } else {
149                         computed = f->f_libm_real(x, y);
150                 }
151
152                 if (f->f_libml_real) {
153                         if (f->f_narg == 1) {
154                                 computedl = f->f_libml_real(xl);
155                         } else {
156                                 computedl = f->f_libml_real(xl, yl);
157                         }
158                 }
159
160                 /* Skip bogus results */
161                 if (fpclassify(computed) == FP_NORMAL &&
162                     fpclassify(exact) == FP_NORMAL) {
163                         myulp = calculp_double(computed, exact);
164                         ULP_UPDATE(u, myulp);
165                 } else {
166                         u->ulp_skipped++;
167                 }
168
169                 if (f->f_libml_real) {
170                         if (fpclassify(computedl) == FP_NORMAL &&
171                             fpclassify(exactl) == FP_NORMAL) {
172                                 myulpl = calculp_long_double(computedl, exactl);
173                                 ULP_UPDATEL(u, myulpl);
174                         }
175                 } else {
176                         u->ulp_skippedl++;
177                 }
178         }
179
180         u->ulp_avg  /= (i - u->ulp_skipped);
181         u->ulp_avgl /= (i - u->ulp_skippedl);
182
183         /* Free resources */
184         mpfr_clears(mp_exact,  mp_x,  mp_y,  (mpfr_ptr)NULL);
185         mpfr_clears(mp_exactl, mp_xl, mp_yl, (mpfr_ptr)NULL);
186 }
187
188 void
189 printulps_double(struct ulp u)
190 {
191         if (u.ulp_max > 9.9 || u.ulp_min > 9.9) {
192                 printf("%-10.4e %-10.e %-10.4e   ",
193                     u.ulp_max, u.ulp_min, u.ulp_avg);
194         } else {
195                 printf("%-10.4f %-10.4f %-10.4f   ",
196                     u.ulp_max, u.ulp_min, u.ulp_avg);
197         }
198         printf("%5u\n", u.ulp_skipped);
199 }
200
201 void
202 printulps_long_double(struct ulp u)
203 {
204         if (u.ulp_maxl > 9.9 || u.ulp_minl > 9.9) {
205                 printf("%-10.4e %-10.4e %-10.4e   ",
206                     (double)u.ulp_maxl,
207                     (double)u.ulp_minl,
208                     (double)u.ulp_avgl);
209         } else {
210                 printf("%-10.4f %-10.4f %-10.4f   ",
211                     (double)u.ulp_maxl,
212                     (double)u.ulp_minl,
213                     (double)u.ulp_avgl);
214         }
215         printf("%5u\n", u.ulp_skippedl);
216 }