modes.c (36468B)
1 /* modes.c - modos Ecuaciones, Tabla, Grafico y Estadistica (Programador usa la 2 * pantalla de Calcular con otro contexto de evaluacion). */ 3 #include <string.h> 4 #include "calc.h" 5 #include "dmath.h" 6 #include "str.h" 7 8 bool edit_key(Calc *c, Expr *e, int key); 9 10 static void say(Calc *c, const char *m) { str_put(c->msg, sizeof c->msg, 0, m); } 11 12 /* evalua la entrada (c->in) con los ajustes actuales */ 13 static bool eval_in(Calc *c, Num *out) 14 { 15 EvalCtx cx = calc_ctx(c); 16 int pos; 17 ex_tidy(&c->in); 18 c->err = calc_eval(&c->in, &cx, out, &pos); 19 if (c->err) { say(c, calc_err_text(c->err)); c->in.cur = pos < c->in.n ? pos : c->in.n; return false; } 20 if (!out->exact) out->v = dec_round_sig(out->v, 15); 21 return true; 22 } 23 24 /* evalua e con X = x (para tabla, grafico y SOLVE) */ 25 static bool eval_at(Calc *c, const Expr *e, Num x, Num *out) 26 { 27 EvalCtx cx = calc_ctx(c); 28 static Num vars[NVARS]; 29 memcpy(vars, c->vars, sizeof vars); 30 vars[VAR_X] = x; 31 cx.vars = vars; 32 int pos; 33 return calc_eval(e, &cx, out, &pos) == CE_OK; 34 } 35 36 /* ============================================================ ecuaciones */ 37 38 static int eq_cols(int type) { return type == EQ_SYS2 ? 3 : type == EQ_SYS3 ? 4 : type == EQ_POLY2 ? 3 : 4; } 39 static int eq_rows(int type) { return type == EQ_SYS2 ? 2 : type == EQ_SYS3 ? 3 : 1; } 40 int eq_cells(int type) { return eq_cols(type) * eq_rows(type); } 41 42 static Cplx cx_real(Num v) { Cplx z = { v, num_int(0) }; return z; } 43 44 /* Gauss con Num: exacto mientras se pueda */ 45 static void solve_linear(Calc *c) 46 { 47 EqnState *q = &c->eq; 48 int n = q->type == EQ_SYS2 ? 2 : 3, m = n + 1; 49 static Num a[3][4]; 50 for (int i = 0; i < n; i++) 51 for (int j = 0; j < m; j++) a[i][j] = q->coef[i * m + j]; 52 for (int col = 0; col < n; col++) { 53 int p = -1; 54 for (int r = col; r < n; r++) if (a[r][col].v.m) { p = r; break; } 55 if (p < 0) { 56 /* singular: ver si es incompatible o indeterminado */ 57 str_put(q->note, sizeof q->note, 0, "Infinitas soluciones o ninguna"); 58 q->nsol = 0; 59 return; 60 } 61 if (p != col) for (int j = 0; j < m; j++) { Num t = a[p][j]; a[p][j] = a[col][j]; a[col][j] = t; } 62 for (int r = 0; r < n; r++) { 63 if (r == col || !a[r][col].v.m) continue; 64 Num f = num_div(a[r][col], a[col][col]); 65 for (int j = col; j < m; j++) a[r][j] = num_sub(a[r][j], num_mul(f, a[col][j])); 66 } 67 } 68 for (int i = 0; i < n; i++) { 69 Num v = num_div(a[i][n], a[i][i]); 70 if (!v.exact) v.v = dec_round_sig(v.v, 15); 71 if (!v.exact && v.v.m && v.v.e < -15) v = num_int(0); 72 q->sol[i] = cx_real(v); 73 } 74 q->nsol = n; 75 q->note[0] = 0; 76 } 77 78 static Num half_num(Num v) { return num_div(v, num_int(2)); } 79 80 /* raices de a x² + b x + c (a != 0) */ 81 static int quad(Num a, Num b, Num c, Cplx *out) 82 { 83 Num D = num_sub(num_mul(b, b), num_mul(num_int(4), num_mul(a, c))); 84 Num two_a = num_mul(num_int(2), a); 85 Num mb = num_neg(b); 86 if (!D.v.neg) { 87 Num s = num_sqrt(D); 88 out[0] = cx_real(num_div(num_add(mb, s), two_a)); 89 out[1] = cx_real(num_div(num_sub(mb, s), two_a)); 90 if (!D.v.m) return 1; 91 return 2; 92 } 93 Num re = num_div(mb, two_a), im = num_abs(num_div(num_sqrt(num_neg(D)), two_a)); 94 out[0].re = re; out[0].im = im; 95 out[1].re = re; out[1].im = num_neg(im); 96 return 2; 97 } 98 99 static Num poly3(const Num *k, Num x) /* a x³ + b x² + c x + d (Horner) */ 100 { 101 Num v = k[0]; 102 for (int i = 1; i < 4; i++) v = num_add(num_mul(v, x), k[i]); 103 return v; 104 } 105 106 /* divisores positivos de |v| (hasta 64, con v chico) */ 107 static int divisors(int64_t v, int64_t *out, int max) 108 { 109 if (v < 0) v = -v; 110 int n = 0; 111 if (!v || v > 1000000000000ll) return 0; 112 for (int64_t d = 1; d * d <= v && n < max - 1; d++) 113 if (v % d == 0) { out[n++] = d; if (d != v / d) out[n++] = v / d; } 114 return n; 115 } 116 117 static void solve_cubic(Calc *c) 118 { 119 EqnState *q = &c->eq; 120 Num *k = q->coef; 121 Num a = k[0]; 122 /* 1) raiz racional si los coeficientes son enteros */ 123 bool ints = true; 124 for (int i = 0; i < 4; i++) if (!(k[i].exact && ex_is_int(&k[i].x))) ints = false; 125 Num root; 126 bool found = false; 127 if (ints && k[3].x.a == 0) { root = num_int(0); found = true; } 128 if (ints && !found) { 129 static int64_t P[64], Q[64]; 130 int np = divisors(k[3].x.a, P, 64), nq = divisors(k[0].x.a, Q, 64); 131 for (int i = 0; i < np && !found; i++) 132 for (int j = 0; j < nq && !found; j++) 133 for (int s = -1; s <= 1 && !found; s += 2) { 134 Ex r; 135 if (!ex_rat(&r, s * P[i], Q[j])) continue; 136 Num x = num_ex(&r); 137 Num v = poly3(k, x); 138 if (v.exact && ex_is_zero(&v.x)) { root = x; found = true; } 139 } 140 } 141 /* 2) si no, una raiz real por biseccion en decimal */ 142 static Num kk[4]; /* static: la pila de la Pico es chica */ 143 if (!found) { 144 for (int i = 0; i < 4; i++) kk[i] = num_to_dec(k[i]); 145 Num A = kk[0]; 146 /* cota de Cauchy */ 147 Dec bound = DEC_ONE; 148 for (int i = 1; i < 4; i++) { 149 Dec r = dec_abs(dec_div(kk[i].v, A.v)); 150 if (dec_cmp(r, bound) > 0) bound = r; 151 } 152 bound = dec_add(bound, DEC_ONE); 153 Num lo = num_dec(dec_neg(bound)), hi = num_dec(bound); 154 bool lo_neg = poly3(kk, lo).v.neg; 155 for (int it = 0; it < 200; it++) { 156 Num mid = num_dec(dec_div_int(dec_add(lo.v, hi.v), 2)); 157 Num v = poly3(kk, mid); 158 if (!v.v.m) { lo = hi = mid; break; } 159 if (v.v.neg == lo_neg) lo = mid; else hi = mid; 160 if (dec_cmp(dec_abs(dec_sub(hi.v, lo.v)), dec_mul(dec_parse("1e-17", 0), dec_add(dec_abs(lo.v), DEC_ONE))) < 0) break; 161 } 162 root = num_dec(dec_round_sig(lo.v, 15)); 163 k = kk; 164 a = A; 165 } 166 /* deflacion: a x³ + b x² + c x + d = (x - r)(a x² + b' x + c') */ 167 Num b2 = num_add(k[1], num_mul(a, root)); 168 Num c2 = num_add(k[2], num_mul(b2, root)); 169 q->sol[0] = cx_real(root); 170 Cplx r2[2]; 171 int n2 = quad(a, b2, c2, r2); 172 q->sol[1] = r2[0]; 173 q->sol[2] = n2 > 1 ? r2[1] : r2[0]; 174 q->nsol = 3; 175 /* sacar el polvo numerico: partes imaginarias o reales casi cero */ 176 for (int i = 0; i < 3; i++) { 177 if (!q->sol[i].im.exact && q->sol[i].im.v.m && q->sol[i].im.v.e < -13) q->sol[i].im = num_int(0); 178 if (!q->sol[i].re.exact) q->sol[i].re.v = dec_round_sig(q->sol[i].re.v, 15); 179 if (!q->sol[i].im.exact) q->sol[i].im.v = dec_round_sig(q->sol[i].im.v, 15); 180 } 181 } 182 183 /* SOLVE: raiz de L - R cerca de X (Newton con derivada numerica, despues secante) */ 184 static void solve_eq(Calc *c) 185 { 186 EqnState *q = &c->eq; 187 const Expr *e = &q->solve_eq; 188 Dec x = num_to_dec(c->vars[VAR_X]).v, fx; 189 Num r; 190 bool ok = false; 191 for (int it = 0; it < 80; it++) { 192 if (!eval_at(c, e, num_dec(x), &r)) break; 193 fx = r.v; 194 if (!fx.m) { ok = true; break; } 195 Dec h = dec_mul(dec_add(dec_abs(x), DEC_ONE), dec_parse("1e-8", 0)); 196 Num r2; 197 if (!eval_at(c, e, num_dec(dec_add(x, h)), &r2)) break; 198 Dec d = dec_div(dec_sub(r2.v, fx), h); 199 if (!d.m || d.err) { x = dec_add(x, dec_mul(h, dec_int(1000))); continue; } 200 Dec step = dec_div(fx, d); 201 x = dec_sub(x, step); 202 if (step.m == 0 || dec_cmp_abs(step, dec_mul(dec_add(dec_abs(x), DEC_ONE), dec_parse("1e-17", 0))) < 0) { ok = true; break; } 203 } 204 if (ok && eval_at(c, e, num_dec(x), &r)) { 205 q->solve_x = num_dec(dec_round_sig(x, 15)); 206 q->solve_lr = num_dec(dec_round_sig(r.v, 15)); 207 c->vars[VAR_X] = q->solve_x; 208 q->nsol = 1; 209 q->note[0] = 0; 210 } else { 211 q->nsol = 0; 212 str_put(q->note, sizeof q->note, 0, "No se encontró solución"); 213 } 214 } 215 216 static void eq_solve(Calc *c) 217 { 218 EqnState *q = &c->eq; 219 q->note[0] = 0; 220 q->shown = 0; 221 switch (q->type) { 222 case EQ_SYS2: case EQ_SYS3: solve_linear(c); break; 223 case EQ_POLY2: { 224 if (!q->coef[0].v.m) { str_put(q->note, sizeof q->note, 0, "a no puede ser 0"); q->nsol = 0; break; } 225 Cplx r[2]; 226 int n = quad(q->coef[0], q->coef[1], q->coef[2], r); 227 q->sol[0] = r[0]; q->sol[1] = r[1]; 228 q->nsol = n; 229 /* vertice: x = -b/2a, y = f(x) */ 230 Num vx = half_num(num_div(num_neg(q->coef[1]), q->coef[0])); 231 q->sol[2].re = vx; 232 q->sol[2].im = num_add(num_mul(num_add(num_mul(q->coef[0], vx), q->coef[1]), vx), q->coef[2]); 233 break; 234 } 235 case EQ_POLY3: 236 if (!q->coef[0].v.m) { str_put(q->note, sizeof q->note, 0, "a no puede ser 0"); q->nsol = 0; break; } 237 solve_cubic(c); 238 break; 239 case EQ_SOLVE: solve_eq(c); break; 240 } 241 q->phase = 1; 242 q->solved = true; 243 } 244 245 static void eq_type_form(Calc *c) 246 { 247 Form *f = &c->form; 248 form_begin(f, "Ecuaciones"); 249 form_button(f, 100 + EQ_SYS2, "1 Sistema de 2 ecuaciones (x, y)"); 250 form_button(f, 100 + EQ_SYS3, "2 Sistema de 3 ecuaciones (x, y, z)"); 251 form_button(f, 100 + EQ_POLY2, "3 Polinomio de grado 2"); 252 form_button(f, 100 + EQ_POLY3, "4 Polinomio de grado 3"); 253 form_button(f, 100 + EQ_SOLVE, "5 Resolver f(X) = g(X) (SOLVE)"); 254 f->cur = c->eq.type; 255 c->form_kind = FORM_EQN_TYPE; 256 c->scr = SCR_FORM; 257 } 258 259 static bool eqn_key(Calc *c, int key) 260 { 261 EqnState *q = &c->eq; 262 if (q->type == EQ_SOLVE && q->phase == 0) { 263 if (key == K_OK) { 264 if (ex_empty(&c->in)) return true; 265 ex_tidy(&c->in); 266 q->solve_eq = c->in; 267 eq_solve(c); 268 return true; 269 } 270 if (key == K_BACK) { if (ex_empty(&c->in)) eq_type_form(c); else ex_clear(&c->in); return true; } 271 return edit_key(c, &c->in, key); 272 } 273 if (q->phase == 1) { 274 int n = q->type == EQ_POLY2 ? 3 : q->nsol; 275 if (key == K_DOWN || key == K_OK || key == K_RIGHT) { if (n) q->shown = (q->shown + 1) % n; } 276 else if (key == K_UP || key == K_LEFT) { if (n) q->shown = (q->shown + n - 1) % n; } 277 else if (key == K_BACK) { 278 q->phase = 0; 279 if (q->type == EQ_SOLVE) { c->in = q->solve_eq; c->in.cur = c->in.n; } 280 } else if (key == K_TAB) { 281 /* S⇔D sobre las soluciones */ 282 for (int i = 0; i < 3; i++) { 283 q->sol[i].re = num_to_dec(q->sol[i].re); 284 q->sol[i].im = num_to_dec(q->sol[i].im); 285 } 286 } 287 return true; 288 } 289 int cells = eq_cells(q->type), cols = eq_cols(q->type); 290 if (q->editing) { 291 if (key == K_OK) { 292 Num v; 293 if (!eval_in(c, &v)) return true; 294 q->coef[q->sel] = v; 295 q->editing = false; 296 ex_clear(&c->in); 297 q->sel = (q->sel + 1) % cells; 298 return true; 299 } 300 if (key == K_BACK) { q->editing = false; ex_clear(&c->in); return true; } 301 if ((key == K_UP || key == K_DOWN) && ex_empty(&c->in)) { q->editing = false; } 302 else return edit_key(c, &c->in, key); 303 } 304 switch (key) { 305 case K_LEFT: q->sel = (q->sel + cells - 1) % cells; return true; 306 case K_RIGHT: q->sel = (q->sel + 1) % cells; return true; 307 case K_UP: q->sel = (q->sel + cells - cols) % cells; return true; 308 case K_DOWN: q->sel = (q->sel + cols) % cells; return true; 309 case K_OK: eq_solve(c); return true; 310 case K_BACK: eq_type_form(c); return true; 311 case K_BKSP: case K_DEL: q->coef[q->sel] = num_int(0); return true; 312 } 313 if (key >= 32 && key < 127) { 314 q->editing = true; 315 ex_clear(&c->in); 316 return edit_key(c, &c->in, key); 317 } 318 return false; 319 } 320 321 /* ============================================================ tabla */ 322 323 static void table_form(Calc *c); 324 325 void table_compute(Calc *c) 326 { 327 TableState *t = &c->tb; 328 t->rows = 0; 329 Num x = t->start; 330 if (!t->step.v.m || t->step.v.neg != dec_sub(t->end.v, t->start.v).neg) { 331 if (dec_cmp(t->end.v, t->start.v) != 0) { say(c, "Paso inválido"); return; } 332 } 333 for (int i = 0; i < TABLE_ROWS; i++) { 334 if (t->step.v.neg ? num_cmp(x, t->end) < 0 : num_cmp(x, t->end) > 0) break; 335 t->x[i] = x; 336 Num v; 337 t->ferr[i] = !eval_at(c, &t->f, x, &v); 338 t->fx[i] = v; 339 if (t->use_g) { t->gerr[i] = !eval_at(c, &t->g, x, &v); t->gx[i] = v; } 340 t->rows++; 341 x = num_add(x, t->step); 342 if (!t->step.v.m) break; 343 } 344 if (t->rows == TABLE_ROWS && num_cmp(x, t->end) <= 0) say(c, "Tabla cortada a 30 filas"); 345 t->top = t->sel = 0; 346 t->ready = true; 347 } 348 349 static void num_field(Calc *c, Form *f, int id, const char *label, Num v) 350 { 351 char b[64]; 352 static Expr e; 353 calc_result_expr(c, &v, true, &e); 354 expr_to_text(&e, b, sizeof b); 355 form_text(f, id, label, b, 30); 356 } 357 358 static void table_form(Calc *c) 359 { 360 Form *f = &c->form; 361 static const char *const NY[] = { "No", "Sí" }; 362 form_begin(f, "Rango de la tabla"); 363 num_field(c, f, 1, "Inicio", c->tb.start); 364 num_field(c, f, 2, "Fin", c->tb.end); 365 num_field(c, f, 3, "Paso", c->tb.step); 366 form_choice(f, 4, "Usar g(x)", NY, 2, c->tb.use_g); 367 form_button(f, 5, "Armar la tabla"); 368 f->cur = 4; 369 c->form_kind = FORM_TABLE; 370 c->scr = SCR_FORM; 371 } 372 373 static bool parse_num(Calc *c, const char *s, Num *out) 374 { 375 static Expr e; 376 if (!expr_from_text(&e, s)) return false; 377 EvalCtx cx = calc_ctx(c); 378 int pos; 379 return calc_eval(&e, &cx, out, &pos) == CE_OK; 380 } 381 382 static bool table_key(Calc *c, int key) 383 { 384 TableState *t = &c->tb; 385 if (t->phase < 2) { 386 Expr *e = t->phase == 0 ? &t->f : &t->g; 387 if (key == K_OK) { 388 *e = c->in; 389 ex_tidy(e); 390 if (t->phase == 0 && t->use_g) { t->phase = 1; c->in = t->g; c->in.cur = c->in.n; } 391 else table_form(c); 392 return true; 393 } 394 if (key == K_DOWN && t->phase == 0 && t->use_g) { t->f = c->in; t->phase = 1; c->in = t->g; c->in.cur = c->in.n; return true; } 395 if (key == K_UP && t->phase == 1) { t->g = c->in; t->phase = 0; c->in = t->f; c->in.cur = c->in.n; return true; } 396 if (key == K_BACK) { ex_clear(&c->in); return true; } 397 return edit_key(c, &c->in, key); 398 } 399 switch (key) { 400 case K_UP: if (t->sel > 0) t->sel--; return true; 401 case K_DOWN: if (t->sel + 1 < t->rows) t->sel++; return true; 402 case K_BACK: case K_OK: t->phase = 0; c->in = t->f; c->in.cur = c->in.n; return true; 403 } 404 return false; 405 } 406 407 /* ============================================================ grafico */ 408 409 static Dec lerp(Dec a, Dec b, int i, int n) { return dec_add(a, dec_div_int(dec_mul_int(dec_sub(b, a), i), n)); } 410 411 /* Dos pasadas (sin guardar los 640 valores: la RAM de la Pico es poca): la primera 412 * busca el rango de y si es automatico, la segunda pasa cada valor a fila de pantalla. */ 413 void graph_compute(Calc *c) 414 { 415 GraphState *g = &c->gr; 416 int W = g->w > 0 ? g->w : 320, H = g->h > 0 ? g->h : 260; 417 if (W > GRAPH_W) W = GRAPH_W; 418 static Prog pf; 419 static Num vars[NVARS]; 420 EvalCtx cx = calc_ctx(c); 421 memcpy(vars, c->vars, sizeof vars); 422 cx.vars = vars; 423 Dec lo = DEC_ZERO, hi = DEC_ZERO, yr = DEC_ONE; 424 bool any = false; 425 for (int pass = g->autoy ? 0 : 1; pass < 2; pass++) { 426 if (pass == 1) { 427 if (g->autoy) { 428 if (!any) { lo = dec_int(-1); hi = dec_int(1); } 429 Dec span = dec_sub(hi, lo); 430 if (!span.m) span = DEC_ONE; 431 Dec m = dec_div_int(span, 10); 432 g->ymin = dec_sub(lo, m); 433 g->ymax = dec_add(hi, m); 434 } 435 yr = dec_sub(g->ymax, g->ymin); 436 } 437 for (int k = 0; k < 2; k++) { 438 const Expr *e = k ? &g->g : &g->f; 439 int16_t *dst = k ? g->gy : g->fy; 440 int pos; 441 bool use = (k == 0 || g->use_g) && !ex_empty(e) && compile(e, &cx, &pf, &pos) == CE_OK; 442 for (int i = 0; i < W; i++) { 443 if (pass == 1) dst[i] = INT16_MIN; 444 if (!use) continue; 445 vars[VAR_X] = num_dec(lerp(g->xmin, g->xmax, i, W - 1)); 446 Num r; 447 if (eval(&pf, &cx, &r, &pos) != CE_OK || r.v.err) continue; 448 Dec y = r.v; 449 if (pass == 0) { 450 if (y.e > 6) continue; /* asintotas: no estirar por ellas */ 451 if (!any || dec_cmp(y, lo) < 0) lo = y; 452 if (!any || dec_cmp(y, hi) > 0) hi = y; 453 any = true; 454 continue; 455 } 456 if (!yr.m) continue; 457 Dec rr = dec_div(dec_mul_int(dec_sub(g->ymax, y), H - 1), yr); 458 int64_t row; 459 if (dec_to_int(dec_round_int(rr), &row) && row > -30000 && row < 30000) dst[i] = (int16_t)row; 460 else dst[i] = rr.neg ? -30000 : 30000; 461 } 462 } 463 } 464 if (g->tx < 0 || g->tx >= W) g->tx = W / 2; 465 g->ready = true; 466 } 467 468 /* valor de la curva en la columna del trazo */ 469 static void trace_value(Calc *c) 470 { 471 GraphState *g = &c->gr; 472 int W = g->w > 0 ? g->w : 320; 473 Num x = num_dec(lerp(g->xmin, g->xmax, g->tx, W - 1)), r; 474 if (!eval_at(c, g->curve ? &g->g : &g->f, x, &r)) r = num_err(); 475 g->tval = r; 476 } 477 478 static void graph_form(Calc *c) 479 { 480 Form *f = &c->form; 481 static const char *const NY[] = { "No", "Sí" }; 482 form_begin(f, "Ventana del gráfico"); 483 num_field(c, f, 1, "x mín", num_dec(c->gr.xmin)); 484 num_field(c, f, 2, "x máx", num_dec(c->gr.xmax)); 485 form_choice(f, 5, "y automática", NY, 2, c->gr.autoy); 486 num_field(c, f, 3, "y mín", num_dec(c->gr.ymin)); 487 num_field(c, f, 4, "y máx", num_dec(c->gr.ymax)); 488 form_choice(f, 6, "Usar g(x)", NY, 2, c->gr.use_g); 489 form_button(f, 7, "Graficar"); 490 f->cur = 6; 491 c->form_kind = FORM_GRAPH; 492 c->scr = SCR_FORM; 493 } 494 495 static void zoom(Calc *c, int num, int den) 496 { 497 GraphState *g = &c->gr; 498 int W = g->w > 0 ? g->w : 320; 499 Dec cxv = lerp(g->xmin, g->xmax, g->tx, W - 1); 500 Dec hw = dec_div_int(dec_mul_int(dec_sub(g->xmax, g->xmin), num), 2 * den); 501 g->xmin = dec_sub(cxv, hw); 502 g->xmax = dec_add(cxv, hw); 503 if (!g->autoy) { 504 Dec cy = dec_div_int(dec_add(g->ymin, g->ymax), 2); 505 Dec hh = dec_div_int(dec_mul_int(dec_sub(g->ymax, g->ymin), num), 2 * den); 506 g->ymin = dec_sub(cy, hh); 507 g->ymax = dec_add(cy, hh); 508 } 509 g->tx = W / 2; 510 g->ready = false; 511 } 512 513 static bool graph_key(Calc *c, int key) 514 { 515 GraphState *g = &c->gr; 516 if (g->phase < 2) { 517 Expr *e = g->phase == 0 ? &g->f : &g->g; 518 if (key == K_OK) { 519 *e = c->in; 520 ex_tidy(e); 521 if (g->phase == 0 && g->use_g) { g->phase = 1; c->in = g->g; c->in.cur = c->in.n; } 522 else { g->phase = 2; g->ready = false; g->curve = 0; g->tx = -1; } 523 return true; 524 } 525 if (key == K_DOWN && g->phase == 0) { g->f = c->in; g->use_g = true; g->phase = 1; c->in = g->g; c->in.cur = c->in.n; return true; } 526 if (key == K_UP && g->phase == 1) { g->g = c->in; g->phase = 0; c->in = g->f; c->in.cur = c->in.n; return true; } 527 if (key == K_BACK) { ex_clear(&c->in); if (g->phase == 1) { g->g = c->in; g->use_g = false; } return true; } 528 return edit_key(c, &c->in, key); 529 } 530 int W = g->w > 0 ? g->w : 320; 531 switch (key) { 532 case K_LEFT: g->tx = g->tx > 0 ? g->tx - 1 : 0; break; 533 case K_RIGHT: g->tx = g->tx + 1 < W ? g->tx + 1 : W - 1; break; 534 case ',': case '<': g->tx = g->tx > 10 ? g->tx - 10 : 0; break; 535 case '.': g->tx = g->tx + 10 < W ? g->tx + 10 : W - 1; break; 536 case K_UP: case K_DOWN: if (g->use_g) g->curve ^= 1; break; 537 case '+': case K_PLUS: zoom(c, 1, 2); break; 538 case '-': case K_MINUS: zoom(c, 2, 1); break; 539 case 'w': case 'W': graph_form(c); return true; 540 case K_OK: case K_BACK: g->phase = 0; c->in = g->f; c->in.cur = c->in.n; return true; 541 default: return false; 542 } 543 trace_value(c); 544 return true; 545 } 546 547 /* ============================================================ estadistica */ 548 549 static const char *const ST_NAMES[ST_NTYPES] = { 550 "1 variable", "y = a + bx", "y = a + bx + cx²", "y = a + b·ln x", "y = a·e^(bx)", "y = a·b^x", 551 "y = a·x^b", "y = a + b/x", 552 }; 553 const char *stat_type_name(int t) { return t >= 0 && t < ST_NTYPES ? ST_NAMES[t] : "?"; } 554 555 static void stat_form(Calc *c) 556 { 557 Form *f = &c->form; 558 form_begin(f, "Estadística: tipo"); 559 for (int i = 0; i < ST_NTYPES; i++) { 560 static char lab[ST_NTYPES][40]; 561 int p = str_int(lab[i], 40, 0, i + 1); 562 p = str_put(lab[i], 40, p, " "); 563 str_put(lab[i], 40, p, ST_NAMES[i]); 564 form_button(f, 200 + i, lab[i]); 565 } 566 f->cur = c->st.type; 567 c->form_kind = FORM_STAT_TYPE; 568 c->scr = SCR_FORM; 569 } 570 571 /* modelo de regresion: y(x), o x(y) si inverse */ 572 Num stat_model(const EvalCtx *cx, Num v, bool inverse) 573 { 574 Num a = cx->ra, b = cx->rb, cc = cx->rc; 575 switch (cx->rtype) { 576 case ST_LIN: return inverse ? num_div(num_sub(v, a), b) : num_add(a, num_mul(b, v)); 577 case ST_QUAD: 578 if (!inverse) return num_add(num_add(a, num_mul(b, v)), num_mul(cc, num_mul(v, v))); 579 { 580 Cplx r[2]; 581 if (quad(cc, b, num_sub(a, v), r) && !r[0].im.v.m) return r[0].re; 582 return num_err(); 583 } 584 case ST_LOG: return inverse ? num_dec(d_exp(dec_div(dec_sub(v.v, a.v), b.v))) : num_add(a, num_mul(b, num_ln(v))); 585 case ST_EXP: return inverse ? num_dec(dec_div(d_ln(dec_div(v.v, a.v)), b.v)) : num_mul(a, num_dec(d_exp(dec_mul(b.v, v.v)))); 586 case ST_ABEXP: return inverse ? num_dec(dec_div(d_ln(dec_div(v.v, a.v)), d_ln(b.v))) : num_mul(a, num_pow(b, v)); 587 case ST_POW: return inverse ? num_dec(d_pow(dec_div(v.v, a.v), dec_div(DEC_ONE, b.v))) : num_mul(a, num_pow(v, b)); 588 case ST_INV: return inverse ? num_div(b, num_sub(v, a)) : num_add(a, num_div(b, v)); 589 } 590 return num_err(); 591 } 592 593 static void res(StatState *s, const char *label, Num v) 594 { 595 if (s->nres >= STAT_RES) return; 596 if (!v.exact) v.v = dec_round_sig(v.v, 15); 597 s->rlabel[s->nres] = label; 598 s->rval[s->nres] = v; 599 s->rerr[s->nres] = num_bad(&v); 600 s->nres++; 601 } 602 603 /* ordena copias de x (insercion: son pocos datos) */ 604 static void sorted(const StatState *s, Num *o) 605 { 606 for (int i = 0; i < s->n; i++) { 607 Num v = s->x[i]; 608 int j = i; 609 while (j > 0 && num_cmp(o[j - 1], v) > 0) { o[j] = o[j - 1]; j--; } 610 o[j] = v; 611 } 612 } 613 614 static Num median(const Num *v, int from, int to) /* [from, to) */ 615 { 616 int n = to - from; 617 if (n <= 0) return num_err(); 618 if (n & 1) return v[from + n / 2]; 619 return half_num(num_add(v[from + n / 2 - 1], v[from + n / 2])); 620 } 621 622 static void stat_compute(Calc *c) 623 { 624 StatState *s = &c->st; 625 s->nres = 0; 626 s->rvalid = false; 627 int n = s->n; 628 if (!n) { say(c, "No hay datos"); return; } 629 static Num N, sx, sx2, sy, sy2, sxy, mx, my, vx, vy, a, b, cc, r; /* static: pila chica */ 630 N = num_int(n); sx = sx2 = sy = sy2 = sxy = num_int(0); 631 for (int i = 0; i < n; i++) { 632 sx = num_add(sx, s->x[i]); 633 sx2 = num_add(sx2, num_mul(s->x[i], s->x[i])); 634 sy = num_add(sy, s->y[i]); 635 sy2 = num_add(sy2, num_mul(s->y[i], s->y[i])); 636 sxy = num_add(sxy, num_mul(s->x[i], s->y[i])); 637 } 638 mx = num_div(sx, N); my = num_div(sy, N); 639 /* sigma² = Σx²/n - x̄² */ 640 vx = num_sub(num_div(sx2, N), num_mul(mx, mx)); 641 vy = num_sub(num_div(sy2, N), num_mul(my, my)); 642 res(s, "n", N); 643 res(s, "x̄", mx); 644 res(s, "Σx", sx); 645 res(s, "Σx²", sx2); 646 res(s, "σx", num_sqrt(vx)); 647 res(s, "sx", n > 1 ? num_sqrt(num_div(num_mul(vx, N), num_int(n - 1))) : num_err()); 648 static Num o[STAT_MAX]; 649 sorted(s, o); 650 res(s, "mín", o[0]); 651 res(s, "Q1", median(o, 0, n / 2)); 652 res(s, "Med", median(o, 0, n)); 653 res(s, "Q3", median(o, (n + 1) / 2, n)); 654 res(s, "máx", o[n - 1]); 655 if (s->type == ST_1VAR) return; 656 res(s, "y\xcc\x84", my); 657 res(s, "Σy", sy); 658 res(s, "Σy²", sy2); 659 res(s, "Σxy", sxy); 660 res(s, "σy", num_sqrt(vy)); 661 /* regresion: se linealiza X o Y segun el modelo y se ajusta Y = A + B X */ 662 int t = s->type; 663 cc = num_int(0); r = num_err(); 664 if (t == ST_QUAD) { 665 /* minimos cuadrados con 3 incognitas: ecuaciones normales */ 666 static Num S[5], T[3]; 667 for (int k = 0; k < 5; k++) S[k] = num_int(0); 668 for (int k = 0; k < 3; k++) T[k] = num_int(0); 669 for (int i = 0; i < n; i++) { 670 Num p = num_int(1); 671 for (int k = 0; k < 5; k++) { 672 S[k] = num_add(S[k], p); 673 if (k < 3) T[k] = num_add(T[k], num_mul(p, s->y[i])); 674 p = num_mul(p, s->x[i]); 675 } 676 } 677 static EqnState save; 678 save = c->eq; 679 c->eq.type = EQ_SYS3; 680 for (int i = 0; i < 3; i++) { 681 for (int j = 0; j < 3; j++) c->eq.coef[i * 4 + j] = S[i + j]; 682 c->eq.coef[i * 4 + 3] = T[i]; 683 } 684 solve_linear(c); 685 bool okq = c->eq.nsol == 3; 686 a = c->eq.sol[0].re; b = c->eq.sol[1].re; cc = c->eq.sol[2].re; 687 c->eq = save; 688 if (!okq) { say(c, "No se puede ajustar"); return; } 689 } else { 690 /* se linealiza al vuelo (sin arreglos: la pila y la RAM de la Pico son chicas) */ 691 static Num sX, sY, sXX, sXY, sYY, X, Y; 692 sX = sY = sXX = sXY = sYY = num_int(0); 693 for (int i = 0; i < n; i++) { 694 X = s->x[i]; Y = s->y[i]; 695 if (t == ST_LOG || t == ST_POW) X = num_ln(s->x[i]); 696 if (t == ST_EXP || t == ST_ABEXP || t == ST_POW) Y = num_ln(s->y[i]); 697 if (t == ST_INV) X = num_div(num_int(1), s->x[i]); 698 if (num_bad(&X) || num_bad(&Y)) { say(c, "Datos fuera del dominio del modelo"); return; } 699 sX = num_add(sX, X); sY = num_add(sY, Y); 700 sXX = num_add(sXX, num_mul(X, X)); sXY = num_add(sXY, num_mul(X, Y)); 701 sYY = num_add(sYY, num_mul(Y, Y)); 702 } 703 Num Sxx = num_sub(sXX, num_div(num_mul(sX, sX), N)); 704 Num Sxy = num_sub(sXY, num_div(num_mul(sX, sY), N)); 705 Num Syy = num_sub(sYY, num_div(num_mul(sY, sY), N)); 706 if (!Sxx.v.m) { say(c, "Todas las x iguales"); return; } 707 b = num_div(Sxy, Sxx); 708 a = num_sub(num_div(sY, N), num_mul(b, num_div(sX, N))); 709 r = Syy.v.m ? num_div(Sxy, num_sqrt(num_mul(Sxx, Syy))) : num_err(); 710 if (t == ST_EXP || t == ST_ABEXP || t == ST_POW) a = num_dec(d_exp(a.v)); 711 if (t == ST_ABEXP) b = num_dec(d_exp(b.v)); 712 } 713 res(s, "a", a); 714 res(s, "b", b); 715 if (t == ST_QUAD) res(s, "c", cc); 716 else res(s, "r", r); 717 s->ra = a; s->rb = b; s->rc = cc; 718 s->rtype = t; 719 s->rvalid = true; 720 } 721 722 static bool stat_key(Calc *c, int key) 723 { 724 StatState *s = &c->st; 725 int cols = s->type == ST_1VAR ? 1 : 2; 726 if (s->phase == 1) { 727 if (key == K_UP && s->rtop > 0) s->rtop--; 728 else if (key == K_DOWN && s->rtop + 1 < s->nres) s->rtop++; 729 else if (key == K_BACK || key == K_OK) s->phase = 0; 730 return true; 731 } 732 if (s->editing) { 733 if (key == K_OK) { 734 Num v; 735 if (!eval_in(c, &v)) return true; 736 if (s->row >= s->n) { 737 if (s->n >= STAT_MAX) { say(c, "Máximo 80 datos"); return true; } 738 s->x[s->n] = num_int(0); s->y[s->n] = num_int(1); 739 if (s->type != ST_1VAR) s->y[s->n] = num_int(0); 740 s->n++; 741 } 742 if (s->col == 0) s->x[s->row] = v; else s->y[s->row] = v; 743 s->editing = false; 744 ex_clear(&c->in); 745 if (cols == 2 && s->col == 0) s->col = 1; 746 else { s->col = 0; s->row++; } 747 c->dirty_save = true; 748 return true; 749 } 750 if (key == K_BACK) { s->editing = false; ex_clear(&c->in); return true; } 751 return edit_key(c, &c->in, key); 752 } 753 switch (key) { 754 case K_UP: if (s->row > 0) s->row--; return true; 755 case K_DOWN: if (s->row < s->n) s->row++; return true; 756 case K_LEFT: if (s->col > 0) s->col--; return true; 757 case K_RIGHT: if (s->col + 1 < cols) s->col++; return true; 758 case K_DEL: case K_BKSP: 759 if (s->row < s->n) { 760 memmove(&s->x[s->row], &s->x[s->row + 1], sizeof s->x[0] * (size_t)(s->n - s->row - 1)); 761 memmove(&s->y[s->row], &s->y[s->row + 1], sizeof s->y[0] * (size_t)(s->n - s->row - 1)); 762 s->n--; 763 c->dirty_save = true; 764 } 765 return true; 766 case K_OK: case '=': stat_compute(c); if (s->nres) { s->phase = 1; s->rtop = 0; } return true; 767 case K_BACK: stat_form(c); return true; 768 } 769 if (key >= 32 && key < 127) { 770 s->editing = true; 771 ex_clear(&c->in); 772 return edit_key(c, &c->in, key); 773 } 774 return false; 775 } 776 777 /* ============================================================ matrices */ 778 779 static MVal *mat_target(Calc *c) 780 { 781 int t = c->mt.target; 782 return t < 4 ? &c->mt.m[t] : &c->mt.v[t - 4]; 783 } 784 785 void mat_pick_form(Calc *c) 786 { 787 static char lab[8][32]; 788 Form *f = &c->form; 789 form_begin(f, "Editar matrices y vectores"); 790 for (int i = 0; i < 8; i++) { 791 const MVal *m = i < 4 ? &c->mt.m[i] : &c->mt.v[i - 4]; 792 int p = str_put(lab[i], 32, 0, i < 4 ? "Mat" : "Vct"); 793 char l[2] = { (char)('A' + i % 4), 0 }; 794 p = str_put(lab[i], 32, p, l); 795 if (m->r) { 796 p = str_put(lab[i], 32, p, " "); 797 if (i < 4) { p = str_int(lab[i], 32, p, m->r); p = str_put(lab[i], 32, p, "×"); } 798 p = str_int(lab[i], 32, p, m->c); 799 if (i >= 4) str_put(lab[i], 32, p, " elementos"); 800 } else str_put(lab[i], 32, p, " (sin definir)"); 801 form_button(f, 300 + i, lab[i]); 802 } 803 f->cur = c->mt.target; 804 c->form_kind = FORM_MAT_PICK; 805 c->scr = SCR_FORM; 806 } 807 808 static void mat_dim_form(Calc *c) 809 { 810 static const char *const N4[] = { "1", "2", "3", "4" }; 811 static const char *const N23[] = { "2", "3" }; 812 Form *f = &c->form; 813 MVal *m = mat_target(c); 814 bool vec = c->mt.target >= 4; 815 form_begin(f, vec ? "Dimensión del vector" : "Dimensión de la matriz"); 816 if (vec) form_choice(f, 2, "Elementos", N23, 2, m->r && m->c == 3 ? 1 : 0); 817 else { 818 form_choice(f, 1, "Filas", N4, 4, m->r ? m->r - 1 : 1); 819 form_choice(f, 2, "Columnas", N4, 4, m->r ? m->c - 1 : 1); 820 } 821 form_button(f, 10, "Editar los valores"); 822 f->cur = f->n - 1; 823 c->form_kind = FORM_MAT_DIM; 824 c->scr = SCR_FORM; 825 } 826 827 static bool matrix_key(Calc *c, int key) 828 { 829 MatState *t = &c->mt; 830 MVal *m = mat_target(c); 831 int cells = m->r * m->c; 832 if (t->editing) { 833 if (key == K_OK) { 834 Num v; 835 if (!eval_in(c, &v)) return true; 836 *mv_at(m, t->sel / m->c, t->sel % m->c) = v; 837 t->editing = false; 838 ex_clear(&c->in); 839 t->sel = (t->sel + 1) % cells; 840 c->dirty_save = true; 841 return true; 842 } 843 if (key == K_BACK) { t->editing = false; ex_clear(&c->in); return true; } 844 return edit_key(c, &c->in, key); 845 } 846 switch (key) { 847 case K_LEFT: t->sel = (t->sel + cells - 1) % cells; return true; 848 case K_RIGHT: t->sel = (t->sel + 1) % cells; return true; 849 case K_UP: t->sel = (t->sel + cells - m->c) % cells; return true; 850 case K_DOWN: t->sel = (t->sel + m->c) % cells; return true; 851 case K_BKSP: case K_DEL: *mv_at(m, t->sel / m->c, t->sel % m->c) = num_int(0); return true; 852 case K_BACK: case K_OK: t->phase = 0; ex_clear(&c->in); return true; 853 case '=': mat_pick_form(c); return true; 854 } 855 if (key >= 32 && key < 127) { 856 t->editing = true; 857 ex_clear(&c->in); 858 return edit_key(c, &c->in, key); 859 } 860 return false; 861 } 862 863 /* ============================================================ despacho */ 864 865 void mode_enter(Calc *c, int mode) 866 { 867 c->mode = mode; 868 ex_clear(&c->in); 869 c->shown = false; 870 c->err = 0; 871 switch (mode) { 872 case MODE_EQN: eq_type_form(c); break; 873 case MODE_TABLE: c->tb.phase = 0; c->in = c->tb.f; c->in.cur = c->in.n; break; 874 case MODE_GRAPH: c->gr.phase = 0; c->in = c->gr.f; c->in.cur = c->in.n; break; 875 case MODE_STAT: stat_form(c); break; 876 case MODE_PROG: say(c, "Tab cambia la base"); break; 877 case MODE_CMPLX: say(c, "i imaginaria < ∠ y polar"); break; 878 case MODE_MATRIX: c->mt.phase = 0; say(c, "= define A-D MatA ` A VctA"); break; 879 } 880 c->dirty_save = true; 881 } 882 883 bool mode_form_done(Calc *c, int id) 884 { 885 Form *f = &c->form; 886 switch (c->form_kind) { 887 case FORM_EQN_TYPE: 888 if (id >= 100 && id < 105) { 889 int t = id - 100; 890 if (t != c->eq.type) { 891 for (int i = 0; i < 12; i++) c->eq.coef[i] = num_int(0); 892 c->eq.sel = 0; 893 } 894 c->eq.type = t; 895 c->eq.phase = 0; 896 c->eq.editing = false; 897 ex_clear(&c->in); 898 if (t == EQ_SOLVE) { c->in = c->eq.solve_eq; c->in.cur = c->in.n; } 899 c->scr = SCR_MAIN; 900 } 901 return true; 902 case FORM_MAT_PICK: 903 if (id >= 300 && id < 308) { c->mt.target = id - 300; mat_dim_form(c); } 904 return true; 905 case FORM_MAT_DIM: { 906 if (id != 10) return true; 907 MVal *m = mat_target(c); 908 bool vec = c->mt.target >= 4; 909 int r = vec ? 1 : form_get(f, 1)->value + 1, cc = form_get(f, 2)->value + (vec ? 2 : 1); 910 if (m->r != r || m->c != cc) { 911 /* conservar lo que entre en la dimension nueva */ 912 static MVal old; 913 old = *m; 914 mv_zero(m, vec ? MV_VEC : MV_MAT, r, cc); 915 for (int i = 0; i < r && i < old.r; i++) 916 for (int j = 0; j < cc && j < old.c; j++) *mv_at(m, i, j) = *mv_get(&old, i, j); 917 } 918 c->mt.phase = 1; 919 c->mt.sel = 0; 920 c->mt.editing = false; 921 c->scr = SCR_MAIN; 922 c->dirty_save = true; 923 return true; 924 } 925 case FORM_STAT_TYPE: 926 if (id >= 200 && id < 200 + ST_NTYPES) { 927 c->st.type = id - 200; 928 c->st.phase = 0; c->st.row = c->st.col = 0; 929 c->scr = SCR_MAIN; 930 } 931 return true; 932 case FORM_TABLE: { 933 if (id != 5) return true; 934 Num a, b, s; 935 if (!parse_num(c, form_get(f, 1)->text, &a) || !parse_num(c, form_get(f, 2)->text, &b) || 936 !parse_num(c, form_get(f, 3)->text, &s)) { form_msg(f, "Valor inválido", true); return true; } 937 c->tb.start = a; c->tb.end = b; c->tb.step = s; 938 c->tb.use_g = form_get(f, 4)->value; 939 table_compute(c); 940 c->tb.phase = 2; 941 c->scr = SCR_MAIN; 942 return true; 943 } 944 case FORM_GRAPH: { 945 if (id != 7) return true; 946 Num v[4]; 947 for (int i = 0; i < 4; i++) 948 if (!parse_num(c, form_get(f, i + 1)->text, &v[i])) { form_msg(f, "Valor inválido", true); return true; } 949 if (dec_cmp(v[0].v, v[1].v) >= 0) { form_msg(f, "x mín tiene que ser menor", true); return true; } 950 c->gr.xmin = v[0].v; c->gr.xmax = v[1].v; 951 c->gr.autoy = form_get(f, 5)->value; 952 if (!c->gr.autoy) { 953 if (dec_cmp(v[2].v, v[3].v) >= 0) { form_msg(f, "y mín tiene que ser menor", true); return true; } 954 c->gr.ymin = v[2].v; c->gr.ymax = v[3].v; 955 } 956 c->gr.use_g = form_get(f, 6)->value; 957 c->gr.ready = false; 958 c->gr.phase = 2; 959 c->gr.tx = -1; 960 c->scr = SCR_MAIN; 961 return true; 962 } 963 } 964 c->scr = SCR_MAIN; 965 return true; 966 } 967 968 bool mode_key(Calc *c, int key) 969 { 970 switch (c->mode) { 971 case MODE_EQN: return eqn_key(c, key); 972 case MODE_TABLE: return table_key(c, key); 973 case MODE_GRAPH: return graph_key(c, key); 974 case MODE_STAT: return stat_key(c, key); 975 case MODE_MATRIX: return matrix_key(c, key); 976 } 977 return edit_key(c, &c->in, key); 978 } 979 980 Expr *mode_target(Calc *c) 981 { 982 if (c->mode == MODE_EQN && c->eq.phase == 0 && c->eq.type != EQ_SOLVE) c->eq.editing = true; 983 if (c->mode == MODE_STAT && c->st.phase == 0) c->st.editing = true; 984 if (c->mode == MODE_MATRIX && c->mt.phase == 1) c->mt.editing = true; 985 return &c->in; 986 } 987 988 /* para los frontends */ 989 void graph_trace(Calc *c) { trace_value(c); }