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

integrate.c

/*
 * Algebraic manipulator integration routines.
 *
 * These are very low level, so they don't do much, mostly polynomials.
 *
 * Copyright (c) 1987-2005 George Gesslein II.
 */

#include "includes.h"

static int  int_sub();
static int  laplace_sub();
static int  inv_laplace_sub();

/*
 * Make variable "v" always raised to a power,
 * unless it is on the right side of a power operator.
 */
make_powers(equation, np, v)
token_type  *equation;  /* pointer to expression */
int         *np;        /* pointer to length of expression */
long        v;          /* variable */
{
      int   i;
      int   level;

      for (i = 0; i < *np;) {
            level = equation[i].level;
            if (equation[i].kind == OPERATOR && equation[i].token.operatr == POWER) {
                  for (i += 2; i < *np && equation[i].level >= level; i += 2)
                        ;
                  continue;
            }
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  if ((i + 1) >= *np || equation[i+1].token.operatr != POWER) {
                        if (*np + 2 > n_tokens) {
                              error_huge();
                        }
                        level++;
                        equation[i].level = level;
                        i++;
                        blt(&equation[i+2], &equation[i], (*np - i) * sizeof(token_type));
                        *np += 2;
                        equation[i].level = level;
                        equation[i].kind = OPERATOR;
                        equation[i].token.operatr = POWER;
                        i++;
                        equation[i].level = level;
                        equation[i].kind = CONSTANT;
                        equation[i].token.constant = 1.0;
                  }
            }
            i++;
      }
}

/*
 * Integration dispatch routine.
 *
 * Return true if successful.
 */
int
int_dispatch(equation, np, v, func)
token_type  *equation;
int         *np;
long        v;
int         (*func)();  /* integration function to call */
{
      int   i, j;

      make_powers(equation, np, v);
      if (level1_plus(equation, *np)) {
            for (j = 0, i = 1;; i += 2) {
                  if (i >= *np || equation[i].level == 1) {
                        if (!(*func)(equation, np, j, i, v))
                              return false;
                        for (i = j + 1;; i += 2) {
                              if (i >= *np) {
                                    return true;
                              }
                              if (equation[i].level == 1) {
                                    j = i + 1;
                                    break;
                              }
                        }
                  }
            }
      } else {
            if (!(*func)(equation, np, 0, *np, v))
                  return false;
      }
      return true;
}

/*
 * Do the actual integration.
 *
 * Return true if successful.
 */
