Skip to content

Instantly share code, notes, and snippets.

@marionette-of-u

marionette-of-u/symbolic_integration.cpp

Last active August 29, 2015 14:26
Show Gist options
  • Save marionette-of-u/ef42366f82346fdf6496 to your computer and use it in GitHub Desktop.
Save marionette-of-u/ef42366f82346fdf6496 to your computer and use it in GitHub Desktop.
Download ZIP
Symbolic Integration vol.1
Raw
symbolic_integration.cpp
#include <iostream>
#include <vector>
#include <functional>
#include <utility>
#include <memory>
#include <stdexcept>
#include <iterator>
#include <algorithm>
#include <limits>
#include <cassert>
#include <ginac/ginac.h>
namespace qz{
using namespace GiNaC;
inline int factor_class(const ex &f){
exmap map;
if(f.match(pow(wild(0), wild(1)), map)){
const ex &g = map[wild(0)];
if(is_a<symbol>(g)){
return 7;
}
numeric gp = ex_to<numeric>(g);
if(imag(gp) != 0){
if(real(gp) != 0){
return 6;
}else{
return 5;
}
}
return 4;
}else{
if(is_a<symbol>(f)){
return 3;
}
numeric fp = ex_to<numeric>(f);
if(imag(fp) != 0){
if(real(fp) != 0){
return 2;
}else{
return 1;
}
}
return 0;
}
};
inline int lexicographical_compare(const ex &l, const ex &r){
if(is_a<add>(l) && is_a<add>(r) || is_a<mul>(l) && is_a<mul>(r)){
auto iter = l.begin(), jter = r.begin();
int result;
for(; ; ){
result = lexicographical_compare(*iter, *jter);
if(result != 0){
return result;
}
++iter, ++jter;
if(iter == l.end() && jter != r.end()){ return -1; }
if(iter != l.end() && jter == r.end()){ return +1; }
if(iter == l.end() && jter == r.end()){ return 0; }
}
return result;
}
if(!is_a<mul>(l) && is_a<mul>(r) || is_a<add>(l) && !is_a<add>(r)){
auto iter = l.begin(), jter = r.begin();
std::size_t n = 0;
int result;
for(; iter != l.end(); ++iter, ++n){
for(; jter != r.end(); ++jter){
result = lexicographical_compare(*iter, *jter);
if(result != 0){ return result; }
}
}
return n == 1 ? -1 : +1;
}
if(is_a<mul>(l) && !is_a<mul>(r) || !is_a<add>(l) && is_a<add>(r)){
auto iter = l.begin(), jter = r.begin();
std::size_t n = 0;
int result;
for(; jter != r.end(); ++jter, ++n){
for(; iter != l.end(); ++iter){
result = lexicographical_compare(*iter, *jter);
if(result != 0){ return result; }
}
}
return n == 1 ? +1 : -1;
}
auto primitive_compare = [](auto &&a, auto &&b){
return a < b ? -1 : a > b ? +1 : 0;
};
if(l == 0 && r == 0){ return 0; }
if(l == 0 && r != 0){ return -1; }
if(l != 0 && r ==0 ){ return +1; }
int l_class = factor_class(l), r_class = factor_class(r);
if(l_class < r_class){
return -1;
}else if(l_class > r_class){
return +1;
}
int result = 0;
switch(l_class){
case 7:
{
exmap lmap, rmap;
l.match(pow(wild(0), wild(1)), lmap);
symbol fp = ex_to<symbol>(lmap[wild(0)]);
r.match(pow(wild(0), wild(1)), rmap);
symbol gp = ex_to<symbol>(rmap[wild(0)]);
result = primitive_compare(fp.get_name(), gp.get_name());
if(result == 0){ result = lexicographical_compare(lmap[wild(1)], rmap[wild(1)]); }
}
break;
case 6:
{
exmap lmap, rmap;
l.match(pow(wild(0), wild(1)), lmap);
numeric fp = ex_to<numeric>(lmap[wild(0)]);
r.match(pow(wild(0), wild(1)), rmap);
numeric gp = ex_to<numeric>(rmap[wild(0)]);
result = primitive_compare(real(fp), real(gp));
if(result == 0){ result = primitive_compare(imag(fp), imag(gp)); }
if(result == 0){ result = lexicographical_compare(lmap[wild(1)], rmap[wild(1)]); }
}
break;
case 5:
{
exmap lmap, rmap;
l.match(pow(wild(0), wild(1)), lmap);
numeric fp = ex_to<numeric>(lmap[wild(0)]);
r.match(pow(wild(0), wild(1)), rmap);
numeric gp = ex_to<numeric>(rmap[wild(0)]);
result = primitive_compare(imag(fp), imag(gp));
if(result == 0){ result = lexicographical_compare(lmap[wild(1)], rmap[wild(1)]); }
}
break;
case 4:
{
exmap lmap, rmap;
l.match(pow(wild(0), wild(1)), lmap);
numeric fp = ex_to<numeric>(lmap[wild(0)]);
r.match(pow(wild(0), wild(1)), rmap);
numeric gp = ex_to<numeric>(rmap[wild(0)]);
result = primitive_compare(real(fp), real(gp));
if(result == 0){ result = lexicographical_compare(lmap[wild(1)], rmap[wild(1)]); }
}
break;
case 3:
{
symbol fp = ex_to<symbol>(l), gp = ex_to<symbol>(r);
result = primitive_compare(ex_to<symbol>(l).get_name(), ex_to<symbol>(r).get_name());
}
break;
case 2:
{
numeric f = ex_to<numeric>(l), g = ex_to<numeric>(r);
result = primitive_compare(real(f), real(g));
if(result == 0){ result = primitive_compare(imag(f), imag(g)); }
}
break;
case 1:
{
numeric f = ex_to<numeric>(l), g = ex_to<numeric>(r);
result = primitive_compare(imag(f), imag(g));
}
break;
case 0:
{
numeric f = ex_to<numeric>(l), g = ex_to<numeric>(r);
result = primitive_compare(real(f), real(g));
}
break;
}
return result;
}
inline const ex lt(const ex &a){
if(is_a<add>(a)){
return *a.begin();
}else{
return a;
}
}
inline std::tuple<ex, ex> poly_divide(const ex &A, const ex &B, const ex &x){
std::tuple<ex, ex> QR;
ex &Q = std::get<0>(QR), &R = std::get<1>(QR);
R = A.expand();
for(; ; ){
if(R == 0){ break; }
int delta = R.degree(x) - B.degree(x);
if(delta < 0){ break; }
ex T = (R.lcoeff(x) / B.lcoeff(x) * pow(x, delta)).expand();
Q = (Q + T).expand();
R = (R - B * T).expand();
}
return QR;
}
inline std::tuple<ex, ex> poly_pseudo_divide(const ex &A, const ex &B, const ex &x){
ex b = B.lcoeff(x);
int N = A.degree(x) - B.degree(x) + 1;
ex Q = 0, R = A.expand();
for(; ; ){
if(R == 0){ break; }
int delta = R.degree(x) - B.degree(x);
if(delta < 0){ break; }
ex T = R.lcoeff(x) * pow(x, delta);
--N;
Q = (b * Q + T).expand();
R = (b * R - T * B).expand();
}
return std::make_tuple((pow(b, N) * Q).expand(), (pow(b, N) * R).expand());
}
inline std::tuple<ex, ex> multivar_poly_divide(const ex &f, const ex &g){
auto check_exp = [](const ex &a){
if(is_a<mul>(a)){
for(auto iter = a.begin(); iter != a.end(); ++iter){
exmap map;
if(iter->match(pow(wild(0), wild(1)), map)){
if(map[wild(1)] < 0){
return false;
}
}
}
return true;
}else{
return true;
}
};
ex r = 0, p = f, q = 0;
ex ltf = lt(f);
while(p != 0){
ex ltp = lt(p), ltg = lt(g);
ex h = (ltp / ltg).expand();
if(check_exp(h)){
q = (q + h).expand();
p = (p - h * g).expand();
}else{
r = (r + ltp).expand();
p = (p - ltp).expand();
}
}
return std::make_tuple(q, r);
}
inline ex euclidean(
const ex &A,
const ex &B,
const ex &x,
std::function<std::tuple<ex, ex>(const ex&, const ex&, const ex&)> euclidean_division = poly_pseudo_divide
){
ex a = A, b = B;
while(b != 0){
auto qr = euclidean_division(a, b, x);
a = b;
b = std::get<1>(qr);
}
return a;
}
inline ex multivar_euclidean(const ex &A, const ex &B){
ex a = A, b = B;
while(b != 0){
auto qr = multivar_poly_divide(a, b);
a = b;
b = std::get<1>(qr);
}
return a;
}
inline std::tuple<ex, ex, ex> extended_euclidean(
const ex &A,
const ex &B,
const ex &x,
std::function<std::tuple<ex, ex>(const ex&, const ex&, const ex&)> euclidean_division = poly_pseudo_divide
){
ex a = A, b = B, a1 = 1, a2 = 0, b1 = 0, b2 = 1;
while(b != 0){
auto qr = euclidean_division(a, b, x);
a = b, b = std::get<1>(qr);
ex r1 = (a1 - std::get<0>(qr) * b1).expand(), r2 = (a2 - std::get<0>(qr) * b2).expand();
a1 = b1, a2 = b2, b1 = r1, b2 = r2;
}
return std::make_tuple(a1, a2, a);
}
inline std::tuple<ex, ex> half_extended_euclidean(
const ex &A,
const ex &B,
const ex &x,
std::function<std::tuple<ex, ex>(const ex&, const ex&, const ex&)> euclidean_division = poly_divide
){
ex a = A, b = B, a1 = 1, b1 = 0;
while(b != 0){
auto qr = euclidean_division(a, b, x);
a = b, b = std::get<1>(qr);
ex r1 = a1 - std::get<0>(qr) * b1;
a1 = b1, b1 = r1.expand();
}
return std::make_tuple(a1, a);
}
inline std::tuple<ex, ex, ex> extended_euclidean_half_full(
const ex &a,
const ex &b,
const ex &x,
std::function<std::tuple<ex, ex>(const ex&, const ex&, const ex&)> euclidean_division = poly_divide
){
auto sg = half_extended_euclidean(a, b, x, euclidean_division);
auto tr = euclidean_division(std::get<1>(sg) - std::get<0>(sg) * a, b, x);
return std::make_tuple(std::get<0>(sg), std::get<0>(tr), std::get<1>(sg));
}
inline std::tuple<ex, ex> extended_euclidean(
const ex &a,
const ex &b,
const ex &c,
const ex &x,
std::function<std::tuple<ex, ex>(const ex&, const ex&, const ex&)> euclidean_division = poly_divide
){
auto stg = extended_euclidean(a, b, x, euclidean_division);
ex &s = std::get<0>(stg), &t = std::get<1>(stg), &g = std::get<2>(stg);
auto qr = euclidean_division(c, g, x);
ex &q = std::get<0>(qr), &r = std::get<1>(qr);
if(r != 0){
throw std::runtime_error("c is not in the ideal generated by a and b.");
}
s = (q * s).expand(), t = (q * t).expand();
if(s != 0 && s.degree(x) >= b.degree(x)){
auto qr1 = euclidean_division(s, b, x);
q = std::get<0>(qr1), r = std::get<1>(qr1);
s = r, t = (t + q * a).expand();
}
return std::make_tuple(s, t);
}
inline ex half_extended_euclidean(const ex &a, const ex &b, const ex &c, const ex &x){
auto sg = half_extended_euclidean(a, b, x);
ex &s = std::get<0>(sg), &g = std::get<1>(sg);
auto qr = poly_divide(c, g, x);
ex &q = std::get<0>(qr), &r = std::get<1>(qr);
if(r != 0){
throw std::runtime_error("c is not in the ideal generated by a and b.");
}
s = (q * s).expand();
if(s != 0 && s.degree(x) >= b.degree(x)){
auto qr1 = poly_divide(s, b, x);
s = std::get<1>(qr1);
}
return s;
}
inline std::tuple<ex, ex> extended_euclidean_half_full(const ex &a, const ex &b, const ex &c, const ex &x){
ex s = half_extended_euclidean(a, b, c, x);
auto tr = poly_divide(c - s * a, b, x);
return std::make_tuple(s, std::get<0>(tr));
}
inline std::vector<ex> partial_fraction_impl(
const ex &a,
std::vector<ex>::const_iterator begin,
std::vector<ex>::const_iterator end,
const ex &x
){
ex d = 1;
for(auto iter = begin; iter != end; ++iter){
d = d * *iter;
}
d = d.expand();
auto a0_r = poly_divide(a, d, x);
ex &a0 = std::get<0>(a0_r), &r = std::get<1>(a0_r);
if(std::distance(begin, end) == 1){
std::vector<ex> v = { a0, r };
return v;
}
ex d2_dn = 1;
for(auto iter = begin + 1; iter != end; ++iter){
d2_dn = (d2_dn * *iter).expand();
}
auto a1_t = extended_euclidean_half_full(d2_dn, *begin, r, x);
ex &a1 = std::get<0>(a1_t), &t = std::get<1>(a1_t);
std::vector<ex> b0_a2_an = partial_fraction_impl(t, begin + 1, end, x);
std::vector<ex> result;
result.reserve(b0_a2_an.size() + 1);
result.push_back((a0 + b0_a2_an.front()).expand());
result.push_back(std::move(a1));
for(auto iter = b0_a2_an.begin() + 1; iter != b0_a2_an.end(); ++iter){
result.push_back(*iter);
}
return result;
}
inline std::vector<ex> partial_fraction(
const ex &a,
const std::vector<ex> &v,
const ex &x
){
return partial_fraction_impl(a, v.begin(), v.end(), x);
}
inline std::vector<ex> complete_partial_fraction(
const ex &a,
const std::vector<std::pair<ex, int>> &de,
const ex &x
){
std::vector<ex> pf_d;
pf_d.reserve(de.size());
for(auto &p : de){
pf_d.push_back(pow(p.first, p.second).expand());
}
std::vector<ex> an = partial_fraction(a, pf_d, x);
std::vector<std::vector<ex>> anm;
anm.resize(an.size() - 1);
for(std::size_t i = 1; i < an.size(); ++i){
for(std::size_t j = de[i - 1].second; ; ){
auto qa = poly_divide(an[i], de[i - 1].first, x);
ex &q = std::get<0>(qa), &a = std::get<1>(qa);
anm[i - 1].push_back(a);
an[i] = q;
if(j == 1){ break; }
--j;
}
an[0] = (an[0] + an[i]).expand();
}
std::vector<ex> result;
std::size_t s = 1;
for(auto &i : anm){ s += i.size(); }
result.reserve(s);
result.push_back(an[0]);
for(auto &i : anm){
for(auto &j : i){
result.push_back(j);
}
}
return result;
}
inline std::pair<ex, std::vector<ex>> subresultant(const ex &A, const ex &B, const ex &x){
std::vector<ex> R;
R.push_back(A);
R.push_back(B);
ex gamma = -1;
std::vector<ex> Delta = { A.degree(x) - B.degree(x) };
std::vector<ex> Beta = { (A.degree(x) - B.degree(x) + 1) & 1 == 1 ? -1 : +1 };
ex *beta = &Beta.back();
std::vector<ex> r_vec;
while(!R.back().is_zero()){
r_vec.push_back(R.back().lcoeff(x));
std::tuple<ex, ex> qr = poly_pseudo_divide(R[R.size() - 2], R[R.size() - 1], x);
ex &rem = std::get<1>(qr);
R.push_back(rem / *beta);
gamma = (pow(-r_vec.back(), Delta.back()) * pow(gamma, 1 - Delta.back())).expand();
Delta.push_back(R[R.size() - 2].degree(x) - R[R.size() - 1].degree(x));
Beta.push_back((-r_vec.back() * pow(gamma, Delta.back())).expand());
beta = &Beta.back();
}
std::size_t k = R.size() - 2;
if(R[k].degree(x) > 0){
return std::make_pair(0, std::move(R));
}
if(R[k - 1].degree(x) == 1){
return std::make_pair(R[k], std::move(R));
}
int s = 1;
ex c = 1;
for(std::size_t j = 1; j < k - 1; ++j){
if(R[j - 1].degree(x) & 1 == 1 && R[j].degree(x) & 1 == 1){
s = -s;
}
c = (c * pow((Beta[j - 1] / pow(r_vec[j - 1], 1 + Delta[j - 1])), R[j].degree(x)) * pow(r_vec[j - 1], R[j - 1].degree(x) - R[j + 1].degree(x)).expand()).expand();
}
return std::make_pair((s * c * pow(R[k], R[k - 1].degree(x))).expand(), std::move(R));
}
template<class Lambda>
inline std::vector<std::pair<int, ex>> decomp_factor(const Lambda &lambda){
std::vector<std::pair<int, ex>> result;
ex f = lambda();
if(is_a<mul>(f)){
for(auto iter = f.begin(); iter != f.end(); ++iter){
exmap map;
if(iter->match(pow(wild(0), wild(1)), map)){
result.push_back(std::make_pair(ex_to<numeric>(map[wild(1)]).to_int(), map[wild(0)]));
}else{
result.push_back(std::make_pair(1, *iter));
}
}
}else{
exmap map;
if(f.match(pow(wild(0), wild(1)), map)){
result.push_back(std::make_pair(ex_to<numeric>(map[wild(1)]).to_int(), map[wild(0)]));
}else{
result.push_back(std::make_pair(1, f));
}
}
return result;
}
inline std::vector<std::pair<int, ex>> factor_vec(const ex &x){
return decomp_factor([&]() -> ex{ return factor(x); });
}
inline ex content(const ex &A, const symbol &x){
if(is_a<add>(A)){
auto iter = A.begin(), jter = iter; ++jter;
ex f = gcd(iter->lcoeff(x), jter->lcoeff(x), nullptr, nullptr);
for(; jter != A.end(); ++jter){
f = gcd(f, jter->lcoeff(x));
}
return f;
}else{
return A.begin()->lcoeff(x);
}
}
inline std::vector<std::pair<int, ex>> squarefree(const ex &a, const symbol &x){
std::vector<std::pair<int, ex>> res;
ex w = a;
ex z = w.diff(x);
ex g = gcd(w, z);
if(g.is_equal(a)){
res.push_back(std::make_pair(1, a));
return res;
}
ex y;
int j = 1;
do{
w = quo(w, g, x);
y = quo(z, g, x);
z = y - w.diff(x);
g = gcd(w, z);
res.push_back(std::make_pair(j++, g));
}while(!z.is_zero());
return res;
}
inline std::tuple<ex, ex> hermite_reduce(const ex &A_, const ex &D, const symbol &x){
ex A = A_;
ex g = 0, D_minus = gcd(D, D.diff(x), nullptr, nullptr).expand(), D_star = quo(D, D_minus, x);
while(D_minus.degree(x) > 0){
ex D_minus2 = gcd(D_minus, D_minus.diff(x), nullptr, nullptr).expand();
ex D_minus_star = quo(D_minus, D_minus2, x);
auto BC = extended_euclidean(quo((-D_star * D_minus.diff(x)).expand(), D_minus, x), D_minus_star, A, x);
ex &B = std::get<0>(BC), &C = std::get<1>(BC);
A = (C - (quo((B.diff(x) * D_star).expand(), D_minus_star, x))).expand();
g = (g + (B / D_minus)).expand();
D_minus = D_minus2;
}
return std::make_tuple(g, A / D_star);
}
inline std::vector<numeric> horner(const std::vector<numeric> &a, std::size_t n, const numeric &x0){
std::vector<numeric> b;
b = a;
for(std::size_t i = 0; i <= n; ++i){
for(std::size_t j = 1; j <= n - i; ++j){
b[j] += x0 * b[j - 1];
}
}
return b;
}
inline std::vector<numeric> aberth(
const ex &e,
const std::vector<numeric> &coeff,
const numeric &max_coeff,
const symbol &x,
const numeric &epsilon = 1.0 / 1024.0
){
std::vector<numeric> z;
std::size_t n = e.degree(x);
z.resize(n);
numeric zc = -coeff[1] / (coeff[0] * n);
std::vector<numeric> a = horner(coeff, n, zc);
ex b, db;
b += abs(a[0]);
for(std::size_t i = 1; i <= n; ++i){
b += -abs(a[i]) * pow(x, i);
}
db = b.diff(x);
numeric xn = 100.0, xm;
for(; ; ){
xm = xn - ex_to<numeric>(b.subs(x == xn) / db.subs(x == xn));
if(abs(xn - xm) <= epsilon){
xn = xm;
break;
}
xn = xm;
}
numeric pi = ex_to<numeric>(Pi.evalf());
for(std::size_t j = 0; j < n; ++j){
numeric theta = (2 * pi / n) * j + pi / (2 * n);
z[j] = zc + cos(theta) + I * sin(theta);
}
return z;
}
inline std::vector<ex> nth_roots(const ex &poly, const symbol &x, const numeric epsilon = 1.0 / 1024.0){
numeric max_coeff = 0;
std::vector<numeric> coeff;
std::size_t n = poly.degree(x);
if(n == 1){
std::vector<ex> result({ (-poly.coeff(x, 0) / poly.coeff(x, 1)).expand() });
return result;
}
if(n == 2){
ex D = pow(pow(poly.coeff(x, 1), 2) - 4 * poly.coeff(x, 2) * poly.coeff(x, 0), ex(1) / ex(2));
std::vector<ex> result({
((-poly.coeff(x, 1) + D) / (2 * poly.coeff(x, 2))).expand(),
((-poly.coeff(x, 1) - D) / (2 * poly.coeff(x, 2))).expand()
});
return result;
}
if(n == 3){
ex A = poly.coeff(x, 3), B = poly.coeff(x, 2), C = poly.coeff(x, 1), D = poly.coeff(x, 0);
ex b = B / A, c = C/ A, d = D / A;
ex p = c - pow(b, 2) / 3, q = d - b * c / 3 + 2 * pow(b, 3) / 27;
symbol t;
ex T = pow(t, 2) + q * t - pow(p, 3) / 27;
std::vector<ex> sol_t = nth_roots(T.expand(), t);
ex u = pow(sol_t[0], ex(1) / ex(3)), v = pow(sol_t[1], ex(1) / ex(3)), w = exp(2 * Pi * I / 3);
std::vector<ex> result({
(u + v - b / 3).expand(),
(u * w + v * pow(w, 2) - b / 3).expand(),
(u * pow(w, 2) + v * w - b / 3).expand(),
});
return result;
}
if(n == 4){
auto cbrt = [](ex a, int n){
ex Pi2 = 2 * Pi;
ex rho = pow(abs(a), ex(1) / ex(3));
ex theta = ((Pi2 * n) + atan2(imag_part(a), real_part(a))) / 3;
return (rho * cos(theta) + rho * sin(theta) * I).expand();
};
ex a = poly.coeff(x, 4), b = poly.coeff(x, 3), c = poly.coeff(x, 2), d = poly.coeff(x, 1), e = poly.coeff(x, 0);
b /= a;
c /= a;
d /= a;
e /= a;
ex b2 = b * b;
ex b3 = b * b2;
ex b4 = b2 * b2;
ex alpha = -3 * b2 / 8 + c;
ex beta = (b3 / 8 - b * c / 2 + d).expand();
ex gamma = -3 *b4 / 256 + b2 * c / 16 - b * d / 4 + e;
ex alpha2 = alpha * alpha;
ex t = -b / 4;
if(beta.is_zero()){
ex rad = pow(alpha2, ex(1) / ex(2));
ex r1 = pow((-alpha + rad) / 2, ex(1) / ex(2));
ex r2 = pow((-alpha - rad) / 2, ex(1) / ex(2));
std::vector<ex> result({
(t + r1).expand(),
(t - r1).expand(),
(t + r2).expand(),
(t - r2).expand()
});
return result;
}else{
ex alpha3 = alpha * alpha2;
ex P = -(alpha2 / 12 + gamma);
ex Q = -alpha3 / 108 + alpha * gamma / 3 - beta * beta / 8;
ex R = -Q / 2 + pow(pow(Q, 2) / 4 + pow(P, 3) / 27, ex(1) / ex(2));
ex U = cbrt(R, 0);
ex y = -5 * alpha / 6 + U;
if(U.is_zero()){
y -= cbrt(Q, 0);
}else{
y -= P / (3 * U);
}
ex W = pow(alpha + 2 * y, ex(1) / ex(2));
ex r1 = pow(-(3 * alpha + 2 * y + 2 * beta / W), ex(1) / ex(2));
ex r2 = pow(-(3 * alpha + 2 * y - 2 * beta / W), ex(1) / ex(2));
std::vector<ex> result({
(t + (+W + r1) / 2).expand(),
(t + (+W - r1) / 2).expand(),
(t + (-W + r2) / 2).expand(),
(t + (-W - r2) / 2).expand()
});
return result;
}
}
coeff.resize(n + 1);
for(std::size_t i = 0; i <= n; ++i){
coeff[i] = ex_to<numeric>(poly.coeff(x, i));
if(max_coeff < abs(coeff[i])){
max_coeff = abs(coeff[i]);
}
}
std::vector<numeric> z = aberth(poly, coeff, max_coeff, x, epsilon), zp;
for(; ; ){
zp = z;
for(std::size_t j = 0; j < n; ++j){
numeric f = ex_to<numeric>(poly.subs(x == z[j]));
numeric df = coeff[0];
for(std::size_t i = 0; i < n; ++i){
if(i != j){
df *= zp[j] - zp[i];
}
}
z[j] = zp[j] - f / df;
}
numeric err = 0;
for(std::size_t j = 0; j < n; ++j){
numeric err_prime = abs(ex_to<numeric>(poly.subs(x == z[j])));
if(err < err_prime){
err = err_prime;
}
}
if(err < epsilon){
break;
}
}
std::vector<ex> ze(z.size());
for(auto &a : z){
ze.push_back(a);
}
return ze;
}
inline ex int_rational_log_part(const ex &A, const ex &D, const symbol &x){
std::vector<std::pair<ex, ex>> QS;
ex sum;
symbol t;
std::pair<ex, std::vector<ex>> R = subresultant(D, (A - t * D.diff(x)).expand(), x);
std::vector<std::pair<int, ex>> Q = squarefree(R.first, t);
for(const std::pair<int, ex> &Qi : Q){
int i = Qi.first;
const ex &q = Qi.second;
if(q.degree(t) <= 0){
continue;
}
if(i == D.degree(x)){
std::vector<ex> roots = nth_roots(q, t);
for(auto &a : roots){
sum += a * log(D.subs(t == a));
}
}else{
ex Si = [&]() -> const ex&{
for(std::size_t m = 1; m < R.second.size() - 1; ++m){
if(R.second[m].degree(x) == i){
return R.second[m];
}
}
throw;
}();
std::vector<std::pair<int, ex>> Aq = squarefree(Si.lcoeff(x), t);
for(const std::pair<int, ex> &Aj : Aq){
int j = Aj.first;
const ex &a = Aj.second;
Si = Si / pow(gcd(a, q, nullptr, nullptr), j);
}
std::vector<ex> roots = nth_roots(q, t);
for(auto &a : roots){
sum += a * log(Si.subs(t == a));
}
}
++i;
}
return sum;
}
class p_integral_s{
public:
ex expression;
symbol x;
inline p_integral_s(){}
inline p_integral_s(const ex &expression, const symbol &x) : expression(expression), x(x){}
inline p_integral_s(const p_integral_s &other) : expression(other.expression), x(other.x){}
~p_integral_s() = default;
};
inline bool operator ==(const p_integral_s &lhs, const p_integral_s &rhs){
return lhs.expression == rhs.expression && lhs.x == rhs.x;
}
inline bool operator <(const p_integral_s &lhs, const p_integral_s &rhs){
int comp = lhs.x.compare(rhs.x);
if(comp != 0){
return comp;
}else{
return lhs.expression.compare(rhs.expression);
}
}
}
namespace GiNaC{
using p_integral = structure<qz::p_integral_s, compare_std_less>;
template<>
inline void p_integral::print(const print_context &c, unsigned level) const{
if(is_a<print_tree>(c)){
inherited::print(c, level);
}
const qz::p_integral_s &p_i = get_struct();
c.s << "∫(" << p_i.expression << ")d" << p_i.x;
}
template<>
inline std::size_t p_integral::nops() const{
return 1;
}
template<>
inline ex p_integral::op(std::size_t i) const{
switch(i){
case 0:
return get_struct().expression;
default:
throw std::range_error("p_integral::op(): no such operand");
}
}
template<>
inline ex &p_integral::let_op(std::size_t i){
ensure_if_modifiable();
switch(i){
case 0:
return get_struct().expression;
default:
throw std::range_error("p_integral::let_op(): no such operand");
}
}
template<>
inline const char *p_integral::get_class_name(){ return "p_integral"; }
}
namespace qz{
inline ex p_integ(const ex &expression, const symbol &x){
if(expression.is_zero()){
return 0;
}
return p_integral(p_integral_s(expression, x));
}
inline ex integrate_rational_function(const ex &f, const symbol &x){
std::tuple<ex, ex> gh = hermite_reduce(f.numer(), f.denom(), x);
std::tuple<ex, ex> QR = poly_divide(std::get<1>(gh).numer(), std::get<1>(gh).denom(), x);
if(std::get<1>(QR).is_zero()){
return std::get<0>(gh) + p_integ(std::get<0>(QR), x);
}else{
return std::get<0>(gh) + p_integ(std::get<0>(QR), x) + int_rational_log_part(std::get<1>(QR), std::get<1>(gh).denom(), x);
}
}
}
int main(){
using namespace GiNaC;
using namespace qz;
symbol t("t"), u("u"), v("v"), x("x"), y("y"), z("z");
ex f = ex(36) / (pow(x, 5) - 2 * pow(x, 4) - 2 * pow(x, 3) + 4 * pow(x, 2) + x - 2);
std::cout << integrate_rational_function(f, x).subs(x == Pi).evalf() << std::endl;
return 0;
}
void lexicographical_compare_test(){
using namespace GiNaC;
using namespace qz;
symbol t("t"), u("u"), v("v"), x("x"), y("y"), z("z");
{
ex e = t;
ex f = t;
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == 0);
}
{
ex e = t;
ex f = u;
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = u;
ex f = t;
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e = 1;
ex f = 1;
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == 0);
}
{
ex e = 2;
ex f = 1;
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e = 1;
ex f = 2;
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = 1 + I;
ex f = 1 + I;
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == 0);
}
{
ex e = 2 + I;
ex f = 1 + I;
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e = 1 + I;
ex f = 2 + I;
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = 1 + 2 * I;
ex f = 1 + I;
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e = 1 + I;
ex f = 1 + 2 * I;
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(x, 2);
ex f = pow(x, 2);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == 0);
}
{
ex e = pow(x, 2);
ex f = pow(y, 2);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(y, 2);
ex f = pow(x, 2);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e = pow(x, 3);
ex f = pow(x, 2);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e = pow(x, 2);
ex f = pow(x, 3);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(x, 3);
ex f = pow(y, 2);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(x, 2);
ex f = pow(y, 3);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(y, 3);
ex f = pow(x, 2);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e = pow(y, 2);
ex f = pow(x, 3);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e = pow(x, pow(x, 2));
ex f = pow(x, pow(x, 2));
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == 0);
}
{
ex e = pow(x, pow(x, 2));
ex f = pow(x, pow(y, 2));
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(x, pow(y, 2));
ex f = pow(x, pow(x, 2));
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e =pow(x, pow(x, 3));
ex f = pow(x, pow(x, 2));
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e = pow(x, pow(x, 2));
ex f = pow(x, pow(x, 3));
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(x, pow(x, 3));
ex f = pow(x, pow(y, 2));
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(x, pow(x, 2));
ex f = pow(x, pow(y, 3));
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(x, pow(y, 3));
ex f = pow(x, pow(x, 2));
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e = pow(x, pow(y, 2));
ex f = pow(x, pow(x, 3));
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e = pow(x, 1 + I);
ex f = pow(x, 1 + I);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == 0);
}
{
ex e = pow(x, 1 + I);
ex f = pow(y, 1 + I);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(x, 1 + 2 * I);
ex f = pow(x, 1 + I);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == +1);
}
{
ex e = pow(x, 1 + I);
ex f = pow(x, 1 + 2 * I);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(x, 1 + I) + pow(y, 1 + I);
ex f = pow(x, 1 + I) + pow(y, 1 + 2 * I);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(y, 1 + I) * pow(x, 1 + I);
ex f = pow(y, 1 + I) * pow(x, 1 + 2 * I);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(y, 1 + I) + pow(x, 1 + I);
ex f = pow(y, 1 + I) * pow(x, 1 + 2 * I);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(x, 1 + I);
ex f = pow(y, 1 + I) * pow(x, 1 + 2 * I);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(x, 1 + I);
ex f = pow(y, 1 + I) + pow(x, 1 + 2 * I);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
{
ex e = pow(x, 3 + I);
ex f = pow(y, 1 + I) + pow(x, 2);
std::cout << e << ", " << f << " : " << lexicographical_compare(e, f) << std::endl;
assert(lexicographical_compare(e, f) == -1);
}
}
// g++ -std=c++14 main.cpp -o test -O3 -lcln -lginac
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment