83cln::cl_N arithmetic_geometric_mean(
const cln::cl_N & a_0,
const cln::cl_N & b_0)
85 cln::cl_N a_old = a_0 * cln::cl_float(1, cln::float_format(
Digits));
86 cln::cl_N b_old = b_0 * cln::cl_float(1, cln::float_format(
Digits));
89 cln::cl_N res = a_old;
94 a_new = (a_old+b_old)/2;
95 b_new =
sqrt(a_old*b_old);
97 if ( (
abs(a_new-b_new) >
abs(a_new+b_new) )
99 ( (
abs(a_new-b_new) ==
abs(a_new+b_new)) && (imagpart(b_new/a_new) <= 0) ) ) {
107 }
while (res != resbuf);
118cln::cl_N agm_helper_second_kind(
const cln::cl_N & a_0,
const cln::cl_N & b_0,
const cln::cl_N & c_0)
120 cln::cl_N a_old = a_0 * cln::cl_float(1, cln::float_format(
Digits));
121 cln::cl_N b_old = b_0 * cln::cl_float(1, cln::float_format(
Digits));
122 cln::cl_N c_old = c_0 * cln::cl_float(1, cln::float_format(
Digits));
126 cln::cl_N res = square(a_old)-square(c_old)/2;
128 cln::cl_N pre = cln::cl_float(1, cln::float_format(
Digits));
132 a_new = (a_old+b_old)/2;
133 b_new =
sqrt(a_old*b_old);
135 if ( (
abs(a_new-b_new) >
abs(a_new+b_new) )
137 ( (
abs(a_new-b_new) ==
abs(a_new+b_new)) && (imagpart(b_new/a_new) <= 0) ) ) {
141 c_new = square(c_old)/4/a_new;
143 res -= pre*square(c_new);
150 }
while (res != resbuf);
169 return EllipticK(k).hold();
174 ex result =
Pi/2/
numeric(arithmetic_geometric_mean(1,kbar));
176 return result.
evalf();
187 return EllipticK(k).evalf();
190 return EllipticK(k).hold();
196 return -EllipticK(k)/k + EllipticE(k)/k/(1-k*k);
208 for (
int i=0; i<(order+1)/2; ++i)
210 ser +=
Pi/2 *
numeric(cln::square(cln::binomial(2*i,i))) *
pow(s/4,2*i);
216 ser +=
pseries(rel, std::move(nseq));
218 return ser.
series(rel, order);
221 if ( (k_pt ==
_ex1) || (k_pt ==
_ex_1) ) {
222 throw std::runtime_error(
"EllipticK_series: don't know how to do the series expansion at this point!");
231 c.
s <<
"\\mathrm{K}(";
243 do_not_evalf_params());
249 return EllipticE(k).hold();
254 ex result =
Pi/2 *
numeric( agm_helper_second_kind(1,kbar,
ex_to<numeric>(k).to_cl_N()) / arithmetic_geometric_mean(1,kbar) );
256 return result.
evalf();
271 return EllipticE(k).evalf();
274 return EllipticE(k).hold();
280 return -EllipticK(k)/k + EllipticE(k)/k;
292 for (
int i=0; i<(order+1)/2; ++i)
294 ser -=
Pi/2 *
numeric(cln::square(cln::binomial(2*i,i)))/(2*i-1) *
pow(s/4,2*i);
300 ser +=
pseries(rel, std::move(nseq));
302 return ser.
series(rel, order);
305 if ( (k_pt ==
_ex1) || (k_pt ==
_ex_1) ) {
306 throw std::runtime_error(
"EllipticE_series: don't know how to do the series expansion at this point!");
315 c.
s <<
"\\mathrm{K}(";
327 do_not_evalf_params());
342cln::cl_N iterated_integral_do_sum(
const std::vector<int> & m,
const std::vector<const integration_kernel *> & kernel,
const cln::cl_N & lambda,
int N_trunc)
344 if ( cln::zerop(lambda) ) {
348 cln::cl_F one = cln::cl_float(1, cln::float_format(
Digits));
350 const int depth = m.size();
356 if ( N_trunc == 0 ) {
358 bool flag_accidental_zero =
false;
367 multi_iterator_ordered_eq<int> i_multi(1,N+1,depth-1);
368 for( i_multi.init(); !i_multi.overflow(); i_multi++) {
370 for (
int j=1; j<depth; j++) {
372 tt = tt * kernel[0]->series_coeff(N-i_multi[depth-2]) / cln::expt(cln::cl_I(i_multi[depth-2]),m[1]);
375 tt = tt * kernel[j-1]->series_coeff(i_multi[depth-j]-i_multi[depth-j-1]) / cln::expt(cln::cl_I(i_multi[depth-j-1]),m[j]);
378 tt = tt * kernel[depth-1]->series_coeff(i_multi[0]);
384 subexpr = kernel[0]->series_coeff(N) * one;
386 flag_accidental_zero = cln::zerop(subexpr);
387 res += cln::expt(lambda, N) / cln::expt(cln::cl_I(N),m[0]) * subexpr;
390 }
while ( (res != resbuf) || flag_accidental_zero );
394 for (
int N=1; N<=N_trunc; N++) {
398 for( i_multi.init(); !i_multi.overflow(); i_multi++) {
400 for (
int j=1; j<depth; j++) {
402 tt = tt * kernel[0]->series_coeff(N-i_multi[depth-2]) / cln::expt(cln::cl_I(i_multi[depth-2]),m[1]);
405 tt = tt * kernel[j-1]->series_coeff(i_multi[depth-j]-i_multi[depth-j-1]) / cln::expt(cln::cl_I(i_multi[depth-j-1]),m[j]);
408 tt = tt * kernel[depth-1]->series_coeff(i_multi[0]);
414 subexpr = kernel[0]->series_coeff(N) * one;
416 res += cln::expt(lambda, N) / cln::expt(cln::cl_I(N),m[0]) * subexpr;
424cln::cl_N iterated_integral_prepare_m_lst(
const std::vector<const integration_kernel *> & kernel_in,
const cln::cl_N & lambda,
int N_trunc)
426 size_t depth = kernel_in.size();
431 std::vector<const integration_kernel *> kernel;
432 kernel.reserve(depth);
436 for (
const auto & it : kernel_in) {
447 cln::cl_N result = iterated_integral_do_sum(m, kernel, lambda, N_trunc);
454cln::cl_N iterated_integral_shuffle(
const std::vector<const integration_kernel *> & kernel,
const cln::cl_N & lambda,
int N_trunc)
456 cln::cl_F one = cln::cl_float(1, cln::float_format(
Digits));
458 const size_t depth = kernel.size();
460 size_t i_trailing = 0;
461 for (
size_t i=0; i<depth; i++) {
467 if ( i_trailing == 0 ) {
468 return cln::expt(cln::log(lambda), depth) / cln::factorial(depth) * one;
471 if ( i_trailing == depth ) {
472 return iterated_integral_prepare_m_lst(kernel, lambda, N_trunc);
476 std::vector<const integration_kernel *> a,b;
477 for (
size_t i=0; i<i_trailing; i++) {
478 a.push_back(kernel[i]);
480 for (
size_t i=i_trailing; i<depth; i++) {
481 b.push_back(kernel[i]);
484 cln::cl_N result = iterated_integral_prepare_m_lst(a, lambda, N_trunc) * cln::expt(cln::log(lambda), depth-i_trailing) / cln::factorial(depth-i_trailing);
486 for( i_shuffle.init(); !i_shuffle.overflow(); i_shuffle++) {
487 std::vector<const integration_kernel *> new_kernel;
488 new_kernel.reserve(depth);
489 for (
size_t i=0; i<depth; i++) {
490 new_kernel.push_back(i_shuffle[i]);
493 result -= iterated_integral_shuffle(new_kernel, lambda, N_trunc);
518 bool flag_not_numeric =
false;
519 for (
const auto & it : k_lst) {
521 flag_not_numeric =
true;
524 if ( flag_not_numeric) {
528 for (
const auto & it : k_lst) {
530 flag_not_numeric =
true;
533 if ( flag_not_numeric) {
542 const size_t depth = k_lst.
nops();
544 std::vector<const integration_kernel *> kernel_vec;
545 kernel_vec.reserve(depth);
547 for (
const auto & it : k_lst) {
551 size_t i_trailing = 0;
552 for (
size_t i=0; i<depth; i++) {
553 if ( !(kernel_vec[i]->has_trailing_zero()) ) {
560 cln::cl_F one = cln::cl_float(1, cln::float_format(
Digits));
569 result = iterated_integral_shuffle(kernel_vec, lambda_cln, N_trunc_int);
572 cln::cl_N
coeff = one;
573 for (
size_t i_plus=depth; i_plus>i_trailing; i_plus--) {
574 coeff =
coeff * kernel_vec[i_plus-1]->series_coeff(0);
575 kernel_vec[i_plus-1] = &L0;
576 if ( i_plus==i_trailing+1 ) {
577 result +=
coeff * iterated_integral_shuffle(kernel_vec, lambda_cln, N_trunc_int);
581 result +=
coeff * iterated_integral_shuffle(kernel_vec, lambda_cln, N_trunc_int);
621 do_not_evalf_params().
628 do_not_evalf_params().
Interface to GiNaC's sums of expressions.
The basic integration kernel with a logarithmic singularity at the origin.
const basic & hold() const
Stop further evaluation.
size_t nops() const override
Number of operands/members.
Exception class thrown by classes which provide their own series expansion to signal that ordinary Ta...
Lightweight wrapper for GiNaC's symbolic objects.
ex series(const ex &r, int order, unsigned options=0) const
Compute the truncated series expansion of an expression.
ex subs(const exmap &m, unsigned options=0) const
bool info(unsigned inf) const
void print(const print_context &c, unsigned level=0) const
Print expression to stream.
ex evalf() const override
Evaluate object numerically.
static unsigned register_new(function_options const &opt)
The class multi_iterator_ordered_eq defines a multi_iterator , such that.
The class multi_iterator_shuffle_prime defines a multi_iterator, which runs over all shuffles of a an...
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Base class for print_contexts.
std::ostream & s
stream to output to
This class holds a extended truncated power series (positive and negative integer powers).
This class holds a relation consisting of two expressions and a logical relation between them.
@ no_pattern
disable pattern matching
Interface to GiNaC's constant types and some special constants.
#define REGISTER_FUNCTION(NAME, OPT)
Interface to GiNaC's initially known functions.
Interface to GiNaC's integration kernels for iterated integrals.
Definition of GiNaC's lst.
Interface to GiNaC's products of expressions.
static ex iterated_integral3_evalf(const ex &kernel_lst, const ex &lambda, const ex &N_trunc)
static ex iterated_integral2_eval(const ex &kernel_lst, const ex &lambda)
const numeric pow(const numeric &x, const numeric &y)
static ex EllipticE_evalf(const ex &k)
static void EllipticK_print_latex(const ex &k, const print_context &c)
container< std::list > lst
static ex EllipticE_series(const ex &k, const relational &rel, int order, unsigned options)
static ex iterated_integral3_eval(const ex &kernel_lst, const ex &lambda, const ex &N_trunc)
static ex EllipticK_series(const ex &k, const relational &rel, int order, unsigned options)
std::vector< expair > epvector
expair-vector
const numeric abs(const numeric &x)
Absolute value.
function iterated_integral(const T1 &kernel_lst, const T2 &lambda)
const numeric sqrt(const numeric &x)
Numeric square root.
attribute_pure const T & ex_to(const ex &e)
Return a reference to the basic-derived class T object embedded in an expression.
static ex iterated_integral_evalf_impl(const ex &kernel_lst, const ex &lambda, const ex &N_trunc)
static void EllipticE_print_latex(const ex &k, const print_context &c)
static ex EllipticE_eval(const ex &k)
const constant Pi("Pi", PiEvalf, "\\pi", domain::positive)
Pi.
static ex EllipticK_eval(const ex &k)
static ex EllipticE_deriv(const ex &k, unsigned deriv_param)
bool is_a(const basic &obj)
Check if obj is a T, including base classes.
static ex EllipticK_evalf(const ex &k)
static ex EllipticK_deriv(const ex &k, unsigned deriv_param)
_numeric_digits Digits
Accuracy in decimal digits.
ex coeff(const ex &thisex, const ex &s, int n=1)
static ex iterated_integral2_evalf(const ex &kernel_lst, const ex &lambda)
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.
Interface to GiNaC's symbolic objects.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...
Utilities for summing over multiple indices.
Interface to GiNaC's wildcard objects.