1

I need to evaluate a product of two real functions, namely $F(x)\cdot G(x)$. The function $F(x)$ is a diverging and obscure member of scipy.special, while $G(x)$ is a gaussian. While the product is a well behaved function and of order ~1, its factors by themselves, for x range that I am interested in, are of order $10^{\pm 500}$. My problem is that each of them overflow pythons float. My go-to move is to try numpy.float128 which would cover range of interest quite nicely, if not for the error:

TypeError: ufunc 'pbdv' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

So, function $F$ doesn't support float128. Is there a well known path to deal with this? Only thing I can think of right now is to express offending function as power series and calculate Cauchy product of sums - would that work?

i_prob_should_know_this
  • 141
  • 3
  • 2
    The typical solution is to work in log space, potentially tracking signs separately if you can't guarantee positivity. Taking the log of a Gaussian is easy, but I have no idea what properties of F might make this easily exploitable. What function is F? – helloworld922 Aug 26 '22 at 14:57
  • Its https://reference.wolfram.com/language/ref/ParabolicCylinderD.html . I did think of that once, but discarded it on a hunch - let me take another look, thanks for suggestion. – i_prob_should_know_this Aug 26 '22 at 15:12
  • Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. – Community Aug 26 '22 at 16:31
  • What is the "diverging and obscure member of scipy.special" specifically? – njuffa Aug 27 '22 at 08:28
  • It is a special function, parabolic cylindrical function (scipy.special.pbdv). More info here: https://reference.wolfram.com/language/ref/ParabolicCylinderD.html – i_prob_should_know_this Aug 27 '22 at 20:17
  • Often libraries of special functions whose values reach infinitesimal or exponentially high ranges will offer versions that are normalized e.g. by exact powers of two, to cope with just such intermediate overflow errors. It's hard to say more without looking into your specifics. – hardmath Aug 28 '22 at 16:01

0 Answers0