static int
int_sub(equation, np, loc, eloc, v)
token_type  *equation;
int         *np;
int         loc, eloc;
long        v;
{
      int         i, j, k;
      int         len;
      int         level, mlevel;
      int         count;
      int         exp_pos;
      complexs    c;
      int         div_flag;

      exp_pos = -1;
      level = min_level(&equation[loc], eloc - loc);
      for (i = loc, count = 0; i < eloc; i += 2) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  count++;
                  if (count > 1)
                        return false;
                  if (i > loc && equation[i].level == equation[i-1].level
                      && equation[i-1].token.operatr == POWER) {
                        exp_pos = i - 1;
                        if (equation[exp_pos].level == level) {
                              continue;
                        } else if (equation[exp_pos].level == level + 1) {
                              for (k = loc + 1; k < eloc; k += 2) {
                                    if (equation[k].level == level) {
                                          switch (equation[k].token.operatr) {
                                          case TIMES:
                                          case DIVIDE:
                                                break;
                                          default:
                                                return false;
                                          }
                                    }
                              }
                              continue;
                        }
                        return false;
                  }
                  if (equation[i].level == level || equation[i].level == (level + 1)) {
                        for (k = loc + 1; k < eloc; k += 2) {
                              if (equation[k].level == level) {
                                    switch (equation[k].token.operatr) {
                                    case DIVIDE:
                                    case TIMES:
                                          continue;
                                    case POWER:
                                          if (k == (i + 1))
                                                continue;
                                    default:
                                          return false;
                                    }
                              }
                        }
                        if (equation[i].level == (level + 1)) {
                              if ((i + 1) < eloc && equation[i+1].level == equation[i].level
                                  && equation[i+1].token.operatr == POWER) {
                                    continue;
                              }
                        } else {
                              continue;
                        }
                  }
                  return false;
            }
      }
      if (exp_pos >= 0) {
            for (i = exp_pos - 2; i > loc && equation[i].level > equation[exp_pos].level; i -= 2)
                  ;
            i++;
            if (!parse_complex(&equation[i], exp_pos - i, &c))
                  return false;
            c = complex_log(c);
            if (i > loc && equation[i-1].token.operatr == DIVIDE) {
                  c.re = -c.re;
                  c.im = -c.im;
            }
            if (*np + 6 > n_tokens) {
                  error_huge();
            }
            for (j = loc; j < eloc; j++)
                  equation[j].level++;
            blt(&equation[eloc+6], &equation[eloc], (*np - eloc) * sizeof(token_type));
            *np += 6;
            k = eloc;
            equation[k].level = level;
            equation[k].kind = OPERATOR;
            equation[k].token.operatr = DIVIDE;
            k++;
            equation[k].level = level + 1;
            equation[k].kind = CONSTANT;
            equation[k].token.constant = c.re;
            k++;
            equation[k].level = level + 1;
            equation[k].kind = OPERATOR;
            equation[k].token.operatr = PLUS;
            k++;
            equation[k].level = level + 2;
            equation[k].kind = CONSTANT;
            equation[k].token.constant = c.im;
            k++;
            equation[k].level = level + 2;
            equation[k].kind = OPERATOR;
            equation[k].token.operatr = TIMES;
            k++;
            equation[k].level = level + 2;
            equation[k].kind = VARIABLE;
            equation[k].token.variable = IMAGINARY;
            return true;
      }
      mlevel = level + 1;
      for (j = loc; j < eloc; j++)
            equation[j].level += 2;
      for (i = loc; i < eloc; i += 2) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  div_flag = (i > loc && equation[i-1].token.operatr == DIVIDE);
                  i++;
                  if (i >= eloc || equation[i].token.operatr != POWER)
                        return false;
                  level = equation[i].level;
                  i++;
                  if (div_flag) {
                        if (equation[i].level == level
                            && equation[i].kind == CONSTANT
                            && equation[i].token.constant == 1.0)
                              return false;
                        if (*np + 2 > n_tokens)
                              error_huge();
                        for (j = i; j < eloc && equation[j].level >= level; j++)
                              equation[j].level++;
                        equation[i-3].token.operatr = TIMES;
                        blt(&equation[i+2], &equation[i], (*np - i) * sizeof(token_type));
                        *np += 2;
                        eloc += 2;
                        equation[i].level = level + 1;
                        equation[i].kind = CONSTANT;
                        equation[i].token.constant = -1.0;
                        equation[i+1].level = level + 1;
                        equation[i+1].kind = OPERATOR;
                        equation[i+1].token.operatr = TIMES;
                  }
                  for (j = i; j < eloc && equation[j].level >= level; j++)
                        equation[j].level++;
                  len = j - i;
                  if (*np + len + 5 > n_tokens)
                        error_huge();
                  blt(&equation[j+2], &equation[j], (*np - j) * sizeof(token_type));
                  *np += 2;
                  eloc += 2;
                  len += 2;
                  level++;
                  equation[j].level = level;
                  equation[j].kind = OPERATOR;
                  equation[j].token.operatr = PLUS;
                  j++;
                  equation[j].level = level;
                  equation[j].kind = CONSTANT;
                  equation[j].token.constant = 1.0;
                  j++;
                  blt(&equation[eloc+len+1], &equation[eloc], (*np - eloc) * sizeof(token_type));
                  *np += len + 1;
                  k = eloc;
                  equation[k].level = mlevel;
                  equation[k].kind = OPERATOR;
                  equation[k].token.operatr = DIVIDE;
                  k++;
                  blt(&equation[k], &equation[i], len * sizeof(token_type));
                  return true;
            }
      }
      if (*np + 2 > n_tokens) {
            error_huge();
      }
      blt(&equation[eloc+2], &equation[eloc], (*np - eloc) * sizeof(token_type));
      *np += 2;
      k = eloc;
      equation[k].level = mlevel;
      equation[k].kind = OPERATOR;
      equation[k].token.operatr = TIMES;
      k++;
      equation[k].level = mlevel;
      equation[k].kind = VARIABLE;
      equation[k].token.variable = v;
      return true;
}

