Re: [Bug-gsl] [bug #43256] gsl_sf_coupling_6j overflows

[Alexey A. Illarionov]

Hi

Currently all coupling routines in GSL are quite primitive and produce
quite poor result. And there are algorithms that could improve
precision calculations dramatically, but... to this moment it is
unknown how to calculate coupling coefficients for large arguments
without involving big integers. At some point we would always deal
with overflows and precision loss. Unless there is a will to make GMP
a dependence of GSL we are stuck with this "feature".

There is a good paper on coupling coefficient algorithms by Tuzun et
al in Comp Phys Comm 112 (1998)

P.S. If someone interested I could transform my code to GSL patch that
calculates coupling_3j (not coupling_6j) using "edge recursive"
algorithm. It uses only double precision numbers and produces quite
accurate values up to sum(j_1+j_2+j_3) around 110 -- within 1e-12
relative accuracy. For comparison current GSL implementation gives
only max(rel_acc) = 1e-4 at the very same sum(j_1+j_2+j_3) = 110

Patrick Alken wrote:
> Thanks for the report. If you are able to fix the issue and send a
> patch that would be great, otherwise it could take a while for one
> of us to look into this.
> 
> On 09/19/2014 06:45 AM, Anders Søndergaard wrote:
>> URL: <http://savannah.gnu.org/bugs/?43256>
>> 
>> Summary: gsl_sf_coupling_6j overflows Project: GNU Scientific
>> Library Submitted by: andersas Submitted on: Fri 19 Sep 2014
>> 12:45:52 PM GMT Category: Runtime error Severity: 3 - Normal 
>> Operating System: Ubuntu linux Status: None Assigned to: None 
>> Open/Closed: Open Release: 1.16 (ubuntu latest) Discussion Lock:
>> Any
>> 
>> _______________________________________________________
>> 
>> Details:
>> 
>> I need to calculate 6-j symbols for somewhat large J. For example
>> the symbol:
>> 
>> (14, 16, 16 78, 62, 76) =
>> 6/1739*sqrt(16971005238954/1556382731177197) = 0.000360286
>> 
>> but this fails because gsl_sf_fact overflows. It is used
>> internally in gsl_sf_coupling_6j.
>> 
>> 
>> I have looked a bit at the code, and I think in these large cases
>> it is possible to work with the logarithm of the factorials and
>> in the end exponentiate.
>> 
>> Other libraries I have tested also gets these symbols wrong (but
>> silently returns the wrong value), except pythons sympy, which
>> calculates the symbols as the square root of a rational number.
>> 
>> I would rather not interface my Fortran program with Python.
>> 
>> 
>> 
>> 
>> 
>> 
>> _______________________________________________________
>> 
>> Reply to this item at:
>> 
>> <http://savannah.gnu.org/bugs/?43256>
>> 
>> _______________________________________________ Message sent
>> via/by Savannah http://savannah.gnu.org/
>> 
>> 
> 
> 


-- 
Regards, Alexey A Illarionov
С уважением, Алексей Александрович Илларионов.