juegos

dist.c (6764B)

---

 /* dist.c - distribuciones (como el menu DIST de la fx-991): normal, binomial, Poisson. */
 #include "dist.h"
 #include "dmath.h"
 
 static Dec D(const char *s) { return dec_parse(s, 0); }
 
 /* ln(n!) para n entero >= 0: producto directo hasta 20, Stirling despues */
 static Dec lnfact(int64_t n)
 {
     if (n < 2) return DEC_ZERO;
     if (n <= 20) {
         Dec p = DEC_ONE;
         for (int64_t i = 2; i <= n; i++) p = dec_mul_int(p, i);
         return d_ln(p);
     }
     Dec x = dec_int(n), lx = d_ln(x);
     /* n ln n - n + ln(2πn)/2 + 1/12n - 1/360n³ + 1/1260n⁵ - 1/1680n⁷ */
     Dec r = dec_sub(dec_mul(x, lx), x);
     r = dec_add(r, dec_div_int(d_ln(dec_mul(dec_mul_int(d_pi(), 2), x)), 2));
     Dec x2 = dec_mul(x, x), xi = dec_div(DEC_ONE, x);
     r = dec_add(r, dec_div_int(xi, 12));
     xi = dec_div(xi, x2); r = dec_sub(r, dec_div_int(xi, 360));
     xi = dec_div(xi, x2); r = dec_add(r, dec_div_int(xi, 1260));
     xi = dec_div(xi, x2); r = dec_sub(r, dec_div_int(xi, 1680));
     return r;
 }
 
 /* erfc(x) para x >= 0 */
 static Dec erfc_pos(Dec x)
 {
     if (dec_cmp(x, D("2.5")) < 0) {
         /* 1 - erf(x), erf por serie: 2/√π Σ (-1)^n x^(2n+1) / (n! (2n+1)) */
         Dec x2 = dec_mul(x, x), t = x, sum = x;
         for (int n = 1; n < 200; n++) {
             t = dec_neg(dec_div_int(dec_mul(t, x2), n));
             Dec term = dec_div_int(t, 2 * n + 1);
             sum = dec_add(sum, term);
             if (!term.m || term.e < sum.e - 20) break;
         }
         Dec erf = dec_div(dec_mul_int(sum, 2), dec_sqrt(d_pi()));
         return dec_sub(DEC_ONE, erf);
     }
     if (dec_cmp(x, D("27")) > 0) return DEC_ZERO;
     /* fraccion continua (Lentz): erfc x = e^(-x²)/√π · 1/(x + (1/2)/(x + 1/(x + (3/2)/(x + ...)))) */
     Dec tiny = D("1e-90"), f = x, C = x, Dd = DEC_ZERO;
     for (int k = 1; k < 400; k++) {
         Dec a = dec_div_int(dec_int(k), 2);
         Dd = dec_add(x, dec_mul(a, Dd));
         if (!Dd.m) Dd = tiny;
         C = dec_add(x, dec_div(a, C));
         if (!C.m) C = tiny;
         Dd = dec_div(DEC_ONE, Dd);
         Dec delta = dec_mul(C, Dd);
         f = dec_mul(f, delta);
         Dec d1 = dec_abs(dec_sub(delta, DEC_ONE));
         if (!d1.m || d1.e < -18) break;
     }
     return dec_div(d_exp(dec_neg(dec_mul(x, x))), dec_mul(f, dec_sqrt(d_pi())));
 }
 
 /* Φ(z) y su complemento, cada uno preciso en su cola */
 static Dec phi(Dec z)
 {
     Dec r2 = dec_sqrt(dec_int(2));
     Dec e = erfc_pos(dec_div(dec_abs(z), r2));
     return z.neg ? dec_div_int(e, 2) : dec_sub(DEC_ONE, dec_div_int(e, 2));
 }
 
 static Dec npdf(Dec z)       /* densidad normal estandar */
 {
     return dec_div(d_exp(dec_neg(dec_div_int(dec_mul(z, z), 2))), dec_sqrt(dec_mul_int(d_pi(), 2)));
 }
 
 static bool sigma_ok(Dec s) { return !s.err && s.m && !s.neg; }
 
 Dec dist_npd(Dec x, Dec s, Dec mu)
 {
     if (!sigma_ok(s)) return dec_err();
     return dec_div(npdf(dec_div(dec_sub(x, mu), s)), s);
 }
 
 Dec dist_ncd(Dec lo, Dec hi, Dec s, Dec mu)
 {
     if (!sigma_ok(s)) return dec_err();
     Dec a = dec_div(dec_sub(lo, mu), s), b = dec_div(dec_sub(hi, mu), s);
     if (dec_cmp(a, b) > 0) { Dec t = a; a = b; b = t; }
     /* en la cola derecha conviene restar complementos */
     Dec r;
     if (!a.neg && a.m) {
         Dec r2 = dec_sqrt(dec_int(2));
         r = dec_div_int(dec_sub(erfc_pos(dec_div(a, r2)), erfc_pos(dec_div(b, r2))), 2);
     } else r = dec_sub(phi(b), phi(a));
     return dec_round_sig(r, 15);
 }
 
 Dec dist_invn(Dec p, Dec s, Dec mu)
 {
     if (!sigma_ok(s) || p.err || p.neg || !p.m || dec_cmp(p, DEC_ONE) >= 0) return dec_err();
     /* primera aproximacion (Abramowitz-Stegun 26.2.23) y despues Newton sobre Φ */
     bool low = dec_cmp(p, D("0.5")) < 0;
     Dec q = low ? p : dec_sub(DEC_ONE, p);
     Dec t = dec_sqrt(dec_mul_int(d_ln(q), -2));
     Dec num = dec_add(D("2.515517"), dec_add(dec_mul(D("0.802853"), t), dec_mul(D("0.010328"), dec_mul(t, t))));
     Dec den = dec_add(DEC_ONE, dec_add(dec_mul(D("1.432788"), t),
                        dec_add(dec_mul(D("0.189269"), dec_mul(t, t)), dec_mul(D("0.001308"), dec_mul(dec_mul(t, t), t)))));
     Dec z = dec_sub(t, dec_div(num, den));
     if (low) z = dec_neg(z);
     for (int i = 0; i < 12; i++) {
         Dec f = dec_sub(phi(z), p), d = npdf(z);
         if (!d.m) break;
         Dec step = dec_div(f, d);
         z = dec_sub(z, step);
         if (!step.m || step.e < -17) break;
     }
     return dec_round_sig(dec_add(mu, dec_mul(s, z)), 15);
 }
 
 /* P(X = k) con log para no desbordar */
 static Dec binom_term(int64_t k, int64_t n, Dec lp, Dec lq)
 {
     Dec l = dec_sub(dec_sub(lnfact(n), lnfact(k)), lnfact(n - k));
     l = dec_add(l, dec_add(dec_mul_int(lp, k), dec_mul_int(lq, n - k)));
     return d_exp(l);
 }
 
 static bool binom_args(Dec x, Dec n, Dec p, int64_t *k, int64_t *N)
 {
     return !p.err && !p.neg && dec_cmp(p, DEC_ONE) <= 0 && dec_to_int(n, N) && *N >= 0 && *N <= 100000 &&
            dec_to_int(x, k);
 }
 
 /* casos borde p = 0 y p = 1 */
 static bool binom_edge(int64_t k, int64_t N, Dec p, bool cum, Dec *out)
 {
     if (p.m && dec_cmp(p, DEC_ONE) != 0) return false;
     int64_t at = p.m ? N : 0;                       /* todo el peso en 0 o en N */
     *out = (cum ? k >= at : k == at) ? DEC_ONE : DEC_ZERO;
     return true;
 }
 
 Dec dist_bpd(Dec x, Dec n, Dec p)
 {
     int64_t k, N;
     if (!binom_args(x, n, p, &k, &N)) return dec_err();
     if (k < 0 || k > N) return DEC_ZERO;
     Dec r;
     if (binom_edge(k, N, p, false, &r)) return r;
     return dec_round_sig(binom_term(k, N, d_ln(p), d_ln(dec_sub(DEC_ONE, p))), 15);
 }
 
 Dec dist_bcd(Dec x, Dec n, Dec p)
 {
     int64_t k, N;
     if (!binom_args(x, n, p, &k, &N)) return dec_err();
     if (k < 0) return DEC_ZERO;
     if (k >= N) return DEC_ONE;
     Dec r;
     if (binom_edge(k, N, p, true, &r)) return r;
     if (k > 20000) return dec_err();
     Dec lp = d_ln(p), lq = d_ln(dec_sub(DEC_ONE, p)), s = DEC_ZERO;
     for (int64_t i = 0; i <= k; i++) s = dec_add(s, binom_term(i, N, lp, lq));
     if (dec_cmp(s, DEC_ONE) > 0) s = DEC_ONE;
     return dec_round_sig(s, 15);
 }
 
 static Dec pois_term(int64_t k, Dec lam, Dec ll)
 {
     return d_exp(dec_sub(dec_sub(dec_mul_int(ll, k), lam), lnfact(k)));
 }
 
 Dec dist_ppd(Dec x, Dec lam)
 {
     int64_t k;
     if (lam.err || lam.neg || !lam.m || !dec_to_int(x, &k)) return dec_err();
     if (k < 0) return DEC_ZERO;
     return dec_round_sig(pois_term(k, lam, d_ln(lam)), 15);
 }
 
 Dec dist_pcd(Dec x, Dec lam)
 {
     int64_t k;
     if (lam.err || lam.neg || !lam.m || !dec_to_int(x, &k)) return dec_err();
     if (k < 0) return DEC_ZERO;
     if (k > 20000) return dec_err();
     Dec ll = d_ln(lam), s = DEC_ZERO;
     for (int64_t i = 0; i <= k; i++) s = dec_add(s, pois_term(i, lam, ll));
     if (dec_cmp(s, DEC_ONE) > 0) s = DEC_ONE;
     return dec_round_sig(s, 15);
 }