/*
 * The integrate command.
 */
int
integrate_cmd(cp)
char  *cp;
{
      int         i, j;
      long        v;
      double            d1, order;
      token_type  *source, *dest;
      int         n1, n2, *nps, *np;
      int         def_flag;

      order = 1.0;
      if (notdefined(cur_equation)) {
            return false;
      }
      if (def_flag = (strcmp_tospace(cp, "definite") == 0)) {
            cp = skip_param(cp);
      }
      i = next_espace();
      if (n_rhs[cur_equation]) {
            source = rhs[cur_equation];
            nps = &n_rhs[cur_equation];
            dest = rhs[i];
            np = &n_rhs[i];
      } else {
            source = lhs[cur_equation];
            nps = &n_lhs[cur_equation];
            dest = lhs[i];
            np = &n_lhs[i];
      }
      if (*cp) {
            cp = parse_var2(&v, cp);
            if (cp == NULL) {
                  return false;
            }
            if (*cp) {
                  order = strtod(cp, &cp);
            }
            if (*cp || order <= 0 || fmod(order, 1.0) != 0.0) {
                  error(_("The order must be a positive integer."));
                  return false;
            }
      } else {
            if (!prompt_var(&v)) {
                  return false;
            }
      }
      partial_flag = false;
      uf_simp(source, nps);
      partial_flag = true;
      factorv(source, nps, v);
      blt(dest, source, *nps * sizeof(token_type));
      n1 = *nps;
      for (d1 = 0; d1 < order; d1++) {
            if (!int_dispatch(dest, &n1, v, int_sub)) {
                  error(_("Integration failed."));
                  return false;
            }
            simp_loop(dest, &n1);
      }
      if (def_flag) {
            my_strlcpy(prompt_str, _("Enter lower bound: "), sizeof(prompt_str));
            if (!get_expr(tlhs, &n_tlhs)) {
                  return false;
            }
            my_strlcpy(prompt_str, _("Enter upper bound: "), sizeof(prompt_str));
            if (!get_expr(trhs, &n_trhs)) {
                  return false;
            }
            blt(scratch, dest, n1 * sizeof(token_type));
            n2 = n1;
            subst_var_with_exp(scratch, &n2, tlhs, n_tlhs, v);
            subst_var_with_exp(dest, &n1, trhs, n_trhs, v);
            if (n1 + 1 + n2 > n_tokens) {
                  error_huge();
            }
            for (j = 0; j < n1; j++) {
                  dest[j].level++;
            }
            for (j = 0; j < n2; j++) {
                  scratch[j].level++;
            }
            dest[n1].kind = OPERATOR;
            dest[n1].level = 1;
            dest[n1].token.operatr = MINUS;
            n1++;
            blt(&dest[n1], scratch, n2 * sizeof(token_type));
            n1 += n2;
      }
      simpa_side(dest, &n1, false, false);
      if (n_rhs[cur_equation]) {
            blt(lhs[i], lhs[cur_equation], n_lhs[cur_equation] * sizeof(token_type));
            n_lhs[i] = n_lhs[cur_equation];
      }
      *np = n1;
      cur_equation = i;
      return_result(cur_equation);
      return true;
}

