Go to the documentation of this file.
44 #include "polynomial/chinrem_gcd.h"
54 #define FAST_COMPARE 1
57 #define USE_REMEMBER 0
62 #define USE_TRIAL_DIVISION 0
70 static int gcd_called = 0;
71 static int sr_gcd_called = 0;
72 static int heur_gcd_called = 0;
73 static int heur_gcd_failed = 0;
76 static 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) {
354 #endif // def DO_GINAC_ASSERT
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);
726 typedef std::pair<ex, ex> ex2;
727 typedef 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));
737 typedef 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++)
969 cont =
gcd(
r.coeff(
x, i), cont,
nullptr,
nullptr,
false);
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);
1081 static 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) {
1191 #endif // def DO_GINAC_ASSERT
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) {
1234 #endif // def DO_GINAC_ASSERT
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);
1415 static ex
gcd_pf_pow(
const ex& a,
const ex& b, ex* ca, ex* cb);
1419 static 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();
1960 size_t yun_max_exponent = yun.empty() ? 0 : ex_to<numeric>(yun.back().coeff).to_int();
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));
1967 for (
size_t k=0;
k<yun.size();
k++) {
1968 if (yun[
k].
coeff == i_exponent)
1969 prod *=
pow(yun[
k].rest, i_exponent-1-j);
1973 cofac.push_back(prod.
expand());
1976 size_t num_factors =
factor.size();
1982 matrix sys(max_denom_deg + 1, num_factors);
1984 for (
int i=0; i<=max_denom_deg; i++) {
1985 for (
size_t j=0; j<num_factors; j++)
1986 sys(i, j) = cofac[j].
coeff(
x, i);
1993 matrix vars(num_factors, 1);
1994 for (
size_t i=0; i<num_factors; i++)
2000 for (
size_t i=0; i<num_factors; i++)
2001 sum += sol(i, 0) /
factor[i];
2003 return red_poly + sum;
2058 auto it = rev_lookup.
find(e_replaced);
2059 if (it != rev_lookup.end())
2063 if (! is_a<numeric>(e_replaced))
2064 for (
auto & it : repl)
2065 if (is_a<power>(it.second) && e_replaced.
is_equal(it.second.op(0))) {
2074 for (
auto & it : repl) {
2077 if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).
is_rational()) {
2082 return dynallocate<power>(it.first, ratio);
2086 ex es = dynallocate<symbol>();
2093 rev_lookup.erase(it.second);
2094 rev_lookup.insert({
exp(e_replaced.
op(0)/Num), es});
2095 repl.erase(it.first);
2096 repl.insert({es,
exp(e_replaced.
op(0)/Num)});
2097 return dynallocate<power>(es, Num);
2102 }
else if (is_a<power>(e_replaced) && !is_a<numeric>(e_replaced.
op(0))
2103 && ! (is_a<symbol>(e_replaced.
op(0))
2104 && is_a<numeric>(e_replaced.
op(1)) && ex_to<numeric>(e_replaced.
op(1)).is_integer())) {
2105 for (
auto & it : repl) {
2107 || (is_a<power>(it.second) && e_replaced.
op(0).
is_equal(it.second.op(0)))) {
2109 if (is_a<power>(it.second))
2110 ratio =
normal(e_replaced.
op(1) / it.second.
op(1));
2112 ratio = e_replaced.
op(1);
2113 if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).
is_rational()) {
2118 return dynallocate<power>(it.first, ratio);
2122 ex es = dynallocate<symbol>();
2129 rev_lookup.erase(it.second);
2130 rev_lookup.insert({
pow(e_replaced.
op(0), e_replaced.
op(1)/Num), es});
2131 repl.erase(it.first);
2132 repl.insert({es,
pow(e_replaced.
op(0), e_replaced.
op(1)/Num)});
2133 return dynallocate<power>(es, Num);
2143 ex es = dynallocate<symbol>();
2145 repl.insert({es, e_replaced});
2146 rev_lookup.insert({e_replaced, es});
2154 ex es = dynallocate<symbol>();
2155 repl.insert(std::make_pair(es, e_replaced));
2156 rev_lookup.insert(std::make_pair(e_replaced, es));
2171 for (
auto & it : repl)
2172 if (it.second.is_equal(e_replaced))
2178 ex es = dynallocate<symbol>();
2179 repl.insert(std::make_pair(es, e_replaced));
2198 int nmod = modifier.
nops();
2200 for (
int imod = nmod; imod < modifier.
nops(); ++imod) {
2202 this_repl.insert(std::make_pair(modifier.
op(imod).
op(0), modifier.
op(imod).
op(1)));
2214 return dynallocate<lst>({*
this,
_ex1});
2238 return dynallocate<lst>({numex,
denom()});
2256 return dynallocate<lst>({num, den});
2260 return dynallocate<lst>({num,
_ex1});
2262 throw(std::overflow_error(
"frac_cancel: division by zero in frac_cancel"));
2270 pre_factor = den_lcm / num_lcm;
2274 if (
gcd(num, den, &cnum, &cden,
false) !=
_ex1) {
2281 if (is_exactly_a<numeric>(den)) {
2290 if (ex_to<numeric>(den.
unit(
x)).is_negative()) {
2299 return dynallocate<lst>({num * pre_factor.
numer(), den * pre_factor.
denom()});
2310 nums.reserve(
seq.size()+1);
2311 dens.reserve(
seq.size()+1);
2312 int nmod = modifier.
nops();
2313 for (
auto & it :
seq) {
2315 nums.push_back(
n.op(0));
2316 dens.push_back(
n.op(1));
2319 nums.push_back(
n.op(0));
2320 dens.push_back(
n.op(1));
2328 auto num_it = nums.begin(), num_itend = nums.end();
2329 auto den_it = dens.begin(), den_itend = dens.end();
2331 for (
int imod = nmod; imod < modifier.
nops(); ++imod) {
2332 while (num_it != num_itend) {
2339 num_it = nums.begin();
2340 den_it = dens.begin();
2343 ex num = *num_it++, den = *den_it++;
2344 while (num_it != num_itend) {
2346 ex next_num = *num_it++, next_den = *den_it++;
2349 while ((den_it != den_itend) && next_den.is_equal(*den_it)) {
2350 next_num += *num_it;
2356 ex co_den1, co_den2;
2357 ex g =
gcd(den, next_den, &co_den1, &co_den2,
false);
2358 num = ((num * co_den2) + (next_num * co_den1)).
expand();
2377 int nmod = modifier.
nops();
2378 for (
auto & it :
seq) {
2380 num.push_back(
n.op(0));
2381 den.push_back(
n.op(1));
2383 n = ex_to<numeric>(
overall_coeff).normal(repl, rev_lookup, modifier);
2384 num.push_back(
n.op(0));
2385 den.push_back(
n.op(1));
2386 auto num_it = num.begin(), num_itend = num.end();
2387 auto den_it = den.begin(), den_itend = den.end();
2388 for (
int imod = nmod; imod < modifier.
nops(); ++imod) {
2389 while (num_it != num_itend) {
2395 num_it = num.begin();
2396 den_it = den.begin();
2400 return frac_cancel(dynallocate<mul>(num), dynallocate<mul>(den));
2411 int nmod = modifier.
nops();
2412 ex n_basis = ex_to<basic>(
basis).
normal(repl, rev_lookup, modifier);
2413 for (
int imod = nmod; imod < modifier.
nops(); ++imod)
2416 nmod = modifier.
nops();
2418 for (
int imod = nmod; imod < modifier.
nops(); ++imod)
2420 n_exponent = n_exponent.
op(0) / n_exponent.
op(1);
2427 return dynallocate<lst>({
pow(n_basis.
op(0), n_exponent),
pow(n_basis.
op(1), n_exponent)});
2432 return dynallocate<lst>({
pow(n_basis.
op(1), -n_exponent),
pow(n_basis.
op(0), -n_exponent)});
2468 for (
auto & it :
seq) {
2491 exmap repl, rev_lookup;
2494 ex e =
bp->normal(repl, rev_lookup, modifier);
2498 if (!repl.empty()) {
2499 for(
int i=0; i < modifier.
nops(); ++i)
2505 return e.
op(0) / e.
op(1);
2516 exmap repl, rev_lookup;
2519 ex e =
bp->normal(repl, rev_lookup, modifier);
2526 for(
int i=0; i < modifier.
nops(); ++i)
2541 exmap repl, rev_lookup;
2544 ex e =
bp->normal(repl, rev_lookup, modifier);
2551 for(
int i=0; i < modifier.
nops(); ++i)
2566 exmap repl, rev_lookup;
2569 ex e =
bp->normal(repl, rev_lookup, modifier);
2576 for(
int i=0; i < modifier.
nops(); ++i)
2599 return bp->to_rational(repl);
2604 return bp->to_polynomial(repl);
2691 if (is_exactly_a<mul>(basis_pref) || is_exactly_a<power>(basis_pref)) {
2708 s.reserve(
seq.size());
2709 for (
auto & it :
seq)
2724 s.reserve(
seq.size());
2725 for (
auto & it :
seq)
2742 if (is_exactly_a<add>(e)) {
2744 size_t num = e.
nops();
2745 exvector terms; terms.reserve(num);
2749 for (
size_t i=0; i<num; i++) {
2752 if (is_exactly_a<add>(
x) || is_exactly_a<mul>(
x) || is_a<power>(
x)) {
2776 for (
size_t i=0; i<num; i++) {
2781 if (is_exactly_a<mul>(t)) {
2782 for (
size_t j=0; j<t.
nops(); j++) {
2785 for (
size_t k=0;
k<t.
nops();
k++) {
2789 v.push_back(t.
op(
k));
2791 t = dynallocate<mul>(v);
2801 return dynallocate<add>(terms);
2803 }
else if (is_exactly_a<mul>(e)) {
2805 size_t num = e.
nops();
2808 for (
size_t i=0; i<num; i++)
2811 return dynallocate<mul>(v);
2813 }
else if (is_exactly_a<power>(e)) {
2814 const ex e_exp(e.
op(1));
2820 return pow(pre_res, e_exp);
2834 if (is_exactly_a<add>(e) || is_exactly_a<mul>(e) || is_exactly_a<power>(e)) {
2854 throw(std::runtime_error(
"resultant(): arguments must be polynomials"));
2856 const int h1 = ee1.
degree(s);
2857 const int l1 = ee1.
ldegree(s);
2858 const int h2 = ee2.
degree(s);
2859 const int l2 = ee2.
ldegree(s);
2861 const int msize = h1 + h2;
2864 for (
int l = h1; l >= l1; --l) {
2866 for (
int k = 0;
k < h2; ++
k)
2869 for (
int l = h2; l >= l2; --l) {
2871 for (
int k = 0;
k < h1; ++
k)
2872 m(
k+h2,
k+h2-l) = e;
2875 return m.determinant();
Interface to GiNaC's initially known functions.
container & append(const ex &b)
Add element at back.
std::vector< expair > epvector
expair-vector
ex expand(unsigned options=0) const
Interface to class signaling failure of operation.
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.
Interface to GiNaC's constant types and some special constants.
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.
bool is_integer() const
True if object is a non-complex integer.
numeric max_coefficient() const override
Implementation ex::max_coefficient().
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...
Interface to GiNaC's products of expressions.
ex to_rational(exmap &repl) const
Rationalization of non-rational functions.
ex normal() const
Normalization of rational functions.
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].
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for pseries.
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for symbols.
pseries(const ex &rel_, const epvector &ops_)
Construct pseries from a vector of coefficients and powers.
ex coeff(const ex &s, int n=1) const
const numeric iquo(const numeric &a, const numeric &b)
Numeric integer quotient.
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a product.
This class holds a relation consisting of two expressions and a logical relation between them.
Makes the interface to the underlying bignum package available.
Function object for map().
numeric max_coefficient() const override
Implementation ex::max_coefficient().
static ex frac_cancel(const ex &n, const ex &d)
Fraction cancellation.
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for a numeric.
ex subs(const exmap &m, unsigned options=0) const
bool is_rational() const
True if object is an exact rational number, may even be complex (denominator may be unity).
Interface to GiNaC's sums of expressions.
function psi(const T1 &p1)
const numeric imag() const
Imaginary part of a number.
ex denom() const
Get denominator of an expression.
sym_desc(const ex &s)
Initialize symbol, leave other variables uninitialized.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
std::vector< ex > exvector
This class holds a two-component object, a basis and and exponent representing exponentiation.
ex sqrfree(const ex &a, const lst &l)
Compute a square-free factorization of a multivariate polynomial in Q[X].
@ evaluated
.eval() has already done its job
size_t max_lcnops
Maximum number of terms of leading coefficient of symbol in both polynomials.
bool is_rational(const numeric &x)
bool is_real() const
True if object is a real integer, rational or float (but not complex).
int degree(const ex &s) const
virtual expair split_ex_to_pair(const ex &e) const
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.
virtual ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const
Default implementation of ex::normal().
ex resultant(const ex &e1, const ex &e2, const ex &s)
Resultant of two expressions e1,e2 with respect to symbol s.
int deg_b
Highest degree of symbol in polynomial "b".
static void add_symbol(const ex &s, sym_desc_vec &v)
virtual size_t nops() const
Number of operands/members.
std::vector< sym_desc > sym_desc_vec
Interface to relations between expressions.
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
bool is_equal(const ex &other) const
ex operator()(const ex &e) override
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 to_rational(exmap &repl) const override
Implementation of ex::to_rational() for symbols.
ex unit(const ex &x) const
Compute unit part (= sign of leading coefficient) of a multivariate polynomial in Q[x].
ex primpart(const ex &x) const
Compute primitive part of a multivariate polynomial in Q[x].
This structure holds information about the highest and lowest degrees in which a symbol appears in tw...
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for powers.
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for powers.
bool operator<(const sym_desc &x) const
Comparison operator for sorting.
@ no_heur_gcd
Usually GiNaC tries heuristic GCD first, because typically it's much faster than anything else.
ex numer_denom() const
Get numerator and denominator of an expression.
const basic & clearflag(unsigned f) const
Clear some status_flags.
Interface to symbolic matrices.
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...
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a sum.
@ no_part_factored
GiNaC tries to avoid expanding expressions when computing GCDs.
static ex gcd_pf_pow_pow(const ex &a, const ex &b, ex *ca, ex *cb)
ex smod(const numeric &xi) const
bool info(unsigned inf) const
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...
virtual ex default_overall_coeff() const
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a numeric.
static void collect_symbols(const ex &e, sym_desc_vec &v)
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].
Lightweight wrapper for GiNaC's symbolic objects.
virtual ex smod(const numeric &xi) const
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition of GiNaC's lst.
const_iterator begin() const noexcept
bool is_zero() const
True if object is zero.
const numeric isqrt(const numeric &x)
Integer numeric square root.
ex normal(const ex &thisex)
numeric integer_content() const override
numeric max_coefficient() const
Return maximum (absolute value) coefficient of a polynomial.
static epvector sqrfree_yun(const ex &a, const symbol &x)
Compute square-free factorization of multivariate polynomial a(x) using Yun's algorithm.
ex numer(const ex &thisex)
bool find(const ex &pattern, exset &found) const
Find all occurrences of a pattern.
const numeric exp(const numeric &x)
Exponential function.
@ hash_calculated
.calchash() has already done its job
ex diff(const symbol &s, unsigned nth=1) const
Compute partial derivative of an expression.
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].
int ldegree(const ex &s) const
ex to_polynomial(exmap &repl) const
const numeric inverse() const
Inverse of a number.
ex expand(const ex &thisex, unsigned options=0)
ex coeff(const ex &s, int n=1) const override
Return coefficient of degree n in object s.
Interface to GiNaC's light-weight expression handles.
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for powers.
const numeric I
Imaginary unit.
virtual ex coeff(const ex &s, int n=1) const
Return coefficient of degree n in object s.
Interface to GiNaC's symbolic objects.
ex op(size_t i) const override
Return operand/member at position i.
static numeric lcm_of_coefficients_denominators(const ex &e)
Compute LCM of denominators of coefficients of a polynomial.
container & remove_first()
Remove first element.
static ex interpolate(const ex &gamma, const numeric &xi, const ex &x, int degree_hint=1)
xi-adic polynomial interpolation
ex subs(const exmap &m, unsigned options=0) const override
Substitute a set of objects by arbitrary expressions.
ex recombine_pair_to_ex(const expair &p) const override
ex collect_common_factors(const ex &e)
Collect common factors in sums.
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 sqrfree_parfrac(const ex &a, const symbol &x)
Compute square-free partial fraction decomposition of rational function a(x).
ex denom(const ex &thisex)
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
int int_length() const
Size in binary notation.
virtual numeric integer_content() const
epvector seq
Vector of {coefficient, power} pairs.
virtual ex thisexpairseq(const epvector &v, const ex &oc, bool do_index_renaming=false) const
Wrapper template for making GiNaC classes out of STL containers.
static bool get_first_symbol(const ex &e, ex &x)
Return pointer to first symbol found in expression.
@ no_pattern
disable pattern matching
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].
virtual ex to_polynomial(exmap &repl) const
std::map< ex, ex, ex_is_less > exmap
Exception thrown by heur_gcd() to signal failure.
ex lcoeff(const ex &s) const
int ldeg_b
Lowest degree of symbol in polynomial "b".
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
ex coeff(const ex &thisex, const ex &s, int n=1)
ex factor(const ex &poly, unsigned options)
Interface function to the outside world.
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for symbols.
ex var
Series variable (holds a symbol)
int degree(const ex &thisex, const ex &s)
ex numer_denom(const ex &thisex)
Interface to GiNaC's ABC.
Interface to sequences of expression pairs.
Function object to be applied by basic::normal().
const numeric pow(const numeric &x, const numeric &y)
Interface to class for extended truncated power series.
int deg_a
Highest degree of symbol in polynomial "a".
numeric integer_content() const override
static ex gcd_pf_pow(const ex &a, const ex &b, ex *ca, ex *cb)
numeric max_coefficient() const override
Implementation ex::max_coefficient().
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for a numeric.
bool is_negative(const numeric &x)
#define is_ex_the_function(OBJ, FUNCNAME)
static ex gcd_pf_mul(const ex &a, const ex &b, ex *ca, ex *cb)
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for expairseqs.
static numeric lcmcoeff(const ex &e, const numeric &l)
ex recombine_pair_to_ex(const expair &p) const override
ex lcm(const ex &a, const ex &b, bool check_args)
Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X].
int ldeg_a
Lowest degree of symbol in polynomial "a".
ex expand(unsigned options=0) const override
Expand expression, i.e.
ex numer() const
Get numerator of an expression.
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,...
int max_deg
Maximum of deg_a and deg_b (Used for sorting)
virtual numeric max_coefficient() const
Implementation ex::max_coefficient().
const numeric abs(const numeric &x)
Absolute value.
const numeric smod(const numeric &a_, const numeric &b_)
Modulus (in symmetric representation).
const numeric real() const
Real part of a number.
bool is_integer(const numeric &x)
size_t nops() const override
Number of operands/members.
const numeric denom() const
Denominator.
virtual ex map(map_function &f) const
Construct new expression by applying the specified function to all sub-expressions (one level only,...
const numeric numer() const
Numerator.
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).
static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
Collect statistical information about symbols in polynomials.
ptr< basic > bp
pointer to basic object managed by this
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].
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].
Interface to GiNaC's overloaded operators.
numeric integer_content() const
Compute the integer content (= GCD of all numeric coefficients) of an expanded polynomial.
numeric integer_content() const override
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for expairseqs.
ex content(const ex &x) const
Compute content part (= unit normal GCD of all coefficients) of a multivariate polynomial in Q[x].
virtual ex recombine_pair_to_ex(const expair &p) const
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
@ use_sr_gcd
By default GiNaC uses modular GCD algorithm.
This file defines several functions that work on univariate and multivariate polynomials and rational...
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
virtual ex to_rational(exmap &repl) const
Default implementation of ex::to_rational().
ex sym
Reference to symbol.
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].
This page is part of the GiNaC
developer's reference. It was generated automatically by doxygen. For
an introduction, see the tutorial.