Re: Numerical instability of j_1(x)=sin(Pi*x)^2+sin(Pi*(gamma(x)+1)/x)^2

* Georgi Guninski [2023-02-01 08:51]:
> The function:
> 
> j_1(x)=sin(Pi*x)^2+sin(Pi*(gamma(x)+1)/x)^2
> 
> has real zeros at the primes for x>1.
[...]
> ? j_1(40.0+20.*I)
> fails in pari, but I don't care about this.
[...]

For the record, the failure is related to

(10:05) gp > sin(10^40 * I)
  ***   at top-level: sin(10^40*I)
  ***                 ^------------
  *** sin: overflow in expo().

The result is about -exp(10^40) / (2*I), but the 'exponent' of 10^40 is
so large that the resulting floating point number can't be represented in PARI.
(Exponents are coded by C longs, so the hardcoded limit for exponents is 2^63;
i.e. floats about 2^(2^63) ...)

Cheers,

    K.B.
--
Pr Karim Belabas, U. Bordeaux, Vice-président en charge du Numérique
Institut de Mathématiques de Bordeaux UMR 5251 - (+33) 05 40 00 29 77
http://www.math.u-bordeaux.fr/~kbelabas/
`