selftest_num.c (8329B)
1 /* selftest_num.c - pruebas de Dec y dmath contra long double (64 bits de mantisa, 2 * ~19 digitos: alcanza como oraculo para 15 digitos). */ 3 #include <math.h> 4 #include <stdio.h> 5 #include <stdlib.h> 6 #include <string.h> 7 #include "selftest.h" 8 #include "dec.h" 9 #include "dmath.h" 10 11 static long double to_ld(Dec d) 12 { 13 char b[64]; 14 dec_to_text(d, b, sizeof b); 15 return strtold(b, 0); 16 } 17 18 static Dec from_ld(long double x) 19 { 20 char b[64]; 21 snprintf(b, sizeof b, "%.20Le", x); 22 return dec_parse(b, 0); 23 } 24 25 /* error relativo en unidades del digito 15 */ 26 static double ulp15(Dec got, long double want) 27 { 28 if (got.err) return 1e9; 29 long double g = to_ld(got); 30 if (want == 0) return fabsl(g) < 1e-30L ? 0 : 1e9; 31 long double rel = fabsl((g - want) / want); 32 return (double)(rel * 1e15L); 33 } 34 35 static void eq(const char *expr, Dec got, const char *want) 36 { 37 char b[64]; 38 dec_to_text(dec_round_sig(got, 15), b, sizeof b); 39 st_check(!strcmp(b, want), "%s = %s (esperaba %s)", expr, b, want); 40 } 41 42 typedef struct { const char *name; Dec (*f)(Dec); long double (*g)(long double); double lo, hi; bool logscale; } Fn1; 43 44 static Dec f_exp(Dec x) { return d_exp(x); } 45 static Dec f_ln(Dec x) { return d_ln(x); } 46 static Dec f_log(Dec x) { return d_log10(x); } 47 static Dec f_sinr(Dec x) { return d_sin(x, ANG_RAD); } 48 static Dec f_cosr(Dec x) { return d_cos(x, ANG_RAD); } 49 static Dec f_tanr(Dec x) { return d_tan(x, ANG_RAD); } 50 static Dec f_sind(Dec x) { return d_sin(x, ANG_DEG); } 51 static Dec f_atan(Dec x) { return d_atan(x, ANG_RAD); } 52 static Dec f_asin(Dec x) { return d_asin(x, ANG_RAD); } 53 static Dec f_acos(Dec x) { return d_acos(x, ANG_RAD); } 54 static Dec f_sqrt(Dec x) { return dec_sqrt(x); } 55 static Dec f_cbrt(Dec x) { return d_root(x, 3); } 56 static Dec f_sinh(Dec x) { return d_sinh(x); } 57 static Dec f_cosh(Dec x) { return d_cosh(x); } 58 static Dec f_tanh(Dec x) { return d_tanh(x); } 59 static Dec f_asinh(Dec x) { return d_asinh(x); } 60 static Dec f_acosh(Dec x) { return d_acosh(x); } 61 static Dec f_atanh(Dec x) { return d_atanh(x); } 62 static Dec f_p10(Dec x) { return d_pow10(x); } 63 static long double g_sind(long double x) { return sinl(fmodl(x, 360.0L) * (3.14159265358979323846264338327950288L / 180)); } 64 static long double g_p10(long double x) { return powl(10.0L, x); } 65 66 static const Fn1 FNS[] = { 67 { "exp", f_exp, expl, -200, 200, false }, 68 { "exp chico", f_exp, expl, -1e-6, 1e-6, false }, 69 { "ln", f_ln, logl, 1e-90, 1e90, true }, 70 { "ln ~1", f_ln, logl, 0.9, 1.1, false }, 71 { "log", f_log, log10l, 1e-90, 1e90, true }, 72 { "sin rad", f_sinr, sinl, -100, 100, false }, 73 { "sin rad chico", f_sinr, sinl, -1e-5, 1e-5, false }, 74 { "cos rad", f_cosr, cosl, -100, 100, false }, 75 { "tan rad", f_tanr, tanl, -1.5, 1.5, false }, 76 { "sin deg", f_sind, g_sind, -720, 720, false }, 77 { "atan", f_atan, atanl, -1e6, 1e6, false }, 78 { "atan chico", f_atan, atanl, -2, 2, false }, 79 { "asin", f_asin, asinl, -1, 1, false }, 80 { "acos", f_acos, acosl, -1, 1, false }, 81 { "sqrt", f_sqrt, sqrtl, 1e-90, 1e90, true }, 82 { "cbrt", f_cbrt, cbrtl, 1e-90, 1e90, true }, 83 { "sinh", f_sinh, sinhl, -50, 50, false }, 84 { "sinh chico", f_sinh, sinhl, -1e-3, 1e-3, false }, 85 { "cosh", f_cosh, coshl, -50, 50, false }, 86 { "tanh", f_tanh, tanhl, -5, 5, false }, 87 { "asinh", f_asinh, asinhl, -1e5, 1e5, false }, 88 { "asinh chico", f_asinh, asinhl, -1e-4, 1e-4, false }, 89 { "acosh", f_acosh, acoshl, 1, 1e5, false }, 90 { "atanh", f_atanh, atanhl, -0.999, 0.999, false }, 91 { "10^x", f_p10, g_p10, -90, 90, false }, 92 }; 93 94 void selftest_num(int n) 95 { 96 /* casos fijos */ 97 eq("0.1+0.2", dec_add(dec_parse("0.1", 0), dec_parse("0.2", 0)), "3E-1"); 98 eq("1/3*3", dec_mul(dec_div(DEC_ONE, dec_int(3)), dec_int(3)), "1"); 99 eq("2/3", dec_div(dec_int(2), dec_int(3)), "6.66666666666667E-1"); 100 eq("1e99*10", dec_mul(dec_parse("1e99", 0), dec_int(10)), "Error"); 101 eq("sqrt 2", dec_sqrt(dec_int(2)), "1.4142135623731"); 102 eq("sqrt 16", dec_sqrt(dec_int(16)), "4"); 103 eq("1-0.9", dec_sub(DEC_ONE, dec_parse("0.9", 0)), "1E-1"); 104 eq("123456789*987654321", dec_mul(dec_int(123456789), dec_int(987654321)), "1.21932631112635E17"); 105 eq("sin 30", d_sin(dec_int(30), ANG_DEG), "5E-1"); 106 eq("cos 90", d_cos(dec_int(90), ANG_DEG), "0"); 107 eq("tan 90", d_tan(dec_int(90), ANG_DEG), "Error"); 108 eq("sin pi", d_sin(d_pi(), ANG_RAD), "0"); 109 eq("tan 45", d_tan(dec_int(45), ANG_DEG), "1"); 110 eq("asin 1 deg", d_asin(DEC_ONE, ANG_DEG), "9E1"); 111 eq("acos -1 deg", d_acos(dec_int(-1), ANG_DEG), "1.8E2"); 112 eq("atan 1 deg", d_atan(DEC_ONE, ANG_DEG), "4.5E1"); 113 eq("ln 1", d_ln(DEC_ONE), "0"); 114 /* verificados con bc -l (scale=40), donde long double no alcanza */ 115 eq("log 1.00003612497617417", d_log10(dec_parse("1.00003612497617417", 0)), "1.56885944379845E-5"); 116 eq("sin -9.42467122311921382", d_sin(dec_parse("-9.42467122311921382", 0), ANG_RAD), "-1.0673764996322E-4"); 117 eq("cos -7.85385408338804454", d_cos(dec_parse("-7.85385408338804454", 0), ANG_RAD), "1.275505860927E-4"); 118 eq("atanh -0.998687747670192153", d_atanh(dec_parse("-0.998687747670192153", 0)), "-3.66425056072682"); 119 eq("log 1000", d_log10(dec_int(1000)), "3"); 120 eq("ln e", d_ln(d_e()), "1"); 121 eq("e^1", d_exp(DEC_ONE), "2.71828182845905"); 122 eq("69!", d_fact(dec_int(69)), "1.71122452428141E98"); 123 eq("70!", d_fact(dec_int(70)), "Error"); 124 eq("10C3", d_ncr(dec_int(10), dec_int(3)), "1.2E2"); 125 eq("10P3", d_npr(dec_int(10), dec_int(3)), "7.2E2"); 126 eq("2^10", d_pow(dec_int(2), dec_int(10)), "1.024E3"); 127 eq("2^0.5", d_pow(dec_int(2), dec_parse("0.5", 0)), "1.4142135623731"); 128 eq("cbrt -27", d_root(dec_int(-27), 3), "-3"); 129 eq("(-8)^0.5", d_pow(dec_int(-8), dec_parse("0.5", 0)), "Error"); 130 eq("fmod 725 360", dec_fmod(dec_int(725), dec_int(360)), "5"); 131 eq("round 2.5", dec_round_int(dec_parse("2.5", 0)), "3"); 132 eq("round -2.5", dec_round_int(dec_parse("-2.5", 0)), "-3"); 133 eq("parse 1.2e-7", dec_parse("1.2e-7", 0), "1.2E-7"); 134 eq("1e-99/10", dec_div(dec_parse("1e-99", 0), dec_int(10)), "1E-100"); 135 { 136 int64_t v = 0; 137 st_check(dec_to_int(dec_parse("123456789012345678", 0), &v) && v == 123456789012345678ll, "to_int 18 digitos"); 138 st_check(!dec_to_int(dec_parse("1.5", 0), &v), "to_int 1.5 no es entero"); 139 } 140 141 /* aritmetica al azar */ 142 srand(12345); 143 double worst = 0; 144 for (int i = 0; i < n * 4; i++) { 145 long double a = (rand() / (long double)RAND_MAX - 0.5L) * powl(10, rand() % 40 - 20); 146 long double b = (rand() / (long double)RAND_MAX - 0.5L) * powl(10, rand() % 40 - 20); 147 Dec A = from_ld(a), B = from_ld(b); 148 a = to_ld(A); b = to_ld(B); 149 double e1 = ulp15(dec_add(A, B), a + b), e2 = ulp15(dec_mul(A, B), a * b), e3 = ulp15(dec_div(A, B), a / b); 150 if (fabsl(a + b) < 1e-3L * fabsl(a)) e1 = 0; /* cancelacion: el oraculo tambien pierde */ 151 double e = e1 > e2 ? e1 : e2; 152 if (e3 > e) e = e3; 153 if (e > worst) worst = e; 154 } 155 st_check(worst < 0.1, "+ * / al azar: peor error %.3g ulp15", worst); 156 157 /* funciones al azar */ 158 for (size_t f = 0; f < sizeof FNS / sizeof FNS[0]; f++) { 159 const Fn1 *F = &FNS[f]; 160 double w = 0; 161 long double wx = 0; 162 for (int i = 0; i < n; i++) { 163 long double u = rand() / (long double)RAND_MAX, x; 164 if (F->logscale) x = expl(logl(F->lo) + u * (logl(F->hi) - logl(F->lo))); 165 else x = F->lo + u * (F->hi - F->lo); 166 Dec X = from_ld(x); 167 x = to_ld(X); 168 /* cerca de 1, x en long double ya trae error relativo grande para ln: log1p(x-1) */ 169 long double want = F->g == logl && fabsl(x - 1) < 0.2L ? log1pl(to_ld(dec_sub(X, DEC_ONE))) : F->g(x); 170 double e = ulp15(F->f(X), want); 171 /* el x pasado a long double ya trae ~5e-20 de error relativo: se descuenta 172 * segun el numero de condicion |x f'(x) / f(x)| */ 173 long double h = 1e-9L, d = (F->g(x * (1 + h)) - F->g(x * (1 - h))) / (2 * h); 174 double cond = want != 0 ? (double)fabsl(d / want) : 0; 175 e -= cond * 1e-4; 176 if (e > w) { w = e; wx = x; } 177 } 178 st_check(w < 0.5, "%-14s peor error %.3g ulp15 (x = %.18Lg)", F->name, w, wx); 179 } 180 }