Re: [Bug-gsl] bugs when using gsl_sf_hyperg_U

[Patrick Alken]

Just wanted to forward this to the list so it is archived

On 03/11/2014 03:42 PM, Feng Gao wrote:
> Hi Ray and Patrick,
>
> Thanks for your help so much thus far. I have included a simple cpp 
> file which demonstrates several cases that will invoke bugs. I put the 
> gsl folder in the same current folder and used
>
> g++ -Wall -o gsl_sf_hyperg_U_bug_report gsl_sf_hyperg_U_bug_report.cpp 
> -lgsl -lgslcblas -lm
>
> to compile.
>
> Thanks again!
>
> Best Regards,
>
> Feng
>
>
> 2014-03-11 15:13 GMT-04:00 Raymond Rogers <address@hidden 
> <mailto:address@hidden>>:
>
>     Hi Feng,
>
>         Take the writing below as coming from somebody who had just a
>     brush with the subject matter; numerical approximatioin of special
>     function.s.  I think that somebody on the GSL mailing list is an
>     expert and could actually help.  I didn't have much luck but that
>     was years ago.  I do suggest giving actual c (or some such) code
>     you  are trying to run.  A simple program.
>     Let me know I if I can help further.
>     
> ------------------------------------------------------------------------------------------
>         No joy in my searh but I did look at the DLMF.  In particular
>     http://dlmf.nist.gov/13.2#vii
>     To change  -b to +b  (and +a)
>     and
>     http://dlmf.nist.gov/13.8
>     for Asymptotic expansions.
>     Unfortunately I didn't find your case where the a,b go to
>     infinity; but I did find an earlier paper:
>     http://www.sciencedirect.com/science/article/pii/S0377042796001331
>     Which directly that covers b-> -inf   equation 3. Unfortunately it
>     solves the problem in terms of Bessel functions :)  So some work
>     has to be done.
>
>     When I worked with large negative a I used a recursive form to
>     walk a up to positive.  One of the relations on
>     http://dlmf.nist.gov/13.3
>     In particular 13.3.8   walks b and leaves a, z  fixed.
>     I did it in a matrix form because I am simpleminded but a straight
>     recursion would work.
>     Be aware that there are two ways to work recursions one explodes
>     errors (the dominant recursion) and the other smooths them out
>     (the minimall one).
>     If you have access to a library or have "Numerical Methods for
>     Special Functions" Gil,Segura, Temme chapter 4 covers recurrence
>     theory.  It's really not very hard.  If needed I could copy that
>     chapter and email it or you might get a preview on Google books
>
>     I did email Pr Temme  (who is a scholar on numerical
>     approximations and Special functions) at:
>     http://homepages.cwi.nl/~nicot/ <http://homepages.cwi.nl/%7Enicot/>
>     He was very helpful and came back in 1 day with an answer to my
>     question about large negative a;  a simple formula that worked
>     with a<-1000 to astounding accuracy.  I imagine asking him for a
>     pointer instead of just asking for the formula would be more
>     diplomatic.  I did go out and buy one of his books after he
>     answered  and read some parts.
>
>
>     Ray
>
>
>
>
>
>     On 03/11/2014 12:15 PM, Feng Gao wrote:
> >     Hi Ray,
> >
> >     Thanks very much for your effort!
> >
> >     Best Regards,
> >
> >     Feng
> >
> >
> >     2014-03-11 11:56 GMT-04:00 Raymond Rogers
> >     <address@hidden <mailto:address@hidden>>:
> >
> >         Okay, let me look around on my computer.  I have to find it;
> >         the original email I had got extinguished by difficult
> >         manuever.  It's hard but you _can_ erase all of your old
> >         emails onn gmail if you really try (or are really dumb!).
> >         Anyhow.... I will check now.
> >         Ray
> >
> >         On 03/11/2014 11:01 AM, Feng Gao wrote:
> > >         Hi Ray,
> > >
> > >         Thanks for your response! These would be helpful indeed.
> > >         Please notice that all the problems I had were for negative
> > >         b. Actually in my case the value of a is fixed at 1 (so I
> > >         don't need to worry about large negative a); the value of b
> > >         is (-inf, 2) and the value of x is > 0. If you have anything
> > >         that works for this scenario, that will be extremely helpful.
> > >
> > >         Best Regards,
> > >
> > >         Feng
> > >
> > >
> > >         2014-03-11 9:03 GMT-04:00 Raymond Rogers
> > >         <address@hidden
> > >         <mailto:address@hidden>>:
> > >
> > >             Hello,
> > >                 I had a similar problem with negative a (Lamm's
> > >             equation) and developed some programs that worked.  
> > >             That was about three years ago and I would have to find
> > >             it.  In addition the program blows up for large negative
> > >             a; but I got in touch with a professional in these
> > >             things and he gave me a asymptotic formula.  I also have
> > >             his book that explains a lot about doing these calculations.
> > >                 Let me know if I can help along this line.  I
> > >             probably have a bunch of links for asymptotic analysis.
> > >
> > >             Ray
> > >
> > >
> > >             On 03/10/2014 10:38 PM, Feng Gao wrote:
> > >
> > >                 Hi,
> > >
> > >                 Sorry for so many emails... I also tried the GSL
> > >                 SHELL software and it
> > >                 seemse that sf.hypergU(1, -500, 0.1) can be
> > >                 calculated correctly. However,
> > >                 I found that any negative non-integers of b for U(a,
> > >                 b, x) cannot be
> > >                 calculated. These negative non-integers should be
> > >                 suitable arguments and
> > >                 the manual states that b is of type double.
> > >
> > >                 Thanks very much for your attention and I will be
> > >                 very grateful if you can
> > >                 resolve these problems.
> > >
> > >                 Best Regards,
> > >
> > >                 Feng
> > >
> > >
> > >                 2014-03-10 21:01 GMT-04:00 Feng Gao
> > >                 <address@hidden <mailto:address@hidden>>:
> > >
> > >                     Sorry for the typo but for example 1) listed in
> > >                     the previous email:
> > >
> > >                     gsl_sf_hyperg_U(1, -500, 0.1) gives the correct
> > >                     value: ~1.99e-3.
> > >
> > >                     The problem happens when calculating
> > >                     gsl_sf_hyperg_U(1, -500, 40): gsl:
> > >                     gsl: hyperg.c:165: ERROR: overflow. The value
> > >                     should be ~0.00184 so I
> > >                     suppose this is a bug.
> > >
> > >                     Again, thanks very much if you can help me solve
> > >                     this problem since I very
> > >                     much need it in my research.
> > >
> > >                     Best Regards,
> > >
> > >                     Feng
> > >
> > >
> > >                     2014-03-10 19:53 GMT-04:00 Feng Gao
> > >                     <address@hidden
> > >                     <mailto:address@hidden>>:
> > >
> > >                     Hi,
> > >
> > >                         I have been using the gsl_sf_hyperg_U
> > >                         function for a while and I found
> > >                         several bugs when calculating normal values.
> > >                         For example,
> > >
> > >                         1) gsl_sf_hyperg_U(1, -500, 0.1) gives
> > >                         incorrect value: 9.3246e+03, which
> > >                         actually should be around 0.19.
> > >
> > >                         2) Many domain error bugs. For example,
> > >                         gsl_sf_hyperg_U(1, -6.67, 1)
> > >                         should have value ~0.128 but will invoke
> > >                         domain error for gsl.
> > >
> > >                         Hope that you can further investigate this
> > >                         and come up with solutions.
> > >                         Thanks!
> > >
> > >                         Best Regards,
> > >
> > >                         Feng
> > >
> > >
> > >
> > >             -- 
> > >             Act IV, Sc. IV
> > >             What is a man,
> > >             If his chief good and profit of his time
> > >             Be to sleep and feed.  Be a beast, no more
> >
> >         -- 
> >         Act IV, Sc. IV
> >         What is a man,
> >         If his chief good and profit of his time
> >         Be to sleep and feed.  Be a beast, no more
>
>     -- 
>     Act IV, Sc. IV
>     What is a man,
>     If his chief good and profit of his time
>     Be to sleep and feed.  Be a beast, no more