juegos

exact.c (10407B)

---

 /* exact.c - aritmetica exacta (a + b·√r)/d · π^p (ver exact.h). */
 #include "exact.h"
 #include "dmath.h"
 
 #define OVF(op) do { if (op) return false; } while (0)
 
 static int64_t iabs(int64_t v) { return v < 0 ? -v : v; }
 
 int64_t i64_gcd(int64_t a, int64_t b)
 {
     uint64_t x = (uint64_t)(a < 0 ? -(uint64_t)a : (uint64_t)a), y = (uint64_t)(b < 0 ? -(uint64_t)b : (uint64_t)b);
     while (y) { uint64_t t = x % y; x = y; y = t; }
     return (int64_t)x;
 }
 
 /* saca los cuadrados de r: r = s²·r' (division de prueba hasta 10^5) */
 static bool squarefree(int64_t *r, int64_t *s)
 {
     *s = 1;
     int64_t v = *r;
     for (int64_t p = 2; p <= 100000 && p * p <= v; p += (p == 2 ? 1 : 2)) {
         int64_t pp = p * p;
         while (v % pp == 0) { v /= pp; OVF(__builtin_mul_overflow(*s, p, s)); }
     }
     *r = v;
     return true;
 }
 
 bool ex_norm(Ex *x)
 {
     if (x->d == 0) return false;
     if (x->d < 0) {
         if (x->d == INT64_MIN || x->a == INT64_MIN || x->b == INT64_MIN) return false;
         x->d = -x->d; x->a = -x->a; x->b = -x->b;
     }
     if (x->b == 0 || x->r == 0) { x->b = 0; x->r = 1; }
     else if (x->r < 0) return false;                      /* raiz de negativo: no es real */
     else {
         int64_t s;
         if (!squarefree(&x->r, &s)) return false;
         OVF(__builtin_mul_overflow(x->b, s, &x->b));
         if (x->r == 1) { OVF(__builtin_add_overflow(x->a, x->b, &x->a)); x->b = 0; }
     }
     if (x->b && x->p) return false;                       /* π·√r: no lo representamos */
     int64_t g = i64_gcd(i64_gcd(x->a, x->b), x->d);
     if (g > 1) { x->a /= g; x->b /= g; x->d /= g; }
     if (x->a == 0 && x->b == 0) { x->d = 1; x->p = 0; }
     return true;
 }
 
 bool ex_rat(Ex *x, int64_t a, int64_t d)
 {
     x->a = a; x->b = 0; x->r = 1; x->d = d; x->p = 0;
     return ex_norm(x);
 }
 
 bool ex_eq(const Ex *x, const Ex *y)
 {
     return x->a == y->a && x->b == y->b && x->d == y->d && x->p == y->p && (x->b == 0 || x->r == y->r);
 }
 
 int ex_sign(const Ex *x)
 {
     Dec v = ex_to_dec(x);
     return !v.m ? 0 : v.neg ? -1 : 1;
 }
 
 bool ex_neg(Ex *o, const Ex *x)
 {
     if (x->a == INT64_MIN || x->b == INT64_MIN) return false;
     *o = *x;
     o->a = -x->a; o->b = -x->b;
     return true;
 }
 
 bool ex_add(Ex *o, const Ex *x, const Ex *y)
 {
     if (ex_is_zero(x)) { *o = *y; return true; }
     if (ex_is_zero(y)) { *o = *x; return true; }
     if (x->p != y->p) return false;
     if (x->b && y->b && x->r != y->r) return false;
     Ex t;
     int64_t u, v;
     t.p = x->p;
     t.r = x->b ? x->r : y->r;
     OVF(__builtin_mul_overflow(x->a, y->d, &u));
     OVF(__builtin_mul_overflow(y->a, x->d, &v));
     OVF(__builtin_add_overflow(u, v, &t.a));
     OVF(__builtin_mul_overflow(x->b, y->d, &u));
     OVF(__builtin_mul_overflow(y->b, x->d, &v));
     OVF(__builtin_add_overflow(u, v, &t.b));
     OVF(__builtin_mul_overflow(x->d, y->d, &t.d));
     if (!ex_norm(&t)) return false;
     *o = t;
     return true;
 }
 
 bool ex_sub(Ex *o, const Ex *x, const Ex *y)
 {
     Ex n;
     return ex_neg(&n, y) && ex_add(o, x, &n);
 }
 
 bool ex_mul(Ex *o, const Ex *x, const Ex *y)
 {
     if (ex_is_zero(x) || ex_is_zero(y)) return ex_rat(o, 0, 1);
     if (x->p + y->p > 1) return false;
     Ex t;
     int64_t u, v;
     t.p = (uint8_t)(x->p + y->p);
     OVF(__builtin_mul_overflow(x->d, y->d, &t.d));
     if (!x->b || !y->b) {
         const Ex *s = x->b ? x : y, *q = x->b ? y : x;      /* q es racional (o π racional) */
         OVF(__builtin_mul_overflow(s->a, q->a, &t.a));
         OVF(__builtin_mul_overflow(s->b, q->a, &t.b));
         t.r = s->r;
     } else if (x->r == y->r) {
         /* (a1 + b1√r)(a2 + b2√r) = a1a2 + b1b2 r + (a1b2 + a2b1)√r */
         int64_t w;
         OVF(__builtin_mul_overflow(x->a, y->a, &u));
         OVF(__builtin_mul_overflow(x->b, y->b, &v));
         OVF(__builtin_mul_overflow(v, x->r, &v));
         OVF(__builtin_add_overflow(u, v, &t.a));
         OVF(__builtin_mul_overflow(x->a, y->b, &u));
         OVF(__builtin_mul_overflow(y->a, x->b, &w));
         OVF(__builtin_add_overflow(u, w, &t.b));
         t.r = x->r;
     } else if (!x->a && !y->a) {                           /* b1√r · b2√s = b1b2 √(rs) */
         t.a = 0;
         OVF(__builtin_mul_overflow(x->b, y->b, &t.b));
         OVF(__builtin_mul_overflow(x->r, y->r, &t.r));
     } else return false;
     if (!ex_norm(&t)) return false;
     *o = t;
     return true;
 }
 
 bool ex_div(Ex *o, const Ex *x, const Ex *y)
 {
     if (ex_is_zero(y)) return false;
     Ex inv, xx = *x;
     if (y->p) {                                            /* (algo·π)/(c·π) */
         if (!x->p) return false;
         xx.p = 0;
     }
     /* 1/((a + b√r)/d) = d(a - b√r)/(a² - b²r) */
     int64_t a2, b2r, den;
     OVF(__builtin_mul_overflow(y->a, y->a, &a2));
     OVF(__builtin_mul_overflow(y->b, y->b, &b2r));
     OVF(__builtin_mul_overflow(b2r, y->r, &b2r));
     OVF(__builtin_sub_overflow(a2, b2r, &den));
     if (den == 0) return false;
     OVF(__builtin_mul_overflow(y->d, y->a, &inv.a));
     OVF(__builtin_mul_overflow(y->d, -y->b, &inv.b));
     inv.r = y->r; inv.d = den; inv.p = 0;
     if (!ex_norm(&inv)) return false;
     return ex_mul(o, &xx, &inv);
 }
 
 bool ex_pow_int(Ex *o, const Ex *x, int64_t n)
 {
     if (n < 0) {
         Ex one, p;
         ex_rat(&one, 1, 1);
         if (n == INT64_MIN || !ex_pow_int(&p, x, -n)) return false;
         return ex_div(o, &one, &p);
     }
     if (n > 200) return false;
     Ex r, b = *x;
     ex_rat(&r, 1, 1);
     while (n) {
         if (n & 1) { if (!ex_mul(&r, &r, &b)) return false; }
         n >>= 1;
         if (n && !ex_mul(&b, &b, &b)) return false;
     }
     *o = r;
     return true;
 }
 
 bool ex_sqrt(Ex *o, const Ex *x)
 {
     if (!ex_is_rat(x) || x->a < 0) return false;
     /* √(a/d) = √(a·d)/d */
     Ex t = { 0, 1, 0, x->d, 0 };
     OVF(__builtin_mul_overflow(x->a, x->d, &t.r));
     if (t.r == 0) return ex_rat(o, 0, 1);
     if (!ex_norm(&t)) return false;
     *o = t;
     return true;
 }
 
 /* raiz n-esima entera de v si existe */
 static bool iroot(int64_t v, int64_t n, int64_t *out)
 {
     if (v < 0) return false;
     if (v < 2) { *out = v; return true; }
     Dec r = d_root(dec_int(v), n);
     int64_t c;
     if (!dec_to_int(dec_round_int(r), &c)) return false;
     for (int64_t k = c - 1; k <= c + 1; k++) {
         if (k < 0) continue;
         int64_t p = 1;
         bool ovf = false;
         for (int64_t i = 0; i < n && !ovf; i++) ovf = __builtin_mul_overflow(p, k, &p);
         if (!ovf && p == v) { *out = k; return true; }
     }
     return false;
 }
 
 bool ex_root(Ex *o, const Ex *x, int64_t n)
 {
     if (n == 2) return ex_sqrt(o, x);
     if (!ex_is_rat(x) || n < 1) return false;
     bool neg = x->a < 0;
     if (neg && !(n & 1)) return false;
     int64_t ra, rd;
     if (!iroot(iabs(x->a), n, &ra) || !iroot(x->d, n, &rd)) return false;
     return ex_rat(o, neg ? -ra : ra, rd);
 }
 
 Dec ex_to_dec(const Ex *x)
 {
     Dec v = dec_int(x->a);
     if (x->b) v = dec_add(v, dec_mul(dec_int(x->b), dec_sqrt(dec_int(x->r))));
     v = dec_div(v, dec_int(x->d));
     if (x->p) v = dec_mul(v, d_pi());
     return v;
 }
 
 /* ------------------------------------------------------------ trigonometria exacta */
 
 /* angulo -> grados (racional). Solo sirve si da entero. */
 static bool to_deg(const Ex *x, int ang, int64_t *deg)
 {
     Ex t;
     if (x->b) return false;
     if (ang == ANG_RAD) {
         if (!x->p && x->a) return false;                   /* radianes sin π: no hay exacto */
         Ex k;
         ex_rat(&k, 180, 1);
         t = *x; t.p = 0;
         if (!ex_mul(&t, &t, &k)) return false;
     } else {
         if (x->p) return false;
         t = *x;
         if (ang == ANG_GRA) { Ex k; ex_rat(&k, 9, 10); if (!ex_mul(&t, &t, &k)) return false; }
     }
     if (t.d != 1) return false;
     *deg = ((t.a % 360) + 360) % 360;
     return true;
 }
 
 static bool from_deg(Ex *o, int64_t deg, int ang)
 {
     if (ang == ANG_DEG) return ex_rat(o, deg, 1);
     if (ang == ANG_GRA) return ex_rat(o, deg * 10, 9);
     if (!ex_rat(o, deg, 180)) return false;
     if (deg) o->p = 1;
     return true;
 }
 
 /* sin de 0..90 grados con forma exacta */
 static bool sin_tab(int64_t t, Ex *o)
 {
     switch (t) {
     case 0: return ex_rat(o, 0, 1);
     case 18: *o = (Ex){ -1, 1, 5, 4, 0 }; return true;    /* (√5 - 1)/4 */
     case 30: return ex_rat(o, 1, 2);
     case 45: *o = (Ex){ 0, 1, 2, 2, 0 }; return true;     /* √2/2 */
     case 54: *o = (Ex){ 1, 1, 5, 4, 0 }; return true;     /* (1 + √5)/4 */
     case 60: *o = (Ex){ 0, 1, 3, 2, 0 }; return true;     /* √3/2 */
     case 90: return ex_rat(o, 1, 1);
     }
     return false;
 }
 
 static bool sin_deg(int64_t t, Ex *o)
 {
     t = ((t % 360) + 360) % 360;
     bool neg = t >= 180;
     if (neg) t -= 180;
     if (t > 90) t = 180 - t;
     if (!sin_tab(t, o)) return false;
     return neg ? ex_neg(o, o) : true;
 }
 
 bool ex_sin(Ex *o, const Ex *x, int ang)
 {
     int64_t d;
     return to_deg(x, ang, &d) && sin_deg(d, o);
 }
 
 bool ex_cos(Ex *o, const Ex *x, int ang)
 {
     int64_t d;
     return to_deg(x, ang, &d) && sin_deg(d + 90, o);
 }
 
 bool ex_tan(Ex *o, const Ex *x, int ang, bool *undef)
 {
     int64_t d;
     Ex s, c;
     *undef = false;
     if (!to_deg(x, ang, &d) || !sin_deg(d, &s) || !sin_deg(d + 90, &c)) return false;
     if (ex_is_zero(&c)) { *undef = true; return false; }
     return ex_div(o, &s, &c);
 }
 
 static const int64_t TAB[] = { 0, 18, 30, 45, 54, 60, 90 };
 
 bool ex_asin(Ex *o, const Ex *x, int ang)
 {
     for (int i = 0; i < 7; i++) {
         Ex s, n;
         sin_tab(TAB[i], &s);
         n = s;
         ex_neg(&n, &s);
         if (ex_eq(&s, x)) return from_deg(o, TAB[i], ang);
         if (ex_eq(&n, x)) return from_deg(o, -TAB[i], ang);
     }
     return false;
 }
 
 bool ex_acos(Ex *o, const Ex *x, int ang)
 {
     for (int i = 0; i < 7; i++) {
         Ex s, n;
         sin_tab(TAB[i], &s);
         n = s;
         ex_neg(&n, &s);
         if (ex_eq(&s, x)) return from_deg(o, 90 - TAB[i], ang);
         if (ex_eq(&n, x)) return from_deg(o, 90 + TAB[i], ang);
     }
     return false;
 }
 
 bool ex_atan(Ex *o, const Ex *x, int ang)
 {
     for (int i = 0; i < 6; i++) {
         Ex s, c, t, n;
         sin_tab(TAB[i], &s);
         sin_deg(TAB[i] + 90, &c);
         if (!ex_div(&t, &s, &c)) continue;
         n = t;
         ex_neg(&n, &t);
         if (ex_eq(&t, x)) return from_deg(o, TAB[i], ang);
         if (ex_eq(&n, x)) return from_deg(o, -TAB[i], ang);
     }
     return false;
 }