1 /*
2 * Mathomatic simplifying routines.
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 /*
28 * MAX_COMPARE_TERMS defines the maximum number of terms
29 * on the same level of parentheses that can be compared.
30 * Most of the stack usage is related this value.
31 * Space for this many pointers is reserved on the stack
32 * with each level of parentheses while comparing two expressions.
33 */
34 #define MAX_COMPARE_TERMS (DEFAULT_N_TOKENS / 6)
35
36 static int org_recurse(token_type *equation, int *np, int loc, int level, int *elocp);
37 static int const_recurse(token_type *equation, int *np, int loc, int level, int iflag);
38 static int compare_recurse(token_type *p1, int n1, int l1, token_type *p2, int n2, int l2, int *diff_signp);
39 static int order_recurse(token_type *equation, int *np, int loc, int level);
40
41 /*
42 * Fix up levels of parentheses in an equation side.
43 * Remove unnecessary parentheses and put similar operators on the same level.
44 *
45 * This is the inner-most loop in Mathomatic, make it fast.
46 */
47 void
organize(equation,np)48 organize(equation, np)
49 token_type *equation; /* equation side pointer */
50 int *np; /* pointer to length of equation side */
51 {
52 #if DEBUG
53 if (equation == NULL || np == NULL) {
54 error_bug("NULL pointer passed to organize().");
55 }
56 #endif
57 if (*np <= 0 || (*np & 1) == 0) {
58 printf("Bad expression size = %d.\n", *np);
59 error_bug("Internal error: organize() called with bad expression size.");
60 }
61 if (*np > n_tokens) {
62 error_bug("Internal error: expression array overflow detected in organize().");
63 }
64 org_recurse(equation, np, 0, 1, NULL);
65 }
66
67 static inline void
org_up_level(bp,ep,level,invert)68 org_up_level(bp, ep, level, invert)
69 token_type *bp, *ep;
70 int level, invert;
71 {
72 if (invert) {
73 for (; bp <= ep; bp++) {
74 bp->level--;
75 if (bp->level == level && bp->kind == OPERATOR) {
76 switch (bp->token.operatr) {
77 case PLUS:
78 bp->token.operatr = MINUS;
79 break;
80 case MINUS:
81 bp->token.operatr = PLUS;
82 break;
83 case TIMES:
84 bp->token.operatr = DIVIDE;
85 break;
86 case DIVIDE:
87 bp->token.operatr = TIMES;
88 break;
89 }
90 }
91 }
92 } else {
93 for (; bp <= ep; bp++) {
94 bp->level--;
95 }
96 }
97 }
98
99 /*
100 * Recurse through every sub-expression in "equation", starting at "loc",
101 * moving up levels to "level" of parentheses.
102 */
103 static int
org_recurse(equation,np,loc,level,elocp)104 org_recurse(equation, np, loc, level, elocp)
105 token_type *equation; /* equation side pointer */
106 int *np; /* pointer to length of equation side */
107 int loc, level, *elocp;
108 {
109 token_type *p1, *bp, *ep;
110 int op, sub_op;
111 int i;
112 int eloc;
113 int sub_eloc;
114 int min1;
115 int invert;
116
117 bp = &equation[loc];
118 ep = &equation[*np];
119 min1 = bp->level;
120 for (p1 = bp + 1; p1 < ep; p1 += 2) {
121 if (p1->level < min1) {
122 if (p1->level < level)
123 break;
124 min1 = p1->level;
125 }
126 }
127 ep = p1;
128 eloc = (ep - equation) - 1;
129 if (elocp)
130 *elocp = eloc;
131 if (eloc == loc) {
132 bp->level = max(level - 1, 1);
133 return 0;
134 }
135 if (min1 > level) {
136 for (p1 = bp; p1 < ep; p1++)
137 p1->level -= min1 - level;
138 }
139 op = 0;
140 for (p1 = bp + 1; p1 < ep; p1 += 2) {
141 if (p1->level == level) {
142 op = p1->token.operatr;
143 break;
144 }
145 }
146 for (i = loc; i <= eloc; i += 2) {
147 if (equation[i].level > level) {
148 sub_op = org_recurse(equation, np, i, level + 1, &sub_eloc);
149 switch (sub_op) {
150 case PLUS:
151 case MINUS:
152 if (op != PLUS && op != MINUS)
153 break;
154 invert = (i - 1 >= loc && equation[i-1].token.operatr == MINUS);
155 org_up_level(&equation[i], &equation[sub_eloc], level, invert);
156 break;
157 case TIMES:
158 case DIVIDE:
159 if (op != TIMES && op != DIVIDE)
160 break;
161 invert = (i - 1 >= loc && equation[i-1].token.operatr == DIVIDE);
162 org_up_level(&equation[i], &equation[sub_eloc], level, invert);
163 break;
164 }
165 i = sub_eloc;
166 }
167 }
168 return op;
169 }
170
171 /*
172 * The quickest, most basic simplification loop.
173 * Just does constant simplification.
174 */
175 void
elim_loop(equation,np)176 elim_loop(equation, np)
177 token_type *equation; /* pointer to the beginning of equation side to simplify */
178 int *np; /* pointer to length of equation side */
179 {
180 if (abort_flag) {
181 /* Control-C pressed, gracefully return to main prompt and leave unsimplified */
182 abort_flag = false;
183 #if DEBUG && !SILENT
184 char *cp, buf[100];
185
186 my_strlcpy(prompt_str, _("Enter debug level, or an empty line to abort the current operation: "), sizeof(prompt_str));
187 if ((cp = get_string(buf, sizeof(buf))) == NULL || *cp == '\0') {
188 longjmp(jmp_save, 13);
189 } else {
190 debug_level = decstrtol(cp, NULL);
191 printf(_("Debug level set to %d.\n"), debug_level);
192 }
193 #else
194 longjmp(jmp_save, 13);
195 #endif
196 }
197 side_debug(6, equation, *np);
198 do {
199 do {
200 do {
201 organize(equation, np);
202 } while (combine_constants(equation, np, true));
203 } while (elim_k(equation, np));
204 } while (simp_pp(equation, np));
205 if (reorder(equation, np)) {
206 do {
207 organize(equation, np);
208 } while (elim_k(equation, np));
209 }
210 side_debug(5, equation, *np);
211 }
212
213 /*
214 * Configurable high level simplify routine.
215 */
216 void
simp_ssub(equation,np,v,d,power_flag,times_flag,fc_level)217 simp_ssub(equation, np, v, d, power_flag, times_flag, fc_level)
218 token_type *equation; /* pointer to the beginning of equation side to simplify */
219 int *np; /* pointer to length of equation side */
220 long v; /* variable to factor, 0L or MATCH_ANY to factor all variables */
221 double d; /* factor expressions raised to the power of this if v */
222 int power_flag; /* factor_power() flag */
223 int times_flag; /* factor_times() flag */
224 int fc_level; /* factor constants code, passed to factor_constants() */
225 {
226 do {
227 do {
228 do {
229 do {
230 do {
231 do {
232 do {
233 do {
234 elim_loop(equation, np);
235 } while (simp2_power(equation, np));
236 } while (times_flag && factor_times(equation, np));
237 } while (elim_sign(equation, np));
238 } while (subtract_itself(equation, np));
239 } while (factor_constants(equation, np, fc_level));
240 } while (factor_divide(equation, np, v, d));
241 } while (factor_plus(equation, np, v, d));
242 } while (power_flag && factor_power(equation, np));
243 }
244
245 /*
246 * Quickly and very basically simplify an equation space.
247 * No factoring is done.
248 */
249 void
simp_equation(n)250 simp_equation(n)
251 int n; /* equation space number to simplify */
252 {
253 if (empty_equation_space(n))
254 return;
255 simp_loop(lhs[n], &n_lhs[n]);
256 if (n_rhs[n] > 0) {
257 simp_loop(rhs[n], &n_rhs[n]);
258 }
259 }
260
261 /*
262 * For quick, mid-range simplification of an equation side.
263 * Trivial factoring is done.
264 */
265 void
mid_simp_side(equation,np)266 mid_simp_side(equation, np)
267 token_type *equation; /* pointer to the beginning of equation side to simplify */
268 int *np; /* pointer to length of equation side */
269 {
270 simp_ssub(equation, np, 0L, 1.0, true, true, 6);
271 }
272
273 /*
274 * For quick, mid-range simplification of an equation space.
275 * Trivial factoring is done.
276 */
277 void
mid_simp_equation(n)278 mid_simp_equation(n)
279 int n; /* equation space number to simplify */
280 {
281 if (empty_equation_space(n))
282 return;
283 mid_simp_side(lhs[n], &n_lhs[n]);
284 if (n_rhs[n] > 0) {
285 mid_simp_side(rhs[n], &n_rhs[n]);
286 }
287 }
288
289 /*
290 * This function is the mid-range simplifier used by the solver.
291 */
292 void
simps_side(equation,np,zsolve)293 simps_side(equation, np, zsolve)
294 token_type *equation; /* pointer to the beginning of equation side to simplify */
295 int *np; /* pointer to length of equation side */
296 int zsolve; /* true for solving for zero */
297 {
298 elim_loop(equation, np);
299 simp_constant_power(equation, np);
300 do {
301 simp_ssub(equation, np, 0L, 0.0, !zsolve, true, 6);
302 } while (super_factor(equation, np, 0));
303 }
304
305 /*
306 * This function is used by the factor command.
307 */
308 void
simpv_side(equation,np,v)309 simpv_side(equation, np, v)
310 token_type *equation; /* pointer to the beginning of equation side to factor */
311 int *np; /* pointer to length of equation side */
312 long v; /* variable to factor, 0 for all variables */
313 {
314 simp_ssub(equation, np, v, 0.0, v == 0, true, 6);
315 }
316
317 /*
318 * For factoring an equation space like the factor command.
319 * Trivial factoring is done.
320 */
321 void
simpv_equation(n,v)322 simpv_equation(n, v)
323 int n; /* equation space number to simplify */
324 long v; /* Mathomatic variable to factor or 0 */
325 {
326 if (empty_equation_space(n))
327 return;
328 simpv_side(lhs[n], &n_lhs[n], v);
329 if (n_rhs[n] > 0) {
330 simpv_side(rhs[n], &n_rhs[n], v);
331 }
332 }
333
334 /*
335 * Factor out and simplify imaginary constants in an equation side.
336 * Do complex number root approximation and other complex number simplification.
337 *
338 * Return true if anything was approximated.
339 */
340 int
factor_imaginary(equation,np)341 factor_imaginary(equation, np)
342 token_type *equation; /* pointer to the beginning of equation side */
343 int *np; /* pointer to equation side length */
344 {
345 int rv;
346
347 rv = approximate_complex_roots(equation, np);
348 factorv(equation, np, IMAGINARY);
349 return rv;
350 }
351
352 /*
353 * Factor out only the specified variable "v" and simplify a little.
354 * Does not do any approximation.
355 *
356 * If v is IMAGINARY, simplify imaginary number division, so that the imaginary unit is
357 * not in the divisor.
358 */
359 void
factorv(equation,np,v)360 factorv(equation, np, v)
361 token_type *equation; /* pointer to the beginning of equation side to factor */
362 int *np; /* pointer to equation side length */
363 long v; /* Mathomatic variable to factor */
364 {
365 do {
366 do {
367 simp_loop(equation, np);
368 } while (factor_plus(equation, np, v, 0.0));
369 } while (v == IMAGINARY && div_imaginary(equation, np));
370 }
371
372 /*
373 * Simplify and approximate for the calculate command and for more successful comparisons.
374 * Includes complex number simplification.
375 */
376 void
calc_simp(equation,np)377 calc_simp(equation, np)
378 token_type *equation;
379 int *np;
380 {
381 approximate_roots = true;
382 subst_constants(equation, np);
383 mid_simp_side(equation, np);
384 factor_imaginary(equation, np);
385 ufactor(equation, np);
386 factor_imaginary(equation, np);
387 uf_simp(equation, np);
388 factor_imaginary(equation, np);
389 mid_simp_side(equation, np);
390 make_simple_fractions(equation, np);
391 uf_tsimp(equation, np);
392 approximate_roots = false;
393 }
394
395 /*
396 * Approximate an equation side for the approximate command.
397 */
398 void
approximate(equation,np)399 approximate(equation, np)
400 token_type *equation;
401 int *np;
402 {
403 if (repeat_flag) {
404 /* do more */
405 calc_simp(equation, np);
406 } else {
407 subst_constants(equation, np);
408 approximate_roots = true;
409 simp_loop(equation, np);
410 factor_imaginary(equation, np);
411 approximate_roots = false;
412 }
413 }
414
415 /*
416 * Try to eliminate imaginary units (i#) from an equation side by converting
417 * expressions like (i#*(b^.5)) to ((-1*b)^.5).
418 *
419 * Return true if any imaginary units where substituted.
420 */
421 int
simp_i(equation,np)422 simp_i(equation, np)
423 token_type *equation;
424 int *np;
425 {
426 int i;
427 int level;
428 int rv = false;
429
430 simp_loop(equation, np);
431 for (i = 0; i < *np; i++) {
432 switch (equation[i].kind) {
433 case VARIABLE:
434 if (equation[i].token.variable == IMAGINARY) {
435 if (*np + 2 > n_tokens) {
436 error_huge();
437 }
438 level = equation[i].level + 1;
439 blt(&equation[i+2], &equation[i], (*np - i) * sizeof(token_type));
440 *np += 2;
441 equation[i].level = level;
442 equation[i].kind = CONSTANT;
443 equation[i].token.constant = -1.0;
444 i++;
445 equation[i].level = level;
446 equation[i].kind = OPERATOR;
447 equation[i].token.operatr = POWER;
448 i++;
449 equation[i].level = level;
450 equation[i].kind = CONSTANT;
451 equation[i].token.constant = 0.5;
452 rv = true;
453 }
454 break;
455 case OPERATOR:
456 #if 0
457 if (equation[i].token.operatr != POWER)
458 continue;
459 level = equation[i].level;
460 /* Skip imaginary units on the right side of any power operator. */
461 for (; i < *np && equation[i].level >= level; i += 2)
462 ;
463 i--;
464 #endif
465 break;
466 case CONSTANT:
467 break;
468 }
469 }
470 do {
471 do {
472 do {
473 do {
474 do {
475 organize(equation, np);
476 } while (combine_constants(equation, np, false));
477 } while (elim_k(equation, np));
478 } while (simp_pp(equation, np));
479 /* The following factor_power() messes up absolute values and complex numbers sometimes: */
480 } while (factor_power(equation, np));
481 } while (factor_times(equation, np));
482 simp_loop(equation, np);
483 return rv;
484 }
485
486 /*
487 * Combine all like denominators.
488 */
489 void
simp_divide(equation,np)490 simp_divide(equation, np)
491 token_type *equation;
492 int *np;
493 {
494 do {
495 do {
496 simp_loop(equation, np);
497 } while (factor_constants(equation, np, 1));
498 } while (factor_divide(equation, np, 0L, 0.0));
499 }
500
501 /*
502 * Combine all like denominators containing variable "v".
503 * Don't call factor_times().
504 *
505 * For beauty simplifier simpb_side() below.
506 */
507 void
simp2_divide(equation,np,v,fc_level)508 simp2_divide(equation, np, v, fc_level)
509 token_type *equation;
510 int *np;
511 long v;
512 int fc_level;
513 {
514 do {
515 do {
516 do {
517 do {
518 do {
519 elim_loop(equation, np);
520 } while (simp2_power(equation, np));
521 } while (elim_sign(equation, np));
522 } while (subtract_itself(equation, np));
523 } while (factor_constants(equation, np, fc_level));
524 } while (factor_divide(equation, np, v, 0.0));
525 }
526
527 /*
528 * Compare function for qsort(3) within simpb_side() below.
529 */
530 static int
simpb_vcmp(p1,p2)531 simpb_vcmp(p1, p2)
532 sort_type *p1, *p2;
533 {
534 if (((p1->v & VAR_MASK) == SIGN) == ((p2->v & VAR_MASK) == SIGN)) {
535 if (p2->count == p1->count) {
536 if (p1->v < p2->v)
537 return -1;
538 if (p1->v == p2->v)
539 return 0;
540 return 1;
541 }
542 return(p2->count - p1->count);
543 } else {
544 if ((p1->v & VAR_MASK) == SIGN) {
545 return -1;
546 } else {
547 return 1;
548 }
549 }
550 }
551
552 /*
553 * Beauty simplifier for equation sides.
554 * This is the neat simplify and variable factor routine, to make a pleasing display.
555 * Factors variables in order: "sign" variables first, then by frequency.
556 */
557 void
simpb_side(equation,np,uf_power_flag,power_flag,fc_level)558 simpb_side(equation, np, uf_power_flag, power_flag, fc_level)
559 token_type *equation; /* pointer to the beginning of equation side */
560 int *np; /* pointer to length of equation side */
561 int uf_power_flag; /* uf_allpower() flag */
562 int power_flag; /* factor_power() flag */
563 int fc_level; /* factor constants code, passed to factor_constants() */
564 {
565 int i;
566 int vc, cnt; /* counts */
567 long v1, last_v; /* Mathomatic variables */
568 sort_type va[MAX_VARS]; /* array of all real variables found in equation side */
569
570 elim_loop(equation, np);
571 if (uf_power_flag) {
572 uf_allpower(equation, np);
573 }
574 last_v = 0;
575 for (vc = 0; vc < ARR_CNT(va);) {
576 cnt = 0;
577 v1 = -1;
578 for (i = 0; i < *np; i += 2) {
579 if (equation[i].kind == VARIABLE && equation[i].token.variable > last_v) {
580 if (v1 == -1 || equation[i].token.variable < v1) {
581 v1 = equation[i].token.variable;
582 cnt = 1;
583 } else if (equation[i].token.variable == v1) {
584 cnt++;
585 }
586 }
587 }
588 if (v1 == -1)
589 break;
590 last_v = v1;
591 if (v1 > IMAGINARY) { /* ignore constant variables */
592 va[vc].v = v1;
593 va[vc].count = cnt;
594 vc++;
595 }
596 }
597 if (vc) {
598 /* sort variable array va[] */
599 qsort((char *) va, vc, sizeof(*va), simpb_vcmp);
600 /* factor division by most frequently occurring variables first */
601 simp2_divide(equation, np, va[0].v, fc_level);
602 for (i = 1; i < vc; i++) {
603 if (factor_divide(equation, np, va[i].v, 0.0))
604 simp2_divide(equation, np, va[i].v, fc_level);
605 }
606 simp2_divide(equation, np, 0L, fc_level);
607 /* factor all sub-expressions in order of most frequently occurring variables */
608 for (i = 0; i < vc; i++) {
609 while (factor_plus(equation, np, va[i].v, 0.0)) {
610 simp2_divide(equation, np, 0L, fc_level);
611 }
612 }
613 }
614 #if 0 /* Messes up support/gcd.in and tests/quartic.in */
615 make_simple_fractions(equation, np);
616 #endif
617 /* make sure equation side is completely factored */
618 while (factor_divide(equation, np, MATCH_ANY, 0.0)) {
619 simp2_divide(equation, np, MATCH_ANY, fc_level);
620 }
621 while (factor_plus(equation, np, MATCH_ANY, 0.0)) {
622 simp2_divide(equation, np, 0L, fc_level);
623 }
624 simp_ssub(equation, np, MATCH_ANY, 0.0, power_flag, true, fc_level);
625 #if 0 /* Seems to complicate and change constants slightly when enabled. */
626 if (power_flag) {
627 make_simple_fractions(equation, np);
628 factor_power(equation, np);
629 simp_ssub(equation, np, MATCH_ANY, 0.0, power_flag, true, fc_level);
630 }
631 #endif
632 }
633
634 /*
635 * Convert expressions with any algebraic fractions into a single simple fraction.
636 * Used by the fraction command.
637 *
638 * Globals tlhs[] and trhs[] are wiped out.
639 */
640 void
simple_frac_side(equation,np)641 simple_frac_side(equation, np)
642 token_type *equation; /* pointer to the beginning of equation side */
643 int *np; /* pointer to length of equation side */
644 {
645 if (*np == 1) {
646 make_simple_fractions(equation, np);
647 fractions_and_group(equation, np);
648 return;
649 }
650 simp_loop(equation, np);
651 poly_factor(equation, np, true);
652 do {
653 do {
654 do {
655 simp_ssub(equation, np, 0L, 0.0, false, true, 5);
656 } while (poly_gcd_simp(equation, np));
657 } while (uf_power(equation, np));
658 } while (super_factor(equation, np, 3));
659 side_debug(2, equation, *np);
660
661 make_simple_fractions(equation, np);
662 uf_tsimp(equation, np);
663 #if 0 /* Not really needed, but makes it end up the same as the simplify command. */
664 make_simple_fractions(equation, np);
665 uf_power(equation, np);
666 integer_root_simp(equation, np);
667 simpb_side(equation, np, true, true, 3);
668 #endif
669 poly_factor(equation, np, true);
670 simpb_side(equation, np, true, false, 2);
671 simpb_side(equation, np, true, false, 2); /* Added for thoroughness, making sure everything is uf_power()ed. */
672 fractions_and_group(equation, np);
673 }
674
675 /*
676 * This is the slow and thorough simplify of the simplify command.
677 * Applies many equivalent algebraic transformations and their inverses (like unfactor and factor),
678 * then does generalized polynomial simplifications.
679 *
680 * Globals tlhs[] and trhs[] are wiped out.
681 */
682 void
simpa_side(equation,np,quick_flag,frac_flag)683 simpa_side(equation, np, quick_flag, frac_flag)
684 token_type *equation; /* pointer to the beginning of equation side to simplify */
685 int *np; /* pointer to length of the equation side */
686 int quick_flag; /* "simplify quick" option, simpler with no (x+1)^2 expansion */
687 int frac_flag; /* "simplify fraction" option, simplify to the ratio of two polynomials */
688 {
689 int i;
690 int flag, poly_flag = true;
691 jmp_buf save_save;
692
693 if (*np == 1) { /* no need to full simplify a single constant or variable */
694 make_simple_fractions(equation, np);
695 simpb_side(equation, np, true, !frac_flag, 2);
696 return;
697 }
698 debug_string(2, "Simplify input:");
699 side_debug(2, equation, *np);
700 simp_loop(equation, np);
701 #if 1
702 do {
703 simp_ssub(equation, np, 0L, 1.0, false, true, 5);
704 } while (uf_power(equation, np));
705 while (factor_power(equation, np)) {
706 simp_loop(equation, np);
707 }
708 #else
709 simpb_side(equation, np, false, true, 5);
710 #endif
711 if (rationalize_denominators) {
712 rationalize(equation, np);
713 }
714 unsimp_power(equation, np);
715 uf_tsimp(equation, np);
716
717 /* Here is the only place in Mathomatic that we do complete modulus (%) simplification: */
718 uf_pplus(equation, np);
719 uf_repeat(equation, np);
720 do {
721 elim_loop(equation, np);
722 } while (mod_simp(equation, np));
723
724 /* Here we try to simplify out unnecessary negative constants and imaginary numbers: */
725 simp_i(equation, np);
726 unsimp_power(equation, np);
727 uf_times(equation, np);
728 simp_ssub(equation, np, 0L, 1.0, true, true, 5);
729 unsimp_power(equation, np);
730 uf_neg_help(equation, np);
731 uf_tsimp(equation, np);
732 do {
733 do {
734 simp_ssub(equation, np, 0L, 1.0, false, true, 6);
735 } while (uf_power(equation, np));
736 } while (!quick_flag && super_factor(equation, np, 2));
737 if (poly_gcd_simp(equation, np)) {
738 simp_ssub(equation, np, 0L, 1.0, false, true, 6);
739 }
740 side_debug(2, equation, *np);
741 unsimp_power(equation, np);
742 uf_times(equation, np);
743 factorv(equation, np, IMAGINARY);
744 uf_pplus(equation, np);
745 simp_ssub(equation, np, 0L, 1.0, true, false, 5);
746 if (poly_gcd_simp(equation, np)) {
747 factorv(equation, np, IMAGINARY);
748 uf_pplus(equation, np);
749 simp_ssub(equation, np, 0L, 1.0, true, false, 5);
750 }
751 uf_times(equation, np);
752 uf_pplus(equation, np);
753 factor_imaginary(equation, np);
754 uf_power(equation, np);
755 do {
756 do {
757 simp_ssub(equation, np, 0L, 1.0, false, true, 6);
758 } while (uf_power(equation, np));
759 } while (!quick_flag && super_factor(equation, np, 2));
760
761 /* Here we do the greatest expansion; if it fails, do less expansion. */
762 partial_flag = frac_flag;
763 n_tlhs = *np;
764 blt(tlhs, equation, n_tlhs * sizeof(token_type));
765 blt(save_save, jmp_save, sizeof(jmp_save));
766 if ((i = setjmp(jmp_save)) != 0) { /* trap errors */
767 blt(jmp_save, save_save, sizeof(jmp_save));
768 if (i == 13) { /* critical error code */
769 longjmp(jmp_save, i);
770 }
771 /* an error occurred, restore the original expression */
772 *np = n_tlhs;
773 blt(equation, tlhs, n_tlhs * sizeof(token_type));
774 if (i == 14) {
775 debug_string(1, "Simplify not expanding fully, due to oversized expression.");
776 } else {
777 debug_string(0, "Simplify not expanding fully, due to some error.");
778 }
779 partial_flag = true; /* expand less */
780 uf_tsimp(equation, np);
781 } else {
782 if (quick_flag) {
783 uf_tsimp(equation, np);
784 } else {
785 /* expand powers of 2 and higher, might result in error_huge() trap */
786 do {
787 uf_power(equation, np);
788 uf_repeat(equation, np);
789 } while (uf_tsimp(equation, np));
790 }
791 blt(jmp_save, save_save, sizeof(jmp_save));
792 }
793 partial_flag = true;
794
795 simpb_side(equation, np, true, true, 2);
796 debug_string(1, "Simplify result before applying polynomial operations:");
797 side_debug(1, equation, *np);
798 for (flag = false;;) {
799 /* divide top and bottom of fractions by any polynomial GCD found */
800 if (poly_gcd_simp(equation, np)) {
801 flag = false;
802 simpb_side(equation, np, false, true, 3);
803 }
804 /* factor polynomials */
805 if (!flag && poly_factor(equation, np, true)) {
806 flag = true;
807 simpb_side(equation, np, false, true, 3);
808 continue;
809 }
810 /* simplify algebraic fractions with polynomial and smart division */
811 if (!frac_flag && div_remainder(equation, np, poly_flag, quick_flag)) {
812 flag = false;
813 simpb_side(equation, np, false, true, 3);
814 continue;
815 }
816 break;
817 }
818 debug_string(2, "Raw simplify result after applying polynomial operations:");
819 side_debug(2, equation, *np);
820 simp_constant_power(equation, np);
821 simp_ssub(equation, np, 0L, 1.0, true, true, 5);
822 unsimp_power(equation, np);
823 make_simple_fractions(equation, np);
824 factor_power(equation, np);
825 uf_tsimp(equation, np);
826 make_simple_fractions(equation, np);
827 uf_power(equation, np);
828 integer_root_simp(equation, np);
829 simpb_side(equation, np, true, true, 3);
830 poly_factor(equation, np, true);
831 simpb_side(equation, np, true, !frac_flag, 2);
832 }
833
834 /*
835 * This routine is used by the simplify command,
836 * and is the slowest and most thorough simplify of all.
837 * It repeats the simplification until the smallest expression is achieved,
838 * if repeat_flag is true.
839 *
840 * Globals tes[], tlhs[], and trhs[] are wiped out.
841 */
842 void
simpa_repeat_side(equation,np,quick_flag,frac_flag)843 simpa_repeat_side(equation, np, quick_flag, frac_flag)
844 token_type *equation; /* pointer to the beginning of equation side to simplify */
845 int *np; /* pointer to length of the equation side */
846 int quick_flag; /* "simplify quick" option, simpler fractions with no (x+1)^2 expansion */
847 int frac_flag; /* "simplify fraction" option, simplify to the ratio of two polynomials */
848 {
849 if (*np <= 0)
850 return;
851 simpa_side(equation, np, quick_flag, frac_flag);
852 if (repeat_flag && *np > 1) {
853 do {
854 n_tes = *np;
855 blt(tes, equation, n_tes * sizeof(token_type));
856 simpa_side(equation, np, quick_flag, frac_flag);
857 } while (*np < n_tes);
858 if (*np != n_tes) {
859 *np = n_tes;
860 blt(equation, tes, n_tes * sizeof(token_type));
861 }
862 }
863 }
864
865 /*
866 * Completely and repeatedly (if repeat_flag) simplify an equation space to the smallest possible size.
867 *
868 * Globals tes[], tlhs[], and trhs[] are clobbered.
869 */
870 void
simpa_repeat(n,quick_flag,frac_flag)871 simpa_repeat(n, quick_flag, frac_flag)
872 int n; /* equation space number to simplify */
873 int quick_flag; /* "simplify quick" option, simpler fractions with no (x+1)^2 expansion */
874 int frac_flag; /* "simplify fraction" option, simplify to the ratio of two polynomials */
875 {
876 if (empty_equation_space(n))
877 return;
878 simpa_repeat_side(lhs[n], &n_lhs[n], quick_flag, frac_flag);
879 if (n_rhs[n] > 0) {
880 simpa_repeat_side(rhs[n], &n_rhs[n], quick_flag, frac_flag);
881 }
882 }
883
884 /*
885 * This routine is used by the fraction command.
886 * It repeats the simplification until the smallest expression is achieved,
887 * if repeat_flag is true.
888 *
889 * Globals tes[], tlhs[], and trhs[] are wiped out.
890 */
891 void
simple_frac_repeat_side(equation,np)892 simple_frac_repeat_side(equation, np)
893 token_type *equation; /* pointer to the beginning of equation side to simplify */
894 int *np; /* pointer to length of the equation side */
895 {
896 if (*np <= 0)
897 return;
898 simple_frac_side(equation, np);
899 if (repeat_flag) {
900 do {
901 n_tes = *np;
902 blt(tes, equation, n_tes * sizeof(token_type));
903 simple_frac_side(equation, np);
904 } while (*np < n_tes);
905 if (*np != n_tes) {
906 *np = n_tes;
907 blt(equation, tes, n_tes * sizeof(token_type));
908 }
909 }
910 }
911
912 /*
913 * Commonly used quick simplify routine that doesn't factor.
914 *
915 * Return true if factor_times() did something.
916 */
917 int
simp_loop(equation,np)918 simp_loop(equation, np)
919 token_type *equation; /* pointer to the beginning of equation side to simplify */
920 int *np; /* pointer to length of equation side */
921 {
922 int i;
923 int rv = false;
924
925 do {
926 do {
927 do {
928 do {
929 elim_loop(equation, np);
930 } while (simp2_power(equation, np));
931 i = factor_times(equation, np);
932 if (i)
933 rv = true;
934 } while (i);
935 } while (elim_sign(equation, np));
936 } while (subtract_itself(equation, np));
937 return rv;
938 }
939
940 /*
941 * Convert (x^n)^m to x^(n*m) when appropriate.
942 * Fixed to work the same as Maxima when symb_flag is false,
943 * otherwise always converts (x^n)^m to x^(n*m), which sometimes simplifies more.
944 *
945 * Earlier versions of Mathomatic always converted (x^n)^m to x^(n*m),
946 * which is the wrong thing to do, because they are not always equivalent.
947 * That is now the job of the "simplify symbolic" command, which sets symb_flag to true.
948 *
949 * Return true if equation side was modified.
950 */
951 int
simp_pp(equation,np)952 simp_pp(equation, np)
953 token_type *equation;
954 int *np;
955 {
956 int i, j, k;
957 int ilevel, jlevel;
958 int modified = false;
959 double numerator, denominator;
960
961 for (i = 1; i < *np; i += 2) {
962 #if DEBUG
963 if (equation[i].kind != OPERATOR) {
964 error_bug("Bug found in simp_pp(), operators are misplaced.");
965 }
966 #endif
967 if (equation[i].token.operatr != POWER)
968 continue;
969 ilevel = equation[i].level;
970 for (j = i + 2; j < *np; j += 2) {
971 jlevel = equation[j].level;
972 if (jlevel == ilevel - 1 && equation[j].token.operatr == POWER) {
973 if (!symb_flag && (equation[i-1].level != ilevel || equation[i-1].kind != CONSTANT || equation[i-1].token.constant < 0)) {
974 if (jlevel == equation[j+1].level && equation[j+1].kind == CONSTANT) {
975 f_to_fraction(equation[j+1].token.constant, &numerator, &denominator);
976 if (fmod(denominator, 2.0) == 0.0) {
977 if ((i + 2) == j && equation[i+1].kind == CONSTANT) {
978 f_to_fraction(equation[i+1].token.constant, &numerator, &denominator);
979 if (fmod(numerator, 2.0) == 0.0) {
980 break;
981 }
982 } else
983 break;
984 }
985 } else {
986 if ((i + 2) == j && equation[i+1].kind == CONSTANT) {
987 f_to_fraction(equation[i+1].token.constant, &numerator, &denominator);
988 if (fmod(numerator, 2.0) == 0.0) {
989 break;
990 }
991 } else
992 break;
993 }
994 }
995 equation[j].token.operatr = TIMES;
996 for (k = j; k < *np && equation[k].level >= jlevel; k++) {
997 equation[k].level += 2;
998 }
999 for (k = i + 1; k < j; k++)
1000 equation[k].level++;
1001 i -= 2;
1002 modified = true;
1003 break;
1004 }
1005 if (jlevel <= ilevel)
1006 break;
1007 }
1008 }
1009 return modified;
1010 }
1011
1012 /*
1013 * Simplify surds like 12^(1/2) to 2*3^(1/2).
1014 *
1015 * Return true if equation side was modified.
1016 */
1017 int
integer_root_simp(equation,np)1018 integer_root_simp(equation, np)
1019 token_type *equation;
1020 int *np;
1021 {
1022 int modified = false;
1023 int i, j;
1024 int level;
1025 double d1, d2, numerator, denominator;
1026
1027 for (i = 1; (i + 3) < *np; i += 2) {
1028 #if DEBUG
1029 if (equation[i].kind != OPERATOR) {
1030 error_bug("Bug found in integer_root_simp(), operators are misplaced.");
1031 }
1032 #endif
1033 if (equation[i].token.operatr == POWER) {
1034 level = equation[i].level;
1035 if (equation[i-1].level == level
1036 && equation[i+1].level == level + 1
1037 && equation[i+2].level == level + 1
1038 && equation[i+3].level == level + 1
1039 && equation[i+2].token.operatr == DIVIDE
1040 && equation[i-1].kind == CONSTANT
1041 && equation[i+1].kind == CONSTANT
1042 && equation[i+3].kind == CONSTANT) {
1043 if ((i + 4) < *np && equation[i+4].level >= level)
1044 continue;
1045 numerator = equation[i+1].token.constant;
1046 if (numerator > 50.0 || numerator < 1.0 || fmod(numerator, 1.0) != 0.0)
1047 continue;
1048 denominator = equation[i+3].token.constant;
1049 if (denominator > 50.0 || denominator < 2.0 || fmod(denominator, 1.0) != 0.0)
1050 continue;
1051 errno = 0;
1052 d2 = pow(equation[i-1].token.constant, numerator);
1053 if (errno) {
1054 continue;
1055 }
1056 if (!factor_one(d2))
1057 continue;
1058 d1 = 1.0;
1059 for (j = 0; j < uno; j++) {
1060 if (unique[j] > 0.0) {
1061 while (ucnt[j] >= denominator) {
1062 d1 *= unique[j];
1063 ucnt[j] -= denominator;
1064 }
1065 }
1066 }
1067 if (d1 == 1.0)
1068 continue;
1069 if (*np + 2 > n_tokens) {
1070 error_huge();
1071 }
1072 equation[i+1].token.constant = 1.0;
1073 equation[i-1].token.constant = multiply_out_unique();
1074 for (j = i - 1; j < (i + 4); j++) {
1075 equation[j].level++;
1076 }
1077 blt(&equation[i+1], &equation[i-1], (*np - (i - 1)) * sizeof(token_type));
1078 *np += 2;
1079 equation[i-1].level = level;
1080 equation[i-1].kind = CONSTANT;
1081 equation[i-1].token.constant = d1;
1082 equation[i].level = level;
1083 equation[i].kind = OPERATOR;
1084 equation[i].token.operatr = TIMES;
1085 modified = true;
1086 i += 4;
1087 }
1088 }
1089 }
1090 return modified;
1091 }
1092
1093 /*
1094 * Simplify constant^(constant*x) to (constant^constant)^x.
1095 * I am not sure if this is a good thing, because it changes the base of the exponent.
1096 * The tests don't require this routine.
1097 *
1098 * Return true if equation side was modified.
1099 */
1100 int
simp_constant_power(equation,np)1101 simp_constant_power(equation, np)
1102 token_type *equation;
1103 int *np;
1104 {
1105 int i, j;
1106 int level;
1107 int modified = false;
1108
1109 for (i = 1; i < *np; i += 2) {
1110 if (equation[i].token.operatr != POWER) {
1111 continue;
1112 }
1113 level = equation[i].level;
1114 if (equation[i-1].level != level || equation[i-1].kind != CONSTANT)
1115 continue;
1116 if (equation[i-1].token.constant < 0 && !symb_flag)
1117 continue;
1118 if (equation[i+1].level != level + 1 || equation[i+1].kind != CONSTANT
1119 || equation[i+1].token.constant == 1.0)
1120 continue;
1121 j = i + 2;
1122 if (j >= *np || equation[j].level != level + 1)
1123 continue;
1124 switch (equation[j].token.operatr) {
1125 case TIMES:
1126 break;
1127 case DIVIDE:
1128 if (*np + 2 > n_tokens) {
1129 error_huge();
1130 }
1131 blt(&equation[j+2], &equation[j], (*np - j) * sizeof(token_type));
1132 *np += 2;
1133 equation[j+1].level = level + 1;
1134 equation[j+1].kind = CONSTANT;
1135 equation[j+1].token.constant = 1.0;
1136 break;
1137 default:
1138 continue;
1139 }
1140 equation[j].level = level;
1141 equation[j].token.operatr = POWER;
1142 equation[i-1].level++;
1143 equation[i].level++;
1144 modified = true;
1145 }
1146 return modified;
1147 }
1148
1149 /*
1150 * Convert x^-y to 1/(x^y).
1151 *
1152 * Return true if equation side was modified.
1153 */
1154 int
simp2_power(equation,np)1155 simp2_power(equation, np)
1156 token_type *equation;
1157 int *np;
1158 {
1159 int i, i1, j, k;
1160 int level;
1161 int op;
1162 int modified = false;
1163
1164 for (i = 1; i < *np; i += 2) {
1165 if (equation[i].token.operatr == POWER) {
1166 level = equation[i].level;
1167 op = 0;
1168 k = -1;
1169 for (j = i + 1; j < *np && equation[j].level >= level; j++) {
1170 if (equation[j].level == level + 1) {
1171 if (equation[j].kind == OPERATOR) {
1172 op = equation[j].token.operatr;
1173 } else if (equation[j].kind == CONSTANT && equation[j].token.constant < 0.0) {
1174 k = j;
1175 }
1176 }
1177 }
1178 if (j - i <= 2 && equation[i+1].kind == CONSTANT && equation[i+1].token.constant < 0.0) {
1179 k = i + 1;
1180 } else if (k < 0)
1181 continue;
1182 switch (op) {
1183 case 0:
1184 case TIMES:
1185 case DIVIDE:
1186 if (*np + 2 > n_tokens) {
1187 error_huge();
1188 }
1189 equation[k].token.constant = -equation[k].token.constant;
1190 for (k = i - 2;; k--) {
1191 if (k < 0 || equation[k].level < level)
1192 break;
1193 }
1194 k++;
1195 for (i1 = k; i1 < j; i1++)
1196 equation[i1].level++;
1197 blt(&equation[k+2], &equation[k], (*np - k) * sizeof(token_type));
1198 *np += 2;
1199 equation[k].level = level;
1200 equation[k].kind = CONSTANT;
1201 equation[k].token.constant = 1.0;
1202 k++;
1203 equation[k].level = level;
1204 equation[k].kind = OPERATOR;
1205 equation[k].token.operatr = DIVIDE;
1206 modified = true;
1207 }
1208 }
1209 }
1210 return modified;
1211 }
1212
1213 /*
1214 * fmod(3) in the standard C math library has some problems.
1215 * Hopefully this fixes them.
1216 */
1217 double
fixed_fmod(k1,k2)1218 fixed_fmod(k1, k2)
1219 double k1, k2;
1220 {
1221 double d;
1222
1223 if (k2 == 0 || !isfinite(k1) || !isfinite(k2) || (fmod(k1, 1.0) == 0 && fmod(k2, 1.0) == 0)) {
1224 k1 = fmod(k1, k2);
1225 } else {
1226 /* fmod() doesn't work well with fractional amounts. */
1227 /* Another way to get the remainder of division: */
1228 k1 = modf(k1 / k2, &d) * k2;
1229 }
1230 return k1;
1231 }
1232
1233 /*
1234 * Combine two or more constants on the same level of parentheses.
1235 * If "iflag" is false, don't produce imaginary numbers.
1236 *
1237 * Return true if equation side was modified.
1238 */
1239 int
combine_constants(equation,np,iflag)1240 combine_constants(equation, np, iflag)
1241 token_type *equation; /* pointer to the beginning of equation side */
1242 int *np; /* pointer to length of equation side */
1243 int iflag; /* produce imaginary numbers flag */
1244 {
1245 return const_recurse(equation, np, 0, 1, iflag);
1246 }
1247
1248 /*
1249 * Do the floating point arithmetic for Mathomatic.
1250 *
1251 * Return true if successful.
1252 * domain_check must be set to false after this.
1253 */
1254 int
calc(op1p,k1p,op2,k2)1255 calc(op1p, k1p, op2, k2)
1256 int *op1p; /* Pointer to operator 1, which comes immediately before operand 1 if this operator exists. */
1257 double *k1p; /* Pointer to operand 1 and where to store the result of the calculation. */
1258 int op2; /* Operator 2; always exists and comes immediately before operand 2, and usually after operand 1. */
1259 double k2; /* Operand 2; ignored for unary operators. */
1260 {
1261 #if _REENTRANT
1262 int sign = 1;
1263 #endif
1264 int op1;
1265 double d, d1, d2;
1266
1267 domain_check = false;
1268 errno = 0;
1269 if (op1p) {
1270 op1 = *op1p;
1271 } else {
1272 op1 = 0;
1273 }
1274 switch (op2) {
1275 case PLUS:
1276 case MINUS:
1277 if (op1 == MINUS)
1278 d = -(*k1p);
1279 else
1280 d = *k1p;
1281 d1 = fabs(d) * epsilon;
1282 if (op2 == PLUS) {
1283 d += k2;
1284 } else {
1285 d -= k2;
1286 }
1287 if (fabs(d) < d1)
1288 d = 0.0;
1289 if (op1 == 0) {
1290 *k1p = d;
1291 } else {
1292 if (d >= 0.0) {
1293 *op1p = PLUS;
1294 *k1p = d;
1295 } else {
1296 *op1p = MINUS;
1297 *k1p = -d;
1298 }
1299 }
1300 break;
1301 case TIMES:
1302 case DIVIDE:
1303 if (op1 == 0)
1304 op1 = TIMES;
1305 if (op1 == op2) {
1306 *k1p *= k2;
1307 } else {
1308 if (op1 == DIVIDE) {
1309 check_divide_by_zero(*k1p);
1310 *k1p = k2 / *k1p;
1311 *op1p = TIMES;
1312 } else if (op2 == DIVIDE) {
1313 check_divide_by_zero(k2);
1314 *k1p = *k1p / k2;
1315 }
1316 }
1317 break;
1318 case IDIVIDE:
1319 check_divide_by_zero(k2);
1320 modf(*k1p / k2, k1p);
1321 break;
1322 case MODULUS:
1323 if (k2 == 0) {
1324 warning(_("Modulo 0 encountered."));
1325 }
1326 *k1p = fixed_fmod(*k1p, k2);
1327 if (modulus_mode && *k1p < 0.0) {
1328 *k1p += fabs(k2); /* make positive */
1329 }
1330 if (modulus_mode == 1 && k2 < 0.0 && *k1p > 0.0) {
1331 *k1p += k2; /* make negative */
1332 }
1333 /* modulus_mode == 0 result same sign as dividend,
1334 modulus_mode == 1 result same sign as divisor,
1335 modulus_mode == 2 result is always positive or zero */
1336 break;
1337 case POWER:
1338 if (*k1p < 0.0 && fmod(k2, 1.0) != 0.0) {
1339 /* it's probably imaginary; pow() will give a domain error, so skip these calculations */
1340 break;
1341 }
1342 domain_check = true;
1343 if (*k1p == 0.0 && k2 == 0.0) {
1344 warning(_("0^0 encountered, might be considered indeterminate."));
1345 d = 1.0; /* 0^0 = 1 */
1346 } else if (*k1p == 0 && k2 < 0.0) {
1347 warning(_("Divide by zero (0 raised to negative power)."));
1348 d = INFINITY;
1349 } else {
1350 d = pow(*k1p, k2);
1351 if (preserve_surds && !approximate_roots) {
1352 if (isfinite(k2) && fmod(k2, 1.0) != 0.0 && f_to_fraction(*k1p, &d1, &d2)) {
1353 if (!f_to_fraction(d, &d1, &d2)) {
1354 domain_check = false;
1355 return false;
1356 }
1357 }
1358 }
1359 }
1360 #if 1
1361 if (errno == ERANGE) { /* preserve overflowed powers rather than aborting with an error message */
1362 domain_check = false;
1363 return false;
1364 }
1365 #endif
1366 check_err();
1367 if (domain_check)
1368 *k1p = d;
1369 break;
1370 case FACTORIAL:
1371 #if NOGAMMA
1372 #warning "Using integer factorial only."
1373 /* calculate integer factorial if gamma() functions don't exist */
1374 if (*k1p > 170.0 || *k1p < 0.0 || fmod(*k1p, 1.0) != 0.0) {
1375 return false;
1376 }
1377 d = 1.0;
1378 for (d1 = 2.0; d1 <= *k1p; d1 += 1.0) {
1379 d *= d1;
1380 }
1381 #else
1382 #if _XOPEN_SOURCE >= 600 || _ISOC99_SOURCE || __USE_ISOC99 || MINGW || __APPLE__ || USE_TGAMMA
1383 /* Use true gamma for factorial if available, it is the nicest gamma() function. */
1384 /* Problem is I don't know how to properly tell if it is available. */
1385 d = tgamma(*k1p + 1.0);
1386 if (errno) { /* don't evaluate if overflow */
1387 return false;
1388 }
1389 #else
1390 #if !_REENTRANT
1391 #warning "Using lgamma(3); Not a problem, but tgamma(3) is more direct, if available."
1392 d = exp(lgamma(*k1p + 1.0)) * signgam;
1393 #else /* use re-entrant version: */
1394 #warning "Using lgamma_r(3); Not a problem, but tgamma(3) is more direct, if available."
1395 d = exp(lgamma_r(*k1p + 1.0, &sign));
1396 d *= sign;
1397 #endif
1398 /* Just compile with -DUSE_TGAMMA, if tgamma(3) is available and not being used. */
1399 if (errno) { /* don't evaluate if overflow */
1400 return false;
1401 }
1402 #endif
1403 #endif
1404 *k1p = d;
1405 break;
1406 default:
1407 return false;
1408 }
1409 return true;
1410 }
1411
1412 static int
const_recurse(equation,np,loc,level,iflag)1413 const_recurse(equation, np, loc, level, iflag)
1414 token_type *equation;
1415 int *np, loc, level, iflag;
1416 {
1417 int loc1, old_loc;
1418 int const_count = 0;
1419 int op;
1420 int modified = false;
1421 double d1, d2, d3, numerator, denominator;
1422 complexs cv, p;
1423
1424 loc1 = old_loc = loc;
1425 for (;; loc++) {
1426 beginning:
1427 if (loc >= *np || equation[loc].level < level) {
1428 if (loc - old_loc == 1) /* decrement the level of parentheses if only one constant left */
1429 equation[old_loc].level = max(level - 1, 1);
1430 return modified;
1431 }
1432 if (equation[loc].level > level) {
1433 modified |= const_recurse(equation, np, loc, level + 1, iflag);
1434 for (; loc < *np && equation[loc].level > level; loc++)
1435 ;
1436 goto beginning;
1437 }
1438 if (equation[loc].kind == CONSTANT) {
1439 if (const_count == 0) {
1440 loc1 = loc;
1441 const_count++;
1442 continue;
1443 }
1444 op = equation[loc-1].token.operatr;
1445 d1 = equation[loc1].token.constant;
1446 d2 = equation[loc].token.constant;
1447 if (calc((loc1 <= old_loc) ? NULL : &equation[loc1-1].token.operatr, &d1, op, d2)) {
1448 if (op == POWER && !domain_check) {
1449 if (!f_to_fraction(d2, &numerator, &denominator)) { /* if irrational power */
1450 if (!iflag || (preserve_surds && !approximate_roots))
1451 return modified;
1452 cv.re = d1;
1453 cv.im = 0.0;
1454 p.re = d2;
1455 p.im = 0.0;
1456 cv = complex_pow(cv, p);
1457 if (*np + 2 > n_tokens) {
1458 error_huge();
1459 }
1460 blt(&equation[loc1+2], &equation[loc1], (*np - loc1) * sizeof(token_type));
1461 *np += 2;
1462 equation[loc1].level = level;
1463 equation[loc1].kind = CONSTANT;
1464 equation[loc1].token.constant = cv.re;
1465 loc1++;
1466 equation[loc1].level = level;
1467 equation[loc1].kind = OPERATOR;
1468 equation[loc1].token.operatr = PLUS;
1469 level++;
1470 equation[loc].level = level;
1471 equation[loc].kind = VARIABLE;
1472 equation[loc].token.variable = IMAGINARY;
1473 loc++;
1474 equation[loc].level = level;
1475 equation[loc].kind = OPERATOR;
1476 equation[loc].token.operatr = TIMES;
1477 loc++;
1478 equation[loc].level = level;
1479 equation[loc].kind = CONSTANT;
1480 equation[loc].token.constant = cv.im;
1481 return true;
1482 }
1483 errno = 0;
1484 d3 = pow(-d1, d2);
1485 check_err();
1486 if (!always_positive(denominator)) {
1487 if (*np + 2 > n_tokens) {
1488 error_huge();
1489 }
1490 blt(&equation[loc1+2], &equation[loc1], (*np - loc1) * sizeof(token_type));
1491 *np += 2;
1492 equation[loc1].level = level + 1;
1493 equation[loc1].kind = CONSTANT;
1494 equation[loc1].token.constant = -d1;
1495 loc1++;
1496 equation[loc1].level = level + 1;
1497 equation[loc1].kind = OPERATOR;
1498 equation[loc1].token.operatr = POWER;
1499 equation[loc].level = level + 1;
1500 equation[loc].kind = CONSTANT;
1501 equation[loc].token.constant = d2;
1502 loc++;
1503 equation[loc].level = level;
1504 equation[loc].kind = OPERATOR;
1505 equation[loc].token.operatr = TIMES;
1506 loc++;
1507 equation[loc].level = level;
1508 equation[loc].kind = CONSTANT;
1509 if (always_positive(numerator)) {
1510 equation[loc].token.constant = 1.0;
1511 } else {
1512 equation[loc].token.constant = -1.0;
1513 }
1514 return true;
1515 }
1516 if (!iflag)
1517 return modified;
1518 if (*np + 2 > n_tokens) {
1519 error_huge();
1520 }
1521 blt(&equation[loc1+2], &equation[loc1], (*np - loc1) * sizeof(token_type));
1522 *np += 2;
1523 if (d2 == 0.5) {
1524 equation[loc1].level = level + 1;
1525 equation[loc1].kind = CONSTANT;
1526 equation[loc1].token.constant = -d1;
1527 loc1++;
1528 equation[loc1].level = level + 1;
1529 equation[loc1].kind = OPERATOR;
1530 equation[loc1].token.operatr = POWER;
1531 equation[loc].level = level + 1;
1532 equation[loc].kind = CONSTANT;
1533 equation[loc].token.constant = d2;
1534 loc++;
1535 equation[loc].level = level;
1536 equation[loc].kind = OPERATOR;
1537 equation[loc].token.operatr = TIMES;
1538 loc++;
1539 equation[loc].level = level;
1540 equation[loc].kind = VARIABLE;
1541 equation[loc].token.variable = IMAGINARY;
1542 } else {
1543 equation[loc1].level = level;
1544 equation[loc1].kind = CONSTANT;
1545 equation[loc1].token.constant = d3;
1546 loc1++;
1547 equation[loc1].level = level;
1548 equation[loc1].kind = OPERATOR;
1549 equation[loc1].token.operatr = TIMES;
1550 level++;
1551 equation[loc].level = level;
1552 equation[loc].kind = VARIABLE;
1553 equation[loc].token.variable = IMAGINARY;
1554 loc++;
1555 equation[loc].level = level;
1556 equation[loc].kind = OPERATOR;
1557 equation[loc].token.operatr = POWER;
1558 loc++;
1559 equation[loc].level = level;
1560 equation[loc].kind = CONSTANT;
1561 equation[loc].token.constant = d2 * 2.0;
1562 }
1563 return true;
1564 } else {
1565 equation[loc1].token.constant = d1;
1566 modified = true;
1567 domain_check = false;
1568 blt(&equation[loc-1], &equation[loc+1], (*np - (loc + 1)) * sizeof(token_type));
1569 *np -= 2;
1570 loc -= 2;
1571 }
1572 } else {
1573 domain_check = false;
1574 }
1575 }
1576 }
1577 }
1578
1579 /*
1580 * Eliminate or fix operations involving constants that can be simplified or normalized.
1581 *
1582 * Return true if equation side was significantly modified.
1583 */
1584 int
elim_k(equation,np)1585 elim_k(equation, np)
1586 token_type *equation; /* equation side pointer */
1587 int *np; /* pointer to length of equation side */
1588 {
1589 token_type *p1, *p2, *p3, *p4;
1590 token_type *ep; /* end pointer */
1591 int modified = false;
1592 int level;
1593 double d, numerator, denominator;
1594 int flag;
1595
1596 for (p1 = equation + 1;;) {
1597 ep = &equation[*np];
1598 if (p1 >= ep)
1599 break;
1600 if (p1->kind != OPERATOR) {
1601 p1++;
1602 continue;
1603 }
1604 level = p1->level;
1605 switch (p1->token.operatr) {
1606 case PLUS:
1607 case MINUS:
1608 /* Fix addition/subtraction of negative, zero, or infinity constants. */
1609 p2 = p1 + 1;
1610 if (p1 + 2 < ep && (p1 + 2)->level == level + 1
1611 && ((p1 + 2)->token.operatr == TIMES || (p1 + 2)->token.operatr == DIVIDE)
1612 && p2->kind == CONSTANT && p2->token.constant < 0.0) {
1613 if (p1->token.operatr == PLUS)
1614 p1->token.operatr = MINUS;
1615 else
1616 p1->token.operatr = PLUS;
1617 p2->token.constant = -(p2->token.constant);
1618 }
1619 if (p2->level == level && p2->kind == CONSTANT) {
1620 if (p2->token.constant < 0.0) {
1621 if (p1->token.operatr == PLUS)
1622 p1->token.operatr = MINUS;
1623 else
1624 p1->token.operatr = PLUS;
1625 p2->token.constant = -p2->token.constant;
1626 }
1627 if (p2->token.constant == 0.0) {
1628 blt(p1, p1 + 2, (char *) ep - (char *) (p1 + 2));
1629 *np -= 2;
1630 modified = true;
1631 continue;
1632 }
1633 }
1634 if ((p1 - 1)->level == level && (p1 - 1)->kind == CONSTANT && isinf((p1 - 1)->token.constant)) {
1635 p2 = p1 - 1;
1636 }
1637 if (p2->level == level && p2->kind == CONSTANT && isinf(p2->token.constant)) {
1638 flag = false;
1639 for (p3 = p1;; p3--) {
1640 if (p3->level < level) {
1641 p3++;
1642 break;
1643 }
1644 if (p3->kind == CONSTANT && p3 != p2 && !isfinite(p3->token.constant)) {
1645 flag = true;
1646 }
1647 if (p3 == equation)
1648 break;
1649 }
1650 for (p4 = p1; p4 < ep && p4->level >= level; p4++) {
1651 if (p4->kind == CONSTANT && p4 != p2 && !isfinite(p4->token.constant)) {
1652 flag = true;
1653 }
1654 }
1655 if (!flag) { /* no other infinities on level */
1656 if (p2 > p3 && (p2 - 1)->token.operatr == MINUS) {
1657 p2->token.constant = -(p2->token.constant);
1658 }
1659 blt(p2 + 1, p4, (char *) ep - (char *) p4);
1660 *np -= p4 - (p2 + 1);
1661 ep = &equation[*np];
1662 blt(p3, p2, (char *) ep - (char *) p2);
1663 *np -= p2 - p3;
1664 return true;
1665 }
1666 }
1667 break;
1668 }
1669 p2 = p1 - 1;
1670 switch (p1->token.operatr) {
1671 case PLUS:
1672 /* remove 0+ */
1673 if (p2->level == level && p2->kind == CONSTANT && p2->token.constant == 0.0) {
1674 blt(p2, p1 + 1, (char *) ep - (char *) (p1 + 1));
1675 *np -= 2;
1676 modified = true;
1677 continue;
1678 }
1679 break;
1680 case MINUS:
1681 /* 0-x to (-1*x) */
1682 if (p2->level == level && p2->kind == CONSTANT && p2->token.constant == 0.0) {
1683 if (p2 == equation || (p2 - 1)->level < level) {
1684 p2->token.constant = -1.0;
1685 p1->token.operatr = TIMES;
1686 binary_parenthesize(equation, *np, p1 - equation);
1687 modified = true;
1688 continue;
1689 }
1690 }
1691 break;
1692 case TIMES:
1693 if (p2->level == level && p2->kind == CONSTANT) {
1694 if (p2->token.constant == 0.0) {
1695 /* Replace 0*x with 0. */
1696 for (p2 = p1 + 2; p2 < ep; p2 += 2) {
1697 if (p2->level < level)
1698 break;
1699 }
1700 blt(p1, p2, (char *) ep - (char *) p2);
1701 *np -= p2 - p1;
1702 modified = true;
1703 continue;
1704 }
1705 if (fabs(p2->token.constant - 1.0) <= epsilon) {
1706 /* Replace 1*x with x. */
1707 blt(p2, p1 + 1, (char *) ep - (char *) (p1 + 1));
1708 *np -= 2;
1709 modified = true;
1710 continue;
1711 }
1712 }
1713 if ((p1 + 1)->level == level && (p1 + 1)->kind == CONSTANT) {
1714 /* Move any constant to the beginning of a multiplicative sub-expression. */
1715 d = (p1 + 1)->token.constant;
1716 for (p2 = p1 - 1; p2 > equation; p2--) {
1717 if ((p2 - 1)->level < level)
1718 break;
1719 }
1720 if (p2->level == level && p2->kind == CONSTANT) {
1721 /* already a constant at the beginning, so don't move anything */
1722 break;
1723 }
1724 blt(p2 + 2, p2, (char *) p1 - (char *) p2);
1725 p2->level = level;
1726 p2->kind = CONSTANT;
1727 p2->token.constant = d;
1728 (p2 + 1)->level = level;
1729 (p2 + 1)->kind = OPERATOR;
1730 (p2 + 1)->token.operatr = TIMES;
1731 if (p2 > equation) {
1732 p1 = p2 - 1;
1733 } else {
1734 p1 = equation + 1;
1735 }
1736 continue;
1737 }
1738 break;
1739 case DIVIDE:
1740 if (p2->level == level && p2->kind == CONSTANT && p2->token.constant == 0.0) {
1741 /* Replace 0/x with 0. */
1742 for (p2 = p1 + 2; p2 < ep; p2 += 2) {
1743 if (p2->level < level)
1744 break;
1745 }
1746 blt(p1, p2, (char *) ep - (char *) p2);
1747 *np -= p2 - p1;
1748 modified = true;
1749 continue;
1750 }
1751 p2 = p1 + 1;
1752 if (p2->level == level && p2->kind == CONSTANT) {
1753 /* Replace division by a constant with times its reciprocal. */
1754 f_to_fraction(p2->token.constant, &numerator, &denominator);
1755 check_divide_by_zero(numerator);
1756 p2->token.constant = denominator / numerator;
1757 p1->token.operatr = TIMES;
1758 continue;
1759 }
1760 if (p2->level == level && p2->kind == VARIABLE && (p2->token.variable & VAR_MASK) == SIGN) {
1761 /* Replace division by a sign variable with times the sign variable. */
1762 p1->token.operatr = TIMES;
1763 continue;
1764 }
1765 break;
1766 case MODULUS:
1767 case IDIVIDE:
1768 if (p2->level == level && p2->kind == CONSTANT && p2->token.constant == 0.0) {
1769 /* Replace 0%x with 0. */
1770 for (p2 = p1 + 2; p2 < ep; p2 += 2) {
1771 if (p2->level < level)
1772 break;
1773 }
1774 blt(p1, p2, (char *) ep - (char *) p2);
1775 *np -= p2 - p1;
1776 modified = true;
1777 continue;
1778 }
1779 if (p1->token.operatr == MODULUS) {
1780 if ((p1 + 1)->level == level && (p1 + 1)->kind == CONSTANT) {
1781 d = fabs((p1 + 1)->token.constant);
1782 if (d > epsilon && fmod(1.0 / d, 1.0) == 0.0) {
1783 for (p2 = p1 - 1; p2 > equation; p2--) {
1784 if ((p2 - 1)->level < level)
1785 break;
1786 }
1787 if (is_integer_expr(p2, p1 - p2)) {
1788 blt(p2, p1 + 1, (char *) ep - (char *) (p1 + 1));
1789 *np -= (p1 + 1) - p2;
1790 p2->token.constant = 0.0;
1791 if (p2 > equation) {
1792 p1 = p2 - 1;
1793 } else {
1794 p1 = equation + 1;
1795 }
1796 modified = true;
1797 continue;
1798 }
1799 }
1800 }
1801 }
1802 break;
1803 case POWER:
1804 if (p2->level == level && p2->kind == CONSTANT) {
1805 if (p2->token.constant == 1.0) {
1806 /* Replace 1^x with 1. */
1807 for (p2 = p1 + 2; p2 < ep; p2 += 2) {
1808 if (p2->level <= level)
1809 break;
1810 }
1811 blt(p1, p2, (char *) ep - (char *) p2);
1812 *np -= p2 - p1;
1813 modified = true;
1814 continue;
1815 }
1816 }
1817 p2 = p1 + 1;
1818 if (p2->level == level && p2->kind == CONSTANT) {
1819 if (p2->token.constant == 0.0) {
1820 /* Replace x^0 with 1. */
1821 for (p2 = p1 - 1; p2 > equation; p2--) {
1822 if ((p2 - 1)->level <= level)
1823 break;
1824 }
1825 blt(p2, p1 + 1, (char *) ep - (char *) (p1 + 1));
1826 *np -= (p1 + 1) - p2;
1827 p2->token.constant = 1.0;
1828 p1 = p2 + 1;
1829 modified = true;
1830 continue;
1831 }
1832 if (fabs(p2->token.constant - 1.0) <= epsilon) {
1833 /* Replace x^1 with x. */
1834 blt(p1, p1 + 2, (char *) ep - (char *) (p1 + 2));
1835 *np -= 2;
1836 modified = true;
1837 continue;
1838 }
1839 }
1840 break;
1841 }
1842 p1 += 2;
1843 }
1844 return modified;
1845 }
1846
1847 /*
1848 * Compare two sub-expressions for equality.
1849 * This is quick and low level and does not simplify or modify the expressions.
1850 *
1851 * If equal, return true with *diff_signp = 0.
1852 * if p1 * -1 == p2, return true with *diff_signp = 1.
1853 * Otherwise return false.
1854 */
1855 int
se_compare(p1,n1,p2,n2,diff_signp)1856 se_compare(p1, n1, p2, n2, diff_signp)
1857 token_type *p1; /* first sub-expression pointer */
1858 int n1; /* first sub-expression length */
1859 token_type *p2; /* second sub-expression pointer */
1860 int n2; /* second sub-expression length */
1861 int *diff_signp; /* different sign flag pointer */
1862 {
1863 int l1, l2;
1864 int rv;
1865 #if DEBUG
1866 int rv_should_be_false = false;
1867
1868 if (n1 < 1 || n2 < 1 || (n1 & 1) != 1 || (n2 & 1) != 1 || diff_signp == NULL || p1 == NULL || p2 == NULL) {
1869 error_bug("Programming error in call to se_compare().");
1870 }
1871 #endif
1872 if (((n1 > n2) ? ((n1 + 1) / (n2 + 1)) : ((n2 + 1) / (n1 + 1))) > 3) {
1873 /* expressions are grossly different in size, no need to compare, they are different */
1874 #if DEBUG
1875 rv_should_be_false = true;
1876 #else
1877 *diff_signp = false;
1878 return false;
1879 #endif
1880 }
1881 /* Find the proper ground levels of parentheses for the two sub-expressions: */
1882 l1 = min_level(p1, n1);
1883 l2 = min_level(p2, n2);
1884 rv = compare_recurse(p1, n1, l1, p2, n2, l2, diff_signp);
1885 #if DEBUG
1886 if (rv && rv_should_be_false) {
1887 error_bug("Expression compare optimization failed in se_compare().");
1888 }
1889 #endif
1890 return rv;
1891 }
1892
1893 /*
1894 * Recursively compare each parenthesized sub-expression.
1895 * This is the most used function in Mathomatic.
1896 * Optimizing or streamlining this would greatly speed things up.
1897 */
1898 static int
compare_recurse(p1,n1,l1,p2,n2,l2,diff_signp)1899 compare_recurse(p1, n1, l1, p2, n2, l2, diff_signp)
1900 token_type *p1;
1901 int n1, l1;
1902 token_type *p2;
1903 int n2, l2, *diff_signp;
1904 {
1905 #define compare_epsilon epsilon
1906
1907 token_type *pv1, *ep1, *ep2;
1908 int i, j;
1909 int len;
1910 int first;
1911 int oc2; /* operand count 2 */
1912 token_type *opa2[MAX_COMPARE_TERMS]; /* operand pointer array 2 */
1913 char used[MAX_COMPARE_TERMS]; /* operand used flag array 2 */
1914 int last_op1, op1 = 0, op2 = 0;
1915 int diff_op = false;
1916 double d1, c1, c2;
1917
1918 *diff_signp = false;
1919 if (n1 == 1 && n2 == 1) {
1920 if (p1->kind != p2->kind) {
1921 return false;
1922 }
1923 switch (p1->kind) {
1924 case VARIABLE:
1925 if (sign_cmp_flag && (p1->token.variable & VAR_MASK) == SIGN) {
1926 return((p2->token.variable & VAR_MASK) == SIGN);
1927 } else {
1928 return(p1->token.variable == p2->token.variable);
1929 }
1930 break;
1931 case CONSTANT:
1932 c1 = p1->token.constant;
1933 c2 = p2->token.constant;
1934 if (c1 == c2) {
1935 return true;
1936 } else if (c1 == -c2) {
1937 *diff_signp = true;
1938 return true;
1939 }
1940 d1 = fabs(c1) * compare_epsilon;
1941 if (fabs(c1 - c2) < d1) {
1942 return true;
1943 }
1944 if (fabs(c1 + c2) < d1) {
1945 *diff_signp = true;
1946 return true;
1947 }
1948 break;
1949 case OPERATOR:
1950 error_bug("Programming error in call to compare_recurse().");
1951 break;
1952 }
1953 return false;
1954 }
1955 #if DEBUG
1956 if (n1 < 1 || n2 < 1 || (n1 & 1) != 1 || (n2 & 1) != 1) {
1957 error_bug("Programming error in call to compare_recurse().");
1958 }
1959 #endif
1960 ep1 = &p1[n1];
1961 ep2 = &p2[n2];
1962 for (pv1 = p1 + 1; pv1 < ep1; pv1 += 2) {
1963 if (pv1->level == l1) {
1964 op1 = pv1->token.operatr;
1965 break;
1966 }
1967 }
1968 for (pv1 = p2 + 1; pv1 < ep2; pv1 += 2) {
1969 if (pv1->level == l2) {
1970 op2 = pv1->token.operatr;
1971 break;
1972 }
1973 }
1974 if (op2 == 0) {
1975 if (op1 != TIMES && op1 != DIVIDE)
1976 return false;
1977 goto no_op2;
1978 }
1979 switch (op1) {
1980 case PLUS:
1981 case MINUS:
1982 if (op2 != PLUS && op2 != MINUS)
1983 diff_op = true;
1984 break;
1985 case 0:
1986 if (op2 != TIMES && op2 != DIVIDE)
1987 return false;
1988 break;
1989 case TIMES:
1990 case DIVIDE:
1991 if (op2 != TIMES && op2 != DIVIDE)
1992 diff_op = true;
1993 break;
1994 default:
1995 if (op2 != op1)
1996 diff_op = true;
1997 break;
1998 }
1999 if (diff_op) {
2000 if (p1->kind == CONSTANT && p1->level == l1 && op1 == TIMES) {
2001 if (fabs(fabs(p1->token.constant) - 1.0) <= compare_epsilon) {
2002 if (!compare_recurse(p1 + 2, n1 - 2, min_level(p1 + 2, n1 - 2), p2, n2, l2, diff_signp)) {
2003 return false;
2004 }
2005 if (p1->token.constant < 0.0) {
2006 *diff_signp ^= true;
2007 }
2008 return true;
2009 }
2010 }
2011 if (p2->kind == CONSTANT && p2->level == l2 && op2 == TIMES) {
2012 if (fabs(fabs(p2->token.constant) - 1.0) <= compare_epsilon) {
2013 if (!compare_recurse(p1, n1, l1, p2 + 2, n2 - 2, min_level(p2 + 2, n2 - 2), diff_signp)) {
2014 return false;
2015 }
2016 if (p2->token.constant < 0.0) {
2017 *diff_signp ^= true;
2018 }
2019 return true;
2020 }
2021 }
2022 return false;
2023 }
2024 no_op2:
2025 opa2[0] = p2;
2026 used[0] = false;
2027 oc2 = 1;
2028 for (pv1 = p2 + 1; pv1 < ep2; pv1 += 2) {
2029 if (pv1->level == l2) {
2030 opa2[oc2] = pv1 + 1;
2031 used[oc2] = false;
2032 if (++oc2 >= ARR_CNT(opa2)) {
2033 debug_string(1, "Expression too big to compare, because MAX_COMPARE_TERMS exceeded.");
2034 return false;
2035 }
2036 }
2037 }
2038 opa2[oc2] = pv1 + 1;
2039 last_op1 = 0;
2040 first = true;
2041 for (pv1 = p1;;) {
2042 for (len = 1; &pv1[len] < ep1; len += 2)
2043 if (pv1[len].level <= l1)
2044 break;
2045 for (i = 0;; i++) {
2046 if (i >= oc2) {
2047 if ((op1 == TIMES || op1 == DIVIDE) && pv1->level == l1 && pv1->kind == CONSTANT) {
2048 if (fabs(fabs(pv1->token.constant) - 1.0) <= compare_epsilon) {
2049 if (pv1->token.constant < 0.0) {
2050 *diff_signp ^= true;
2051 }
2052 break;
2053 }
2054 }
2055 return false;
2056 }
2057 if (used[i]) {
2058 continue;
2059 }
2060 switch (op1) {
2061 case PLUS:
2062 case MINUS:
2063 break;
2064 case 0:
2065 case TIMES:
2066 case DIVIDE:
2067 if ((last_op1 == DIVIDE) != ((i != 0)
2068 && ((opa2[i] - 1)->token.operatr == DIVIDE)))
2069 continue;
2070 break;
2071 default:
2072 if ((last_op1 == 0) != (i == 0))
2073 return false;
2074 break;
2075 }
2076 if (compare_recurse(pv1, len, (pv1->level <= l1) ? l1 : (l1 + 1),
2077 opa2[i], (opa2[i+1] - opa2[i]) - 1, (opa2[i]->level <= l2) ? l2 : (l2 + 1), &j)) {
2078 switch (op1) {
2079 case 0:
2080 case TIMES:
2081 case DIVIDE:
2082 *diff_signp ^= j;
2083 break;
2084 case PLUS:
2085 case MINUS:
2086 if (last_op1 == MINUS)
2087 j = !j;
2088 if (i != 0 && (opa2[i] - 1)->token.operatr == MINUS)
2089 j = !j;
2090 if (!first) {
2091 if (*diff_signp != j)
2092 continue;
2093 } else {
2094 *diff_signp = j;
2095 first = false;
2096 }
2097 break;
2098 /* case POWER:
2099 Someday fix so that (x-y)^6 compares identical to (y-x)^6.
2100 Also (x-y)^5 should compare identical with different sign to (y-x)^5.
2101 */
2102 default:
2103 if (j)
2104 continue;
2105 break;
2106 }
2107 used[i] = true;
2108 break;
2109 }
2110 }
2111 pv1 += len;
2112 if (pv1 >= ep1)
2113 break;
2114 last_op1 = pv1->token.operatr;
2115 pv1++;
2116 }
2117 for (i = 0; i < oc2; i++) {
2118 if (!used[i]) {
2119 if ((op2 == TIMES || op2 == DIVIDE) && opa2[i]->level == l2 && opa2[i]->kind == CONSTANT) {
2120 if (fabs(fabs(opa2[i]->token.constant) - 1.0) <= compare_epsilon) {
2121 if (opa2[i]->token.constant < 0.0) {
2122 *diff_signp ^= true;
2123 }
2124 continue;
2125 }
2126 }
2127 return false;
2128 }
2129 }
2130 return true;
2131 }
2132
2133 /*
2134 * Take out meaningless "sign" variables and negative constants.
2135 * Simplify imaginary numbers raised to the power of some constants.
2136 *
2137 * Return true if equation side was modified.
2138 */
2139 int
elim_sign(equation,np)2140 elim_sign(equation, np)
2141 token_type *equation;
2142 int *np;
2143 {
2144 int i, j, k;
2145 int level;
2146 int modified = false;
2147 int op;
2148 double d;
2149 double numerator, denominator;
2150
2151 for (i = 1; i < *np; i += 2) {
2152 #if DEBUG
2153 if (equation[i].kind != OPERATOR) {
2154 error_bug("Error in elim_sign().");
2155 }
2156 #endif
2157 level = equation[i].level;
2158 if (equation[i+1].kind == CONSTANT
2159 && equation[i].token.operatr == POWER
2160 && ((equation[i+1].level == level)
2161 || (equation[i+1].level == level + 1))) {
2162 if (equation[i+1].level == level + 1) {
2163 if (i + 3 >= *np || equation[i+2].token.operatr != TIMES)
2164 continue;
2165 for (k = i + 2; k < *np && equation[k].level >= (level + 1); k += 2)
2166 ;
2167 if (k <= i + 2)
2168 continue;
2169 if (!is_integer_expr(&equation[i+3], k - (i + 3)))
2170 continue;
2171 }
2172 f_to_fraction(equation[i+1].token.constant, &numerator, &denominator);
2173 if (always_positive(numerator)) {
2174 if (equation[i-1].level == level
2175 && equation[i-1].kind == VARIABLE
2176 && equation[i-1].token.variable == IMAGINARY) {
2177 equation[i-1].kind = CONSTANT;
2178 equation[i-1].token.constant = -1.0;
2179 equation[i+1].token.constant /= 2.0;
2180 modified = true;
2181 continue;
2182 }
2183 op = 0;
2184 for (j = i - 1; (j >= 0) && equation[j].level >= level; j--) {
2185 if (equation[j].level <= level + 1) {
2186 if (equation[j].kind == OPERATOR) {
2187 op = equation[j].token.operatr;
2188 break;
2189 }
2190 }
2191 }
2192 switch (op) {
2193 case 0:
2194 case TIMES:
2195 case DIVIDE:
2196 for (j = i - 1; (j >= 0) && equation[j].level >= level; j--) {
2197 if (equation[j].level <= level + 1) {
2198 if (equation[j].kind == VARIABLE
2199 && (equation[j].token.variable & VAR_MASK) == SIGN) {
2200 equation[j].kind = CONSTANT;
2201 equation[j].token.constant = 1.0;
2202 modified = true;
2203 } else if (equation[j].kind == CONSTANT
2204 && equation[j].token.constant < 0.0) {
2205 equation[j].token.constant = -equation[j].token.constant;
2206 modified = true;
2207 }
2208 }
2209 }
2210 }
2211 } else {
2212 if (equation[i-1].level == level && equation[i-1].kind == VARIABLE) {
2213 if (equation[i-1].token.variable == IMAGINARY && equation[i+1].level == level) {
2214 d = fmod(equation[i+1].token.constant, 4.0);
2215 if (d == 1.0) {
2216 equation[i].token.operatr = TIMES;
2217 equation[i+1].token.constant = 1.0;
2218 modified = true;
2219 } else if (d == 3.0) {
2220 equation[i].token.operatr = TIMES;
2221 equation[i+1].token.constant = -1.0;
2222 modified = true;
2223 }
2224 } else if ((equation[i-1].token.variable & VAR_MASK) == SIGN) {
2225 if (fmod(denominator, 2.0) == 1.0) {
2226 /* odd root, make root (denominator) 1 */
2227 numerator = fmod(numerator, 2.0);
2228 /* all sign^2 now translate to 1 */
2229 if (numerator != equation[i+1].token.constant) {
2230 equation[i+1].token.constant = numerator;
2231 modified = true;
2232 }
2233 }
2234 }
2235 }
2236 }
2237 }
2238 }
2239 return modified;
2240 }
2241
2242 /*
2243 * Do complex number division and move the imaginary number from the denominator to the numerator.
2244 * This is done by multiplying the numerator and denominator by
2245 * the conjugate of the complex number in the denominator.
2246 * This often complicates expressions, so use with care.
2247 *
2248 * Includes simple division (1/i# -> -1*i#) conversions now.
2249 *
2250 * Return true if equation side was modified.
2251 */
2252 int
div_imaginary(equation,np)2253 div_imaginary(equation, np)
2254 token_type *equation;
2255 int *np;
2256 {
2257 int i, j, k;
2258 int n;
2259 int level;
2260 int ilevel;
2261 int op;
2262 int iloc, biloc, eiloc;
2263 int eloc;
2264 int modified = false;
2265
2266 for (i = 1; i < *np; i += 2) {
2267 #if DEBUG
2268 if (equation[i].kind != OPERATOR) {
2269 error_bug("Error in div_imaginary().");
2270 }
2271 #endif
2272 if (equation[i].token.operatr == DIVIDE) {
2273 level = equation[i].level;
2274 #if 1 /* This necessary code makes complex number calculations give different (maybe wrong) results. */
2275 /* It converts x/i to x*-1*i. Might be better to do this only when symb_flag is true. */
2276 if (equation[i+1].level == level
2277 && equation[i+1].kind == VARIABLE
2278 && equation[i+1].token.variable == IMAGINARY) {
2279 if (*np + 2 > n_tokens) {
2280 error_huge();
2281 }
2282 blt(&equation[i+2], &equation[i], (*np - i) * sizeof(token_type));
2283 *np += 2;
2284 equation[i].level = level;
2285 equation[i].kind = OPERATOR;
2286 equation[i].token.operatr = TIMES;
2287 i++;
2288 equation[i].level = level;
2289 equation[i].kind = CONSTANT;
2290 equation[i].token.constant = -1.0;
2291 i++;
2292 equation[i].level = level;
2293 equation[i].kind = OPERATOR;
2294 equation[i].token.operatr = TIMES;
2295 modified = true;
2296 continue;
2297 }
2298 #endif
2299 op = 0;
2300 iloc = biloc = eiloc = -1;
2301 k = i;
2302 for (j = i + 1; j < *np && equation[j].level > level; j++) {
2303 if (equation[j].kind == OPERATOR && equation[j].level == level + 1) {
2304 op = equation[j].token.operatr;
2305 k = j;
2306 if (iloc >= 0 && eiloc < 0) {
2307 eiloc = j;
2308 }
2309 } else if (equation[j].kind == VARIABLE && equation[j].token.variable == IMAGINARY) {
2310 if (iloc >= 0) {
2311 op = 0;
2312 break;
2313 }
2314 iloc = j;
2315 biloc = k + 1;
2316 }
2317 }
2318 eloc = j;
2319 if (iloc >= 0 && eiloc < 0) {
2320 eiloc = j;
2321 }
2322 if (iloc < 0 || (op != PLUS && op != MINUS))
2323 continue;
2324 ilevel = equation[iloc].level;
2325 if (ilevel != level + 1) {
2326 if (ilevel != level + 2) {
2327 continue;
2328 }
2329 if (iloc > biloc && equation[iloc-1].token.operatr != TIMES)
2330 continue;
2331 if (iloc + 1 < eiloc) {
2332 switch (equation[iloc+1].token.operatr) {
2333 case TIMES:
2334 case DIVIDE:
2335 break;
2336 default:
2337 continue;
2338 }
2339 }
2340 }
2341 if ((eloc - (i + 1)) + 5 + (eiloc - biloc) + *np + 2 > n_tokens)
2342 error_huge();
2343 n = eloc - (i + 1);
2344 blt(scratch, &equation[i+1], n * sizeof(token_type));
2345 scratch[iloc-(i+1)].kind = CONSTANT;
2346 scratch[iloc-(i+1)].token.constant = 0.0;
2347 for (j = 0; j < n; j++)
2348 scratch[j].level += 2;
2349 scratch[n].level = level + 2;
2350 scratch[n].kind = OPERATOR;
2351 scratch[n].token.operatr = POWER;
2352 n++;
2353 scratch[n].level = level + 2;
2354 scratch[n].kind = CONSTANT;
2355 scratch[n].token.constant = 2.0;
2356 n++;
2357 scratch[n].level = level + 1;
2358 scratch[n].kind = OPERATOR;
2359 scratch[n].token.operatr = PLUS;
2360 n++;
2361 blt(&scratch[n], &equation[biloc], (eiloc - biloc) * sizeof(token_type));
2362 j = n;
2363 n += (eiloc - biloc);
2364 for (k = j; k < n; k++)
2365 scratch[k].level += 2;
2366 scratch[n].level = level + 2;
2367 scratch[n].kind = OPERATOR;
2368 scratch[n].token.operatr = POWER;
2369 n++;
2370 scratch[n].level = level + 2;
2371 scratch[n].kind = CONSTANT;
2372 scratch[n].token.constant = 2.0;
2373 n++;
2374 scratch[j+(iloc-biloc)].kind = CONSTANT;
2375 scratch[j+(iloc-biloc)].token.constant = 1.0;
2376 blt(&equation[iloc+2], &equation[iloc], (*np - iloc) * sizeof(token_type));
2377 *np += 2;
2378 ilevel++;
2379 equation[iloc].level = ilevel;
2380 equation[iloc].kind = CONSTANT;
2381 equation[iloc].token.constant = -1.0;
2382 iloc++;
2383 equation[iloc].level = ilevel;
2384 equation[iloc].kind = OPERATOR;
2385 equation[iloc].token.operatr = TIMES;
2386 iloc++;
2387 equation[iloc].level = ilevel;
2388 blt(&equation[i+1+n], &equation[i], (*np - i) * sizeof(token_type));
2389 *np += n + 1;
2390 blt(&equation[i+1], scratch, n * sizeof(token_type));
2391 i += n + 1;
2392 equation[i].token.operatr = TIMES;
2393 modified = true;
2394 }
2395 }
2396 return modified;
2397 }
2398
2399 /*
2400 * Convert 1/m*(-1*b+y) to (y-b)/m in equation side.
2401 *
2402 * Returns true if equation side was modified.
2403 * "elim_k()" should be called after this if it returns true.
2404 *
2405 * This routine also checks the validity of the equation side.
2406 */
2407 int
reorder(equation,np)2408 reorder(equation, np)
2409 token_type *equation;
2410 int *np;
2411 {
2412 return(order_recurse(equation, np, 0, 1));
2413 }
2414
2415 static void
swap(equation,np,level,i1,i2)2416 swap(equation, np, level, i1, i2)
2417 token_type *equation;
2418 int *np;
2419 int level;
2420 int i1, i2;
2421 {
2422 int e1, e2;
2423 int n1, n2;
2424
2425 for (e1 = i1 + 1;; e1 += 2) {
2426 if (e1 >= *np || equation[e1].level <= level)
2427 break;
2428 }
2429 for (e2 = i2 + 1;; e2 += 2) {
2430 if (e2 >= *np || equation[e2].level <= level)
2431 break;
2432 }
2433 n1 = e1 - i1;
2434 n2 = e2 - i2;
2435 blt(scratch, &equation[i1], (e2 - i1) * sizeof(token_type));
2436 if ((i1 + n2) != e1) {
2437 blt(&equation[i1+n2], &equation[e1], (i2 - e1) * sizeof(token_type));
2438 }
2439 blt(&equation[i1], &scratch[i2-i1], n2 * sizeof(token_type));
2440 blt(&equation[e2-n1], scratch, n1 * sizeof(token_type));
2441 }
2442
2443 /*
2444 * Recursively check and possibly reorder a "level" sub-expression
2445 * which starts at "loc" in an equation side.
2446 */
2447 static int
order_recurse(equation,np,loc,level)2448 order_recurse(equation, np, loc, level)
2449 token_type *equation;
2450 int *np, loc, level;
2451 {
2452 int i, j, k, n;
2453 int op = 0;
2454 int modified = false;
2455
2456 /* check passed sub-expression for validity: */
2457 if ((loc & 1) != 0) {
2458 goto corrupt;
2459 }
2460 for (i = loc; i < *np;) {
2461 if (equation[i].level < level) {
2462 if (equation[i].kind != OPERATOR || equation[i].level < 1) {
2463 goto corrupt;
2464 }
2465 break;
2466 }
2467 if (equation[i].level > level) {
2468 modified |= order_recurse(equation, np, i, level + 1);
2469 i++;
2470 for (; i < *np && equation[i].level > level; i++)
2471 ;
2472 continue;
2473 } else {
2474 if (equation[i].kind == OPERATOR) {
2475 if ((i & 1) == 0 || equation[i].token.operatr == 0) {
2476 goto corrupt;
2477 }
2478 if (op) {
2479 switch (equation[i].token.operatr) {
2480 case PLUS:
2481 case MINUS:
2482 if (op == PLUS || op == MINUS)
2483 break;
2484 goto corrupt;
2485 case TIMES:
2486 case DIVIDE:
2487 if (op == TIMES || op == DIVIDE)
2488 break;
2489 default:
2490 goto corrupt;
2491 }
2492 } else {
2493 op = equation[i].token.operatr;
2494 }
2495 } else {
2496 if ((i & 1) != 0)
2497 goto corrupt;
2498 }
2499 }
2500 i++;
2501 }
2502 if ((i & 1) == 0) { /* "i" should be at the end of the current sub-expression */
2503 goto corrupt;
2504 }
2505 /* reorder sub-expression if needed: */
2506 switch (op) {
2507 case PLUS:
2508 case MINUS:
2509 if (equation[loc].kind == CONSTANT && equation[loc].token.constant < 0.0) {
2510 if (equation[loc].level == level || (equation[loc+1].level == level + 1
2511 && (equation[loc+1].token.operatr == TIMES || equation[loc+1].token.operatr == DIVIDE))) {
2512 for (j = loc + 1; j < i; j += 2) {
2513 if (equation[j].level == level && equation[j].token.operatr == PLUS) {
2514 swap(equation, np, level, loc, j + 1);
2515 modified = true;
2516 break;
2517 }
2518 }
2519 }
2520 }
2521 break;
2522 case TIMES:
2523 case DIVIDE:
2524 for (j = loc + 1;; j += 2) {
2525 if (j >= i)
2526 return modified;
2527 if (equation[j].level == level && equation[j].token.operatr == DIVIDE)
2528 break;
2529 }
2530 for (k = j + 2; k < i;) {
2531 if (equation[k].level == level && equation[k].token.operatr == TIMES) {
2532 for (n = k + 2; n < i && equation[n].level > level; n += 2)
2533 ;
2534 n -= k;
2535 blt(scratch, &equation[k], n * sizeof(token_type));
2536 blt(&equation[j+n], &equation[j], (k - j) * sizeof(token_type));
2537 blt(&equation[j], scratch, n * sizeof(token_type));
2538 j += n;
2539 k += n;
2540 modified = true;
2541 continue;
2542 }
2543 k += 2;
2544 }
2545 break;
2546 }
2547 return modified;
2548 corrupt:
2549 error_bug("Internal representation of expression is corrupt!");
2550 return modified; /* not reached */
2551 }
2552
2553 /*
2554 * Try to rationalize the denominator of algebraic fractions.
2555 * Only works with a square root in the denominator
2556 * and sometimes helps with simplification.
2557 *
2558 * Returns true if something was done.
2559 * Unfactoring needs to be done immediately, if this returns true.
2560 */
2561 int
rationalize(equation,np)2562 rationalize(equation, np)
2563 token_type *equation;
2564 int *np;
2565 {
2566 int i, j, k;
2567 int i1, k1;
2568 int div_level;
2569 int end_loc, neg_one_loc = -1;
2570 int flag;
2571 int modified = false;
2572 int count;
2573
2574 for (i = 1;; i += 2) {
2575 be_thorough:
2576 if (i >= *np)
2577 break;
2578 #if DEBUG
2579 if (equation[i].kind != OPERATOR) {
2580 error_bug("Bug in rationalize().");
2581 }
2582 #endif
2583 if (equation[i].token.operatr == DIVIDE) {
2584 div_level = equation[i].level;
2585 count = 0;
2586 j = -1;
2587 for (end_loc = i + 2; end_loc < *np && equation[end_loc].level > div_level; end_loc += 2) {
2588 if (equation[end_loc].level == (div_level + 1)) {
2589 count++;
2590 if (j < 0) {
2591 j = end_loc;
2592 }
2593 }
2594 }
2595 if (j < 0)
2596 continue;
2597 switch (equation[j].token.operatr) {
2598 case PLUS:
2599 case MINUS:
2600 break;
2601 default:
2602 continue;
2603 }
2604 i1 = i;
2605 do_again:
2606 flag = false;
2607 for (k = j - 2; k > i1; k -= 2) {
2608 if (equation[k].level == (div_level + 2)) {
2609 switch (equation[k].token.operatr) {
2610 case TIMES:
2611 case DIVIDE:
2612 flag = 1;
2613 break;
2614 case POWER:
2615 flag = 2;
2616 break;
2617 }
2618 break;
2619 }
2620 }
2621 if (flag) {
2622 for (k = j - 2; k > i1; k -= 2) {
2623 if ((equation[k].level == (div_level + 2)
2624 || (equation[k].level == (div_level + 3) && flag == 1))
2625 && equation[k].token.operatr == POWER
2626 && equation[k].level == equation[k+1].level
2627 && equation[k+1].kind == CONSTANT
2628 && fmod(equation[k+1].token.constant, 1.0) == 0.5) {
2629 for (k1 = i + 2; k1 < end_loc; k1 += 2) {
2630 if (equation[k1].token.operatr == POWER
2631 && equation[k1].level == equation[k1+1].level
2632 && equation[k1+1].kind == CONSTANT
2633 && fmod(equation[k1+1].token.constant, 1.0) == 0.5) {
2634 /* make sure we will actually be eliminating ^.5 from the denominator: */
2635 if (k != k1 && !(equation[k1].level == (div_level + 2) && count == 1)) {
2636 i += 2;
2637 goto be_thorough;
2638 }
2639 /* avoid rationalizing absolute values: */
2640 if (equation[k1-1].level == (equation[k1].level + 1)
2641 && equation[k1-2].level == equation[k1-1].level
2642 && equation[k1-1].kind == CONSTANT
2643 && equation[k1-2].token.operatr == POWER) {
2644 i += 2;
2645 goto be_thorough;
2646 }
2647 }
2648 }
2649 neg_one_loc = i1 + 1;
2650 k = i1 - i;
2651 blt(scratch, &equation[i+1], k * sizeof(token_type));
2652 scratch[k].level = div_level + 2;
2653 scratch[k].kind = CONSTANT;
2654 scratch[k].token.constant = -1.0;
2655 k++;
2656 scratch[k].level = div_level + 2;
2657 scratch[k].kind = OPERATOR;
2658 scratch[k].token.operatr = TIMES;
2659 k++;
2660 blt(&scratch[k], &equation[neg_one_loc], (end_loc - neg_one_loc) * sizeof(token_type));
2661 for (k1 = 0; k1 < (j - neg_one_loc); k1++, k++)
2662 scratch[k].level++;
2663 k = end_loc - (i + 1) + 2;
2664 if (*np + 2 * (k + 1) > n_tokens) {
2665 error_huge();
2666 }
2667 blt(&equation[end_loc+(2*(k+1))], &equation[end_loc], (*np - end_loc) * sizeof(token_type));
2668 *np += 2 * (k + 1);
2669 k1 = end_loc;
2670 equation[k1].level = div_level;
2671 equation[k1].kind = OPERATOR;
2672 equation[k1].token.operatr = TIMES;
2673 k1++;
2674 blt(&equation[k1], scratch, k * sizeof(token_type));
2675 k1 += k;
2676 equation[k1].level = div_level;
2677 equation[k1].kind = OPERATOR;
2678 equation[k1].token.operatr = DIVIDE;
2679 k1++;
2680 blt(&equation[k1], scratch, k * sizeof(token_type));
2681 k1 += k;
2682 debug_string(1, "Square roots in denominator rationalized:");
2683 side_debug(1, &equation[i+1], k1 - (i + 1));
2684 i = k1;
2685 modified = true;
2686 goto be_thorough;
2687 }
2688 }
2689 }
2690 if (j >= end_loc)
2691 continue;
2692 i1 = j;
2693 for (j += 2; j < end_loc; j += 2) {
2694 if (equation[j].level == (div_level + 1)) {
2695 break;
2696 }
2697 }
2698 goto do_again;
2699 }
2700 }
2701 return modified;
2702 }
2703