1
$\begingroup$

Cross posted on StackOverflow

I am trying to use the gsl library for calculating 1F1. I have some C code.

The following works and matches Mathematica's results for small real arguments z but not large arguments:

/* file 1F1.c */
#include <stdio.h>
#include <stdlib.h>

#include <gsl/gsl_sf.h>
#include <gsl/gsl_math.h>

double _1F1( double b,
            double c,
            double z)
{
  double Fa = gsl_sf_hyperg_1F1(b, c, z);
  return Fa;
}

int main()
{

  double b = 100000.0;
  double c = 14.0;
  double z = 0.9;
  double result;   
  
  result = _1F1(b, c, z);
  printf("1F1(%.1f, %.1f, %.1f): %e\n", b, c, z, result);

  /* result = 3.2042*10^236 in Mathematica*/

  z = 10*z;

  result = _1F1(b, c, z);
  printf("1F1(%.1f, %.1f, %.1f): %e\n", b, c, z, result);       

  /* result = 9.52444189237289*10^794 in Mathematica*/

  return 0;
}

I get:

$    gcc -std=c17 -Wall -pedantic 1F1.c -lm -lgsl -lcblas; ./a.out
1F1(100000.0, 14.0, 0.9): 3.204195e+236
1F1(100000.0, 14.0, 9.0): -nan

So, how do I get around and evaluate for large z, as Mathematica seems to be able to do?

I also tried using the Cephes library version and got into the same problem.

I guess I am trying to understand why Mathematica does not have this problem, i.e. what approximation formula for large arguments they are using here. Is there something else I can use here?

$\endgroup$
16
  • $\begingroup$ I don't think you can even represent the result as a double, even if the routine internally were to use some higher precision. The largest double is approximately $1.8 \times 10^308$. Mathematica maybe uses some kind of arbitrary precision arithmetic or just quadruple precision. $\endgroup$
    – Tyberius
    Nov 11, 2023 at 0:36
  • 1
    $\begingroup$ I think @Tyberius is on the right track. The 1F1 function grows exponentially. Just because Mathematica gives you some number does not mean that you should expect it to be correct or accurate. $\endgroup$ Nov 11, 2023 at 2:38
  • 2
    $\begingroup$ @user3236841 A better approach is to find an implementation to gives you the log of 1F1. That's a more reasonable function for which you can reason about correctness without having to worry about floating point overflow. $\endgroup$ Nov 11, 2023 at 21:48
  • 2
    $\begingroup$ Abramowicz-Stegun is the historically first place I would look. Today that is all collected in a website that is called something like the NIST handbook of special functions or similar. $\endgroup$ Nov 12, 2023 at 4:16
  • 1
    $\begingroup$ I checked your example in mpmath and obtained the same result: python >>> import mpmath >>> from mpmath import * >>> hyp1f1(100000.0, 14.0, 9.0) mpf('9.5244418923728841e+794') $\endgroup$
    – nicoguaro
    Nov 13, 2023 at 16:01

1 Answer 1

2
$\begingroup$

So, how do I get around and evaluate for large z, as Mathematica seems to be able to do?

I think that they use arbitrary precision arithmetic (see this).

If I check this computation in Maxima (I tried it on my phone, though)

hypergeometric([100000.0], [14.0], 9.0)

I get the following

$$9.524441892372884_B \times 10^{794}$$

If I try it on mpmath

>>> import mpmath
>>> from mpmath import *
>>> hyp1f1(100000.0, 14.0, 9.0)

I get the following

mpf('9.5244418923728841e+794')

This is the same result as in SymPy since it uses mpmath internally.

Is there something else I can use here?

You would need an implementation of special functions and arbitrary precision at the same time. I do not know of C packages for that off the top of my head, though.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.