GiNaC 1.8.6
normal.cpp
Go to the documentation of this file.
1
8/*
9 * GiNaC Copyright (C) 1999-2023 Johannes Gutenberg University Mainz, Germany
10 *
11 * This program is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation; either version 2 of the License, or
14 * (at your option) any later version.
15 *
16 * This program is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with this program; if not, write to the Free Software
23 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
25
26#include "normal.h"
27#include "basic.h"
28#include "ex.h"
29#include "add.h"
30#include "constant.h"
31#include "expairseq.h"
32#include "fail.h"
33#include "inifcns.h"
34#include "lst.h"
35#include "mul.h"
36#include "numeric.h"
37#include "power.h"
38#include "relational.h"
39#include "operators.h"
40#include "matrix.h"
41#include "pseries.h"
42#include "symbol.h"
43#include "utils.h"
44#include "polynomial/chinrem_gcd.h"
45
46#include <algorithm>
47#include <map>
48
49namespace GiNaC {
50
51// If comparing expressions (ex::compare()) is fast, you can set this to 1.
52// Some routines like quo(), rem() and gcd() will then return a quick answer
53// when they are called with two identical arguments.
54#define FAST_COMPARE 1
55
56// Set this if you want divide_in_z() to use remembering
57#define USE_REMEMBER 0
58
59// Set this if you want divide_in_z() to use trial division followed by
60// polynomial interpolation (always slower except for completely dense
61// polynomials)
62#define USE_TRIAL_DIVISION 0
63
64// Set this to enable some statistical output for the GCD routines
65#define STATISTICS 0
66
67
68#if STATISTICS
69// Statistics variables
70static int gcd_called = 0;
71static int sr_gcd_called = 0;
72static int heur_gcd_called = 0;
73static int heur_gcd_failed = 0;
74
75// Print statistics at end of program
76static struct _stat_print {
77 _stat_print() {}
78 ~_stat_print() {
79 std::cout << "gcd() called " << gcd_called << " times\n";
80 std::cout << "sr_gcd() called " << sr_gcd_called << " times\n";
81 std::cout << "heur_gcd() called " << heur_gcd_called << " times\n";
82 std::cout << "heur_gcd() failed " << heur_gcd_failed << " times\n";
83 }
84} stat_print;
85#endif
86
87
95static bool get_first_symbol(const ex &e, ex &x)
96{
97 if (is_a<symbol>(e)) {
98 x = e;
99 return true;
100 } else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
101 for (size_t i=0; i<e.nops(); i++)
102 if (get_first_symbol(e.op(i), x))
103 return true;
104 } else if (is_exactly_a<power>(e)) {
105 if (get_first_symbol(e.op(0), x))
106 return true;
107 }
108 return false;
109}
110
111
112/*
113 * Statistical information about symbols in polynomials
114 */
115
122struct sym_desc {
124 sym_desc(const ex& s)
125 : sym(s), deg_a(0), deg_b(0), ldeg_a(0), ldeg_b(0), max_deg(0), max_lcnops(0)
126 { }
127
130
132 int deg_a;
133
135 int deg_b;
136
139
142
145
148
150 bool operator<(const sym_desc &x) const
151 {
152 if (max_deg == x.max_deg)
153 return max_lcnops < x.max_lcnops;
154 else
155 return max_deg < x.max_deg;
156 }
157};
158
159// Vector of sym_desc structures
160typedef std::vector<sym_desc> sym_desc_vec;
161
162// Add symbol the sym_desc_vec (used internally by get_symbol_stats())
163static void add_symbol(const ex &s, sym_desc_vec &v)
164{
165 for (auto & it : v)
166 if (it.sym.is_equal(s)) // If it's already in there, don't add it a second time
167 return;
168
169 v.push_back(sym_desc(s));
170}
171
172// Collect all symbols of an expression (used internally by get_symbol_stats())
173static void collect_symbols(const ex &e, sym_desc_vec &v)
174{
175 if (is_a<symbol>(e)) {
176 add_symbol(e, v);
177 } else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
178 for (size_t i=0; i<e.nops(); i++)
179 collect_symbols(e.op(i), v);
180 } else if (is_exactly_a<power>(e)) {
181 collect_symbols(e.op(0), v);
182 }
183}
184
197static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
198{
199 collect_symbols(a, v);
200 collect_symbols(b, v);
201 for (auto & it : v) {
202 int deg_a = a.degree(it.sym);
203 int deg_b = b.degree(it.sym);
204 it.deg_a = deg_a;
205 it.deg_b = deg_b;
206 it.max_deg = std::max(deg_a, deg_b);
207 it.max_lcnops = std::max(a.lcoeff(it.sym).nops(), b.lcoeff(it.sym).nops());
208 it.ldeg_a = a.ldegree(it.sym);
209 it.ldeg_b = b.ldegree(it.sym);
210 }
211 std::sort(v.begin(), v.end());
212
213#if 0
214 std::clog << "Symbols:\n";
215 auto it = v.begin(), itend = v.end();
216 while (it != itend) {
217 std::clog << " " << it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << std::endl;
218 std::clog << " lcoeff_a=" << a.lcoeff(it->sym) << ", lcoeff_b=" << b.lcoeff(it->sym) << std::endl;
219 ++it;
220 }
221#endif
222}
223
224
225/*
226 * Computation of LCM of denominators of coefficients of a polynomial
227 */
228
229// Compute LCM of denominators of coefficients by going through the
230// expression recursively (used internally by lcm_of_coefficients_denominators())
231static numeric lcmcoeff(const ex &e, const numeric &l)
232{
234 return lcm(ex_to<numeric>(e).denom(), l);
235 else if (is_exactly_a<add>(e)) {
236 numeric c = *_num1_p;
237 for (size_t i=0; i<e.nops(); i++)
238 c = lcmcoeff(e.op(i), c);
239 return lcm(c, l);
240 } else if (is_exactly_a<mul>(e)) {
241 numeric c = *_num1_p;
242 for (size_t i=0; i<e.nops(); i++)
243 c *= lcmcoeff(e.op(i), *_num1_p);
244 return lcm(c, l);
245 } else if (is_exactly_a<power>(e)) {
246 if (is_a<symbol>(e.op(0)))
247 return l;
248 else
249 return pow(lcmcoeff(e.op(0), l), ex_to<numeric>(e.op(1)));
250 }
251 return l;
252}
253
262{
263 return lcmcoeff(e, *_num1_p);
264}
265
271static ex multiply_lcm(const ex &e, const numeric &lcm)
272{
273 if (lcm.is_equal(*_num1_p))
274 // e * 1 -> e;
275 return e;
276
277 if (is_exactly_a<mul>(e)) {
278 // (a*b*...)*lcm -> (a*lcma)*(b*lcmb)*...*(lcm/(lcma*lcmb*...))
279 size_t num = e.nops();
280 exvector v;
281 v.reserve(num + 1);
282 numeric lcm_accum = *_num1_p;
283 for (size_t i=0; i<num; i++) {
284 numeric op_lcm = lcmcoeff(e.op(i), *_num1_p);
285 v.push_back(multiply_lcm(e.op(i), op_lcm));
286 lcm_accum *= op_lcm;
287 }
288 v.push_back(lcm / lcm_accum);
289 return dynallocate<mul>(v);
290 } else if (is_exactly_a<add>(e)) {
291 // (a+b+...)*lcm -> a*lcm+b*lcm+...
292 size_t num = e.nops();
293 exvector v;
294 v.reserve(num);
295 for (size_t i=0; i<num; i++)
296 v.push_back(multiply_lcm(e.op(i), lcm));
297 return dynallocate<add>(v);
298 } else if (is_exactly_a<power>(e)) {
299 if (!is_a<symbol>(e.op(0))) {
300 // (b^e)*lcm -> (b*lcm^(1/e))^e if lcm^(1/e) ∈ ℚ (i.e. not a float)
301 // but not for symbolic b, as evaluation would undo this again
302 numeric root_of_lcm = lcm.power(ex_to<numeric>(e.op(1)).inverse());
303 if (root_of_lcm.is_rational())
304 return pow(multiply_lcm(e.op(0), root_of_lcm), e.op(1));
305 }
306 }
307 // can't recurse down into e
308 return dynallocate<mul>(e, lcm);
309}
310
311
319{
320 return bp->integer_content();
321}
322
324{
325 return *_num1_p;
326}
327
329{
330 return abs(*this);
331}
332
334{
335 numeric c = *_num0_p, l = *_num1_p;
336 for (auto & it : seq) {
337 GINAC_ASSERT(!is_exactly_a<numeric>(it.rest));
338 GINAC_ASSERT(is_exactly_a<numeric>(it.coeff));
339 c = gcd(ex_to<numeric>(it.coeff).numer(), c);
340 l = lcm(ex_to<numeric>(it.coeff).denom(), l);
341 }
342 GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
343 c = gcd(ex_to<numeric>(overall_coeff).numer(), c);
344 l = lcm(ex_to<numeric>(overall_coeff).denom(), l);
345 return c/l;
346}
347
349{
350#ifdef DO_GINAC_ASSERT
351 for (auto & it : seq) {
352 GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(it)));
353 }
354#endif // def DO_GINAC_ASSERT
355 GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
356 return abs(ex_to<numeric>(overall_coeff));
357}
358
359
360/*
361 * Polynomial quotients and remainders
362 */
363
373ex quo(const ex &a, const ex &b, const ex &x, bool check_args)
374{
375 if (b.is_zero())
376 throw(std::overflow_error("quo: division by zero"));
377 if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b))
378 return a / b;
379#if FAST_COMPARE
380 if (a.is_equal(b))
381 return _ex1;
382#endif
384 throw(std::invalid_argument("quo: arguments must be polynomials over the rationals"));
385
386 // Polynomial long division
387 ex r = a.expand();
388 if (r.is_zero())
389 return r;
390 int bdeg = b.degree(x);
391 int rdeg = r.degree(x);
392 ex blcoeff = b.expand().coeff(x, bdeg);
393 bool blcoeff_is_numeric = is_exactly_a<numeric>(blcoeff);
394 exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
395 while (rdeg >= bdeg) {
396 ex term, rcoeff = r.coeff(x, rdeg);
397 if (blcoeff_is_numeric)
398 term = rcoeff / blcoeff;
399 else {
400 if (!divide(rcoeff, blcoeff, term, false))
401 return dynallocate<fail>();
402 }
403 term *= pow(x, rdeg - bdeg);
404 v.push_back(term);
405 r -= (term * b).expand();
406 if (r.is_zero())
407 break;
408 rdeg = r.degree(x);
409 }
410 return dynallocate<add>(v);
411}
412
413
423ex rem(const ex &a, const ex &b, const ex &x, bool check_args)
424{
425 if (b.is_zero())
426 throw(std::overflow_error("rem: division by zero"));
427 if (is_exactly_a<numeric>(a)) {
428 if (is_exactly_a<numeric>(b))
429 return _ex0;
430 else
431 return a;
432 }
433#if FAST_COMPARE
434 if (a.is_equal(b))
435 return _ex0;
436#endif
438 throw(std::invalid_argument("rem: arguments must be polynomials over the rationals"));
439
440 // Polynomial long division
441 ex r = a.expand();
442 if (r.is_zero())
443 return r;
444 int bdeg = b.degree(x);
445 int rdeg = r.degree(x);
446 ex blcoeff = b.expand().coeff(x, bdeg);
447 bool blcoeff_is_numeric = is_exactly_a<numeric>(blcoeff);
448 while (rdeg >= bdeg) {
449 ex term, rcoeff = r.coeff(x, rdeg);
450 if (blcoeff_is_numeric)
451 term = rcoeff / blcoeff;
452 else {
453 if (!divide(rcoeff, blcoeff, term, false))
454 return dynallocate<fail>();
455 }
456 term *= pow(x, rdeg - bdeg);
457 r -= (term * b).expand();
458 if (r.is_zero())
459 break;
460 rdeg = r.degree(x);
461 }
462 return r;
463}
464
465
472ex decomp_rational(const ex &a, const ex &x)
473{
474 ex nd = numer_denom(a);
475 ex numer = nd.op(0), denom = nd.op(1);
476 ex q = quo(numer, denom, x);
477 if (is_exactly_a<fail>(q))
478 return a;
479 else
480 return q + rem(numer, denom, x) / denom;
481}
482
483
492ex prem(const ex &a, const ex &b, const ex &x, bool check_args)
493{
494 if (b.is_zero())
495 throw(std::overflow_error("prem: division by zero"));
496 if (is_exactly_a<numeric>(a)) {
497 if (is_exactly_a<numeric>(b))
498 return _ex0;
499 else
500 return b;
501 }
503 throw(std::invalid_argument("prem: arguments must be polynomials over the rationals"));
504
505 // Polynomial long division
506 ex r = a.expand();
507 ex eb = b.expand();
508 int rdeg = r.degree(x);
509 int bdeg = eb.degree(x);
510 ex blcoeff;
511 if (bdeg <= rdeg) {
512 blcoeff = eb.coeff(x, bdeg);
513 if (bdeg == 0)
514 eb = _ex0;
515 else
516 eb -= blcoeff * pow(x, bdeg);
517 } else
518 blcoeff = _ex1;
519
520 int delta = rdeg - bdeg + 1, i = 0;
521 while (rdeg >= bdeg && !r.is_zero()) {
522 ex rlcoeff = r.coeff(x, rdeg);
523 ex term = (pow(x, rdeg - bdeg) * eb * rlcoeff).expand();
524 if (rdeg == 0)
525 r = _ex0;
526 else
527 r -= rlcoeff * pow(x, rdeg);
528 r = (blcoeff * r).expand() - term;
529 rdeg = r.degree(x);
530 i++;
531 }
532 return pow(blcoeff, delta - i) * r;
533}
534
535
544ex sprem(const ex &a, const ex &b, const ex &x, bool check_args)
545{
546 if (b.is_zero())
547 throw(std::overflow_error("prem: division by zero"));
548 if (is_exactly_a<numeric>(a)) {
549 if (is_exactly_a<numeric>(b))
550 return _ex0;
551 else
552 return b;
553 }
555 throw(std::invalid_argument("prem: arguments must be polynomials over the rationals"));
556
557 // Polynomial long division
558 ex r = a.expand();
559 ex eb = b.expand();
560 int rdeg = r.degree(x);
561 int bdeg = eb.degree(x);
562 ex blcoeff;
563 if (bdeg <= rdeg) {
564 blcoeff = eb.coeff(x, bdeg);
565 if (bdeg == 0)
566 eb = _ex0;
567 else
568 eb -= blcoeff * pow(x, bdeg);
569 } else
570 blcoeff = _ex1;
571
572 while (rdeg >= bdeg && !r.is_zero()) {
573 ex rlcoeff = r.coeff(x, rdeg);
574 ex term = (pow(x, rdeg - bdeg) * eb * rlcoeff).expand();
575 if (rdeg == 0)
576 r = _ex0;
577 else
578 r -= rlcoeff * pow(x, rdeg);
579 r = (blcoeff * r).expand() - term;
580 rdeg = r.degree(x);
581 }
582 return r;
583}
584
585
595bool divide(const ex &a, const ex &b, ex &q, bool check_args)
596{
597 if (b.is_zero())
598 throw(std::overflow_error("divide: division by zero"));
599 if (a.is_zero()) {
600 q = _ex0;
601 return true;
602 }
603 if (is_exactly_a<numeric>(b)) {
604 q = a / b;
605 return true;
606 } else if (is_exactly_a<numeric>(a))
607 return false;
608#if FAST_COMPARE
609 if (a.is_equal(b)) {
610 q = _ex1;
611 return true;
612 }
613#endif
614 if (check_args && (!a.info(info_flags::rational_polynomial) ||
616 throw(std::invalid_argument("divide: arguments must be polynomials over the rationals"));
617
618 // Find first symbol
619 ex x;
620 if (!get_first_symbol(a, x) && !get_first_symbol(b, x))
621 throw(std::invalid_argument("invalid expression in divide()"));
622
623 // Try to avoid expanding partially factored expressions.
624 if (is_exactly_a<mul>(b)) {
625 // Divide sequentially by each term
626 ex rem_new, rem_old = a;
627 for (size_t i=0; i < b.nops(); i++) {
628 if (! divide(rem_old, b.op(i), rem_new, false))
629 return false;
630 rem_old = rem_new;
631 }
632 q = rem_new;
633 return true;
634 } else if (is_exactly_a<power>(b)) {
635 const ex& bb(b.op(0));
636 int exp_b = ex_to<numeric>(b.op(1)).to_int();
637 ex rem_new, rem_old = a;
638 for (int i=exp_b; i>0; i--) {
639 if (! divide(rem_old, bb, rem_new, false))
640 return false;
641 rem_old = rem_new;
642 }
643 q = rem_new;
644 return true;
645 }
646
647 if (is_exactly_a<mul>(a)) {
648 // Divide sequentially each term. If some term in a is divisible
649 // by b we are done... and if not, we can't really say anything.
650 size_t i;
651 ex rem_i;
652 bool divisible_p = false;
653 for (i=0; i < a.nops(); ++i) {
654 if (divide(a.op(i), b, rem_i, false)) {
655 divisible_p = true;
656 break;
657 }
658 }
659 if (divisible_p) {
660 exvector resv;
661 resv.reserve(a.nops());
662 for (size_t j=0; j < a.nops(); j++) {
663 if (j==i)
664 resv.push_back(rem_i);
665 else
666 resv.push_back(a.op(j));
667 }
668 q = dynallocate<mul>(resv);
669 return true;
670 }
671 } else if (is_exactly_a<power>(a)) {
672 // The base itself might be divisible by b, in that case we don't
673 // need to expand a
674 const ex& ab(a.op(0));
675 int a_exp = ex_to<numeric>(a.op(1)).to_int();
676 ex rem_i;
677 if (divide(ab, b, rem_i, false)) {
678 q = rem_i * pow(ab, a_exp - 1);
679 return true;
680 }
681// code below is commented-out because it leads to a significant slowdown
682// for (int i=2; i < a_exp; i++) {
683// if (divide(power(ab, i), b, rem_i, false)) {
684// q = rem_i*power(ab, a_exp - i);
685// return true;
686// }
687// } // ... so we *really* need to expand expression.
688 }
689
690 // Polynomial long division (recursive)
691 ex r = a.expand();
692 if (r.is_zero()) {
693 q = _ex0;
694 return true;
695 }
696 int bdeg = b.degree(x);
697 int rdeg = r.degree(x);
698 ex blcoeff = b.expand().coeff(x, bdeg);
699 bool blcoeff_is_numeric = is_exactly_a<numeric>(blcoeff);
700 exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
701 while (rdeg >= bdeg) {
702 ex term, rcoeff = r.coeff(x, rdeg);
703 if (blcoeff_is_numeric)
704 term = rcoeff / blcoeff;
705 else
706 if (!divide(rcoeff, blcoeff, term, false))
707 return false;
708 term *= pow(x, rdeg - bdeg);
709 v.push_back(term);
710 r -= (term * b).expand();
711 if (r.is_zero()) {
712 q = dynallocate<add>(v);
713 return true;
714 }
715 rdeg = r.degree(x);
716 }
717 return false;
718}
719
720
721#if USE_REMEMBER
722/*
723 * Remembering
724 */
725
726typedef std::pair<ex, ex> ex2;
727typedef std::pair<ex, bool> exbool;
728
729struct ex2_less {
730 bool operator() (const ex2 &p, const ex2 &q) const
731 {
732 int cmp = p.first.compare(q.first);
733 return ((cmp<0) || (!(cmp>0) && p.second.compare(q.second)<0));
734 }
735};
736
737typedef std::map<ex2, exbool, ex2_less> ex2_exbool_remember;
738#endif
739
740
757static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_iterator var)
758{
759 q = _ex0;
760 if (b.is_zero())
761 throw(std::overflow_error("divide_in_z: division by zero"));
762 if (b.is_equal(_ex1)) {
763 q = a;
764 return true;
765 }
766 if (is_exactly_a<numeric>(a)) {
767 if (is_exactly_a<numeric>(b)) {
768 q = a / b;
769 return q.info(info_flags::integer);
770 } else
771 return false;
772 }
773#if FAST_COMPARE
774 if (a.is_equal(b)) {
775 q = _ex1;
776 return true;
777 }
778#endif
779
780#if USE_REMEMBER
781 // Remembering
782 static ex2_exbool_remember dr_remember;
783 ex2_exbool_remember::const_iterator remembered = dr_remember.find(ex2(a, b));
784 if (remembered != dr_remember.end()) {
785 q = remembered->second.first;
786 return remembered->second.second;
787 }
788#endif
789
790 if (is_exactly_a<power>(b)) {
791 const ex& bb(b.op(0));
792 ex qbar = a;
793 int exp_b = ex_to<numeric>(b.op(1)).to_int();
794 for (int i=exp_b; i>0; i--) {
795 if (!divide_in_z(qbar, bb, q, var))
796 return false;
797 qbar = q;
798 }
799 return true;
800 }
801
802 if (is_exactly_a<mul>(b)) {
803 ex qbar = a;
804 for (const auto & it : b) {
805 sym_desc_vec sym_stats;
806 get_symbol_stats(a, it, sym_stats);
807 if (!divide_in_z(qbar, it, q, sym_stats.begin()))
808 return false;
809
810 qbar = q;
811 }
812 return true;
813 }
814
815 // Main symbol
816 const ex &x = var->sym;
817
818 // Compare degrees
819 int adeg = a.degree(x), bdeg = b.degree(x);
820 if (bdeg > adeg)
821 return false;
822
823#if USE_TRIAL_DIVISION
824
825 // Trial division with polynomial interpolation
826 int i, k;
827
828 // Compute values at evaluation points 0..adeg
829 vector<numeric> alpha; alpha.reserve(adeg + 1);
830 exvector u; u.reserve(adeg + 1);
831 numeric point = *_num0_p;
832 ex c;
833 for (i=0; i<=adeg; i++) {
834 ex bs = b.subs(x == point, subs_options::no_pattern);
835 while (bs.is_zero()) {
836 point += *_num1_p;
837 bs = b.subs(x == point, subs_options::no_pattern);
838 }
839 if (!divide_in_z(a.subs(x == point, subs_options::no_pattern), bs, c, var+1))
840 return false;
841 alpha.push_back(point);
842 u.push_back(c);
843 point += *_num1_p;
844 }
845
846 // Compute inverses
847 vector<numeric> rcp; rcp.reserve(adeg + 1);
848 rcp.push_back(*_num0_p);
849 for (k=1; k<=adeg; k++) {
850 numeric product = alpha[k] - alpha[0];
851 for (i=1; i<k; i++)
852 product *= alpha[k] - alpha[i];
853 rcp.push_back(product.inverse());
854 }
855
856 // Compute Newton coefficients
857 exvector v; v.reserve(adeg + 1);
858 v.push_back(u[0]);
859 for (k=1; k<=adeg; k++) {
860 ex temp = v[k - 1];
861 for (i=k-2; i>=0; i--)
862 temp = temp * (alpha[k] - alpha[i]) + v[i];
863 v.push_back((u[k] - temp) * rcp[k]);
864 }
865
866 // Convert from Newton form to standard form
867 c = v[adeg];
868 for (k=adeg-1; k>=0; k--)
869 c = c * (x - alpha[k]) + v[k];
870
871 if (c.degree(x) == (adeg - bdeg)) {
872 q = c.expand();
873 return true;
874 } else
875 return false;
876
877#else
878
879 // Polynomial long division (recursive)
880 ex r = a.expand();
881 if (r.is_zero())
882 return true;
883 int rdeg = adeg;
884 ex eb = b.expand();
885 ex blcoeff = eb.coeff(x, bdeg);
886 exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
887 while (rdeg >= bdeg) {
888 ex term, rcoeff = r.coeff(x, rdeg);
889 if (!divide_in_z(rcoeff, blcoeff, term, var+1))
890 break;
891 term = (term * pow(x, rdeg - bdeg)).expand();
892 v.push_back(term);
893 r -= (term * eb).expand();
894 if (r.is_zero()) {
895 q = dynallocate<add>(v);
896#if USE_REMEMBER
897 dr_remember[ex2(a, b)] = exbool(q, true);
898#endif
899 return true;
900 }
901 rdeg = r.degree(x);
902 }
903#if USE_REMEMBER
904 dr_remember[ex2(a, b)] = exbool(q, false);
905#endif
906 return false;
907
908#endif
909}
910
911
912/*
913 * Separation of unit part, content part and primitive part of polynomials
914 */
915
923ex ex::unit(const ex &x) const
924{
925 ex c = expand().lcoeff(x);
926 if (is_exactly_a<numeric>(c))
927 return c.info(info_flags::negative) ?_ex_1 : _ex1;
928 else {
929 ex y;
930 if (get_first_symbol(c, y))
931 return c.unit(y);
932 else
933 throw(std::invalid_argument("invalid expression in unit()"));
934 }
935}
936
937
945ex ex::content(const ex &x) const
946{
947 if (is_exactly_a<numeric>(*this))
948 return info(info_flags::negative) ? -*this : *this;
949
950 ex e = expand();
951 if (e.is_zero())
952 return _ex0;
953
954 // First, divide out the integer content (which we can calculate very efficiently).
955 // If the leading coefficient of the quotient is an integer, we are done.
956 ex c = e.integer_content();
957 ex r = e / c;
958 int deg = r.degree(x);
959 ex lcoeff = r.coeff(x, deg);
961 return c;
962
963 // GCD of all coefficients
964 int ldeg = r.ldegree(x);
965 if (deg == ldeg)
966 return lcoeff * c / lcoeff.unit(x);
967 ex cont = _ex0;
968 for (int i=ldeg; i<=deg; i++)
969 cont = gcd(r.coeff(x, i), cont, nullptr, nullptr, false);
970 return cont * c;
971}
972
973
981ex ex::primpart(const ex &x) const
982{
983 // We need to compute the unit and content anyway, so call unitcontprim()
984 ex u, c, p;
985 unitcontprim(x, u, c, p);
986 return p;
987}
988
989
997ex ex::primpart(const ex &x, const ex &c) const
998{
999 if (is_zero() || c.is_zero())
1000 return _ex0;
1001 if (is_exactly_a<numeric>(*this))
1002 return _ex1;
1003
1004 // Divide by unit and content to get primitive part
1005 ex u = unit(x);
1006 if (is_exactly_a<numeric>(c))
1007 return *this / (c * u);
1008 else
1009 return quo(*this, c * u, x, false);
1010}
1011
1012
1022void ex::unitcontprim(const ex &x, ex &u, ex &c, ex &p) const
1023{
1024 // Quick check for zero (avoid expanding)
1025 if (is_zero()) {
1026 u = _ex1;
1027 c = p = _ex0;
1028 return;
1029 }
1030
1031 // Special case: input is a number
1032 if (is_exactly_a<numeric>(*this)) {
1034 u = _ex_1;
1035 c = abs(ex_to<numeric>(*this));
1036 } else {
1037 u = _ex1;
1038 c = *this;
1039 }
1040 p = _ex1;
1041 return;
1042 }
1043
1044 // Expand input polynomial
1045 ex e = expand();
1046 if (e.is_zero()) {
1047 u = _ex1;
1048 c = p = _ex0;
1049 return;
1050 }
1051
1052 // Compute unit and content
1053 u = unit(x);
1054 c = content(x);
1055
1056 // Divide by unit and content to get primitive part
1057 if (c.is_zero()) {
1058 p = _ex0;
1059 return;
1060 }
1061 if (is_exactly_a<numeric>(c))
1062 p = *this / (c * u);
1063 else
1064 p = quo(e, c * u, x, false);
1065}
1066
1067
1068/*
1069 * GCD of multivariate polynomials
1070 */
1071
1081static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
1082{
1083#if STATISTICS
1084 sr_gcd_called++;
1085#endif
1086
1087 // The first symbol is our main variable
1088 const ex &x = var->sym;
1089
1090 // Sort c and d so that c has higher degree
1091 ex c, d;
1092 int adeg = a.degree(x), bdeg = b.degree(x);
1093 int cdeg, ddeg;
1094 if (adeg >= bdeg) {
1095 c = a;
1096 d = b;
1097 cdeg = adeg;
1098 ddeg = bdeg;
1099 } else {
1100 c = b;
1101 d = a;
1102 cdeg = bdeg;
1103 ddeg = adeg;
1104 }
1105
1106 // Remove content from c and d, to be attached to GCD later
1107 ex cont_c = c.content(x);
1108 ex cont_d = d.content(x);
1109 ex gamma = gcd(cont_c, cont_d, nullptr, nullptr, false);
1110 if (ddeg == 0)
1111 return gamma;
1112 c = c.primpart(x, cont_c);
1113 d = d.primpart(x, cont_d);
1114
1115 // First element of subresultant sequence
1116 ex r = _ex0, ri = _ex1, psi = _ex1;
1117 int delta = cdeg - ddeg;
1118
1119 for (;;) {
1120
1121 // Calculate polynomial pseudo-remainder
1122 r = prem(c, d, x, false);
1123 if (r.is_zero())
1124 return gamma * d.primpart(x);
1125
1126 c = d;
1127 cdeg = ddeg;
1128 if (!divide_in_z(r, ri * pow(psi, delta), d, var))
1129 throw(std::runtime_error("invalid expression in sr_gcd(), division failed"));
1130 ddeg = d.degree(x);
1131 if (ddeg == 0) {
1132 if (is_exactly_a<numeric>(r))
1133 return gamma;
1134 else
1135 return gamma * r.primpart(x);
1136 }
1137
1138 // Next element of subresultant sequence
1139 ri = c.expand().lcoeff(x);
1140 if (delta == 1)
1141 psi = ri;
1142 else if (delta)
1143 divide_in_z(pow(ri, delta), pow(psi, delta-1), psi, var+1);
1144 delta = cdeg - ddeg;
1145 }
1146}
1147
1148
1155{
1156 return bp->max_coefficient();
1157}
1158
1162{
1163 return *_num1_p;
1164}
1165
1167{
1168 return abs(*this);
1169}
1170
1172{
1173 GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
1174 numeric cur_max = abs(ex_to<numeric>(overall_coeff));
1175 for (auto & it : seq) {
1176 numeric a;
1177 GINAC_ASSERT(!is_exactly_a<numeric>(it.rest));
1178 a = abs(ex_to<numeric>(it.coeff));
1179 if (a > cur_max)
1180 cur_max = a;
1181 }
1182 return cur_max;
1183}
1184
1186{
1187#ifdef DO_GINAC_ASSERT
1188 for (auto & it : seq) {
1189 GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(it)));
1190 }
1191#endif // def DO_GINAC_ASSERT
1192 GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
1193 return abs(ex_to<numeric>(overall_coeff));
1194}
1195
1196
1203ex basic::smod(const numeric &xi) const
1204{
1205 return *this;
1206}
1207
1208ex numeric::smod(const numeric &xi) const
1209{
1210 return GiNaC::smod(*this, xi);
1211}
1212
1213ex add::smod(const numeric &xi) const
1214{
1215 epvector newseq;
1216 newseq.reserve(seq.size()+1);
1217 for (auto & it : seq) {
1218 GINAC_ASSERT(!is_exactly_a<numeric>(it.rest));
1219 numeric coeff = GiNaC::smod(ex_to<numeric>(it.coeff), xi);
1220 if (!coeff.is_zero())
1221 newseq.push_back(expair(it.rest, coeff));
1222 }
1223 GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
1224 numeric coeff = GiNaC::smod(ex_to<numeric>(overall_coeff), xi);
1225 return dynallocate<add>(std::move(newseq), coeff);
1226}
1227
1228ex mul::smod(const numeric &xi) const
1229{
1230#ifdef DO_GINAC_ASSERT
1231 for (auto & it : seq) {
1232 GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(it)));
1233 }
1234#endif // def DO_GINAC_ASSERT
1235 mul & mulcopy = dynallocate<mul>(*this);
1236 GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
1237 mulcopy.overall_coeff = GiNaC::smod(ex_to<numeric>(overall_coeff),xi);
1240 return mulcopy;
1241}
1242
1243
1245static ex interpolate(const ex &gamma, const numeric &xi, const ex &x, int degree_hint = 1)
1246{
1247 exvector g; g.reserve(degree_hint);
1248 ex e = gamma;
1249 numeric rxi = xi.inverse();
1250 for (int i=0; !e.is_zero(); i++) {
1251 ex gi = e.smod(xi);
1252 g.push_back(gi * pow(x, i));
1253 e = (e - gi) * rxi;
1254 }
1255 return dynallocate<add>(g);
1256}
1257
1260
1277static bool heur_gcd_z(ex& res, const ex &a, const ex &b, ex *ca, ex *cb,
1278 sym_desc_vec::const_iterator var)
1279{
1280#if STATISTICS
1281 heur_gcd_called++;
1282#endif
1283
1284 // Algorithm only works for non-vanishing input polynomials
1285 if (a.is_zero() || b.is_zero())
1286 return false;
1287
1288 // GCD of two numeric values -> CLN
1289 if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b)) {
1290 numeric g = gcd(ex_to<numeric>(a), ex_to<numeric>(b));
1291 if (ca)
1292 *ca = ex_to<numeric>(a) / g;
1293 if (cb)
1294 *cb = ex_to<numeric>(b) / g;
1295 res = g;
1296 return true;
1297 }
1298
1299 // The first symbol is our main variable
1300 const ex &x = var->sym;
1301
1302 // Remove integer content
1304 numeric rgc = gc.inverse();
1305 ex p = a * rgc;
1306 ex q = b * rgc;
1307 int maxdeg = std::max(p.degree(x), q.degree(x));
1308
1309 // Find evaluation point
1310 numeric mp = p.max_coefficient();
1311 numeric mq = q.max_coefficient();
1312 numeric xi;
1313 if (mp > mq)
1314 xi = mq * (*_num2_p) + (*_num2_p);
1315 else
1316 xi = mp * (*_num2_p) + (*_num2_p);
1317
1318 // 6 tries maximum
1319 for (int t=0; t<6; t++) {
1320 if (xi.int_length() * maxdeg > 100000) {
1321 throw gcdheu_failed();
1322 }
1323
1324 // Apply evaluation homomorphism and calculate GCD
1325 ex cp, cq;
1326 ex gamma;
1327 bool found = heur_gcd_z(gamma,
1330 &cp, &cq, var+1);
1331 if (found) {
1332 gamma = gamma.expand();
1333 // Reconstruct polynomial from GCD of mapped polynomials
1334 ex g = interpolate(gamma, xi, x, maxdeg);
1335
1336 // Remove integer content
1337 g /= g.integer_content();
1338
1339 // If the calculated polynomial divides both p and q, this is the GCD
1340 ex dummy;
1341 if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) {
1342 g *= gc;
1343 res = g;
1344 return true;
1345 }
1346 }
1347
1348 // Next evaluation point
1349 xi = iquo(xi * isqrt(isqrt(xi)) * numeric(73794), numeric(27011));
1350 }
1351 return false;
1352}
1353
1371static bool heur_gcd(ex& res, const ex& a, const ex& b, ex *ca, ex *cb,
1372 sym_desc_vec::const_iterator var)
1373{
1376 try {
1377 return heur_gcd_z(res, a, b, ca, cb, var);
1378 } catch (gcdheu_failed) {
1379 return false;
1380 }
1381 }
1382
1383 // convert polynomials to Z[X]
1385 const numeric ab_lcm = lcmcoeff(b, a_lcm);
1386
1387 const ex ai = a*ab_lcm;
1388 const ex bi = b*ab_lcm;
1390 throw std::logic_error("heur_gcd: not an integer polynomial [1]");
1391
1393 throw std::logic_error("heur_gcd: not an integer polynomial [2]");
1394
1395 bool found = false;
1396 try {
1397 found = heur_gcd_z(res, ai, bi, ca, cb, var);
1398 } catch (gcdheu_failed) {
1399 return false;
1400 }
1401
1402 // GCD is not unique, it's defined up to a unit (i.e. invertible
1403 // element). If the coefficient ring is a field, every its element is
1404 // invertible, so one can multiply the polynomial GCD with any element
1405 // of the coefficient field. We use this ambiguity to make cofactors
1406 // integer polynomials.
1407 if (found)
1408 res /= ab_lcm;
1409 return found;
1410}
1411
1412
1413// gcd helper to handle partially factored polynomials (to avoid expanding
1414// large expressions). At least one of the arguments should be a power.
1415static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb);
1416
1417// gcd helper to handle partially factored polynomials (to avoid expanding
1418// large expressions). At least one of the arguments should be a product.
1419static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb);
1420
1432ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned options)
1433{
1434#if STATISTICS
1435 gcd_called++;
1436#endif
1437
1438 // GCD of numerics -> CLN
1439 if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b)) {
1440 numeric g = gcd(ex_to<numeric>(a), ex_to<numeric>(b));
1441 if (ca || cb) {
1442 if (g.is_zero()) {
1443 if (ca)
1444 *ca = _ex0;
1445 if (cb)
1446 *cb = _ex0;
1447 } else {
1448 if (ca)
1449 *ca = ex_to<numeric>(a) / g;
1450 if (cb)
1451 *cb = ex_to<numeric>(b) / g;
1452 }
1453 }
1454 return g;
1455 }
1456
1457 // Check arguments
1459 throw(std::invalid_argument("gcd: arguments must be polynomials over the rationals"));
1460 }
1461
1462 // Partially factored cases (to avoid expanding large expressions)
1464 if (is_exactly_a<mul>(a) || is_exactly_a<mul>(b))
1465 return gcd_pf_mul(a, b, ca, cb);
1466#if FAST_COMPARE
1467 if (is_exactly_a<power>(a) || is_exactly_a<power>(b))
1468 return gcd_pf_pow(a, b, ca, cb);
1469#endif
1470 }
1471
1472 // Some trivial cases
1473 ex aex = a.expand();
1474 if (aex.is_zero()) {
1475 if (ca)
1476 *ca = _ex0;
1477 if (cb)
1478 *cb = _ex1;
1479 return b;
1480 }
1481 ex bex = b.expand();
1482 if (bex.is_zero()) {
1483 if (ca)
1484 *ca = _ex1;
1485 if (cb)
1486 *cb = _ex0;
1487 return a;
1488 }
1489 if (aex.is_equal(_ex1) || bex.is_equal(_ex1)) {
1490 if (ca)
1491 *ca = a;
1492 if (cb)
1493 *cb = b;
1494 return _ex1;
1495 }
1496#if FAST_COMPARE
1497 if (a.is_equal(b)) {
1498 if (ca)
1499 *ca = _ex1;
1500 if (cb)
1501 *cb = _ex1;
1502 return a;
1503 }
1504#endif
1505
1506 if (is_a<symbol>(aex)) {
1507 if (! bex.subs(aex==_ex0, subs_options::no_pattern).is_zero()) {
1508 if (ca)
1509 *ca = a;
1510 if (cb)
1511 *cb = b;
1512 return _ex1;
1513 }
1514 }
1515
1516 if (is_a<symbol>(bex)) {
1517 if (! aex.subs(bex==_ex0, subs_options::no_pattern).is_zero()) {
1518 if (ca)
1519 *ca = a;
1520 if (cb)
1521 *cb = b;
1522 return _ex1;
1523 }
1524 }
1525
1526 if (is_exactly_a<numeric>(aex)) {
1527 numeric bcont = bex.integer_content();
1528 numeric g = gcd(ex_to<numeric>(aex), bcont);
1529 if (ca)
1530 *ca = ex_to<numeric>(aex)/g;
1531 if (cb)
1532 *cb = bex/g;
1533 return g;
1534 }
1535
1536 if (is_exactly_a<numeric>(bex)) {
1537 numeric acont = aex.integer_content();
1538 numeric g = gcd(ex_to<numeric>(bex), acont);
1539 if (ca)
1540 *ca = aex/g;
1541 if (cb)
1542 *cb = ex_to<numeric>(bex)/g;
1543 return g;
1544 }
1545
1546 // Gather symbol statistics
1547 sym_desc_vec sym_stats;
1548 get_symbol_stats(a, b, sym_stats);
1549
1550 // The symbol with least degree which is contained in both polynomials
1551 // is our main variable
1552 auto vari = sym_stats.begin();
1553 while ((vari != sym_stats.end()) &&
1554 (((vari->ldeg_b == 0) && (vari->deg_b == 0)) ||
1555 ((vari->ldeg_a == 0) && (vari->deg_a == 0))))
1556 vari++;
1557
1558 // No common symbols at all, just return 1:
1559 if (vari == sym_stats.end()) {
1560 // N.B: keep cofactors factored
1561 if (ca)
1562 *ca = a;
1563 if (cb)
1564 *cb = b;
1565 return _ex1;
1566 }
1567 // move symbol contained only in one of the polynomials to the end:
1568 rotate(sym_stats.begin(), vari, sym_stats.end());
1569
1570 sym_desc_vec::const_iterator var = sym_stats.begin();
1571 const ex &x = var->sym;
1572
1573 // Cancel trivial common factor
1574 int ldeg_a = var->ldeg_a;
1575 int ldeg_b = var->ldeg_b;
1576 int min_ldeg = std::min(ldeg_a,ldeg_b);
1577 if (min_ldeg > 0) {
1578 ex common = pow(x, min_ldeg);
1579 return gcd((aex / common).expand(), (bex / common).expand(), ca, cb, false) * common;
1580 }
1581
1582 // Try to eliminate variables
1583 if (var->deg_a == 0 && var->deg_b != 0 ) {
1584 ex bex_u, bex_c, bex_p;
1585 bex.unitcontprim(x, bex_u, bex_c, bex_p);
1586 ex g = gcd(aex, bex_c, ca, cb, false);
1587 if (cb)
1588 *cb *= bex_u * bex_p;
1589 return g;
1590 } else if (var->deg_b == 0 && var->deg_a != 0) {
1591 ex aex_u, aex_c, aex_p;
1592 aex.unitcontprim(x, aex_u, aex_c, aex_p);
1593 ex g = gcd(aex_c, bex, ca, cb, false);
1594 if (ca)
1595 *ca *= aex_u * aex_p;
1596 return g;
1597 }
1598
1599 // Try heuristic algorithm first, fall back to PRS if that failed
1600 ex g;
1602 bool found = heur_gcd(g, aex, bex, ca, cb, var);
1603 if (found) {
1604 // heur_gcd have already computed cofactors...
1605 if (g.is_equal(_ex1)) {
1606 // ... but we want to keep them factored if possible.
1607 if (ca)
1608 *ca = a;
1609 if (cb)
1610 *cb = b;
1611 }
1612 return g;
1613 }
1614#if STATISTICS
1615 else {
1616 heur_gcd_failed++;
1617 }
1618#endif
1619 }
1621 g = sr_gcd(aex, bex, var);
1622 } else {
1623 exvector vars;
1624 for (std::size_t n = sym_stats.size(); n-- != 0; )
1625 vars.push_back(sym_stats[n].sym);
1626 g = chinrem_gcd(aex, bex, vars);
1627 }
1628
1629 if (g.is_equal(_ex1)) {
1630 // Keep cofactors factored if possible
1631 if (ca)
1632 *ca = a;
1633 if (cb)
1634 *cb = b;
1635 } else {
1636 if (ca)
1637 divide(aex, g, *ca, false);
1638 if (cb)
1639 divide(bex, g, *cb, false);
1640 }
1641 return g;
1642}
1643
1644// gcd helper to handle partially factored polynomials (to avoid expanding
1645// large expressions). Both arguments should be powers.
1646static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb)
1647{
1648 ex p = a.op(0);
1649 const ex& exp_a = a.op(1);
1650 ex pb = b.op(0);
1651 const ex& exp_b = b.op(1);
1652
1653 // a = p^n, b = p^m, gcd = p^min(n, m)
1654 if (p.is_equal(pb)) {
1655 if (exp_a < exp_b) {
1656 if (ca)
1657 *ca = _ex1;
1658 if (cb)
1659 *cb = pow(p, exp_b - exp_a);
1660 return pow(p, exp_a);
1661 } else {
1662 if (ca)
1663 *ca = pow(p, exp_a - exp_b);
1664 if (cb)
1665 *cb = _ex1;
1666 return pow(p, exp_b);
1667 }
1668 }
1669
1670 ex p_co, pb_co;
1671 ex p_gcd = gcd(p, pb, &p_co, &pb_co, false);
1672 // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==> gcd(a,b) = 1
1673 if (p_gcd.is_equal(_ex1)) {
1674 if (ca)
1675 *ca = a;
1676 if (cb)
1677 *cb = b;
1678 return _ex1;
1679 }
1680
1681 // there are common factors:
1682 // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==>
1683 // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m
1684 if (exp_a < exp_b) {
1685 ex pg = gcd(pow(p_co, exp_a), pow(p_gcd, exp_b-exp_a)*pow(pb_co, exp_b), ca, cb, false);
1686 return pow(p_gcd, exp_a)*pg;
1687 } else {
1688 ex pg = gcd(pow(p_gcd, exp_a - exp_b)*pow(p_co, exp_a), pow(pb_co, exp_b), ca, cb, false);
1689 return pow(p_gcd, exp_b)*pg;
1690 }
1691}
1692
1693static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb)
1694{
1695 if (is_exactly_a<power>(a) && is_exactly_a<power>(b))
1696 return gcd_pf_pow_pow(a, b, ca, cb);
1697
1698 if (is_exactly_a<power>(b) && (! is_exactly_a<power>(a)))
1699 return gcd_pf_pow(b, a, cb, ca);
1700
1701 GINAC_ASSERT(is_exactly_a<power>(a));
1702
1703 ex p = a.op(0);
1704 const ex& exp_a = a.op(1);
1705 if (p.is_equal(b)) {
1706 // a = p^n, b = p, gcd = p
1707 if (ca)
1708 *ca = pow(p, exp_a - 1);
1709 if (cb)
1710 *cb = _ex1;
1711 return p;
1712 }
1713 if (is_a<symbol>(p)) {
1714 // Cancel trivial common factor
1715 int ldeg_a = ex_to<numeric>(exp_a).to_int();
1716 int ldeg_b = b.ldegree(p);
1717 int min_ldeg = std::min(ldeg_a, ldeg_b);
1718 if (min_ldeg > 0) {
1719 ex common = pow(p, min_ldeg);
1720 return gcd(pow(p, ldeg_a - min_ldeg), (b / common).expand(), ca, cb, false) * common;
1721 }
1722 }
1723
1724 ex p_co, bpart_co;
1725 ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
1726
1727 if (p_gcd.is_equal(_ex1)) {
1728 // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
1729 if (ca)
1730 *ca = a;
1731 if (cb)
1732 *cb = b;
1733 return _ex1;
1734 }
1735 // a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
1736 ex rg = gcd(pow(p_gcd, exp_a-1)*pow(p_co, exp_a), bpart_co, ca, cb, false);
1737 return p_gcd*rg;
1738}
1739
1740static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb)
1741{
1742 if (is_exactly_a<mul>(a) && is_exactly_a<mul>(b)
1743 && (b.nops() > a.nops()))
1744 return gcd_pf_mul(b, a, cb, ca);
1745
1746 if (is_exactly_a<mul>(b) && (!is_exactly_a<mul>(a)))
1747 return gcd_pf_mul(b, a, cb, ca);
1748
1749 GINAC_ASSERT(is_exactly_a<mul>(a));
1750 size_t num = a.nops();
1751 exvector g; g.reserve(num);
1752 exvector acc_ca; acc_ca.reserve(num);
1753 ex part_b = b;
1754 for (size_t i=0; i<num; i++) {
1755 ex part_ca, part_cb;
1756 g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, false));
1757 acc_ca.push_back(part_ca);
1758 part_b = part_cb;
1759 }
1760 if (ca)
1761 *ca = dynallocate<mul>(acc_ca);
1762 if (cb)
1763 *cb = part_b;
1764 return dynallocate<mul>(g);
1765}
1766
1774ex lcm(const ex &a, const ex &b, bool check_args)
1775{
1776 if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b))
1777 return lcm(ex_to<numeric>(a), ex_to<numeric>(b));
1779 throw(std::invalid_argument("lcm: arguments must be polynomials over the rationals"));
1780
1781 ex ca, cb;
1782 ex g = gcd(a, b, &ca, &cb, false);
1783 return ca * cb * g;
1784}
1785
1786
1787/*
1788 * Square-free factorization
1789 */
1790
1798static epvector sqrfree_yun(const ex &a, const symbol &x)
1799{
1800 ex w = a;
1801 ex z = w.diff(x);
1802 ex g = gcd(w, z);
1803 if (g.is_zero()) {
1804 // manifest zero or hidden zero
1805 return {};
1806 }
1807 if (g.is_equal(_ex1)) {
1808 // w(x) and w'(x) share no factors: w(x) is square-free
1809 return {expair(a, _ex1)};
1810 }
1811
1813 ex i = 0; // exponent
1814 do {
1815 w = quo(w, g, x);
1816 if (w.is_zero()) {
1817 // hidden zero
1818 break;
1819 }
1820 z = quo(z, g, x) - w.diff(x);
1821 i += 1;
1822 if (w.is_equal(x)) {
1823 // shortcut for x^n with n ∈ ℕ
1824 i += quo(z, w.diff(x), x);
1825 factors.push_back(expair(w, i));
1826 break;
1827 }
1828 g = gcd(w, z);
1829 if (!g.is_equal(_ex1)) {
1830 factors.push_back(expair(g, i));
1831 }
1832 } while (!z.is_zero());
1833
1834 // correct for lost factor
1835 // (being based on GCDs, Yun's algorithm only finds factors up to a unit)
1836 const ex lost_factor = quo(a, mul{factors}, x);
1837 if (lost_factor.is_equal(_ex1)) {
1838 // trivial lost factor
1839 return factors;
1840 }
1841 if (!factors.empty() && factors[0].coeff.is_equal(1)) {
1842 // multiply factor^1 with lost_factor
1843 factors[0].rest *= lost_factor;
1844 return factors;
1845 }
1846 // no factor^1: prepend lost_factor^1 to the results
1847 epvector results = {expair(lost_factor, 1)};
1848 std::move(factors.begin(), factors.end(), std::back_inserter(results));
1849 return results;
1850}
1851
1852
1888ex sqrfree(const ex &a, const lst &l)
1889{
1890 if (is_exactly_a<numeric>(a) ||
1891 is_a<symbol>(a)) // shortcuts
1892 return a;
1893
1894 // If no lst of variables to factorize in was specified we have to
1895 // invent one now. Maybe one can optimize here by reversing the order
1896 // or so, I don't know.
1897 lst args;
1898 if (l.nops()==0) {
1899 sym_desc_vec sdv;
1900 get_symbol_stats(a, _ex0, sdv);
1901 for (auto & it : sdv)
1902 args.append(it.sym);
1903 } else {
1904 args = l;
1905 }
1906
1907 // Find the symbol to factor in at this stage
1908 if (!is_a<symbol>(args.op(0)))
1909 throw (std::runtime_error("sqrfree(): invalid factorization variable"));
1910 const symbol &x = ex_to<symbol>(args.op(0));
1911
1912 // convert the argument from something in Q[X] to something in Z[X]
1914 const ex tmp = multiply_lcm(a, lcm);
1915
1916 // find the factors
1917 epvector factors = sqrfree_yun(tmp, x);
1918 if (factors.empty()) {
1919 // the polynomial was a hidden zero
1920 return _ex0;
1921 }
1922
1923 // remove symbol x and proceed recursively with the remaining symbols
1924 args.remove_first();
1925
1926 // recurse down the factors in remaining variables
1927 if (args.nops()>0) {
1928 for (auto & it : factors)
1929 it.rest = sqrfree(it.rest, args);
1930 }
1931
1932 // Done with recursion, now construct the final result
1933 ex result = mul(factors);
1934
1935 // Put in the rational overall factor again and return
1936 return result * lcm.inverse();
1937}
1938
1939
1947ex sqrfree_parfrac(const ex & a, const symbol & x)
1948{
1949 // Find numerator and denominator
1950 ex nd = numer_denom(a);
1951 ex numer = nd.op(0), denom = nd.op(1);
1952//std::clog << "numer = " << numer << ", denom = " << denom << std::endl;
1953
1954 // Convert N(x)/D(x) -> Q(x) + R(x)/D(x), so degree(R) < degree(D)
1955 ex red_poly = quo(numer, denom, x), red_numer = rem(numer, denom, x).expand();
1956//std::clog << "red_poly = " << red_poly << ", red_numer = " << red_numer << std::endl;
1957
1958 // Factorize denominator and compute cofactors
1959 epvector yun = sqrfree_yun(denom, x);
1960 exvector factor, cofac;
1961 size_t dim = 0;
1962 for (size_t i=0; i<yun.size(); i++) {
1963 numeric i_exponent = ex_to<numeric>(yun[i].coeff);
1964 for (size_t j=0; j<i_exponent; j++) {
1965 factor.push_back(pow(yun[i].rest, j+1));
1966 dim += degree(yun[i].rest, x);
1967 ex prod = _ex1;
1968 for (size_t k=0; k<yun.size(); k++) {
1969 if (yun[k].coeff == i_exponent)
1970 prod *= pow(yun[k].rest, i_exponent-1-j);
1971 else
1972 prod *= pow(yun[k].rest, yun[k].coeff);
1973 }
1974 cofac.push_back(prod.expand());
1975 }
1976 }
1977//std::clog << "factors : " << exprseq(factor) << std::endl;
1978//std::clog << "cofactors: " << exprseq(cofac) << std::endl;
1979
1980 // Construct linear system for decomposition
1981 matrix sys(dim, dim);
1982 matrix rhs(dim, 1);
1983 matrix vars(dim, 1);
1984 for (size_t i=0, n=0, f=0; i<yun.size(); i++) {
1985 size_t i_expo = to_int(ex_to<numeric>(yun[i].coeff));
1986 for (size_t j=0; j<i_expo; j++) {
1987 for (size_t k=0; k<size_t(degree(yun[i].rest, x)); k++) {
1988 GINAC_ASSERT(n < dim && f < factor.size());
1989
1990 // column n of coefficient matrix
1991 for (size_t r=0; r+k<dim; r++) {
1992 sys(r+k, n) = cofac[f].coeff(x, r);
1993 }
1994
1995 // element n of right hand side vector
1996 rhs(n, 0) = red_numer.coeff(x, n);
1997
1998 // element n of free variables vector
1999 vars(n, 0) = symbol();
2000
2001 n++;
2002 }
2003 f++;
2004 }
2005 }
2006//std::clog << "coeffs: " << sys << std::endl;
2007//std::clog << "rhs : " << rhs << std::endl;
2008
2009 // Solve resulting linear system and sum up decomposed fractions
2010 matrix sol = sys.solve(vars, rhs);
2011//std::clog << "sol : " << sol << std::endl;
2012 ex sum = red_poly;
2013 for (size_t i=0, n=0, f=0; i<yun.size(); i++) {
2014 size_t i_expo = to_int(ex_to<numeric>(yun[i].coeff));
2015 for (size_t j=0; j<i_expo; j++) {
2016 ex frac_numer = 0;
2017 for (size_t k=0; k<size_t(degree(yun[i].rest, x)); k++) {
2018 GINAC_ASSERT(n < dim && f < factor.size());
2019 frac_numer += sol(n, 0) * pow(x, k);
2020 n++;
2021 }
2022 sum += frac_numer / factor[f];
2023
2024 f++;
2025 }
2026 }
2027
2028 return sum;
2029}
2030
2031
2032/*
2033 * Normal form of rational functions
2034 */
2035
2036/*
2037 * Note: The internal normal() functions (= basic::normal() and overloaded
2038 * functions) all return lists of the form {numerator, denominator}. This
2039 * is to get around mul::eval()'s automatic expansion of numeric coefficients.
2040 * E.g. (a+b)/3 is automatically converted to a/3+b/3 but we want to keep
2041 * the information that (a+b) is the numerator and 3 is the denominator.
2042 */
2043
2044
2077static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup, lst & modifier)
2078{
2079 // Since the repl contains replaced expressions we should search for them
2080 ex e_replaced = e.subs(repl, subs_options::no_pattern);
2081
2082 // Expression already replaced? Then return the assigned symbol
2083 auto it = rev_lookup.find(e_replaced);
2084 if (it != rev_lookup.end())
2085 return it->second;
2086
2087 // The expression can be the base of substituted power, which requires a more careful search
2088 if (! is_a<numeric>(e_replaced))
2089 for (auto & it : repl)
2090 if (is_a<power>(it.second) && e_replaced.is_equal(it.second.op(0))) {
2091 ex degree = pow(it.second.op(1), _ex_1);
2092 if (is_a<numeric>(degree) && ex_to<numeric>(degree).is_integer())
2093 return pow(it.first, degree);
2094 }
2095
2096 // We treat powers and the exponent functions differently because
2097 // they can be rationalised more efficiently
2098 if (is_a<function>(e_replaced) && is_ex_the_function(e_replaced, exp)) {
2099 for (auto & it : repl) {
2100 if (is_a<function>(it.second) && is_ex_the_function(it.second, exp)) {
2101 ex ratio = normal(e_replaced.op(0) / it.second.op(0));
2102 if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).is_rational()) {
2103 // Different exponents can be treated as powers of the same basic equation
2104 if (ex_to<numeric>(ratio).is_integer()) {
2105 // If ratio is an integer then this is simply the power of the existing symbol.
2106 // std::clog << e_replaced << " is a " << ratio << " power of " << it.first << std::endl;
2107 return dynallocate<power>(it.first, ratio);
2108 } else {
2109 // otherwise we need to give the replacement pattern to change
2110 // the previous expression...
2111 ex es = dynallocate<symbol>();
2112 ex Num = numer(ratio);
2113 modifier.append(it.first == power(es, denom(ratio)));
2114 // std::clog << e_replaced << " is power " << Num << " and "
2115 // << it.first << " is power " << denom(ratio) << " of the common base "
2116 // << exp(e_replaced.op(0)/Num) << std::endl;
2117 // ... and modify the replacement tables
2118 rev_lookup.erase(it.second);
2119 rev_lookup.insert({exp(e_replaced.op(0)/Num), es});
2120 repl.erase(it.first);
2121 repl.insert({es, exp(e_replaced.op(0)/Num)});
2122 return dynallocate<power>(es, Num);
2123 }
2124 }
2125 }
2126 }
2127 } else if (is_a<power>(e_replaced) && !is_a<numeric>(e_replaced.op(0)) // We do not replace simple monomials like x^3 or sqrt(2)
2128 && ! (is_a<symbol>(e_replaced.op(0))
2129 && is_a<numeric>(e_replaced.op(1)) && ex_to<numeric>(e_replaced.op(1)).is_integer())) {
2130 for (auto & it : repl) {
2131 if (e_replaced.op(0).is_equal(it.second) // The base is an allocated symbol or base of power
2132 || (is_a<power>(it.second) && e_replaced.op(0).is_equal(it.second.op(0)))) {
2133 ex ratio; // We bind together two above cases
2134 if (is_a<power>(it.second))
2135 ratio = normal(e_replaced.op(1) / it.second.op(1));
2136 else
2137 ratio = e_replaced.op(1);
2138 if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).is_rational()) {
2139 // Different powers can be treated as powers of the same basic equation
2140 if (ex_to<numeric>(ratio).is_integer()) {
2141 // If ratio is an integer then this is simply the power of the existing symbol.
2142 //std::clog << e_replaced << " is a " << ratio << " power of " << it.first << std::endl;
2143 return dynallocate<power>(it.first, ratio);
2144 } else {
2145 // otherwise we need to give the replacement pattern to change
2146 // the previous expression...
2147 ex es = dynallocate<symbol>();
2148 ex Num = numer(ratio);
2149 modifier.append(it.first == power(es, denom(ratio)));
2150 //std::clog << e_replaced << " is power " << Num << " and "
2151 // << it.first << " is power " << denom(ratio) << " of the common base "
2152 // << pow(e_replaced.op(0), e_replaced.op(1)/Num) << std::endl;
2153 // ... and modify the replacement tables
2154 rev_lookup.erase(it.second);
2155 rev_lookup.insert({pow(e_replaced.op(0), e_replaced.op(1)/Num), es});
2156 repl.erase(it.first);
2157 repl.insert({es, pow(e_replaced.op(0), e_replaced.op(1)/Num)});
2158 return dynallocate<power>(es, Num);
2159 }
2160 }
2161 }
2162 }
2163 // There is no existing substitution, thus we are creating a new one.
2164 // This needs to be done separately to treat possible occurrences of
2165 // b = e_replaced.op(0) elsewhere in the expression as pow(b, 1).
2166 ex degree = pow(e_replaced.op(1), _ex_1);
2167 if (is_a<numeric>(degree) && ex_to<numeric>(degree).is_integer()) {
2168 ex es = dynallocate<symbol>();
2169 modifier.append(e_replaced.op(0) == power(es, degree));
2170 repl.insert({es, e_replaced});
2171 rev_lookup.insert({e_replaced, es});
2172 return es;
2173 }
2174 }
2175
2176 // Otherwise create new symbol and add to list, taking care that the
2177 // replacement expression doesn't itself contain symbols from repl,
2178 // because subs() is not recursive
2179 ex es = dynallocate<symbol>();
2180 repl.insert(std::make_pair(es, e_replaced));
2181 rev_lookup.insert(std::make_pair(e_replaced, es));
2182 return es;
2183}
2184
2190static ex replace_with_symbol(const ex & e, exmap & repl)
2191{
2192 // Since the repl contains replaced expressions we should search for them
2193 ex e_replaced = e.subs(repl, subs_options::no_pattern);
2194
2195 // Expression already replaced? Then return the assigned symbol
2196 for (auto & it : repl)
2197 if (it.second.is_equal(e_replaced))
2198 return it.first;
2199
2200 // Otherwise create new symbol and add to list, taking care that the
2201 // replacement expression doesn't itself contain symbols from repl,
2202 // because subs() is not recursive
2203 ex es = dynallocate<symbol>();
2204 repl.insert(std::make_pair(es, e_replaced));
2205 return es;
2206}
2207
2208
2211 ex operator()(const ex & e) override { return normal(e); }
2212};
2213
2217ex basic::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2218{
2219 if (nops() == 0)
2220 return dynallocate<lst>({replace_with_symbol(*this, repl, rev_lookup, modifier), _ex1});
2221
2222 normal_map_function map_normal;
2223 size_t nmod = modifier.nops(); // To watch new modifiers to the replacement list
2224 ex result = replace_with_symbol(map(map_normal), repl, rev_lookup, modifier);
2225 for (size_t imod = nmod; imod < modifier.nops(); ++imod) {
2226 exmap this_repl;
2227 this_repl.insert(std::make_pair(modifier.op(imod).op(0), modifier.op(imod).op(1)));
2228 result = result.subs(this_repl, subs_options::no_pattern);
2229 }
2230
2231 // Sometimes we may obtain negative powers, they need to be placed to denominator
2232 if (is_a<power>(result) && result.op(1).info(info_flags::negative))
2233 return dynallocate<lst>({_ex1, power(result.op(0), -result.op(1))});
2234 else
2235 return dynallocate<lst>({result, _ex1});
2236}
2237
2238
2241ex symbol::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2242{
2243 return dynallocate<lst>({*this, _ex1});
2244}
2245
2246
2251ex numeric::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2252{
2253 numeric num = numer();
2254 ex numex = num;
2255
2256 if (num.is_real()) {
2257 if (!num.is_integer())
2258 numex = replace_with_symbol(numex, repl, rev_lookup, modifier);
2259 } else { // complex
2260 numeric re = num.real(), im = num.imag();
2261 ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl, rev_lookup, modifier);
2262 ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl, rev_lookup, modifier);
2263 numex = re_ex + im_ex * replace_with_symbol(I, repl, rev_lookup, modifier);
2264 }
2265
2266 // Denominator is always a real integer (see numeric::denom())
2267 return dynallocate<lst>({numex, denom()});
2268}
2269
2270
2275static ex frac_cancel(const ex &n, const ex &d)
2276{
2277 ex num = n;
2278 ex den = d;
2279 numeric pre_factor = *_num1_p;
2280
2281//std::clog << "frac_cancel num = " << num << ", den = " << den << std::endl;
2282
2283 // Handle trivial case where denominator is 1
2284 if (den.is_equal(_ex1))
2285 return dynallocate<lst>({num, den});
2286
2287 // Handle special cases where numerator or denominator is 0
2288 if (num.is_zero())
2289 return dynallocate<lst>({num, _ex1});
2290 if (den.expand().is_zero())
2291 throw(std::overflow_error("frac_cancel: division by zero in frac_cancel"));
2292
2293 // Bring numerator and denominator to Z[X] by multiplying with
2294 // LCM of all coefficients' denominators
2297 num = multiply_lcm(num, num_lcm);
2298 den = multiply_lcm(den, den_lcm);
2299 pre_factor = den_lcm / num_lcm;
2300
2301 // Cancel GCD from numerator and denominator
2302 ex cnum, cden;
2303 if (gcd(num, den, &cnum, &cden, false) != _ex1) {
2304 num = cnum;
2305 den = cden;
2306 }
2307
2308 // Make denominator unit normal (i.e. coefficient of first symbol
2309 // as defined by get_first_symbol() is made positive)
2310 if (is_exactly_a<numeric>(den)) {
2311 if (ex_to<numeric>(den).is_negative()) {
2312 num *= _ex_1;
2313 den *= _ex_1;
2314 }
2315 } else {
2316 ex x;
2317 if (get_first_symbol(den, x)) {
2318 GINAC_ASSERT(is_exactly_a<numeric>(den.unit(x)));
2319 if (ex_to<numeric>(den.unit(x)).is_negative()) {
2320 num *= _ex_1;
2321 den *= _ex_1;
2322 }
2323 }
2324 }
2325
2326 // Return result as list
2327//std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << std::endl;
2328 return dynallocate<lst>({num * pre_factor.numer(), den * pre_factor.denom()});
2329}
2330
2331
2335ex add::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2336{
2337 // Normalize children and split each one into numerator and denominator
2338 exvector nums, dens;
2339 nums.reserve(seq.size()+1);
2340 dens.reserve(seq.size()+1);
2341 size_t nmod = modifier.nops(); // To watch new modifiers to the replacement list
2342 for (auto & it : seq) {
2343 ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, modifier);
2344 nums.push_back(n.op(0));
2345 dens.push_back(n.op(1));
2346 }
2347 ex n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, modifier);
2348 nums.push_back(n.op(0));
2349 dens.push_back(n.op(1));
2350 GINAC_ASSERT(nums.size() == dens.size());
2351
2352 // Now, nums is a vector of all numerators and dens is a vector of
2353 // all denominators
2354//std::clog << "add::normal uses " << nums.size() << " summands:\n";
2355
2356 // Add fractions sequentially
2357 auto num_it = nums.begin(), num_itend = nums.end();
2358 auto den_it = dens.begin(), den_itend = dens.end();
2359//std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl;
2360 for (size_t imod = nmod; imod < modifier.nops(); ++imod) {
2361 while (num_it != num_itend) {
2362 *num_it = num_it->subs(modifier.op(imod), subs_options::no_pattern);
2363 ++num_it;
2364 *den_it = den_it->subs(modifier.op(imod), subs_options::no_pattern);
2365 ++den_it;
2366 }
2367 // Reset iterators for the next round
2368 num_it = nums.begin();
2369 den_it = dens.begin();
2370 }
2371
2372 ex num = *num_it++, den = *den_it++;
2373 while (num_it != num_itend) {
2374//std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl;
2375 ex next_num = *num_it++, next_den = *den_it++;
2376
2377 // Trivially add sequences of fractions with identical denominators
2378 while ((den_it != den_itend) && next_den.is_equal(*den_it)) {
2379 next_num += *num_it;
2380 num_it++; den_it++;
2381 }
2382
2383 // Addition of two fractions, taking advantage of the fact that
2384 // the heuristic GCD algorithm computes the cofactors at no extra cost
2385 ex co_den1, co_den2;
2386 ex g = gcd(den, next_den, &co_den1, &co_den2, false);
2387 num = ((num * co_den2) + (next_num * co_den1)).expand();
2388 den *= co_den2; // this is the lcm(den, next_den)
2389 }
2390//std::clog << " common denominator = " << den << std::endl;
2391
2392 // Cancel common factors from num/den
2393 return frac_cancel(num, den);
2394}
2395
2396
2400ex mul::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2401{
2402 // Normalize children, separate into numerator and denominator
2403 exvector num; num.reserve(seq.size());
2404 exvector den; den.reserve(seq.size());
2405 ex n;
2406 size_t nmod = modifier.nops(); // To watch new modifiers to the replacement list
2407 for (auto & it : seq) {
2408 n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, modifier);
2409 num.push_back(n.op(0));
2410 den.push_back(n.op(1));
2411 }
2412 n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, modifier);
2413 num.push_back(n.op(0));
2414 den.push_back(n.op(1));
2415 auto num_it = num.begin(), num_itend = num.end();
2416 auto den_it = den.begin();
2417 for (size_t imod = nmod; imod < modifier.nops(); ++imod) {
2418 while (num_it != num_itend) {
2419 *num_it = num_it->subs(modifier.op(imod), subs_options::no_pattern);
2420 ++num_it;
2421 *den_it = den_it->subs(modifier.op(imod), subs_options::no_pattern);
2422 ++den_it;
2423 }
2424 num_it = num.begin();
2425 den_it = den.begin();
2426 }
2427
2428 // Perform fraction cancellation
2429 return frac_cancel(dynallocate<mul>(num), dynallocate<mul>(den));
2430}
2431
2432
2437ex power::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2438{
2439 // Normalize basis and exponent (exponent gets reassembled)
2440 size_t nmod = modifier.nops(); // To watch new modifiers to the replacement list
2441 ex n_basis = ex_to<basic>(basis).normal(repl, rev_lookup, modifier);
2442 for (size_t imod = nmod; imod < modifier.nops(); ++imod)
2443 n_basis = n_basis.subs(modifier.op(imod), subs_options::no_pattern);
2444
2445 nmod = modifier.nops();
2446 ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup, modifier);
2447 for (size_t imod = nmod; imod < modifier.nops(); ++imod)
2448 n_exponent = n_exponent.subs(modifier.op(imod), subs_options::no_pattern);
2449 n_exponent = n_exponent.op(0) / n_exponent.op(1);
2450
2451 if (n_exponent.info(info_flags::integer)) {
2452
2453 if (n_exponent.info(info_flags::positive)) {
2454
2455 // (a/b)^n -> {a^n, b^n}
2456 return dynallocate<lst>({pow(n_basis.op(0), n_exponent), pow(n_basis.op(1), n_exponent)});
2457
2458 } else if (n_exponent.info(info_flags::negative)) {
2459
2460 // (a/b)^-n -> {b^n, a^n}
2461 return dynallocate<lst>({pow(n_basis.op(1), -n_exponent), pow(n_basis.op(0), -n_exponent)});
2462 }
2463
2464 } else {
2465
2466 if (n_exponent.info(info_flags::positive)) {
2467
2468 // (a/b)^x -> {sym((a/b)^x), 1}
2469 return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup, modifier), _ex1});
2470
2471 } else if (n_exponent.info(info_flags::negative)) {
2472
2473 if (n_basis.op(1).is_equal(_ex1)) {
2474
2475 // a^-x -> {1, sym(a^x)}
2476 return dynallocate<lst>({_ex1, replace_with_symbol(pow(n_basis.op(0), -n_exponent), repl, rev_lookup, modifier)});
2477
2478 } else {
2479
2480 // (a/b)^-x -> {sym((b/a)^x), 1}
2481 return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(1) / n_basis.op(0), -n_exponent), repl, rev_lookup, modifier), _ex1});
2482 }
2483 }
2484 }
2485
2486 // (a/b)^x -> {sym((a/b)^x, 1}
2487 return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup, modifier), _ex1});
2488}
2489
2490
2494ex pseries::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2495{
2496 epvector newseq;
2497 for (auto & it : seq) {
2498 ex restexp = it.rest.normal();
2499 if (!restexp.is_zero())
2500 newseq.push_back(expair(restexp, it.coeff));
2501 }
2502 ex n = pseries(relational(var,point), std::move(newseq));
2503 return dynallocate<lst>({replace_with_symbol(n, repl, rev_lookup, modifier), _ex1});
2504}
2505
2506
2519{
2520 exmap repl, rev_lookup;
2521 lst modifier;
2522
2523 ex e = bp->normal(repl, rev_lookup, modifier);
2524 GINAC_ASSERT(is_a<lst>(e));
2525
2526 // Re-insert replaced symbols
2527 if (!repl.empty()) {
2528 for(size_t i=0; i < modifier.nops(); ++i)
2529 e = e.subs(modifier.op(i), subs_options::no_pattern);
2530 e = e.subs(repl, subs_options::no_pattern);
2531 }
2532
2533 // Convert {numerator, denominator} form back to fraction
2534 return e.op(0) / e.op(1);
2535}
2536
2544{
2545 exmap repl, rev_lookup;
2546 lst modifier;
2547
2548 ex e = bp->normal(repl, rev_lookup, modifier);
2549 GINAC_ASSERT(is_a<lst>(e));
2550
2551 // Re-insert replaced symbols
2552 if (repl.empty())
2553 return e.op(0);
2554 else {
2555 for(size_t i=0; i < modifier.nops(); ++i)
2556 e = e.subs(modifier.op(i), subs_options::no_pattern);
2557
2558 return e.op(0).subs(repl, subs_options::no_pattern);
2559 }
2560}
2561
2569{
2570 exmap repl, rev_lookup;
2571 lst modifier;
2572
2573 ex e = bp->normal(repl, rev_lookup, modifier);
2574 GINAC_ASSERT(is_a<lst>(e));
2575
2576 // Re-insert replaced symbols
2577 if (repl.empty())
2578 return e.op(1);
2579 else {
2580 for(size_t i=0; i < modifier.nops(); ++i)
2581 e = e.subs(modifier.op(i), subs_options::no_pattern);
2582
2583 return e.op(1).subs(repl, subs_options::no_pattern);
2584 }
2585}
2586
2594{
2595 exmap repl, rev_lookup;
2596 lst modifier;
2597
2598 ex e = bp->normal(repl, rev_lookup, modifier);
2599 GINAC_ASSERT(is_a<lst>(e));
2600
2601 // Re-insert replaced symbols
2602 if (repl.empty())
2603 return e;
2604 else {
2605 for(size_t i=0; i < modifier.nops(); ++i)
2606 e = e.subs(modifier.op(i), subs_options::no_pattern);
2607
2608 return e.subs(repl, subs_options::no_pattern);
2609 }
2610}
2611
2612
2627{
2628 return bp->to_rational(repl);
2629}
2630
2632{
2633 return bp->to_polynomial(repl);
2634}
2635
2639{
2640 return replace_with_symbol(*this, repl);
2641}
2642
2644{
2645 return replace_with_symbol(*this, repl);
2646}
2647
2648
2652{
2653 return *this;
2654}
2655
2659{
2660 return *this;
2661}
2662
2663
2668{
2669 if (is_real()) {
2670 if (!is_rational())
2671 return replace_with_symbol(*this, repl);
2672 } else { // complex
2673 numeric re = real();
2674 numeric im = imag();
2675 ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl);
2676 ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl);
2677 return re_ex + im_ex * replace_with_symbol(I, repl);
2678 }
2679 return *this;
2680}
2681
2686{
2687 if (is_real()) {
2688 if (!is_integer())
2689 return replace_with_symbol(*this, repl);
2690 } else { // complex
2691 numeric re = real();
2692 numeric im = imag();
2693 ex re_ex = re.is_integer() ? re : replace_with_symbol(re, repl);
2694 ex im_ex = im.is_integer() ? im : replace_with_symbol(im, repl);
2695 return re_ex + im_ex * replace_with_symbol(I, repl);
2696 }
2697 return *this;
2698}
2699
2700
2704{
2706 return pow(basis.to_rational(repl), exponent);
2707 else
2708 return replace_with_symbol(*this, repl);
2709}
2710
2714{
2716 return pow(basis.to_polynomial(repl), exponent);
2718 {
2719 ex basis_pref = collect_common_factors(basis);
2720 if (is_exactly_a<mul>(basis_pref) || is_exactly_a<power>(basis_pref)) {
2721 // (A*B)^n will be automagically transformed to A^n*B^n
2722 ex t = pow(basis_pref, exponent);
2723 return t.to_polynomial(repl);
2724 }
2725 else
2726 return pow(replace_with_symbol(pow(basis, _ex_1), repl), -exponent);
2727 }
2728 else
2729 return replace_with_symbol(*this, repl);
2730}
2731
2732
2735{
2736 epvector s;
2737 s.reserve(seq.size());
2738 for (auto & it : seq)
2739 s.push_back(split_ex_to_pair(recombine_pair_to_ex(it).to_rational(repl)));
2740
2741 ex oc = overall_coeff.to_rational(repl);
2742 if (oc.info(info_flags::numeric))
2743 return thisexpairseq(std::move(s), overall_coeff);
2744 else
2745 s.push_back(expair(oc, _ex1));
2746 return thisexpairseq(std::move(s), default_overall_coeff());
2747}
2748
2751{
2752 epvector s;
2753 s.reserve(seq.size());
2754 for (auto & it : seq)
2755 s.push_back(split_ex_to_pair(recombine_pair_to_ex(it).to_polynomial(repl)));
2756
2757 ex oc = overall_coeff.to_polynomial(repl);
2758 if (oc.info(info_flags::numeric))
2759 return thisexpairseq(std::move(s), overall_coeff);
2760 else
2761 s.push_back(expair(oc, _ex1));
2762 return thisexpairseq(std::move(s), default_overall_coeff());
2763}
2764
2765
2769static ex find_common_factor(const ex & e, ex & factor, exmap & repl)
2770{
2771 if (is_exactly_a<add>(e)) {
2772
2773 size_t num = e.nops();
2774 exvector terms; terms.reserve(num);
2775 ex gc;
2776
2777 // Find the common GCD
2778 for (size_t i=0; i<num; i++) {
2779 ex x = e.op(i).to_polynomial(repl);
2780
2781 if (is_exactly_a<add>(x) || is_exactly_a<mul>(x) || is_a<power>(x)) {
2782 ex f = 1;
2783 x = find_common_factor(x, f, repl);
2784 x *= f;
2785 }
2786
2787 if (gc.is_zero())
2788 gc = x;
2789 else
2790 gc = gcd(gc, x);
2791
2792 terms.push_back(x);
2793 }
2794
2795 if (gc.is_equal(_ex1))
2796 return e;
2797
2798 if (gc.is_zero())
2799 return _ex0;
2800
2801 // The GCD is the factor we pull out
2802 factor *= gc;
2803
2804 // Now divide all terms by the GCD
2805 for (size_t i=0; i<num; i++) {
2806 ex x;
2807
2808 // Try to avoid divide() because it expands the polynomial
2809 ex &t = terms[i];
2810 if (is_exactly_a<mul>(t)) {
2811 for (size_t j=0; j<t.nops(); j++) {
2812 if (t.op(j).is_equal(gc)) {
2813 exvector v; v.reserve(t.nops());
2814 for (size_t k=0; k<t.nops(); k++) {
2815 if (k == j)
2816 v.push_back(_ex1);
2817 else
2818 v.push_back(t.op(k));
2819 }
2820 t = dynallocate<mul>(v);
2821 goto term_done;
2822 }
2823 }
2824 }
2825
2826 divide(t, gc, x);
2827 t = x;
2828term_done: ;
2829 }
2830 return dynallocate<add>(terms);
2831
2832 } else if (is_exactly_a<mul>(e)) {
2833
2834 size_t num = e.nops();
2835 exvector v; v.reserve(num);
2836
2837 for (size_t i=0; i<num; i++)
2838 v.push_back(find_common_factor(e.op(i), factor, repl));
2839
2840 return dynallocate<mul>(v);
2841
2842 } else if (is_exactly_a<power>(e)) {
2843 const ex e_exp(e.op(1));
2844 if (e_exp.info(info_flags::integer)) {
2845 ex eb = e.op(0).to_polynomial(repl);
2846 ex factor_local(_ex1);
2847 ex pre_res = find_common_factor(eb, factor_local, repl);
2848 factor *= pow(factor_local, e_exp);
2849 return pow(pre_res, e_exp);
2850
2851 } else
2852 return e.to_polynomial(repl);
2853
2854 } else
2855 return e;
2856}
2857
2858
2862{
2863 if (is_exactly_a<add>(e) || is_exactly_a<mul>(e) || is_exactly_a<power>(e)) {
2864
2865 exmap repl;
2866 ex factor = 1;
2867 ex r = find_common_factor(e, factor, repl);
2869
2870 } else
2871 return e;
2872}
2873
2874
2877ex resultant(const ex & e1, const ex & e2, const ex & s)
2878{
2879 const ex ee1 = e1.expand();
2880 const ex ee2 = e2.expand();
2881 if (!ee1.info(info_flags::polynomial) ||
2883 throw(std::runtime_error("resultant(): arguments must be polynomials"));
2884
2885 const int h1 = ee1.degree(s);
2886 const int l1 = ee1.ldegree(s);
2887 const int h2 = ee2.degree(s);
2888 const int l2 = ee2.ldegree(s);
2889
2890 const int msize = h1 + h2;
2891 matrix m(msize, msize);
2892
2893 for (int l = h1; l >= l1; --l) {
2894 const ex e = ee1.coeff(s, l);
2895 for (int k = 0; k < h2; ++k)
2896 m(k, k+h1-l) = e;
2897 }
2898 for (int l = h2; l >= l2; --l) {
2899 const ex e = ee2.coeff(s, l);
2900 for (int k = 0; k < h1; ++k)
2901 m(k+h2, k+h2-l) = e;
2902 }
2903
2904 return m.determinant();
2905}
2906
2907
2908} // namespace GiNaC
Interface to GiNaC's sums of expressions.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
Interface to GiNaC's ABC.
ex coeff(const ex &s, int n=1) const override
Return coefficient of degree n in object s.
Definition: add.cpp:294
numeric integer_content() const override
Definition: normal.cpp:333
numeric max_coefficient() const override
Implementation ex::max_coefficient().
Definition: normal.cpp:1171
ex expand(unsigned options=0) const override
Expand expression, i.e.
Definition: add.cpp:572
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a sum.
Definition: normal.cpp:2335
ex recombine_pair_to_ex(const expair &p) const override
Form an ex out of an expair, using the corresponding semantics.
Definition: add.cpp:564
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition: normal.cpp:1213
virtual size_t nops() const
Number of operands/members.
Definition: basic.cpp:229
const basic & clearflag(unsigned f) const
Clear some status_flags.
Definition: basic.h:291
virtual numeric integer_content() const
Definition: normal.cpp:323
virtual numeric max_coefficient() const
Implementation ex::max_coefficient().
Definition: normal.cpp:1161
virtual ex to_polynomial(exmap &repl) const
Definition: normal.cpp:2643
virtual ex to_rational(exmap &repl) const
Default implementation of ex::to_rational().
Definition: normal.cpp:2638
virtual ex coeff(const ex &s, int n=1) const
Return coefficient of degree n in object s.
Definition: basic.cpp:337
virtual ex smod(const numeric &xi) const
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition: normal.cpp:1203
virtual ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const
Default implementation of ex::normal().
Definition: normal.cpp:2217
virtual ex map(map_function &f) const
Construct new expression by applying the specified function to all sub-expressions (one level only,...
Definition: basic.cpp:294
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
size_t nops() const override
Number of operands/members.
Definition: container.h:118
ex op(size_t i) const override
Return operand/member at position i.
Definition: container.h:295
container & remove_first()
Remove first element.
Definition: container.h:400
container & append(const ex &b)
Add element at back.
Definition: container.h:391
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
ex to_rational(exmap &repl) const
Rationalization of non-rational functions.
Definition: normal.cpp:2626
ex unit(const ex &x) const
Compute unit part (= sign of leading coefficient) of a multivariate polynomial in Q[x].
Definition: normal.cpp:923
numeric integer_content() const
Compute the integer content (= GCD of all numeric coefficients) of an expanded polynomial.
Definition: normal.cpp:318
ex primpart(const ex &x) const
Compute primitive part of a multivariate polynomial in Q[x].
Definition: normal.cpp:981
ex lcoeff(const ex &s) const
Definition: ex.h:176
const_iterator begin() const noexcept
Definition: ex.h:662
bool find(const ex &pattern, exset &found) const
Find all occurrences of a pattern.
Definition: ex.cpp:105
ex diff(const symbol &s, unsigned nth=1) const
Compute partial derivative of an expression.
Definition: ex.cpp:86
ex expand(unsigned options=0) const
Definition: ex.cpp:73
ex numer_denom() const
Get numerator and denominator of an expression.
Definition: normal.cpp:2593
bool is_equal(const ex &other) const
Definition: ex.h:345
ex normal() const
Normalization of rational functions.
Definition: normal.cpp:2518
int degree(const ex &s) const
Definition: ex.h:173
ptr< basic > bp
pointer to basic object managed by this
Definition: ex.h:251
numeric max_coefficient() const
Return maximum (absolute value) coefficient of a polynomial.
Definition: normal.cpp:1154
ex to_polynomial(exmap &repl) const
Definition: normal.cpp:2631
void unitcontprim(const ex &x, ex &u, ex &c, ex &p) const
Compute unit part, content part, and primitive part of a multivariate polynomial in Q[x].
Definition: normal.cpp:1022
size_t nops() const
Definition: ex.h:135
ex smod(const numeric &xi) const
Definition: ex.h:202
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:841
bool info(unsigned inf) const
Definition: ex.h:132
ex denom() const
Get denominator of an expression.
Definition: normal.cpp:2568
bool is_zero() const
Definition: ex.h:213
ex op(size_t i) const
Definition: ex.h:136
ex content(const ex &x) const
Compute content part (= unit normal GCD of all coefficients) of a multivariate polynomial in Q[x].
Definition: normal.cpp:945
int ldegree(const ex &s) const
Definition: ex.h:174
ex numer() const
Get numerator of an expression.
Definition: normal.cpp:2543
ex coeff(const ex &s, int n=1) const
Definition: ex.h:175
A pair of expressions.
Definition: expair.h:38
epvector seq
Definition: expairseq.h:127
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for expairseqs.
Definition: normal.cpp:2750
virtual ex default_overall_coeff() const
Definition: expairseq.cpp:574
virtual ex thisexpairseq(const epvector &v, const ex &oc, bool do_index_renaming=false) const
Create an object of this type.
Definition: expairseq.cpp:489
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for expairseqs.
Definition: normal.cpp:2734
virtual ex recombine_pair_to_ex(const expair &p) const
Form an ex out of an expair, using the corresponding semantics.
Definition: expairseq.cpp:559
virtual expair split_ex_to_pair(const ex &e) const
Form an expair from an ex, using the corresponding semantics.
Definition: expairseq.cpp:532
Exception thrown by heur_gcd() to signal failure.
Definition: normal.cpp:1259
@ integer_polynomial
Definition: flags.h:256
@ rational_polynomial
Definition: flags.h:258
Symbolic matrices.
Definition: matrix.h:38
matrix solve(const matrix &vars, const matrix &rhs, unsigned algo=solve_algo::automatic) const
Solve a linear system consisting of a m x n matrix and a m x p right hand side by applying an elimina...
Definition: matrix.cpp:995
Product of expressions.
Definition: mul.h:32
ex recombine_pair_to_ex(const expair &p) const override
Form an ex out of an expair, using the corresponding semantics.
Definition: mul.cpp:999
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition: normal.cpp:1228
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a product.
Definition: normal.cpp:2400
numeric integer_content() const override
Definition: normal.cpp:348
numeric max_coefficient() const override
Implementation ex::max_coefficient().
Definition: normal.cpp:1185
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Definition: numeric.h:82
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for a numeric.
Definition: normal.cpp:2667
numeric integer_content() const override
Definition: normal.cpp:328
bool is_rational() const
True if object is an exact rational number, may even be complex (denominator may be unity).
Definition: numeric.cpp:1201
const numeric real() const
Real part of a number.
Definition: numeric.cpp:1339
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for a numeric.
Definition: normal.cpp:2685
bool is_real() const
True if object is a real integer, rational or float (but not complex).
Definition: numeric.cpp:1208
const numeric numer() const
Numerator.
Definition: numeric.cpp:1356
bool is_integer() const
True if object is a non-complex integer.
Definition: numeric.cpp:1154
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a numeric.
Definition: normal.cpp:2251
const numeric denom() const
Denominator.
Definition: numeric.cpp:1387
const numeric imag() const
Imaginary part of a number.
Definition: numeric.cpp:1346
numeric max_coefficient() const override
Implementation ex::max_coefficient().
Definition: normal.cpp:1166
int int_length() const
Size in binary notation.
Definition: numeric.cpp:1418
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition: normal.cpp:1208
const numeric inverse() const
Inverse of a number.
Definition: numeric.cpp:1053
bool is_zero() const
True if object is zero.
Definition: numeric.cpp:1129
This class holds a two-component object, a basis and and exponent representing exponentiation.
Definition: power.h:39
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for powers.
Definition: normal.cpp:2437
ex basis
Definition: power.h:105
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for powers.
Definition: normal.cpp:2703
ex exponent
Definition: power.h:106
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for powers.
Definition: normal.cpp:2713
ex var
Series variable (holds a symbol)
Definition: pseries.h:118
pseries(const ex &rel_, const epvector &ops_)
Construct pseries from a vector of coefficients and powers.
Definition: pseries.cpp:71
epvector seq
Vector of {coefficient, power} pairs.
Definition: pseries.h:115
ex point
Expansion point.
Definition: pseries.h:121
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for pseries.
Definition: normal.cpp:2494
This class holds a relation consisting of two expressions and a logical relation between them.
Definition: relational.h:35
@ evaluated
.eval() has already done its job
Definition: flags.h:203
@ hash_calculated
.calchash() has already done its job
Definition: flags.h:205
@ no_pattern
disable pattern matching
Definition: flags.h:51
Basic CAS symbol.
Definition: symbol.h:39
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for symbols.
Definition: normal.cpp:2651
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for symbols.
Definition: normal.cpp:2241
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for symbols.
Definition: normal.cpp:2658
Interface to GiNaC's constant types and some special constants.
Interface to GiNaC's light-weight expression handles.
Interface to sequences of expression pairs.
upvec factors
Definition: factor.cpp:1430
unsigned options
Definition: factor.cpp:2475
vector< int > k
Definition: factor.cpp:1435
size_t n
Definition: factor.cpp:1432
size_t c
Definition: factor.cpp:757
ex x
Definition: factor.cpp:1610
size_t r
Definition: factor.cpp:757
mvec m
Definition: factor.cpp:758
ex cont
Definition: factor.cpp:2191
Interface to class signaling failure of operation.
#define is_ex_the_function(OBJ, FUNCNAME)
Definition: function.h:765
Interface to GiNaC's initially known functions.
Definition of GiNaC's lst.
Interface to symbolic matrices.
Interface to GiNaC's products of expressions.
Definition: add.cpp:38
const numeric I
Imaginary unit.
Definition: numeric.cpp:1433
ex denom(const ex &thisex)
Definition: ex.h:763
const numeric pow(const numeric &x, const numeric &y)
Definition: numeric.h:251
static ex replace_with_symbol(const ex &e, exmap &repl, exmap &rev_lookup, lst &modifier)
Create a symbol for replacing the expression "e" (or return a previously assigned symbol).
Definition: normal.cpp:2077
std::map< ex, ex, ex_is_less > exmap
Definition: basic.h:50
ex sqrfree(const ex &a, const lst &l)
Compute a square-free factorization of a multivariate polynomial in Q[X].
Definition: normal.cpp:1888
std::vector< expair > epvector
expair-vector
Definition: expairseq.h:33
ex resultant(const ex &e1, const ex &e2, const ex &s)
Resultant of two expressions e1,e2 with respect to symbol s.
Definition: normal.cpp:2877
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2320
static ex interpolate(const ex &gamma, const numeric &xi, const ex &x, int degree_hint=1)
xi-adic polynomial interpolation
Definition: normal.cpp:1245
static void add_symbol(const ex &s, sym_desc_vec &v)
Definition: normal.cpp:163
bool is_negative(const numeric &x)
Definition: numeric.h:269
const ex _ex1
Definition: utils.cpp:385
bool is_rational(const numeric &x)
Definition: numeric.h:290
ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned options)
Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X) and b(X) in Z[X].
Definition: normal.cpp:1432
static ex find_common_factor(const ex &e, ex &factor, exmap &repl)
Remove the common factor in the terms of a sum 'e' by calculating the GCD, and multiply it into the e...
Definition: normal.cpp:2769
function psi(const T1 &p1)
Definition: inifcns.h:165
ex rhs(const ex &thisex)
Definition: ex.h:832
ex prem(const ex &a, const ex &b, const ex &x, bool check_args)
Pseudo-remainder of polynomials a(x) and b(x) in Q[x].
Definition: normal.cpp:492
static ex gcd_pf_pow_pow(const ex &a, const ex &b, ex *ca, ex *cb)
Definition: normal.cpp:1646
static epvector sqrfree_yun(const ex &a, const symbol &x)
Compute square-free factorization of multivariate polynomial a(x) using Yun's algorithm.
Definition: normal.cpp:1798
const numeric exp(const numeric &x)
Exponential function.
Definition: numeric.cpp:1439
ex quo(const ex &a, const ex &b, const ex &x, bool check_args)
Quotient q(x) of polynomials a(x) and b(x) in Q[x].
Definition: normal.cpp:373
static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
Compute GCD of multivariate polynomials using the subresultant PRS algorithm.
Definition: normal.cpp:1081
static bool get_first_symbol(const ex &e, ex &x)
Return pointer to first symbol found in expression.
Definition: normal.cpp:95
const ex _ex_1
Definition: utils.cpp:352
const numeric smod(const numeric &a_, const numeric &b_)
Modulus (in symmetric representation).
Definition: numeric.cpp:2346
int degree(const ex &thisex, const ex &s)
Definition: ex.h:751
static ex gcd_pf_mul(const ex &a, const ex &b, ex *ca, ex *cb)
Definition: normal.cpp:1740
ex sqrfree_parfrac(const ex &a, const symbol &x)
Compute square-free partial fraction decomposition of rational function a(x).
Definition: normal.cpp:1947
ex lcm(const ex &a, const ex &b, bool check_args)
Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X].
Definition: normal.cpp:1774
const numeric iquo(const numeric &a, const numeric &b)
Numeric integer quotient.
Definition: numeric.cpp:2409
static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_iterator var)
Exact polynomial division of a(X) by b(X) in Z[X].
Definition: normal.cpp:757
const numeric isqrt(const numeric &x)
Integer numeric square root.
Definition: numeric.cpp:2487
ex normal(const ex &thisex)
Definition: ex.h:769
ex decomp_rational(const ex &a, const ex &x)
Decompose rational function a(x)=N(x)/D(x) into P(x)+n(x)/D(x) with degree(n, x) < degree(D,...
Definition: normal.cpp:472
static bool heur_gcd(ex &res, const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const_iterator var)
Compute GCD of multivariate polynomials using the heuristic GCD algorithm.
Definition: normal.cpp:1371
const numeric * _num1_p
Definition: utils.cpp:384
ex coeff(const ex &thisex, const ex &s, int n=1)
Definition: ex.h:757
static ex multiply_lcm(const ex &e, const numeric &lcm)
Bring polynomial from Q[X] to Z[X] by multiplying in the previously determined LCM of the coefficient...
Definition: normal.cpp:271
static bool heur_gcd_z(ex &res, const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const_iterator var)
Compute GCD of multivariate polynomials using the heuristic GCD algorithm.
Definition: normal.cpp:1277
ex numer(const ex &thisex)
Definition: ex.h:760
ex factor(const ex &poly, unsigned options)
Interface function to the outside world.
Definition: factor.cpp:2576
int to_int(const numeric &x)
Definition: numeric.h:302
ex collect_common_factors(const ex &e)
Collect common factors in sums.
Definition: normal.cpp:2861
static numeric lcm_of_coefficients_denominators(const ex &e)
Compute LCM of denominators of coefficients of a polynomial.
Definition: normal.cpp:261
bool is_integer(const numeric &x)
Definition: numeric.h:272
static numeric lcmcoeff(const ex &e, const numeric &l)
Definition: normal.cpp:231
static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
Collect statistical information about symbols in polynomials.
Definition: normal.cpp:197
ex rem(const ex &a, const ex &b, const ex &x, bool check_args)
Remainder r(x) of polynomials a(x) and b(x) in Q[x].
Definition: normal.cpp:423
ex sprem(const ex &a, const ex &b, const ex &x, bool check_args)
Sparse pseudo-remainder of polynomials a(x) and b(x) in Q[x].
Definition: normal.cpp:544
std::vector< sym_desc > sym_desc_vec
Definition: normal.cpp:160
static ex frac_cancel(const ex &n, const ex &d)
Fraction cancellation.
Definition: normal.cpp:2275
bool divide(const ex &a, const ex &b, ex &q, bool check_args)
Exact polynomial division of a(X) by b(X) in Q[X].
Definition: normal.cpp:595
static void collect_symbols(const ex &e, sym_desc_vec &v)
Definition: normal.cpp:173
const ex _ex0
Definition: utils.cpp:369
std::vector< ex > exvector
Definition: basic.h:48
ex numer_denom(const ex &thisex)
Definition: ex.h:766
static ex gcd_pf_pow(const ex &a, const ex &b, ex *ca, ex *cb)
Definition: normal.cpp:1693
const numeric * _num0_p
Definition: utils.cpp:367
ex expand(const ex &thisex, unsigned options=0)
Definition: ex.h:730
This file defines several functions that work on univariate and multivariate polynomials and rational...
Makes the interface to the underlying bignum package available.
Interface to GiNaC's overloaded operators.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
Interface to class for extended truncated power series.
Interface to relations between expressions.
@ use_sr_gcd
By default GiNaC uses modular GCD algorithm.
Definition: normal.h:61
@ no_part_factored
GiNaC tries to avoid expanding expressions when computing GCDs.
Definition: normal.h:55
@ no_heur_gcd
Usually GiNaC tries heuristic GCD first, because typically it's much faster than anything else.
Definition: normal.h:47
Function object for map().
Definition: basic.h:85
Function object to be applied by basic::normal().
Definition: normal.cpp:2210
ex operator()(const ex &e) override
Definition: normal.cpp:2211
This structure holds information about the highest and lowest degrees in which a symbol appears in tw...
Definition: normal.cpp:122
int max_deg
Maximum of deg_a and deg_b (Used for sorting)
Definition: normal.cpp:144
int ldeg_a
Lowest degree of symbol in polynomial "a".
Definition: normal.cpp:138
ex sym
Reference to symbol.
Definition: normal.cpp:129
int ldeg_b
Lowest degree of symbol in polynomial "b".
Definition: normal.cpp:141
bool operator<(const sym_desc &x) const
Comparison operator for sorting.
Definition: normal.cpp:150
int deg_b
Highest degree of symbol in polynomial "b".
Definition: normal.cpp:135
int deg_a
Highest degree of symbol in polynomial "a".
Definition: normal.cpp:132
size_t max_lcnops
Maximum number of terms of leading coefficient of symbol in both polynomials.
Definition: normal.cpp:147
sym_desc(const ex &s)
Initialize symbol, leave other variables uninitialized.
Definition: normal.cpp:124
Interface to GiNaC's symbolic objects.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...

This page is part of the GiNaC developer's reference. It was generated automatically by doxygen. For an introduction, see the tutorial.