static int
laplace_sub(equation, np, loc, eloc, v)
token_type  *equation;
int         *np;
int         loc, eloc;
long        v;
{
      int         i, j, k;
      int         len;
      int         level, level2, mlevel;
      int         count;
      int         exp_pos;
      complexs    c;

      exp_pos = -1;
      level = min_level(&equation[loc], eloc - loc);
      for (i = loc, count = 0; i < eloc; i += 2) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  count++;
                  if (count > 1)
                        return false;
                  if (i > loc && equation[i].level == equation[i-1].level
                      && equation[i-1].token.operatr == POWER) {
                        if (equation[i-2].level != equation[i-1].level) {
                              return false;
                        }
                        exp_pos = i - 2;
                        if (equation[i].level == level) {
                              continue;
                        } else if (equation[i].level == level + 1) {
                              for (k = loc + 1; k < eloc; k += 2) {
                                    if (equation[k].level == level) {
                                          switch (equation[k].token.operatr) {
                                          case TIMES:
                                          case DIVIDE:
                                                break;
                                          default:
                                                return false;
                                          }
                                    }
                              }
                              continue;
                        }
                        return false;
                  }
                  if (equation[i].level == level || equation[i].level == (level + 1)) {
                        for (k = loc + 1; k < eloc; k += 2) {
                              if (equation[k].level == level) {
                                    switch (equation[k].token.operatr) {
                                    case DIVIDE:
                                          if (k == (i - 1))
                                                return false;
                                    case TIMES:
                                          continue;
                                    case POWER:
                                          if (k == (i + 1))
                                                continue;
                                    default:
                                          return false;
                                    }
                              }
                        }
                        if (equation[i].level == (level + 1)) {
                              if ((i + 1) < eloc && equation[i+1].level == equation[i].level
                                  && equation[i+1].token.operatr == POWER) {
                                    continue;
                              }
                        } else {
                              continue;
                        }
                  }
                  return false;
            }
      }
      if (exp_pos >= 0) {
            if (*np + 6 > n_tokens) {
                  error_huge();
            }
            if (!parse_complex(&equation[exp_pos], 1, &c))
                  return false;
            c = complex_log(c);
            if (exp_pos > loc && equation[exp_pos-1].token.operatr == DIVIDE) {
                  equation[exp_pos-1].token.operatr = TIMES;
                  c.re = -c.re;
                  c.im = -c.im;
            }
            equation[exp_pos].kind = CONSTANT;
            equation[exp_pos].token.constant = 1.0;
            exp_pos++;
            level2 = equation[exp_pos].level + 1;
            equation[exp_pos].token.operatr = DIVIDE;
            exp_pos++;
            equation[exp_pos].level = level2;
            exp_pos++;
            blt(&equation[exp_pos+6], &equation[exp_pos], (*np - exp_pos) * sizeof(token_type));
            *np += 6;
            k = exp_pos;
            equation[k].level = level2;
            equation[k].kind = OPERATOR;
            equation[k].token.operatr = MINUS;
            k++;
            equation[k].level = level2 + 1;
            equation[k].kind = CONSTANT;
            equation[k].token.constant = c.re;
            k++;
            equation[k].level = level2 + 1;
            equation[k].kind = OPERATOR;
            equation[k].token.operatr = PLUS;
            k++;
            equation[k].level = level2 + 2;
            equation[k].kind = CONSTANT;
            equation[k].token.constant = c.im;
            k++;
            equation[k].level = level2 + 2;
            equation[k].kind = OPERATOR;
            equation[k].token.operatr = TIMES;
            k++;
            equation[k].level = level2 + 2;
            equation[k].kind = VARIABLE;
            equation[k].token.variable = IMAGINARY;
            return true;
      }
      mlevel = level + 1;
      for (j = loc; j < eloc; j++)
            equation[j].level += 2;
      for (i = loc; i < eloc; i += 2) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  i++;
                  if (i >= eloc || equation[i].token.operatr != POWER)
                        return false;
                  level = equation[i].level;
                  i++;
                  for (j = i; j < eloc && equation[j].level >= level; j++)
                        equation[j].level++;
                  len = j - i;
                  if (*np + len + 7 > n_tokens)
                        error_huge();
                  blt(&equation[j+4], &equation[j], (*np - j) * sizeof(token_type));
                  *np += 4;
                  eloc += 4;
                  level++;
                  equation[j].level = level;
                  equation[j].kind = OPERATOR;
                  equation[j].token.operatr = PLUS;
                  j++;
                  equation[j].level = level;
                  equation[j].kind = CONSTANT;
                  equation[j].token.constant = 1.0;
                  j++;
                  for (k = i; k < j; k++)
                        equation[k].level++;
                  equation[j].level = level;
                  equation[j].kind = OPERATOR;
                  equation[j].token.operatr = TIMES;
                  j++;
                  equation[j].level = level;
                  equation[j].kind = CONSTANT;
                  equation[j].token.constant = -1.0;
                  j++;
                  blt(&equation[eloc+len+3], &equation[eloc], (*np - eloc) * sizeof(token_type));
                  *np += len + 3;
                  k = eloc;
                  equation[k].level = mlevel;
                  equation[k].kind = OPERATOR;
                  equation[k].token.operatr = TIMES;
                  k++;
                  blt(&equation[k], &equation[i], len * sizeof(token_type));
                  k += len;
                  equation[k].level = mlevel + 1;
                  equation[k].kind = OPERATOR;
                  equation[k].token.operatr = FACTORIAL;
                  k++;
                  equation[k].level = mlevel + 1;
                  equation[k].kind = CONSTANT;
                  equation[k].token.constant = 0.0;
                  k++;
                  return true;
            }
      }
      if (*np + 2 > n_tokens) {
            error_huge();
      }
      blt(&equation[eloc+2], &equation[eloc], (*np - eloc) * sizeof(token_type));
      *np += 2;
      k = eloc;
      equation[k].level = mlevel;
      equation[k].kind = OPERATOR;
      equation[k].token.operatr = DIVIDE;
      k++;
      equation[k].level = mlevel;
      equation[k].kind = VARIABLE;
      equation[k].token.variable = v;
      return true;
}

