/*
* Mathomatic symbolic solve routines.
*
* Copyright (C) 1987-2007 George Gesslein II.
*/
#include "includes.h"
#define MAX_RAISE_POWER 20 /* Maximum number of times to increase power in solve function. */
static int flip();
static int g_of_f();
static int quad_solve();
static int increase();
static int repeat_count;
/*
* Solve using equation spaces.
*
* Return true if successful.
*/
int
solve(want, have)
int want; /* equation number containing what to solve for */
int have; /* equation number to solve */
{
int rv;
if (n_lhs[want] > 0) {
rv = solve_sub(lhs[want], n_lhs[want], lhs[have], &n_lhs[have], rhs[have], &n_rhs[have]);
} else {
rv = solve_sub(rhs[want], n_rhs[want], rhs[have], &n_rhs[have], lhs[have], &n_lhs[have]);
}
n_lhs[want] = 0;
n_rhs[want] = 0;
#if !SILENT
if (rv <= 0) {
printf(_("Solve failed.\n"));
}
#endif
return(rv > 0);
}
/*
* Main Mathomatic symbolic solve routine.
*
* This works by moving everything containing the variable to solve for
* to the LHS (via transposition), then moving everything not containing the variable to the
* RHS. Many tricks are used, and this routine works very well.
*
* tlhs and trhs are used by this routine.
*
* Returns a positive integer if successful, with the result placed in the passed LHS and RHS.
* Returns 1 for normal success.
* Returns 2 if successful and a solution was zero and removed.
* Returns 0 on failure. Might succeed at a numeric solve.
* Returns -1 if the equation is an identity.
* Returns -2 if unsolvable in all realms.
*/
int
solve_sub(wantp, wantn, leftp, leftnp, rightp, rightnp)
token_type *wantp; /* expression to solve for */
int wantn; /* size of expression to solve for (currently should always be 1) */
token_type *leftp; /* LHS of equation */
int *leftnp; /* size of LHS */
token_type *rightp; /* RHS of equation */
int *rightnp; /* size of RHS */
{
int i, j;
int found, found_count;
int worked;
int diff_sign;
int op, op_kind;
token_type *p1, *b1, *ep;
long v;
int need_flip;
int uf_flag = false; /* unfactor flag */
int qtries = 0;
int zflag; /* true if RHS is currently zero */
int zsolve; /* true if we are solving for zero */
int inc_count = 0;
int zero_solved = false;
double numerator, denominator;
int success = 1;
repeat_count = 0;
/* copy the equation to temporary storage where it will be manipulated */
n_tlhs = *leftnp;
blt(tlhs, leftp, n_tlhs * sizeof(*leftp));
n_trhs = *rightnp;
blt(trhs, rightp, n_trhs * sizeof(*rightp));
if (wantn != 1) {
error(_("This program will only solve for a single variable or for zero."));
return false;
}
if (n_tlhs <= 0 || n_trhs <= 0) {
error(_("Please enter an equation or type \"help\"."));
return false;
}
if (wantp->kind == VARIABLE) {
v = wantp->token.variable;
if (!found_var(trhs, n_trhs, v) && !found_var(tlhs, n_tlhs, v)) {
error(_("Variable not found."));
return false;
}
zsolve = false;
} else {
v = 0;
if (wantp->kind != CONSTANT || wantp->token.constant != 0.0) {
error(_("This program will only solve for a single variable or for zero."));
return false;
}
zsolve = true;
}
uf_power(tlhs, &n_tlhs);
uf_power(trhs, &n_trhs);
simp_again:
/* Make sure equation is a bit simplified. */
list_tdebug(2);
simps_side(tlhs, &n_tlhs, zsolve);
if (uf_flag) {
simp_loop(trhs, &n_trhs);
uf_simp(trhs, &n_trhs);
factorv(trhs, &n_trhs, v);
} else {
simps_side(trhs, &n_trhs, zsolve);
}
list_tdebug(1);
no_simp:
/* First selectively move subexpressions from the RHS to the LHS. */
op = 0;
ep = &trhs[n_trhs];
if (zsolve) {
for (b1 = p1 = trhs; p1 < ep; p1++) {
if (p1->level == 1 && p1->kind == OPERATOR) {
op = p1->token.operatr;
b1 = p1 + 1;
if (op == DIVIDE) {
if (!g_of_f(op, b1, trhs, &n_trhs, tlhs, &n_tlhs))
return false;
goto simp_again;
}
}
}
} else {
for (b1 = p1 = trhs; p1 < ep; p1++) {
if (p1->kind == VARIABLE && v == p1->token.variable) {
if (op == 0) {
for (p1++;; p1++) {
if (p1 >= ep) {
op = PLUS;
break;
}
if (p1->level == 1 && p1->kind == OPERATOR) {
switch (p1->token.operatr) {
case TIMES:
case DIVIDE:
op = TIMES;
break;
case PLUS:
case MINUS:
op = PLUS;
break;
default:
op = p1->token.operatr;
break;
}
break;
}
}
}
switch (op) {
case TIMES:
case DIVIDE:
case POWER:
b1 = trhs;
op = PLUS;
for (p1 = b1; p1 < ep; p1++)
p1->level++;
break;
}
if (!g_of_f(op, b1, trhs, &n_trhs, tlhs, &n_tlhs))
return false;
goto simp_again;
} else if (p1->level == 1 && p1->kind == OPERATOR) {
op = p1->token.operatr;
b1 = p1 + 1;
}
}
}
if (uf_flag) {
simps_side(trhs, &n_trhs, zsolve);
}
left_again:
worked = true;
uf_flag = false;
see_work:
/* See if we have solved the equation. */
if (se_compare(wantp, wantn, tlhs, n_tlhs, &diff_sign)) {
if (!diff_sign) {
if (!zsolve) {
ep = &trhs[n_trhs];
for (p1 = trhs; p1 < ep; p1 += 2) {
if (p1->kind == VARIABLE && p1->token.variable == v) {
return false;
}
}
} else {
zero_simp:
list_tdebug(4);
uf_power(trhs, &n_trhs);
do {
do {
simp_ssub(trhs, &n_trhs, 0L, 0.0, false, true, 0);
} while (uf_power(trhs, &n_trhs));
} while (super_factor(trhs, &n_trhs, 1));
list_tdebug(3);
ep = &trhs[n_trhs];
op = 0;
for (p1 = &trhs[1]; p1 < ep; p1 += 2) {
if (p1->level == 1) {
op = p1->token.operatr;
if (op == DIVIDE) {
goto no_simp;
}
if (op != TIMES) {
break;
}
}
}
switch (op) {
case TIMES:
for (p1 = trhs; p1 < ep;) {
b1 = p1;
for (;; p1++) {
if (p1 >= ep || (p1->kind == OPERATOR && p1->level == 1)) {
blt(b1 + 1, p1, (char *) ep - (char *) p1);
n_trhs -= p1 - (b1 + 1);
b1->kind = CONSTANT;
b1->token.constant = 1.0;
goto zero_simp;
}
if (p1->kind != CONSTANT && p1->kind != OPERATOR
&& (p1->kind != VARIABLE || (p1->token.variable & VAR_MASK) > SIGN)) {
break;
}
}
p1 = b1;
for (p1++; p1 < ep && p1->level > 1; p1 += 2)
;
p1++;
}
break;
case POWER:
if ((p1 + 1)->level == 1
&& (p1 + 1)->kind == CONSTANT
&& (p1 + 1)->token.constant > 0.0) {
n_trhs -= 2;
goto zero_simp;
}
break;
}
}
if (zsolve) {
debug_string(1, _("Solve for zero completed:"));
} else {
debug_string(1, _("Solve completed:"));
}
list_tdebug(1);
blt(leftp, tlhs, n_tlhs * sizeof(*leftp));
*leftnp = n_tlhs;
blt(rightp, trhs, n_trhs * sizeof(*rightp));
*rightnp = n_trhs;
return success;
}
}
/* Move what we don't want in the LHS to the RHS. */
found_count = 0;
need_flip = 0;
found = 0;
op = 0;
ep = &tlhs[n_tlhs];
for (b1 = p1 = tlhs;; p1++) {
if (p1 >= ep || (p1->level == 1 && p1->kind == OPERATOR)) {
if (!found) {
if ((p1 < ep || found_count || zsolve || n_tlhs > 1 || tlhs[0].kind != CONSTANT)
&& (p1 - b1 != 1 || b1->kind != CONSTANT || b1->token.constant != 1.0
|| p1 >= ep || p1->token.operatr != DIVIDE)) {
if (op == 0) {
for (;; p1++) {
if (p1 >= ep) {
op = PLUS;
break;
}
if (p1->level == 1 && p1->kind == OPERATOR) {
switch (p1->token.operatr) {
case TIMES:
case DIVIDE:
op = TIMES;
break;
case PLUS:
case MINUS:
op = PLUS;
break;
default:
op = p1->token.operatr;
break;
}
break;
}
}
}
if (zsolve) {
if (p1 < ep) {
switch (op) {
case PLUS:
case MINUS:
case DIVIDE:
break;
default:
goto fin1;
}
} else {
if (op != DIVIDE) {
b1 = tlhs;
op = PLUS;
for (p1 = b1; p1 < ep; p1++)
p1->level++;
}
}
}
if (!g_of_f(op, b1, tlhs, &n_tlhs, trhs, &n_trhs))
return false;
list_tdebug(2);
if (uf_flag) {
simp_loop(tlhs, &n_tlhs);
} else {
simps_side(tlhs, &n_tlhs, zsolve);
}
simps_side(trhs, &n_trhs, zsolve);
list_tdebug(1);
goto see_work;
}
} else if (op == DIVIDE) {
need_flip += found;
}
if (p1 >= ep) {
if (found_count == 0) { /* if solve variable no longer in LHS */
for (i = 0; i < n_trhs; i += 2) {
if (trhs[i].kind == VARIABLE && trhs[i].token.variable == v) {
/* solve variable in RHS */
return false;
}
}
calc_simp(tlhs, &n_tlhs);
calc_simp(trhs, &n_trhs);
if (se_compare(tlhs, n_tlhs, trhs, n_trhs, &diff_sign) && !diff_sign) {
error(_("This equation is an identity."));
#if !SILENT
printf(_("That is, the LHS is identical to the RHS.\n"));
#endif
return -1;
}
found = false;
for (i = 0; i < n_tlhs; i += 2) {
if (tlhs[i].kind == VARIABLE
&& tlhs[i].token.variable > IMAGINARY) {
found = true;
break;
}
}
for (i = 0; i < n_trhs; i += 2) {
if (trhs[i].kind == VARIABLE
&& trhs[i].token.variable > IMAGINARY) {
found = true;
break;
}
}
if (found) {
error(_("This equation is independent of the solve variable."));
} else {
error(_("There are no possible values for the solve variable."));
}
return -2;
}
zflag = (n_trhs == 1 && trhs[0].kind == CONSTANT && trhs[0].token.constant == 0.0);
if (zflag) {
/* overwrite -0.0 */
trhs[0].token.constant = 0.0;
}
if (need_flip >= found_count) {
if (!flip(tlhs, &n_tlhs, trhs, &n_trhs))
return false;
list_tdebug(2);
simps_side(tlhs, &n_tlhs, zsolve);
simps_side(trhs, &n_trhs, zsolve);
list_tdebug(1);
goto left_again;
}
if (worked && !uf_flag) {
worked = false;
debug_string(1, _("Unfactoring..."));
partial_flag = false;
uf_simp(tlhs, &n_tlhs);
partial_flag = true;
factorv(tlhs, &n_tlhs, v);
list_tdebug(1);
uf_flag = true;
goto see_work;
}
if (uf_flag) {
simps_side(tlhs, &n_tlhs, zsolve);
list_tdebug(1);
uf_flag = false;
goto see_work;
}
op = 0;
b1 = tlhs;
for (i = 1; i < n_tlhs; i += 2) {
if (tlhs[i].level == 1) {
op_kind = tlhs[i].token.operatr;
if (op_kind == TIMES || op_kind == DIVIDE) {
if (op == 0) {
op = TIMES;
}
} else {
op = op_kind;
break;
}
if (zflag) {
if (op_kind == DIVIDE
|| (tlhs[i+1].kind == VARIABLE && tlhs[i+1].token.variable == v
&& (tlhs[i+1].level == 1
|| (tlhs[i+1].level == 2 && tlhs[i+2].token.operatr == POWER
&& tlhs[i+3].level == 2 && tlhs[i+3].kind == CONSTANT && tlhs[i+3].token.constant > 0.0)))) {
op = op_kind;
b1 = &tlhs[i+1];
if (op_kind == DIVIDE)
break;
}
} else {
if (op_kind == DIVIDE) {
for (j = i + 2; j < n_tlhs && tlhs[j].level > 1; j += 2) {
if (tlhs[j].level == 2) {
op_kind = tlhs[j].token.operatr;
if (op_kind == PLUS || op_kind == MINUS) {
op = DIVIDE;
b1 = &tlhs[i+1];
}
break;
}
}
}
}
}
}
if ((zflag && zero_solved && op == TIMES
&& b1[0].kind == VARIABLE && b1[0].token.variable == v
&& (b1[0].level == 1 || (b1[0].level == 2 && b1[1].token.operatr == POWER
&& b1[2].level == 2 && b1[2].kind == CONSTANT && b1[2].token.constant > 0.0)))
|| op == DIVIDE) {
if (op == TIMES) {
qtries = 0; /* might be quadratic after removing solution */
success = 2;
#if !SILENT
printf(_("Removing possible solution: \""));
list_proc(b1, 1, false);
printf(" = 0\".\n");
#endif
} else {
debug_string(1, _("Juggling..."));
uf_flag = true;
}
if (!g_of_f(op, b1, tlhs, &n_tlhs, trhs, &n_trhs))
return false;
goto simp_again;
}
b1 = NULL;
for (i = 1; i < n_tlhs; i += 2) {
if (tlhs[i].token.operatr == POWER
&& tlhs[i+1].level == tlhs[i].level
&& tlhs[i+1].kind == CONSTANT
&& fabs(tlhs[i+1].token.constant) < 1.0) {
if (!f_to_fraction(tlhs[i+1].token.constant, &numerator, &denominator)
|| fabs(numerator) != 1.0 || denominator < 2.0) {
continue;
}
for (j = i - 1; j >= 0 && tlhs[j].level >= tlhs[i].level; j--) {
if (tlhs[j].kind == VARIABLE && tlhs[j].token.variable == v) {
if (b1) {
if (fabs(b1->token.constant) < fabs(tlhs[i+1].token.constant)) {
b1 = &tlhs[i+1];
}
} else {
b1 = &tlhs[i+1];
}
break;
}
}
}
}
if (b1 && zero_solved) {
inc_count++;
if (inc_count > MAX_RAISE_POWER)
return false;
#if !SILENT
printf(_("Raising both sides to the power of %g and unfactoring...\n"), 1.0 / b1->token.constant);
#endif
zero_solved = false;
qtries = 0;
if (!increase(b1->token.constant, v)) {
error(_("Failed to raise both sides to a power."));
return false;
}
uf_flag = true;
goto simp_again;
}
if (qtries) {
return false;
}
debug_string(1, _("Solving for zero..."));
*leftnp = n_tlhs;
blt(leftp, tlhs, n_tlhs * sizeof(*leftp));
*rightnp = n_trhs;
blt(rightp, trhs, n_trhs * sizeof(*rightp));
if (solve_sub(&zero_token, 1, leftp, leftnp, rightp, rightnp) <= 0)
return false;
if (zero_solved) {
qtries++;
}
zero_solved = true;
if (quad_solve(v)) {
goto left_again;
} else {
goto simp_again;
}
} else {
fin1:
found = 0;
op = p1->token.operatr;
b1 = p1 + 1;
}
} else if (p1->kind == VARIABLE) {
if (v == p1->token.variable) {
found_count++;
found++;
}
}
}
}
/*
* Isolate the expression containing variable "v" raised to the power of "d",
* then raise both sides of the equation to the power of 1/d.
*
* Return true if successful.
*/
static int
increase(d, v)
double d;
long v;
{
int flag, foundp, found2;
int len1, len2;
int op;
token_type *b1, *p1, *p2;
token_type *ep;
partial_flag = false;
ufactor(tlhs, &n_tlhs);
partial_flag = true;
simp_ssub(tlhs, &n_tlhs, v, d, true, false, 2);
simp_ssub(tlhs, &n_tlhs, 0L, 0.0, true, true, 2);
lonely:
ep = &tlhs[n_tlhs];
len1 = 0;
len2 = 0;
foundp = false;
for (p1 = &tlhs[1];; p1 += 2) {
if (p1 >= ep) {
return false;
}
if (p1->level == 1) {
break;
}
if (p1->token.operatr == POWER
&& (p1 + 1)->level == p1->level
&& (p1 + 1)->kind == CONSTANT
&& (p1 + 1)->token.constant == d) {
flag = false;
for (b1 = p1 - 1;; b1--) {
if (b1->level < p1->level) {
b1++;
break;
}
if (b1->kind == VARIABLE && b1->token.variable == v) {
flag = true;
}
if (b1 == tlhs)
break;
}
if (flag) {
foundp = true;
len1 = max(len1, p1 - b1);
}
}
}
found2 = false;
for (p2 = p1 + 2;; p2 += 2) {
if (p2 >= ep) {
break;
}
if (p2->token.operatr == POWER
&& (p2 + 1)->level == p2->level
&& (p2 + 1)->kind == CONSTANT
&& (p2 + 1)->token.constant == d) {
flag = false;
for (b1 = p2 - 1;; b1--) {
if (b1->level < p2->level) {
b1++;
break;
}
if (b1->kind == VARIABLE && b1->token.variable == v) {
flag = true;
}
if (b1 == tlhs)
break;
}
if (flag) {
found2 = true;
len2 = max(len2, p2 - b1);
}
}
}
if (foundp && found2) {
if (len2 > len1)
foundp = false;
}
b1 = p1 + 1;
op = p1->token.operatr;
if (op == POWER && b1->level == 1 && b1->kind == CONSTANT && b1->token.constant == d) {
return(g_of_f(POWER, b1, tlhs, &n_tlhs, trhs, &n_trhs));
}
if (!foundp) {
b1 = tlhs;
if (p1 - b1 == 1 && p1->token.operatr == DIVIDE
&& b1->kind == CONSTANT && b1->token.constant == 1.0) {
if (!flip(tlhs, &n_tlhs, trhs, &n_trhs))
return false;
goto end;
}
switch (p1->token.operatr) {
case TIMES:
case DIVIDE:
op = TIMES;
break;
case PLUS:
case MINUS:
op = PLUS;
break;
default:
op = p1->token.operatr;
break;
}
}
if (!g_of_f(op, b1, tlhs, &n_tlhs, trhs, &n_trhs))
return false;
end:
list_tdebug(2);
simp_loop(tlhs, &n_tlhs);
simp_loop(trhs, &n_trhs);
list_tdebug(1);
goto lonely;
}
/*
* Quadratic and biquadratic solve routine.
* Solves any equation of the form "0 = ax^(2n)+bx^n+c" for "x^n",
* where "x" is an expression containing the solve variable,
* and "n" is a constant. Uses the quadratic formula.
*
* The equation to solve is in tlhs and trhs, it must already be solved for zero.
*
* Return true if successful, with solved equation in tlhs and trhs.
*/
static int
quad_solve(v)
long v; /* solve variable */
{
int i, j, k;
token_type *p1, *p2, *ep;
token_type *x1p = NULL, *x2p;
token_type *a1p = NULL, *a2p = NULL, *a2ep = NULL;
token_type *b1p, *b2p, *b2ep;
token_type *x1tp, *a1tp;
token_type x1_storage[100];
int op, op2, opx1, opx2;
int found, diff_sign;
int len, alen, blen, aloc, nx1;
double high_power = 0.0;
uf_simp(trhs, &n_trhs);
while (factor_plus(trhs, &n_trhs, v, 0.0)) {
simp_loop(trhs, &n_trhs);
}
list_tdebug(1);
found = false;
op = 0;
ep = &trhs[n_trhs];
for (x1tp = p1 = trhs;; p1++) {
if (p1 >= ep || (p1->level == 1 && p1->kind == OPERATOR)) {
if (p1 < ep) {
switch (p1->token.operatr) {
case PLUS:
case MINUS:
break;
default:
return false;
}
}
if (op == TIMES || op == DIVIDE) {
found = false;
op2 = 0;
for (a1tp = p2 = x1tp;; p2++) {
if (p2 >= p1)
break;
if (p2->level == 2) {
if (p2->kind == OPERATOR) {
x1tp = p2 + 1;
op2 = p2->token.operatr;
found = false;
}
} else {
if (p2->kind == OPERATOR) {
if (p2->level == 3 && p2->token.operatr == POWER) {
if (found && (op2 == TIMES || op2 == 0)
&& (p2 + 1)->level == 3
&& (p2 + 1)->kind == CONSTANT
&& (p2 + 1)->token.constant > high_power) {
high_power = (p2 + 1)->token.constant;
x1p = x1tp;
a1p = a1tp;
a2p = p2 + 2;
a2ep = p1;
}
}
} else if (p2->kind == VARIABLE && p2->token.variable == v) {
found = true;
}
}
}
} else if (op == POWER && found && (p1 - 1)->level == 2
&& (p1 - 1)->kind == CONSTANT && (p1 - 1)->token.constant > high_power) {
high_power = (p1 - 1)->token.constant;
a1p = x1p = x1tp;
a2p = p1;
a2ep = a2p;
}
if (p1 >= ep) {
break;
}
}
if (p1->level == 1) {
if (p1->kind == OPERATOR) {
op = 0;
x1tp = p1 + 1;
found = false;
}
} else {
if (p1->kind == OPERATOR) {
if (p1->level == 2) {
op = p1->token.operatr;
}
} else if (op == 0 && p1->kind == VARIABLE && p1->token.variable == v) {
found = true;
}
}
}
if (high_power == 0.0)
return false;
#if !SILENT
list_var(v, 0);
printf(_("Equation is a degree %.*g polynomial in %s.\n"), precision, high_power, var_str);
#endif
if (a1p > trhs && (a1p - 1)->token.operatr == MINUS)
opx1 = MINUS;
else
opx1 = PLUS;
if (high_power == 2.0) {
nx1 = (a2p - x1p) - 2;
if (nx1 > ARR_CNT(x1_storage))
return false;
blt(x1_storage, x1p, nx1 * sizeof(token_type));
} else {
nx1 = (a2p - x1p);
if (nx1 > ARR_CNT(x1_storage))
return false;
blt(x1_storage, x1p, nx1 * sizeof(token_type));
x1_storage[nx1-1].token.constant /= 2.0;
}
opx2 = 0;
op = 0;
for (x2p = p1 = trhs;; p1++) {
if (p1 >= ep || (p1->level == 1 && p1->kind == OPERATOR)) {
if (se_compare(x1_storage, nx1, x2p, p1 - x2p, &diff_sign)) {
b1p = x2p;
b2p = p1;
b2ep = b2p;
break;
}
if (op == TIMES || op == DIVIDE) {
op2 = 0;
for (b1p = p2 = x2p;; p2++) {
if (p2 >= p1 || (p2->level == 2 && p2->kind == OPERATOR)) {
if (op2 == 0 || op2 == TIMES) {
if (se_compare(x1_storage, nx1, x2p, p2 - x2p, &diff_sign)) {
b2p = p2;
b2ep = p1;
goto big_bbreak;
}
}
if (p2 >= p1)
break;
}
if (p2->level == 2 && p2->kind == OPERATOR) {
x2p = p2 + 1;
op2 = p2->token.operatr;
}
}
}
if (p1 >= ep)
return false;
}
if (p1->level == 1) {
if (p1->kind == OPERATOR) {
op = 0;
opx2 = p1->token.operatr;
x2p = p1 + 1;
}
} else {
if (p1->kind == OPERATOR && p1->level == 2) {
op = p1->token.operatr;
}
}
}
big_bbreak:
switch (opx2) {
case 0:
opx2 = PLUS;
case PLUS:
if (diff_sign)
opx2 = MINUS;
break;
case MINUS:
if (diff_sign)
opx2 = PLUS;
break;
default:
return false;
}
blt(scratch, b1p, (char *) x2p - (char *) b1p);
len = x2p - b1p;
scratch[len].level = 7;
scratch[len].kind = CONSTANT;
if (opx2 == MINUS)
scratch[len].token.constant = -1.0;
else
scratch[len].token.constant = 1.0;
len++;
blt(&scratch[len], b2p, (char *) b2ep - (char *) b2p);
len += (b2ep - b2p);
blen = len;
j = min_level(scratch, len);
j = 7 - j;
for (i = 0; i < len; i++)
scratch[i].level += j;
scratch[len].level = 6;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = POWER;
len++;
scratch[len].level = 6;
scratch[len].kind = CONSTANT;
scratch[len].token.constant = 2.0;
len++;
scratch[len].level = 5;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = MINUS;
len++;
scratch[len].level = 6;
scratch[len].kind = CONSTANT;
scratch[len].token.constant = 4.0;
len++;
scratch[len].level = 6;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = TIMES;
len++;
aloc = len;
blt(&scratch[len], a1p, (char *) x1p - (char *) a1p);
len += (x1p - a1p);
scratch[len].level = 7;
scratch[len].kind = CONSTANT;
if (opx1 == MINUS)
scratch[len].token.constant = -1.0;
else
scratch[len].token.constant = 1.0;
len++;
blt(&scratch[len], a2p, (char *) a2ep - (char *) a2p);
len += (a2ep - a2p);
alen = len - aloc;
j = min_level(&scratch[aloc], len - aloc);
j = 7 - j;
for (i = aloc; i < len; i++)
scratch[i].level += j;
scratch[len].level = 6;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = TIMES;
len++;
k = len;
scratch[len] = zero_token;
len++;
for (p2 = p1 = trhs;; p1++) {
if (p1 >= ep || (p1->level == 1 && p1->kind == OPERATOR)) {
if (!((p2 <= x1p && p1 > x1p) || (p2 <= x2p && p1 > x2p))) {
if (p2 == trhs) {
scratch[len].level = 1;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = PLUS;
len++;
}
blt(&scratch[len], p2, (char *) p1 - (char *) p2);
len += (p1 - p2);
}
if (p1 >= ep)
break;
else
p2 = p1;
}
}
for (i = k; i < len; i++)
scratch[i].level += 6;
scratch[len].level = 4;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = POWER;
len++;
scratch[len].level = 4;
scratch[len].kind = CONSTANT;
scratch[len].token.constant = 0.5;
len++;
scratch[len].level = 3;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = TIMES;
len++;
scratch[len].level = 3;
scratch[len].kind = VARIABLE;
next_sign(&scratch[len].token.variable);
len++;
scratch[len].level = 2;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = MINUS;
len++;
if (len + blen + 3 + alen > n_tokens) {
error_huge();
}
blt(&scratch[len], scratch, blen * sizeof(token_type));
len += blen;
scratch[len].level = 1;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = DIVIDE;
len++;
scratch[len].level = 2;
scratch[len].kind = CONSTANT;
scratch[len].token.constant = 2.0;
len++;
scratch[len].level = 2;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = TIMES;
len++;
blt(&scratch[len], &scratch[aloc], alen * sizeof(token_type));
len += alen;
if (found_var(scratch, len, v))
return false;
debug_string(1, _("Plugging equation into the quadratic formula:"));
blt(tlhs, x1_storage, nx1 * sizeof(token_type));
n_tlhs = nx1;
simp_loop(tlhs, &n_tlhs);
blt(trhs, scratch, len * sizeof(token_type));
n_trhs = len;
simp_loop(trhs, &n_trhs);
list_tdebug(2);
uf_simp(trhs, &n_trhs);
simps_side(trhs, &n_trhs, false);
list_tdebug(1);
#if !SILENT
if (high_power == 2.0) {
printf(_("Equation was quadratic.\n"));
} else if (high_power == 4.0) {
printf(_("Equation was biquadratic.\n"));
} else {
printf(_("Equation was solved with the quadratic formula.\n"));
}
#endif
return true;
}
/*
* This is the heart of Mathomatic solving:
* It applies an identical mathematical operation to both sides of an equation.
*
* Apply the inverse of the operation "op" followed by expression "operandp",
* which is somewhere in "side1p", to both sides of an equation,
* which is "side1p" and "side2p".
*
* Return true unless something is wrong.
*/
static int
g_of_f(op, operandp, side1p, side1np, side2p, side2np)
int op; /* current operator */
token_type *operandp; /* operand pointer */
token_type *side1p; /* equation side pointer */
int *side1np; /* pointer to the length of "side1p" */
token_type *side2p; /* equation side pointer */
int *side2np; /* pointer to the length of "side2p" */
{
token_type *p1, *p2, *ep;
int oldn, operandn;
double numerator, denominator;
static int prev_n1, prev_n2;
double d1, d2;
oldn = *side1np;
ep = &side1p[oldn];
if (*side1np == prev_n1 && *side2np == prev_n2) {
if (++repeat_count >= 5) {
debug_string(1, _("Infinite loop aborted in solve routine."));
return false;
}
} else {
prev_n1 = *side1np;
prev_n2 = *side2np;
repeat_count = 0;
}
switch (op) {
case PLUS:
case MINUS:
case TIMES:
case DIVIDE:
case POWER:
case MODULUS:
break;
default:
return false;
}
for (p1 = operandp + 1; p1 < ep; p1 += 2) {
if (p1->level == 1) {
switch (p1->token.operatr) {
case FACTORIAL:
op = PLUS;
continue;
case MODULUS:
operandp = p1 + 1;
continue;
}
break;
}
}
operandn = p1 - operandp;
if (op == POWER && operandp == side1p) {
if (!get_constant(side2p, *side2np, &d1))
return false;
if (!get_constant(operandp, operandn, &d2))
return false;
debug_string(1, _("Taking logarithm of both sides:"));
*side2np = 1;
side2p->level = 1;
side2p->kind = CONSTANT;
errno = 0;
side2p->token.constant = log(d1) / log(d2);
check_err();
blt(side1p, p1 + 1, (*side1np - (operandn + 1)) * sizeof(token_type));
*side1np -= operandn + 1;
return true;
}
#if !SILENT
if (debug_level > 0) {
switch (op) {
case PLUS:
printf(_("Subtracting"));
break;
case MINUS:
printf(_("Adding"));
break;
case TIMES:
printf(_("Dividing both sides of the equation by"));
break;
case DIVIDE:
printf(_("Multiplying both sides of the equation by"));
break;
case POWER:
printf(_("Raising both sides of the equation to the power of"));
break;
case MODULUS:
printf(_("Applying inverse modulus of"));
break;
}
printf(" \"");
if (op == POWER)
printf("1/(");
list_proc(operandp, operandn, false);
switch (op) {
case PLUS:
printf(_("\" from both sides of the equation:\n"));
break;
case MINUS:
case MODULUS:
printf(_("\" to both sides of the equation:\n"));
break;
case POWER:
printf(")");
default:
printf("\":\n");
break;
}
}
#endif
if (*side1np + operandn + 3 > n_tokens || *side2np + operandn + 5 > n_tokens) {
error_huge();
}
for (p2 = side1p; p2 < ep; p2++)
p2->level++;
ep = &side2p[*side2np];
for (p2 = side2p; p2 < ep; p2++)
p2->level++;
p2 = &side1p[oldn];
switch (op) {
case MODULUS:
p2->level = 1;
p2->kind = OPERATOR;
p2->token.operatr = PLUS;
p2++;
p2->level = 2;
p2->kind = VARIABLE;
parse_var(&p2->token.variable, "integer");
p2++;
p2->level = 2;
p2->kind = OPERATOR;
p2->token.operatr = TIMES;
p2++;
blt(p2, operandp, (char *) p1 - (char *) operandp);
*side1np += 3 + operandn;
break;
case POWER:
p2->level = 1;
p2->kind = OPERATOR;
p2->token.operatr = POWER;
p2++;
p2->level = 2;
p2->kind = CONSTANT;
p2->token.constant = 1.0;
p2++;
p2->level = 2;
p2->kind = OPERATOR;
p2->token.operatr = DIVIDE;
p2++;
blt(p2, operandp, (char *) p1 - (char *) operandp);
*side1np += 3 + operandn;
break;
case TIMES:
p2->level = 1;
p2->kind = OPERATOR;
p2->token.operatr = DIVIDE;
p2++;
blt(p2, operandp, (char *) p1 - (char *) operandp);
*side1np += 1 + operandn;
break;
case DIVIDE:
p2->level = 1;
p2->kind = OPERATOR;
p2->token.operatr = TIMES;
p2++;
blt(p2, operandp, (char *) p1 - (char *) operandp);
*side1np += 1 + operandn;
break;
case PLUS:
p2->level = 1;
p2->kind = OPERATOR;
p2->token.operatr = MINUS;
p2++;
blt(p2, operandp, (char *) p1 - (char *) operandp);
*side1np += 1 + operandn;
break;
case MINUS:
p2->level = 1;
p2->kind = OPERATOR;
p2->token.operatr = PLUS;
p2++;
blt(p2, operandp, (char *) p1 - (char *) operandp);
*side1np += 1 + operandn;
break;
}
blt(&side2p[*side2np], &side1p[oldn], (*side1np - oldn) * sizeof(*side1p));
*side2np += *side1np - oldn;
if (op == POWER && operandn == 1 && operandp->kind == CONSTANT) {
f_to_fraction(operandp->token.constant, &numerator, &denominator);
if (always_positive(numerator)) {
ep = &side2p[*side2np];
for (p2 = side2p; p2 < ep; p2++)
p2->level++;
p2->level = 1;
p2->kind = OPERATOR;
p2->token.operatr = TIMES;
p2++;
p2->level = 1;
p2->kind = VARIABLE;
next_sign(&p2->token.variable);
*side2np += 2;
}
}
if (op == POWER || op == MODULUS) {
*side1np = (operandp - 1) - side1p;
}
return true;
}
/*
* Take the reciprocal of both equation sides.
*
* Return true if successful.
*/
static int
flip(side1p, side1np, side2p, side2np)
token_type *side1p; /* equation side pointer */
int *side1np; /* pointer to equation side length */
token_type *side2p;
int *side2np;
{
token_type *p1, *ep;
debug_string(1, _("Taking the reciprocal of both sides of the equation..."));
if (*side1np + 2 > n_tokens || *side2np + 2 > n_tokens) {
error_huge();
}
ep = &side1p[*side1np];
for (p1 = side1p; p1 < ep; p1++)
p1->level++;
ep = &side2p[*side2np];
for (p1 = side2p; p1 < ep; p1++)
p1->level++;
blt(side1p + 2, side1p, *side1np * sizeof(*side1p));
*side1np += 2;
blt(side2p + 2, side2p, *side2np * sizeof(*side2p));
*side2np += 2;
*side1p = one_token;
side1p++;
side1p->level = 1;
side1p->kind = OPERATOR;
side1p->token.operatr = DIVIDE;
*side2p = one_token;
side2p++;
side2p->level = 1;
side2p->kind = OPERATOR;
side2p->token.operatr = DIVIDE;
return true;
}
syntax highlighted by Code2HTML, v. 0.9.1