dmath.c (13979B)
1 /* dmath.c - funciones trascendentes en decimal (ver dmath.h). 2 * 3 * Reduccion de argumento en Dec con constantes partidas en dos (Cody-Waite) y 4 * series de Taylor sobre argumentos chicos. Los 3 digitos de guarda de Dec 5 * absorben el error; el resultado se redondea a 15 al mostrarlo o guardarlo. */ 6 #include "dmath.h" 7 8 enum { 9 K_PI, K_PI_LO, K_PI2, K_PI2_LO, K_LN10, K_LN10_LO, K_LN2, K_LN2_LO, K_E, K_SQRT3, 10 K_PI6, K_P1, K_P2, K_P3, K_PI180, K_PI200, K_LN10I, K_COUNT 11 }; 12 static const char *const KSTR[K_COUNT] = { 13 "3.14159265358979323", "8.46264338327950288e-18", 14 "1.57079632679489661", "9.23132169163975144e-18", 15 "2.30258509299404", "5.68401799145468436e-15", /* hi de 15 digitos: k*hi es exacto */ 16 "0.693147180559945", "3.09417232121458177e-16", 17 "2.71828182845904523", "1.73205080756887729", 18 "0.523598775598298873", 19 "1.5707963", "2.6794897e-8", "-3.80768678308360249e-16", /* pi/2 en tres partes */ 20 "0.0174532925199432957", "0.0157079632679489661", "0.434294481903251827", 21 }; 22 static Dec K[K_COUNT]; 23 static bool k_ready; 24 25 static const Dec *kon(void) 26 { 27 if (!k_ready) { 28 for (int i = 0; i < K_COUNT; i++) K[i] = dec_parse(KSTR[i], 0); 29 k_ready = true; 30 } 31 return K; 32 } 33 34 Dec d_pi(void) { return kon()[K_PI]; } 35 Dec d_e(void) { return kon()[K_E]; } 36 37 /* el termino ya no cambia la suma (con 2 digitos de margen) */ 38 static bool tiny(Dec term, Dec sum) { return !term.m || (sum.m && term.e < sum.e - 20); } 39 40 static Dec half(Dec x) { return dec_div_int(x, 2); } 41 42 /* ------------------------------------------------------------ exp / ln */ 43 44 Dec d_exp(Dec x) 45 { 46 const Dec *k = kon(); 47 if (x.err) return x; 48 if (dec_cmp(x, dec_int(231)) > 0) return dec_err(); 49 if (dec_cmp(x, dec_int(-232)) < 0) return DEC_ZERO; 50 int64_t n; 51 dec_to_int(dec_round_int(dec_div(x, k[K_LN10])), &n); 52 Dec r = dec_sub(dec_sub(x, dec_mul_int(k[K_LN10], n)), dec_mul_int(k[K_LN10_LO], n)); 53 /* exp(r), |r| <= ~1.16: Taylor */ 54 Dec sum = DEC_ONE, term = DEC_ONE; 55 for (int i = 1; i < 40; i++) { 56 term = dec_div_int(dec_mul(term, r), i); 57 sum = dec_add(sum, term); 58 if (tiny(term, sum)) break; 59 } 60 return dec_scale10(sum, (int)n); 61 } 62 63 /* 2*(s + s^3/3 + s^5/5 + ...) = ln((1+s)/(1-s)), para |s| <= 1/3 */ 64 static Dec atanh_series(Dec s) 65 { 66 Dec s2 = dec_mul(s, s), p = s, sum = s; 67 for (int i = 3; i < 200; i += 2) { 68 p = dec_mul(p, s2); 69 Dec t = dec_div_int(p, i); 70 sum = dec_add(sum, t); 71 if (tiny(t, sum)) break; 72 } 73 return dec_add(sum, sum); 74 } 75 76 Dec d_log1p(Dec x) 77 { 78 if (x.err) return x; 79 Dec y = dec_add(DEC_ONE, x); 80 if (y.neg || !y.m) return dec_err(); 81 /* cerca de 0: s = x/(2+x) sale directo de x, sin perder digitos */ 82 if (dec_cmp_abs(x, dec_parse("0.5", 0)) <= 0) 83 return atanh_series(dec_div(x, dec_add(dec_int(2), x))); 84 return d_ln(y); 85 } 86 87 Dec d_ln(Dec x) 88 { 89 const Dec *k = kon(); 90 if (x.err || x.neg || !x.m) return dec_err(); 91 /* entre 0.5 y 2: directo, preciso cerca de 1 */ 92 if (dec_cmp(x, dec_parse("0.5", 0)) >= 0 && dec_cmp(x, dec_int(2)) <= 0) 93 return atanh_series(dec_div(dec_sub(x, DEC_ONE), dec_add(x, DEC_ONE))); 94 int E = x.e, j = 0; 95 Dec y = x; 96 y.e = 0; /* y en [1, 10) */ 97 Dec lim = dec_parse("1.41421356", 0); 98 while (dec_cmp(y, lim) > 0) { y = half(y); j++; } 99 Dec s = atanh_series(dec_div(dec_sub(y, DEC_ONE), dec_add(y, DEC_ONE))); 100 Dec hi = dec_add(dec_mul_int(k[K_LN10], E), dec_mul_int(k[K_LN2], j)); 101 Dec lo = dec_add(dec_mul_int(k[K_LN10_LO], E), dec_mul_int(k[K_LN2_LO], j)); 102 return dec_add(hi, dec_add(s, lo)); 103 } 104 105 Dec d_log10(Dec x) 106 { 107 if (x.err || x.neg || !x.m) return dec_err(); 108 if (x.m == 100000000000000000ull) return dec_int(x.e); /* potencia de 10: exacto */ 109 return dec_mul(d_ln(x), kon()[K_LN10I]); 110 } 111 112 Dec d_pow10(Dec x) 113 { 114 int64_t n; 115 if (x.err) return x; 116 if (dec_to_int(x, &n)) { 117 if (n > DEC_EMAX) return dec_err(); 118 if (n < -DEC_EMAX - 1) return DEC_ZERO; 119 return dec_scale10(DEC_ONE, (int)n); 120 } 121 const Dec *k = kon(); 122 return d_exp(dec_add(dec_mul(x, k[K_LN10]), dec_mul(x, k[K_LN10_LO]))); 123 } 124 125 Dec d_root(Dec x, int64_t n) 126 { 127 if (x.err || n == 0) return dec_err(); 128 if (n < 0) return dec_div(DEC_ONE, d_root(x, -n)); 129 if (n == 1 || !x.m) return x; 130 if (n == 2) return dec_sqrt(x); 131 bool neg = x.neg; 132 if (neg && !(n & 1)) return dec_err(); 133 Dec a = dec_abs(x); 134 Dec r = d_exp(dec_div_int(d_ln(a), n)); 135 /* un paso de Newton: r -= (r^n - a) / (n r^(n-1)) */ 136 Dec rn1 = dec_pow_int(r, n - 1); 137 Dec f = dec_sub(dec_mul(rn1, r), a); 138 r = dec_sub(r, dec_div(f, dec_mul_int(rn1, n))); 139 if (neg) r = dec_neg(r); 140 return r; 141 } 142 143 Dec d_pow(Dec x, Dec y) 144 { 145 int64_t n; 146 if (x.err || y.err) return dec_err(); 147 if (dec_to_int(y, &n) && n > -100000 && n < 100000) { 148 if (!x.m && n <= 0) return dec_err(); 149 return dec_pow_int(x, n); 150 } 151 if (!x.m) return y.neg || !y.m ? dec_err() : DEC_ZERO; 152 if (x.neg) return dec_err(); 153 return d_exp(dec_mul(y, d_ln(x))); 154 } 155 156 /* ------------------------------------------------------------ trigonometria */ 157 158 static Dec sin_k(Dec r) /* |r| <= pi/4 */ 159 { 160 Dec r2 = dec_mul(r, r), t = r, sum = r; 161 for (int i = 2; i < 60; i += 2) { 162 t = dec_neg(dec_div_int(dec_mul(t, r2), (int64_t)i * (i + 1))); 163 sum = dec_add(sum, t); 164 if (tiny(t, sum)) break; 165 } 166 return sum; 167 } 168 169 static Dec cos_k(Dec r) 170 { 171 Dec r2 = dec_mul(r, r), t = DEC_ONE, sum = DEC_ONE; 172 for (int i = 1; i < 60; i += 2) { 173 t = dec_neg(dec_div_int(dec_mul(t, r2), (int64_t)i * (i + 1))); 174 sum = dec_add(sum, t); 175 if (tiny(t, sum)) break; 176 } 177 return sum; 178 } 179 180 /* x = q*(pi/2) + r con |r| <= pi/4. Devuelve q mod 4 (o -1 si esta fuera de rango). 181 * zero = r es exactamente 0 (multiplos exactos en grados/gradianes). */ 182 static int reduce(Dec x, int ang, Dec *r, bool *zero) 183 { 184 const Dec *k = kon(); 185 int64_t q; 186 *zero = false; 187 if (ang == ANG_RAD) { 188 if (x.e >= 10) return -1; 189 dec_to_int(dec_round_int(dec_div(x, k[K_PI2])), &q); 190 *r = dec_sub(dec_sub(dec_sub(x, dec_mul_int(k[K_P1], q)), dec_mul_int(k[K_P2], q)), 191 dec_mul_int(k[K_P3], q)); 192 } else { 193 Dec full = dec_int(ang == ANG_DEG ? 360 : 400), quarter = dec_int(ang == ANG_DEG ? 90 : 100); 194 if (x.e >= 12) return -1; 195 Dec y = dec_fmod(x, full); /* exacto en decimal */ 196 dec_to_int(dec_round_int(dec_div(y, quarter)), &q); 197 Dec d = dec_sub(y, dec_mul_int(quarter, q)); 198 *zero = !d.m; 199 *r = dec_mul(d, k[ang == ANG_DEG ? K_PI180 : K_PI200]); 200 } 201 return (int)(((q % 4) + 4) % 4); 202 } 203 204 static Dec trig(Dec x, int ang, int which) /* 0 sin, 1 cos, 2 tan */ 205 { 206 if (x.err) return x; 207 Dec r; 208 bool zero; 209 int q = reduce(x, ang, &r, &zero); 210 if (q < 0) return dec_err(); 211 Dec s, c; 212 if (zero) { s = DEC_ZERO; c = DEC_ONE; } 213 else { s = sin_k(r); c = cos_k(r); } 214 Dec S, C; 215 switch (q) { 216 case 0: S = s; C = c; break; 217 case 1: S = c; C = dec_neg(s); break; 218 case 2: S = dec_neg(s); C = dec_neg(c); break; 219 default: S = dec_neg(c); C = s; break; 220 } 221 Dec res; 222 if (which == 0) res = S; 223 else if (which == 1) res = C; 224 else { 225 if (!C.m) return dec_err(); 226 res = dec_div(S, C); 227 } 228 /* sin(pi) en radianes da ~1e-18: es cero (como la Casio) */ 229 if (res.m && ang == ANG_RAD) { 230 int lim = -16 + (x.e > 0 ? x.e : 0); 231 if (res.e < lim) res = DEC_ZERO; 232 } 233 return res; 234 } 235 236 Dec d_sin(Dec x, int ang) { return trig(x, ang, 0); } 237 Dec d_cos(Dec x, int ang) { return trig(x, ang, 1); } 238 Dec d_tan(Dec x, int ang) { return trig(x, ang, 2); } 239 240 Dec d_to_rad(Dec x, int ang) 241 { 242 if (ang == ANG_RAD) return x; 243 return dec_mul(x, kon()[ang == ANG_DEG ? K_PI180 : K_PI200]); 244 } 245 246 Dec d_from_rad(Dec x, int ang) 247 { 248 if (ang == ANG_RAD || x.err) return x; 249 const Dec *k = kon(); 250 /* dividir por pi con 36 digitos: x / (pi_hi + pi_lo) ~ (x/pi_hi) * (1 - pi_lo/pi_hi) */ 251 Dec q = dec_div(x, k[K_PI]); 252 q = dec_sub(q, dec_mul(q, dec_div(k[K_PI_LO], k[K_PI]))); 253 return dec_mul_int(q, ang == ANG_DEG ? 180 : 200); 254 } 255 256 Dec d_ang_conv(Dec x, int from, int to) 257 { 258 if (from == to) return x; 259 if (from == ANG_DEG && to == ANG_GRA) return dec_div_int(dec_mul_int(x, 10), 9); 260 if (from == ANG_GRA && to == ANG_DEG) return dec_div_int(dec_mul_int(x, 9), 10); 261 return d_from_rad(d_to_rad(x, from), to); 262 } 263 264 static Dec atan_k(Dec t) /* |t| <= 2 - sqrt(3) */ 265 { 266 Dec t2 = dec_mul(t, t), p = t, sum = t; 267 for (int i = 3; i < 200; i += 2) { 268 p = dec_neg(dec_mul(p, t2)); 269 Dec term = dec_div_int(p, i); 270 sum = dec_add(sum, term); 271 if (tiny(term, sum)) break; 272 } 273 return sum; 274 } 275 276 /* atan en radianes */ 277 static Dec atan_r(Dec x) 278 { 279 const Dec *k = kon(); 280 if (x.err) return x; 281 bool neg = x.neg; 282 Dec t = dec_abs(x), add = DEC_ZERO; 283 bool inv = false; 284 if (dec_cmp(t, DEC_ONE) > 0) { t = dec_div(DEC_ONE, t); inv = true; } 285 if (dec_cmp(t, dec_parse("0.267949192431122706", 0)) > 0) { 286 /* atan t = pi/6 + atan((sqrt3 t - 1)/(sqrt3 + t)) */ 287 t = dec_div(dec_sub(dec_mul(k[K_SQRT3], t), DEC_ONE), dec_add(k[K_SQRT3], t)); 288 add = k[K_PI6]; 289 } 290 Dec r = dec_add(add, atan_k(t)); 291 if (inv) r = dec_add(dec_sub(k[K_PI2], r), k[K_PI2_LO]); 292 return neg ? dec_neg(r) : r; 293 } 294 295 Dec d_atan(Dec x, int ang) { return d_from_rad(atan_r(x), ang); } 296 297 Dec d_asin(Dec x, int ang) 298 { 299 if (x.err) return x; 300 int c = dec_cmp_abs(x, DEC_ONE); 301 if (c > 0) return dec_err(); 302 Dec r; 303 if (c == 0) { r = kon()[K_PI2]; if (x.neg) r = dec_neg(r); } 304 else { 305 Dec d = dec_sqrt(dec_mul(dec_sub(DEC_ONE, x), dec_add(DEC_ONE, x))); 306 r = atan_r(dec_div(x, d)); 307 } 308 if (ang != ANG_RAD && c == 0) return dec_int(x.neg ? (ang == ANG_DEG ? -90 : -100) : (ang == ANG_DEG ? 90 : 100)); 309 return d_from_rad(r, ang); 310 } 311 312 Dec d_acos(Dec x, int ang) 313 { 314 if (x.err) return x; 315 int c = dec_cmp_abs(x, DEC_ONE); 316 if (c > 0) return dec_err(); 317 if (c == 0 && x.neg) 318 return ang == ANG_RAD ? d_pi() : dec_int(ang == ANG_DEG ? 180 : 200); 319 if (c == 0) return DEC_ZERO; 320 Dec r = atan_r(dec_sqrt(dec_div(dec_sub(DEC_ONE, x), dec_add(DEC_ONE, x)))); 321 return d_from_rad(dec_add(r, r), ang); 322 } 323 324 Dec d_atan2(Dec y, Dec x, int ang) 325 { 326 const Dec *k = kon(); 327 if (x.err || y.err) return dec_err(); 328 if (!x.m && !y.m) return dec_err(); 329 Dec r; 330 if (!x.m) r = y.neg ? dec_neg(k[K_PI2]) : k[K_PI2]; 331 else { 332 r = atan_r(dec_div(y, x)); 333 if (x.neg) r = y.neg ? dec_sub(r, k[K_PI]) : dec_add(r, k[K_PI]); 334 } 335 return d_from_rad(r, ang); 336 } 337 338 /* ------------------------------------------------------------ hiperbolicas */ 339 340 Dec d_sinh(Dec x) 341 { 342 if (x.err) return x; 343 if (dec_cmp_abs(x, DEC_ONE) < 0) { /* serie: sin cancelacion cerca de 0 */ 344 Dec x2 = dec_mul(x, x), t = x, sum = x; 345 for (int i = 2; i < 60; i += 2) { 346 t = dec_div_int(dec_mul(t, x2), (int64_t)i * (i + 1)); 347 sum = dec_add(sum, t); 348 if (tiny(t, sum)) break; 349 } 350 return sum; 351 } 352 Dec e = d_exp(x); 353 if (e.err) return e; 354 return half(dec_sub(e, dec_div(DEC_ONE, e))); 355 } 356 357 Dec d_cosh(Dec x) 358 { 359 if (x.err) return x; 360 Dec e = d_exp(dec_abs(x)); 361 if (e.err) return e; 362 return half(dec_add(e, dec_div(DEC_ONE, e))); 363 } 364 365 Dec d_tanh(Dec x) 366 { 367 if (x.err) return x; 368 if (x.e >= 2) return x.neg ? dec_neg(DEC_ONE) : DEC_ONE; /* |x| >= 10: 1 en 18 digitos... casi */ 369 Dec s = d_sinh(x), c = d_cosh(x); 370 return dec_div(s, c); 371 } 372 373 Dec d_asinh(Dec x) 374 { 375 if (x.err) return x; 376 bool neg = x.neg; 377 Dec a = dec_abs(x); 378 Dec r; 379 if (a.e >= 9) r = dec_add(d_ln(a), d_ln(dec_int(2))); /* x grande: ln(2x) */ 380 else { 381 /* asinh a = log1p(a + a^2/(1 + sqrt(1 + a^2))) */ 382 Dec a2 = dec_mul(a, a); 383 r = d_log1p(dec_add(a, dec_div(a2, dec_add(DEC_ONE, dec_sqrt(dec_add(DEC_ONE, a2)))))); 384 } 385 return neg ? dec_neg(r) : r; 386 } 387 388 Dec d_acosh(Dec x) 389 { 390 if (x.err) return x; 391 if (dec_cmp(x, DEC_ONE) < 0) return dec_err(); 392 if (x.e >= 9) return dec_add(d_ln(x), d_ln(dec_int(2))); 393 Dec u = dec_sub(x, DEC_ONE); /* log1p(u + sqrt(u(u+2))) */ 394 return d_log1p(dec_add(u, dec_sqrt(dec_mul(u, dec_add(u, dec_int(2)))))); 395 } 396 397 Dec d_atanh(Dec x) 398 { 399 if (x.err) return x; 400 if (dec_cmp_abs(x, DEC_ONE) >= 0) return dec_err(); 401 /* atanh x = log1p(2x/(1-x)) / 2, con x >= 0 (para x negativo 1+u cancelaria) */ 402 Dec a = dec_abs(x); 403 Dec r = half(d_log1p(dec_div(dec_add(a, a), dec_sub(DEC_ONE, a)))); 404 return x.neg ? dec_neg(r) : r; 405 } 406 407 /* ------------------------------------------------------------ combinatoria */ 408 409 Dec d_fact(Dec x) 410 { 411 int64_t n; 412 if (x.err || !dec_to_int(x, &n) || n < 0 || n > 69) return dec_err(); 413 Dec r = DEC_ONE; 414 for (int64_t i = 2; i <= n; i++) r = dec_mul_int(r, i); 415 return r; 416 } 417 418 static bool nr_args(Dec n, Dec r, int64_t *pn, int64_t *pr) 419 { 420 return dec_to_int(n, pn) && dec_to_int(r, pr) && *pn >= 0 && *pr >= 0 && *pr <= *pn 421 && *pn < 10000000000ll; 422 } 423 424 Dec d_npr(Dec n, Dec r) 425 { 426 int64_t a, b; 427 if (n.err || r.err || !nr_args(n, r, &a, &b)) return dec_err(); 428 Dec res = DEC_ONE; 429 for (int64_t i = 0; i < b; i++) { 430 res = dec_mul_int(res, a - i); 431 if (res.err) return res; 432 } 433 return res; 434 } 435 436 Dec d_ncr(Dec n, Dec r) 437 { 438 int64_t a, b; 439 if (n.err || r.err || !nr_args(n, r, &a, &b)) return dec_err(); 440 if (b > a - b) b = a - b; 441 Dec res = DEC_ONE; 442 for (int64_t i = 1; i <= b; i++) { 443 res = dec_div_int(dec_mul_int(res, a - b + i), i); 444 if (res.err) return res; 445 } 446 return res; 447 }