83 cln::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);
118 cln::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();
172 cln::cl_N kbar =
sqrt(1-square(ex_to<numeric>(
k).to_cl_N()));
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)
216 ser +=
pseries(rel, std::move(nseq));
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();
252 cln::cl_N kbar =
sqrt(1-square(ex_to<numeric>(
k).to_cl_N()));
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)
300 ser +=
pseries(rel, std::move(nseq));
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());
342 cln::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++) {
397 multi_iterator_ordered_eq<int> i_multi(1,N+1,depth-1);
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;
424 cln::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) {
437 if ( is_a<basic_log_kernel>(*it) ) {
442 kernel.push_back( &ex_to<integration_kernel>(*it) );
447 cln::cl_N result = iterated_integral_do_sum(
m, kernel, lambda, N_trunc);
454 cln::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++) {
462 if ( !(is_a<basic_log_kernel>(*(kernel[i]))) ) {
467 if ( i_trailing == 0 ) {
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);
485 multi_iterator_shuffle_prime<const integration_kernel *> i_shuffle(a,b);
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);
516 lst k_lst = ex_to<lst>(kernel_lst);
518 bool flag_not_numeric =
false;
519 for (
const auto & it : k_lst) {
520 if ( !is_a<integration_kernel>(it) ) {
521 flag_not_numeric =
true;
524 if ( flag_not_numeric) {
528 for (
const auto & it : k_lst) {
529 if ( !(ex_to<integration_kernel>(it).is_numeric()) ) {
530 flag_not_numeric =
true;
533 if ( flag_not_numeric) {
539 int N_trunc_int = ex_to<numeric>(N_trunc).
to_int();
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) {
548 kernel_vec.push_back( &ex_to<integration_kernel>(it) );
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));
561 cln::cl_N lambda_cln = ex_to<numeric>(lambda.
evalf()).to_cl_N();
565 if ( is_a<basic_log_kernel>(*(kernel_vec[depth-1])) ) {
569 result = iterated_integral_shuffle(kernel_vec, lambda_cln, N_trunc_int);
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);
580 if ( !(is_a<basic_log_kernel>(*(kernel_vec[i_plus-2]))) ) {
581 result +=
coeff * iterated_integral_shuffle(kernel_vec, lambda_cln, N_trunc_int);
621 do_not_evalf_params().
628 do_not_evalf_params().
Exception class thrown by classes which provide their own series expansion to signal that ordinary Ta...
Interface to GiNaC's symbolic exponentiation (basis^exponent).
static ex iterated_integral_evalf_impl(const ex &kernel_lst, const ex &lambda, const ex &N_trunc)
Interface to GiNaC's symbolic objects.
int to_int(const numeric &x)
ex subs(const exmap &m, unsigned options=0) const
static ex EllipticK_eval(const ex &k)
const basic & hold() const
Stop further evaluation.
static ex EllipticK_evalf(const ex &k)
static ex iterated_integral2_evalf(const ex &kernel_lst, const ex &lambda)
static ex EllipticK_deriv(const ex &k, unsigned deriv_param)
Interface to GiNaC's constant types and some special constants.
static ex EllipticK_series(const ex &k, const relational &rel, int order, unsigned options)
This class holds a relation consisting of two expressions and a logical relation between them...
Interface to GiNaC's sums of expressions.
Interface to GiNaC's products of expressions.
Interface to GiNaC's wildcard objects.
bool info(unsigned inf) const
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...
static ex EllipticE_deriv(const ex &k, unsigned deriv_param)
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
static ex EllipticE_eval(const ex &k)
static void EllipticK_print_latex(const ex &k, const print_context &c)
REGISTER_FUNCTION(conjugate_function, eval_func(conjugate_eval). evalf_func(conjugate_evalf). expl_derivative_func(conjugate_expl_derivative). info_func(conjugate_info). print_func< print_latex >(conjugate_print_latex). conjugate_func(conjugate_conjugate). real_part_func(conjugate_real_part). imag_part_func(conjugate_imag_part). set_name("conjugate","conjugate"))
ex evalf() const override
Evaluate object numerically.
This class holds a extended truncated power series (positive and negative integer powers)...
_numeric_digits Digits
Accuracy in decimal digits.
static void EllipticE_print_latex(const ex &k, const print_context &c)
const numeric log(const numeric &x)
Natural logarithm.
static ex EllipticE_evalf(const ex &k)
static ex EllipticE_series(const ex &k, const relational &rel, int order, unsigned options)
const numeric abs(const numeric &x)
Absolute value.
The basic integration kernel with a logarithmic singularity at the origin.
Interface to GiNaC's overloaded operators.
static ex iterated_integral3_evalf(const ex &kernel_lst, const ex &lambda, const ex &N_trunc)
const numeric sqrt(const numeric &x)
Numeric square root.
const numeric binomial(const numeric &n, const numeric &k)
The Binomial coefficients.
Base class for print_contexts.
ex series(const ex &r, int order, unsigned options=0) const
Compute the truncated series expansion of an expression.
function iterated_integral(const T1 &kernel_lst, const T2 &lambda)
Definition of GiNaC's lst.
const numeric pow(const numeric &x, const numeric &y)
Lightweight wrapper for GiNaC's symbolic objects.
const constant Pi("Pi", PiEvalf, "\i", domain::positive)
Pi.
Interface to GiNaC's initially known functions.
ex coeff(const ex &thisex, const ex &s, int n=1)
Utilities for summing over multiple indices.
Interface to class for extended truncated power series.
Wrapper template for making GiNaC classes out of STL containers.
static ex iterated_integral3_eval(const ex &kernel_lst, const ex &lambda, const ex &N_trunc)
Interface to relations between expressions.
Makes the interface to the underlying bignum package available.
ex evalf() const override
Evaluate object numerically.
static ex iterated_integral2_eval(const ex &kernel_lst, const ex &lambda)
std::vector< expair > epvector
expair-vector
const numeric factorial(const numeric &n)
Factorial combinatorial function.
Interface to GiNaC's integration kernels for iterated integrals.
static unsigned register_new(function_options const &opt)