static int
inv_laplace_sub(equation, np, loc, eloc, v)
token_type  *equation;
int         *np;
int         loc, eloc;
long        v;
{
      int   i, j, k;
      int   len;
      int   level, mlevel;
      int   count;

      count = 0;
      for (i = loc; i < eloc; i += 2) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  count++;
            }
      }
      if (count > 1)
            return false;
      mlevel = min_level(&equation[loc], eloc - loc) + 1;
      for (j = loc; j < eloc; j++)
            equation[j].level += 2;
      for (i = loc; i < eloc; i += 2) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  i++;
                  if (i >= eloc || equation[i].token.operatr != POWER)
                        return false;
                  if ((i - 2) <= loc || equation[i-2].token.operatr != DIVIDE)
                        return false;
                  level = equation[i].level;
                  i++;
                  for (j = i; j < eloc && equation[j].level >= level; j++)
                        equation[j].level++;
                  len = j - i;
                  if (*np + len + 7 > n_tokens)
                        error_huge();
                  equation[i-3].token.operatr = TIMES;
                  blt(&equation[j+2], &equation[j], (*np - j) * sizeof(token_type));
                  *np += 2;
                  eloc += 2;
                  len += 2;
                  level++;
                  equation[j].level = level;
                  equation[j].kind = OPERATOR;
                  equation[j].token.operatr = MINUS;
                  j++;
                  equation[j].level = level;
                  equation[j].kind = CONSTANT;
                  equation[j].token.constant = 1.0;
                  j++;
                  blt(&equation[eloc+len+3], &equation[eloc], (*np - eloc) * sizeof(token_type));
                  *np += len + 3;
                  k = eloc;
                  equation[k].level = mlevel;
                  equation[k].kind = OPERATOR;
                  equation[k].token.operatr = DIVIDE;
                  k++;
                  blt(&equation[k], &equation[i], len * sizeof(token_type));
                  k += len;
                  equation[k].level = mlevel + 1;
                  equation[k].kind = OPERATOR;
                  equation[k].token.operatr = FACTORIAL;
                  k++;
                  equation[k].level = mlevel + 1;
                  equation[k].kind = CONSTANT;
                  equation[k].token.constant = 0.0;
                  k++;
                  return true;
            }
      }
      return false;
}

