1 /*
2 * Mathomatic symbolic differentiation routines and related commands.
3 *
4 * Copyright (C) 1987-2012 George Gesslein II.
5
6 This library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the License, or (at your option) any later version.
10
11 This library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
15
16 You should have received a copy of the GNU Lesser General Public
17 License along with this library; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19
20 The chief copyright holder can be contacted at gesslein@mathomatic.org, or
21 George Gesslein II, P.O. Box 224, Lansing, NY 14882-0224 USA.
22
23 */
24
25 #include "includes.h"
26
27 static int d_recurse(token_type *equation, int *np, int loc, int level, long v);
28
29 /*
30 * Compute the derivative of an equation side, with respect to variable "v",
31 * using the fast, rule-based transform method.
32 * This is done by recursively applying the proper rule of differentiation
33 * for each operator encountered.
34 *
35 * Return true if successful.
36 * The result must be simplified by the caller.
37 */
38 int
differentiate(equation,np,v)39 differentiate(equation, np, v)
40 token_type *equation; /* pointer to source and destination equation side */
41 int *np; /* pointer to the length of the equation side */
42 long v; /* differentiation variable */
43 {
44 int i;
45
46 organize(equation, np);
47 /* First put every times and divide on a level by itself. */
48 for (i = 1; i < *np; i += 2) {
49 switch (equation[i].token.operatr) {
50 case TIMES:
51 case DIVIDE:
52 binary_parenthesize(equation, *np, i);
53 }
54 }
55 return d_recurse(equation, np, 0, 1, v);
56 }
57
58 /*
59 * Recursive differentiation routine.
60 *
61 * Symbolically differentiate expression in "equation"
62 * (which is a standard equation side) starting at "loc".
63 * The current level of parentheses is "level" and
64 * do the differentiation with respect to variable "v".
65 *
66 * Return true if successful.
67 * Return false if it is beyond this program's capabilities or an error was encountered.
68 */
69 static int
d_recurse(equation,np,loc,level,v)70 d_recurse(equation, np, loc, level, v)
71 token_type *equation;
72 int *np, loc, level;
73 long v;
74 {
75 int i, j;
76 int n;
77 int op;
78 int oploc, endloc;
79 complexs c;
80
81 if (equation[loc].level < level) {
82 /* First differentiate if it is a single variable or constant. */
83 /* If it is the specified variable, change it to the constant 1, */
84 /* otherwise change it to the constant 0. */
85 if (equation[loc].kind == VARIABLE
86 && ((v == MATCH_ANY && (equation[loc].token.variable & VAR_MASK) > SIGN)
87 || equation[loc].token.variable == v)) {
88 equation[loc].kind = CONSTANT;
89 equation[loc].token.constant = 1.0;
90 } else {
91 equation[loc].kind = CONSTANT;
92 equation[loc].token.constant = 0.0;
93 }
94 return true;
95 }
96 for (op = 0, oploc = endloc = loc + 1; endloc < *np && equation[endloc].level >= level; endloc += 2) {
97 if (equation[endloc].level == level) {
98 switch (op) {
99 case 0:
100 case PLUS:
101 case MINUS:
102 break;
103 default:
104 /* Oops. More than one operator on the same level in this expression. */
105 error_bug("Internal error in d_recurse(): differentiating with unparenthesized operators is not allowed.");
106 return false;
107 }
108 op = equation[endloc].token.operatr;
109 oploc = endloc;
110 }
111 }
112 switch (op) {
113 case 0:
114 case PLUS:
115 case MINUS:
116 break;
117 case TIMES:
118 goto d_times;
119 case DIVIDE:
120 goto d_divide;
121 case POWER:
122 goto d_power;
123 default:
124 /* Differentiate an unsupported operator. */
125 /* This is possible if the expression doesn't contain the specified variable. */
126 /* In that case, the expression is replaced with "0", otherwise return false. */
127 for (i = loc; i < endloc; i += 2) {
128 if (equation[i].kind == VARIABLE
129 && ((v == MATCH_ANY && (equation[i].token.variable & VAR_MASK) > SIGN)
130 || equation[i].token.variable == v)) {
131 return false;
132 }
133 }
134 blt(&equation[loc+1], &equation[endloc], (*np - endloc) * sizeof(token_type));
135 *np -= (endloc - (loc + 1));
136 equation[loc].level = level;
137 equation[loc].kind = CONSTANT;
138 equation[loc].token.constant = 0.0;
139 return true;
140 }
141 /* Differentiate PLUS and MINUS operators. */
142 /* Use addition rule: d(u+v) = d(u) + d(v), */
143 /* where "d()" is the derivative function */
144 /* and "u" and "v" are expressions. */
145 for (i = loc; i < *np && equation[i].level >= level;) {
146 if (equation[i].kind != OPERATOR) {
147 if (!d_recurse(equation, np, i, level + 1, v))
148 return false;
149 i++;
150 for (; i < *np && equation[i].level > level; i += 2)
151 ;
152 continue;
153 }
154 i++;
155 }
156 return true;
157 d_times:
158 /* Differentiate TIMES operator. */
159 /* Use product rule: d(u*v) = u*d(v) + v*d(u). */
160 if (*np + 1 + (endloc - loc) > n_tokens) {
161 error_huge();
162 }
163 for (i = loc; i < endloc; i++)
164 equation[i].level++;
165 blt(&equation[endloc+1], &equation[loc], (*np - loc) * sizeof(token_type));
166 *np += 1 + (endloc - loc);
167 equation[endloc].level = level;
168 equation[endloc].kind = OPERATOR;
169 equation[endloc].token.operatr = PLUS;
170 if (!d_recurse(equation, np, endloc + (oploc - loc) + 2, level + 2, v))
171 return false;
172 return(d_recurse(equation, np, loc, level + 2, v));
173 d_divide:
174 /* Differentiate DIVIDE operator. */
175 /* Use quotient rule: d(u/v) = (v*d(u) - u*d(v))/v^2. */
176 if (*np + 3 + (endloc - loc) + (endloc - oploc) > n_tokens) {
177 error_huge();
178 }
179 for (i = loc; i < endloc; i++)
180 equation[i].level += 2;
181 equation[oploc].token.operatr = TIMES;
182 j = 1 + (endloc - loc);
183 blt(&equation[endloc+1], &equation[loc], (*np - loc) * sizeof(token_type));
184 *np += j;
185 equation[endloc].level = level + 1;
186 equation[endloc].kind = OPERATOR;
187 equation[endloc].token.operatr = MINUS;
188 j += endloc;
189 blt(&equation[j+2+(endloc-oploc)], &equation[j], (*np - j) * sizeof(token_type));
190 *np += 2 + (endloc - oploc);
191 equation[j].level = level;
192 equation[j].kind = OPERATOR;
193 equation[j].token.operatr = DIVIDE;
194 blt(&equation[j+1], &equation[oploc+1], (endloc - (oploc + 1)) * sizeof(token_type));
195 j += endloc - oploc;
196 equation[j].level = level + 1;
197 equation[j].kind = OPERATOR;
198 equation[j].token.operatr = POWER;
199 j++;
200 equation[j].level = level + 1;
201 equation[j].kind = CONSTANT;
202 equation[j].token.constant = 2.0;
203 if (!d_recurse(equation, np, endloc + (oploc - loc) + 2, level + 3, v))
204 return false;
205 return(d_recurse(equation, np, loc, level + 3, v));
206 d_power:
207 /* Differentiate POWER operator. */
208 /* Since we don't have symbolic logarithms, do all we can without them. */
209 for (i = oploc; i < endloc; i++) {
210 if (equation[i].kind == VARIABLE
211 && ((v == MATCH_ANY && (equation[i].token.variable & VAR_MASK) > SIGN)
212 || equation[i].token.variable == v)) {
213 if (parse_complex(&equation[loc], oploc - loc, &c)) {
214 c = complex_log(c);
215 n = (endloc - oploc) + 6;
216 if (*np + n > n_tokens) {
217 error_huge();
218 }
219 blt(&equation[endloc+n], &equation[endloc], (*np - endloc) * sizeof(token_type));
220 *np += n;
221 n = endloc;
222 equation[n].level = level;
223 equation[n].kind = OPERATOR;
224 equation[n].token.operatr = TIMES;
225 n++;
226 equation[n].level = level + 1;
227 equation[n].kind = CONSTANT;
228 equation[n].token.constant = c.re;
229 n++;
230 equation[n].level = level + 1;
231 equation[n].kind = OPERATOR;
232 equation[n].token.operatr = PLUS;
233 n++;
234 equation[n].level = level + 2;
235 equation[n].kind = CONSTANT;
236 equation[n].token.constant = c.im;
237 n++;
238 equation[n].level = level + 2;
239 equation[n].kind = OPERATOR;
240 equation[n].token.operatr = TIMES;
241 n++;
242 equation[n].level = level + 2;
243 equation[n].kind = VARIABLE;
244 equation[n].token.variable = IMAGINARY;
245 n++;
246 equation[n].level = level;
247 equation[n].kind = OPERATOR;
248 equation[n].token.operatr = TIMES;
249 n++;
250 blt(&equation[n], &equation[oploc+1], (endloc - (oploc + 1)) * sizeof(token_type));
251 for (i = loc; i < endloc; i++) {
252 equation[i].level++;
253 }
254 return(d_recurse(equation, np, n, level + 1, v));
255 }
256 return false;
257 }
258 }
259 blt(scratch, &equation[oploc+1], (endloc - (oploc + 1)) * sizeof(token_type));
260 n = endloc - (oploc + 1);
261 scratch[n].level = level;
262 scratch[n].kind = OPERATOR;
263 scratch[n].token.operatr = TIMES;
264 n++;
265 if (n + (endloc - loc) + 2 > n_tokens) {
266 error_huge();
267 }
268 blt(&scratch[n], &equation[loc], (endloc - loc) * sizeof(token_type));
269 i = n;
270 n += oploc + 1 - loc;
271 for (; i < n; i++)
272 scratch[i].level++;
273 n += endloc - (oploc + 1);
274 for (; i < n; i++)
275 scratch[i].level += 2;
276 scratch[n].level = level + 2;
277 scratch[n].kind = OPERATOR;
278 scratch[n].token.operatr = MINUS;
279 n++;
280 scratch[n].level = level + 2;
281 scratch[n].kind = CONSTANT;
282 scratch[n].token.constant = 1.0;
283 n++;
284 if (n + (oploc - loc) + 1 > n_tokens) {
285 error_huge();
286 }
287 scratch[n].level = level;
288 scratch[n].kind = OPERATOR;
289 scratch[n].token.operatr = TIMES;
290 n++;
291 j = n;
292 blt(&scratch[n], &equation[loc], (oploc - loc) * sizeof(token_type));
293 n += oploc - loc;
294 if (*np - (endloc - loc) + n > n_tokens) {
295 error_huge();
296 }
297 blt(&equation[loc+n], &equation[endloc], (*np - endloc) * sizeof(token_type));
298 *np += loc + n - endloc;
299 blt(&equation[loc], scratch, n * sizeof(token_type));
300 return(d_recurse(equation, np, loc + j, level + 1, v));
301 }
302
303 /*
304 * The derivative command.
305 */
306 int
derivative_cmd(cp)307 derivative_cmd(cp)
308 char *cp;
309 {
310 int i, len;
311 long v = 0; /* Mathomatic variable */
312 long l1, order = 1;
313 token_type *source, *dest;
314 int n1, *nps, *np;
315 int simplify_flag = true, solved;
316
317 if (current_not_defined()) {
318 return false;
319 }
320 solved = solved_equation(cur_equation);
321 if (strcmp_tospace(cp, "nosimplify") == 0) {
322 simplify_flag = false;
323 cp = skip_param(cp);
324 }
325 i = next_espace();
326 if (n_rhs[cur_equation]) {
327 if (!solved) {
328 warning(_("Not a solved equation. Only the RHS will be differentiated."));
329 }
330 source = rhs[cur_equation];
331 nps = &n_rhs[cur_equation];
332 dest = rhs[i];
333 np = &n_rhs[i];
334 } else {
335 source = lhs[cur_equation];
336 nps = &n_lhs[cur_equation];
337 dest = lhs[i];
338 np = &n_lhs[i];
339 }
340 /* parse the command line or prompt: */
341 if (*cp) {
342 if (is_all(cp)) {
343 cp = skip_param(cp);
344 v = MATCH_ANY;
345 } else {
346 if (isvarchar(*cp)) {
347 cp = parse_var2(&v, cp);
348 if (cp == NULL) {
349 return false;
350 }
351 }
352 }
353 if (*cp) {
354 order = decstrtol(cp, &cp);
355 }
356 if (order <= 0) {
357 error(_("The order must be a positive integer."));
358 return false;
359 }
360 if (extra_characters(cp))
361 return false;
362 }
363 if (no_vars(source, *nps, &v)) {
364 warning(_("Current expression contains no variables; the derivative will be zero."));
365 } else {
366 if (v == 0) {
367 if (!prompt_var(&v)) {
368 return false;
369 }
370 }
371 if (v && v != MATCH_ANY && !found_var(source, *nps, v)) {
372 warning(_("Specified variable not found; the derivative will be zero."));
373 }
374 }
375 if (v == 0) {
376 error(_("No differentiation variable specified."));
377 return false;
378 }
379 #if !SILENT
380 list_var(v, 0);
381 if (n_rhs[cur_equation]) {
382 fprintf(gfp, _("Differentiating the RHS with respect to %s"), var_str);
383 } else {
384 fprintf(gfp, _("Differentiating with respect to %s"), var_str);
385 }
386 if (order != 1) {
387 fprintf(gfp, _(" %ld times"), order);
388 }
389 if (simplify_flag) {
390 fprintf(gfp, _(" and simplifying"));
391 } else {
392 fprintf(gfp, _(" and not simplifying"));
393 }
394 fprintf(gfp, "...\n");
395 #endif
396 blt(dest, source, *nps * sizeof(token_type));
397 n1 = *nps;
398 /* do the actual differentiating and simplifying: */
399 for (l1 = 0; l1 < order; l1++) {
400 if (order != 1) {
401 if (n1 == 1 && dest[0].kind == CONSTANT && dest[0].token.constant == 0.0) {
402 fprintf(gfp, _("0 reached after %ld derivatives taken.\n"), l1);
403 order = l1;
404 break;
405 }
406 }
407 if (!differentiate(dest, &n1, v)) {
408 error(_("Differentiation failed."));
409 return false;
410 }
411 if (simplify_flag) {
412 simpa_repeat_side(dest, &n1, true, false);
413 } else {
414 elim_loop(dest, &n1);
415 }
416 }
417 *np = n1;
418 if (n_rhs[cur_equation]) {
419 blt(lhs[i], lhs[cur_equation], n_lhs[cur_equation] * sizeof(token_type));
420 n_lhs[i] = n_lhs[cur_equation];
421 if (solved && isvarchar('\'')) {
422 len = list_var(lhs[i][0].token.variable, 0);
423 for (l1 = 0; l1 < order && len < (MAX_VAR_LEN - 1); l1++) {
424 var_str[len++] = '\'';
425 }
426 var_str[len] = '\0';
427 if (l1 == order) {
428 parse_var(&lhs[i][0].token.variable, var_str);
429 }
430 }
431 }
432 cur_equation = i;
433 return return_result(cur_equation);
434 }
435
436 /*
437 * The extrema command.
438 */
439 int
extrema_cmd(cp)440 extrema_cmd(cp)
441 char *cp;
442 {
443 int i;
444 long v = 0; /* Mathomatic variable */
445 long l1, order = 1;
446 token_type want;
447 token_type *source;
448 int n;
449
450 if (current_not_defined()) {
451 return false;
452 }
453 i = next_espace();
454 if (n_rhs[cur_equation]) {
455 if (!solved_equation(cur_equation)) {
456 error(_("The current equation is not solved for a variable."));
457 return false;
458 }
459 source = rhs[cur_equation];
460 n = n_rhs[cur_equation];
461 } else {
462 source = lhs[cur_equation];
463 n = n_lhs[cur_equation];
464 }
465 if (*cp) {
466 if (isvarchar(*cp)) {
467 cp = parse_var2(&v, cp);
468 if (cp == NULL) {
469 return false;
470 }
471 }
472 if (*cp) {
473 order = decstrtol(cp, &cp);
474 }
475 if (order <= 0) {
476 error(_("The order must be a positive integer."));
477 return false;
478 }
479 if (extra_characters(cp))
480 return false;
481 }
482 show_usage = false;
483 if (no_vars(source, n, &v)) {
484 error(_("Current expression contains no variables."));
485 return false;
486 }
487 if (v == 0) {
488 if (!prompt_var(&v)) {
489 return false;
490 }
491 }
492 if (!found_var(source, n, v)) {
493 error(_("Specified variable not found; the derivative would be zero."));
494 return false;
495 }
496 blt(rhs[i], source, n * sizeof(token_type));
497 /* take derivatives with respect to the specified variable and simplify: */
498 for (l1 = 0; l1 < order; l1++) {
499 if (!differentiate(rhs[i], &n, v)) {
500 error(_("Differentiation failed."));
501 return false;
502 }
503 simpa_repeat_side(rhs[i], &n, true, false);
504 }
505 if (!found_var(rhs[i], n, v)) {
506 error(_("There are no solutions."));
507 return false;
508 }
509 n_rhs[i] = n;
510 /* set equal to zero: */
511 n_lhs[i] = 1;
512 lhs[i][0] = zero_token;
513 cur_equation = i;
514 /* lastly, solve for the specified variable and simplify: */
515 want.level = 1;
516 want.kind = VARIABLE;
517 want.token.variable = v;
518 if (solve_sub(&want, 1, lhs[i], &n_lhs[i], rhs[i], &n_rhs[i]) <= 0) {
519 error(_("Solve failed."));
520 return false;
521 }
522 simpa_repeat_side(rhs[i], &n_rhs[i], false, false);
523 return return_result(cur_equation);
524 }
525
526 /*
527 * The taylor command.
528 */
529 int
taylor_cmd(cp)530 taylor_cmd(cp)
531 char *cp;
532 {
533 long v = 0; /* Mathomatic variable */
534 int i, j, k, i1;
535 int level;
536 long l1, n, order = -1L;
537 double d;
538 char *cp_start, *cp1 = NULL, buf[MAX_CMD_LEN];
539 int our;
540 int our_nlhs, our_nrhs;
541 token_type *ep, *source, *dest;
542 int n1, *nps, *np;
543 int simplify_flag = true;
544
545 cp_start = cp;
546 if (current_not_defined()) {
547 return false;
548 }
549 if (strcmp_tospace(cp, "nosimplify") == 0) {
550 simplify_flag = false;
551 cp = skip_param(cp);
552 }
553 i = next_espace();
554 blt(lhs[i], lhs[cur_equation], n_lhs[cur_equation] * sizeof(token_type));
555 n_lhs[i] = n_lhs[cur_equation];
556 n_rhs[i] = 0;
557 our = alloc_next_espace();
558 n_lhs[i] = 0;
559 if (our < 0) {
560 error(_("Out of free equation spaces."));
561 show_usage = false;
562 return false;
563 }
564 if (n_rhs[cur_equation]) {
565 source = rhs[cur_equation];
566 nps = &n_rhs[cur_equation];
567 dest = rhs[i];
568 np = &n_rhs[i];
569 } else {
570 source = lhs[cur_equation];
571 nps = &n_lhs[cur_equation];
572 dest = lhs[i];
573 np = &n_lhs[i];
574 }
575 if (*cp && isvarchar(*cp)) {
576 cp = parse_var2(&v, cp);
577 if (cp == NULL) {
578 return false;
579 }
580 }
581 if (*cp) {
582 order = decstrtol(cp, &cp1);
583 if (cp1 != skip_param(cp) || order < 0) {
584 error(_("Positive integer required for order."));
585 return false;
586 }
587 cp = cp1;
588 }
589 show_usage = false;
590 no_vars(source, *nps, &v);
591 if (v == 0) {
592 if (!prompt_var(&v)) {
593 return false;
594 }
595 }
596 if (!found_var(source, *nps, v)) {
597 warning(_("Specified differentiation variable not found; the derivative will be 0."));
598 }
599 blt(rhs[our], source, *nps * sizeof(token_type));
600 our_nrhs = *nps;
601 /* Simplify and take the first derivative: */
602 uf_simp(rhs[our], &our_nrhs);
603 if (!differentiate(rhs[our], &our_nrhs, v)) {
604 error(_("Differentiation failed."));
605 return false;
606 }
607 if (*cp) {
608 input_column += (cp - cp_start);
609 cp = parse_expr(lhs[our], &our_nlhs, cp, true);
610 if (cp == NULL || extra_characters(cp) || our_nlhs <= 0) {
611 show_usage = true;
612 return false;
613 }
614 } else {
615 #if !SILENT
616 list_var(v, 0);
617 printf(_("Taylor series expansion around %s = point.\n"), var_str);
618 #endif
619 my_strlcpy(prompt_str, _("Enter point (an expression; usually 0): "), sizeof(prompt_str));
620 if (!get_expr(lhs[our], &our_nlhs)) {
621 return false;
622 }
623 }
624 if (order < 0) {
625 my_strlcpy(prompt_str, _("Enter order (number of derivatives to take): "), sizeof(prompt_str));
626 if ((cp1 = get_string(buf, sizeof(buf))) == NULL)
627 return false;
628 if (*cp1) {
629 cp = NULL;
630 order = decstrtol(cp1, &cp);
631 if (cp == NULL || *cp || order < 0) {
632 error(_("Positive integer required for order."));
633 return false;
634 }
635 } else {
636 order = LONG_MAX - 1;
637 #if !SILENT
638 printf(_("Derivatives will be taken until they reach zero...\n"));
639 #endif
640 }
641 }
642 #if !SILENT
643 fprintf(gfp, _("Taylor series"));
644 if (n_rhs[cur_equation]) {
645 fprintf(gfp, _(" of the RHS"));
646 }
647 list_var(v, 0);
648 fprintf(gfp, _(" with respect to %s"), var_str);
649 if (simplify_flag) {
650 fprintf(gfp, _(", simplified"));
651 } else {
652 fprintf(gfp, _(", not simplified"));
653 }
654 fprintf(gfp, "...\n");
655 #endif
656 n = 0;
657 i1 = 0;
658 blt(dest, source, *nps * sizeof(token_type));
659 n1 = *nps;
660 loop_again:
661 for (k = i1; k < n1; k += 2) {
662 if (dest[k].kind == VARIABLE && dest[k].token.variable == v) {
663 level = dest[k].level;
664 if ((n1 + our_nlhs - 1) > n_tokens)
665 error_huge();
666 blt(&dest[k+our_nlhs], &dest[k+1], (n1 - (k + 1)) * sizeof(token_type));
667 n1 += our_nlhs - 1;
668 j = k;
669 blt(&dest[k], lhs[our], our_nlhs * sizeof(token_type));
670 k += our_nlhs;
671 for (; j < k; j++)
672 dest[j].level += level;
673 k--;
674 }
675 }
676 if ((n1 + our_nlhs + 7) > n_tokens)
677 error_huge();
678 for (k = i1; k < n1; k++)
679 dest[k].level++;
680 ep = &dest[n1];
681 ep->level = 1;
682 ep->kind = OPERATOR;
683 ep->token.operatr = TIMES;
684 ep++;
685 ep->level = 3;
686 ep->kind = VARIABLE;
687 ep->token.variable = v;
688 ep++;
689 ep->level = 3;
690 ep->kind = OPERATOR;
691 ep->token.operatr = MINUS;
692 n1 += 3;
693 j = n1;
694 blt(&dest[n1], lhs[our], our_nlhs * sizeof(token_type));
695 n1 += our_nlhs;
696 for (; j < n1; j++)
697 dest[j].level += 3;
698 ep = &dest[n1];
699 ep->level = 2;
700 ep->kind = OPERATOR;
701 ep->token.operatr = POWER;
702 ep++;
703 ep->level = 2;
704 ep->kind = CONSTANT;
705 ep->token.constant = n;
706 ep++;
707 ep->level = 1;
708 ep->kind = OPERATOR;
709 ep->token.operatr = DIVIDE;
710 ep++;
711 for (d = 1.0, l1 = 2; l1 <= n; l1++)
712 d *= l1;
713 ep->level = 1;
714 ep->kind = CONSTANT;
715 ep->token.constant = d;
716 n1 += 4;
717 for (; i1 < n1; i1++)
718 dest[i1].level++;
719 if (simplify_flag) {
720 uf_simp(dest, &n1);
721 }
722 side_debug(1, dest, n1);
723 if (exp_contains_infinity(dest, n1)) {
724 error(_("Result invalid because it contains infinity or NaN."));
725 return false;
726 }
727 if (n < order) {
728 if (n > 0) {
729 if (!differentiate(rhs[our], &our_nrhs, v)) {
730 error(_("Differentiation failed."));
731 return false;
732 }
733 }
734 /* symb_flag = symblify; */
735 simpa_repeat_side(rhs[our], &our_nrhs, true, false /* was true */);
736 /* symb_flag = false; */
737 if (our_nrhs != 1 || rhs[our][0].kind != CONSTANT || rhs[our][0].token.constant != 0.0) {
738 i1 = n1;
739 if ((i1 + 1 + our_nrhs) > n_tokens)
740 error_huge();
741 for (j = 0; j < i1; j++)
742 dest[j].level++;
743 dest[i1].level = 1;
744 dest[i1].kind = OPERATOR;
745 dest[i1].token.operatr = PLUS;
746 i1++;
747 blt(&dest[i1], rhs[our], our_nrhs * sizeof(token_type));
748 n1 = i1 + our_nrhs;
749 n++;
750 goto loop_again;
751 }
752 }
753 #if !SILENT
754 fprintf(gfp, _("%ld non-zero derivative%s applied.\n"), n, (n == 1) ? "" : "s");
755 #endif
756 if (n_rhs[cur_equation]) {
757 n_lhs[i] = n_lhs[cur_equation];
758 }
759 *np = n1;
760 cur_equation = i;
761 return return_result(cur_equation);
762 }
763
764 /*
765 * The limit command.
766 */
767 int
limit_cmd(cp)768 limit_cmd(cp)
769 char *cp;
770 {
771 int i;
772 long v = 0; /* Mathomatic variable */
773 token_type solved_v, want;
774 char *cp_start;
775
776 cp_start = cp;
777 if (current_not_defined()) {
778 return false;
779 }
780 i = next_espace();
781 if (n_rhs[cur_equation] == 0) {
782 /* make expression into an equation: */
783 blt(rhs[cur_equation], lhs[cur_equation], n_lhs[cur_equation] * sizeof(token_type));
784 n_rhs[cur_equation] = n_lhs[cur_equation];
785 n_lhs[cur_equation] = 1;
786 lhs[cur_equation][0].level = 1;
787 lhs[cur_equation][0].kind = VARIABLE;
788 parse_var(&lhs[cur_equation][0].token.variable, "limit");
789 }
790 if (!solved_equation(cur_equation)) {
791 error(_("The current equation is not solved for a variable."));
792 return false;
793 }
794 solved_v = lhs[cur_equation][0];
795 /* parse the command line or prompt: */
796 if (*cp) {
797 cp = parse_var2(&v, cp);
798 if (cp == NULL) {
799 return false;
800 }
801 }
802 show_usage = false;
803 if (no_vars(rhs[cur_equation], n_rhs[cur_equation], &v)) {
804 warning(_("Current expression contains no variables; that is the answer."));
805 return return_result(cur_equation);
806 }
807 if (v == 0) {
808 if (!prompt_var(&v)) {
809 return false;
810 }
811 }
812 if (!found_var(rhs[cur_equation], n_rhs[cur_equation], v)) {
813 warning(_("Limit variable not found; answer is original expression."));
814 return return_result(cur_equation);
815 }
816 if (*cp == '=') {
817 cp = skip_space(cp + 1);
818 }
819 if (*cp) {
820 input_column += (cp - cp_start);
821 cp = parse_expr(tes, &n_tes, cp, true);
822 if (cp == NULL || extra_characters(cp) || n_tes <= 0) {
823 show_usage = true;
824 return false;
825 }
826 } else {
827 list_var(v, 0);
828 snprintf(prompt_str, sizeof(prompt_str), _("as %s goes to: "), var_str);
829 if (!get_expr(tes, &n_tes)) {
830 return false;
831 }
832 }
833 simp_loop(tes, &n_tes);
834 #if !SILENT
835 list_var(v, 0);
836 fprintf(gfp, _("Taking the limit as %s goes to "), var_str);
837 list_proc(tes, n_tes, false);
838 fprintf(gfp, "\n");
839 #endif
840 /* copy the current equation to a new equation space, then simplify and work on the copy: */
841 copy_espace(cur_equation, i);
842 simpa_side(rhs[i], &n_rhs[i], false, false);
843
844 /* see if the limit expression is positive infinity: */
845 if (n_tes == 1 && tes[0].kind == CONSTANT && tes[0].token.constant == INFINITY) {
846 /* To take the limit to positive infinity, */
847 /* replace infinity with zero and replace the limit variable with its reciprocal: */
848 n_tes = 1;
849 tes[0] = zero_token;
850 tlhs[0] = one_token;
851 tlhs[1].level = 1;
852 tlhs[1].kind = OPERATOR;
853 tlhs[1].token.operatr = DIVIDE;
854 tlhs[2].level = 1;
855 tlhs[2].kind = VARIABLE;
856 tlhs[2].token.variable = v;
857 n_tlhs = 3;
858 subst_var_with_exp(rhs[i], &n_rhs[i], tlhs, n_tlhs, v);
859 }
860
861 /* General limit taking, solve for the limit variable: */
862 debug_string(0, _("Solving..."));
863 want.level = 1;
864 want.kind = VARIABLE;
865 want.token.variable = v;
866 if (solve_sub(&want, 1, lhs[i], &n_lhs[i], rhs[i], &n_rhs[i]) <= 0) {
867 error(_("Can't take the limit because solve failed."));
868 return false;
869 }
870 /* replace the limit variable (LHS) with the limit expression: */
871 blt(lhs[i], tes, n_tes * sizeof(token_type));
872 n_lhs[i] = n_tes;
873 /* simplify the RHS: */
874 symb_flag = symblify;
875 simpa_side(rhs[i], &n_rhs[i], false, false);
876 symb_flag = false;
877 if (exp_contains_nan(rhs[i], n_rhs[i])) {
878 error(_("Unable to take limit; result contains NaN (Not a Number)."));
879 return false;
880 }
881 /* solve back for the original variable: */
882 if (solve_sub(&solved_v, 1, lhs[i], &n_lhs[i], rhs[i], &n_rhs[i]) <= 0) {
883 error(_("Can't take the limit because solve failed."));
884 return false;
885 }
886 /* simplify before returning the result: */
887 simpa_side(rhs[i], &n_rhs[i], false, false);
888 if (exp_contains_nan(rhs[i], n_rhs[i])) {
889 error(_("Unable to take limit; result contains NaN (Not a Number)."));
890 return false;
891 }
892 return return_result(i);
893 }
894