num.c (6726B)
1 /* num.c - operaciones sobre Num: primero intenta exacto, si no, decimal. */ 2 #include "num.h" 3 #include "dmath.h" 4 5 Num num_dec(Dec v) { Num n; n.v = v; n.exact = 0; n.x.a = 0; n.x.b = 0; n.x.r = 1; n.x.d = 1; n.x.p = 0; return n; } 6 Num num_err(void) { return num_dec(dec_err()); } 7 8 Num num_ex(const Ex *x) 9 { 10 Num n; 11 n.x = *x; 12 n.v = ex_to_dec(x); 13 n.exact = !n.v.err; 14 return n; 15 } 16 17 Num num_int(int64_t v) 18 { 19 Ex x; 20 ex_rat(&x, v, 1); 21 return num_ex(&x); 22 } 23 24 Num num_to_dec(Num a) { a.exact = 0; return a; } 25 26 Num num_lit(const char *s, int exp10) 27 { 28 /* mantisa entera y cantidad de decimales */ 29 int64_t m = 0; 30 int frac = 0, nd = 0; 31 bool dot = false, ok = true; 32 for (const char *p = s; *p; p++) { 33 if (*p == '.') { dot = true; continue; } 34 if (*p < '0' || *p > '9') break; 35 if (!m && *p == '0') { if (dot) frac++; continue; } 36 if (++nd > 18) { ok = false; break; } 37 m = m * 10 + (*p - '0'); 38 if (dot) frac++; 39 } 40 /* valor decimal */ 41 char buf[64]; 42 int k = 0; 43 for (const char *p = s; *p && k < 56; p++) buf[k++] = *p; 44 buf[k] = 0; 45 Dec v = dec_scale10(dec_parse(buf, 0), exp10); 46 if (!ok || v.err) return num_dec(v); 47 int e = exp10 - frac; 48 Ex x; 49 int64_t pw = 1; 50 for (int i = 0; i < (e < 0 ? -e : e); i++) 51 if (__builtin_mul_overflow(pw, 10, &pw)) return num_dec(v); 52 if (e >= 0) { 53 int64_t a; 54 if (__builtin_mul_overflow(m, pw, &a)) return num_dec(v); 55 ex_rat(&x, a, 1); 56 } else ex_rat(&x, m, pw); 57 return num_ex(&x); 58 } 59 60 static Num both(Dec v, bool ok, const Ex *x) 61 { 62 if (ok) { 63 Num n = num_ex(x); 64 if (!n.v.err) return n; 65 } 66 return num_dec(v); 67 } 68 69 Num num_add(Num a, Num b) 70 { 71 Ex x; 72 bool ok = a.exact && b.exact && ex_add(&x, &a.x, &b.x); 73 return both(dec_add(a.v, b.v), ok, &x); 74 } 75 76 Num num_sub(Num a, Num b) 77 { 78 Ex x; 79 bool ok = a.exact && b.exact && ex_sub(&x, &a.x, &b.x); 80 return both(dec_sub(a.v, b.v), ok, &x); 81 } 82 83 Num num_mul(Num a, Num b) 84 { 85 Ex x; 86 bool ok = a.exact && b.exact && ex_mul(&x, &a.x, &b.x); 87 return both(dec_mul(a.v, b.v), ok, &x); 88 } 89 90 Num num_div(Num a, Num b) 91 { 92 if (!b.v.m) return num_err(); 93 Ex x; 94 bool ok = a.exact && b.exact && ex_div(&x, &a.x, &b.x); 95 return both(dec_div(a.v, b.v), ok, &x); 96 } 97 98 Num num_neg(Num a) 99 { 100 Ex x; 101 bool ok = a.exact && ex_neg(&x, &a.x); 102 return both(dec_neg(a.v), ok, &x); 103 } 104 105 Num num_abs(Num a) { return a.v.neg ? num_neg(a) : a; } 106 107 Num num_sqrt(Num a) 108 { 109 if (a.v.neg) return num_err(); 110 Ex x; 111 bool ok = a.exact && ex_sqrt(&x, &a.x); 112 return both(dec_sqrt(a.v), ok, &x); 113 } 114 115 Num num_root(Num idx, Num a) 116 { 117 int64_t n; 118 if (!dec_to_int(idx.v, &n) || n == 0) { 119 /* indice no entero: a^(1/idx) */ 120 return num_pow(a, num_div(num_int(1), idx)); 121 } 122 Ex x; 123 bool ok = a.exact && n > 0 && ex_root(&x, &a.x, n); 124 return both(d_root(a.v, n), ok, &x); 125 } 126 127 Num num_pow(Num a, Num b) 128 { 129 int64_t n; 130 Ex x; 131 if (b.exact && ex_is_int(&b.x) && dec_to_int(b.v, &n)) { 132 if (!a.v.m && n <= 0) return num_err(); 133 bool ok = a.exact && ex_pow_int(&x, &a.x, n); 134 return both(dec_pow_int(a.v, n), ok, &x); 135 } 136 /* exponente racional p/q: raiz q-esima exacta a la potencia p (o negativo con q impar) */ 137 if (b.exact && ex_is_rat(&b.x) && b.x.d <= 1000) { 138 Num r = num_root(num_int(b.x.d), a); 139 if (num_bad(&r)) return r; 140 return num_pow(r, num_int(b.x.a)); 141 } 142 return num_dec(d_pow(a.v, b.v)); 143 } 144 145 Num num_fact(Num a) 146 { 147 Dec v = d_fact(a.v); 148 int64_t n; 149 if (!v.err && dec_to_int(v, &n) && v.e < 18) return num_int(n); 150 return num_dec(v); 151 } 152 153 static Num int_result(Dec v) 154 { 155 int64_t n; 156 if (!v.err && v.e < 18 && dec_to_int(v, &n)) return num_int(n); 157 return num_dec(v); 158 } 159 160 Num num_npr(Num n, Num r) { return int_result(d_npr(n.v, r.v)); } 161 Num num_ncr(Num n, Num r) { return int_result(d_ncr(n.v, r.v)); } 162 163 /* log_b(x) exacto si x = b^k (racionales, k entero chico) */ 164 static bool log_exact(const Num *b, const Num *x, int64_t *k) 165 { 166 if (!b->exact || !x->exact || !ex_is_rat(&b->x) || !ex_is_rat(&x->x)) return false; 167 if (b->x.a <= 0 || x->x.a <= 0 || ex_eq(&b->x, &(Ex){ 1, 0, 1, 1, 0 })) return false; 168 for (int s = -1; s <= 1; s += 2) 169 for (int64_t i = 0; i <= 64; i++) { 170 Ex p; 171 if (!ex_pow_int(&p, &b->x, s * i)) break; 172 if (ex_eq(&p, &x->x)) { *k = s * i; return true; } 173 } 174 return false; 175 } 176 177 Num num_logb(Num b, Num x) 178 { 179 int64_t k; 180 if (log_exact(&b, &x, &k)) return num_int(k); 181 Dec lb = d_ln(b.v), lx = d_ln(x.v); 182 if (lb.err || lx.err || !lb.m) return num_err(); 183 return num_dec(dec_div(lx, lb)); 184 } 185 186 Num num_log10(Num x) 187 { 188 int64_t k; 189 Num ten = num_int(10); 190 if (log_exact(&ten, &x, &k)) return num_int(k); 191 return num_dec(d_log10(x.v)); 192 } 193 194 Num num_ln(Num x) { return num_dec(d_ln(x.v)); } 195 196 Num num_sin(Num x, int ang) 197 { 198 Ex r; 199 bool ok = x.exact && ex_sin(&r, &x.x, ang); 200 return both(d_sin(x.v, ang), ok, &r); 201 } 202 203 Num num_cos(Num x, int ang) 204 { 205 Ex r; 206 bool ok = x.exact && ex_cos(&r, &x.x, ang); 207 return both(d_cos(x.v, ang), ok, &r); 208 } 209 210 Num num_tan(Num x, int ang) 211 { 212 Ex r; 213 bool undef = false; 214 bool ok = x.exact && ex_tan(&r, &x.x, ang, &undef); 215 if (undef) return num_err(); 216 return both(d_tan(x.v, ang), ok, &r); 217 } 218 219 Num num_asin(Num x, int ang) 220 { 221 Ex r; 222 bool ok = x.exact && ex_asin(&r, &x.x, ang); 223 return both(d_asin(x.v, ang), ok, &r); 224 } 225 226 Num num_acos(Num x, int ang) 227 { 228 Ex r; 229 bool ok = x.exact && ex_acos(&r, &x.x, ang); 230 return both(d_acos(x.v, ang), ok, &r); 231 } 232 233 Num num_atan(Num x, int ang) 234 { 235 Ex r; 236 bool ok = x.exact && ex_atan(&r, &x.x, ang); 237 return both(d_atan(x.v, ang), ok, &r); 238 } 239 240 Num num_int_part(Num x) 241 { 242 if (x.exact && ex_is_rat(&x.x)) return num_int(x.x.a / x.x.d); 243 return int_result(dec_trunc(x.v)); 244 } 245 246 Num num_floor(Num x) 247 { 248 if (x.exact && ex_is_rat(&x.x)) { 249 int64_t q = x.x.a / x.x.d; 250 if (x.x.a % x.x.d && x.x.a < 0) q--; 251 return num_int(q); 252 } 253 return int_result(dec_floor(x.v)); 254 } 255 256 Num num_gcd(Num a, Num b) 257 { 258 int64_t x, y; 259 if (!dec_to_int(a.v, &x) || !dec_to_int(b.v, &y)) return num_err(); 260 return num_int(i64_gcd(x, y)); 261 } 262 263 Num num_lcm(Num a, Num b) 264 { 265 int64_t x, y, g, l; 266 if (!dec_to_int(a.v, &x) || !dec_to_int(b.v, &y)) return num_err(); 267 g = i64_gcd(x, y); 268 if (!g) return num_int(0); 269 if (__builtin_mul_overflow(x / g, y, &l)) return num_dec(dec_mul(dec_int(x / g), dec_int(y))); 270 return num_int(l < 0 ? -l : l); 271 } 272 273 int num_cmp(Num a, Num b) { return dec_cmp(a.v, b.v); }