GRFreire/mathomatic
1
0
Fork 0
mirror of https://github.com/mfillpot/mathomatic.git synced 2026-01-08 04:29:39 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity
master
mathomatic/integrate.c

996 lines
24 KiB
C
Raw Permalink Blame History

/*
* Mathomatic integration routines and commands.
*
* Copyright (C) 1987-2012 George Gesslein II.
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
The chief copyright holder can be contacted at gesslein@mathomatic.org, or
George Gesslein II, P.O. Box 224, Lansing, NY 14882-0224 USA.
*/
#include "includes.h"
static int integrate_sub(token_type *equation, int *np, int loc, int eloc, long v);
static int laplace_sub(token_type *equation, int *np, int loc, int eloc, long v);
static int inv_laplace_sub(token_type *equation, int *np, int loc, int eloc, long v);
static int constant_var_number = 1; /* makes unique numbers for the constant of integration */
/*
* Make variable "v" always raised to a power,
* unless it is on the right side of a power operator.
*/
void
make_powers(equation, np, v)
token_type *equation; /* pointer to beginning of equation side */
int *np; /* pointer to length of equation side */
long v; /* Mathomatic 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 for polynomials.
* Handles the level 1 additive operators,
* sending each polynomial term to the specified integration function.
*
* Return true if successful.
*/
int
int_dispatch(equation, np, v, func)
token_type *equation; /* pointer to beginning of equation side to integrate */
int *np; /* pointer to length of equation side */
long v; /* integration variable */
int (*func)(token_type *equation, int *np, int loc, int eloc, long v); /* integration function to call for each term */
{
int i, j;
make_powers(equation, np, v);
for (j = 0, i = 1;; i += 2) {
if (i >= *np) {
return((*func)(equation, np, j, i, v));
}
if (equation[i].level == 1
&& (equation[i].token.operatr == PLUS || equation[i].token.operatr == MINUS)) {
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;
}
}
}
}
return true;
}
/*
* Do the actual integration of a polynomial term.
*
* Return true if successful.
*/
static int
integrate_sub(equation, np, loc, eloc, v)
token_type *equation; /* pointer to beginning of equation side */
int *np; /* pointer to length of equation side */
int loc; /* beginning location of term */
int eloc; /* end location of term */
long v; /* variable of integration */
{
int i, j, k;
int len;
int level, vlevel, mlevel;
int count;
int div_flag;
level = min_level(&equation[loc], eloc - loc);
/* determine if the term is a polynomial term in "v" */
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;
vlevel = equation[i].level;
if (vlevel == level || vlevel == (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 (vlevel == (level + 1)) {
if ((i + 1) < eloc && equation[i+1].level == vlevel
&& equation[i+1].token.operatr == POWER) {
continue;
}
} else {
continue;
}
}
return false;
}
}
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;
blt(&equation[eloc+len+1], &equation[eloc], (*np - eloc) * sizeof(token_type));
*np += len + 1;
equation[eloc].level = mlevel;
equation[eloc].kind = OPERATOR;
equation[eloc].token.operatr = DIVIDE;
blt(&equation[eloc+1], &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;
equation[eloc].level = mlevel;
equation[eloc].kind = OPERATOR;
equation[eloc].token.operatr = TIMES;
eloc++;
equation[eloc].level = mlevel;
equation[eloc].kind = VARIABLE;
equation[eloc].token.variable = v;
return true;
}
/*
* The integrate command.
*/
int
integrate_cmd(cp)
char *cp;
{
int i, j;
int len;
long v = 0;
token_type *source, *dest;
int n1, n2, *nps, *np;
int definite_flag = false, constant_flag = false, solved;
double integrate_order = 1.0;
char var_name_buf[MAX_VAR_LEN], *cp_start;
long l1;
cp_start = cp;
if (current_not_defined()) {
return false;
}
n_tlhs = 0;
n_trhs = 0;
solved = solved_equation(cur_equation);
i = next_espace();
for (;; cp = skip_param(cp)) {
if (strcmp_tospace(cp, "definite") == 0) {
definite_flag = true;
continue;
}
if (strcmp_tospace(cp, "constant") == 0) {
constant_flag = true;
continue;
}
break;
}
if (constant_flag && definite_flag) {
error(_("Conflicting options given."));
return false;
}
if (n_rhs[cur_equation]) {
if (!solved) {
warning(_("Not a solved equation."));
}
debug_string(0, _("Only the RHS will be transformed."));
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) {
if (isvarchar(*cp)) {
cp = parse_var2(&v, cp);
if (cp == NULL) {
return false;
}
}
if (*cp) {
integrate_order = strtod(cp, &cp);
}
if (!isfinite(integrate_order) || integrate_order <= 0 || fmod(integrate_order, 1.0) != 0.0) {
error(_("The order must be a positive integer."));
return false;
}
}
if (*cp) {
cp = skip_comma_space(cp);
input_column += (cp - cp_start);
cp = parse_expr(tlhs, &n_tlhs, cp, false);
if (cp == NULL || n_tlhs <= 0) {
return false;
}
}
if (*cp) {
cp_start = cp;
cp = skip_comma_space(cp);
input_column += (cp - cp_start);
cp = parse_expr(trhs, &n_trhs, cp, false);
if (cp == NULL || extra_characters(cp) || n_trhs <= 0) {
return false;
}
}
show_usage = false;
if (v == 0) {
if (!prompt_var(&v)) {
return false;
}
}
#if !SILENT
list_var(v, 0);
if (n_rhs[cur_equation]) {
fprintf(gfp, _("Integrating the RHS with respect to %s"), var_str);
} else {
fprintf(gfp, _("Integrating with respect to %s"), var_str);
}
if (integrate_order != 1.0) {
fprintf(gfp, _(" %.*g times"), precision, integrate_order);
}
fprintf(gfp, _(" and simplifying...\n"));
#endif
partial_flag = false;
uf_simp(source, nps);
partial_flag = true;
factorv(source, nps, v);
blt(dest, source, *nps * sizeof(token_type));
n1 = *nps;
for (l1 = 0; l1 < integrate_order; l1++) {
if (!int_dispatch(dest, &n1, v, integrate_sub)) {
error(_("Integration failed, not a polynomial."));
return false;
}
if (constant_flag) {
if (n1 + 2 > n_tokens) {
error_huge();
}
for (j = 0; j < n1; j++) {
dest[j].level++;
}
dest[n1].kind = OPERATOR;
dest[n1].level = 1;
dest[n1].token.operatr = PLUS;
n1++;
dest[n1].kind = VARIABLE;
dest[n1].level = 1;
snprintf(var_name_buf, sizeof(var_name_buf), "C_%d", constant_var_number);
if (parse_var(&dest[n1].token.variable, var_name_buf) == NULL) {
return false;
}
n1++;
constant_var_number++;
if (constant_var_number < 0) {
constant_var_number = 1;
}
}
simp_loop(dest, &n1);
}
if (definite_flag) {
if (n_tlhs == 0) {
my_strlcpy(prompt_str, _("Enter lower bound: "), sizeof(prompt_str));
if (!get_expr(tlhs, &n_tlhs)) {
return false;
}
}
if (n_trhs == 0) {
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);
*np = 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];
if (solved && isvarchar('\'')) {
len = list_var(lhs[i][0].token.variable, 0);
for (l1 = 0; l1 < integrate_order && len > 0 && var_str[len-1] == '\''; l1++) {
var_str[--len] = '\0';
}
parse_var(&lhs[i][0].token.variable, var_str);
}
}
cur_equation = i;
return return_result(cur_equation);
}
/*
* Do the actual Laplace transformation of a polynomial term.
*
* Return true if successful.
*/
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, mlevel;
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;
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;
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 = 1.0;
return true;
}
}
if (*np + 2 > n_tokens) {
error_huge();
}
blt(&equation[eloc+2], &equation[eloc], (*np - eloc) * sizeof(token_type));
*np += 2;
equation[eloc].level = mlevel;
equation[eloc].kind = OPERATOR;
equation[eloc].token.operatr = DIVIDE;
eloc++;
equation[eloc].level = mlevel;
equation[eloc].kind = VARIABLE;
equation[eloc].token.variable = v;
return true;
}
/*
* Do the actual inverse Laplace transformation of a polynomial term.
*
* Return true if successful.
*/
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;
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;
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 = 1.0;
return true;
}
}
return false;
}
/*
* The laplace command.
*/
int
laplace_cmd(cp)
char *cp;
{
int i;
long v = 0;
int inverse_flag, solved;
token_type *source, *dest;
int n1, *nps, *np;
if (current_not_defined()) {
return false;
}
solved = solved_equation(cur_equation);
i = next_espace();
if (n_rhs[cur_equation]) {
if (!solved) {
warning(_("Not a solved equation."));
}
#if !SILENT
fprintf(gfp, _("Only the RHS will be transformed.\n"));
#endif
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;
}
if (extra_characters(cp)) {
return false;
}
}
show_usage = false;
if (v == 0) {
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;
if (inverse_flag) {
if (!poly_in_v(dest, n1, v, true) || !int_dispatch(dest, &n1, v, inv_laplace_sub)) {
error(_("Inverse Laplace transformation failed."));
return false;
}
} else {
if (!poly_in_v(dest, n1, v, false) || !int_dispatch(dest, &n1, v, laplace_sub)) {
error(_("Laplace transformation failed, not a polynomial."));
return false;
}
}
#if 1
simp_loop(dest, &n1);
#else
simpa_side(dest, &n1, false, false);
#endif
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 return_result(cur_equation);
}
/*
* Numerical integrate command.
*/
int
nintegrate_cmd(cp)
char *cp;
{
long v = 0; /* Mathomatic variable */
int i, j, k, i1, i2;
int len;
int level;
int iterations = 1000; /* must be even */
int first_size = 0;
int trap_flag, singularity, solved;
token_type *ep, *source, *dest;
int n1, *nps, *np;
char *cp_start;
cp_start = cp;
if (current_not_defined()) {
return false;
}
n_tlhs = 0;
n_trhs = 0;
solved = solved_equation(cur_equation);
i = next_espace();
if (n_rhs[cur_equation]) {
if (!solved) {
warning(_("Not a solved 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) {
cp = parse_var2(&v, cp);
if (cp == NULL) {
return false;
}
if (*cp) {
iterations = decstrtol(cp, &cp);
}
if (iterations <= 0 || (iterations % 2) != 0) {
error(_("Number of partitions must be a positive, even integer."));
return false;
}
}
if (*cp) {
input_column += (cp - cp_start);
cp = parse_expr(tlhs, &n_tlhs, cp, false);
if (cp == NULL || n_tlhs <= 0) {
return false;
}
}
if (*cp) {
cp_start = cp;
cp = skip_comma_space(cp);
input_column += (cp - cp_start);
cp = parse_expr(trhs, &n_trhs, cp, false);
if (cp == NULL || extra_characters(cp) || n_trhs <= 0) {
return false;
}
}
show_usage = false;
if (v == 0) {
if (!prompt_var(&v)) {
return false;
}
}
#if !SILENT
list_var(v, 0);
if (n_rhs[cur_equation]) {
fprintf(gfp, _("Numerically integrating the RHS with respect to %s...\n"), var_str);
} else {
fprintf(gfp, _("Numerically integrating with respect to %s...\n"), var_str);
}
#endif
singularity = false;
for (j = 1; j < *nps; j += 2) {
if (source[j].token.operatr == DIVIDE) {
for (k = j + 1; k < *nps && source[k].level >= source[j].level; k++) {
if (source[k].kind == VARIABLE && source[k].token.variable == v) {
singularity = true;
}
if (source[k].kind == OPERATOR && source[k].level == source[j].level)
break;
}
}
}
if (singularity) {
warning(_("Singularity detected, result of numerical integration might be wrong."));
}
if (n_tlhs == 0) {
my_strlcpy(prompt_str, _("Enter lower bound: "), sizeof(prompt_str));
if (!get_expr(tlhs, &n_tlhs)) {
return false;
}
}
subst_constants(tlhs, &n_tlhs);
simp_loop(tlhs, &n_tlhs);
if (exp_contains_infinity(tlhs, n_tlhs)) {
error(_("Not computable because: Lower bound contains infinity or NaN."));
return false;
}
if (n_trhs == 0) {
my_strlcpy(prompt_str, _("Enter upper bound: "), sizeof(prompt_str));
if (!get_expr(trhs, &n_trhs)) {
return false;
}
}
subst_constants(trhs, &n_trhs);
simp_loop(trhs, &n_trhs);
if (exp_contains_infinity(trhs, n_trhs)) {
error(_("Not computable because: Upper bound contains infinity or NaN."));
return false;
}
if ((n_tlhs + n_trhs + 3) > n_tokens) {
error_huge();
}
#if !SILENT
fprintf(gfp, _("Approximating the definite integral\n"));
if (trap_flag) {
fprintf(gfp, _("using the trapezoid method (%d partitions)...\n"), iterations);
} else {
fprintf(gfp, _("using Simpson's rule (%d partitions)...\n"), iterations);
}
#endif
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] = zero_token;
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;
i2 = n_tlhs + 2 + n_trhs;
if ((n1 + i2) > n_tokens)
error_huge();
blt(&dest[k+1+i2], &dest[k+1], (n1 - (k + 1)) * sizeof(token_type));
n1 += i2;
i2 = k;
blt(&dest[k], tlhs, n_tlhs * sizeof(token_type));
k += n_tlhs;
level++;
for (; i2 < k; i2++) {
dest[i2].level += level;
}
ep = &dest[k];
ep->level = level;
ep->kind = OPERATOR;
ep->token.operatr = PLUS;
ep++;
level++;
ep->level = level;
ep->kind = CONSTANT;
ep->token.constant = j;
ep++;
ep->level = level;
ep->kind = OPERATOR;
ep->token.operatr = TIMES;
k += 3;
i2 = k;
blt(&dest[k], trhs, n_trhs * sizeof(token_type));
k += n_trhs;
for (; i2 < k; i2++) {
dest[i2].level += level;
}
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;
}
/* simplify and approximate the partial result quickly: */
approximate_roots = true;
elim_loop(dest, &n1);
ufactor(dest, &n1);
simp_divide(dest, &n1);
factor_imaginary(dest, &n1);
approximate_roots = false;
side_debug(1, dest, n1);
if (exp_contains_infinity(dest, n1)) {
error(_("Integration failed because result contains infinity or NaN (a singularity)."));
return false;
}
/* detect an ever growing result: */
switch (j) {
case 0:
break;
case 1:
first_size = n1;
if (first_size < 4)
first_size = 4;
break;
default:
if ((n1 / 8) >= first_size) {
error(_("Result growing, integration failed."));
return false;
}
break;
}
}
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++;
/* simplify and approximate the result even more: */
approximate_roots = true;
do {
elim_loop(dest, &n1);
ufactor(dest, &n1);
simp_divide(dest, &n1);
} while (factor_imaginary(dest, &n1));
approximate_roots = false;
#if !SILENT
fprintf(gfp, _("Numerical integration successful:\n"));
#endif
*np = 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];
if (solved && isvarchar('\'')) {
len = list_var(lhs[i][0].token.variable, 0);
if (len > 0 && var_str[len-1] == '\'') {
var_str[--len] = '\0';
}
parse_var(&lhs[i][0].token.variable, var_str);
}
}
return return_result(i);
}
Reference in New Issue View Git Blame Copy Permalink