From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 23797 invoked by alias); 1 Oct 2002 18:40:49 -0000 Mailing-List: contact gsl-discuss-help@sources.redhat.com; run by ezmlm Precedence: bulk List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sources.redhat.com Received: (qmail 23790 invoked from network); 1 Oct 2002 18:40:49 -0000 Received: from unknown (HELO mail.phy.duke.edu) (152.3.182.2) by sources.redhat.com with SMTP; 1 Oct 2002 18:40:49 -0000 Received: from ganesh.phy.duke.edu (ganesh.phy.duke.edu [152.3.182.51]) by mail.phy.duke.edu (Postfix) with ESMTP id EE7BD30374 for ; Tue, 1 Oct 2002 14:40:48 -0400 (EDT) Received: by ganesh.phy.duke.edu (Postfix, from userid 1337) id 68BF0520A2; Tue, 1 Oct 2002 14:40:48 -0400 (EDT) Date: Tue, 01 Oct 2002 11:40:00 -0000 From: "Robert G. Brown" To: gsl-discuss@sources.redhat.com Subject: Laguerre Polynomial overflow... Message-ID: MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII X-SW-Source: 2002-q4/txt/msg00000.txt.bz2 Dear GSL listpersons: I'm using gsl-1.0-1 prebuilt in Red Hat 7.2 (and have been using various gsl routines from roughly 0.7-something on). Currently I'm testing an ODE-based solution to the time-independent Schrodinger equation by solving the O(3) simple harmonic oscillator problem, the radial eigensolution to which is basically a laguerre polynomial times a few other simple functions. To compare the ODE solution to the exact solution, I've been evaluating ... status = gsl_sf_laguerre_n_e(nsho, lsho+0.5 , rsho*rsho, &lgp); if(status != GSL_SUCCESS) return(0.0); ret = sqrt(2.0*fac(nsho)/fac(lsho + nsho + 1.0))*pow(rsho,(double)lsho)* lgp.val*exp(-rsho*rsho*0.5)/norm; fprintf(stderr,"sho3d( %2d, %2d, %2d, %.5e) = %.5e\n",nsho,lsho,msho,rsho,ret! fprintf(stderr,"gsl_sf_laguerre_n( %2d, %2d, %.5e) = %.5e\n",nsho,lsho,rsho,lgp.val); return(ret); which works fine for nsho = 0 or 1, but when I hit nsho > 2, the routine runs for some values of lsho and rsho but then dies on an overflow beyond a certain point: sho3d( 2, 0, 0, 2.21000e+00) = 7.38523e-02 gsl_sf_laguerre_n( 2, 0, 4.88410e+00) = 1.59197e+00 sho3d( 2, 0, 0, 2.22000e+00) = 7.70712e-02 gsl_sf_laguerre_n( 2, 0, 4.92840e+00) = 1.69856e+00 sho3d( 2, 0, 0, 2.23000e+00) = 8.02146e-02 gsl_sf_laguerre_n( 2, 0, 4.97290e+00) = 1.80762e+00 gsl: laguerre.c:96: ERROR: overflow Abort The interesting thing is that (as you can see) the solution is extremely pedestrian here -- n is two, \ell is 0, and the argument is around 5. The returned value is order unity. I can see no good reason for an overflow to be generated externally, and presume that it is occuring during the operation of e.g. the recursion relation that generates the result internally. I also seem to generate no GSL error condition -- the code crashes even though I try to trap the overflow and continue. For what it is worth, the solution is dead on up to that point, validated somewhat backwards by the nearly exact coincidence of this solution with the one I'm obtaining by explicit ODE solution with e.g. rk45 (also from gsl elsewhere in the code). If nobody on this list has encountered this or can suggest any fix beyond going into the source and figuring out where and why the overflow occurs and fixing it, I'm happy enough to do the latter, but given the hassle of building from source and the possibility that it has already been fixed in a more recent version, I thought I'd ask first. rgb Robert G. Brown http://www.phy.duke.edu/~rgb/ Duke University Dept. of Physics, Box 90305 Durham, N.C. 27708-0305 Phone: 1-919-660-2567 Fax: 919-660-2525 email:rgb@phy.duke.edu