/*
 * The laplace command.
 */
int
laplace_cmd(cp)
char  *cp;
{
      int         i;
      long        v;
      int         inverse_flag;
      token_type  *source, *dest;
      int         n1, *nps, *np;

      if (notdefined(cur_equation)) {
            return false;
      }
      i = next_espace();
      if (n_rhs[cur_equation]) {
            source = rhs[cur_equation];
            nps = &n_rhs[cur_equation];
            dest = rhs[i];
            np = &n_rhs[i];
      } else {
            source = lhs[cur_equation];
            nps = &n_lhs[cur_equation];
            dest = lhs[i];
            np = &n_lhs[i];
      }
      inverse_flag = (strcmp_tospace(cp, "inverse") == 0);
      if (inverse_flag) {
            cp = skip_param(cp);
      }
      if (*cp) {
            cp = parse_var2(&v, cp);
            if (cp == NULL) {
                  return false;
            }
      } else {
            if (!prompt_var(&v)) {
                  return false;
            }
      }
      if (extra_garbage(cp)) {
            return false;
      }
      partial_flag = false;
      uf_simp(source, nps);
      partial_flag = true;
      factorv(source, nps, v);
      blt(dest, source, *nps * sizeof(token_type));
      n1 = *nps;
      if (inverse_flag) {
            if (!poly_in_v(dest, n1, v, true)
                || !int_dispatch(dest, &n1, v, inv_laplace_sub)) {
                  error(_("Inverse Laplace failed."));
                  return false;
            }
      } else {
            if (!int_dispatch(dest, &n1, v, laplace_sub)) {
                  error(_("Laplace failed."));
                  return false;
            }
      }
      simp_loop(dest, &n1);
      if (n_rhs[cur_equation]) {
            blt(lhs[i], lhs[cur_equation], n_lhs[cur_equation] * sizeof(token_type));
            n_lhs[i] = n_lhs[cur_equation];
      }
      *np = n1;
      cur_equation = i;
      return_result(cur_equation);
      return true;
}

#if   !BASICS
/*
 * Numerical integrate command.
 */
