1 /*
2 * Mathomatic unfactorizing (expanding) 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 static int unf_sub(token_type *equation, int *np, int b1, int loc, int e1, int level, int ii);
28
29 /*
30 * Unfactor times and divide only (products of sums) and simplify.
31 *
32 * Return true if equation side was unfactored.
33 */
34 int
uf_tsimp(equation,np)35 uf_tsimp(equation, np)
36 token_type *equation;
37 int *np;
38 {
39 int rv;
40
41 rv = uf_times(equation, np);
42 simp_loop(equation, np);
43 while (uf_times(equation, np)) {
44 rv = true;
45 simp_loop(equation, np);
46 }
47 return rv;
48 }
49
50 /*
51 * Unfactor power only.
52 * (a * b)^c -> a^c * b^c
53 *
54 * Return true if equation side was unfactored.
55 */
56 int
uf_power(equation,np)57 uf_power(equation, np)
58 token_type *equation;
59 int *np;
60 {
61 int count = -1;
62
63 do {
64 organize(equation, np);
65 if (++count > 0)
66 break;
67 } while (sub_ufactor(equation, np, 2));
68 return count;
69 }
70
71 /*
72 * Unfactor power only.
73 * a^(b + c) -> a^b * a^c
74 *
75 * Return true if equation side was unfactored.
76 */
77 int
uf_pplus(equation,np)78 uf_pplus(equation, np)
79 token_type *equation;
80 int *np;
81 {
82 int count = -1;
83
84 do {
85 organize(equation, np);
86 if (++count > 0)
87 break;
88 } while (sub_ufactor(equation, np, 4));
89 return count;
90 }
91
92 /*
93 * Unfactor all power operators.
94 * Same as a call to uf_pplus() and uf_power(), only faster.
95 */
96 void
uf_allpower(equation,np)97 uf_allpower(equation, np)
98 token_type *equation;
99 int *np;
100 {
101 do {
102 organize(equation, np);
103 } while (sub_ufactor(equation, np, 0));
104 }
105
106 /*
107 * Unfactor power only if it will help with expansion and exponent is <= 100.
108 * (a + 1)^2 -> (a + 1)*(a + 1)
109 * Afterwards we also sneak in simplifying division by irrational constants,
110 * because it works best here.
111 *
112 * uf_times() is usually called after this to complete the expansion.
113 */
114 void
uf_repeat(equation,np)115 uf_repeat(equation, np)
116 token_type *equation;
117 int *np;
118 {
119 int count = -1;
120
121 do {
122 organize(equation, np);
123 if (++count > 0)
124 break;
125 } while (sub_ufactor(equation, np, 6));
126 patch_root_div(equation, np);
127 }
128
129 /*
130 * Unfactor power only.
131 * a^2 -> a*a
132 * Useful for removing all integer powers.
133 */
134 void
uf_repeat_always(equation,np)135 uf_repeat_always(equation, np)
136 token_type *equation;
137 int *np;
138 {
139 int count = -1;
140
141 do {
142 organize(equation, np);
143 if (++count > 0)
144 break;
145 } while (sub_ufactor(equation, np, 8));
146 }
147
148 /*
149 * Totally unfactor equation side and simplify.
150 */
151 void
uf_simp(equation,np)152 uf_simp(equation, np)
153 token_type *equation; /* pointer to beginning of equation side */
154 int *np; /* pointer to length of equation side */
155 {
156 uf_tsimp(equation, np);
157 uf_power(equation, np);
158 uf_repeat(equation, np);
159 uf_tsimp(equation, np);
160 }
161
162 /*
163 * Unfactor equation side and simplify.
164 * Don't call uf_repeat().
165 */
166 void
uf_simp_no_repeat(equation,np)167 uf_simp_no_repeat(equation, np)
168 token_type *equation;
169 int *np;
170 {
171 uf_power(equation, np);
172 uf_tsimp(equation, np);
173 }
174
175 /*
176 * Totally unfactor equation side with no simplification.
177 */
178 int
ufactor(equation,np)179 ufactor(equation, np)
180 token_type *equation;
181 int *np;
182 {
183 int rv;
184
185 uf_repeat(equation, np);
186 rv = uf_times(equation, np);
187 uf_allpower(equation, np);
188 return rv;
189 }
190
191 /*
192 * Increase the level of numerators by 2, so that the divide operator is not unfactored.
193 */
194 static void
no_divide(equation,np)195 no_divide(equation, np)
196 token_type *equation;
197 int *np;
198 {
199 int i, j;
200 int level;
201
202 for (i = 1; i < *np; i += 2) {
203 if (equation[i].token.operatr == DIVIDE) {
204 level = equation[i].level;
205 for (j = i - 1; j >= 0; j--) {
206 if (equation[j].level < level)
207 break;
208 equation[j].level += 2;
209 }
210 }
211 }
212 }
213
214 /*
215 * Unfactor times and divide only (products of sums like (a+b)*(c+d)).
216 * (a + b)*c -> a*c + b*c
217 * If partial_flag is true, don't expand (a+b)/(c+d) nor (a+b)/c.
218 *
219 * Return true if equation side was unfactored.
220 */
221 int
uf_times(equation,np)222 uf_times(equation, np)
223 token_type *equation;
224 int *np;
225 {
226 int i;
227 int rv = false;
228
229 do {
230 organize(equation, np);
231 if (reorder(equation, np)) {
232 organize(equation, np);
233 }
234 group_proc(equation, np);
235 if (partial_flag) {
236 no_divide(equation, np);
237 }
238 i = sub_ufactor(equation, np, 1);
239 rv |= i;
240 } while (i);
241 organize(equation, np);
242 return rv;
243 }
244
245 /*
246 * General equation side algebraic expansion routine.
247 * Expands everything.
248 * The type of expansion to be done is indicated by the code in "ii".
249 *
250 * Return true if equation side was modified.
251 */
252 int
sub_ufactor(equation,np,ii)253 sub_ufactor(equation, np, ii)
254 token_type *equation;
255 int *np;
256 int ii;
257 {
258 int modified = false;
259 int i;
260 int b1, e1;
261 int level;
262
263 for (i = 1; i < *np; i += 2) {
264 switch (equation[i].token.operatr) {
265 case TIMES:
266 case DIVIDE:
267 if (ii == 1)
268 break;
269 else
270 continue;
271 case POWER:
272 if ((ii & 1) == 0) /* even number codes only */
273 break;
274 default:
275 continue;
276 }
277 level = equation[i].level;
278 for (b1 = i - 2; b1 >= 0; b1 -= 2)
279 if (equation[b1].level < level)
280 break;
281 b1++;
282 for (e1 = i + 2; e1 < *np; e1 += 2) {
283 if (equation[e1].level < level)
284 break;
285 }
286 if (unf_sub(equation, np, b1, i, e1, level, ii)) {
287 modified = true;
288 i = b1 - 1;
289 continue;
290 }
291 }
292 return modified;
293 }
294
295 static int
unf_sub(equation,np,b1,loc,e1,level,ii)296 unf_sub(equation, np, b1, loc, e1, level, ii)
297 token_type *equation;
298 int *np;
299 int b1, loc, e1, level;
300 int ii;
301 {
302 int i, j, k;
303 int b2, eb1, be1;
304 int len;
305 double d1, d2;
306
307 switch (equation[loc].token.operatr) {
308 case TIMES:
309 case DIVIDE:
310 if (ii != 1)
311 break;
312 for (i = b1 + 1; i < e1; i += 2) {
313 if (equation[i].level == level + 1) {
314 switch (equation[i].token.operatr) {
315 case PLUS:
316 case MINUS:
317 break;
318 default:
319 continue;
320 }
321 for (b2 = i - 2; b2 >= b1; b2 -= 2)
322 if (equation[b2].level <= level)
323 break;
324 b2++;
325 eb1 = b2;
326 for (be1 = i + 2; be1 < e1; be1 += 2)
327 if (equation[be1].level <= level)
328 break;
329 if (eb1 > b1 && equation[eb1-1].token.operatr == DIVIDE) {
330 i = be1 - 2;
331 continue;
332 }
333 len = 0;
334 u_again:
335 if (len + (eb1 - b1) + (i - b2) + (e1 - be1) + 1 > n_tokens) {
336 error_huge();
337 }
338 blt(&scratch[len], &equation[b1], (eb1 - b1) * sizeof(token_type));
339 j = len;
340 len += (eb1 - b1);
341 for (; j < len; j++)
342 scratch[j].level++;
343 blt(&scratch[len], &equation[b2], (i - b2) * sizeof(token_type));
344 len += (i - b2);
345 blt(&scratch[len], &equation[be1], (e1 - be1) * sizeof(token_type));
346 j = len;
347 len += (e1 - be1);
348 for (; j < len; j++)
349 scratch[j].level++;
350 if (i < be1) {
351 scratch[len] = equation[i];
352 scratch[len].level--;
353 len++;
354 b2 = i + 1;
355 for (i += 2; i < be1; i += 2) {
356 if (equation[i].level == (level + 1))
357 break;
358 }
359 goto u_again;
360 } else {
361 if (*np - (e1 - b1) + len > n_tokens) {
362 error_huge();
363 }
364 blt(&equation[b1+len], &equation[e1], (*np - e1) * sizeof(token_type));
365 *np += len - (e1 - b1);
366 blt(&equation[b1], scratch, len * sizeof(token_type));
367 return true;
368 }
369 }
370 }
371 break;
372 case POWER:
373 #if !ALWAYS_UNFACTOR_POWER
374 /* avoid absolute values: */
375 if ((loc + 3) < *np && equation[loc+1].level == level && equation[loc+1].kind == CONSTANT
376 && equation[loc+2].level == (level - 1) && equation[loc+2].token.operatr == POWER
377 && equation[loc+3].kind == CONSTANT && ((equation[loc+3].level == (level - 1))
378 || ((loc + 5) < *np && equation[loc+3].level == level
379 && equation[loc+4].level == level && equation[loc+4].token.operatr == DIVIDE
380 && equation[loc+5].level == level && equation[loc+5].kind == CONSTANT
381 && ((loc + 6) >= *np || equation[loc+6].level < level)))) {
382 return false;
383 }
384 #endif
385 if (ii == 2 || ii == 0) {
386 for (i = b1 + 1; i < loc; i += 2) {
387 if (equation[i].level != level + 1)
388 continue;
389 switch (equation[i].token.operatr) {
390 case TIMES:
391 case DIVIDE:
392 break;
393 default:
394 goto do_pplus;
395 }
396 b2 = b1;
397 len = 0;
398 u1_again:
399 if (len + (i - b2) + (e1 - loc) + 1 > n_tokens) {
400 error_huge();
401 }
402 blt(&scratch[len], &equation[b2], (i - b2) * sizeof(token_type));
403 len += (i - b2);
404 blt(&scratch[len], &equation[loc], (e1 - loc) * sizeof(token_type));
405 j = len;
406 len += (e1 - loc);
407 for (; j < len; j++)
408 scratch[j].level++;
409 if (i < loc) {
410 scratch[len] = equation[i];
411 scratch[len].level--;
412 len++;
413 b2 = i + 1;
414 for (i += 2; i < loc; i += 2) {
415 if (equation[i].level == (level + 1))
416 break;
417 }
418 goto u1_again;
419 } else {
420 if (*np - (e1 - b1) + len > n_tokens) {
421 error_huge();
422 }
423 blt(&equation[b1+len], &equation[e1], (*np - e1) * sizeof(token_type));
424 *np += len - (e1 - b1);
425 blt(&equation[b1], scratch, len * sizeof(token_type));
426 return true;
427 }
428 }
429 }
430 do_pplus:
431 if (ii == 4 || ii == 0) {
432 for (i = loc + 2; i < e1; i += 2) {
433 if (equation[i].level != level + 1)
434 continue;
435 switch (equation[i].token.operatr) {
436 case PLUS:
437 case MINUS:
438 break;
439 default:
440 goto do_repeat;
441 }
442 b2 = loc + 1;
443 len = 0;
444 u2_again:
445 if (len + (loc - b1) + (i - b2) + 2 > n_tokens) {
446 error_huge();
447 }
448 j = len;
449 blt(&scratch[len], &equation[b1], (loc + 1 - b1) * sizeof(token_type));
450 len += (loc + 1 - b1);
451 for (; j < len; j++)
452 scratch[j].level++;
453 blt(&scratch[len], &equation[b2], (i - b2) * sizeof(token_type));
454 len += (i - b2);
455 if (i < e1) {
456 scratch[len].level = level;
457 scratch[len].kind = OPERATOR;
458 if (equation[i].token.operatr == PLUS)
459 scratch[len].token.operatr = TIMES;
460 else
461 scratch[len].token.operatr = DIVIDE;
462 len++;
463 b2 = i + 1;
464 for (i += 2; i < e1; i += 2) {
465 if (equation[i].level == (level + 1))
466 break;
467 }
468 goto u2_again;
469 } else {
470 if (*np - (e1 - b1) + len > n_tokens) {
471 error_huge();
472 }
473 blt(&equation[b1+len], &equation[e1], (*np - e1) * sizeof(token_type));
474 *np += len - (e1 - b1);
475 blt(&equation[b1], scratch, len * sizeof(token_type));
476 return true;
477 }
478 }
479 }
480 do_repeat:
481 /* unfactor power: (a + 1)^2 -> (a + 1)*(a + 1) */
482 if (ii != 6 && ii != 8)
483 break;
484 i = loc;
485 if (equation[i+1].level != level || equation[i+1].kind != CONSTANT)
486 break;
487 d1 = equation[i+1].token.constant;
488 if (!isfinite(d1) || d1 <= 1.0)
489 break;
490 if (ii != 8) { /* if true, only do useful expansions */
491 if (d1 > 100.0) /* ignore exponents over 100 */
492 break;
493 if ((i - b1) == 1 && equation[b1].kind != CONSTANT)
494 break;
495 if ((i - b1) > 1 && d1 > 2.0 && fmod(d1, 1.0) != 0.0)
496 break;
497 }
498 d1 = ceil(d1) - 1.0;
499 d2 = d1 * ((i - b1) + 1.0);
500 if ((*np + d2) > (n_tokens - 10)) {
501 break; /* too big to expand, do nothing */
502 }
503 j = d1;
504 k = d2;
505 blt(&equation[e1+k], &equation[e1], (*np - e1) * sizeof(token_type));
506 *np += k;
507 equation[i+1].token.constant -= d1;
508 k = e1;
509 while (j-- > 0) {
510 equation[k].level = level;
511 equation[k].kind = OPERATOR;
512 equation[k].token.operatr = TIMES;
513 blt(&equation[k+1], &equation[b1], (i - b1) * sizeof(token_type));
514 k += (i - b1) + 1;
515 }
516 if (equation[i+1].token.constant == 1.0) {
517 blt(&equation[i], &equation[e1], (*np - e1) * sizeof(token_type));
518 *np -= (e1 - i);
519 } else {
520 for (j = b1; j < e1; j++)
521 equation[j].level++;
522 }
523 return true;
524 }
525 return false;
526 }
527
528 static int
usp_sub(equation,np,i)529 usp_sub(equation, np, i)
530 token_type *equation;
531 int *np, i;
532 {
533 int level;
534 int j;
535
536 level = equation[i].level;
537 for (j = i - 2;; j -= 2) {
538 if (j < 0)
539 return false;
540 if (equation[j].level < level) {
541 if (equation[j].level == (level - 1) && equation[j].token.operatr == DIVIDE) {
542 break;
543 } else
544 return false;
545 }
546 }
547 if ((*np + 2) > n_tokens) {
548 error_huge();
549 }
550 equation[j].token.operatr = TIMES;
551 for (j = i + 1;; j++) {
552 if (j >= *np || equation[j].level < level)
553 break;
554 equation[j].level++;
555 }
556 i++;
557 blt(&equation[i+2], &equation[i], (*np - i) * sizeof(token_type));
558 *np += 2;
559 equation[i].level = level + 1;
560 equation[i].kind = CONSTANT;
561 equation[i].token.constant = -1.0;
562 i++;
563 equation[i].level = level + 1;
564 equation[i].kind = OPERATOR;
565 equation[i].token.operatr = TIMES;
566 return true;
567 }
568
569 /*
570 * Convert a/(x^y) to a*x^(-1*y).
571 * If y is a constant, don't do.
572 *
573 * Return true if equation side is modified.
574 */
575 int
unsimp_power(equation,np)576 unsimp_power(equation, np)
577 token_type *equation;
578 int *np;
579 {
580 int modified = false;
581 int i;
582
583 for (i = 1; i < *np; i += 2) {
584 if (equation[i].token.operatr == POWER) {
585 if (equation[i].level == equation[i+1].level
586 && equation[i+1].kind == CONSTANT)
587 continue;
588 modified |= usp_sub(equation, np, i);
589 }
590 }
591 return modified;
592 }
593
594 #if 0 /* the following is not currently used */
595 /*
596 * Convert a/(x^y) to a*(1/x)^y.
597 * If y is a constant, don't do.
598 *
599 * Return true if equation side is modified.
600 */
601 int
602 unsimp2_power(equation, np)
603 token_type *equation;
604 int *np;
605 {
606 int modified = false;
607 int i;
608
609 for (i = 1; i < *np; i += 2) {
610 if (equation[i].token.operatr == POWER) {
611 modified |= usp2_sub(equation, np, i);
612 }
613 }
614 return modified;
615 }
616
617 int
618 usp2_sub(equation, np, i)
619 token_type *equation;
620 int *np, i;
621 {
622 int level;
623 int j, k;
624
625 level = equation[i].level;
626 if (equation[i+1].level == level && equation[i+1].kind == CONSTANT) {
627 return false;
628 }
629 for (j = i - 2;; j -= 2) {
630 if (j < 0)
631 return false;
632 if (equation[j].level < level) {
633 if (equation[j].level == (level - 1) && equation[j].token.operatr == DIVIDE) {
634 break;
635 } else
636 return false;
637 }
638 }
639 if ((*np + 2) > n_tokens) {
640 error_huge();
641 }
642 equation[j].token.operatr = TIMES;
643 j++;
644 for (k = j; k < i; k++) {
645 equation[k].level++;
646 }
647 blt(&equation[j+2], &equation[j], (*np - j) * sizeof(token_type));
648 *np += 2;
649 equation[j].level = level + 1;
650 equation[j].kind = CONSTANT;
651 equation[j].token.constant = 1.0;
652 j++;
653 equation[j].level = level + 1;
654 equation[j].kind = OPERATOR;
655 equation[j].token.operatr = DIVIDE;
656 return true;
657 }
658 #endif
659
660 /*
661 * Convert anything times a negative constant to a positive constant
662 * divided by -1.
663 * When uf_times() is done after this, the additive denominators will be
664 * attempted to be negated, possibly getting rid of unneeded times -1.
665 */
666 void
uf_neg_help(equation,np)667 uf_neg_help(equation, np)
668 token_type *equation;
669 int *np;
670 {
671 int i;
672 int level;
673
674 for (i = 0; i < *np - 1; i += 2) {
675 if (equation[i].kind == CONSTANT && equation[i].token.constant < 0.0) {
676 level = equation[i].level;
677 if (equation[i+1].level == level) {
678 switch (equation[i+1].token.operatr) {
679 case TIMES:
680 case DIVIDE:
681 if ((*np + 2) > n_tokens) {
682 error_huge();
683 }
684 blt(&equation[i+3], &equation[i+1], (*np - (i + 1)) * sizeof(token_type));
685 *np += 2;
686 equation[i].token.constant = -equation[i].token.constant;
687 i++;
688 equation[i].level = level;
689 equation[i].kind = OPERATOR;
690 equation[i].token.operatr = DIVIDE;
691 i++;
692 equation[i].level = level;
693 equation[i].kind = CONSTANT;
694 equation[i].token.constant = -1.0;
695 break;
696 }
697 }
698 }
699 }
700 }
701
702 /*
703 * Simplify division by irrational constants (roots like 2^.5) by
704 * converting 1/(k1^k2) to 1/(k1^(k2-1))/k1 if k1 is integer,
705 * otherwise convert 1/(k1^k2) to 1*((1/k1)^k2).
706 *
707 * Return true if equation side was modified.
708 */
709 int
patch_root_div(equation,np)710 patch_root_div(equation, np)
711 token_type *equation;
712 int *np;
713 {
714 int i;
715 int level;
716 int modified = false;
717
718 for (i = 1; i < *np - 2; i += 2) {
719 if (equation[i].token.operatr == DIVIDE) {
720 level = equation[i].level + 1;
721 if (equation[i+2].token.operatr == POWER
722 && equation[i+2].level == level
723 && equation[i+1].kind == CONSTANT
724 && equation[i+3].level == level
725 && equation[i+3].kind == CONSTANT) {
726 if (fmod(equation[i+1].token.constant, 1.0) == 0.0) {
727 if (!rationalize_denominators
728 || !isfinite(equation[i+3].token.constant)
729 || equation[i+3].token.constant <= 0.0
730 || equation[i+3].token.constant >= 1.0) {
731 continue;
732 }
733 if (*np + 2 > n_tokens)
734 error_huge();
735 equation[i+3].token.constant -= 1.0;
736 blt(&equation[i+2], &equation[i], (*np - i) * sizeof(token_type));
737 *np += 2;
738 i++;
739 equation[i].level = level - 1;
740 equation[i].kind = CONSTANT;
741 equation[i].token.constant = equation[i+2].token.constant;
742 i++;
743 } else {
744 equation[i].token.operatr = TIMES;
745 equation[i+1].token.constant = 1.0 / equation[i+1].token.constant;
746 }
747 modified = true;
748 }
749 }
750 }
751 return modified;
752 }
753