44#include "polynomial/chinrem_gcd.h"
62#define USE_TRIAL_DIVISION 0
70static int gcd_called = 0;
71static int sr_gcd_called = 0;
72static int heur_gcd_called = 0;
73static int heur_gcd_failed = 0;
76static struct _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";
97 if (is_a<symbol>(e)) {
100 }
else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
101 for (
size_t i=0; i<e.
nops(); i++)
104 }
else if (is_exactly_a<power>(e)) {
166 if (it.sym.is_equal(s))
175 if (is_a<symbol>(e)) {
177 }
else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
178 for (
size_t i=0; i<e.
nops(); i++)
180 }
else if (is_exactly_a<power>(e)) {
201 for (
auto & it : v) {
202 int deg_a = a.
degree(it.sym);
203 int deg_b = b.
degree(it.sym);
206 it.max_deg = std::max(deg_a, deg_b);
211 std::sort(v.begin(), v.end());
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;
234 return lcm(ex_to<numeric>(e).
denom(), l);
235 else if (is_exactly_a<add>(e)) {
237 for (
size_t i=0; i<e.
nops(); i++)
240 }
else if (is_exactly_a<mul>(e)) {
242 for (
size_t i=0; i<e.
nops(); i++)
245 }
else if (is_exactly_a<power>(e)) {
246 if (is_a<symbol>(e.
op(0)))
277 if (is_exactly_a<mul>(e)) {
279 size_t num = e.
nops();
283 for (
size_t i=0; i<num; i++) {
288 v.push_back(
lcm / lcm_accum);
289 return dynallocate<mul>(v);
290 }
else if (is_exactly_a<add>(e)) {
292 size_t num = e.
nops();
295 for (
size_t i=0; i<num; i++)
297 return dynallocate<add>(v);
298 }
else if (is_exactly_a<power>(e)) {
299 if (!is_a<symbol>(e.
op(0))) {
302 numeric root_of_lcm =
lcm.power(ex_to<numeric>(e.
op(1)).inverse());
308 return dynallocate<mul>(e,
lcm);
320 return bp->integer_content();
336 for (
auto & it :
seq) {
339 c =
gcd(ex_to<numeric>(it.coeff).numer(),
c);
340 l =
lcm(ex_to<numeric>(it.coeff).denom(), l);
350#ifdef DO_GINAC_ASSERT
351 for (
auto & it :
seq) {
376 throw(std::overflow_error(
"quo: division by zero"));
377 if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b))
384 throw(std::invalid_argument(
"quo: arguments must be polynomials over the rationals"));
391 int rdeg =
r.degree(
x);
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;
400 if (!
divide(rcoeff, blcoeff, term,
false))
401 return dynallocate<fail>();
403 term *=
pow(
x, rdeg - bdeg);
410 return dynallocate<add>(v);
426 throw(std::overflow_error(
"rem: division by zero"));
427 if (is_exactly_a<numeric>(a)) {
428 if (is_exactly_a<numeric>(b))
438 throw(std::invalid_argument(
"rem: arguments must be polynomials over the rationals"));
445 int rdeg =
r.degree(
x);
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;
453 if (!
divide(rcoeff, blcoeff, term,
false))
454 return dynallocate<fail>();
456 term *=
pow(
x, rdeg - bdeg);
477 if (is_exactly_a<fail>(q))
495 throw(std::overflow_error(
"prem: division by zero"));
496 if (is_exactly_a<numeric>(a)) {
497 if (is_exactly_a<numeric>(b))
503 throw(std::invalid_argument(
"prem: arguments must be polynomials over the rationals"));
508 int rdeg =
r.degree(
x);
512 blcoeff = eb.
coeff(
x, bdeg);
516 eb -= blcoeff *
pow(
x, bdeg);
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();
527 r -= rlcoeff *
pow(
x, rdeg);
532 return pow(blcoeff, delta - i) *
r;
547 throw(std::overflow_error(
"prem: division by zero"));
548 if (is_exactly_a<numeric>(a)) {
549 if (is_exactly_a<numeric>(b))
555 throw(std::invalid_argument(
"prem: arguments must be polynomials over the rationals"));
560 int rdeg =
r.degree(
x);
564 blcoeff = eb.
coeff(
x, bdeg);
568 eb -= blcoeff *
pow(
x, bdeg);
572 while (rdeg >= bdeg && !
r.is_zero()) {
573 ex rlcoeff =
r.coeff(
x, rdeg);
574 ex term = (
pow(
x, rdeg - bdeg) * eb * rlcoeff).
expand();
578 r -= rlcoeff *
pow(
x, rdeg);
598 throw(std::overflow_error(
"divide: division by zero"));
603 if (is_exactly_a<numeric>(b)) {
606 }
else if (is_exactly_a<numeric>(a))
616 throw(std::invalid_argument(
"divide: arguments must be polynomials over the rationals"));
621 throw(std::invalid_argument(
"invalid expression in divide()"));
624 if (is_exactly_a<mul>(b)) {
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))
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))
647 if (is_exactly_a<mul>(a)) {
652 bool divisible_p =
false;
653 for (i=0; i < a.
nops(); ++i) {
654 if (
divide(a.
op(i), b, rem_i,
false)) {
661 resv.reserve(a.
nops());
662 for (
size_t j=0; j < a.
nops(); j++) {
664 resv.push_back(rem_i);
666 resv.push_back(a.
op(j));
668 q = dynallocate<mul>(resv);
671 }
else if (is_exactly_a<power>(a)) {
674 const ex& ab(a.
op(0));
675 int a_exp = ex_to<numeric>(a.
op(1)).to_int();
677 if (
divide(ab, b, rem_i,
false)) {
678 q = rem_i *
pow(ab, a_exp - 1);
697 int rdeg =
r.degree(
x);
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;
706 if (!
divide(rcoeff, blcoeff, term,
false))
708 term *=
pow(
x, rdeg - bdeg);
712 q = dynallocate<add>(v);
726typedef std::pair<ex, ex> ex2;
727typedef std::pair<ex, bool> exbool;
730 bool operator() (
const ex2 &p,
const ex2 &q)
const
732 int cmp = p.first.compare(q.first);
733 return ((cmp<0) || (!(cmp>0) && p.second.compare(q.second)<0));
737typedef std::map<ex2, exbool, ex2_less> ex2_exbool_remember;
761 throw(std::overflow_error(
"divide_in_z: division by zero"));
766 if (is_exactly_a<numeric>(a)) {
767 if (is_exactly_a<numeric>(b)) {
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;
790 if (is_exactly_a<power>(b)) {
791 const ex& bb(b.
op(0));
793 int exp_b = ex_to<numeric>(b.
op(1)).to_int();
794 for (
int i=exp_b; i>0; i--) {
802 if (is_exactly_a<mul>(b)) {
804 for (
const auto & it : b) {
816 const ex &
x = var->sym;
823#if USE_TRIAL_DIVISION
829 vector<numeric> alpha; alpha.reserve(adeg + 1);
833 for (i=0; i<=adeg; i++) {
841 alpha.push_back(point);
847 vector<numeric> rcp; rcp.reserve(adeg + 1);
849 for (
k=1;
k<=adeg;
k++) {
850 numeric product = alpha[
k] - alpha[0];
852 product *= alpha[
k] - alpha[i];
853 rcp.push_back(product.
inverse());
859 for (
k=1;
k<=adeg;
k++) {
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]);
868 for (
k=adeg-1;
k>=0;
k--)
869 c =
c * (
x - alpha[
k]) + v[
k];
871 if (
c.degree(
x) == (adeg - bdeg)) {
886 exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
887 while (rdeg >= bdeg) {
888 ex term, rcoeff =
r.coeff(
x, rdeg);
895 q = dynallocate<add>(v);
897 dr_remember[ex2(a, b)] = exbool(q,
true);
904 dr_remember[ex2(a, b)] = exbool(q,
false);
926 if (is_exactly_a<numeric>(
c))
933 throw(std::invalid_argument(
"invalid expression in unit()"));
947 if (is_exactly_a<numeric>(*
this))
958 int deg =
r.degree(
x);
964 int ldeg =
r.ldegree(
x);
968 for (
int i=ldeg; i<=deg; i++)
1001 if (is_exactly_a<numeric>(*
this))
1006 if (is_exactly_a<numeric>(
c))
1007 return *
this / (
c * u);
1009 return quo(*
this,
c * u,
x,
false);
1032 if (is_exactly_a<numeric>(*
this)) {
1035 c =
abs(ex_to<numeric>(*
this));
1061 if (is_exactly_a<numeric>(
c))
1062 p = *
this / (
c * u);
1064 p =
quo(e,
c * u,
x,
false);
1081static ex sr_gcd(
const ex &a,
const ex &b, sym_desc_vec::const_iterator var)
1088 const ex &
x = var->sym;
1107 ex cont_c =
c.content(
x);
1109 ex gamma =
gcd(cont_c, cont_d,
nullptr,
nullptr,
false);
1112 c =
c.primpart(
x, cont_c);
1117 int delta = cdeg - ddeg;
1129 throw(std::runtime_error(
"invalid expression in sr_gcd(), division failed"));
1132 if (is_exactly_a<numeric>(
r))
1135 return gamma *
r.primpart(
x);
1139 ri =
c.expand().lcoeff(
x);
1144 delta = cdeg - ddeg;
1156 return bp->max_coefficient();
1175 for (
auto & it :
seq) {
1178 a =
abs(ex_to<numeric>(it.coeff));
1187#ifdef DO_GINAC_ASSERT
1188 for (
auto & it :
seq) {
1216 newseq.reserve(
seq.size()+1);
1217 for (
auto & it :
seq) {
1225 return dynallocate<add>(std::move(newseq),
coeff);
1230#ifdef DO_GINAC_ASSERT
1231 for (
auto & it :
seq) {
1235 mul & mulcopy = dynallocate<mul>(*
this);
1247 exvector g; g.reserve(degree_hint);
1250 for (
int i=0; !e.
is_zero(); i++) {
1252 g.push_back(gi *
pow(
x, i));
1255 return dynallocate<add>(g);
1278 sym_desc_vec::const_iterator var)
1289 if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b)) {
1290 numeric g =
gcd(ex_to<numeric>(a), ex_to<numeric>(b));
1292 *ca = ex_to<numeric>(a) / g;
1294 *cb = ex_to<numeric>(b) / g;
1300 const ex &
x = var->sym;
1314 xi = mq * (*_num2_p) + (*_num2_p);
1316 xi = mp * (*_num2_p) + (*_num2_p);
1319 for (
int t=0; t<6; t++) {
1372 sym_desc_vec::const_iterator var)
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]");
1393 throw std::logic_error(
"heur_gcd: not an integer polynomial [2]");
1397 found =
heur_gcd_z(res, ai, bi, ca, cb, var);
1415static ex
gcd_pf_pow(
const ex& a,
const ex& b, ex* ca, ex* cb);
1419static ex
gcd_pf_mul(
const ex& a,
const ex& b, ex* ca, ex* cb);
1439 if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b)) {
1440 numeric g =
gcd(ex_to<numeric>(a), ex_to<numeric>(b));
1449 *ca = ex_to<numeric>(a) / g;
1451 *cb = ex_to<numeric>(b) / g;
1459 throw(std::invalid_argument(
"gcd: arguments must be polynomials over the rationals"));
1464 if (is_exactly_a<mul>(a) || is_exactly_a<mul>(b))
1467 if (is_exactly_a<power>(a) || is_exactly_a<power>(b))
1506 if (is_a<symbol>(aex)) {
1516 if (is_a<symbol>(bex)) {
1526 if (is_exactly_a<numeric>(aex)) {
1530 *ca = ex_to<numeric>(aex)/g;
1536 if (is_exactly_a<numeric>(bex)) {
1542 *cb = ex_to<numeric>(bex)/g;
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))))
1559 if (vari == sym_stats.end()) {
1568 rotate(sym_stats.begin(), vari, sym_stats.end());
1570 sym_desc_vec::const_iterator var = sym_stats.begin();
1571 const ex &
x = var->sym;
1574 int ldeg_a = var->ldeg_a;
1575 int ldeg_b = var->ldeg_b;
1576 int min_ldeg = std::min(ldeg_a,ldeg_b);
1578 ex common =
pow(
x, min_ldeg);
1579 return gcd((aex / common).
expand(), (bex / common).
expand(), ca, cb,
false) * common;
1583 if (var->deg_a == 0 && var->deg_b != 0 ) {
1584 ex bex_u, bex_c, bex_p;
1586 ex g =
gcd(aex, bex_c, ca, cb,
false);
1588 *cb *= bex_u * bex_p;
1590 }
else if (var->deg_b == 0 && var->deg_a != 0) {
1591 ex aex_u, aex_c, aex_p;
1593 ex g =
gcd(aex_c, bex, ca, cb,
false);
1595 *ca *= aex_u * aex_p;
1602 bool found =
heur_gcd(g, aex, bex, ca, cb, var);
1621 g =
sr_gcd(aex, bex, var);
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);
1637 divide(aex, g, *ca,
false);
1639 divide(bex, g, *cb,
false);
1649 const ex& exp_a = a.
op(1);
1651 const ex& exp_b = b.
op(1);
1655 if (exp_a < exp_b) {
1659 *cb =
pow(p, exp_b - exp_a);
1660 return pow(p, exp_a);
1663 *ca =
pow(p, exp_a - exp_b);
1666 return pow(p, exp_b);
1671 ex p_gcd =
gcd(p, pb, &p_co, &pb_co,
false);
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;
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;
1695 if (is_exactly_a<power>(a) && is_exactly_a<power>(b))
1698 if (is_exactly_a<power>(b) && (! is_exactly_a<power>(a)))
1704 const ex& exp_a = a.
op(1);
1708 *ca =
pow(p, exp_a - 1);
1713 if (is_a<symbol>(p)) {
1715 int ldeg_a = ex_to<numeric>(exp_a).to_int();
1717 int min_ldeg = std::min(ldeg_a, ldeg_b);
1719 ex common =
pow(p, min_ldeg);
1720 return gcd(
pow(p, ldeg_a - min_ldeg), (b / common).
expand(), ca, cb,
false) * common;
1725 ex p_gcd =
gcd(p, b, &p_co, &bpart_co,
false);
1736 ex rg =
gcd(
pow(p_gcd, exp_a-1)*
pow(p_co, exp_a), bpart_co, ca, cb,
false);
1742 if (is_exactly_a<mul>(a) && is_exactly_a<mul>(b)
1746 if (is_exactly_a<mul>(b) && (!is_exactly_a<mul>(a)))
1750 size_t num = a.
nops();
1752 exvector acc_ca; acc_ca.reserve(num);
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);
1761 *ca = dynallocate<mul>(acc_ca);
1764 return dynallocate<mul>(g);
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"));
1782 ex g =
gcd(a, b, &ca, &cb,
false);
1843 factors[0].rest *= lost_factor;
1848 std::move(
factors.begin(),
factors.end(), std::back_inserter(results));
1890 if (is_exactly_a<numeric>(a) ||
1901 for (
auto & it : sdv)
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));
1927 if (args.
nops()>0) {
1929 it.rest =
sqrfree(it.rest, args);
1936 return result *
lcm.inverse();
1962 for (
size_t i=0; i<yun.size(); i++) {
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);
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);
1974 cofac.push_back(prod.
expand());
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++) {
1991 for (
size_t r=0;
r+
k<dim;
r++) {
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++) {
2017 for (
size_t k=0;
k<size_t(
degree(yun[i].rest,
x));
k++) {
2019 frac_numer += sol(
n, 0) *
pow(
x,
k);
2022 sum += frac_numer /
factor[f];
2083 auto it = rev_lookup.
find(e_replaced);
2084 if (it != rev_lookup.end())
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))) {
2099 for (
auto & it : repl) {
2102 if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).
is_rational()) {
2107 return dynallocate<power>(it.first, ratio);
2111 ex es = dynallocate<symbol>();
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);
2127 }
else if (is_a<power>(e_replaced) && !is_a<numeric>(e_replaced.
op(0))
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) {
2132 || (is_a<power>(it.second) && e_replaced.
op(0).
is_equal(it.second.op(0)))) {
2134 if (is_a<power>(it.second))
2135 ratio =
normal(e_replaced.
op(1) / it.second.
op(1));
2137 ratio = e_replaced.
op(1);
2138 if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).
is_rational()) {
2143 return dynallocate<power>(it.first, ratio);
2147 ex es = dynallocate<symbol>();
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);
2168 ex es = dynallocate<symbol>();
2170 repl.insert({es, e_replaced});
2171 rev_lookup.insert({e_replaced, es});
2179 ex es = dynallocate<symbol>();
2180 repl.insert(std::make_pair(es, e_replaced));
2181 rev_lookup.insert(std::make_pair(e_replaced, es));
2196 for (
auto & it : repl)
2197 if (it.second.is_equal(e_replaced))
2203 ex es = dynallocate<symbol>();
2204 repl.insert(std::make_pair(es, e_replaced));
2223 size_t nmod = modifier.
nops();
2225 for (
size_t imod = nmod; imod < modifier.
nops(); ++imod) {
2227 this_repl.insert(std::make_pair(modifier.
op(imod).
op(0), modifier.
op(imod).
op(1)));
2233 return dynallocate<lst>({
_ex1,
power(result.
op(0), -result.
op(1))});
2235 return dynallocate<lst>({result,
_ex1});
2243 return dynallocate<lst>({*
this,
_ex1});
2267 return dynallocate<lst>({numex,
denom()});
2285 return dynallocate<lst>({num, den});
2289 return dynallocate<lst>({num,
_ex1});
2291 throw(std::overflow_error(
"frac_cancel: division by zero in frac_cancel"));
2299 pre_factor = den_lcm / num_lcm;
2303 if (
gcd(num, den, &cnum, &cden,
false) !=
_ex1) {
2310 if (is_exactly_a<numeric>(den)) {
2319 if (ex_to<numeric>(den.
unit(
x)).is_negative()) {
2328 return dynallocate<lst>({num * pre_factor.
numer(), den * pre_factor.
denom()});
2339 nums.reserve(
seq.size()+1);
2340 dens.reserve(
seq.size()+1);
2341 size_t nmod = modifier.
nops();
2342 for (
auto & it :
seq) {
2344 nums.push_back(
n.op(0));
2345 dens.push_back(
n.op(1));
2348 nums.push_back(
n.op(0));
2349 dens.push_back(
n.op(1));
2357 auto num_it = nums.begin(), num_itend = nums.end();
2358 auto den_it = dens.begin(), den_itend = dens.end();
2360 for (
size_t imod = nmod; imod < modifier.
nops(); ++imod) {
2361 while (num_it != num_itend) {
2368 num_it = nums.begin();
2369 den_it = dens.begin();
2372 ex num = *num_it++, den = *den_it++;
2373 while (num_it != num_itend) {
2375 ex next_num = *num_it++, next_den = *den_it++;
2378 while ((den_it != den_itend) && next_den.is_equal(*den_it)) {
2379 next_num += *num_it;
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();
2406 size_t nmod = modifier.
nops();
2407 for (
auto & it :
seq) {
2409 num.push_back(
n.op(0));
2410 den.push_back(
n.op(1));
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) {
2424 num_it = num.begin();
2425 den_it = den.begin();
2429 return frac_cancel(dynallocate<mul>(num), dynallocate<mul>(den));
2440 size_t nmod = modifier.
nops();
2441 ex n_basis = ex_to<basic>(
basis).
normal(repl, rev_lookup, modifier);
2442 for (
size_t imod = nmod; imod < modifier.
nops(); ++imod)
2445 nmod = modifier.
nops();
2447 for (
size_t imod = nmod; imod < modifier.
nops(); ++imod)
2449 n_exponent = n_exponent.
op(0) / n_exponent.
op(1);
2456 return dynallocate<lst>({
pow(n_basis.
op(0), n_exponent),
pow(n_basis.
op(1), n_exponent)});
2461 return dynallocate<lst>({
pow(n_basis.
op(1), -n_exponent),
pow(n_basis.
op(0), -n_exponent)});
2497 for (
auto & it :
seq) {
2520 exmap repl, rev_lookup;
2523 ex e =
bp->normal(repl, rev_lookup, modifier);
2527 if (!repl.empty()) {
2528 for(
size_t i=0; i < modifier.
nops(); ++i)
2534 return e.
op(0) / e.
op(1);
2545 exmap repl, rev_lookup;
2548 ex e =
bp->normal(repl, rev_lookup, modifier);
2555 for(
size_t i=0; i < modifier.
nops(); ++i)
2570 exmap repl, rev_lookup;
2573 ex e =
bp->normal(repl, rev_lookup, modifier);
2580 for(
size_t i=0; i < modifier.
nops(); ++i)
2595 exmap repl, rev_lookup;
2598 ex e =
bp->normal(repl, rev_lookup, modifier);
2605 for(
size_t i=0; i < modifier.
nops(); ++i)
2628 return bp->to_rational(repl);
2633 return bp->to_polynomial(repl);
2720 if (is_exactly_a<mul>(basis_pref) || is_exactly_a<power>(basis_pref)) {
2737 s.reserve(
seq.size());
2738 for (
auto & it :
seq)
2753 s.reserve(
seq.size());
2754 for (
auto & it :
seq)
2771 if (is_exactly_a<add>(e)) {
2773 size_t num = e.
nops();
2774 exvector terms; terms.reserve(num);
2778 for (
size_t i=0; i<num; i++) {
2781 if (is_exactly_a<add>(
x) || is_exactly_a<mul>(
x) || is_a<power>(
x)) {
2805 for (
size_t i=0; i<num; i++) {
2810 if (is_exactly_a<mul>(t)) {
2811 for (
size_t j=0; j<t.
nops(); j++) {
2814 for (
size_t k=0;
k<t.
nops();
k++) {
2818 v.push_back(t.
op(
k));
2820 t = dynallocate<mul>(v);
2830 return dynallocate<add>(terms);
2832 }
else if (is_exactly_a<mul>(e)) {
2834 size_t num = e.
nops();
2837 for (
size_t i=0; i<num; i++)
2840 return dynallocate<mul>(v);
2842 }
else if (is_exactly_a<power>(e)) {
2843 const ex e_exp(e.
op(1));
2849 return pow(pre_res, e_exp);
2863 if (is_exactly_a<add>(e) || is_exactly_a<mul>(e) || is_exactly_a<power>(e)) {
2883 throw(std::runtime_error(
"resultant(): arguments must be polynomials"));
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);
2890 const int msize = h1 + h2;
2893 for (
int l = h1; l >= l1; --l) {
2895 for (
int k = 0;
k < h2; ++
k)
2898 for (
int l = h2; l >= l2; --l) {
2900 for (
int k = 0;
k < h1; ++
k)
2901 m(
k+h2,
k+h2-l) = e;
2904 return m.determinant();
Interface to GiNaC's sums of expressions.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Interface to GiNaC's ABC.
ex coeff(const ex &s, int n=1) const override
Return coefficient of degree n in object s.
numeric integer_content() const override
numeric max_coefficient() const override
Implementation ex::max_coefficient().
ex expand(unsigned options=0) const override
Expand expression, i.e.
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a sum.
ex recombine_pair_to_ex(const expair &p) const override
Form an ex out of an expair, using the corresponding semantics.
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
virtual size_t nops() const
Number of operands/members.
const basic & clearflag(unsigned f) const
Clear some status_flags.
virtual numeric integer_content() const
virtual numeric max_coefficient() const
Implementation ex::max_coefficient().
virtual ex to_polynomial(exmap &repl) const
virtual ex to_rational(exmap &repl) const
Default implementation of ex::to_rational().
virtual ex coeff(const ex &s, int n=1) const
Return coefficient of degree n in object s.
virtual ex smod(const numeric &xi) const
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
virtual ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const
Default implementation of ex::normal().
virtual ex map(map_function &f) const
Construct new expression by applying the specified function to all sub-expressions (one level only,...
Wrapper template for making GiNaC classes out of STL containers.
size_t nops() const override
Number of operands/members.
ex op(size_t i) const override
Return operand/member at position i.
container & remove_first()
Remove first element.
container & append(const ex &b)
Add element at back.
Lightweight wrapper for GiNaC's symbolic objects.
ex to_rational(exmap &repl) const
Rationalization of non-rational functions.
ex unit(const ex &x) const
Compute unit part (= sign of leading coefficient) of a multivariate polynomial in Q[x].
numeric integer_content() const
Compute the integer content (= GCD of all numeric coefficients) of an expanded polynomial.
ex primpart(const ex &x) const
Compute primitive part of a multivariate polynomial in Q[x].
ex lcoeff(const ex &s) const
const_iterator begin() const noexcept
bool find(const ex &pattern, exset &found) const
Find all occurrences of a pattern.
ex diff(const symbol &s, unsigned nth=1) const
Compute partial derivative of an expression.
ex expand(unsigned options=0) const
ex numer_denom() const
Get numerator and denominator of an expression.
bool is_equal(const ex &other) const
ex normal() const
Normalization of rational functions.
int degree(const ex &s) const
ptr< basic > bp
pointer to basic object managed by this
numeric max_coefficient() const
Return maximum (absolute value) coefficient of a polynomial.
ex to_polynomial(exmap &repl) const
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].
ex smod(const numeric &xi) const
ex subs(const exmap &m, unsigned options=0) const
bool info(unsigned inf) const
ex denom() const
Get denominator of an expression.
ex content(const ex &x) const
Compute content part (= unit normal GCD of all coefficients) of a multivariate polynomial in Q[x].
int ldegree(const ex &s) const
ex numer() const
Get numerator of an expression.
ex coeff(const ex &s, int n=1) const
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for expairseqs.
virtual ex default_overall_coeff() const
virtual ex thisexpairseq(const epvector &v, const ex &oc, bool do_index_renaming=false) const
Create an object of this type.
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for expairseqs.
virtual ex recombine_pair_to_ex(const expair &p) const
Form an ex out of an expair, using the corresponding semantics.
virtual expair split_ex_to_pair(const ex &e) const
Form an expair from an ex, using the corresponding semantics.
Exception thrown by heur_gcd() to signal failure.
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...
ex recombine_pair_to_ex(const expair &p) const override
Form an ex out of an expair, using the corresponding semantics.
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a product.
numeric integer_content() const override
numeric max_coefficient() const override
Implementation ex::max_coefficient().
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for a numeric.
numeric integer_content() const override
bool is_rational() const
True if object is an exact rational number, may even be complex (denominator may be unity).
const numeric real() const
Real part of a number.
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for a numeric.
bool is_real() const
True if object is a real integer, rational or float (but not complex).
const numeric numer() const
Numerator.
bool is_integer() const
True if object is a non-complex integer.
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a numeric.
const numeric denom() const
Denominator.
const numeric imag() const
Imaginary part of a number.
numeric max_coefficient() const override
Implementation ex::max_coefficient().
int int_length() const
Size in binary notation.
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
const numeric inverse() const
Inverse of a number.
bool is_zero() const
True if object is zero.
This class holds a two-component object, a basis and and exponent representing exponentiation.
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for powers.
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for powers.
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for powers.
ex var
Series variable (holds a symbol)
pseries(const ex &rel_, const epvector &ops_)
Construct pseries from a vector of coefficients and powers.
epvector seq
Vector of {coefficient, power} pairs.
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for pseries.
This class holds a relation consisting of two expressions and a logical relation between them.
@ evaluated
.eval() has already done its job
@ hash_calculated
.calchash() has already done its job
@ no_pattern
disable pattern matching
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for symbols.
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for symbols.
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for symbols.
Interface to GiNaC's constant types and some special constants.
Interface to GiNaC's light-weight expression handles.
Interface to sequences of expression pairs.
Interface to class signaling failure of operation.
#define is_ex_the_function(OBJ, FUNCNAME)
Interface to GiNaC's initially known functions.
Definition of GiNaC's lst.
Interface to symbolic matrices.
Interface to GiNaC's products of expressions.
const numeric I
Imaginary unit.
ex denom(const ex &thisex)
const numeric pow(const numeric &x, const numeric &y)
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).
std::map< ex, ex, ex_is_less > exmap
ex sqrfree(const ex &a, const lst &l)
Compute a square-free factorization of a multivariate polynomial in Q[X].
std::vector< expair > epvector
expair-vector
ex resultant(const ex &e1, const ex &e2, const ex &s)
Resultant of two expressions e1,e2 with respect to symbol s.
const numeric abs(const numeric &x)
Absolute value.
static ex interpolate(const ex &gamma, const numeric &xi, const ex &x, int degree_hint=1)
xi-adic polynomial interpolation
static void add_symbol(const ex &s, sym_desc_vec &v)
bool is_negative(const numeric &x)
bool is_rational(const numeric &x)
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].
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...
function psi(const T1 &p1)
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].
static ex gcd_pf_pow_pow(const ex &a, const ex &b, ex *ca, ex *cb)
static epvector sqrfree_yun(const ex &a, const symbol &x)
Compute square-free factorization of multivariate polynomial a(x) using Yun's algorithm.
const numeric exp(const numeric &x)
Exponential function.
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].
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.
static bool get_first_symbol(const ex &e, ex &x)
Return pointer to first symbol found in expression.
const numeric smod(const numeric &a_, const numeric &b_)
Modulus (in symmetric representation).
int degree(const ex &thisex, const ex &s)
static ex gcd_pf_mul(const ex &a, const ex &b, ex *ca, ex *cb)
ex sqrfree_parfrac(const ex &a, const symbol &x)
Compute square-free partial fraction decomposition of rational function a(x).
ex lcm(const ex &a, const ex &b, bool check_args)
Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X].
const numeric iquo(const numeric &a, const numeric &b)
Numeric integer quotient.
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].
const numeric isqrt(const numeric &x)
Integer numeric square root.
ex normal(const ex &thisex)
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,...
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.
ex coeff(const ex &thisex, const ex &s, int n=1)
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...
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.
ex numer(const ex &thisex)
ex factor(const ex &poly, unsigned options)
Interface function to the outside world.
int to_int(const numeric &x)
ex collect_common_factors(const ex &e)
Collect common factors in sums.
static numeric lcm_of_coefficients_denominators(const ex &e)
Compute LCM of denominators of coefficients of a polynomial.
bool is_integer(const numeric &x)
static numeric lcmcoeff(const ex &e, const numeric &l)
static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
Collect statistical information about symbols in polynomials.
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].
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].
std::vector< sym_desc > sym_desc_vec
static ex frac_cancel(const ex &n, const ex &d)
Fraction cancellation.
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].
static void collect_symbols(const ex &e, sym_desc_vec &v)
std::vector< ex > exvector
ex numer_denom(const ex &thisex)
static ex gcd_pf_pow(const ex &a, const ex &b, ex *ca, ex *cb)
ex expand(const ex &thisex, unsigned options=0)
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.
@ no_part_factored
GiNaC tries to avoid expanding expressions when computing GCDs.
@ no_heur_gcd
Usually GiNaC tries heuristic GCD first, because typically it's much faster than anything else.
Function object for map().
Function object to be applied by basic::normal().
ex operator()(const ex &e) override
This structure holds information about the highest and lowest degrees in which a symbol appears in tw...
int max_deg
Maximum of deg_a and deg_b (Used for sorting)
int ldeg_a
Lowest degree of symbol in polynomial "a".
ex sym
Reference to symbol.
int ldeg_b
Lowest degree of symbol in polynomial "b".
bool operator<(const sym_desc &x) const
Comparison operator for sorting.
int deg_b
Highest degree of symbol in polynomial "b".
int deg_a
Highest degree of symbol in polynomial "a".
size_t max_lcnops
Maximum number of terms of leading coefficient of symbol in both polynomials.
sym_desc(const ex &s)
Initialize symbol, leave other variables uninitialized.
Interface to GiNaC's symbolic objects.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...