Logo Search packages:      
Sourcecode: mathomatic version File versions  Download package

solve.c

/*
 * Algebraic manipulator solve routines.
 *
 * Copyright (c) 1987-2005 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) {
            printf(_("Solve failed.\n"));
      }
#endif
      return rv;
}

/*
 * Main Mathomatic 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.
 *
 * Returns true if successful, with the result placed in the passed LHS and RHS.
 * Returns 2 if a solution was zero and removed.
 * Returns false on failure.
 */
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 */
int         *leftnp;    /* size of LHS */
token_type  *rightp;    /* RHS */
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, want0;
      long        v;
      int         need_flip;
      int         uf_flag;
      int         qtries;
      int         zflag;
      int         zsolve;
      int         inc_count;
      int         zero_solved;
      double            numerator, denominator;
      int         success;

      success = 1;
      repeat_count = 0;
      qtries = 0;
      zero_solved = false;
      uf_flag = false;
      inc_count = 0;
      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."));
            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;
      }
#if   true
      uf_power(tlhs, &n_tlhs);
      uf_power(trhs, &n_trhs);
#endif
simp_again:
      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:
      /* Move 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;
                                    }
                              }
                        }
                        if (op == TIMES || op == DIVIDE || op == POWER) {
                              b1 = trhs;
                              op = PLUS;
                              for (p1 = b1; p1 < ep; p1++)
                                    p1->level++;
                        }
                        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 = 1;
      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, true));
                        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 left equation side is identical to the right equation side.\n"));
#endif
                                    return false;
                              }
                              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 false;
                        }
                        zflag = (n_trhs == 1 && trhs[0].kind == CONSTANT && trhs[0].token.constant == 0.0);
                        if (need_flip >= found_count) {
                              if (zflag) {
                                    /* overwrite -0.0 */
                                    trhs[0].token.constant = 0.0;
                              }
                              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--;
                              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(_("increase() function failed."));
                                    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));
                        want0.level = 1;
                        want0.kind = CONSTANT;
                        want0.token.constant = 0.0;
                        if (!solve_sub(&want0, 1, leftp, leftnp, rightp, rightnp))
                              return false;
                        if (zero_solved) {
                              qtries++;
                        } else {
                              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",
 * where "x" is an expression containing passed variable "v",
 * and "n" is a constant.
 *
 * Return true if successful.
 */
static int
quad_solve(v)
long  v;
{
      int         i, j, k;
      token_type  *p1, *p2, *ep;
      token_type  *x1p, *x2p;
      token_type  *a1p, *a2p, *a2ep;
      token_type  *b1p, *b2p, *b2ep;
      token_type  *x1tp, *a1tp;
      token_type  x1_storage[100];
      int         nx1;
      int         op, op2, opx1, opx2;
      int         found;
      int         diff_sign;
      int         len;
      int         alen, blen;
      int         aloc;
      double            high_power;

      uf_simp(trhs, &n_trhs);
      factorv(trhs, &n_trhs, v);
      list_tdebug(1);

      found = false;
      op = 0;
      high_power = 0.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) {
                        if (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) {
                        if (high_power == 0.0)
                              return false;
                        else
                              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) {
                        if (p1->kind == VARIABLE && p1->token.variable == v)
                              found = true;
                  }
            }
      }
      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) {
                                    if (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) {
                        if (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, (int) ((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, (int) ((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, (int) ((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, (int) ((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].level = 1;
      scratch[len].kind = CONSTANT;
      scratch[len].token.constant = 0.0;
      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, (int) ((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 a degree %g quadratic.\n"), high_power);
      }
#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" */
{
      int         i;
      token_type  *p1, *p2, *ep;
      int         oldn;
      int         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;
      case FACTORIAL:
            i = oldn - 2;
            if (i >= 0 && side1p[i].level == 1 && side1p[i].token.operatr == op) {
                  i = *side2np - 2;
                  if (i >= 0 && side2p[i].level == 1 && side2p[i].token.operatr == op) {
                        debug_string(1, _("Removing factorial from both sides of the equation:"));
                        *side1np -= 2;
                        *side2np -= 2;
                        return true;
                  }
            }
      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:
                        if (op != MODULUS) {
                              printf(_("Something is wrong with the solve routine!\n"));
                        } else {
                              operandp = p1 + 1;
                        }
                        continue;
                  }
                  break;
            }
      }
      operandn = p1 - operandp;
      if (op == MODULUS) {
            if (*side2np == 1 && side2p->kind == CONSTANT) {
                  if (side2p->token.constant < 0.0
                      || (operandn == 1 && operandp->kind == CONSTANT
                      && side2p->token.constant >= fabs(operandp->token.constant))) {
                        error(_("There are no solutions."));
                        return false;
                  }
            }
      }
      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;
            p2->token.variable = V_INTEGER;
            p2++;
            p2->level = 2;
            p2->kind = OPERATOR;
            p2->token.operatr = TIMES;
            p2++;
            blt(p2, operandp, (int) ((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, (int) ((char *) p1 - (char *) operandp));
            *side1np += 3 + operandn;
            break;
      case TIMES:
            p2->level = 1;
            p2->kind = OPERATOR;
            p2->token.operatr = DIVIDE;
            p2++;
            blt(p2, operandp, (int) ((char *) p1 - (char *) operandp));
            *side1np += 1 + operandn;
            break;
      case DIVIDE:
            p2->level = 1;
            p2->kind = OPERATOR;
            p2->token.operatr = TIMES;
            p2++;
            blt(p2, operandp, (int) ((char *) p1 - (char *) operandp));
            *side1np += 1 + operandn;
            break;
      case PLUS:
            p2->level = 1;
            p2->kind = OPERATOR;
            p2->token.operatr = MINUS;
            p2++;
            blt(p2, operandp, (int) ((char *) p1 - (char *) operandp));
            *side1np += 1 + operandn;
            break;
      case MINUS:
            p2->level = 1;
            p2->kind = OPERATOR;
            p2->token.operatr = PLUS;
            p2++;
            blt(p2, operandp, (int) ((char *) p1 - (char *) operandp));
            *side1np += 1 + operandn;
            break;
      }
      blt(&side2p[*side2np], &side1p[oldn], (int) (*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->level = 1;
      side1p->kind = CONSTANT;
      side1p->token.constant = 1.0;
      side1p++;
      side1p->level = 1;
      side1p->kind = OPERATOR;
      side1p->token.operatr = DIVIDE;

      side2p->level = 1;
      side2p->kind = CONSTANT;
      side2p->token.constant = 1.0;
      side2p++;
      side2p->level = 1;
      side2p->kind = OPERATOR;
      side2p->token.operatr = DIVIDE;
      return true;
}

Generated by  Doxygen 1.6.0   Back to index