How to implement 0F1 and 1F1 for double precision

[Ondřej Čertík]

Hi Fredrik,

I need an efficient (in terms of speed) Fortran implementation of the
0F1 function,
I only need double precision, but it should be accurate as much as it
could (1e-15 or 1e-14
for all parameters would be nice). Since you have a lot of experience with this,
I'll ask here.

Here is my Fortran implementation of a related function that can be
expressed using 1F1:

https://gist.github.com/3810032

I just use the direct series for 1F1(a, x) as long as roughly "x < a"
and recursion relations
otherwise. Since I only need this for integer "a", this works for this
particular case.

I wanted to just reuse some other codes for this, for example the
Fortran implementation
in scipy, but it simply fails (silently!) for larger values, see for
example my bug report [1].
Mpmath works great, but it's slow.

I would like to have an implementation in Fortran that is:

a) accurate for all parameters, 1e-14 or 1e-15 is enough (probably
even 1e-13 is fine)
b) loudly fails if convergence cannot be achieved
c) as fast as possible

(I only need real parameters and "x" for now, also because the real arithmetic
is much faster than the complex one there need to be two implementations anyway,
one for real, one for complex.)

The SciPy's implementation fails b). Mpmath fails c). My
implementation for some special
values of 1F1 above satisfies all a), b), c). But I would like to have
general functions
for both 1F1 and 0F1 that work for all parameter values. Which means,
that even though
the 1F1 and 0F1 series work for all "x" (radius of convergence is
infinite), in practice there
are cancellation problems. So I looked at how mpmath does it, and it
switches to asymptotic
series for large "x" (if magz >= 8). But it's using hyp2f0 via the
identity [2], which has radius of convergence 0.
The beginning of function hyp2f0 is:

def _hyp2f0(ctx, a_s, b_s, z, **kwargs):
(a, atype), (b, btype) = a_s
# We want to try aggressively to use the asymptotic expansion,
# and fall back only when absolutely necessary
try:
kwargsb = kwargs.copy()
kwargsb['maxterms'] = kwargsb.get('maxterms', ctx.prec)
return ctx.hypsum(2, 0, (atype,btype), [a,b], z, **kwargsb)
except ctx.NoConvergence:
if kwargs.get('force_series'):
raise

but hypsum() seems to be just summing up the hypergeometric series for 2F0 here,
which has radius of convergence 0.
The force_series=True, so it will fail if it doesn't converge.

But I don't understand the logic here --- how can the series for 2F0 converge?
I noticed that ctx_fp.py just has a very simple implementation, but
ctx_mp.py is quite complex, it seems to be using some convergence accelerators?

Is the relation [2] the best general way to calculate 0F1 for large
values of "x" using
the asymptotic series?

What is the best method for calculating 1F1 for large "x"? Mpmath is using
the function hypercomb() from hypergeometric.py, but I am not quite sure what
it does, nor what exactly formula is implemented in the _hyp1f1 function.

Btw, I've put up my notes on hypergeometric series at [3], it shows
how to convert all basic
functions that I was aware of into hypergeometric series.

Thanks for any help or tips.

Ondrej

[1] http://mail.scipy.org/pipermail/scipy-user/2012-September/033189.html

[2] http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric0F1/06/02/03/0004/

[3] http://theoretical-physics.net/dev/src/math/hyper.html

[Fredrik Johansson]

Hi Ondrej

On Mon, Oct 1, 2012 at 9:38 AM, Ondřej Čertík <ondrej...@gmail.com> wrote:
> Hi Fredrik,

>
> but hypsum() seems to be just summing up the hypergeometric series for 2F0 here,
> which has radius of convergence 0.
> The force_series=True, so it will fail if it doesn't converge.

> But I don't understand the logic here --- how can the series for 2F0 converge?
> I noticed that ctx_fp.py just has a very simple implementation, but
> ctx_mp.py is quite complex, it seems to be using some convergence accelerators?

The 2F0 series "converges" when the terms become small enough. This is
just a heuristic. The complexity in hypsum is mostly intended to take
care of increasing the precision when parameters are close to negative
integers and/or there is significant cancellation. Convergence
acceleration is only used for 3F2, 4F3, etc. on the unit circle.

> Is the relation [2] the best general way to calculate 0F1 for large
> values of "x" using
> the asymptotic series?

More or less. The Wolfram Functions site also gives a trigonometric
version, which might be more convenient.

