1 /*
2 Qalculate (library)
3
4 Copyright (C) 2003-2007, 2008, 2016-2019 Hanna Knutsson (hanna.knutsson@protonmail.com)
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10 */
11
12 #include "support.h"
13
14 #include "MathStructure.h"
15 #include "Calculator.h"
16 #include "BuiltinFunctions.h"
17 #include "Number.h"
18 #include "Function.h"
19 #include "Variable.h"
20 #include "Unit.h"
21 #include "Prefix.h"
22
23 #if HAVE_UNORDERED_MAP
24 # include <unordered_map>
25 using std::unordered_map;
26 #elif defined(__GNUC__)
27
28 # ifndef __has_include
29 # define __has_include(x) 0
30 # endif
31
32 # if (defined(__clang__) && __has_include(<tr1/unordered_map>)) || (__GNUC__ >= 4 && __GNUC_MINOR__ >= 3)
33 # include <tr1/unordered_map>
34 namespace Sgi = std;
35 # define unordered_map std::tr1::unordered_map
36 # else
37 # if __GNUC__ < 3
38 # include <hash_map.h>
39 namespace Sgi { using ::hash_map; }; // inherit globals
40 # else
41 # include <ext/hash_map>
42 # if __GNUC__ == 3 && __GNUC_MINOR__ == 0
43 namespace Sgi = std; // GCC 3.0
44 # else
45 namespace Sgi = ::__gnu_cxx; // GCC 3.1 and later
46 # endif
47 # endif
48 # define unordered_map Sgi::hash_map
49 # endif
50 #else // ... there are other compilers, right?
51 namespace Sgi = std;
52 # define unordered_map Sgi::hash_map
53 #endif
54
55 #include "MathStructure-support.h"
56
57 using std::string;
58 using std::cout;
59 using std::vector;
60 using std::endl;
61
sr_gcd(const MathStructure & m1,const MathStructure & m2,MathStructure & mgcd,const sym_desc_vec & sym_stats,size_t var_i,const EvaluationOptions & eo)62 bool sr_gcd(const MathStructure &m1, const MathStructure &m2, MathStructure &mgcd, const sym_desc_vec &sym_stats, size_t var_i, const EvaluationOptions &eo) {
63
64 if(var_i >= sym_stats.size()) return false;
65 const MathStructure &xvar = sym_stats[var_i].sym;
66
67 MathStructure c, d;
68 Number adeg = m1.degree(xvar);
69 Number bdeg = m2.degree(xvar);
70 Number cdeg, ddeg;
71 if(adeg.isGreaterThanOrEqualTo(bdeg)) {
72 c = m1;
73 d = m2;
74 cdeg = adeg;
75 ddeg = bdeg;
76 } else {
77 c = m2;
78 d = m1;
79 cdeg = bdeg;
80 ddeg = adeg;
81 }
82
83 MathStructure cont_c, cont_d;
84 c.polynomialContent(xvar, cont_c, eo);
85 d.polynomialContent(xvar, cont_d, eo);
86
87 MathStructure gamma;
88
89 if(!MathStructure::gcd(cont_c, cont_d, gamma, eo, NULL, NULL, false)) return false;
90 mgcd = gamma;
91 if(ddeg.isZero()) {
92 return true;
93 }
94
95 MathStructure prim_c, prim_d;
96 c.polynomialPrimpart(xvar, cont_c, prim_c, eo);
97 d.polynomialPrimpart(xvar, cont_d, prim_d, eo);
98 c = prim_c;
99 d = prim_d;
100
101 MathStructure r;
102 MathStructure ri(1, 1, 0);
103 MathStructure psi(1, 1, 0);
104 Number delta(cdeg);
105 delta -= ddeg;
106
107 while(true) {
108
109 if(CALCULATOR->aborted()) return false;
110
111 if(!prem(c, d, xvar, r, eo, false)) return false;
112
113 if(r.isZero()) {
114 mgcd = gamma;
115 MathStructure mprim;
116 d.polynomialPrimpart(xvar, mprim, eo);
117 if(CALCULATOR->aborted()) return false;
118 mgcd.calculateMultiply(mprim, eo);
119 return true;
120 }
121
122 c = d;
123 cdeg = ddeg;
124
125 MathStructure psi_pow(psi);
126 psi_pow.calculateRaise(delta, eo);
127 ri.calculateMultiply(psi_pow, eo);
128
129 if(!divide_in_z(r, ri, d, sym_stats, var_i, eo)) {
130 return false;
131 }
132
133 ddeg = d.degree(xvar);
134 if(ddeg.isZero()) {
135 if(r.isNumber()) {
136 mgcd = gamma;
137 } else {
138 r.polynomialPrimpart(xvar, mgcd, eo);
139 if(CALCULATOR->aborted()) return false;
140 mgcd.calculateMultiply(gamma, eo);
141 }
142 return true;
143 }
144
145 c.lcoefficient(xvar, ri);
146 if(delta.isOne()) {
147 psi = ri;
148 } else if(!delta.isZero()) {
149 MathStructure ri_pow(ri);
150 ri_pow.calculateRaise(delta, eo);
151 MathStructure psi_pow(psi);
152 delta--;
153 psi_pow.calculateRaise(delta, eo);
154 divide_in_z(ri_pow, psi_pow, psi, sym_stats, var_i + 1, eo);
155 }
156 delta = cdeg;
157 delta -= ddeg;
158
159 }
160
161 return false;
162
163 }
164
165
interpolate(const MathStructure & gamma,const Number & xi,const MathStructure & xvar,MathStructure & minterp,const EvaluationOptions & eo)166 bool interpolate(const MathStructure &gamma, const Number &xi, const MathStructure &xvar, MathStructure &minterp, const EvaluationOptions &eo) {
167 MathStructure e(gamma);
168 Number rxi(xi);
169 rxi.recip();
170 minterp.clear();
171 for(long int i = 0; !e.isZero(); i++) {
172 if(CALCULATOR->aborted()) return false;
173 MathStructure gi;
174 polynomial_smod(e, xi, gi, eo);
175 if(minterp.isZero() && !gi.isZero()) {
176 minterp = gi;
177 if(i != 0) {
178 if(minterp.isOne()) {
179 minterp = xvar;
180 if(i != 1) minterp.raise(i);
181 } else {
182 minterp.multiply(xvar, true);
183 if(i != 1) minterp[minterp.size() - 1].raise(i);
184 minterp.calculateMultiplyLast(eo);
185 }
186 }
187 } else if(!gi.isZero()) {
188 minterp.add(gi, true);
189 if(i != 0) {
190 if(minterp[minterp.size() - 1].isOne()) {
191 minterp[minterp.size() - 1] = xvar;
192 if(i != 1) minterp[minterp.size() - 1].raise(i);
193 } else {
194 minterp[minterp.size() - 1].multiply(xvar, true);
195 if(i != 1) minterp[minterp.size() - 1][minterp[minterp.size() - 1].size() - 1].raise(i);
196 minterp[minterp.size() - 1].calculateMultiplyLast(eo);
197 }
198 }
199 }
200 if(!gi.isZero()) e.calculateSubtract(gi, eo);
201 e.calculateMultiply(rxi, eo);
202 }
203 minterp.calculatesub(eo, eo, false);
204 return true;
205 }
206
heur_gcd(const MathStructure & m1,const MathStructure & m2,MathStructure & mgcd,const EvaluationOptions & eo,MathStructure * ca,MathStructure * cb,sym_desc_vec & sym_stats,size_t var_i)207 bool heur_gcd(const MathStructure &m1, const MathStructure &m2, MathStructure &mgcd, const EvaluationOptions &eo, MathStructure *ca, MathStructure *cb, sym_desc_vec &sym_stats, size_t var_i) {
208
209 if(m1.isZero() || m2.isZero()) return false;
210
211 if(m1.isNumber() && m2.isNumber()) {
212 mgcd = m1;
213 if(!mgcd.number().gcd(m2.number())) mgcd.set(1, 1, 0);
214 if(ca) {
215 *ca = m1;
216 ca->number() /= mgcd.number();
217 }
218 if(cb) {
219 *cb = m2;
220 cb->number() /= mgcd.number();
221 }
222 return true;
223 }
224
225 if(var_i >= sym_stats.size()) return false;
226 const MathStructure &xvar = sym_stats[var_i].sym;
227
228 Number nr_gc;
229 integer_content(m1, nr_gc);
230 Number nr_rgc;
231 integer_content(m2, nr_rgc);
232 nr_gc.gcd(nr_rgc);
233 nr_rgc = nr_gc;
234 nr_rgc.recip();
235 MathStructure p(m1);
236 p.calculateMultiply(nr_rgc, eo);
237 MathStructure q(m2);
238 q.calculateMultiply(nr_rgc, eo);
239 Number maxdeg(p.degree(xvar));
240 Number maxdeg2(q.degree(xvar));
241 if(maxdeg2.isGreaterThan(maxdeg)) maxdeg = maxdeg2;
242 Number mp(p.maxCoefficient());
243 Number mq(q.maxCoefficient());
244 Number xi;
245 if(mp.isGreaterThan(mq)) {
246 xi = mq;
247 } else {
248 xi = mp;
249 }
250 xi *= 2;
251 xi += 2;
252
253 for(int t = 0; t < 6; t++) {
254
255 if(CALCULATOR->aborted()) return false;
256
257 if(!xi.isInteger() || (maxdeg * xi.integerLength()).isGreaterThan(100000L)) {
258 return false;
259 }
260
261 MathStructure cp, cq;
262 MathStructure gamma;
263 MathStructure psub(p);
264 psub.calculateReplace(xvar, xi, eo, true);
265 MathStructure qsub(q);
266 qsub.calculateReplace(xvar, xi, eo, true);
267
268 if(heur_gcd(psub, qsub, gamma, eo, &cp, &cq, sym_stats, var_i + 1)) {
269
270 if(!interpolate(gamma, xi, xvar, mgcd, eo)) return false;
271
272 Number ig;
273 integer_content(mgcd, ig);
274 ig.recip();
275 mgcd.calculateMultiply(ig, eo);
276
277 MathStructure dummy;
278
279 if(divide_in_z(p, mgcd, ca ? *ca : dummy, sym_stats, var_i, eo) && divide_in_z(q, mgcd, cb ? *cb : dummy, sym_stats, var_i, eo)) {
280 mgcd.calculateMultiply(nr_gc, eo);
281 return true;
282 }
283 }
284
285 Number xi2(xi);
286 xi2.isqrt();
287 xi2.isqrt();
288 xi *= xi2;
289 xi *= 73794L;
290 xi.iquo(27011L);
291
292 }
293
294 return false;
295
296 }
297
lcm(const MathStructure & m1,const MathStructure & m2,MathStructure & mlcm,const EvaluationOptions & eo,bool check_args)298 bool MathStructure::lcm(const MathStructure &m1, const MathStructure &m2, MathStructure &mlcm, const EvaluationOptions &eo, bool check_args) {
299 if(m1.isNumber() && m2.isNumber()) {
300 mlcm = m1;
301 return mlcm.number().lcm(m2.number());
302 }
303 if(check_args && (!m1.isRationalPolynomial() || !m2.isRationalPolynomial())) {
304 return false;
305 }
306 MathStructure ca, cb;
307 if(!MathStructure::gcd(m1, m2, mlcm, eo, &ca, &cb, false)) return false;
308 mlcm.calculateMultiply(ca, eo);
309 mlcm.calculateMultiply(cb, eo);
310 return true;
311 }
312
has_noninteger_coefficient(const MathStructure & mstruct)313 bool has_noninteger_coefficient(const MathStructure &mstruct) {
314 if(mstruct.isNumber() && mstruct.number().isRational() && !mstruct.number().isInteger()) return true;
315 if(mstruct.isFunction() || mstruct.isPower()) return false;
316 for(size_t i = 0; i < mstruct.size(); i++) {
317 if(has_noninteger_coefficient(mstruct[i])) return true;
318 }
319 return false;
320 }
321
gcd(const MathStructure & m1,const MathStructure & m2,MathStructure & mresult,const EvaluationOptions & eo2,MathStructure * ca,MathStructure * cb,bool check_args)322 bool MathStructure::gcd(const MathStructure &m1, const MathStructure &m2, MathStructure &mresult, const EvaluationOptions &eo2, MathStructure *ca, MathStructure *cb, bool check_args) {
323
324 EvaluationOptions eo = eo2;
325 eo.keep_zero_units = false;
326
327 if(ca) *ca = m1;
328 if(cb) *cb = m2;
329 mresult.set(1, 1, 0);
330
331 if(CALCULATOR->aborted()) return false;
332
333 if(m1.isOne() || m2.isOne()) {
334 return true;
335 }
336 if(m1.isNumber() && m2.isNumber()) {
337 mresult = m1;
338 if(!mresult.number().gcd(m2.number())) mresult.set(1, 1, 0);
339 if(ca || cb) {
340 if(mresult.isZero()) {
341 if(ca) ca->clear();
342 if(cb) cb->clear();
343 } else {
344 if(ca) {
345 *ca = m1;
346 ca->number() /= mresult.number();
347 }
348 if(cb) {
349 *cb = m2;
350 cb->number() /= mresult.number();
351 }
352 }
353 }
354 return true;
355 }
356
357 if(m1 == m2) {
358 if(ca) ca->set(1, 1, 0);
359 if(cb) cb->set(1, 1, 0);
360 mresult = m1;
361 return true;
362 }
363 if(m1.isZero()) {
364 if(ca) ca->clear();
365 if(cb) cb->set(1, 1, 0);
366 mresult = m2;
367 return true;
368 }
369 if(m2.isZero()) {
370 if(ca) ca->set(1, 1, 0);
371 if(cb) cb->clear();
372 mresult = m1;
373 return true;
374 }
375
376 if(check_args && (!m1.isRationalPolynomial() || !m2.isRationalPolynomial())) {
377 return false;
378 }
379
380 if(has_noninteger_coefficient(m1) || has_noninteger_coefficient(m2)) {
381 Number nlcm1;
382 lcm_of_coefficients_denominators(m1, nlcm1);
383 Number nlcm2;
384 lcm_of_coefficients_denominators(m2, nlcm2);
385 nlcm1.lcm(nlcm2);
386 if(!nlcm1.isOne()) {
387 MathStructure mtmp1, mtmp2;
388 multiply_lcm(m1, nlcm1, mtmp1, eo);
389 multiply_lcm(m2, nlcm1, mtmp2, eo);
390 if(mtmp1.equals(m1, true, true) || mtmp2.equals(m2, true, true)) return false;
391 return gcd(mtmp1, mtmp2, mresult, eo, ca, cb, false);
392 }
393 }
394
395 if(m1.isMultiplication()) {
396 if(m2.isMultiplication() && m2.size() > m1.size())
397 goto factored_2;
398 factored_1:
399 mresult.clear();
400 mresult.setType(STRUCT_MULTIPLICATION);
401 MathStructure acc_ca;
402 acc_ca.setType(STRUCT_MULTIPLICATION);
403 MathStructure part_2(m2);
404 MathStructure part_ca, part_cb;
405 for(size_t i = 0; i < m1.size(); i++) {
406 mresult.addChild(m_zero);
407 MathStructure::gcd(m1[i], part_2, mresult[i], eo, &part_ca, &part_cb, false);
408 if(CALCULATOR->aborted()) return false;
409 acc_ca.addChild(part_ca);
410 part_2 = part_cb;
411 }
412 mresult.calculatesub(eo, eo, false);
413 if(ca) {
414 *ca = acc_ca;
415 ca->calculatesub(eo, eo, false);
416 }
417 if(cb) *cb = part_2;
418 return true;
419 } else if(m2.isMultiplication()) {
420 if(m1.isMultiplication() && m1.size() > m2.size())
421 goto factored_1;
422 factored_2:
423 mresult.clear();
424 mresult.setType(STRUCT_MULTIPLICATION);
425 MathStructure acc_cb;
426 acc_cb.setType(STRUCT_MULTIPLICATION);
427 MathStructure part_1(m1);
428 MathStructure part_ca, part_cb;
429 for(size_t i = 0; i < m2.size(); i++) {
430 mresult.addChild(m_zero);
431 MathStructure::gcd(part_1, m2[i], mresult[i], eo, &part_ca, &part_cb, false);
432 if(CALCULATOR->aborted()) return false;
433 acc_cb.addChild(part_cb);
434 part_1 = part_ca;
435 }
436 mresult.calculatesub(eo, eo, false);
437 if(ca) *ca = part_1;
438 if(cb) {
439 *cb = acc_cb;
440 cb->calculatesub(eo, eo, false);
441 }
442 return true;
443 }
444
445 if(m1.isPower()) {
446 MathStructure p(m1[0]);
447 if(m2.isPower()) {
448 if(m1[0] == m2[0]) {
449 if(m1[1].number().isLessThan(m2[1].number())) {
450 if(ca) ca->set(1, 1, 0);
451 if(cb) {
452 *cb = m2;
453 (*cb)[1].number() -= m1[1].number();
454 POWER_CLEAN((*cb))
455 }
456 mresult = m1;
457 return true;
458 } else {
459 if(ca) {
460 *ca = m1;
461 (*ca)[1].number() -= m2[1].number();
462 POWER_CLEAN((*ca))
463 }
464 if(cb) cb->set(1, 1, 0);
465 mresult = m2;
466 return true;
467 }
468 } else {
469 MathStructure p_co, pb_co;
470 MathStructure p_gcd;
471 if(!MathStructure::gcd(m1[0], m2[0], p_gcd, eo, &p_co, &pb_co, false)) return false;
472 if(p_gcd.isOne()) {
473 return true;
474 } else {
475 if(m1[1].number().isLessThan(m2[1].number())) {
476 mresult = p_gcd;
477 mresult.calculateRaise(m1[1], eo);
478 mresult.multiply(m_zero, true);
479 p_co.calculateRaise(m1[1], eo);
480 pb_co.calculateRaise(m2[1], eo);
481 p_gcd.raise(m2[1]);
482 p_gcd[1].number() -= m1[1].number();
483 p_gcd.calculateRaiseExponent(eo);
484 p_gcd.calculateMultiply(pb_co, eo);
485 bool b = MathStructure::gcd(p_co, p_gcd, mresult[mresult.size() - 1], eo, ca, cb, false);
486 mresult.calculateMultiplyLast(eo);
487 return b;
488 } else {
489 mresult = p_gcd;
490 mresult.calculateRaise(m2[1], eo);
491 mresult.multiply(m_zero, true);
492 p_co.calculateRaise(m1[1], eo);
493 pb_co.calculateRaise(m2[1], eo);
494 p_gcd.raise(m1[1]);
495 p_gcd[1].number() -= m2[1].number();
496 p_gcd.calculateRaiseExponent(eo);
497 p_gcd.calculateMultiply(p_co, eo);
498 bool b = MathStructure::gcd(p_gcd, pb_co, mresult[mresult.size() - 1], eo, ca, cb, false);
499 mresult.calculateMultiplyLast(eo);
500 return b;
501 }
502 }
503 }
504 } else {
505 if(m1[0] == m2) {
506 if(ca) {
507 *ca = m1;
508 (*ca)[1].number()--;
509 POWER_CLEAN((*ca))
510 }
511 if(cb) cb->set(1, 1, 0);
512 mresult = m2;
513 return true;
514 }
515 MathStructure p_co, bpart_co;
516 MathStructure p_gcd;
517 if(!MathStructure::gcd(m1[0], m2, p_gcd, eo, &p_co, &bpart_co, false)) return false;
518 if(p_gcd.isOne()) {
519 return true;
520 } else {
521 mresult = p_gcd;
522 mresult.multiply(m_zero, true);
523 p_co.calculateRaise(m1[1], eo);
524 p_gcd.raise(m1[1]);
525 p_gcd[1].number()--;
526 p_gcd.calculateRaiseExponent(eo);
527 p_gcd.calculateMultiply(p_co, eo);
528 bool b = MathStructure::gcd(p_gcd, bpart_co, mresult[mresult.size() - 1], eo, ca, cb, false);
529 mresult.calculateMultiplyLast(eo);
530 return b;
531 }
532 }
533 } else if(m2.isPower()) {
534 if(m2[0] == m1) {
535 if(ca) ca->set(1, 1, 0);
536 if(cb) {
537 *cb = m2;
538 (*cb)[1].number()--;
539 POWER_CLEAN((*cb))
540 }
541 mresult = m1;
542 return true;
543 }
544 MathStructure p_co, apart_co;
545 MathStructure p_gcd;
546 if(!MathStructure::gcd(m1, m2[0], p_gcd, eo, &apart_co, &p_co, false)) return false;
547 if(p_gcd.isOne()) {
548 return true;
549 } else {
550 mresult = p_gcd;
551 mresult.multiply(m_zero, true);
552 p_co.calculateRaise(m2[1], eo);
553 p_gcd.raise(m2[1]);
554 p_gcd[1].number()--;
555 p_gcd.calculateRaiseExponent(eo);
556 p_gcd.calculateMultiply(p_co, eo);
557 bool b = MathStructure::gcd(apart_co, p_gcd, mresult[mresult.size() - 1], eo, ca, cb, false);
558 mresult.calculateMultiplyLast(eo);
559 return b;
560 }
561 }
562 if(IS_A_SYMBOL(m1) || m1.isUnit()) {
563 MathStructure bex(m2);
564 bex.calculateReplace(m1, m_zero, eo);
565 if(!bex.isZero()) {
566 return true;
567 }
568 }
569 if(IS_A_SYMBOL(m2) || m2.isUnit()) {
570 MathStructure aex(m1);
571 aex.calculateReplace(m2, m_zero, eo);
572 if(!aex.isZero()) {
573 return true;
574 }
575 }
576
577 sym_desc_vec sym_stats;
578 get_symbol_stats(m1, m2, sym_stats);
579
580 if(sym_stats.empty()) return false;
581
582 size_t var_i = 0;
583 const MathStructure &xvar = sym_stats[var_i].sym;
584
585 /*for(size_t i = 0; i< sym_stats.size(); i++) {
586 if(sym_stats[i].sym.size() > 0) {
587 MathStructure m1b(m1), m2b(m2);
588 vector<UnknownVariable*> vars;
589 for(; i < sym_stats.size(); i++) {
590 if(sym_stats[i].sym.size() > 0) {
591 UnknownVariable *var = new UnknownVariable("", format_and_print(sym_stats[i].sym));
592 m1b.replace(sym_stats[i].sym, var);
593 m2b.replace(sym_stats[i].sym, var);
594 vars.push_back(var);
595 }
596 }
597 bool b = gcd(m1b, m2, mresult, eo2, ca, cb, false);
598 size_t i2 = 0;
599 for(i = 0; i < sym_stats.size(); i++) {
600 if(sym_stats[i].sym.size() > 0) {
601 mresult.replace(vars[i2], sym_stats[i].sym);
602 if(ca) ca->replace(vars[i2], sym_stats[i].sym);
603 if(cb) cb->replace(vars[i2], sym_stats[i].sym);
604 vars[i2]->destroy();
605 i2++;
606 }
607 }
608 return b;
609 }
610 }*/
611
612 Number ldeg_a(sym_stats[var_i].ldeg_a);
613 Number ldeg_b(sym_stats[var_i].ldeg_b);
614 Number min_ldeg;
615 if(ldeg_a.isLessThan(ldeg_b)) {
616 min_ldeg = ldeg_a;
617 } else {
618 min_ldeg = ldeg_b;
619 }
620
621 if(min_ldeg.isPositive()) {
622 MathStructure aex(m1), bex(m2);
623 MathStructure common(xvar);
624 if(!min_ldeg.isOne()) common.raise(min_ldeg);
625 aex.calculateDivide(common, eo);
626 bex.calculateDivide(common, eo);
627 MathStructure::gcd(aex, bex, mresult, eo, ca, cb, false);
628 mresult.calculateMultiply(common, eo);
629 return true;
630 }
631
632 if(sym_stats[var_i].deg_a.isZero()) {
633 if(cb) {
634 MathStructure c, mprim;
635 int u;
636 m2.polynomialUnitContentPrimpart(xvar, u, c, mprim, eo);
637 if(!m2.equals(c, true, true)) MathStructure::gcd(m1, c, mresult, eo, ca, cb, false);
638 if(CALCULATOR->aborted()) {
639 if(ca) *ca = m1;
640 if(cb) *cb = m2;
641 mresult.set(1, 1, 0);
642 return false;
643 }
644 if(u == -1) {
645 cb->calculateNegate(eo);
646 }
647 cb->calculateMultiply(mprim, eo);
648 } else {
649 MathStructure c;
650 m2.polynomialContent(xvar, c, eo);
651 if(!m2.equals(c, true, true)) MathStructure::gcd(m1, c, mresult, eo, ca, cb, false);
652 }
653 return true;
654 } else if(sym_stats[var_i].deg_b.isZero()) {
655 if(ca) {
656 MathStructure c, mprim;
657 int u;
658 m1.polynomialUnitContentPrimpart(xvar, u, c, mprim, eo);
659 if(!m1.equals(c, true, true)) MathStructure::gcd(c, m2, mresult, eo, ca, cb, false);
660 if(CALCULATOR->aborted()) {
661 if(ca) *ca = m1;
662 if(cb) *cb = m2;
663 mresult.set(1, 1, 0);
664 return false;
665 }
666 if(u == -1) {
667 ca->calculateNegate(eo);
668 }
669 ca->calculateMultiply(mprim, eo);
670 } else {
671 MathStructure c;
672 m1.polynomialContent(xvar, c, eo);
673 if(!m1.equals(c, true, true)) MathStructure::gcd(c, m2, mresult, eo, ca, cb, false);
674 }
675 return true;
676 }
677
678 if(!heur_gcd(m1, m2, mresult, eo, ca, cb, sym_stats, var_i)) {
679 if(!sr_gcd(m1, m2, mresult, sym_stats, var_i, eo)) {
680 if(ca) *ca = m1;
681 if(cb) *cb = m2;
682 mresult.set(1, 1, 0);
683 return false;
684 }
685 if(!mresult.isOne()) {
686 if(ca) {
687 if(!MathStructure::polynomialDivide(m1, mresult, *ca, eo, false)) {
688 if(ca) *ca = m1;
689 if(cb) *cb = m2;
690 mresult.set(1, 1, 0);
691 return false;
692 }
693 }
694 if(cb) {
695 if(!MathStructure::polynomialDivide(m2, mresult, *cb, eo, false)) {
696 if(ca) *ca = m1;
697 if(cb) *cb = m2;
698 mresult.set(1, 1, 0);
699 return false;
700 }
701 }
702 }
703 if(CALCULATOR->aborted()) {
704 if(ca) *ca = m1;
705 if(cb) *cb = m2;
706 mresult.set(1, 1, 0);
707 return false;
708 }
709 }
710
711 return true;
712 }
713
714
dos_count_points(const MathStructure & m,bool b_unknown)715 Number dos_count_points(const MathStructure &m, bool b_unknown) {
716 if(m.isPower() && (!b_unknown || m[0].containsUnknowns())) {
717 Number nr = dos_count_points(m[0], b_unknown);
718 if(!nr.isZero()) {
719 if(m[0].isNumber() && m[0].number().isNonZero()) {
720 nr *= m[1].number();
721 if(nr.isNegative()) nr.negate();
722 } else {
723 nr *= 2;
724 }
725 return nr;
726 }
727 } else if(m.size() == 0 && ((b_unknown && m.isUnknown()) || (!b_unknown && !m.isNumber()))) {
728 return nr_one;
729 }
730 Number nr;
731 for(size_t i = 0; i < m.size(); i++) {
732 nr += dos_count_points(m[i], b_unknown);
733 }
734 return nr;
735 }
736
do_simplification(MathStructure & mstruct,const EvaluationOptions & eo,bool combine_divisions,bool only_gcd,bool combine_only,bool recursive,bool limit_size,int i_run)737 bool do_simplification(MathStructure &mstruct, const EvaluationOptions &eo, bool combine_divisions, bool only_gcd, bool combine_only, bool recursive, bool limit_size, int i_run) {
738
739 if(!eo.expand || !eo.assume_denominators_nonzero) return false;
740
741 if(recursive) {
742 bool b = false;
743 for(size_t i = 0; i < mstruct.size(); i++) {
744 b = do_simplification(mstruct[i], eo, combine_divisions, only_gcd || (!mstruct.isComparison() && !mstruct.isLogicalAnd() && !mstruct.isLogicalOr()), combine_only, true, limit_size) || b;
745 if(CALCULATOR->aborted()) return b;
746 }
747 if(b) mstruct.calculatesub(eo, eo, false);
748 return do_simplification(mstruct, eo, combine_divisions, only_gcd, combine_only, false, limit_size) || b;
749 }
750 if(mstruct.isPower() && mstruct[1].isNumber() && mstruct[1].number().isRational() && !mstruct[1].number().isInteger() && mstruct[0].isAddition() && mstruct[0].isRationalPolynomial()) {
751 MathStructure msqrfree(mstruct[0]);
752 if(sqrfree(msqrfree, eo) && msqrfree.isPower() && msqrfree.calculateRaise(mstruct[1], eo)) {
753 mstruct = msqrfree;
754 return true;
755 }
756 } else if(mstruct.isFunction() && mstruct.function()->id() == FUNCTION_ID_ROOT && VALID_ROOT(mstruct) && mstruct[0].isAddition() && mstruct[0].isRationalPolynomial()) {
757 MathStructure msqrfree(mstruct[0]);
758 if(sqrfree(msqrfree, eo) && msqrfree.isPower() && msqrfree[1].isInteger() && msqrfree[1].number().isPositive()) {
759 if(msqrfree[1] == mstruct[1]) {
760 if(msqrfree[1].number().isEven()) {
761 if(!msqrfree[0].representsReal(true)) return false;
762 msqrfree.delChild(2);
763 msqrfree.setType(STRUCT_FUNCTION);
764 msqrfree.setFunctionId(FUNCTION_ID_ABS);
765 mstruct = msqrfree;
766 } else {
767 mstruct = msqrfree[0];
768 }
769 return true;
770 } else if(msqrfree[1].number().isIntegerDivisible(mstruct[1].number())) {
771 if(msqrfree[1].number().isEven()) {
772 if(!msqrfree[0].representsReal(true)) return false;
773 msqrfree[0].transform(STRUCT_FUNCTION);
774 msqrfree[0].setFunctionId(FUNCTION_ID_ABS);
775 }
776 msqrfree[1].number().divide(mstruct[1].number());
777 mstruct = msqrfree;
778 mstruct.calculatesub(eo, eo, false);
779 return true;
780 } else if(mstruct[1].number().isIntegerDivisible(msqrfree[1].number())) {
781 if(msqrfree[1].number().isEven()) {
782 if(!msqrfree[0].representsReal(true)) return false;
783 msqrfree[0].transform(STRUCT_FUNCTION);
784 msqrfree[0].setFunctionId(FUNCTION_ID_ABS);
785 }
786 Number new_root(mstruct[1].number());
787 new_root.divide(msqrfree[1].number());
788 mstruct[0] = msqrfree[0];
789 mstruct[1] = new_root;
790 return true;
791 }
792 }
793 }
794 if(!mstruct.isAddition() && !mstruct.isMultiplication()) return false;
795
796 if(combine_divisions) {
797
798 MathStructure divs, nums, numleft, mleft;
799
800 EvaluationOptions eo2 = eo;
801 eo2.do_polynomial_division = false;
802 eo2.keep_zero_units = false;
803
804 if(!mstruct.isAddition()) mstruct.transform(STRUCT_ADDITION);
805
806 // find division by polynomial
807 for(size_t i = 0; i < mstruct.size(); i++) {
808 if(CALCULATOR->aborted()) return false;
809 bool b = false;
810 if(mstruct[i].isMultiplication()) {
811 MathStructure div, num(1, 1, 0);
812 bool b_num = false;
813 for(size_t i2 = 0; i2 < mstruct[i].size(); i2++) {
814 if(mstruct[i][i2].isPower() && mstruct[i][i2][1].isInteger() && mstruct[i][i2][1].number().isNegative()) {
815 bool b_rat = mstruct[i][i2][0].isRationalPolynomial();
816 if(!b_rat && mstruct[i][i2][0].isAddition() && mstruct[i][i2][0].size() > 1) {
817 b_rat = true;
818 for(size_t i3 = 0; i3 < mstruct[i][i2][0].size(); i3++) {
819 if(!mstruct[i][i2][0][i3].representsZero(true) && !mstruct[i][i2][0][i3].isRationalPolynomial()) {
820 b_rat = false;
821 break;
822 }
823 }
824 }
825 if(!b_rat) {
826 div.clear();
827 break;
828 }
829 bool b_minone = mstruct[i][i2][1].isMinusOne();
830 if(b_minone) {
831 if(div.isZero()) div = mstruct[i][i2][0];
832 else div.multiply(mstruct[i][i2][0], true);
833 } else {
834 mstruct[i][i2][1].number().negate();
835 if(div.isZero()) div = mstruct[i][i2];
836 else div.multiply(mstruct[i][i2], true);
837 mstruct[i][i2][1].number().negate();
838 }
839 } else if(mstruct[i][i2].isRationalPolynomial() || mstruct[i][i2].representsZero(true)) {
840 if(!b_num) {b_num = true; num = mstruct[i][i2];}
841 else num.multiply(mstruct[i][i2], true);
842 } else {
843 div.clear();
844 break;
845 }
846 }
847 if(!div.isZero()) {
848 bool b_found = false;
849 for(size_t i3 = 0; i3 < divs.size(); i3++) {
850 if(divs[i3] == div) {
851 if(!num.representsZero(true)) {
852 if(num.isAddition()) {
853 for(size_t i4 = 0; i4 < num.size(); i4++) {
854 nums[i3].add(num[i4], true);
855 }
856 } else {
857 nums[i3].add(num, true);
858 }
859 }
860 b_found = true;
861 b = true;
862 break;
863 }
864 }
865 if(!b_found && (eo.assume_denominators_nonzero || div.representsNonZero(true)) && !div.representsZero(true)) {
866 if(!num.representsZero(true)) {
867 divs.addChild(div);
868 nums.addChild(num);
869 }
870 b = true;
871 }
872 }
873 } else if(mstruct[i].isPower() && mstruct[i][1].isInteger() && mstruct[i][1].number().isNegative()) {
874 bool b_rat = mstruct[i][0].isRationalPolynomial();
875 if(!b_rat && mstruct[i][0].isAddition() && mstruct[i][0].size() > 1) {
876 b_rat = true;
877 for(size_t i2 = 0; i2 < mstruct[i][0].size(); i2++) {
878 if(!mstruct[i][0][i2].representsZero(true) && !mstruct[i][0].isRationalPolynomial()) {
879 b_rat = false;
880 break;
881 }
882 }
883 }
884 if(b_rat) {
885 bool b_minone = mstruct[i][1].isMinusOne();
886 if(!b_minone) mstruct[i][1].number().negate();
887 bool b_found = false;
888 for(size_t i3 = 0; i3 < divs.size(); i3++) {
889 if((b_minone && divs[i3] == mstruct[i][0]) || (!b_minone && divs[i3] == mstruct[i])) {
890 nums[i3].add(m_one, true);
891 b_found = true;
892 b = true;
893 break;
894 }
895 }
896 if(!b_found && (eo.assume_denominators_nonzero || mstruct[i][0].representsNonZero(true)) && !mstruct[i][0].representsZero(true)) {
897 if(b_minone) divs.addChild(mstruct[i][0]);
898 else divs.addChild(mstruct[i]);
899 nums.addChild(m_one);
900 b = true;
901 }
902 if(!b_minone) mstruct[i][1].number().negate();
903 }
904 }
905 if(!b) {
906 if(i_run <= 1 && mstruct[i].isRationalPolynomial()) numleft.addChild(mstruct[i]);
907 else mleft.addChild(mstruct[i]);
908 }
909 }
910
911 for(size_t i = 0; i < divs.size(); i++) {
912 if(divs[i].isAddition()) {
913 for(size_t i2 = 0; i2 < divs[i].size();) {
914 if(divs[i][i2].representsZero(true)) {
915 divs[i].delChild(i2 + 1);
916 } else i2++;
917 }
918 if(divs[i].size() == 1) divs[i].setToChild(1);
919 else if(divs[i].size() == 0) divs[i].clear();
920 }
921 }
922
923 if(divs.size() == 0) {
924 if(mstruct.size() == 1) mstruct.setToChild(1);
925 return false;
926 }
927 bool b_ret = false;
928 if(divs.size() > 1 || numleft.size() > 0) b_ret = true;
929
930 /*divs.setType(STRUCT_VECTOR);
931 nums.setType(STRUCT_VECTOR);
932 numleft.setType(STRUCT_VECTOR);
933 mleft.setType(STRUCT_VECTOR);
934 cout << nums << ":" << divs << ":" << numleft << ":" << mleft << endl;*/
935
936 if(i_run > 1) {
937 for(size_t i = 0; i < divs.size();) {
938 bool b = true;
939 if(!divs[i].isRationalPolynomial() || !nums[i].isRationalPolynomial()) {
940 if(mstruct.size() == 1) mstruct.setToChild(1);
941 return false;
942 }
943 if(!combine_only && divs[i].isAddition() && nums[i].isAddition()) {
944 MathStructure ca, cb, mgcd;
945 if(MathStructure::gcd(nums[i], divs[i], mgcd, eo2, &ca, &cb, false) && !mgcd.isOne() && (!cb.isAddition() || !ca.isAddition() || ca.size() + cb.size() <= nums[i].size() + divs[i].size())) {
946 if(mgcd.representsNonZero(true) || !eo.warn_about_denominators_assumed_nonzero || warn_about_denominators_assumed_nonzero(mgcd, eo)) {
947 if(cb.isOne()) {
948 numleft.addChild(ca);
949 nums.delChild(i + 1);
950 divs.delChild(i + 1);
951 b = false;
952 } else {
953 nums[i] = ca;
954 divs[i] = cb;
955 }
956 b_ret = true;
957 }
958 }
959 }
960 if(CALCULATOR->aborted()) {if(mstruct.size() == 1) mstruct.setToChild(1); return false;}
961 if(b) i++;
962 }
963
964 } else {
965 while(divs.size() > 0) {
966 bool b = true;
967 if(!divs[0].isRationalPolynomial() || !nums[0].isRationalPolynomial()) {
968 if(mstruct.size() == 1) mstruct.setToChild(1);
969 return false;
970 }
971 if(!combine_only && divs[0].isAddition() && nums[0].isAddition()) {
972 MathStructure ca, cb, mgcd;
973 if(MathStructure::gcd(nums[0], divs[0], mgcd, eo2, &ca, &cb, false) && !mgcd.isOne() && (!cb.isAddition() || !ca.isAddition() || ca.size() + cb.size() <= nums[0].size() + divs[0].size())) {
974 if(mgcd.representsNonZero(true) || !eo.warn_about_denominators_assumed_nonzero || warn_about_denominators_assumed_nonzero(mgcd, eo)) {
975 if(cb.isOne()) {
976 numleft.addChild(ca);
977 nums.delChild(1);
978 divs.delChild(1);
979 b = false;
980 } else {
981 nums[0] = ca;
982 divs[0] = cb;
983 }
984 b_ret = true;
985 }
986 }
987 }
988 if(CALCULATOR->aborted()) {if(mstruct.size() == 1) mstruct.setToChild(1); return false;}
989 if(b && divs.size() > 1) {
990 if(!combine_only && divs[1].isAddition() && nums[1].isAddition()) {
991 MathStructure ca, cb, mgcd;
992 EvaluationOptions eo3 = eo2;
993 eo3.transform_trigonometric_functions = false;
994 if(MathStructure::gcd(nums[1], divs[1], mgcd, eo3, &ca, &cb, false) && !mgcd.isOne() && (!cb.isAddition() || !ca.isAddition() || ca.size() + cb.size() <= nums[1].size() + divs[1].size())) {
995 if(mgcd.representsNonZero(true) || !eo.warn_about_denominators_assumed_nonzero || warn_about_denominators_assumed_nonzero(mgcd, eo)) {
996 if(cb.isOne()) {
997 numleft.addChild(ca);
998 nums.delChild(2);
999 divs.delChild(2);
1000 b = false;
1001 } else {
1002 nums[1] = ca;
1003 divs[1] = cb;
1004 }
1005 }
1006 }
1007 }
1008 if(CALCULATOR->aborted()) {if(mstruct.size() == 1) mstruct.setToChild(1); return false;}
1009 if(b) {
1010 MathStructure ca, cb, mgcd;
1011 b = MathStructure::gcd(divs[0], divs[1], mgcd, eo2, &ca, &cb, false) && !mgcd.isOne();
1012 if(CALCULATOR->aborted()) {if(mstruct.size() == 1) mstruct.setToChild(1); return false;}
1013 bool b_merge = true;
1014 if(b) {
1015 if(limit_size && ((cb.isAddition() && ((divs[0].isAddition() && divs[0].size() * cb.size() > 200) || (nums[0].isAddition() && nums[0].size() * cb.size() > 200))) || (ca.isAddition() && nums[1].isAddition() && nums[1].size() * ca.size() > 200))) {
1016 b_merge = false;
1017 } else {
1018 if(!cb.isOne()) {
1019 divs[0].calculateMultiply(cb, eo2);
1020 nums[0].calculateMultiply(cb, eo2);
1021 }
1022 if(!ca.isOne()) {
1023 nums[1].calculateMultiply(ca, eo2);
1024 }
1025 }
1026 } else {
1027 if(limit_size && ((divs[1].isAddition() && ((divs[0].isAddition() && divs[0].size() * divs[1].size() > 200) || (nums[0].isAddition() && nums[0].size() * divs[1].size() > 200))) || (divs[0].isAddition() && nums[1].isAddition() && nums[1].size() * divs[0].size() > 200))) {
1028 b_merge = false;
1029 } else {
1030 nums[0].calculateMultiply(divs[1], eo2);
1031 nums[1].calculateMultiply(divs[0], eo2);
1032 divs[0].calculateMultiply(divs[1], eo2);
1033 }
1034 }
1035 if(b_merge) {
1036 nums[0].calculateAdd(nums[1], eo2);
1037 nums.delChild(2);
1038 divs.delChild(2);
1039 } else {
1040 size_t size_1 = 2, size_2 = 2;
1041 if(nums[0].isAddition()) size_1 += nums[0].size() - 1;
1042 if(divs[0].isAddition()) size_1 += divs[0].size() - 1;
1043 if(nums[1].isAddition()) size_2 += nums[1].size() - 1;
1044 if(divs[1].isAddition()) size_2 += divs[1].size() - 1;
1045 if(size_1 > size_2) {
1046 nums[0].calculateDivide(divs[0], eo);
1047 mleft.addChild(nums[0]);
1048 nums.delChild(1);
1049 divs.delChild(1);
1050 } else {
1051 nums[1].calculateDivide(divs[1], eo);
1052 mleft.addChild(nums[1]);
1053 nums.delChild(2);
1054 divs.delChild(2);
1055 }
1056 }
1057 }
1058 } else if(b && numleft.size() > 0) {
1059 if(limit_size && divs[0].isAddition() && numleft.size() > 0 && divs[0].size() * numleft.size() > 200) break;
1060 if(numleft.size() == 1) numleft.setToChild(1);
1061 else if(numleft.size() > 1) numleft.setType(STRUCT_ADDITION);
1062 numleft.calculateMultiply(divs[0], eo2);
1063 nums[0].calculateAdd(numleft, eo2);
1064 numleft.clear();
1065 } else if(b) break;
1066 }
1067
1068 if(CALCULATOR->aborted()) {if(mstruct.size() == 1) mstruct.setToChild(1); return false;}
1069 while(!combine_only && !only_gcd && divs.size() > 0 && divs[0].isAddition() && !nums[0].isNumber()) {
1070 bool b_unknown = divs[0].containsUnknowns();
1071 if(!b_unknown || nums[0].containsUnknowns()) {
1072 MathStructure mcomps;
1073 mcomps.resizeVector(nums[0].isAddition() ? nums[0].size() : 1, m_zero);
1074 if(nums[0].isAddition()) {
1075 for(size_t i = 0; i < nums[0].size(); i++) {
1076 if((b_unknown && nums[0][i].containsUnknowns()) || (!b_unknown && !nums[0][i].isNumber())) {
1077 mcomps[i].setType(STRUCT_MULTIPLICATION);
1078 if(nums[0][i].isMultiplication()) {
1079 for(size_t i2 = 0; i2 < nums[0][i].size(); i2++) {
1080 if((b_unknown && nums[0][i][i2].containsUnknowns()) || (!b_unknown && !nums[0][i][i2].isNumber())) {
1081 nums[0][i][i2].ref();
1082 mcomps[i].addChild_nocopy(&nums[0][i][i2]);
1083 }
1084 }
1085 } else {
1086 nums[0][i].ref();
1087 mcomps[i].addChild_nocopy(&nums[0][i]);
1088 }
1089 }
1090 }
1091 } else {
1092 mcomps[0].setType(STRUCT_MULTIPLICATION);
1093 if(nums[0].isMultiplication()) {
1094 for(size_t i2 = 0; i2 < nums[0].size(); i2++) {
1095 if((b_unknown && nums[0][i2].containsUnknowns()) || (!b_unknown && !nums[0][i2].isNumber())) {
1096 nums[0][i2].ref();
1097 mcomps[0].addChild_nocopy(&nums[0][i2]);
1098 }
1099 }
1100 } else {
1101 nums[0].ref();
1102 mcomps[0].addChild_nocopy(&nums[0]);
1103 }
1104 }
1105 bool b_found = false;
1106 unordered_map<size_t, size_t> matches;
1107 for(size_t i = 0; i < divs[0].size(); i++) {
1108 if((b_unknown && divs[0][i].containsUnknowns()) || (!b_unknown && !divs[0][i].isNumber())) {
1109 MathStructure mcomp1;
1110 mcomp1.setType(STRUCT_MULTIPLICATION);
1111 if(divs[0][i].isMultiplication()) {
1112 for(size_t i2 = 0; i2 < divs[0][i].size(); i2++) {
1113 if((b_unknown && divs[0][i][i2].containsUnknowns()) || (!b_unknown && !divs[0][i][i2].isNumber())) {
1114 divs[0][i][i2].ref();
1115 mcomp1.addChild_nocopy(&divs[0][i][i2]);
1116 }
1117 }
1118 } else {
1119 divs[0][i].ref();
1120 mcomp1.addChild_nocopy(&divs[0][i]);
1121 }
1122 b_found = false;
1123 for(size_t i2 = 0; i2 < mcomps.size(); i2++) {
1124 if(mcomps[i2] == mcomp1) {
1125 matches[i] = i2;
1126 b_found = true;
1127 break;
1128 }
1129 }
1130 if(!b_found) break;
1131 }
1132 }
1133 if(b_found) {
1134 b_found = false;
1135 size_t i_selected = 0;
1136 if(nums[0].isAddition()) {
1137 Number i_points;
1138 for(size_t i = 0; i < divs[0].size(); i++) {
1139 if((b_unknown && divs[0][i].containsUnknowns()) || (!b_unknown && !divs[0][i].isNumber())) {
1140 Number new_points = dos_count_points(divs[0][i], b_unknown);
1141 if(new_points > i_points) {i_points = new_points; i_selected = i;}
1142 }
1143 }
1144 }
1145 MathStructure mquo, mrem;
1146 if(polynomial_long_division(nums[0].isAddition() ? nums[0][matches[i_selected]] : nums[0], divs[0], m_zero, mquo, mrem, eo2, false) && !mquo.isZero() && mrem != (nums[0].isAddition() ? nums[0][matches[i_selected]] : nums[0])) {
1147 if(CALCULATOR->aborted()) return false;
1148 if((nums[0].isAddition() && nums[0].size() > 1) || !mrem.isZero() || divs[0].representsNonZero(true) || (eo.warn_about_denominators_assumed_nonzero && !warn_about_denominators_assumed_nonzero(divs[0], eo))) {
1149 mleft.addChild(mquo);
1150 if(nums[0].isAddition()) {
1151 nums[0].delChild(matches[i_selected] + 1, true);
1152 if(!mrem.isZero()) {
1153 nums[0].calculateAdd(mrem, eo2);
1154 }
1155 } else {
1156 if(mrem.isZero()) {
1157 divs.clear();
1158 } else {
1159 nums[0] = mrem;
1160 }
1161 }
1162 b_ret = true;
1163 b_found = true;
1164 }
1165 }
1166 }
1167 if(!b_found) break;
1168 } else {
1169 break;
1170 }
1171 }
1172
1173 if(CALCULATOR->aborted()) {if(mstruct.size() == 1) mstruct.setToChild(1); return false;}
1174 if(!combine_only && !only_gcd && divs.size() > 0 && divs[0].isAddition()) {
1175 bool b_unknown = divs[0].containsUnknowns() || nums[0].containsUnknowns();
1176 MathStructure mquo, mrem;
1177 vector<MathStructure> symsd, symsn;
1178 collect_symbols(nums[0], symsn);
1179 if(!symsn.empty()) {
1180 collect_symbols(divs[0], symsd);
1181 }
1182 for(size_t i = 0; i < symsd.size();) {
1183 bool b_found = false;
1184 if(!b_unknown || symsd[i].containsUnknowns()) {
1185 for(size_t i2 = 0; i2 < symsn.size(); i2++) {
1186 if(symsn[i2] == symsd[i]) {
1187 b_found = (nums[0].degree(symsd[i]) >= divs[0].degree(symsd[i]));
1188 break;
1189 }
1190 }
1191 }
1192 if(b_found) i++;
1193 else symsd.erase(symsd.begin() + i);
1194 }
1195 for(size_t i = 0; i < symsd.size(); i++) {
1196 if(polynomial_long_division(nums[0], divs[0], symsd[i], mquo, mrem, eo2, false) && !mquo.isZero() && mrem != nums[0]) {
1197 if(CALCULATOR->aborted()) {if(mstruct.size() == 1) mstruct.setToChild(1); return false;}
1198 if(!mrem.isZero() || divs[0].representsNonZero(true) || (eo.warn_about_denominators_assumed_nonzero && !warn_about_denominators_assumed_nonzero(divs[0], eo))) {
1199 if(mrem.isZero()) {
1200 mleft.addChild(mquo);
1201 divs.clear();
1202 b_ret = true;
1203 break;
1204 } else {
1205 if(mquo.isAddition()) {
1206 for(size_t i = 0; i < mquo.size(); ) {
1207 MathStructure mtest(mquo[i]);
1208 mtest.calculateMultiply(divs[0], eo2);
1209 mtest.calculateAdd(mrem, eo2);
1210 Number point = dos_count_points(mtest, b_unknown);
1211 point -= dos_count_points(mquo[i], b_unknown);
1212 point -= dos_count_points(mrem, b_unknown);
1213 if(point <= 0) {
1214 mquo.delChild(i + 1);
1215 mrem = mtest;
1216 } else {
1217 i++;
1218 }
1219 }
1220 if(mquo.size() == 0) mquo.clear();
1221 else if(mquo.size() == 1) mquo.setToChild(1);
1222 }
1223 if(!mquo.isZero()) {
1224 Number point = dos_count_points(nums[0], b_unknown);
1225 point -= dos_count_points(mquo, b_unknown);
1226 point -= dos_count_points(mrem, b_unknown);
1227 if(point > 0) {
1228 mleft.addChild(mquo);
1229 MathStructure ca, cb, mgcd;
1230 if(mrem.isAddition() && MathStructure::gcd(mrem, divs[0], mgcd, eo2, &ca, &cb, false) && !mgcd.isOne() && (!cb.isAddition() || !ca.isAddition() || ca.size() + cb.size() <= mrem.size() + divs[0].size())) {
1231 mrem = ca;
1232 divs[0] = cb;
1233 }
1234 mrem.calculateDivide(divs[0], eo2);
1235 mleft.addChild(mrem);
1236 divs.clear();
1237 b_ret = true;
1238 break;
1239 }
1240 }
1241 }
1242 }
1243 }
1244 }
1245 }
1246 }
1247
1248 if(!b_ret) {if(mstruct.size() == 1) mstruct.setToChild(1); return false;}
1249
1250 mstruct.clear(true);
1251 for(size_t i = 0; i < divs.size() && (i == 0 || i_run > 1); i++) {
1252 if(i == 0) mstruct = nums[0];
1253 else mstruct.add(nums[i], i > 1);
1254 if(only_gcd || combine_only) {
1255 divs[i].inverse();
1256 if(i == 0) mstruct.multiply(divs[i], true);
1257 else mstruct.last().multiply(divs[i], true);
1258 } else {
1259 if(i == 0) mstruct.calculateDivide(divs[i], eo);
1260 else mstruct.last().calculateDivide(divs[i], eo);
1261 if(i == 0 && i_run <= 1 && mstruct.isAddition()) {
1262 do_simplification(mstruct, eo, true, false, false, false, limit_size, 2);
1263 }
1264 }
1265 }
1266 for(size_t i = 0; i < mleft.size(); i++) {
1267 if(i == 0 && mstruct.isZero()) mstruct = mleft[i];
1268 else mstruct.calculateAdd(mleft[i], eo);
1269 }
1270 for(size_t i = 0; i < numleft.size(); i++) {
1271 if(i == 0 && mstruct.isZero()) mstruct = numleft[i];
1272 else mstruct.calculateAdd(numleft[i], eo);
1273 }
1274
1275 return true;
1276 }
1277
1278 MathStructure divs;
1279 // find division by polynomial
1280 for(size_t i = 0; i < mstruct.size(); i++) {
1281 if(mstruct[i].isMultiplication()) {
1282 for(size_t i2 = 0; i2 < mstruct[i].size(); i2++) {
1283 if(mstruct[i][i2].isPower() && (mstruct[i][i2][0].isAddition() || combine_divisions) && mstruct[i][i2][1].isMinusOne() && mstruct[i][i2][0].isRationalPolynomial()) {
1284 bool b_found = false;
1285 for(size_t i3 = 0; i3 < divs.size(); i3++) {
1286 if(divs[i3] == mstruct[i][i2][0]) {
1287 divs[i3].number()++;
1288 b_found = true;
1289 break;
1290 }
1291 }
1292 if(!b_found) divs.addChild(mstruct[i][i2][0]);
1293 break;
1294 }
1295 }
1296 } else if(mstruct[i].isPower() && (mstruct[i][0].isAddition() || combine_divisions) && mstruct[i][1].isMinusOne() && mstruct[i][0].isRationalPolynomial()) {
1297 bool b_found = false;
1298 for(size_t i3 = 0; i3 < divs.size(); i3++) {
1299 if(divs[i3] == mstruct[i][0]) {
1300 divs[i3].number()++;
1301 b_found = true;
1302 break;
1303 }
1304 }
1305 if(!b_found) divs.addChild(mstruct[i][0]);
1306 }
1307 }
1308
1309 // check if denominators is zero
1310 for(size_t i = 0; i < divs.size(); ) {
1311 if((!combine_divisions && divs[i].number().isZero()) || (!eo.assume_denominators_nonzero && !divs[i].representsNonZero(true)) || divs[i].representsZero(true)) divs.delChild(i + 1);
1312 else i++;
1313 }
1314
1315 if(divs.size() == 0) return false;
1316
1317 // combine numerators with same denominator
1318 MathStructure nums, numleft, mleft;
1319 nums.resizeVector(divs.size(), m_zero);
1320 for(size_t i = 0; i < mstruct.size(); i++) {
1321 bool b = false;
1322 if(mstruct[i].isMultiplication()) {
1323 for(size_t i2 = 0; i2 < mstruct[i].size(); i2++) {
1324 if(mstruct[i][i2].isPower() && (mstruct[i][i2][0].isAddition() || combine_divisions) && mstruct[i][i2][1].isMinusOne() && mstruct[i][i2][0].isRationalPolynomial()) {
1325 for(size_t i3 = 0; i3 < divs.size(); i3++) {
1326 if(divs[i3] == mstruct[i][i2][0]) {
1327 if(mstruct[i].size() == 1) nums[i3].addChild(m_one);
1328 else if(mstruct[i].size() == 2) nums[i3].addChild(mstruct[i][i2 == 0 ? 1 : 0]);
1329 else {nums[i3].addChild(mstruct[i]); nums[i3][nums[i3].size() - 1].delChild(i2 + 1);}
1330 b = true;
1331 break;
1332 }
1333 }
1334 break;
1335 }
1336 }
1337 } else if(mstruct[i].isPower() && (mstruct[i][0].isAddition() || combine_divisions) && mstruct[i][1].isMinusOne() && mstruct[i][0].isRationalPolynomial()) {
1338 for(size_t i3 = 0; i3 < divs.size(); i3++) {
1339 if(divs[i3] == mstruct[i][0]) {
1340 nums[i3].addChild(m_one);
1341 b = true;
1342 break;
1343 }
1344 }
1345 }
1346 if(!b && combine_divisions) {
1347 if(mstruct[i].isRationalPolynomial()) numleft.addChild(mstruct[i]);
1348 else mleft.addChild(mstruct[i]);
1349 }
1350 }
1351
1352 EvaluationOptions eo2 = eo;
1353 eo2.do_polynomial_division = false;
1354 eo2.keep_zero_units = false;
1355
1356 // do polynomial division; give points to incomplete division
1357 vector<long int> points(nums.size(), -1);
1358 bool b = false, b_ready_candidate = false;
1359 for(size_t i = 0; i < divs.size(); i++) {
1360 if(CALCULATOR->aborted()) return false;
1361 if(nums[i].size() > 1 && divs[i].isAddition()) {
1362 nums[i].setType(STRUCT_ADDITION);
1363 MathStructure xvar;
1364 get_first_symbol(nums[i], xvar);
1365 if(nums[i].isRationalPolynomial() && nums[i].degree(xvar).isLessThan(100) && divs[i].degree(xvar).isLessThan(100)) {
1366 MathStructure mquo, mrem, ca, cb;
1367 if(MathStructure::gcd(nums[i], divs[i], mquo, eo2, &ca, &cb, false) && !mquo.isOne()) {
1368 if(!mquo.representsNonZero(true) && eo.warn_about_denominators_assumed_nonzero && !warn_about_denominators_assumed_nonzero(mquo, eo)) {
1369 nums[i].clear();
1370 } else {
1371 points[i] = 1;
1372 b_ready_candidate = true;
1373 b = true;
1374 }
1375 if(ca.isOne()) {
1376 if(cb.isOne()) {
1377 nums[i].set(1, 1, 0, true);
1378 } else {
1379 nums[i] = cb;
1380 nums[i].inverse();
1381 }
1382 } else if(cb.isOne()) {
1383 nums[i] = ca;
1384 } else {
1385 nums[i] = ca;
1386 nums[i].calculateDivide(cb, eo);
1387 }
1388 } else if(!only_gcd && polynomial_long_division(nums[i], divs[i], xvar, mquo, mrem, eo2, false) && !mquo.isZero() && mrem != nums[i]) {
1389 if(mrem.isZero() && !divs[i].representsNonZero(true) && eo.warn_about_denominators_assumed_nonzero && !warn_about_denominators_assumed_nonzero(divs[i], eo)) {
1390 nums[i].clear();
1391 } else {
1392 long int point = 1;
1393 if(!mrem.isZero()) point = nums[i].size();
1394 nums[i].set(mquo);
1395 if(!mrem.isZero()) {
1396 if(mquo.isAddition()) point -= mquo.size();
1397 else point--;
1398 if(mrem.isAddition()) point -= mrem.size();
1399 else point--;
1400 mrem.calculateDivide(divs[i], eo2);
1401 nums[i].calculateAdd(mrem, eo2);
1402 }
1403 b = true;
1404 points[i] = point;
1405 if(point >= 0) {b_ready_candidate = true;}
1406 else if(b_ready_candidate) nums[i].clear();
1407 }
1408 } else {
1409 nums[i].clear();
1410 }
1411 }
1412 } else {
1413 nums[i].clear();
1414 }
1415 }
1416
1417 if(!b) return false;
1418
1419 if(b_ready_candidate) {
1420 // remove polynomial divisions that inrease complexity
1421 for(size_t i = 0; i < nums.size(); i++) {
1422 if(!nums[i].isZero() && points[i] < 0) nums[i].clear();
1423 }
1424 } else {
1425 // no simplying polynomial division found; see if result can be combined with other terms
1426 b = false;
1427 for(size_t i = 0; i < nums.size(); i++) {
1428 if(!nums[i].isZero()) {
1429 if(b) {
1430 nums[i].clear();
1431 } else {
1432 long int point = points[i];
1433 for(size_t i2 = 0; i2 < nums[i].size(); i2++) {
1434 bool b2 = false;
1435 if(!nums[i][i2].contains(divs[i], false, false, false)) {
1436 MathStructure mtest1(mstruct), mtest2(nums[i][i2]);
1437 for(size_t i3 = 0; i3 < mtest1.size(); i3++) {
1438 if(!mtest1[i3].contains(divs[i], false, false, false)) {
1439 int ret = mtest1[i3].merge_addition(mtest2, eo);
1440 if(ret > 0) {
1441 b2 = true;
1442 point++;
1443 if(mtest1[i3].isZero()) point++;
1444 break;
1445 }
1446 if(ret == 0) ret = mtest2.merge_addition(mtest1[i3], eo);
1447 if(ret > 0) {
1448 b2 = true;
1449 point++;
1450 if(mtest2.isZero()) point++;
1451 break;
1452 }
1453 }
1454 }
1455 if(b2) break;
1456 }
1457 if(point >= 0) break;
1458 }
1459 if(point >= 0) b = true;
1460 else nums[i].clear();
1461 }
1462 }
1463 }
1464 }
1465
1466 if(!b) return false;
1467 // replace terms with polynomial division result
1468 for(size_t i = 0; i < mstruct.size(); ) {
1469 bool b_del = false;
1470 if(mstruct[i].isMultiplication()) {
1471 for(size_t i2 = 0; i2 < mstruct[i].size(); i2++) {
1472 if(mstruct[i][i2].isPower() && mstruct[i][i2][0].isAddition() && mstruct[i][i2][1].isMinusOne() && mstruct[i][i2][0].isRationalPolynomial()) {
1473 for(size_t i3 = 0; i3 < divs.size(); i3++) {
1474 if(divs[i3] == mstruct[i][i2][0]) {
1475 b_del = !nums[i3].isZero();
1476 break;
1477 }
1478 }
1479 break;
1480 }
1481 }
1482 } else if(mstruct[i].isPower() && mstruct[i][0].isAddition() && mstruct[i][1].isMinusOne() && mstruct[i][0].isRationalPolynomial()) {
1483 for(size_t i3 = 0; i3 < divs.size(); i3++) {
1484 if(divs[i3] == mstruct[i][0]) {
1485 b_del = !nums[i3].isZero();
1486 break;
1487 }
1488 }
1489 }
1490 if(b_del) mstruct.delChild(i + 1);
1491 else i++;
1492 }
1493 for(size_t i = 0; i < nums.size(); ) {
1494 if(nums[i].isZero()) {
1495 nums.delChild(i + 1);
1496 } else {
1497 nums[i].evalSort();
1498 i++;
1499 }
1500 }
1501 if(mstruct.size() == 0 && nums.size() == 1) {
1502 mstruct.set(nums[0]);
1503 } else {
1504 for(size_t i = 0; i < nums.size(); i++) {
1505 mstruct.addChild(nums[i]);
1506 }
1507 mstruct.evalSort();
1508 }
1509 return true;
1510 }
1511
1512