diff options
| -rw-r--r-- | yacas.c | 350 |
1 files changed, 243 insertions, 107 deletions
@@ -8,6 +8,9 @@ #include <stdint.h> #include <stdlib.h> #include <string.h> +#include <time.h> +#include <unistd.h> + #ifdef USE_READLINE #include <readline/readline.h> #include <readline/history.h> @@ -28,6 +31,8 @@ struct expression { struct expression *operand, *operand2; } specfunc; /* T_SPECFUNC */ }; + + int refcount; }; enum { OP_POW = 0, OP_MUL, OP_DIV, OP_ADD, OP_SUB, OP_DIFF, OP_ASSIGN, OP_EQUALS, OP_LAST, /* not a real operator: */ OP_LPAREN }; @@ -36,7 +41,7 @@ enum { FN_LOG = 0, FN_EXP, FN_SIN, FN_COS, FN_TAN, FN_CSC, FN_SEC, FN_COT, FN_AS enum { RIGHT, LEFT }; const char *op_names[] = { - "^", "*", "/", "+", "-", "diff", "=", "==" + "^", "*", "/", "+", "-", "diff", "=", "==", "none", "(" }; const char *type_names[] = { @@ -58,7 +63,9 @@ const char *fn_names[] = { bool is_constant(struct expression *expr); bool is_constvar(const char *name); bool is_var(const char *name); +bool check_boolean(const char *varname); double eval_numexpr(struct expression *expr); +double eval_numexpr_and_free(struct expression *expr); struct expression *diff(struct expression *expr, const char *var); struct expression *dup_expr(struct expression *expr); struct expression *eval_expr(struct expression *expr); @@ -68,6 +75,7 @@ void setvar(const char *name, struct expression *value, bool constant); bool simp_progress = false; bool memdiag = false; +bool interactive = true; void __attribute__((noreturn)) fatal(const char *msg) { @@ -94,7 +102,7 @@ struct expression *new_var(const char *varname) new->varname = strdup(varname); memusage += sizeof(*new); check_mem(); - if(memdiag) + if(check_boolean("memdiag")) printf("new var %p = %s\n", new, varname); return new; } @@ -106,7 +114,7 @@ struct expression *new_const(double x) new->constant = x; memusage += sizeof(*new); check_mem(); - if(memdiag) + if(check_boolean("memdiag")) printf("new const %p = %g\n", new, x); return new; } @@ -120,8 +128,8 @@ struct expression *new_op(int op, struct expression *l, struct expression *r) new->subexpr.right = r; memusage += sizeof(*new); check_mem(); - if(memdiag) - printf("new op %p\n", new); + if(check_boolean("memdiag")) + printf("new op %p %s\n", new, op_names[op]); return new; } @@ -132,8 +140,8 @@ struct expression *new_specfunc(int fn, struct expression *expr) new->specfunc.fn = fn; new->specfunc.operand = expr; memusage += sizeof(*new); - if(memdiag) - printf("new specfunc %p\n", new); + if(check_boolean("memdiag")) + printf("new specfunc %p %s\n", new, fn_names[fn]); return new; } @@ -192,7 +200,7 @@ void free_expr(struct expression *expr) break; } memusage -= sizeof(struct expression); - if(memdiag) + if(check_boolean("memdiag")) printf("freeing %p\n", expr); free(expr); } @@ -255,6 +263,7 @@ bool expr_references(struct expression *expr, const char *varname) /* simple stuff, doesn't always work */ bool expr_equals(struct expression *l, struct expression *r) { + /* first the low-hanging fruit */ if(l == r) return true; @@ -265,24 +274,41 @@ bool expr_equals(struct expression *l, struct expression *r) l = simp_and_free_max(l); r = simp_and_free_max(r); + /* must come after simplification */ if(is_constant(l) && is_constant(r)) - return eval_numexpr(l) == eval_numexpr(r); + return eval_numexpr_and_free(l) == eval_numexpr_and_free(r); /* not foolproof */ if(l->type != r->type) + { + free_expr(l); + free_expr(r); return false; + } + + bool ret = false; switch(l->type) { case T_SUBEXPR: - /* todo */ - return false; + ret = false; + free_expr(l); + free_expr(r); + return ret; case T_SPECFUNC: - return l->specfunc.fn == r->specfunc.fn && expr_equals(l->specfunc.operand, r->specfunc.operand); + ret = l->specfunc.fn == r->specfunc.fn && expr_equals(l->specfunc.operand, r->specfunc.operand); + free_expr(l); + free_expr(r); + return ret; case T_VAR: /* symbolic comparison */ if(!is_var(l->varname) && !is_var(r->varname)) - return !strcmp(l->varname, r->varname); + { + ret = !strcmp(l->varname, r->varname); + free_expr(l); + free_expr(r); + return ret; + } break; default: break; @@ -291,14 +317,20 @@ bool expr_equals(struct expression *l, struct expression *r) /* try evaluating them and recurse */ l = eval_and_free(l); r = eval_and_free(r); - return expr_equals(l, r); + + ret = expr_equals(l, r); + + free_expr(l); + free_expr(r); + + return ret; } -#define VARMAP_SIZE 10 +#define VARMAP_SIZE 20 /* chained hash map */ struct variable { - const char *name; + char *name; /* dynamic */ struct expression *value; bool constant; struct variable *next; @@ -380,6 +412,7 @@ bool is_var(const char *name) uint32_t idx = var_hash(name) % VARMAP_SIZE; struct variable *i = varmap[idx]; + while(i) { if(!strcmp(i->name, name)) @@ -405,6 +438,47 @@ struct expression *getvar(const char *name) return NULL; } +bool delvar(const char *name) +{ + uint32_t idx = var_hash(name) % VARMAP_SIZE; + + struct variable *i = varmap[idx], *last = NULL; + + while(i) + { + if(!strcmp(i->name, name)) + { + if(last) + last->next = i->next; + else + varmap[idx] = NULL; + free_expr(i->value); + free(i->name); + free(i); + return true; + } + last = i; + i = i->next; + } + + return false; +} + +/* check whether a variable is nonzero without any intermediate functions */ +bool check_boolean(const char *varname) +{ + struct expression *expr = getvar(varname); + if(!expr) + return false; + switch(expr->type) + { + case T_CONST: + return expr->constant == 1; + default: + return false; + } +} + bool is_constant(struct expression *expr) { switch(expr->type) @@ -516,9 +590,45 @@ struct expression *eval_op(int op, struct expression *lexpr, struct expression * { if(op == OP_DIFF) { - if(rexpr->type != T_VAR) + /* allow for syntax in the form f diff x = 2 for evaluating + * the derivative at x = 2 */ + switch(rexpr->type) + { + case T_VAR: + return diff_and_free(eval_expr(lexpr), rexpr->varname); + case T_SUBEXPR: + if(rexpr->subexpr.op == OP_ASSIGN && rexpr->subexpr.left->type == T_VAR && is_constant(rexpr->subexpr.right)) + { + bool have_old = is_var(rexpr->subexpr.left->varname); + + /* set variable (duplicate because setvar frees it */ + struct expression *old; + + if(have_old) + { + old = dup_expr(getvar(rexpr->subexpr.left->varname)); + /* delete so it's not a constant */ + delvar(rexpr->subexpr.left->varname); + } + + /* get derivative */ + struct expression *deriv = diff_and_free(eval_expr(lexpr), rexpr->subexpr.left->varname); + + setvar(rexpr->subexpr.left->varname, dup_expr(rexpr->subexpr.right), false); + + double val = eval_numexpr_and_free(deriv); + + /* don't care about const */ + if(have_old) + setvar(rexpr->subexpr.left->varname, old, false); + /* otherwise keep it around */ + + return new_const(val); + } + /* fall through */ + default: fatal("diff operator must take variable as second operand"); - return diff(eval_expr(lexpr), rexpr->varname); + } } else if(op == OP_ASSIGN) { @@ -615,6 +725,13 @@ double eval_numexpr(struct expression *expr) return x; } +double eval_numexpr_and_free(struct expression *expr) +{ + double ret = eval_numexpr(expr); + free_expr(expr); + return ret; +} + /* returns an entirely new expression object */ struct expression *simp(struct expression *expr) { @@ -786,10 +903,10 @@ struct expression *diff(struct expression *expr, const char *var) return new_const(1); else { - if(is_var(expr->varname)) - return diff(eval_expr(expr), var); //if(is_constant(expr)) //return new_const(0); + if(is_var(expr->varname)) + return diff(eval_expr(expr), var); return new_op(OP_DIFF, dup_expr(expr), new_var(var)); } case T_SUBEXPR: @@ -797,61 +914,61 @@ struct expression *diff(struct expression *expr, const char *var) { case OP_ADD: case OP_SUB: - return new_op(expr->subexpr.op, simp_and_free(diff(expr->subexpr.left, var)), simp_and_free(diff(expr->subexpr.right, var))); + return simp_and_free(new_op(expr->subexpr.op, diff(expr->subexpr.left, var), diff(expr->subexpr.right, var))); case OP_MUL: /* f'g + fg' */ - return new_op(OP_ADD, - simp_and_free(new_op(OP_MUL, simp_and_free(diff(expr->subexpr.left, var)), dup_expr(expr->subexpr.right))), - simp_and_free(new_op(OP_MUL, dup_expr(expr->subexpr.left), simp_and_free(diff(expr->subexpr.right, var))))); + return simp_and_free(new_op(OP_ADD, + new_op(OP_MUL, diff(expr->subexpr.left, var), dup_expr(expr->subexpr.right)), + new_op(OP_MUL, dup_expr(expr->subexpr.left), diff(expr->subexpr.right, var)))); case OP_DIV: - return new_op(OP_DIV, - new_op(OP_SUB, - new_op(OP_MUL, simp_and_free(diff(expr->subexpr.left, var)), dup_expr(expr->subexpr.right)), - new_op(OP_MUL, dup_expr(expr->subexpr.left), simp_and_free(diff(expr->subexpr.right, var)))), - new_op(OP_POW, dup_expr(expr->subexpr.right), new_const(2))); + return simp_and_free(new_op(OP_DIV, + new_op(OP_SUB, + new_op(OP_MUL, diff(expr->subexpr.left, var), dup_expr(expr->subexpr.right)), + new_op(OP_MUL, dup_expr(expr->subexpr.left), diff(expr->subexpr.right, var))), + new_op(OP_POW, dup_expr(expr->subexpr.right), new_const(2)))); case OP_POW: /* just n * f(x) ^ (n - 1) * f'(x) */ -#if 0 if(is_constant(expr->subexpr.right)) { - return new_op(OP_MUL, - new_const(eval_numexpr(expr->subexpr.right)), - new_op(OP_MUL, - simp_and_free(new_op(OP_POW, + return simp_and_free(new_op(OP_MUL, + new_const(eval_numexpr(expr->subexpr.right)), + new_op(OP_MUL, + new_op(OP_POW, dup_expr(expr->subexpr.left), new_op(OP_SUB, dup_expr(expr->subexpr.right), - new_const(1)))), - simp_and_free(diff(expr->subexpr.left, var)))); + new_const(1))), + diff(expr->subexpr.left, var)))); } else -#endif { /* f^(g-1) * (gf' + fg'ln(f)) */ - return new_op(OP_MUL, - new_op(OP_POW, - dup_expr(expr->subexpr.left), - new_op(OP_SUB, - dup_expr(expr->subexpr.right), - new_const(1))), - new_op(OP_ADD, - new_op(OP_MUL, - dup_expr(expr->subexpr.right), - simp_and_free(diff(expr->subexpr.left, var))), - new_op(OP_MUL, - dup_expr(expr->subexpr.left), - new_op(OP_MUL, - simp_and_free(diff(expr->subexpr.right, var)), - new_specfunc(FN_LOG, dup_expr(expr->subexpr.left)))))); + return simp_and_free(new_op(OP_MUL, + new_op(OP_POW, + dup_expr(expr->subexpr.left), + new_op(OP_SUB, + dup_expr(expr->subexpr.right), + new_const(1))), + new_op(OP_ADD, + new_op(OP_MUL, + dup_expr(expr->subexpr.right), + diff(expr->subexpr.left, var)), + new_op(OP_MUL, + dup_expr(expr->subexpr.left), + new_op(OP_MUL, + diff(expr->subexpr.right, var), + new_specfunc(FN_LOG, dup_expr(expr->subexpr.left))))))); } case OP_DIFF: - return diff_and_free(diff(expr->subexpr.left, expr->subexpr.right->varname), var); + //if(expr->subexpr.left->type == T_VAR && !is_var(expr->subexpr.left->varname)) + return simp_and_free(new_op(OP_DIFF, dup_expr(expr), new_var(var))); + /*return diff_and_free(diff(expr->subexpr.left, expr->subexpr.right->varname), var);*/ case OP_ASSIGN: /* todo: how to make this logical? */ - break; + fatal("cannot differentiate assignment"); case OP_EQUALS: + return simp_and_free(new_op(OP_EQUALS, diff(expr->subexpr.left, var), diff(expr->subexpr.right, var))); /* todo: how to make this logical? */ - break; default: fatal("fall through"); } @@ -859,73 +976,73 @@ struct expression *diff(struct expression *expr, const char *var) switch(expr->specfunc.fn) { case FN_LOG: - return new_op(OP_DIV, - simp_and_free(diff(expr->specfunc.operand, var)), - dup_expr(expr->specfunc.operand)); + return simp_and_free(new_op(OP_DIV, + diff(expr->specfunc.operand, var), + dup_expr(expr->specfunc.operand))); case FN_EXP: - return new_op(OP_MUL, - simp_and_free(diff(expr->specfunc.operand, var)), - new_specfunc(FN_EXP, dup_expr(expr->specfunc.operand))); + return simp_and_free(new_op(OP_MUL, + diff(expr->specfunc.operand, var), + new_specfunc(FN_EXP, dup_expr(expr->specfunc.operand)))); case FN_SIN: - return new_op(OP_MUL, + return simp_and_free(new_op(OP_MUL, new_specfunc(FN_COS, dup_expr(expr->specfunc.operand)), - simp_and_free(diff(expr->specfunc.operand, var))); + diff(expr->specfunc.operand, var))); case FN_COS: - return new_op(OP_MUL, + return simp_and_free(new_op(OP_MUL, new_const(-1), new_op(OP_MUL, new_specfunc(FN_SIN, dup_expr(expr->specfunc.operand)), - simp_and_free(diff(expr->specfunc.operand, var)))); + diff(expr->specfunc.operand, var)))); case FN_TAN: - return new_op(OP_DIV, - simp_and_free(diff(expr->specfunc.operand, var)), + return simp_and_free(new_op(OP_DIV, + diff(expr->specfunc.operand, var), new_op(OP_POW, new_specfunc(FN_COS, dup_expr(expr->specfunc.operand)), - new_const(2))); + new_const(2)))); case FN_CSC: - return new_op(OP_MUL, + return simp_and_free(new_op(OP_MUL, new_const(-1), new_op(OP_MUL, - simp_and_free(diff(expr->specfunc.operand, var)), + diff(expr->specfunc.operand, var), new_op(OP_POW, new_op(OP_SUB, new_const(1), new_op(OP_POW, dup_expr(expr->specfunc.operand), new_const(2))), - new_const(-.5)))); + new_const(-.5))))); case FN_SEC: - return new_op(OP_MUL, - simp_and_free(diff(expr->specfunc.operand, var)), + return simp_and_free(new_op(OP_MUL, + diff(expr->specfunc.operand, var), new_op(OP_POW, new_op(OP_SUB, new_const(1), new_op(OP_POW, dup_expr(expr->specfunc.operand), new_const(2))), - new_const(-.5))); + new_const(-.5)))); case FN_COT: - return new_op(OP_MUL, + return simp_and_free(new_op(OP_MUL, new_const(-1), new_op(OP_MUL, - simp_and_free(diff(expr->specfunc.operand, var)), + diff(expr->specfunc.operand, var), new_op(OP_POW, new_specfunc(FN_CSC, dup_expr(expr->specfunc.operand)), - new_const(2)))); + new_const(2))))); case FN_ASIN: - return new_op(OP_MUL, - simp_and_free(diff(expr->specfunc.operand, var)), + return simp_and_free(new_op(OP_MUL, + diff(expr->specfunc.operand, var), new_op(OP_POW, new_op(OP_SUB, new_const(1), new_op(OP_POW, dup_expr(expr->specfunc.operand), new_const(2))), - new_const(-.5))); + new_const(-.5)))); case FN_ACOS: - return new_op(OP_MUL, - simp_and_free(diff(expr->specfunc.operand, var)), + return simp_and_free(new_op(OP_MUL, + diff(expr->specfunc.operand, var), new_op(OP_MUL, new_const(-1), new_op(OP_POW, @@ -934,20 +1051,20 @@ struct expression *diff(struct expression *expr, const char *var) new_op(OP_POW, dup_expr(expr->specfunc.operand), new_const(2))), - new_const(-.5)))); + new_const(-.5))))); case FN_ATAN: - return new_op(OP_DIV, - simp_and_free(diff(expr->specfunc.operand, var)), + return simp_and_free(new_op(OP_DIV, + diff(expr->specfunc.operand, var), new_op(OP_ADD, new_const(1), new_op(OP_POW, dup_expr(expr->specfunc.operand), - new_const(2)))); + new_const(2))))); case FN_ACSC: - return new_op(OP_MUL, + return simp_and_free(new_op(OP_MUL, new_const(-1), new_op(OP_MUL, - simp_and_free(diff(expr->specfunc.operand, var)), + diff(expr->specfunc.operand, var), new_op(OP_DIV, new_op(OP_POW, new_op(OP_SUB, @@ -956,10 +1073,10 @@ struct expression *diff(struct expression *expr, const char *var) new_const(2)), new_const(1)), new_const(-.5)), - new_specfunc(FN_ABS, dup_expr(expr->specfunc.operand))))); + new_specfunc(FN_ABS, dup_expr(expr->specfunc.operand)))))); case FN_ASEC: - return new_op(OP_MUL, - simp_and_free(diff(expr->specfunc.operand, var)), + return simp_and_free(new_op(OP_MUL, + diff(expr->specfunc.operand, var), new_op(OP_DIV, new_op(OP_POW, new_op(OP_SUB, @@ -968,32 +1085,32 @@ struct expression *diff(struct expression *expr, const char *var) new_const(2)), new_const(1)), new_const(-.5)), - new_specfunc(FN_ABS, dup_expr(expr->specfunc.operand)))); + new_specfunc(FN_ABS, dup_expr(expr->specfunc.operand))))); case FN_ACOT: - return new_op(OP_MUL, + return simp_and_free(new_op(OP_MUL, new_const(-1), new_op(OP_DIV, - simp_and_free(diff(expr->specfunc.operand, var)), + diff(expr->specfunc.operand, var), new_op(OP_ADD, new_const(1), new_op(OP_POW, dup_expr(expr->specfunc.operand), - new_const(2))))); + new_const(2)))))); case FN_ABS: - return new_op(OP_MUL, - simp_and_free(diff(expr->specfunc.operand, var)), + return simp_and_free(new_op(OP_MUL, + diff(expr->specfunc.operand, var), new_op(OP_DIV, new_specfunc(FN_ABS, dup_expr(expr->specfunc.operand)), - dup_expr(expr->specfunc.operand))); + dup_expr(expr->specfunc.operand)))); case FN_NEGATE: - return new_specfunc(FN_NEGATE, simp_and_free(diff(expr->specfunc.operand, var))); + return simp_and_free(new_specfunc(FN_NEGATE, diff(expr->specfunc.operand, var))); case FN_SQRT: - return new_op(OP_DIV, - simp_and_free(diff(expr->specfunc.operand, var)), + return simp_and_free(new_op(OP_DIV, + diff(expr->specfunc.operand, var), new_op(OP_MUL, new_const(2), new_specfunc(FN_SQRT, - dup_expr(expr->specfunc.operand)))); + dup_expr(expr->specfunc.operand))))); default: fatal("fall through"); } @@ -1302,25 +1419,33 @@ struct expression *parse_expr(const char *infix) int main(int argc, char *argv[]) { + interactive = isatty(STDIN_FILENO); + setvar("pi", new_const(M_PI), true); setvar("e", new_const(M_E), true); + setvar("memdiag", new_const(0), false); + while(1) { peakmem = 0; char *line = NULL; #ifndef USE_READLINE - printf("yaCAS> "); - fflush(stdout) + if(interactive) + printf("yaCAS> "); + fflush(stdout); size_t n = 0; - size_t len = getline(&line, &n, stdin); + ssize_t len = getline(&line, &n, stdin); /* remove newline */ if(line[len - 1] == '\n') line[len-- - 1] = '\0'; + + if(len < 0) + return 0; #else - line = readline("yaCAS> "); + line = readline(interactive ? "yaCAS> " : ""); if(line && *line) add_history(line); #endif @@ -1335,6 +1460,8 @@ int main(int argc, char *argv[]) continue; } + clock_t start = clock(); + struct expression *top = parse_expr(line); free(line); @@ -1343,6 +1470,7 @@ int main(int argc, char *argv[]) top = eval_and_free(top); top = simp_and_free_max(top); print_expr(top); + printf("(%s)\n", type_names[top->type]); setvar(".", dup_expr(top), true); @@ -1350,9 +1478,17 @@ int main(int argc, char *argv[]) free_expr(top); } else + { printf("malformed expression\n"); + } + + clock_t diff = clock() - start; + if(diff > CLOCKS_PER_SEC / 10) + printf("time: %d.%.03d sec\n", diff / CLOCKS_PER_SEC, diff % CLOCKS_PER_SEC / (CLOCKS_PER_SEC / 1000)); - if(memdiag) + if(eval_numexpr_and_free(new_var("memdiag")) == 2) + { printf("memory: %d peak: %d\n", memusage, peakmem); + } } } |