11

This is a root-finding problem for an analytic function over a single dimension. One standard technique is to approximate $F(k)$ as a polynomial $F(k) \approx c_0 + c_1 k + c_2 k^2 + \cdots$ using a Chebyshev approximation, and to compute the roots of the polynomial semi-analytically, e.g. by setting up the companion matrix and computing its eigenvalues. (...


10

There is a prime missing both in your code and in the expression in your question. In the original paper, the expression is: $$ \log(x/x_0) + (x-x_0) \frac{\sum^\prime_{0\leq j\leq n}p_j T_j^*(x/6)}{\sum^\prime_{0\leq j\leq n}q_j T_j^*(x/6)}. $$ This primed sum is very common and standard when dealing with Chebyshev series: it means that the first term of ...


10

Doing one-off best rational approximations can often be accomplished by "manual" iterations of the Remez algorithm: interpolate a rational approximation with (relative or absolute) alternating sign errors at an initial guess for interpolation points, locate one (or more) points where the actual error exceeds that of the interpolation points and pivot (...


9

The integral in question is also known as the Boys function, after the British chemist Samuel Francis Boys who introduced its use in the early 1950s. A few years ago, I needed to compute this function in double precision, as fast as possible but accurately. I managed to achieve a relative error on the order of $10^{-15}$ across the entire input domain. It ...


9

The inverse Langevin function $\mathcal{L}^{-1}(x)$ is an odd function. Therefore one needs to consider only approximation on the interval $[0, 1]$; the negative half-plane is treated via symmetry about the origin. Because of the singularity at unity, a polynomial approximation is not a fruitful direction of the investigation, although some authors have ...


8

I implemented a complex polygamma function in the Julia standard library (https://github.com/JuliaLang/julia/pull/7125). It is possible to call this from C++ by linking to the libjulia library. The easiest way to do this is to call Julia's @cfunction construct to get a C/C++ function pointer to the compiled function. Then you can call it normally like any ...


8

Cube roots are not nearly as important as square roots (e.g., for normalizing vectors), so that might be why they are discussed much less. In general, if you apply Newton's method to $x^\alpha-\beta$, you get the iteration $$ (1-\alpha^{-1})x+\frac\beta\alpha x^{1-\alpha}, $$ so you just need to pick the equation so that $\alpha$ is negative integer, it's ...


8

I will complement @Richard Zhang 's answer (+1) with a python implementation of his suggested approach. The MATLAB package Chebfun has been partially ported in python. Actually there are two versions available: chebpy and pychebfun. Here's an implementation of the root finding procedure with pychebfun (the approach is similar with chebpy) from scipy....


7

The appropriate and fastest library depends on several things. Which Bessel functions (only J, Y & Hankel or modified Bessel functions I & K too), for which types of arguments (real or complex, integer, fractional or general order)? Amos's libraries are written in Fortran-77 (there are Fortran-90 coverted versions of TOMS 644 on a mirror of Alan ...


7

First of all, from the first paragraph of your attempts at a solution, I assume that the $z_j$ are non-negative? In that case, the integrand has no real problematic points (it's monotonous, decreasing). The only difficulty is the infinite integration range. And then it all depends on the accuracy you want... Do you want machine epsilon precision (or there ...


7

You try to be too literal: These functions are defined using an infinite sum, but of course computers can not execute infinite sums. Rather, they need to truncate any sums after a finite number of terms so that they can return a result after a finite amount of time. But even then, functions like yours might be defined this way, but that's not how they are ...


6

After some searching I found two formulas to solve the accuracy problem. The two modifications are In range $|k| \ge 40$ the (asymptotic) series from http://functions.wolfram.com/08.01.06.0010.02 is used (it is valid for $|k|>1,\,$ but convergence is suboptimal for small $k$): \begin{multline*} E(z)=\sqrt{-z} +\frac{\ln(-z)}{\sqrt{-z}}\sum_{n=0}^\infty \...


6

I hope I understood the question correctly. They try to compute exactly the same thing, so they really are equivalent. I'll use Chebyshev polynomials because they are easy to analyze. Given a function $f(x)$ on $[-1,1]$, the spectral interpolant is the truncation of $$ \begin{aligned} f(x) &= \sum_{n\geq0} \bar a_n T_n(x), \\ \bar a_n &= \frac{1+[n&...


6

You can avoid overflow by rewriting the sum appropriately. Say you have a sum $$ X = \sum_{i=1}^{n} x_i, $$ where some $x_i$'s are positive and large enough to cause overflow, but $\log x_i$ are of reasonable magnitude. If you find the greatest $\log x_M \geq \log x_i$, then you can rewrite the sum as $$ \log X = \log \sum_i e^{\log x_i} = \log x_M + \log \...


6

The problem is that np.cos(t) and np.sqrt(t) generate arrays with the length of t, whereas the second row ([0,1]) maintains the same size. To use np.vectorize with your function, you have to define the output type, and np.vectorize isn't really meant as a decorator except for the simplest cases. In this way however you can generate the function with the ...


6

You're probably better served writing a C wrapper for the Fortran implementation you linked to: Colavecchia, F. D., Gasaneo, G., "f1: a code to compute Appell's F1 hypergeometric function", Computer Physics Communications, Volume 157, Issue 1, p. 32-38 (2004), found at http://cpc.cs.qub.ac.uk/summaries/ADSJ. The R package appell wraps that implementation of ...


6

For the Hankel transform, one can classify the methods into four major groups: Numerical quadrature-based. Fourier-based ones. Asymptotic expansion of Bessel into sines and cosines. Projection-slice methods. The following paper gives a nice overview of these methods (types of methods): M. J. Cree and P. J. Bones, "Algorithms to numerically evaluate the ...


5

As you can see in this plot of $\log|f(\lambda)|$, $$ f(\lambda) = J_{\lambda-1}(1) - 2J_\lambda(1) -J_{\lambda+1}(1), $$ the roots $\lambda_k$ are really regular, and are approximately equal to $-k$ (starting from $k\geq0$, $\lambda_0=1.23219$ is an exception). So the way to get the $k$-th eigenvalue ($k\geq1$), is to bracket the root in $[-k-\frac12,-k+\...


5

Looking at what the python scipy library does for its special functions, the polygamma is found by returning the digamma if the zeroth derivative is requested, otherwise return $(-1)^{n+1}\Gamma(n+1)\zeta(n+1,z)$ where $\Gamma$ is the gamma function and $\zeta$ the two argument Riemann zeta function. Assuming that this identity holds for the complex numbers, ...


5

Step 1: analyze the recurrence I implemented the recurrence in Python, using basic numpy and scipy.special for the erf function. Why? Because it is simple and cheap (for reasonable $p$). import numpy as np import scipy.special as sp def F1(z): return np.sqrt(np.pi/2)*(1.0+sp.erf(z/np.sqrt(2))) def F2(z): return F1(z)*z + np.exp(-z**2/2.0) def F(...


5

Complementing the answers, I would like to say that directly discretizing the differential equation might work really well. I used that approach in a previous question. I used Finite Differences and used an eigensolver to find the roots. That, besides a perturbation analysis followed by Newton's method. The suggestion there was to use Chebyshev ...


5

There is no need for numerical computation here. First, $T(q)$ is a well-known function, the logarithmic integral. Repeated integration by parts gives an asymptotic expansion $$\mathrm{Li}(q) = \frac{q}{\log q}\sum_{k=0}^{K-1} \frac{k!}{\log^k q} + O\left(\frac{q}{\log^{K+1}q}\right).$$ There's also a fairly rapidly convergent representation due to ...


4

Philips in 1988 proved the following relationship: If $f(x)$ is an infinitely differentiable function defined on the interval $[-1,+1]$ and its Legendre expansion is given by $f(x) = \sum_{n=0}^{\infty}a_{n}P_{n}(x)$, then the Legendre coefficients $a_{n}^{(q)}$ of the $q$-th derivative of $f(x)$ are given by $$a_{n}^{(q)} = \frac{(2n+1)}{2^{q-2}(q-1)!} \...


4

This should be obvious, but the first thing you should try is to use the normalized incomplete gamma function $Q(a,x)=\frac{\Gamma(a,x)}{\Gamma(a)}$ for the two summands in the numerator of your fraction (for the first you have to add one recurrence step for $k+1$). In Boost $Q(a,x)$ is called gamma_q.


4

One possibility is to wrap the arbitrary precision interval implementation in Arb. This will not be as fast as a dedicated double precision implementation, but it might still be fast enough. Note 1: this code requires Arb version 2.8.0 or later. #include "acb_hypgeom.h" void hyp1f1ix(double * re, double * im, double a, double b, double x) { long prec; ...


4

I cannot think of any useful trigonometric identity that could help evaluating \begin{equation} a = \tan ( f \tan^{-1} g) \end{equation} so I would try a series expansion \begin{equation} fg + \frac{1}{3}f(f^2-1)g^3 + \frac{1}{15}f(2 f^4 - 5f^2 + 3) g^5 + \ldots \end{equation} but I doubt that if focus is on accuracy and robustness you can do much better ...


4

I think the main issue that makes this tricky is that trying to evaluate the power series for the hypergeometric function is numerically unstable due to catastrophic cancellation for one negative parameter and function value $\ll1$. Otherwise, having real parameters, and a real argument in $[0,1]$ should be fine. If you look at the DLMF entry for ...


4

This is a classical problem in numerical methods research: evaluating the zeros of special functions. Many years of research have gone into devising efficient methods. The canonical starting point for anything related to this is the classic book by Abramowicz and Stegun. Scanned versions of the book can easily be found on the internet (I'm not sure if these ...


4

You write $S(q)$ and $T(q)$ as integrals, but it is easier to think of them as solutions of ODEs: $$ S'(q) = \sin^2\left(\frac{π\Gamma(q)}{2q}\right) $$ with initial conditions $$ S(2) = 0, $$ and similarly for $T(q)$. You can then use any of the common ODE integrators in matlab, mathematica, maple, ..., to solve and plot the solutions so that you can ...


3

You can find the recurrence relations for the first derivative (with respect to the variable) in the Digital Library of Mathematical Functions (DLMF) in Chapter 14, sections 14.10 and 14.11. You'll have to take care of the stability of the recurrence, though. You should consult the references mentioned in this chapter (see final section on software and the ...


Only top voted, non community-wiki answers of a minimum length are eligible