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 #include <algorithm>
23 #include "MathStructure-support.h"
24
25 using std::string;
26 using std::cout;
27 using std::vector;
28 using std::endl;
29
operator <(const sym_desc & x) const30 bool sym_desc::operator<(const sym_desc &x) const {
31 if (max_deg == x.max_deg) return max_lcnops < x.max_lcnops;
32 else return max_deg.isLessThan(x.max_deg);
33 }
34
polynomialUnit(const MathStructure & xvar) const35 int MathStructure::polynomialUnit(const MathStructure &xvar) const {
36 MathStructure coeff;
37 lcoefficient(xvar, coeff);
38 if(coeff.hasNegativeSign()) return -1;
39 return 1;
40 }
41
integer_content(const MathStructure & mpoly,Number & icontent)42 void integer_content(const MathStructure &mpoly, Number &icontent) {
43 if(mpoly.isNumber()) {
44 icontent = mpoly.number();
45 icontent.abs();
46 } else if(mpoly.isAddition()) {
47 icontent.clear();
48 Number l(1, 1);
49 for(size_t i = 0; i < mpoly.size(); i++) {
50 if(mpoly[i].isNumber()) {
51 if(!icontent.isOne()) {
52 Number c = icontent;
53 icontent = mpoly[i].number().numerator();
54 icontent.gcd(c);
55 }
56 Number l2 = l;
57 l = mpoly[i].number().denominator();
58 l.lcm(l2);
59 } else if(mpoly[i].isMultiplication()) {
60 if(!icontent.isOne()) {
61 Number c = icontent;
62 icontent = mpoly[i].overallCoefficient().numerator();
63 icontent.gcd(c);
64 }
65 Number l2 = l;
66 l = mpoly[i].overallCoefficient().denominator();
67 l.lcm(l2);
68 } else {
69 icontent.set(1, 1, 0);
70 }
71 }
72 icontent /= l;
73 } else if(mpoly.isMultiplication()) {
74 icontent = mpoly.overallCoefficient();
75 icontent.abs();
76 } else {
77 icontent.set(1, 1, 0);
78 }
79 }
polynomialContent(const MathStructure & xvar,MathStructure & mcontent,const EvaluationOptions & eo) const80 void MathStructure::polynomialContent(const MathStructure &xvar, MathStructure &mcontent, const EvaluationOptions &eo) const {
81 if(isZero()) {
82 mcontent.clear();
83 return;
84 }
85 if(isNumber()) {
86 mcontent = *this;
87 mcontent.number().setNegative(false);
88 return;
89 }
90
91 MathStructure c;
92 integer_content(*this, c.number());
93 MathStructure r(*this);
94 if(!c.isOne()) r.calculateDivide(c, eo);
95 MathStructure lcoeff;
96 r.lcoefficient(xvar, lcoeff);
97 if(lcoeff.isInteger()) {
98 mcontent = c;
99 return;
100 }
101 Number deg(r.degree(xvar));
102 Number ldeg(r.ldegree(xvar));
103 if(deg == ldeg) {
104 mcontent = lcoeff;
105 if(lcoeff.polynomialUnit(xvar) == -1) {
106 c.number().negate();
107 }
108 mcontent.calculateMultiply(c, eo);
109 return;
110 }
111 mcontent.clear();
112 MathStructure mtmp, coeff;
113 for(Number i(ldeg); i.isLessThanOrEqualTo(deg); i++) {
114 coefficient(xvar, i, coeff);
115 mtmp = mcontent;
116 if(!MathStructure::gcd(coeff, mtmp, mcontent, eo, NULL, NULL, false)) mcontent.set(1, 1, 0);
117 if(mcontent.isOne()) break;
118 }
119
120 if(!c.isOne()) mcontent.calculateMultiply(c, eo);
121
122 }
123
polynomialPrimpart(const MathStructure & xvar,MathStructure & mprim,const EvaluationOptions & eo) const124 void MathStructure::polynomialPrimpart(const MathStructure &xvar, MathStructure &mprim, const EvaluationOptions &eo) const {
125 if(isZero()) {
126 mprim.clear();
127 return;
128 }
129 if(isNumber()) {
130 mprim.set(1, 1, 0);
131 return;
132 }
133
134 MathStructure c;
135 polynomialContent(xvar, c, eo);
136 if(c.isZero()) {
137 mprim.clear();
138 return;
139 }
140 bool b = (polynomialUnit(xvar) == -1);
141 if(c.isNumber()) {
142 if(b) c.number().negate();
143 mprim = *this;
144 mprim.calculateDivide(c, eo);
145 return;
146 }
147 if(b) c.calculateNegate(eo);
148 MathStructure::polynomialQuotient(*this, c, xvar, mprim, eo, false);
149 }
polynomialUnitContentPrimpart(const MathStructure & xvar,int & munit,MathStructure & mcontent,MathStructure & mprim,const EvaluationOptions & eo) const150 void MathStructure::polynomialUnitContentPrimpart(const MathStructure &xvar, int &munit, MathStructure &mcontent, MathStructure &mprim, const EvaluationOptions &eo) const {
151
152 if(isZero()) {
153 munit = 1;
154 mcontent.clear();
155 mprim.clear();
156 return;
157 }
158
159 if(isNumber()) {
160 if(o_number.isNegative()) {
161 munit = -1;
162 mcontent = *this;
163 mcontent.number().abs();
164 } else {
165 munit = 1;
166 mcontent = *this;
167 }
168 mprim.set(1, 1, 0);
169 return;
170 }
171
172 munit = polynomialUnit(xvar);
173
174 polynomialContent(xvar, mcontent, eo);
175
176 if(mcontent.isZero()) {
177 mprim.clear();
178 return;
179 }
180 if(mcontent.isNumber()) {
181 mprim = *this;
182 if(munit == -1) {
183 Number c(mcontent.number());
184 c.negate();
185 mprim.calculateDivide(c, eo);
186 } else if(!mcontent.isOne()) {
187 mprim.calculateDivide(mcontent, eo);
188 }
189 } else {
190 if(munit == -1) {
191 MathStructure c(mcontent);
192 c.calculateNegate(eo);
193 MathStructure::polynomialQuotient(*this, c, xvar, mprim, eo, false);
194 } else {
195 MathStructure::polynomialQuotient(*this, mcontent, xvar, mprim, eo, false);
196 }
197 }
198 }
199
200
polynomialPrimpart(const MathStructure & xvar,const MathStructure & c,MathStructure & mprim,const EvaluationOptions & eo) const201 void MathStructure::polynomialPrimpart(const MathStructure &xvar, const MathStructure &c, MathStructure &mprim, const EvaluationOptions &eo) const {
202 if(isZero() || c.isZero()) {
203 mprim.clear();
204 return;
205 }
206 if(isNumber()) {
207 mprim.set(1, 1, 0);
208 return;
209 }
210 bool b = (polynomialUnit(xvar) == -1);
211 if(c.isNumber()) {
212 MathStructure cn(c);
213 if(b) cn.number().negate();
214 mprim = *this;
215 mprim.calculateDivide(cn, eo);
216 return;
217 }
218 if(b) {
219 MathStructure cn(c);
220 cn.calculateNegate(eo);
221 MathStructure::polynomialQuotient(*this, cn, xvar, mprim, eo, false);
222 } else {
223 MathStructure::polynomialQuotient(*this, c, xvar, mprim, eo, false);
224 }
225 }
226
degree(const MathStructure & xvar) const227 const Number& MathStructure::degree(const MathStructure &xvar) const {
228 const Number *c = NULL;
229 const MathStructure *mcur = NULL;
230 for(size_t i = 0; ; i++) {
231 if(isAddition()) {
232 if(i >= SIZE) break;
233 mcur = &CHILD(i);
234 } else {
235 mcur = this;
236 }
237 if((*mcur) == xvar) {
238 if(!c) {
239 c = &nr_one;
240 }
241 } else if(mcur->isPower() && (*mcur)[0] == xvar && (*mcur)[1].isNumber()) {
242 if(!c || c->isLessThan((*mcur)[1].number())) {
243 c = &(*mcur)[1].number();
244 }
245 } else if(mcur->isMultiplication()) {
246 for(size_t i2 = 0; i2 < mcur->size(); i2++) {
247 if((*mcur)[i2] == xvar) {
248 if(!c) {
249 c = &nr_one;
250 }
251 } else if((*mcur)[i2].isPower() && (*mcur)[i2][0] == xvar && (*mcur)[i2][1].isNumber()) {
252 if(!c || c->isLessThan((*mcur)[i2][1].number())) {
253 c = &(*mcur)[i2][1].number();
254 }
255 }
256 }
257 }
258 if(!isAddition()) break;
259 }
260 if(!c) return nr_zero;
261 return *c;
262 }
ldegree(const MathStructure & xvar) const263 const Number& MathStructure::ldegree(const MathStructure &xvar) const {
264 const Number *c = NULL;
265 const MathStructure *mcur = NULL;
266 for(size_t i = 0; ; i++) {
267 if(isAddition()) {
268 if(i >= SIZE) break;
269 mcur = &CHILD(i);
270 } else {
271 mcur = this;
272 }
273 if((*mcur) == xvar) {
274 c = &nr_one;
275 } else if(mcur->isPower() && (*mcur)[0] == xvar && (*mcur)[1].isNumber()) {
276 if(!c || c->isGreaterThan((*mcur)[1].number())) {
277 c = &(*mcur)[1].number();
278 }
279 } else if(mcur->isMultiplication()) {
280 bool b = false;
281 for(size_t i2 = 0; i2 < mcur->size(); i2++) {
282 if((*mcur)[i2] == xvar) {
283 c = &nr_one;
284 b = true;
285 } else if((*mcur)[i2].isPower() && (*mcur)[i2][0] == xvar && (*mcur)[i2][1].isNumber()) {
286 if(!c || c->isGreaterThan((*mcur)[i2][1].number())) {
287 c = &(*mcur)[i2][1].number();
288 }
289 b = true;
290 }
291 }
292 if(!b) return nr_zero;
293 } else {
294 return nr_zero;
295 }
296 if(!isAddition()) break;
297 }
298 if(!c) return nr_zero;
299 return *c;
300 }
lcoefficient(const MathStructure & xvar,MathStructure & mcoeff) const301 void MathStructure::lcoefficient(const MathStructure &xvar, MathStructure &mcoeff) const {
302 coefficient(xvar, degree(xvar), mcoeff);
303 }
tcoefficient(const MathStructure & xvar,MathStructure & mcoeff) const304 void MathStructure::tcoefficient(const MathStructure &xvar, MathStructure &mcoeff) const {
305 coefficient(xvar, ldegree(xvar), mcoeff);
306 }
coefficient(const MathStructure & xvar,const Number & pownr,MathStructure & mcoeff) const307 void MathStructure::coefficient(const MathStructure &xvar, const Number &pownr, MathStructure &mcoeff) const {
308 const MathStructure *mcur = NULL;
309 mcoeff.clear();
310 for(size_t i = 0; ; i++) {
311 if(isAddition()) {
312 if(i >= SIZE) break;
313 mcur = &CHILD(i);
314 } else {
315 mcur = this;
316 }
317 if((*mcur) == xvar) {
318 if(pownr.isOne()) {
319 if(mcoeff.isZero()) mcoeff.set(1, 1, 0);
320 else mcoeff.add(m_one, true);
321 }
322 } else if(mcur->isPower() && (*mcur)[0] == xvar && (*mcur)[1].isNumber()) {
323 if((*mcur)[1].number() == pownr) {
324 if(mcoeff.isZero()) mcoeff.set(1, 1, 0);
325 else mcoeff.add(m_one, true);
326 }
327 } else if(mcur->isMultiplication()) {
328 bool b = false;
329 for(size_t i2 = 0; i2 < mcur->size(); i2++) {
330 bool b2 = false;
331 if((*mcur)[i2] == xvar) {
332 b = true;
333 if(pownr.isOne()) b2 = true;
334 } else if((*mcur)[i2].isPower() && (*mcur)[i2][0] == xvar && (*mcur)[i2][1].isNumber()) {
335 b = true;
336 if((*mcur)[i2][1].number() == pownr) b2 = true;
337 }
338 if(b2) {
339 if(mcoeff.isZero()) {
340 if(mcur->size() == 1) {
341 mcoeff.set(1, 1, 0);
342 } else {
343 for(size_t i3 = 0; i3 < mcur->size(); i3++) {
344 if(i3 != i2) {
345 if(mcoeff.isZero()) mcoeff = (*mcur)[i3];
346 else mcoeff.multiply((*mcur)[i3], true);
347 }
348 }
349 }
350 } else if(mcur->size() == 1) {
351 mcoeff.add(m_one, true);
352 } else {
353 mcoeff.add(m_zero, true);
354 for(size_t i3 = 0; i3 < mcur->size(); i3++) {
355 if(i3 != i2) {
356 if(mcoeff[mcoeff.size() - 1].isZero()) mcoeff[mcoeff.size() - 1] = (*mcur)[i3];
357 else mcoeff[mcoeff.size() - 1].multiply((*mcur)[i3], true);
358 }
359 }
360 }
361 break;
362 }
363 }
364 if(!b && pownr.isZero()) {
365 if(mcoeff.isZero()) mcoeff = *mcur;
366 else mcoeff.add(*mcur, true);
367 }
368 } else if(pownr.isZero()) {
369 if(mcoeff.isZero()) mcoeff = *mcur;
370 else mcoeff.add(*mcur, true);
371 }
372 if(!isAddition()) break;
373 }
374 mcoeff.evalSort();
375 }
376
polynomialQuotient(const MathStructure & mnum,const MathStructure & mden,const MathStructure & xvar,MathStructure & mquotient,const EvaluationOptions & eo,bool check_args)377 bool MathStructure::polynomialQuotient(const MathStructure &mnum, const MathStructure &mden, const MathStructure &xvar, MathStructure &mquotient, const EvaluationOptions &eo, bool check_args) {
378 mquotient.clear();
379 if(CALCULATOR->aborted()) return false;
380 if(mden.isZero()) {
381 //division by zero
382 return false;
383 }
384 if(mnum.isZero()) {
385 mquotient.clear();
386 return true;
387 }
388 if(mden.isNumber() && mnum.isNumber()) {
389 mquotient = mnum;
390 if(IS_REAL(mden) && IS_REAL(mnum)) {
391 mquotient.number() /= mden.number();
392 } else {
393 mquotient.calculateDivide(mden, eo);
394 }
395 return true;
396 }
397 if(mnum == mden) {
398 mquotient.set(1, 1, 0);
399 return true;
400 }
401 if(check_args && (!mnum.isRationalPolynomial() || !mden.isRationalPolynomial())) {
402 return false;
403 }
404
405 EvaluationOptions eo2 = eo;
406 eo2.keep_zero_units = false;
407
408 Number numdeg = mnum.degree(xvar);
409 Number dendeg = mden.degree(xvar);
410 MathStructure dencoeff;
411 mden.coefficient(xvar, dendeg, dencoeff);
412 MathStructure mrem(mnum);
413 for(size_t i = 0; numdeg.isGreaterThanOrEqualTo(dendeg); i++) {
414 if(i > 1000 || CALCULATOR->aborted()) {
415 return false;
416 }
417 MathStructure numcoeff;
418 mrem.coefficient(xvar, numdeg, numcoeff);
419 numdeg -= dendeg;
420 if(numcoeff == dencoeff) {
421 if(numdeg.isZero()) {
422 numcoeff.set(1, 1, 0);
423 } else {
424 numcoeff = xvar;
425 if(!numdeg.isOne()) {
426 numcoeff.raise(numdeg);
427 }
428 }
429 } else {
430 if(dencoeff.isNumber()) {
431 if(numcoeff.isNumber()) {
432 numcoeff.number() /= dencoeff.number();
433 } else {
434 numcoeff.calculateDivide(dencoeff, eo2);
435 }
436 } else {
437 MathStructure mcopy(numcoeff);
438 if(!MathStructure::polynomialDivide(mcopy, dencoeff, numcoeff, eo2, false)) {
439 return false;
440 }
441 }
442 if(!numdeg.isZero() && !numcoeff.isZero()) {
443 if(numcoeff.isOne()) {
444 numcoeff = xvar;
445 if(!numdeg.isOne()) {
446 numcoeff.raise(numdeg);
447 }
448 } else {
449 numcoeff.multiply(xvar, true);
450 if(!numdeg.isOne()) {
451 numcoeff[numcoeff.size() - 1].raise(numdeg);
452 }
453 numcoeff.calculateMultiplyLast(eo2);
454 }
455 }
456 }
457 if(mquotient.isZero()) mquotient = numcoeff;
458 else mquotient.calculateAdd(numcoeff, eo2);
459 numcoeff.calculateMultiply(mden, eo2);
460 mrem.calculateSubtract(numcoeff, eo2);
461 if(mrem.isZero()) break;
462 numdeg = mrem.degree(xvar);
463 }
464 return true;
465
466 }
add_symbol(const MathStructure & mpoly,sym_desc_vec & v)467 void add_symbol(const MathStructure &mpoly, sym_desc_vec &v) {
468 sym_desc_vec::const_iterator it = v.begin(), itend = v.end();
469 while (it != itend) {
470 if(it->sym == mpoly) return;
471 ++it;
472 }
473 sym_desc d;
474 d.sym = mpoly;
475 v.push_back(d);
476 }
477
collect_symbols(const MathStructure & mpoly,sym_desc_vec & v)478 void collect_symbols(const MathStructure &mpoly, sym_desc_vec &v) {
479 if(IS_A_SYMBOL(mpoly) || mpoly.isUnit()) {
480 add_symbol(mpoly, v);
481 } else if(mpoly.isAddition() || mpoly.isMultiplication()) {
482 for(size_t i = 0; i < mpoly.size(); i++) collect_symbols(mpoly[i], v);
483 } else if(mpoly.isPower()) {
484 collect_symbols(mpoly[0], v);
485 }
486 }
487
add_symbol(const MathStructure & mpoly,vector<MathStructure> & v)488 void add_symbol(const MathStructure &mpoly, vector<MathStructure> &v) {
489 vector<MathStructure>::const_iterator it = v.begin(), itend = v.end();
490 while (it != itend) {
491 if(*it == mpoly) return;
492 ++it;
493 }
494 v.push_back(mpoly);
495 }
496
collect_symbols(const MathStructure & mpoly,vector<MathStructure> & v)497 void collect_symbols(const MathStructure &mpoly, vector<MathStructure> &v) {
498 if(IS_A_SYMBOL(mpoly)) {
499 add_symbol(mpoly, v);
500 } else if(mpoly.isAddition() || mpoly.isMultiplication()) {
501 for(size_t i = 0; i < mpoly.size(); i++) collect_symbols(mpoly[i], v);
502 } else if(mpoly.isPower()) {
503 collect_symbols(mpoly[0], v);
504 }
505 }
506
get_symbol_stats(const MathStructure & m1,const MathStructure & m2,sym_desc_vec & v)507 void get_symbol_stats(const MathStructure &m1, const MathStructure &m2, sym_desc_vec &v) {
508 collect_symbols(m1, v);
509 collect_symbols(m2, v);
510 sym_desc_vec::iterator it = v.begin(), itend = v.end();
511 while(it != itend) {
512 it->deg_a = m1.degree(it->sym);
513 it->deg_b = m2.degree(it->sym);
514 if(it->deg_a.isGreaterThan(it->deg_b)) {
515 it->max_deg = it->deg_a;
516 } else {
517 it->max_deg = it->deg_b;
518 }
519 it->ldeg_a = m1.ldegree(it->sym);
520 it->ldeg_b = m2.ldegree(it->sym);
521 MathStructure mcoeff;
522 m1.lcoefficient(it->sym, mcoeff);
523 it->max_lcnops = mcoeff.size();
524 m2.lcoefficient(it->sym, mcoeff);
525 if(mcoeff.size() > it->max_lcnops) it->max_lcnops = mcoeff.size();
526 ++it;
527 }
528 std::sort(v.begin(), v.end());
529
530 }
531
get_first_symbol(const MathStructure & mpoly,MathStructure & xvar)532 bool get_first_symbol(const MathStructure &mpoly, MathStructure &xvar) {
533 if(IS_A_SYMBOL(mpoly) || mpoly.isUnit()) {
534 xvar = mpoly;
535 return true;
536 } else if(mpoly.isAddition() || mpoly.isMultiplication()) {
537 for(size_t i = 0; i < mpoly.size(); i++) {
538 if(get_first_symbol(mpoly[i], xvar)) return true;
539 }
540 } else if(mpoly.isPower()) {
541 return get_first_symbol(mpoly[0], xvar);
542 }
543 return false;
544 }
545
polynomialDivide(const MathStructure & mnum,const MathStructure & mden,MathStructure & mquotient,const EvaluationOptions & eo,bool check_args)546 bool MathStructure::polynomialDivide(const MathStructure &mnum, const MathStructure &mden, MathStructure &mquotient, const EvaluationOptions &eo, bool check_args) {
547
548 mquotient.clear();
549
550 if(CALCULATOR->aborted()) return false;
551
552 if(mden.isZero()) {
553 //division by zero
554 return false;
555 }
556 if(mnum.isZero()) {
557 mquotient.clear();
558 return true;
559 }
560 if(mden.isNumber()) {
561 mquotient = mnum;
562 if(mnum.isNumber()) {
563 mquotient.number() /= mden.number();
564 } else {
565 mquotient.calculateDivide(mden, eo);
566 }
567 return true;
568 } else if(mnum.isNumber()) {
569 return false;
570 }
571
572 if(mnum == mden) {
573 mquotient.set(1, 1, 0);
574 return true;
575 }
576
577 if(check_args && (!mnum.isRationalPolynomial() || !mden.isRationalPolynomial())) {
578 return false;
579 }
580
581
582 MathStructure xvar;
583 if(!get_first_symbol(mnum, xvar) && !get_first_symbol(mden, xvar)) return false;
584
585 EvaluationOptions eo2 = eo;
586 eo2.keep_zero_units = false;
587
588 Number numdeg = mnum.degree(xvar);
589 Number dendeg = mden.degree(xvar);
590 MathStructure dencoeff;
591 mden.coefficient(xvar, dendeg, dencoeff);
592
593 MathStructure mrem(mnum);
594
595 for(size_t i = 0; numdeg.isGreaterThanOrEqualTo(dendeg); i++) {
596 if(i > 1000 || CALCULATOR->aborted()) {
597 return false;
598 }
599 MathStructure numcoeff;
600 mrem.coefficient(xvar, numdeg, numcoeff);
601 numdeg -= dendeg;
602 if(numcoeff == dencoeff) {
603 if(numdeg.isZero()) {
604 numcoeff.set(1, 1, 0);
605 } else {
606 numcoeff = xvar;
607 if(!numdeg.isOne()) {
608 numcoeff.raise(numdeg);
609 }
610 }
611 } else {
612 if(dencoeff.isNumber()) {
613 if(numcoeff.isNumber()) {
614 numcoeff.number() /= dencoeff.number();
615 } else {
616 numcoeff.calculateDivide(dencoeff, eo2);
617 }
618 } else {
619 MathStructure mcopy(numcoeff);
620 if((mcopy.equals(mnum, true, true) && dencoeff.equals(mden, true, true)) || !MathStructure::polynomialDivide(mcopy, dencoeff, numcoeff, eo2, false)) {
621 return false;
622 }
623 }
624 if(!numdeg.isZero() && !numcoeff.isZero()) {
625 if(numcoeff.isOne()) {
626 numcoeff = xvar;
627 if(!numdeg.isOne()) {
628 numcoeff.raise(numdeg);
629 }
630 } else {
631 numcoeff.multiply(xvar, true);
632 if(!numdeg.isOne()) {
633 numcoeff[numcoeff.size() - 1].raise(numdeg);
634 }
635 numcoeff.calculateMultiplyLast(eo2);
636 }
637 }
638 }
639 if(mquotient.isZero()) mquotient = numcoeff;
640 else mquotient.add(numcoeff, true);
641 if(numcoeff.isAddition() && mden.isAddition() && numcoeff.size() * mden.size() >= (eo.expand == -1 ? 50 : 500)) return false;
642 numcoeff.calculateMultiply(mden, eo2);
643 mrem.calculateSubtract(numcoeff, eo2);
644 if(mrem.isZero()) {
645 return true;
646 }
647 numdeg = mrem.degree(xvar);
648 }
649 return false;
650 }
651
divide_in_z(const MathStructure & mnum,const MathStructure & mden,MathStructure & mquotient,const sym_desc_vec & sym_stats,size_t var_i,const EvaluationOptions & eo)652 bool divide_in_z(const MathStructure &mnum, const MathStructure &mden, MathStructure &mquotient, const sym_desc_vec &sym_stats, size_t var_i, const EvaluationOptions &eo) {
653
654 mquotient.clear();
655 if(mden.isZero()) return false;
656 if(mnum.isZero()) return true;
657 if(mden.isOne()) {
658 mquotient = mnum;
659 return true;
660 }
661 if(mnum.isNumber()) {
662 if(!mden.isNumber()) {
663 return false;
664 }
665 mquotient = mnum;
666 return mquotient.number().divide(mden.number()) && mquotient.isInteger();
667 }
668 if(mnum == mden) {
669 mquotient.set(1, 1, 0);
670 return true;
671 }
672
673 if(mden.isPower()) {
674 MathStructure qbar(mnum);
675 for(Number ni(mden[1].number()); ni.isPositive(); ni--) {
676 if(!divide_in_z(qbar, mden[0], mquotient, sym_stats, var_i, eo)) return false;
677 qbar = mquotient;
678 }
679 return true;
680 }
681
682 if(mden.isMultiplication()) {
683 MathStructure qbar(mnum);
684 for(size_t i = 0; i < mden.size(); i++) {
685 sym_desc_vec sym_stats2;
686 get_symbol_stats(mnum, mden[i], sym_stats2);
687 if(!divide_in_z(qbar, mden[i], mquotient, sym_stats2, 0, eo)) return false;
688 qbar = mquotient;
689 }
690 return true;
691 }
692
693 if(var_i >= sym_stats.size()) return false;
694 const MathStructure &xvar = sym_stats[var_i].sym;
695
696 Number numdeg = mnum.degree(xvar);
697 Number dendeg = mden.degree(xvar);
698 if(dendeg.isGreaterThan(numdeg)) return false;
699 MathStructure dencoeff;
700 MathStructure mrem(mnum);
701 mden.coefficient(xvar, dendeg, dencoeff);
702
703 while(numdeg.isGreaterThanOrEqualTo(dendeg)) {
704 if(CALCULATOR->aborted()) return false;
705 MathStructure numcoeff;
706 mrem.coefficient(xvar, numdeg, numcoeff);
707 MathStructure term;
708 if(!divide_in_z(numcoeff, dencoeff, term, sym_stats, var_i + 1, eo)) break;
709 numdeg -= dendeg;
710 if(!numdeg.isZero() && !term.isZero()) {
711 if(term.isOne()) {
712 term = xvar;
713 if(!numdeg.isOne()) {
714 term.raise(numdeg);
715 }
716 } else {
717 term.multiply(xvar, true);
718 if(!numdeg.isOne()) {
719 term[term.size() - 1].raise(numdeg);
720 }
721 term.calculateMultiplyLast(eo);
722 }
723 }
724 if(mquotient.isZero()) {
725 mquotient = term;
726 } else {
727 mquotient.calculateAdd(term, eo);
728 }
729 if(term.isAddition() && mden.isAddition() && term.size() * mden.size() >= (eo.expand == -1 ? 50 : 500)) return false;
730 term.calculateMultiply(mden, eo);
731 mrem.calculateSubtract(term, eo);
732 if(mrem.isZero()) {
733 return true;
734 }
735 numdeg = mrem.degree(xvar);
736
737 }
738 return false;
739 }
740
prem(const MathStructure & mnum,const MathStructure & mden,const MathStructure & xvar,MathStructure & mrem,const EvaluationOptions & eo,bool check_args)741 bool prem(const MathStructure &mnum, const MathStructure &mden, const MathStructure &xvar, MathStructure &mrem, const EvaluationOptions &eo, bool check_args) {
742
743 mrem.clear();
744 if(mden.isZero()) {
745 //division by zero
746 return false;
747 }
748 if(mnum.isNumber()) {
749 if(!mden.isNumber()) {
750 mrem = mden;
751 }
752 return true;
753 }
754 if(check_args && (!mnum.isRationalPolynomial() || !mden.isRationalPolynomial())) {
755 return false;
756 }
757
758 mrem = mnum;
759 MathStructure eb(mden);
760 Number rdeg = mrem.degree(xvar);
761 Number bdeg = eb.degree(xvar);
762 MathStructure blcoeff;
763 if(bdeg.isLessThanOrEqualTo(rdeg)) {
764 eb.coefficient(xvar, bdeg, blcoeff);
765 if(bdeg == 0) {
766 eb.clear();
767 } else {
768 MathStructure mpow(xvar);
769 mpow.raise(bdeg);
770 mpow.calculateRaiseExponent(eo);
771 //POWER_CLEAN(mpow)
772 mpow.calculateMultiply(blcoeff, eo);
773 eb.calculateSubtract(mpow, eo);
774 }
775 } else {
776 blcoeff.set(1, 1, 0);
777 }
778
779 Number delta(rdeg);
780 delta -= bdeg;
781 delta++;
782 int i = 0;
783 while(rdeg.isGreaterThanOrEqualTo(bdeg) && !mrem.isZero()) {
784 if(CALCULATOR->aborted() || delta < i / 10) {mrem.clear(); return false;}
785 MathStructure rlcoeff;
786 mrem.coefficient(xvar, rdeg, rlcoeff);
787 MathStructure term(xvar);
788 term.raise(rdeg);
789 term[1].number() -= bdeg;
790 term.calculateRaiseExponent(eo);
791 //POWER_CLEAN(term)
792 term.calculateMultiply(rlcoeff, eo);
793 term.calculateMultiply(eb, eo);
794 if(rdeg == 0) {
795 mrem = term;
796 mrem.calculateNegate(eo);
797 } else {
798 if(!rdeg.isZero()) {
799 rlcoeff.multiply(xvar, true);
800 if(!rdeg.isOne()) {
801 rlcoeff[rlcoeff.size() - 1].raise(rdeg);
802 rlcoeff[rlcoeff.size() - 1].calculateRaiseExponent(eo);
803 }
804 rlcoeff.calculateMultiplyLast(eo);
805 }
806 mrem.calculateSubtract(rlcoeff, eo);
807 if(mrem.isAddition() && blcoeff.isAddition() && mrem.size() * blcoeff.size() >= (eo.expand == -1 ? 50 : 500)) {mrem.clear(); return false;}
808 mrem.calculateMultiply(blcoeff, eo);
809 mrem.calculateSubtract(term, eo);
810 }
811 rdeg = mrem.degree(xvar);
812 i++;
813 }
814 delta -= i;
815 blcoeff.raise(delta);
816 blcoeff.calculateRaiseExponent(eo);
817 mrem.calculateMultiply(blcoeff, eo);
818 return true;
819 }
820
maxCoefficient()821 Number MathStructure::maxCoefficient() {
822 if(isNumber()) {
823 Number nr(o_number);
824 nr.abs();
825 return nr;
826 } else if(isAddition()) {
827 Number cur_max(overallCoefficient());
828 cur_max.abs();
829 for(size_t i = 0; i < SIZE; i++) {
830 Number a(CHILD(i).overallCoefficient());
831 a.abs();
832 if(a.isGreaterThan(cur_max)) cur_max = a;
833 }
834 return cur_max;
835 } else if(isMultiplication()) {
836 Number nr(overallCoefficient());
837 nr.abs();
838 return nr;
839 } else {
840 return nr_one;
841 }
842 }
polynomial_smod(const MathStructure & mpoly,const Number & xi,MathStructure & msmod,const EvaluationOptions & eo,MathStructure * mparent,size_t index_smod)843 void polynomial_smod(const MathStructure &mpoly, const Number &xi, MathStructure &msmod, const EvaluationOptions &eo, MathStructure *mparent, size_t index_smod) {
844 if(mpoly.isNumber()) {
845 msmod = mpoly;
846 msmod.number().smod(xi);
847 } else if(mpoly.isAddition()) {
848 msmod.clear();
849 msmod.setType(STRUCT_ADDITION);
850 msmod.resizeVector(mpoly.size(), m_zero);
851 for(size_t i = 0; i < mpoly.size(); i++) {
852 polynomial_smod(mpoly[i], xi, msmod[i], eo, &msmod, i);
853 }
854 msmod.calculatesub(eo, eo, false, mparent, index_smod);
855 } else if(mpoly.isMultiplication()) {
856 msmod = mpoly;
857 if(msmod.size() > 0 && msmod[0].isNumber()) {
858 if(!msmod[0].number().smod(xi) || msmod[0].isZero()) {
859 msmod.clear();
860 }
861 }
862 } else {
863 msmod = mpoly;
864 }
865 }
contains_zerointerval_multiplier(MathStructure & mstruct)866 bool contains_zerointerval_multiplier(MathStructure &mstruct) {
867 if(mstruct.isAddition()) {
868 for(size_t i = 0; i < mstruct.size(); i++) {
869 if(contains_zerointerval_multiplier(mstruct[i])) return true;
870 }
871 } else if(mstruct.isMultiplication()) {
872 for(size_t i = 0; i < mstruct.size(); i++) {
873 if(mstruct[i].isNumber() && !mstruct[i].number().isNonZero()) return true;
874 }
875 } else if(mstruct.isNumber() && !mstruct.number().isNonZero()) {
876 return true;
877 }
878 return false;
879 }
880
polynomial_long_division(const MathStructure & mnum,const MathStructure & mden,const MathStructure & xvar_pre,MathStructure & mquotient,MathStructure & mrem,const EvaluationOptions & eo,bool check_args,bool for_newtonraphson)881 bool polynomial_long_division(const MathStructure &mnum, const MathStructure &mden, const MathStructure &xvar_pre, MathStructure &mquotient, MathStructure &mrem, const EvaluationOptions &eo, bool check_args, bool for_newtonraphson) {
882 mquotient.clear();
883 mrem.clear();
884
885 if(CALCULATOR->aborted()) return false;
886 if(mden.isZero()) {
887 //division by zero
888 mrem.set(mnum);
889 return false;
890 }
891 if(mnum.isZero()) {
892 mquotient.clear();
893 return true;
894 }
895 if(mden.isNumber() && mnum.isNumber()) {
896 mquotient = mnum;
897 if(IS_REAL(mden) && IS_REAL(mnum)) {
898 mquotient.number() /= mden.number();
899 } else {
900 mquotient.calculateDivide(mden, eo);
901 }
902 return true;
903 } else if(mnum.isNumber()) {
904 mrem.set(mnum);
905 return false;
906 }
907 if(mnum == mden) {
908 mquotient.set(1, 1, 0);
909 return true;
910 }
911
912 mrem.set(mnum);
913
914 if(check_args && (!mnum.isRationalPolynomial(true, true) || !mden.isRationalPolynomial(true, true))) {
915 return false;
916 }
917
918 MathStructure xvar(xvar_pre);
919 if(xvar.isZero()) {
920 vector<MathStructure> symsd, symsn;
921 collect_symbols(mnum, symsn);
922 if(symsn.empty()) return false;
923 collect_symbols(mden, symsd);
924 if(symsd.empty()) return false;
925 for(size_t i = 0; i < symsd.size();) {
926 bool b_found = false;
927 for(size_t i2 = 0; i2 < symsn.size(); i2++) {
928 if(symsn[i2] == symsd[i]) {
929 b_found = (mnum.degree(symsd[i]) >= mden.degree(symsd[i]));
930 break;
931 }
932 }
933 if(b_found) i++;
934 else symsd.erase(symsd.begin() + i);
935 }
936 for(size_t i = 0; i < symsd.size(); i++) {
937 MathStructure mquotient2, mrem2;
938 if(polynomial_long_division(mrem, mden, symsd[i], mquotient2, mrem2, eo, false, for_newtonraphson) && !mquotient2.isZero()) {
939 mrem = mrem2;
940 if(mquotient.isZero()) mquotient = mquotient2;
941 else mquotient.calculateAdd(mquotient2, eo);
942 }
943 if(CALCULATOR->aborted()) {
944 mrem = mnum;
945 mquotient.clear();
946 return false;
947 }
948 }
949 return true;
950 }
951
952 EvaluationOptions eo2 = eo;
953 eo2.keep_zero_units = false;
954 eo2.do_polynomial_division = false;
955
956 Number numdeg = mnum.degree(xvar);
957 Number dendeg = mden.degree(xvar);
958 MathStructure dencoeff;
959 mden.coefficient(xvar, dendeg, dencoeff);
960
961 for(size_t i = 0; numdeg.isGreaterThanOrEqualTo(dendeg); i++) {
962 if(i > 10000 || CALCULATOR->aborted()) {
963 mrem = mnum;
964 mquotient.clear();
965 return false;
966 }
967 MathStructure numcoeff;
968 mrem.coefficient(xvar, numdeg, numcoeff);
969 numdeg -= dendeg;
970 if(numcoeff == dencoeff) {
971 if(numdeg.isZero()) {
972 numcoeff.set(1, 1, 0);
973 } else {
974 numcoeff = xvar;
975 if(!numdeg.isOne()) {
976 numcoeff.raise(numdeg);
977 }
978 }
979 } else {
980 if(dencoeff.isNumber()) {
981 if(numcoeff.isNumber()) {
982 numcoeff.number() /= dencoeff.number();
983 } else {
984 numcoeff.calculateDivide(dencoeff, eo2);
985 }
986 } else {
987 if(for_newtonraphson) {mrem = mnum; mquotient.clear(); return false;}
988 MathStructure mcopy(numcoeff);
989 if(!MathStructure::polynomialDivide(mcopy, dencoeff, numcoeff, eo2, check_args)) {
990 if(CALCULATOR->aborted()) {mrem = mnum; mquotient.clear(); return false;}
991 break;
992 }
993 }
994 if(!numdeg.isZero() && !numcoeff.isZero()) {
995 if(numcoeff.isOne()) {
996 numcoeff = xvar;
997 if(!numdeg.isOne()) {
998 numcoeff.raise(numdeg);
999 }
1000 } else {
1001 numcoeff.multiply(xvar, true);
1002 if(!numdeg.isOne()) {
1003 numcoeff[numcoeff.size() - 1].raise(numdeg);
1004 }
1005 numcoeff.calculateMultiplyLast(eo2);
1006 }
1007 }
1008 }
1009 if(mquotient.isZero()) mquotient = numcoeff;
1010 else mquotient.calculateAdd(numcoeff, eo2);
1011 if(numcoeff.isAddition() && mden.isAddition() && numcoeff.size() * mden.size() >= (eo.expand == -1 ? 50 : 500)) {mrem = mnum; mquotient.clear(); return false;}
1012 numcoeff.calculateMultiply(mden, eo2);
1013 mrem.calculateSubtract(numcoeff, eo2);
1014 if(contains_zerointerval_multiplier(mquotient)) {mrem = mnum; mquotient.clear(); return false;}
1015 if(mrem.isZero() || (for_newtonraphson && mrem.isNumber())) break;
1016 if(contains_zerointerval_multiplier(mrem)) {mrem = mnum; mquotient.clear(); return false;}
1017 numdeg = mrem.degree(xvar);
1018 }
1019 return true;
1020
1021 }
isRationalPolynomial(bool allow_non_rational_coefficient,bool allow_interval_coefficient) const1022 bool MathStructure::isRationalPolynomial(bool allow_non_rational_coefficient, bool allow_interval_coefficient) const {
1023 switch(m_type) {
1024 case STRUCT_NUMBER: {
1025 if(allow_interval_coefficient) return o_number.isReal() && o_number.isNonZero();
1026 if(allow_non_rational_coefficient) return o_number.isReal() && !o_number.isInterval() && o_number.isNonZero();
1027 return o_number.isRational() && !o_number.isZero();
1028 }
1029 case STRUCT_MULTIPLICATION: {
1030 for(size_t i = 0; i < SIZE; i++) {
1031 if(CHILD(i).isAddition() || CHILD(i).isMultiplication() || !CHILD(i).isRationalPolynomial(allow_non_rational_coefficient, allow_interval_coefficient)) {
1032 return false;
1033 }
1034 }
1035 return true;
1036 }
1037 case STRUCT_ADDITION: {
1038 for(size_t i = 0; i < SIZE; i++) {
1039 if(CHILD(i).isAddition() || !CHILD(i).isRationalPolynomial(allow_non_rational_coefficient, allow_interval_coefficient)) {
1040 return false;
1041 }
1042 }
1043 return true;
1044 }
1045 case STRUCT_POWER: {
1046 return CHILD(1).isInteger() && CHILD(1).number().isNonNegative() && !CHILD(1).number().isOne() && CHILD(1).number() < 1000 && !CHILD(0).isNumber() && !CHILD(0).isMultiplication() && !CHILD(0).isAddition() && !CHILD(0).isPower() && CHILD(0).isRationalPolynomial(allow_non_rational_coefficient, allow_interval_coefficient);
1047 }
1048 case STRUCT_FUNCTION: {
1049 if(o_function->id() == FUNCTION_ID_UNCERTAINTY || o_function->id() == FUNCTION_ID_INTERVAL || containsInterval() || containsInfinity()) return false;
1050 }
1051 case STRUCT_UNIT: {}
1052 case STRUCT_VARIABLE: {}
1053 case STRUCT_SYMBOLIC: {
1054 return representsNonMatrix() && !representsUndefined(true, true);
1055 }
1056 default: {}
1057 }
1058 return false;
1059 }
overallCoefficient() const1060 const Number &MathStructure::overallCoefficient() const {
1061 switch(m_type) {
1062 case STRUCT_NUMBER: {
1063 return o_number;
1064 }
1065 case STRUCT_MULTIPLICATION: {
1066 for(size_t i = 0; i < SIZE; i++) {
1067 if(CHILD(i).isNumber()) {
1068 return CHILD(i).number();
1069 }
1070 }
1071 return nr_one;
1072 }
1073 case STRUCT_ADDITION: {
1074 for(size_t i = 0; i < SIZE; i++) {
1075 if(CHILD(i).isNumber()) {
1076 return CHILD(i).number();
1077 }
1078 }
1079 return nr_zero;
1080 }
1081 case STRUCT_POWER: {
1082 return nr_zero;
1083 }
1084 default: {}
1085 }
1086 return nr_zero;
1087 }
1088
1089