int
nintegrate_cmd(cp)
char  *cp;
{
      long        v;
      int         i, j, k;
      int         i1, j1;
      int         level;
      int         iterations;
      int         first_size;
      token_type  *ep;
      int         trap_flag;
      int         singularity;
      token_type  *source, *dest;
      int         n1, *nps, *np;

      iterations = 1000;      /* must be even */
      first_size = 0;
      if (notdefined(cur_equation)) {
            return false;
      }
      i = next_espace();
      if (n_rhs[cur_equation]) {
            source = rhs[cur_equation];
            nps = &n_rhs[cur_equation];
            dest = rhs[i];
            np = &n_rhs[i];
      } else {
            source = lhs[cur_equation];
            nps = &n_lhs[cur_equation];
            dest = lhs[i];
            np = &n_lhs[i];
      }
      trap_flag = (strncasecmp(cp, "trap", 4) == 0);
      if (trap_flag) {
            cp = skip_param(cp);
      }
      if (*cp == '\0') {
            if (!prompt_var(&v))
                  return false;
      } else {
            cp = parse_var2(&v, cp);
            if (cp == NULL) {
                  return false;
            }
            if (*cp) {
                  iterations = decstrtol(cp, &cp);
            }
            if (*cp || iterations <= 0 || (iterations % 2) != 0) {
                  error(_("Number of partitions must be a positive, even integer."));
                  return false;
            }
      }
      singularity = false;
      for (j = 1; j < *nps; j += 2) {
            if (source[j].token.operatr == DIVIDE) {
                  level = source[j].level;
                  for (k = j + 1; k < *nps && source[k].level >= level; k++) {
                        if (source[k].kind == VARIABLE
                            && source[k].token.variable == v) {
                              singularity = true;
                        }
                  }
            }
      }
      if (singularity) {
            printf(_("Warning: Singularity detected, result of numerical integration may be wrong.\n"));
      }
      my_strlcpy(prompt_str, _("Enter lower bound: "), sizeof(prompt_str));
      if (!get_expr(tlhs, &n_tlhs)) {
            return false;
      }
      subst_constants(tlhs, &n_tlhs);
      for (j = 0; j < n_tlhs; j += 2) {
            if (tlhs[j].kind == CONSTANT && !isfinite(tlhs[j].token.constant)) {
                  error(_("Error: Bound contains infinity."));
                  return false;
            }
      }
      my_strlcpy(prompt_str, _("Enter upper bound: "), sizeof(prompt_str));
      if (!get_expr(trhs, &n_trhs)) {
            return false;
      }
      subst_constants(trhs, &n_trhs);
      for (j = 0; j < n_trhs; j += 2) {
            if (trhs[j].kind == CONSTANT && !isfinite(trhs[j].token.constant)) {
                  error(_("Error: Bound contains infinity."));
                  return false;
            }
      }
      if ((n_tlhs + 1 + n_trhs + 2) > n_tokens) {
            error_huge();
      }
      if (trap_flag) {
            debug_string(0, _("Approximating the definite integral using the trapezoidal method..."));
      } else {
            debug_string(0, _("Approximating the definite integral using Simpson's rule..."));
      }
      subst_constants(source, nps);
      simp_loop(source, nps);
      for (j = 0; j < n_trhs; j++)
            trhs[j].level += 2;
      trhs[n_trhs].level = 2;
      trhs[n_trhs].kind = OPERATOR;
      trhs[n_trhs].token.operatr = MINUS;
      n_trhs++;
      j = n_trhs;
      blt(&trhs[n_trhs], tlhs, n_tlhs * sizeof(token_type));
      n_trhs += n_tlhs;
      for (; j < n_trhs; j++)
            trhs[j].level += 2;
      trhs[n_trhs].level = 1;
      trhs[n_trhs].kind = OPERATOR;
      trhs[n_trhs].token.operatr = DIVIDE;
      n_trhs++;
      trhs[n_trhs].level = 1;
      trhs[n_trhs].kind = CONSTANT;
      trhs[n_trhs].token.constant = iterations;
      n_trhs++;
      simp_loop(trhs, &n_trhs);
      dest[0].level = 1;
      dest[0].kind = CONSTANT;
      dest[0].token.constant = 0.0;
      n1 = 1;
      for (j = 0; j <= iterations; j++) {
            if ((n1 + 1 + *nps) > n_tokens)
                  error_huge();
            for (k = 0; k < n1; k++)
                  dest[k].level++;
            ep = &dest[n1];
            ep->level = 1;
            ep->kind = OPERATOR;
            ep->token.operatr = PLUS;
            n1++;
            i1 = n1;
            blt(&dest[i1], source, *nps * sizeof(token_type));
            n1 += *nps;
            for (k = i1; k < n1; k++)
                  dest[k].level += 2;
            for (k = i1; k < n1; k += 2) {
                  if (dest[k].kind == VARIABLE
                      && dest[k].token.variable == v) {
                        level = dest[k].level;
                        j1 = n_tlhs + 2 + n_trhs;
                        if ((n1 + j1) > n_tokens)
                              error_huge();
                        blt(&dest[k+1+j1], &dest[k+1], (n1 - (k + 1)) * sizeof(token_type));
                        n1 += j1;
                        j1 = k;
                        blt(&dest[k], tlhs, n_tlhs * sizeof(token_type));
                        k += n_tlhs;
                        for (; j1 < k; j1++)
                              dest[j1].level += level + 1;
                        ep = &dest[k];
                        ep->level = level + 1;
                        ep->kind = OPERATOR;
                        ep->token.operatr = PLUS;
                        ep++;
                        ep->level = level + 2;
                        ep->kind = CONSTANT;
                        ep->token.constant = j;
                        ep++;
                        ep->level = level + 2;
                        ep->kind = OPERATOR;
                        ep->token.operatr = TIMES;
                        k += 3;
                        j1 = k;
                        blt(&dest[k], trhs, n_trhs * sizeof(token_type));
                        k += n_trhs;
                        for (; j1 < k; j1++)
                              dest[j1].level += level + 2;
                        k--;
                  }
            }
            if (j > 0 && j < iterations) {
                  if ((n1 + 2) > n_tokens)
                        error_huge();
                  ep = &dest[n1];
                  ep->level = 2;
                  ep->kind = OPERATOR;
                  ep->token.operatr = TIMES;
                  ep++;
                  ep->level = 2;
                  ep->kind = CONSTANT;
                  if (trap_flag) {
                        ep->token.constant = 2.0;
                  } else {
                        if ((j & 1) == 1) {
                              ep->token.constant = 4.0;
                        } else {
                              ep->token.constant = 2.0;
                        }
                  }
                  n1 += 2;
            }
            in_calc_cmd = true;
            elim_loop(dest, &n1);
            ufactor(dest, &n1);
            simp_divide(dest, &n1);
            in_calc_cmd = false;
            for (k = 0; k < n1; k += 2) {
                  if (dest[k].kind == CONSTANT && !isfinite(dest[k].token.constant)) {
                        error(_("Integration failed because result contains infinity or nan."));
                        return false;
                  }
            }
            if (j > 0) {
                  if (j == 1) {
                        first_size = n1;
                        if (first_size < 4)
                              first_size = 4;
                  } else {
                        if ((n1 / 8) >= first_size) {
                              error(_("Integration failed."));
                              return false;
                        }
                  }
            }
      }
      if ((n1 + 3 + n_trhs) > n_tokens)
            error_huge();
      for (k = 0; k < n1; k++)
            dest[k].level++;
      ep = &dest[n1];
      ep->level = 1;
      ep->kind = OPERATOR;
      ep->token.operatr = DIVIDE;
      ep++;
      ep->level = 1;
      ep->kind = CONSTANT;
      if (trap_flag) {
            ep->token.constant = 2.0;
      } else {
            ep->token.constant = 3.0;
      }
      ep++;
      ep->level = 1;
      ep->kind = OPERATOR;
      ep->token.operatr = TIMES;
      n1 += 3;
      k = n1;
      blt(&dest[k], trhs, n_trhs * sizeof(token_type));
      n1 += n_trhs;
      for (; k < n1; k++)
            dest[k].level++;
      elim_loop(dest, &n1);
      ufactor(dest, &n1);
      simp_divide(dest, &n1);
      debug_string(0, _("Integration successful."));
      if (n_rhs[cur_equation]) {
            blt(lhs[i], lhs[cur_equation], n_lhs[cur_equation] * sizeof(token_type));
            n_lhs[i] = n_lhs[cur_equation];
      }
      *np = n1;
      cur_equation = i;
      return_result(cur_equation);
      return true;
}
#endif

Generated by  Doxygen 1.6.0   Back to index