ulps: Fix long lines
[mathlib.git] / ulps / ulp_real.c
CommitLineData
32adbb0c
SK
1#define _XOPEN_SOURCE 600
2
a95d68a2
SK
3#include <assert.h>
4#include <float.h>
5#include <math.h>
6#include <stdio.h>
8402c1f4 7#include <string.h>
a95d68a2
SK
8
9#include <gmp.h>
10#include <mpfr.h>
11
12#include "gen.h"
13#include "subr_random.h"
14#include "ulp.h"
15
3ec3219c 16double
8727eb1a 17calculp_double(double computed, double exact)
a95d68a2
SK
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
3ec3219c 31long double
a9068f41 32calculp_long_double(long double computedl, long double exactl)
8727eb1a 33{
a9068f41 34 long double xbeforel, xafterl;
35 long double ulpsl;
8727eb1a 36
a9068f41 37 xbeforel = nextafterl(exactl, -INFINITY);
38 xafterl = nextafterl(exactl, +INFINITY);
39 assert(xafterl != xbeforel);
8727eb1a 40
a9068f41 41 ulpsl = fabsl( (exactl - computedl) / (xafterl - xbeforel) );
8727eb1a 42
a9068f41 43 return (ulpsl);
8727eb1a 44}
45
5226af1e
SK
46
47
4a43ae5a 48static void
49populate_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{
e03f9bdc 53 const mpfr_rnd_t tonearest = GMP_RNDN;
4a43ae5a 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);
b28edd96 66 } while (!f->f_u.fp1(tx) || !f->f_u.fp1(txl));
4a43ae5a 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);
b28edd96 74 } while (!f->f_u.fp2(tx, ty) || !f->f_u.fp2(txl, tyl));
4a43ae5a 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);
e03f9bdc 89 mpfr_set_d (mp_y, ty, tonearest);
4a43ae5a 90 mpfr_set_ld(mp_xl, txl, tonearest);
91 mpfr_set_ld(mp_yl, tyl, tonearest);
92}
93
ca3e8900 94void
4565678c 95getfunctionulp_real(const struct fentry *f, struct ulp *u)
a95d68a2 96{
b942df82 97 const mpfr_rnd_t tonearest = GMP_RNDN;
94545ea5 98 mpfr_t mp_exactl, mp_xl, mp_yl;
99 mpfr_t mp_exact, mp_x, mp_y;
a9068f41 100 long double xl, yl, computedl, exactl, myulpl;
101 double x, y, computed, exact, myulp;
a95d68a2 102 size_t i;
b942df82 103
a95d68a2 104 /* Initialize high precision variables */
05ff8068 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);
a95d68a2 107
8727eb1a 108 ULP_INIT(u);
a95d68a2
SK
109
110 for (i = 0; i < NITERATIONS; i++) {
a69590ef 111 if (i % 100 == 0) {
a9068f41 112 printf("Percentage complete: %2.2f\r",
113 (100.0 * i)/NITERATIONS);
a69590ef 114 fflush(stdout);
115 }
0d3e05c6 116
a95d68a2 117 /* Generate random arguments */
4a43ae5a 118 populate_vars(f, &x, &y, &xl, &yl, mp_x, mp_y, mp_xl, mp_yl);
a95d68a2 119
a9068f41 120 /*
121 * Ready to call the mpfr*()
94545ea5 122 * The double version ALWAYS exist!
a9068f41 123 */
a95d68a2 124 if(f->f_narg == 1) {
94545ea5 125 f->f_mpfr(mp_exact, mp_x, tonearest);
a95d68a2 126 } else {
05ff8068 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);
a95d68a2 132 }
94545ea5 133 exact = mpfr_get_d(mp_exact, tonearest);
a95d68a2 134
a9068f41 135 /* We can't tell the same for long double functions though */
24d04885 136 if (f->f_libml_real) {
a9068f41 137 if(f->f_narg == 1) {
94545ea5 138 f->f_mpfr(mp_exactl, mp_xl, tonearest);
a9068f41 139 } else {
94545ea5 140 f->f_mpfr(mp_exactl, mp_xl, mp_yl, tonearest);
a9068f41 141 }
94545ea5 142 exactl = mpfr_get_ld(mp_exactl, tonearest);
a9068f41 143 }
144
a95d68a2
SK
145 /* Ready to call the libm*() */
146 if (f->f_narg == 1) {
24d04885 147 computed = f->f_libm_real(x);
a95d68a2 148 } else {
24d04885 149 computed = f->f_libm_real(x, y);
a95d68a2
SK
150 }
151
24d04885 152 if (f->f_libml_real) {
a9068f41 153 if (f->f_narg == 1) {
24d04885 154 computedl = f->f_libml_real(xl);
a9068f41 155 } else {
24d04885 156 computedl = f->f_libml_real(xl, yl);
a9068f41 157 }
a95d68a2
SK
158 }
159
a9068f41 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);
94545ea5 165 } else {
166 u->ulp_skipped++;
a9068f41 167 }
94545ea5 168
24d04885 169 if (f->f_libml_real) {
a9068f41 170 if (fpclassify(computedl) == FP_NORMAL &&
171 fpclassify(exactl) == FP_NORMAL) {
172 myulpl = calculp_long_double(computedl, exactl);
fe583431 173 ULP_UPDATEL(u, myulpl);
a9068f41 174 }
94545ea5 175 } else {
176 u->ulp_skippedl++;
a9068f41 177 }
a95d68a2 178 }
a95d68a2 179
94545ea5 180 u->ulp_avg /= (i - u->ulp_skipped);
181 u->ulp_avgl /= (i - u->ulp_skippedl);
182
a95d68a2 183 /* Free resources */
4a43ae5a 184 mpfr_clears(mp_exact, mp_x, mp_y, (mpfr_ptr)NULL);
185 mpfr_clears(mp_exactl, mp_xl, mp_yl, (mpfr_ptr)NULL);
a95d68a2 186}
2e1d545e 187
188void
d1e1ecce 189printulps_double(struct ulp u)
2e1d545e 190{
05ff8068 191 if (u.ulp_max > 9.9 || u.ulp_min > 9.9) {
192 printf("%-10.4e %-10.e %-10.4e ",
2e1d545e 193 u.ulp_max, u.ulp_min, u.ulp_avg);
194 } else {
e75639fd 195 printf("%-10.4f %-10.4f %-10.4f ",
2e1d545e 196 u.ulp_max, u.ulp_min, u.ulp_avg);
197 }
31a95504 198 printf("%5u\n", u.ulp_skipped);
2e1d545e 199}
200
201void
202printulps_long_double(struct ulp u)
203{
05ff8068 204 if (u.ulp_maxl > 9.9 || u.ulp_minl > 9.9) {
e75639fd 205 printf("%-10.4e %-10.4e %-10.4e ",
ca3e8900
SK
206 (double)u.ulp_maxl,
207 (double)u.ulp_minl,
208 (double)u.ulp_avgl);
2e1d545e 209 } else {
e75639fd 210 printf("%-10.4f %-10.4f %-10.4f ",
ca3e8900
SK
211 (double)u.ulp_maxl,
212 (double)u.ulp_minl,
213 (double)u.ulp_avgl);
2e1d545e 214 }
31a95504 215 printf("%5u\n", u.ulp_skippedl);
2e1d545e 216}