eval.c (22874B)
1 /* eval.c - evalua el programa RPN con una pila estatica de Num. */ 2 #include "eval.h" 3 #include "dmath.h" 4 #include "dist.h" 5 6 #define STACK 64 7 #define CALC_DEPTH 3 /* ∫ dentro de ∫ dentro de ∫ */ 8 9 static Num st[STACK]; /* static: la pila de la Pico es chica */ 10 static int depth; 11 12 static int run(const Prog *p, int from, int to, const EvalCtx *cx, int base, Num *out, int *pos); 13 14 /* evalua el cuerpo [from, to) con X = x */ 15 static bool body_at(const Prog *p, int from, int to, const EvalCtx *cx, int base, Num x, Num *out) 16 { 17 static Num vars[CALC_DEPTH][NVARS]; 18 if (depth >= CALC_DEPTH) return false; 19 EvalCtx sub = *cx; 20 Num *v = vars[depth]; 21 for (int i = 0; i < NVARS; i++) v[i] = cx->vars ? cx->vars[i] : num_int(0); 22 v[VAR_X] = x; 23 sub.vars = v; 24 depth++; 25 int pos, err = run(p, from, to, &sub, base, out, &pos); 26 depth--; 27 return err == CE_OK && !num_bad(out); 28 } 29 30 static Dec body_dec(const Prog *p, int from, int to, const EvalCtx *cx, int base, Dec x, bool *ok) 31 { 32 Num r; 33 if (!body_at(p, from, to, cx, base, num_dec(x), &r)) { *ok = false; return DEC_ZERO; } 34 return r.v; 35 } 36 37 /* Gauss-Kronrod 7-15 (nodos y pesos con 18 cifras) */ 38 static const char *const GK_X[8] = { 39 "0.991455371120812639", "0.949107912342758525", "0.864864423359769073", "0.741531185599394440", 40 "0.586087235467691130", "0.405845151377397167", "0.207784955007898468", "0" }; 41 static const char *const GK_WK[8] = { 42 "0.0229353220105292250", "0.0630920926299785533", "0.104790010322250184", "0.140653259715525919", 43 "0.169004726639267903", "0.190350578064785410", "0.204432940075298892", "0.209482141084727828" }; 44 static const char *const GK_WG[4] = { 45 "0.129484966168869693", "0.279705391489276668", "0.381830050505118945", "0.417959183673469388" }; 46 47 static Dec gk15(const Prog *p, int from, int to, const EvalCtx *cx, int base, Dec a, Dec b, Dec *err, bool *ok) 48 { 49 static Dec X[8], WK[8], WG[4]; 50 static bool init; 51 if (!init) { 52 for (int i = 0; i < 8; i++) { X[i] = dec_parse(GK_X[i], 0); WK[i] = dec_parse(GK_WK[i], 0); } 53 for (int i = 0; i < 4; i++) WG[i] = dec_parse(GK_WG[i], 0); 54 init = true; 55 } 56 Dec c = dec_div_int(dec_add(a, b), 2), h = dec_div_int(dec_sub(b, a), 2); 57 Dec fc = body_dec(p, from, to, cx, base, c, ok); 58 Dec k = dec_mul(WK[7], fc), g = dec_mul(WG[3], fc); 59 for (int i = 0; i < 7 && *ok; i++) { 60 Dec dx = dec_mul(h, X[i]); 61 Dec f1 = body_dec(p, from, to, cx, base, dec_sub(c, dx), ok); 62 Dec f2 = body_dec(p, from, to, cx, base, dec_add(c, dx), ok); 63 Dec s = dec_add(f1, f2); 64 k = dec_add(k, dec_mul(WK[i], s)); 65 if (i & 1) g = dec_add(g, dec_mul(WG[i / 2], s)); 66 } 67 k = dec_mul(k, h); 68 g = dec_mul(g, h); 69 *err = dec_abs(dec_sub(k, g)); 70 return k; 71 } 72 73 /* adaptativo: parte el intervalo con mas error hasta llegar a la tolerancia */ 74 static Num integrate(const Prog *p, int from, int to, const EvalCtx *cx, int base, Num lo, Num hi) 75 { 76 enum { SEGS = 24 }; 77 static Dec sa[CALC_DEPTH][SEGS], sb[CALC_DEPTH][SEGS], sv[CALC_DEPTH][SEGS], se[CALC_DEPTH][SEGS]; 78 int d = depth < CALC_DEPTH ? depth : CALC_DEPTH - 1, n = 1; 79 bool ok = true; 80 Dec *A = sa[d], *B = sb[d], *V = sv[d], *Er = se[d]; 81 if (!dec_cmp(lo.v, hi.v)) return num_int(0); 82 A[0] = lo.v; B[0] = hi.v; 83 V[0] = gk15(p, from, to, cx, base, A[0], B[0], &Er[0], &ok); 84 for (;;) { 85 Dec tot = DEC_ZERO, terr = DEC_ZERO; 86 int w = 0; 87 for (int i = 0; i < n; i++) { 88 tot = dec_add(tot, V[i]); 89 terr = dec_add(terr, Er[i]); 90 if (dec_cmp(Er[i], Er[w]) > 0) w = i; 91 } 92 if (!ok) return num_err(); 93 Dec tol = dec_mul(dec_add(dec_abs(tot), dec_parse("1e-12", 0)), dec_parse("1e-11", 0)); 94 if (dec_cmp(terr, tol) <= 0 || n >= SEGS) { 95 Num r = num_dec(dec_round_sig(tot, 12)); 96 if (r.v.m && r.v.e < dec_abs(tot).e - 11) r = num_int(0); 97 return r; 98 } 99 Dec m = dec_div_int(dec_add(A[w], B[w]), 2); 100 A[n] = m; B[n] = B[w]; 101 B[w] = m; 102 V[w] = gk15(p, from, to, cx, base, A[w], B[w], &Er[w], &ok); 103 V[n] = gk15(p, from, to, cx, base, A[n], B[n], &Er[n], &ok); 104 n++; 105 } 106 } 107 108 /* derivada: diferencia central con extrapolacion de Richardson */ 109 static Num derive(const Prog *p, int from, int to, const EvalCtx *cx, int base, Num at) 110 { 111 bool ok = true; 112 Dec x = at.v; 113 Dec h = dec_mul(dec_add(dec_abs(x), DEC_ONE), dec_parse("0.01", 0)); 114 Dec D[4]; 115 for (int k = 0; k < 4 && ok; k++) { 116 Dec f1 = body_dec(p, from, to, cx, base, dec_add(x, h), &ok); 117 Dec f0 = body_dec(p, from, to, cx, base, dec_sub(x, h), &ok); 118 D[k] = dec_div(dec_sub(f1, f0), dec_add(h, h)); 119 h = dec_div_int(h, 2); 120 } 121 if (!ok) return num_err(); 122 /* Richardson: el error de la diferencia central va como h², h⁴, h⁶ */ 123 for (int j = 1; j < 4; j++) { 124 int64_t f = (int64_t)1 << (2 * j); 125 for (int k = 3; k >= j; k--) D[k] = dec_div_int(dec_sub(dec_mul_int(D[k], f), D[k - 1]), f - 1); 126 } 127 Dec r = dec_round_sig(D[3], 10); 128 if (r.m && r.e < x.e - 12 && r.e < -12) r = DEC_ZERO; 129 return num_dec(r); 130 } 131 132 /* Σ y Π: X entero de lo a hi (exacto si los terminos lo son) */ 133 static Num series(const Prog *p, int from, int to, const EvalCtx *cx, int base, Num lo, Num hi, bool prod) 134 { 135 int64_t a, b; 136 if (!dec_to_int(lo.v, &a) || !dec_to_int(hi.v, &b) || b < a || b - a > 100000) return num_err(); 137 Num acc = num_int(prod ? 1 : 0), t; 138 for (int64_t x = a; x <= b; x++) { 139 if (!body_at(p, from, to, cx, base, num_int(x), &t)) return num_err(); 140 acc = prod ? num_mul(acc, t) : num_add(acc, t); 141 if (num_bad(&acc)) return acc; 142 } 143 return acc; 144 } 145 146 /* modo programador: entero con signo de 'bits' bits */ 147 static int64_t wrap(int64_t v, int bits) 148 { 149 if (bits >= 64) return v; 150 uint64_t m = (1ull << bits) - 1, u = (uint64_t)v & m; 151 if (u >> (bits - 1)) u |= ~m; 152 return (int64_t)u; 153 } 154 155 static Num prog_bin(int op, Num a, Num b, int bits, bool *ok) 156 { 157 int64_t x, y, r = 0; 158 *ok = dec_to_int(a.v, &x) && dec_to_int(b.v, &y); 159 if (!*ok) return num_err(); 160 uint64_t ux = (uint64_t)x, uy = (uint64_t)y; 161 switch (op) { 162 case OP_ADD: r = (int64_t)(ux + uy); break; 163 case OP_SUB: r = (int64_t)(ux - uy); break; 164 case OP_MUL: r = (int64_t)(ux * uy); break; 165 case OP_DIV: 166 if (!y) { *ok = false; return num_err(); } 167 r = x / y; /* hacia cero, como la Casio */ 168 break; 169 case OP_AND: r = x & y; break; 170 case OP_OR: r = x | y; break; 171 case OP_XOR: r = x ^ y; break; 172 case OP_XNOR: r = ~(x ^ y); break; 173 } 174 return num_int(wrap(r, bits)); 175 } 176 177 int eval(const Prog *p, const EvalCtx *cx, Num *out, int *pos) 178 { 179 depth = 0; 180 return run(p, 0, p->n, cx, 0, out, pos); 181 } 182 183 /* una operacion con argumentos reales a[0..nargs-1] */ 184 static Num op_real(const Prog *p, const Code *c, Num *a, const EvalCtx *cx, int sp) 185 { 186 Num r; 187 switch (c->op) { 188 case OP_NUM: r = p->k[c->arg]; break; 189 case OP_VAR: r = cx->vars ? cx->vars[c->arg] : num_int(0); break; 190 case OP_ANS: r = cx->ans; break; 191 case OP_PREANS: r = cx->preans; break; 192 case OP_PI: { Ex x = { 1, 0, 1, 1, 1 }; r = num_ex(&x); break; } 193 case OP_E: r = num_dec(d_e()); break; 194 case OP_RAN: { 195 uint32_t v = cx->rand ? cx->rand() % 1000 : 0; 196 Ex x; 197 ex_rat(&x, v, 1000); 198 r = num_ex(&x); 199 break; 200 } 201 case OP_ADD: case OP_SUB: case OP_MUL: case OP_DIV: 202 case OP_AND: case OP_OR: case OP_XOR: case OP_XNOR: 203 if (cx->prog) { 204 bool ok; 205 r = prog_bin(c->op, a[0], a[1], cx->bits, &ok); 206 if (!ok) r = num_err(); 207 } else if (c->op == OP_ADD) r = num_add(a[0], a[1]); 208 else if (c->op == OP_SUB) r = num_sub(a[0], a[1]); 209 else if (c->op == OP_MUL) r = num_mul(a[0], a[1]); 210 else if (c->op == OP_DIV) r = num_div(a[0], a[1]); 211 else r = num_err(); 212 break; 213 case OP_NEG: 214 if (cx->prog) { int64_t v; dec_to_int(a[0].v, &v); r = num_int(wrap((int64_t)(0 - (uint64_t)v), cx->bits)); } 215 else r = num_neg(a[0]); 216 break; 217 case OP_POW: r = num_pow(a[0], a[1]); break; 218 case OP_SQRT: r = num_sqrt(a[0]); break; 219 case OP_ROOT: r = num_root(a[0], a[1]); break; 220 case OP_LOGB: r = num_logb(a[0], a[1]); break; 221 case OP_ABS: r = num_abs(a[0]); break; 222 case OP_MIXED: { 223 /* a b/c: el signo del entero manda */ 224 Num f = num_div(a[1], a[2]); 225 r = a[0].v.neg ? num_sub(a[0], f) : num_add(a[0], f); 226 break; 227 } 228 case OP_FACT: r = num_fact(a[0]); break; 229 case OP_PCT: r = num_div(a[0], num_int(100)); break; 230 case OP_INV: r = num_div(num_int(1), a[0]); break; 231 case OP_NPR: r = num_npr(a[0], a[1]); break; 232 case OP_NCR: r = num_ncr(a[0], a[1]); break; 233 case OP_EQ: r = num_sub(a[0], a[1]); break; 234 case OP_INTEG: r = integrate(p, c->arg, c->arg + c->len, cx, sp, a[0], a[1]); break; 235 case OP_DERIV: r = derive(p, c->arg, c->arg + c->len, cx, sp, a[0]); break; 236 case OP_SUM: case OP_PROD: r = series(p, c->arg, c->arg + c->len, cx, sp, a[0], a[1], c->op == OP_PROD); break; 237 case OP_FN: { 238 int ang = cx->ang; 239 Num x = a[0]; 240 switch (c->arg) { 241 case FN_SIN: r = num_sin(x, ang); break; 242 case FN_COS: r = num_cos(x, ang); break; 243 case FN_TAN: r = num_tan(x, ang); break; 244 case FN_ASIN: r = num_asin(x, ang); break; 245 case FN_ACOS: r = num_acos(x, ang); break; 246 case FN_ATAN: r = num_atan(x, ang); break; 247 case FN_SINH: r = num_dec(d_sinh(x.v)); break; 248 case FN_COSH: r = num_dec(d_cosh(x.v)); break; 249 case FN_TANH: r = num_dec(d_tanh(x.v)); break; 250 case FN_ASINH: r = num_dec(d_asinh(x.v)); break; 251 case FN_ACOSH: r = num_dec(d_acosh(x.v)); break; 252 case FN_ATANH: r = num_dec(d_atanh(x.v)); break; 253 case FN_LN: r = num_ln(x); break; 254 case FN_LOG: r = num_log10(x); break; 255 case FN_INT: r = num_int_part(x); break; 256 case FN_INTG: r = num_floor(x); break; 257 case FN_RND: r = x.exact ? x : num_dec(dec_round_sig(x.v, 10)); break; 258 case FN_GCD: r = num_gcd(a[0], a[1]); break; 259 case FN_LCM: r = num_lcm(a[0], a[1]); break; 260 case FN_POL: r = num_dec(dec_sqrt(dec_add(dec_mul(a[0].v, a[0].v), dec_mul(a[1].v, a[1].v)))); break; 261 case FN_REC: r = num_mul(a[0], num_cos(a[1], ang)); break; 262 case FN_NOT: { int64_t v; if (!dec_to_int(x.v, &v)) r = num_err(); else r = num_int(wrap(~v, cx->bits ? cx->bits : 64)); break; } 263 case FN_NEG: { int64_t v; if (!dec_to_int(x.v, &v)) r = num_err(); else r = num_int(wrap((int64_t)(0 - (uint64_t)v), cx->bits ? cx->bits : 64)); break; } 264 case FN_NPD: r = num_dec(dist_npd(a[0].v, a[1].v, a[2].v)); break; 265 case FN_NCD: r = num_dec(dist_ncd(a[0].v, a[1].v, a[2].v, a[3].v)); break; 266 case FN_INVN: r = num_dec(dist_invn(a[0].v, a[1].v, a[2].v)); break; 267 case FN_BPD: r = num_dec(dist_bpd(a[0].v, a[1].v, a[2].v)); break; 268 case FN_BCD: r = num_dec(dist_bcd(a[0].v, a[1].v, a[2].v)); break; 269 case FN_PPD: r = num_dec(dist_ppd(a[0].v, a[1].v)); break; 270 case FN_PCD: r = num_dec(dist_pcd(a[0].v, a[1].v)); break; 271 case FN_ARG: r = x.v.neg ? (ang == ANG_RAD ? num_ex(&(Ex){ 1, 0, 1, 1, 1 }) : num_int(ang == ANG_DEG ? 180 : 200)) : num_int(0); break; 272 case FN_CONJ: case FN_RE: r = x; break; 273 case FN_IM: r = num_int(0); break; 274 case FN_XHAT: case FN_YHAT: { 275 extern Num stat_model(const EvalCtx *cx, Num v, bool inverse); 276 r = cx->rvalid ? stat_model(cx, x, c->arg == FN_XHAT) : num_err(); 277 break; 278 } 279 default: r = num_err(); break; 280 } 281 break; 282 } 283 default: r = num_err(); break; 284 } 285 return r; 286 } 287 288 static int run(const Prog *p, int from, int to, const EvalCtx *cx, int base, Num *out, int *pos) 289 { 290 int sp = base; 291 for (int i = from; i < to; i++) { 292 const Code *c = &p->c[i]; 293 int na = c->nargs; 294 if (c->op == OP_SKIP) { i += c->len; continue; } 295 if (sp - base < na) { if (pos) *pos = c->pos; return CE_SYNTAX; } 296 Num r = op_real(p, c, &st[sp - na], cx, sp); 297 if (num_bad(&r)) { if (pos) *pos = c->pos; return CE_MATH; } 298 sp -= na; 299 if (sp >= STACK) { if (pos) *pos = c->pos; return CE_STACK; } 300 st[sp++] = r; 301 } 302 if (sp != base + 1) { if (pos) *pos = 0; return CE_SYNTAX; } 303 *out = st[base]; 304 return CE_OK; 305 } 306 307 /* ------------------------------------------------------------ complejos */ 308 309 static Num sti[STACK]; /* partes imaginarias, paralelas a st */ 310 311 static bool is0(const Num *v) { return !v->v.m; } 312 313 static void c_mul(Num a, Num b, Num c, Num d, Num *re, Num *im) 314 { 315 *re = num_sub(num_mul(a, c), num_mul(b, d)); 316 *im = num_add(num_mul(a, d), num_mul(b, c)); 317 } 318 319 static bool c_div(Num a, Num b, Num c, Num d, Num *re, Num *im) 320 { 321 Num den = num_add(num_mul(c, c), num_mul(d, d)); 322 if (!den.v.m) return false; 323 *re = num_div(num_add(num_mul(a, c), num_mul(b, d)), den); 324 *im = num_div(num_sub(num_mul(b, c), num_mul(a, d)), den); 325 return true; 326 } 327 328 static void c_sqrt(Num a, Num b, Num *re, Num *im) 329 { 330 if (is0(&b)) { 331 if (!a.v.neg) { *re = num_sqrt(a); *im = num_int(0); } 332 else { *re = num_int(0); *im = num_sqrt(num_neg(a)); } 333 return; 334 } 335 Num m = num_sqrt(num_add(num_mul(a, a), num_mul(b, b))); 336 *re = num_sqrt(num_div(num_add(m, a), num_int(2))); 337 *im = num_sqrt(num_div(num_sub(m, a), num_int(2))); 338 if (b.v.neg) *im = num_neg(*im); 339 } 340 341 /* media vuelta en la unidad de angulo */ 342 static Num half_turn(int ang) 343 { 344 if (ang == ANG_RAD) { Ex p = { 1, 0, 1, 1, 1 }; return num_ex(&p); } 345 return num_int(ang == ANG_DEG ? 180 : 200); 346 } 347 348 Num cplx_arg(Num a, Num b, int ang) 349 { 350 if (is0(&a) && is0(&b)) return num_err(); 351 if (is0(&a)) { Num q = num_div(half_turn(ang), num_int(2)); return b.v.neg ? num_neg(q) : q; } 352 Num t = num_atan(num_div(b, a), ang); 353 if (a.v.neg) t = b.v.neg ? num_sub(t, half_turn(ang)) : num_add(t, half_turn(ang)); 354 return t; 355 } 356 357 /* z^w general: exp(w ln z), en decimal */ 358 static bool c_pow_gen(Num a, Num b, Num c, Num d, int ang, Num *re, Num *im) 359 { 360 if (is0(&a) && is0(&b)) return false; 361 Dec lr = d_ln(dec_sqrt(dec_add(dec_mul(a.v, a.v), dec_mul(b.v, b.v)))); 362 Dec th = d_atan2(b.v, a.v, ANG_RAD); 363 /* (c + di)(lr + i th) = (c lr - d th) + i(c th + d lr) */ 364 Dec x = dec_sub(dec_mul(c.v, lr), dec_mul(d.v, th)), y = dec_add(dec_mul(c.v, th), dec_mul(d.v, lr)); 365 Dec ex = d_exp(x); 366 (void)ang; 367 *re = num_dec(dec_mul(ex, d_cos(y, ANG_RAD))); 368 *im = num_dec(dec_mul(ex, d_sin(y, ANG_RAD))); 369 return !re->v.err && !im->v.err; 370 } 371 372 static bool c_pow(Num a, Num b, Num c, Num d, int ang, Num *re, Num *im) 373 { 374 int64_t n; 375 if (is0(&d) && c.exact && ex_is_int(&c.x) && dec_to_int(c.v, &n) && n >= -1000 && n <= 1000) { 376 if (is0(&b)) { 377 /* real: si la base es negativa y el exponente entero, es real */ 378 *re = num_pow(a, c); *im = num_int(0); 379 return !num_bad(re); 380 } 381 bool inv = n < 0; 382 uint64_t k = inv ? (uint64_t)-n : (uint64_t)n; 383 Num rr = num_int(1), ri = num_int(0), br = a, bi = b; 384 while (k) { 385 if (k & 1) c_mul(rr, ri, br, bi, &rr, &ri); 386 k >>= 1; 387 if (k) c_mul(br, bi, br, bi, &br, &bi); 388 } 389 if (inv) return c_div(num_int(1), num_int(0), rr, ri, re, im); 390 *re = rr; *im = ri; 391 return true; 392 } 393 if (is0(&d) && is0(&b) && !a.v.neg) { *re = num_pow(a, c); *im = num_int(0); return !num_bad(re); } 394 /* raiz cuadrada de negativos y similares: c = 1/2 */ 395 if (is0(&d) && c.exact && ex_is_rat(&c.x) && c.x.a == 1 && c.x.d == 2) { c_sqrt(a, b, re, im); return true; } 396 return c_pow_gen(a, b, c, d, ang, re, im); 397 } 398 399 /* redondeo final: 15 cifras y fuera el polvo numerico */ 400 static void tidy(Num *re, Num *im) 401 { 402 if (!re->exact) re->v = dec_round_sig(re->v, 15); 403 if (!im->exact) im->v = dec_round_sig(im->v, 15); 404 if (re->v.m && im->v.m) { 405 if (!im->exact && im->v.e < re->v.e - 14) *im = num_int(0); 406 else if (!re->exact && re->v.e < im->v.e - 14) *re = num_int(0); 407 } 408 } 409 410 static int run_c(const Prog *p, const EvalCtx *cx, Num *ore, Num *oim, int *pos) 411 { 412 int sp = 0; 413 for (int i = 0; i < p->n; i++) { 414 const Code *c = &p->c[i]; 415 int na = c->nargs; 416 if (c->op == OP_SKIP) { i += c->len; continue; } 417 if (sp < na) { if (pos) *pos = c->pos; return CE_SYNTAX; } 418 Num *a = &st[sp - na], *b = &sti[sp - na]; 419 Num re = num_int(0), im = num_int(0); 420 bool ok = true, done = true; 421 switch (c->op) { 422 case OP_IMAG: im = num_int(1); break; 423 case OP_ANS: re = cx->ans; im = cx->ans_im; break; 424 case OP_VAR: re = cx->vars ? cx->vars[c->arg] : num_int(0); im = cx->vars_im ? cx->vars_im[c->arg] : num_int(0); break; 425 case OP_ADD: re = num_add(a[0], a[1]); im = num_add(b[0], b[1]); break; 426 case OP_SUB: case OP_EQ: re = num_sub(a[0], a[1]); im = num_sub(b[0], b[1]); break; 427 case OP_MUL: c_mul(a[0], b[0], a[1], b[1], &re, &im); break; 428 case OP_DIV: ok = c_div(a[0], b[0], a[1], b[1], &re, &im); break; 429 case OP_INV: ok = c_div(num_int(1), num_int(0), a[0], b[0], &re, &im); break; 430 case OP_NEG: re = num_neg(a[0]); im = num_neg(b[0]); break; 431 case OP_PCT: re = num_div(a[0], num_int(100)); im = num_div(b[0], num_int(100)); break; 432 case OP_SQRT: c_sqrt(a[0], b[0], &re, &im); break; 433 case OP_POW: ok = c_pow(a[0], b[0], a[1], b[1], cx->ang, &re, &im); break; 434 case OP_ABS: re = num_sqrt(num_add(num_mul(a[0], a[0]), num_mul(b[0], b[0]))); break; 435 case OP_POLAR: 436 if (!is0(&b[0]) || !is0(&b[1])) { ok = false; break; } 437 re = num_mul(a[0], num_cos(a[1], cx->ang)); 438 im = num_mul(a[0], num_sin(a[1], cx->ang)); 439 break; 440 case OP_FN: 441 switch (c->arg) { 442 case FN_ARG: re = cplx_arg(a[0], b[0], cx->ang); break; 443 case FN_CONJ: re = a[0]; im = num_neg(b[0]); break; 444 case FN_RE: re = a[0]; break; 445 case FN_IM: re = b[0]; break; 446 default: done = false; 447 } 448 break; 449 default: done = false; 450 } 451 if (!done) { 452 /* el resto solo con argumentos reales */ 453 for (int k = 0; k < na; k++) if (!is0(&b[k])) { if (pos) *pos = c->pos; return CE_MATH; } 454 re = op_real(p, c, a, cx, sp); 455 } 456 if (!ok || num_bad(&re) || num_bad(&im)) { if (pos) *pos = c->pos; return CE_MATH; } 457 sp -= na; 458 if (sp >= STACK) { if (pos) *pos = c->pos; return CE_STACK; } 459 st[sp] = re; sti[sp] = im; 460 sp++; 461 } 462 if (sp != 1) { if (pos) *pos = 0; return CE_SYNTAX; } 463 *ore = st[0]; *oim = sti[0]; 464 tidy(ore, oim); 465 return CE_OK; 466 } 467 468 int calc_eval_c(const Expr *e, const EvalCtx *cx, Num *re, Num *im, int *pos) 469 { 470 static Prog p; 471 int err = compile(e, cx, &p, pos); 472 if (err) return err; 473 depth = 0; 474 return run_c(&p, cx, re, im, pos); 475 } 476 477 /* ------------------------------------------------------------ matrices */ 478 479 #define MSTACK 6 480 static MVal msk[MSTACK]; /* static: cada valor ocupa ~1 KB */ 481 482 static bool all_scalar(const MVal *a, int n) 483 { 484 for (int i = 0; i < n; i++) if (a[i].kind != MV_SCALAR) return false; 485 return true; 486 } 487 488 static bool mat_op(const Prog *p, const Code *c, MVal *a, const EvalCtx *cx, int sp, MVal *r) 489 { 490 int na = c->nargs; 491 if (c->op == OP_MAT || c->op == OP_VCT) { 492 const MVal *src = c->op == OP_MAT ? cx->mats : cx->vcts; 493 if (!src || !src[c->arg].r) return false; /* sin definir */ 494 *r = src[c->arg]; 495 return true; 496 } 497 if (c->op == OP_ANS && cx->ans_mat && cx->mats) { 498 *r = cx->mats[4].r ? cx->mats[4] : cx->vcts[4]; 499 return true; 500 } 501 if (all_scalar(a, na) && !(c->op == OP_FN && c->arg == FN_IDEN)) { 502 static Num args[4]; 503 for (int i = 0; i < na && i < 4; i++) args[i] = a[i].a[0]; 504 Num v = op_real(p, c, args, cx, sp); 505 if (num_bad(&v)) return false; 506 mv_scalar(r, v); 507 return true; 508 } 509 Num v; 510 switch (c->op) { 511 case OP_ADD: return mv_add(r, &a[0], &a[1], false); 512 case OP_SUB: return mv_add(r, &a[0], &a[1], true); 513 case OP_MUL: return mv_mul(r, &a[0], &a[1]); 514 case OP_DIV: return a[1].kind == MV_SCALAR && mv_scale(r, &a[0], a[1].a[0], true); 515 case OP_NEG: return mv_neg(r, &a[0]); 516 case OP_INV: return mv_inv(r, &a[0]); 517 case OP_POW: { 518 int64_t n; 519 return a[1].kind == MV_SCALAR && dec_to_int(a[1].a[0].v, &n) && mv_pow(r, &a[0], n); 520 } 521 case OP_ABS: 522 if (a[0].kind == MV_VEC) { if (!mv_norm(&v, &a[0])) return false; mv_scalar(r, v); return true; } 523 return mv_abs(r, &a[0]); 524 case OP_FN: 525 switch (c->arg) { 526 case FN_DET: if (!mv_det(&v, &a[0])) return false; mv_scalar(r, v); return true; 527 case FN_TRN: return mv_trn(r, &a[0]); 528 case FN_IDEN: { 529 int64_t n; 530 return a[0].kind == MV_SCALAR && dec_to_int(a[0].a[0].v, &n) && mv_ident(r, (int)n); 531 } 532 case FN_DOT: if (!mv_dot(&v, &a[0], &a[1])) return false; mv_scalar(r, v); return true; 533 case FN_CROSS: return mv_cross(r, &a[0], &a[1]); 534 case FN_VANG: if (!mv_angle(&v, &a[0], &a[1], cx->ang)) return false; mv_scalar(r, v); return true; 535 case FN_UNITV: return mv_unit(r, &a[0]); 536 } 537 return false; 538 } 539 return false; 540 } 541 542 int calc_eval_m(const Expr *e, const EvalCtx *cx, MVal *out, int *pos) 543 { 544 static Prog p; 545 int err = compile(e, cx, &p, pos); 546 if (err) return err; 547 depth = 0; 548 int sp = 0; 549 for (int i = 0; i < p.n; i++) { 550 const Code *c = &p.c[i]; 551 int na = c->nargs; 552 if (c->op == OP_SKIP) { i += c->len; continue; } 553 if (sp < na) { if (pos) *pos = c->pos; return CE_SYNTAX; } 554 static MVal r; 555 if (!mat_op(&p, c, &msk[sp - na], cx, sp, &r)) { if (pos) *pos = c->pos; return CE_MATH; } 556 sp -= na; 557 if (sp >= MSTACK) { if (pos) *pos = c->pos; return CE_STACK; } 558 msk[sp++] = r; 559 } 560 if (sp != 1) { if (pos) *pos = 0; return CE_SYNTAX; } 561 *out = msk[0]; 562 if (out->kind == MV_SCALAR && !out->a[0].exact) out->a[0].v = dec_round_sig(out->a[0].v, 15); 563 return CE_OK; 564 } 565 566 int calc_eval(const Expr *e, const EvalCtx *cx, Num *out, int *pos) 567 { 568 static Prog p; 569 int err = compile(e, cx, &p, pos); 570 if (err) return err; 571 return eval(&p, cx, out, pos); 572 }