1 /*
2 * Group and combine algebraic fractions for Mathomatic.
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 sf_recurse(token_type *equation, int *np, int loc, int level, int start_flag);
28 static int sf_sub(token_type *equation, int *np, int loc, int i1, int n1, int i2, int n2, int level, int start_flag);
29
30 static void
group_recurse(equation,np,loc,level)31 group_recurse(equation, np, loc, level)
32 token_type *equation; /* equation side pointer */
33 int *np; /* pointer to length of equation side */
34 int loc; /* starting location within equation side */
35 int level; /* current level of parentheses within the sub-expression starting at "loc" */
36 {
37 int i;
38 int len;
39 int di = -1, edi;
40 int group_flag = false;
41 int e1;
42
43 for (i = loc; i < *np && equation[i].level >= level;) {
44 if (equation[i].level > level) {
45 group_recurse(equation, np, i, level + 1);
46 i++;
47 for (; i < *np && equation[i].level > level; i += 2)
48 ;
49 continue;
50 }
51 i++;
52 }
53 edi = e1 = i;
54 for (i = loc + 1; i < e1; i += 2) {
55 if (equation[i].level == level) {
56 if (equation[i].token.operatr == DIVIDE) {
57 if (di < 0) {
58 di = i;
59 continue;
60 }
61 group_flag = true;
62 for (len = i + 2; len < e1; len += 2) {
63 if (equation[len].level == level && equation[len].token.operatr != DIVIDE)
64 break;
65 }
66 len -= i;
67 if (edi == e1) {
68 i += len;
69 edi = i;
70 continue;
71 }
72 blt(scratch, &equation[i], len * sizeof(token_type));
73 blt(&equation[di+len], &equation[di], (i - di) * sizeof(token_type));
74 blt(&equation[di], scratch, len * sizeof(token_type));
75 edi += len;
76 i += len - 2;
77 } else {
78 if (di >= 0 && edi == e1)
79 edi = i;
80 }
81 }
82 }
83 if (group_flag) {
84 for (i = di + 1; i < edi; i++) {
85 if (equation[i].level == level && equation[i].kind == OPERATOR) {
86 #if DEBUG
87 if (equation[i].token.operatr != DIVIDE) {
88 error_bug("Bug in group_recurse().");
89 }
90 #endif
91 equation[i].token.operatr = TIMES;
92 }
93 equation[i].level++;
94 }
95 }
96 }
97
98 /*
99 * Group denominators of algebraic fractions together in an equation side.
100 * Grouping here means converting "a/b/c/d*e" to "a*e/(b*c*d)" or "a/(b*c*d)*e".
101 * Not guaranteed to put the grouped divisors last, reorder() puts divisors last.
102 */
103 void
group_proc(equation,np)104 group_proc(equation, np)
105 token_type *equation; /* equation side pointer */
106 int *np; /* pointer to length of equation side */
107 {
108 group_recurse(equation, np, 0, 1);
109 }
110
111 /*
112 * Make equation side ready for display.
113 * Basically simplify, then convert non-integer constants in an equation side to fractions,
114 * when exactly equal to simple fractions.
115 * Also groups algebraic fraction denominators with group_proc() above.
116 *
117 * Return true if any fractions were created.
118 */
119 int
fractions_and_group(equation,np)120 fractions_and_group(equation, np)
121 token_type *equation;
122 int *np;
123 {
124 int rv = false;
125
126 elim_loop(equation, np);
127 if (fractions_display) {
128 rv = make_fractions(equation, np);
129 }
130 group_proc(equation, np);
131 return rv;
132 }
133
134 /*
135 * This function is the guts of the display command.
136 * Makes an equation space ready for display.
137 *
138 * Return true if any fractions were created.
139 */
140 int
make_fractions_and_group(n)141 make_fractions_and_group(n)
142 int n; /* equation space number */
143 {
144 int rv = false;
145
146 if (empty_equation_space(n))
147 return false;
148 if (fractions_and_group(lhs[n], &n_lhs[n]))
149 rv = true;
150 if (n_rhs[n] > 0) {
151 if (fractions_and_group(rhs[n], &n_rhs[n]))
152 rv = true;
153 }
154 return rv;
155 }
156
157 /*
158 * Efficiently combine algebraic fractions added together
159 * by putting all terms over a common denominator.
160 * This means converting "(a/b)+(c/d)+f" directly to "(a*d+c*b+b*d*f)/b/d".
161 * The resulting expression is always equivalent to the original expression.
162 *
163 * If start_flag is 0, only combine denominators to convert complex fractions to simple fractions.
164 * Level one addition of fractions will be unchanged.
165 * Can make an expression over-complicated.
166 *
167 * If start_flag is 1, always combine denominators no matter what they are.
168 * This can easily make an expression large and complicated,
169 * but always results in a single, simple fraction.
170 * Used when solving for zero.
171 *
172 * If start_flag is 2, combine denominators
173 * and reduce the result by removing any polynomial GCD between them.
174 * Note that this wipes out globals tlhs[] and trhs[].
175 * This usually results in the simplest and most correct result.
176 *
177 * If start_flag is 3: same as start_flag = 2,
178 * except absolute value and imaginary denominators are combined, too.
179 * This simplifies all algebraic fractions into a single simple fraction.
180 * Removing the polynomial GCD usually prevents making a more complicated
181 * (or larger) algebraic fraction.
182 * Used by the fraction command.
183 *
184 * Return true if the equation side was modified.
185 */
186 int
super_factor(equation,np,start_flag)187 super_factor(equation, np, start_flag)
188 token_type *equation; /* pointer to the beginning of the equation side to process */
189 int *np; /* pointer to the length of the equation side */
190 int start_flag;
191 {
192 int rv;
193
194 group_proc(equation, np);
195 rv = sf_recurse(equation, np, 0, 1, start_flag);
196 organize(equation, np);
197 return rv;
198 }
199
200 static int
sf_recurse(equation,np,loc,level,start_flag)201 sf_recurse(equation, np, loc, level, start_flag)
202 token_type *equation;
203 int *np, loc, level, start_flag;
204 {
205 int modified = false;
206 int i, j, k;
207 int op;
208 int len1, len2;
209
210 if (!start_flag) {
211 for (i = loc + 1; i < *np && equation[i].level >= level; i += 2) {
212 if (equation[i].level == level && equation[i].token.operatr == DIVIDE) {
213 start_flag = true;
214 break;
215 }
216 }
217 }
218 op = 0;
219 for (i = loc; i < *np && equation[i].level >= level;) {
220 if (equation[i].level > level) {
221 modified |= sf_recurse(equation, np, i, level + 1, start_flag);
222 i++;
223 for (; i < *np && equation[i].level > level; i += 2)
224 ;
225 continue;
226 } else if (equation[i].kind == OPERATOR) {
227 op = equation[i].token.operatr;
228 }
229 i++;
230 }
231 if (modified || !start_flag)
232 return modified;
233 switch (op) {
234 case PLUS:
235 case MINUS:
236 break;
237 default:
238 return modified;
239 }
240 sf_again:
241 i = loc;
242 for (k = i + 1; k < *np && equation[k].level > level; k += 2)
243 ;
244 len1 = k - i;
245 for (j = i + len1 + 1; j < *np && equation[j-1].level >= level; j += len2 + 1) {
246 for (k = j + 1; k < *np && equation[k].level > level; k += 2)
247 ;
248 len2 = k - j;
249 #if 0
250 side_debug(0, &equation[i], len1);
251 side_debug(0, &equation[j], len2);
252 #endif
253 if (sf_sub(equation, np, loc, i, len1, j, len2, level + 1, start_flag)) {
254 #if 0
255 int junk;
256 printf("start_flag = %d\n", start_flag);
257 debug_string(0, "super_factor() successful:");
258 for (junk = 1; (loc + junk) < *np && equation[loc+junk].level > level; junk += 2)
259 ;
260 side_debug(0, &equation[loc], junk);
261 #endif
262 modified = true;
263 goto sf_again;
264 }
265 }
266 return modified;
267 }
268
269 static int
sf_sub(equation,np,loc,i1,n1,i2,n2,level,start_flag)270 sf_sub(equation, np, loc, i1, n1, i2, n2, level, start_flag)
271 token_type *equation;
272 int *np, loc, i1, n1, i2, n2, level, start_flag;
273 {
274 int i, j, k;
275 int b1, b2;
276 int len;
277 int e1, e2;
278 int op1, op2;
279 token_type *p1, *p2;
280 int np1, np2;
281 int div_flag1 = false, div_flag2 = false;
282 int rv;
283
284 e1 = i1 + n1;
285 e2 = i2 + n2;
286 op2 = equation[i2-1].token.operatr;
287 if (i1 <= loc) {
288 op1 = PLUS;
289 } else {
290 op1 = equation[i1-1].token.operatr;
291 }
292 for (i = i1 + 1; i < e1; i += 2) {
293 if (equation[i].level == level && equation[i].token.operatr == DIVIDE) {
294 div_flag1 = true;
295 break;
296 }
297 }
298 b1 = i + 1;
299 if (div_flag1) {
300 for (i += 2; i < e1; i += 2) {
301 if (equation[i].level <= level)
302 break;
303 }
304 }
305 for (j = i2 + 1; j < e2; j += 2) {
306 if (equation[j].level == level && equation[j].token.operatr == DIVIDE) {
307 div_flag2 = true;
308 break;
309 }
310 }
311 b2 = j + 1;
312 if (div_flag2) {
313 for (j += 2; j < e2; j += 2) {
314 if (equation[j].level <= level)
315 break;
316 }
317 }
318 if (!div_flag1 && !div_flag2)
319 return false;
320 #if 1
321 if (start_flag <= 2 && start_flag != 1) {
322 #if 0 /* Leave as 0; 1 will not simplify many imaginary number expressions, for example the tangent expansion. */
323 if (div_flag1 && found_var(&equation[b1], i - b1, IMAGINARY)) {
324 return false;
325 }
326 if (div_flag2 && found_var(&equation[b2], j - b2, IMAGINARY)) {
327 return false;
328 }
329 #endif
330 if (div_flag1) {
331 if (exp_is_absolute(&equation[b1], i - b1))
332 return false;
333 }
334 if (div_flag2) {
335 if (exp_is_absolute(&equation[b2], j - b2))
336 return false;
337 }
338 }
339 #endif
340 if (start_flag >= 2 && div_flag1 && div_flag2) {
341 #if DEBUG
342 debug_string(1, "Trying to find a polynomial GCD between 2 denominators in sf_sub():");
343 side_debug(1, &equation[b1], i - b1);
344 side_debug(1, &equation[b2], j - b2);
345 #endif
346 if ((rv = poly2_gcd(&equation[b1], i - b1, &equation[b2], j - b2, 0L, true)) > 0) {
347 p1 = tlhs;
348 np1 = n_tlhs;
349 p2 = trhs;
350 np2 = n_trhs;
351 goto do_gcd_super;
352 }
353 if (rv == 0 && poly2_gcd(&equation[b2], j - b2, &equation[b1], i - b1, 0L, true) > 0) {
354 p1 = trhs;
355 np1 = n_trhs;
356 p2 = tlhs;
357 np2 = n_tlhs;
358 goto do_gcd_super;
359 }
360 #if DEBUG
361 debug_string(1, "Done; polynomial GCD not found.");
362 #endif
363 }
364 if (n1 + n2 + (i - b1) + (j - b2) + 8 > n_tokens) {
365 error_huge();
366 }
367 if (!div_flag1) {
368 for (k = i1; k < e1; k++)
369 equation[k].level++;
370 }
371 if (!div_flag2) {
372 for (k = i2; k < e2; k++)
373 equation[k].level++;
374 }
375 len = (b1 - i1) - 1;
376 blt(scratch, &equation[i1], len * sizeof(token_type));
377 if (op1 == MINUS) {
378 scratch[len].level = level;
379 scratch[len].kind = OPERATOR;
380 scratch[len].token.operatr = TIMES;
381 len++;
382 scratch[len].level = level;
383 scratch[len].kind = CONSTANT;
384 scratch[len].token.constant = -1.0;
385 len++;
386 }
387 if (div_flag1) {
388 blt(&scratch[len], &equation[i], (e1 - i) * sizeof(token_type));
389 len += e1 - i;
390 }
391 if (div_flag2) {
392 scratch[len].level = level;
393 scratch[len].kind = OPERATOR;
394 scratch[len].token.operatr = TIMES;
395 len++;
396 blt(&scratch[len], &equation[b2], (j - b2) * sizeof(token_type));
397 len += j - b2;
398 }
399 for (k = 0; k < len; k++)
400 scratch[k].level += 2;
401 scratch[len].level = level + 1;
402 scratch[len].kind = OPERATOR;
403 scratch[len].token.operatr = op2;
404 len++;
405 k = len;
406 blt(&scratch[len], &equation[i2], (b2 - i2 - 1) * sizeof(token_type));
407 len += b2 - i2 - 1;
408 if (div_flag2) {
409 blt(&scratch[len], &equation[j], (e2 - j) * sizeof(token_type));
410 len += e2 - j;
411 }
412 if (div_flag1) {
413 scratch[len].level = level;
414 scratch[len].kind = OPERATOR;
415 scratch[len].token.operatr = TIMES;
416 len++;
417 blt(&scratch[len], &equation[b1], (i - b1) * sizeof(token_type));
418 len += i - b1;
419 }
420 for (; k < len; k++)
421 scratch[k].level += 2;
422 scratch[len].level = level;
423 scratch[len].kind = OPERATOR;
424 scratch[len].token.operatr = DIVIDE;
425 len++;
426 k = len;
427 if (div_flag1) {
428 blt(&scratch[len], &equation[b1], (i - b1) * sizeof(token_type));
429 len += i - b1;
430 }
431 if (div_flag1 && div_flag2) {
432 scratch[len].level = level;
433 scratch[len].kind = OPERATOR;
434 scratch[len].token.operatr = TIMES;
435 len++;
436 }
437 if (div_flag2) {
438 blt(&scratch[len], &equation[b2], (j - b2) * sizeof(token_type));
439 len += j - b2;
440 }
441 for (; k < len; k++)
442 scratch[k].level++;
443 end_mess:
444 if (*np + len - n1 - (n2 + 1) > n_tokens) {
445 error_huge();
446 }
447 if (op1 == MINUS) {
448 equation[i1-1].token.operatr = PLUS;
449 }
450 blt(&equation[i2-1], &equation[e2], (*np - e2) * sizeof(token_type));
451 *np -= n2 + 1;
452 blt(&equation[i1+len], &equation[e1], (*np - e1) * sizeof(token_type));
453 *np += len - n1;
454 blt(&equation[i1], scratch, len * sizeof(token_type));
455 return true;
456
457 do_gcd_super:
458 #if DEBUG
459 debug_string(1, "Found a polynomial GCD between 2 denominators in sf_sub().");
460 #endif
461
462 if (5 - i1 + e1 + (2*np2) + b2 - i2 + e2 - j + np1 > n_tokens) {
463 error_huge();
464 }
465 for (k = 0; k < np1; k++) {
466 p1[k].level += level;
467 }
468 for (k = 0; k < np2; k++) {
469 p2[k].level += level;
470 }
471 len = (b1 - i1) - 1;
472 blt(scratch, &equation[i1], len * sizeof(token_type));
473 if (op1 == MINUS) {
474 scratch[len].level = level;
475 scratch[len].kind = OPERATOR;
476 scratch[len].token.operatr = TIMES;
477 len++;
478 scratch[len].level = level;
479 scratch[len].kind = CONSTANT;
480 scratch[len].token.constant = -1.0;
481 len++;
482 }
483 /* if (div_flag1) { */
484 blt(&scratch[len], &equation[i], (e1 - i) * sizeof(token_type));
485 len += e1 - i;
486 /* } */
487 /* if (div_flag2) { */
488 scratch[len].level = level;
489 scratch[len].kind = OPERATOR;
490 scratch[len].token.operatr = TIMES;
491 len++;
492 blt(&scratch[len], p2, np2 * sizeof(token_type));
493 len += np2;
494 /* } */
495 for (k = 0; k < len; k++)
496 scratch[k].level += 2;
497 scratch[len].level = level + 1;
498 scratch[len].kind = OPERATOR;
499 scratch[len].token.operatr = op2;
500 len++;
501 k = len;
502 blt(&scratch[len], &equation[i2], (b2 - i2 - 1) * sizeof(token_type));
503 len += b2 - i2 - 1;
504 /* if (div_flag2) { */
505 blt(&scratch[len], &equation[j], (e2 - j) * sizeof(token_type));
506 len += e2 - j;
507 /* } */
508 /* if (div_flag1) { */
509 scratch[len].level = level;
510 scratch[len].kind = OPERATOR;
511 scratch[len].token.operatr = TIMES;
512 len++;
513 blt(&scratch[len], p1, np1 * sizeof(token_type));
514 len += np1;
515 /* } */
516 for (; k < len; k++)
517 scratch[k].level += 2;
518 scratch[len].level = level;
519 scratch[len].kind = OPERATOR;
520 scratch[len].token.operatr = DIVIDE;
521 len++;
522 k = len;
523 /* if (div_flag1) { */
524 blt(&scratch[len], &equation[b1], (i - b1) * sizeof(token_type));
525 len += (i - b1);
526 /* } */
527 /* if (div_flag1 && div_flag2) { */
528 scratch[len].level = level;
529 scratch[len].kind = OPERATOR;
530 scratch[len].token.operatr = TIMES;
531 len++;
532 /* } */
533 /* if (div_flag2) { */
534 blt(&scratch[len], p2, np2 * sizeof(token_type));
535 len += np2;
536 /* } */
537 for (; k < len; k++)
538 scratch[k].level++;
539 goto end_mess;
540 }
541