> What is the best method for calculating 1F1 for large "x"? Mpmath is using
> the function hypercomb() from hypergeometric.py, but I am not quite sure what
> it does, nor what exactly formula is implemented in the _hyp1f1 function.

It's just using this (very standard) asymptotic series:
http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F1/06/02/0005/

hypercomb() is documented here:
http://mpmath.googlecode.com/svn/trunk/doc/build/functions/hypergeometric.html#hypercomb

In general, you want to use the convergent series for small x and the
asymptotic series for large x. What to do in between when the
convergent series has too much cancellation is a very difficult
problem, unless you accept the overhead of switching to multiprecision
arithmetic. The DLMF (http://dlmf.nist.gov/13.29) suggests using
something called "exponentially-improved expansions". I have not tried
using them, so I can't say how well they work, but they are probably
well worth investigating in your case.

Fredrik

[Ondřej Čertík]

On Mon, Oct 1, 2012 at 2:06 AM, Fredrik Johansson
<fredrik....@gmail.com> wrote:
> Hi Ondrej
>
> On Mon, Oct 1, 2012 at 9:38 AM, Ondřej Čertík <ondrej...@gmail.com> wrote:
>> Hi Fredrik,
>>
>> but hypsum() seems to be just summing up the hypergeometric series for 2F0 here,
>> which has radius of convergence 0.
>> The force_series=True, so it will fail if it doesn't converge.
>
>> But I don't understand the logic here --- how can the series for 2F0 converge?
>> I noticed that ctx_fp.py just has a very simple implementation, but
>> ctx_mp.py is quite complex, it seems to be using some convergence accelerators?
>
> The 2F0 series "converges" when the terms become small enough. This is
> just a heuristic.

Ah I see -- that's why it is called "asymptotic series":

http://mathworld.wolfram.com/AsymptoticSeries.html

so what happens is that for a fixed "x", the series first (typically) converges
a bit, and then starts to diverge. The idea is that in the limit x->0
(in our case),
the series (cut after n terms) always gives the correct asymptotic
expression/limit.

So if "x" is small, we might get lucky to get accurate answer within
our precision.

But I would expect that if you want many digits, eventually you will
not be able to do it
with such a series.

> The complexity in hypsum is mostly intended to take
> care of increasing the precision when parameters are close to negative
> integers and/or there is significant cancellation. Convergence
> acceleration is only used for 3F2, 4F3, etc. on the unit circle.

I see. I saw your blog post:

http://fredrikj.net/blog/2010/07/euler-maclaurin-summation-of-hypergeometric-series/

describing it.

>
>> Is the relation [2] the best general way to calculate 0F1 for large
>> values of "x" using
>> the asymptotic series?
>
> More or less. The Wolfram Functions site also gives a trigonometric
> version, which might be more convenient.
>
>> What is the best method for calculating 1F1 for large "x"? Mpmath is using
>> the function hypercomb() from hypergeometric.py, but I am not quite sure what
>> it does, nor what exactly formula is implemented in the _hyp1f1 function.
>
> It's just using this (very standard) asymptotic series:
> http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F1/06/02/0005/
>
> hypercomb() is documented here:
> http://mpmath.googlecode.com/svn/trunk/doc/build/functions/hypergeometric.html#hypercomb

Ah I see, thanks.

>
> In general, you want to use the convergent series for small x and the
> asymptotic series for large x. What to do in between when the
> convergent series has too much cancellation is a very difficult
> problem, unless you accept the overhead of switching to multiprecision
> arithmetic. The DLMF (http://dlmf.nist.gov/13.29) suggests using
> something called "exponentially-improved expansions". I have not tried
> using them, so I can't say how well they work, but they are probably
> well worth investigating in your case.

I've seen those too, but then I went into mpmath and you didn't use them,
but now I can see, that you just increase the precision and that helps in many
cases. In Fortran, one can easily switch to quadruple precision. It's
just a lot slower.
So maybe for the difficult cases, switching to quadruple might do the job.

Thanks for the tips. So as I can see it, I first need to implement 2F0 for small
arguments using the asymptotic series and check that it converges.
Then implement 1F1 and 0F1 for small "x" directly and in terms
of 2F0 for large "x". For 0F1 this should work. For 1F1 I might hit
some cancellation problems, so I'll have to see --- so probably implementing
hypercomb() will be needed.

I can see, that hypsum is calling for example the MPF_hypsum in Sage.
This function works at a set precision, correct? And returns the "magn",
which you use later to determine if you need to rerun it with higher precision.
I am still having problems understanding how you detect, that too much
cancellation has occurred in ctx_mp.py:

zv, have_complex, magnitude = summator(coeffs, v, prec, wp, \
epsshift, mag_dict, **kwargs)
cancel = -magnitude
jumps_resolved = True
if extraprec < max_total_jump:
for n in mag_dict.values():
if (n is None) or (n < prec):
jumps_resolved = False
break
accurate = (cancel < extraprec-25-5 or not accurate_small)
if jumps_resolved:
if accurate:
break

the magnitude is returned from Sage (for example) by the following code:

MPF_set_fixed(a, SRE, wp, prec, ROUND_N)
if mpz_sgn(SRE):
magn = mpz_get_si(a.exp) + <long>mpz_bitcount(a.man)
else:
magn = -wp

what exactly is it? The SRE seems to be the final sum, which is
assigned to "a" (the return value)
and magn is somehow the number of bits of the number, but I don't
understand how this can
detect any cancellation. How can I detect cancellation
using just double precision floating point numbers?

I noticed that cancellation typically occurs when the individual terms
become (temporarily) big,
which reduces accuracy. So maybe I can just check the max absolute
size of the term,
relative to the final answer. If the final answer is let's say 0.1,
but it contains terms of order 1e6
that cancel out then I am in trouble. If on the other hand the final
answer is 1e10, then I should be fine.

Ondrej

[Fredrik Johansson]

On Mon, Oct 1, 2012 at 9:21 PM, Ondřej Čertík <ondrej...@gmail.com> wrote:
> On Mon, Oct 1, 2012 at 2:06 AM, Fredrik Johansson
> <fredrik....@gmail.com> wrote:
>> Hi Ondrej
>>
>> On Mon, Oct 1, 2012 at 9:38 AM, Ondřej Čertík <ondrej...@gmail.com> wrote:
>>> Hi Fredrik,
>>>
>>> but hypsum() seems to be just summing up the hypergeometric series for 2F0 here,
>>> which has radius of convergence 0.
>>> The force_series=True, so it will fail if it doesn't converge.
>>
>>> But I don't understand the logic here --- how can the series for 2F0 converge?
>>> I noticed that ctx_fp.py just has a very simple implementation, but
>>> ctx_mp.py is quite complex, it seems to be using some convergence accelerators?
>>
>> The 2F0 series "converges" when the terms become small enough. This is
>> just a heuristic.
>
> Ah I see -- that's why it is called "asymptotic series":
>
> http://mathworld.wolfram.com/AsymptoticSeries.html
>
> so what happens is that for a fixed "x", the series first (typically) converges
> a bit, and then starts to diverge. The idea is that in the limit x->0
> (in our case),
> the series (cut after n terms) always gives the correct asymptotic
> expression/limit.

Yes. (Assuming you mean |x|->inf, where x is the original variable of
the 0F1/1F1.)

> So if "x" is small, we might get lucky to get accurate answer within
> our precision.
>
> But I would expect that if you want many digits, eventually you will
> not be able to do it
> with such a series.

It still works with many digits, if x increases proportionally...

I haven't used them because increasing the precision works well
enough, and it would be much more complicated to have three algorithms
rather than just two. But they ought to be more useful if you only
have double precision, and selecting between algorithms is easier when
the precision is fixed.

> Thanks for the tips. So as I can see it, I first need to implement 2F0 for small
> arguments using the asymptotic series and check that it converges.
> Then implement 1F1 and 0F1 for small "x" directly and in terms
> of 2F0 for large "x". For 0F1 this should work. For 1F1 I might hit
> some cancellation problems, so I'll have to see --- so probably implementing
> hypercomb() will be needed.

hypercomb() is not very useful in double precision since it depends on
increasing the precision.

The cancellation estimation is specific to the use of fixed-point
arithmetic, which itself is an ugly hack that happens to be fast for
doing arbitrary-precision arithmetic with Python... I wouldn't look at
it for too long if I were you :-)

> I noticed that cancellation typically occurs when the individual terms
> become (temporarily) big,
> which reduces accuracy. So maybe I can just check the max absolute
> size of the term,
> relative to the final answer. If the final answer is let's say 0.1,
> but it contains terms of order 1e6
> that cancel out then I am in trouble. If on the other hand the final
> answer is 1e10, then I should be fine.

Yes, to measure the cancellation with floating point numbers, you can
simply store the (absolute value) largest term, and compare that to
the final result.

Fredrik