1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065 |
- /* specfunc/hyperg_1F1.c
- *
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 3 of the License, or (at
- * your option) any later version.
- *
- * This program is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- */
- /* Author: G. Jungman */
- #include "gsl__config.h"
- #include "gsl_math.h"
- #include "gsl_errno.h"
- #include "gsl_sf_elementary.h"
- #include "gsl_sf_exp.h"
- #include "gsl_sf_bessel.h"
- #include "gsl_sf_gamma.h"
- #include "gsl_sf_laguerre.h"
- #include "gsl_sf_hyperg.h"
- #include "gsl_specfunc__error.h"
- #include "gsl_specfunc__hyperg.h"
- #define _1F1_INT_THRESHOLD (100.0*GSL_DBL_EPSILON)
- /* Asymptotic result for 1F1(a, b, x) x -> -Infinity.
- * Assumes b-a != neg integer and b != neg integer.
- */
- static
- int
- hyperg_1F1_asymp_negx(const double a, const double b, const double x,
- gsl_sf_result * result)
- {
- gsl_sf_result lg_b;
- gsl_sf_result lg_bma;
- double sgn_b;
- double sgn_bma;
- int stat_b = gsl_sf_lngamma_sgn_e(b, &lg_b, &sgn_b);
- int stat_bma = gsl_sf_lngamma_sgn_e(b-a, &lg_bma, &sgn_bma);
- if(stat_b == GSL_SUCCESS && stat_bma == GSL_SUCCESS) {
- gsl_sf_result F;
- int stat_F = gsl_sf_hyperg_2F0_series_e(a, 1.0+a-b, -1.0/x, -1, &F);
- if(F.val != 0) {
- double ln_term_val = a*log(-x);
- double ln_term_err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + fabs(ln_term_val));
- double ln_pre_val = lg_b.val - lg_bma.val - ln_term_val;
- double ln_pre_err = lg_b.err + lg_bma.err + ln_term_err;
- int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
- sgn_bma*sgn_b*F.val, F.err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_F);
- }
- else {
- result->val = 0.0;
- result->err = 0.0;
- return stat_F;
- }
- }
- else {
- DOMAIN_ERROR(result);
- }
- }
- /* Asymptotic result for 1F1(a, b, x) x -> +Infinity
- * Assumes b != neg integer and a != neg integer
- */
- static
- int
- hyperg_1F1_asymp_posx(const double a, const double b, const double x,
- gsl_sf_result * result)
- {
- gsl_sf_result lg_b;
- gsl_sf_result lg_a;
- double sgn_b;
- double sgn_a;
- int stat_b = gsl_sf_lngamma_sgn_e(b, &lg_b, &sgn_b);
- int stat_a = gsl_sf_lngamma_sgn_e(a, &lg_a, &sgn_a);
- if(stat_a == GSL_SUCCESS && stat_b == GSL_SUCCESS) {
- gsl_sf_result F;
- int stat_F = gsl_sf_hyperg_2F0_series_e(b-a, 1.0-a, 1.0/x, -1, &F);
- if(stat_F == GSL_SUCCESS && F.val != 0) {
- double lnx = log(x);
- double ln_term_val = (a-b)*lnx;
- double ln_term_err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + fabs(b)) * fabs(lnx)
- + 2.0 * GSL_DBL_EPSILON * fabs(a-b);
- double ln_pre_val = lg_b.val - lg_a.val + ln_term_val + x;
- double ln_pre_err = lg_b.err + lg_a.err + ln_term_err + 2.0 * GSL_DBL_EPSILON * fabs(x);
- int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
- sgn_a*sgn_b*F.val, F.err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_F);
- }
- else {
- result->val = 0.0;
- result->err = 0.0;
- return stat_F;
- }
- }
- else {
- DOMAIN_ERROR(result);
- }
- }
- /* Asymptotic result from Slater 4.3.7
- *
- * To get the general series, write M(a,b,x) as
- *
- * M(a,b,x)=sum ((a)_n/(b)_n) (x^n / n!)
- *
- * and expand (b)_n in inverse powers of b as follows
- *
- * -log(1/(b)_n) = sum_(k=0)^(n-1) log(b+k)
- * = n log(b) + sum_(k=0)^(n-1) log(1+k/b)
- *
- * Do a taylor expansion of the log in 1/b and sum the resulting terms
- * using the standard algebraic formulas for finite sums of powers of
- * k. This should then give
- *
- * M(a,b,x) = sum_(n=0)^(inf) (a_n/n!) (x/b)^n * (1 - n(n-1)/(2b)
- * + (n-1)n(n+1)(3n-2)/(24b^2) + ...
- *
- * which can be summed explicitly. The trick for summing it is to take
- * derivatives of sum_(i=0)^(inf) a_n*y^n/n! = (1-y)^(-a);
- *
- * [BJG 16/01/2007]
- */
- static
- int
- hyperg_1F1_largebx(const double a, const double b, const double x, gsl_sf_result * result)
- {
- double y = x/b;
- double f = exp(-a*log1p(-y));
- double t1 = -((a*(a+1.0))/(2*b))*pow((y/(1.0-y)),2.0);
- double t2 = (1/(24*b*b))*((a*(a+1)*y*y)/pow(1-y,4))*(12+8*(2*a+1)*y+(3*a*a-a-2)*y*y);
- double t3 = (-1/(48*b*b*b*pow(1-y,6)))*a*((a + 1)*((y*((a + 1)*(a*(y*(y*((y*(a - 2) + 16)*(a - 1)) + 72)) + 96)) + 24)*pow(y, 2)));
- result->val = f * (1 + t1 + t2 + t3);
- result->err = 2*fabs(f*t3) + 2*GSL_DBL_EPSILON*fabs(result->val);
- return GSL_SUCCESS;
- }
-
- /* Asymptotic result for x < 2b-4a, 2b-4a large.
- * [Abramowitz+Stegun, 13.5.21]
- *
- * assumes 0 <= x/(2b-4a) <= 1
- */
- static
- int
- hyperg_1F1_large2bm4a(const double a, const double b, const double x, gsl_sf_result * result)
- {
- double eta = 2.0*b - 4.0*a;
- double cos2th = x/eta;
- double sin2th = 1.0 - cos2th;
- double th = acos(sqrt(cos2th));
- double pre_h = 0.25*M_PI*M_PI*eta*eta*cos2th*sin2th;
- gsl_sf_result lg_b;
- int stat_lg = gsl_sf_lngamma_e(b, &lg_b);
- double t1 = 0.5*(1.0-b)*log(0.25*x*eta);
- double t2 = 0.25*log(pre_h);
- double lnpre_val = lg_b.val + 0.5*x + t1 - t2;
- double lnpre_err = lg_b.err + 2.0 * GSL_DBL_EPSILON * (fabs(0.5*x) + fabs(t1) + fabs(t2));
- #if SMALL_ANGLE
- const double eps = asin(sqrt(cos2th)); /* theta = pi/2 - eps */
- double s1 = (fmod(a, 1.0) == 0.0) ? 0.0 : sin(a*M_PI);
- double eta_reduc = (fmod(eta + 1, 4.0) == 0.0) ? 0.0 : fmod(eta + 1, 8.0);
- double phi1 = 0.25*eta_reduc*M_PI;
- double phi2 = 0.25*eta*(2*eps + sin(2.0*eps));
- double s2 = sin(phi1 - phi2);
- #else
- double s1 = sin(a*M_PI);
- double s2 = sin(0.25*eta*(2.0*th - sin(2.0*th)) + 0.25*M_PI);
- #endif
- double ser_val = s1 + s2;
- double ser_err = 2.0 * GSL_DBL_EPSILON * (fabs(s1) + fabs(s2));
- int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
- ser_val, ser_err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_lg);
- }
- /* Luke's rational approximation.
- * See [Luke, Algorithms for the Computation of Mathematical Functions, p.182]
- *
- * Like the case of the 2F1 rational approximations, these are
- * probably guaranteed to converge for x < 0, barring gross
- * numerical instability in the pre-asymptotic regime.
- */
- static
- int
- hyperg_1F1_luke(const double a, const double c, const double xin,
- gsl_sf_result * result)
- {
- const double RECUR_BIG = 1.0e+50;
- const int nmax = 5000;
- int n = 3;
- const double x = -xin;
- const double x3 = x*x*x;
- const double t0 = a/c;
- const double t1 = (a+1.0)/(2.0*c);
- const double t2 = (a+2.0)/(2.0*(c+1.0));
- double F = 1.0;
- double prec;
- double Bnm3 = 1.0; /* B0 */
- double Bnm2 = 1.0 + t1 * x; /* B1 */
- double Bnm1 = 1.0 + t2 * x * (1.0 + t1/3.0 * x); /* B2 */
-
- double Anm3 = 1.0; /* A0 */
- double Anm2 = Bnm2 - t0 * x; /* A1 */
- double Anm1 = Bnm1 - t0*(1.0 + t2*x)*x + t0 * t1 * (c/(c+1.0)) * x*x; /* A2 */
- while(1) {
- double npam1 = n + a - 1;
- double npcm1 = n + c - 1;
- double npam2 = n + a - 2;
- double npcm2 = n + c - 2;
- double tnm1 = 2*n - 1;
- double tnm3 = 2*n - 3;
- double tnm5 = 2*n - 5;
- double F1 = (n-a-2) / (2*tnm3*npcm1);
- double F2 = (n+a)*npam1 / (4*tnm1*tnm3*npcm2*npcm1);
- double F3 = -npam2*npam1*(n-a-2) / (8*tnm3*tnm3*tnm5*(n+c-3)*npcm2*npcm1);
- double E = -npam1*(n-c-1) / (2*tnm3*npcm2*npcm1);
- double An = (1.0+F1*x)*Anm1 + (E + F2*x)*x*Anm2 + F3*x3*Anm3;
- double Bn = (1.0+F1*x)*Bnm1 + (E + F2*x)*x*Bnm2 + F3*x3*Bnm3;
- double r = An/Bn;
- prec = fabs((F - r)/F);
- F = r;
- if(prec < GSL_DBL_EPSILON || n > nmax) break;
- if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
- An /= RECUR_BIG;
- Bn /= RECUR_BIG;
- Anm1 /= RECUR_BIG;
- Bnm1 /= RECUR_BIG;
- Anm2 /= RECUR_BIG;
- Bnm2 /= RECUR_BIG;
- Anm3 /= RECUR_BIG;
- Bnm3 /= RECUR_BIG;
- }
- else if(fabs(An) < 1.0/RECUR_BIG || fabs(Bn) < 1.0/RECUR_BIG) {
- An *= RECUR_BIG;
- Bn *= RECUR_BIG;
- Anm1 *= RECUR_BIG;
- Bnm1 *= RECUR_BIG;
- Anm2 *= RECUR_BIG;
- Bnm2 *= RECUR_BIG;
- Anm3 *= RECUR_BIG;
- Bnm3 *= RECUR_BIG;
- }
- n++;
- Bnm3 = Bnm2;
- Bnm2 = Bnm1;
- Bnm1 = Bn;
- Anm3 = Anm2;
- Anm2 = Anm1;
- Anm1 = An;
- }
- result->val = F;
- result->err = 2.0 * fabs(F * prec);
- result->err += 2.0 * GSL_DBL_EPSILON * (n-1.0) * fabs(F);
- return GSL_SUCCESS;
- }
- /* Series for 1F1(1,b,x)
- * b > 0
- */
- static
- int
- hyperg_1F1_1_series(const double b, const double x, gsl_sf_result * result)
- {
- double sum_val = 1.0;
- double sum_err = 0.0;
- double term = 1.0;
- double n = 1.0;
- while(fabs(term/sum_val) > 0.25*GSL_DBL_EPSILON) {
- term *= x/(b+n-1);
- sum_val += term;
- sum_err += 8.0*GSL_DBL_EPSILON*fabs(term) + GSL_DBL_EPSILON*fabs(sum_val);
- n += 1.0;
- }
- result->val = sum_val;
- result->err = sum_err;
- result->err += 2.0 * fabs(term);
- return GSL_SUCCESS;
- }
- /* 1F1(1,b,x)
- * b >= 1, b integer
- */
- static
- int
- hyperg_1F1_1_int(const int b, const double x, gsl_sf_result * result)
- {
- if(b < 1) {
- DOMAIN_ERROR(result);
- }
- else if(b == 1) {
- return gsl_sf_exp_e(x, result);
- }
- else if(b == 2) {
- return gsl_sf_exprel_e(x, result);
- }
- else if(b == 3) {
- return gsl_sf_exprel_2_e(x, result);
- }
- else {
- return gsl_sf_exprel_n_e(b-1, x, result);
- }
- }
- /* 1F1(1,b,x)
- * b >=1, b real
- *
- * checked OK: [GJ] Thu Oct 1 16:46:35 MDT 1998
- */
- static
- int
- hyperg_1F1_1(const double b, const double x, gsl_sf_result * result)
- {
- double ax = fabs(x);
- double ib = floor(b + 0.1);
- if(b < 1.0) {
- DOMAIN_ERROR(result);
- }
- else if(b == 1.0) {
- return gsl_sf_exp_e(x, result);
- }
- else if(b >= 1.4*ax) {
- return hyperg_1F1_1_series(b, x, result);
- }
- else if(fabs(b - ib) < _1F1_INT_THRESHOLD && ib < INT_MAX) {
- return hyperg_1F1_1_int((int)ib, x, result);
- }
- else if(x > 0.0) {
- if(x > 100.0 && b < 0.75*x) {
- return hyperg_1F1_asymp_posx(1.0, b, x, result);
- }
- else if(b < 1.0e+05) {
- /* Recurse backward on b, from a
- * chosen offset point. For x > 0,
- * which holds here, this should
- * be a stable direction.
- */
- const double off = ceil(1.4*x-b) + 1.0;
- double bp = b + off;
- gsl_sf_result M;
- int stat_s = hyperg_1F1_1_series(bp, x, &M);
- const double err_rat = M.err / fabs(M.val);
- while(bp > b+0.1) {
- /* M(1,b-1) = x/(b-1) M(1,b) + 1 */
- bp -= 1.0;
- M.val = 1.0 + x/bp * M.val;
- }
- result->val = M.val;
- result->err = err_rat * fabs(M.val);
- result->err += 2.0 * GSL_DBL_EPSILON * (fabs(off)+1.0) * fabs(M.val);
- return stat_s;
- } else if (fabs(x) < fabs(b) && fabs(x) < sqrt(fabs(b)) * fabs(b-x)) {
- return hyperg_1F1_largebx(1.0, b, x, result);
- } else if (fabs(x) > fabs(b)) {
- return hyperg_1F1_1_series(b, x, result);
- } else {
- return hyperg_1F1_large2bm4a(1.0, b, x, result);
- }
- }
- else {
- /* x <= 0 and b not large compared to |x|
- */
- if(ax < 10.0 && b < 10.0) {
- return hyperg_1F1_1_series(b, x, result);
- }
- else if(ax >= 100.0 && GSL_MAX_DBL(fabs(2.0-b),1.0) < 0.99*ax) {
- return hyperg_1F1_asymp_negx(1.0, b, x, result);
- }
- else {
- return hyperg_1F1_luke(1.0, b, x, result);
- }
- }
- }
- /* 1F1(a,b,x)/Gamma(b) for b->0
- * [limit of Abramowitz+Stegun 13.3.7]
- */
- static
- int
- hyperg_1F1_renorm_b0(const double a, const double x, gsl_sf_result * result)
- {
- double eta = a*x;
- if(eta > 0.0) {
- double root_eta = sqrt(eta);
- gsl_sf_result I1_scaled;
- int stat_I = gsl_sf_bessel_I1_scaled_e(2.0*root_eta, &I1_scaled);
- if(I1_scaled.val <= 0.0) {
- result->val = 0.0;
- result->err = 0.0;
- return GSL_ERROR_SELECT_2(stat_I, GSL_EDOM);
- }
- else {
- /* Note that 13.3.7 contains higher terms which are zeroth order
- in b. These make a non-negligible contribution to the sum.
- With the first correction term, the I1 above is replaced by
- I1 + (2/3)*a*(x/(4a))**(3/2)*I2(2*root_eta). We will add
- this as part of the result and error estimate. */
- const double corr1 =(2.0/3.0)*a*pow(x/(4.0*a),1.5)*gsl_sf_bessel_In_scaled(2, 2.0*root_eta)
- ;
- const double lnr_val = 0.5*x + 0.5*log(eta) + fabs(2.0*root_eta) + log(I1_scaled.val+corr1);
- const double lnr_err = GSL_DBL_EPSILON * (1.5*fabs(x) + 1.0) + fabs((I1_scaled.err+corr1)/I1_scaled.val);
- return gsl_sf_exp_err_e(lnr_val, lnr_err, result);
- }
- }
- else if(eta == 0.0) {
- result->val = 0.0;
- result->err = 0.0;
- return GSL_SUCCESS;
- }
- else {
- /* eta < 0 */
- double root_eta = sqrt(-eta);
- gsl_sf_result J1;
- int stat_J = gsl_sf_bessel_J1_e(2.0*root_eta, &J1);
- if(J1.val <= 0.0) {
- result->val = 0.0;
- result->err = 0.0;
- return GSL_ERROR_SELECT_2(stat_J, GSL_EDOM);
- }
- else {
- const double t1 = 0.5*x;
- const double t2 = 0.5*log(-eta);
- const double t3 = fabs(x);
- const double t4 = log(J1.val);
- const double lnr_val = t1 + t2 + t3 + t4;
- const double lnr_err = GSL_DBL_EPSILON * (1.5*fabs(x) + 1.0) + fabs(J1.err/J1.val);
- gsl_sf_result ex;
- int stat_e = gsl_sf_exp_err_e(lnr_val, lnr_err, &ex);
- result->val = -ex.val;
- result->err = ex.err;
- return stat_e;
- }
- }
-
- }
- /* 1F1'(a,b,x)/1F1(a,b,x)
- * Uses Gautschi's version of the CF.
- * [Gautschi, Math. Comp. 31, 994 (1977)]
- *
- * Supposedly this suffers from the "anomalous convergence"
- * problem when b < x. I have seen anomalous convergence
- * in several of the continued fractions associated with
- * 1F1(a,b,x). This particular CF formulation seems stable
- * for b > x. However, it does display a painful artifact
- * of the anomalous convergence; the convergence plateaus
- * unless b >>> x. For example, even for b=1000, x=1, this
- * method locks onto a ratio which is only good to about
- * 4 digits. Apparently the rest of the digits are hiding
- * way out on the plateau, but finite-precision lossage
- * means you will never get them.
- */
- #if 0
- static
- int
- hyperg_1F1_CF1_p(const double a, const double b, const double x, double * result)
- {
- const double RECUR_BIG = GSL_SQRT_DBL_MAX;
- const int maxiter = 5000;
- int n = 1;
- double Anm2 = 1.0;
- double Bnm2 = 0.0;
- double Anm1 = 0.0;
- double Bnm1 = 1.0;
- double a1 = 1.0;
- double b1 = 1.0;
- double An = b1*Anm1 + a1*Anm2;
- double Bn = b1*Bnm1 + a1*Bnm2;
- double an, bn;
- double fn = An/Bn;
- while(n < maxiter) {
- double old_fn;
- double del;
- n++;
- Anm2 = Anm1;
- Bnm2 = Bnm1;
- Anm1 = An;
- Bnm1 = Bn;
- an = (a+n)*x/((b-x+n-1)*(b-x+n));
- bn = 1.0;
- An = bn*Anm1 + an*Anm2;
- Bn = bn*Bnm1 + an*Bnm2;
- if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
- An /= RECUR_BIG;
- Bn /= RECUR_BIG;
- Anm1 /= RECUR_BIG;
- Bnm1 /= RECUR_BIG;
- Anm2 /= RECUR_BIG;
- Bnm2 /= RECUR_BIG;
- }
- old_fn = fn;
- fn = An/Bn;
- del = old_fn/fn;
-
- if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
- }
- *result = a/(b-x) * fn;
- if(n == maxiter)
- GSL_ERROR ("error", GSL_EMAXITER);
- else
- return GSL_SUCCESS;
- }
- #endif /* 0 */
- /* 1F1'(a,b,x)/1F1(a,b,x)
- * Uses Gautschi's series transformation of the
- * continued fraction. This is apparently the best
- * method for getting this ratio in the stable region.
- * The convergence is monotone and supergeometric
- * when b > x.
- * Assumes a >= -1.
- */
- static
- int
- hyperg_1F1_CF1_p_ser(const double a, const double b, const double x, double * result)
- {
- if(a == 0.0) {
- *result = 0.0;
- return GSL_SUCCESS;
- }
- else {
- const int maxiter = 5000;
- double sum = 1.0;
- double pk = 1.0;
- double rhok = 0.0;
- int k;
- for(k=1; k<maxiter; k++) {
- double ak = (a + k)*x/((b-x+k-1.0)*(b-x+k));
- rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0+rhok));
- pk *= rhok;
- sum += pk;
- if(fabs(pk/sum) < 2.0*GSL_DBL_EPSILON) break;
- }
- *result = a/(b-x) * sum;
- if(k == maxiter)
- GSL_ERROR ("error", GSL_EMAXITER);
- else
- return GSL_SUCCESS;
- }
- }
- /* 1F1(a+1,b,x)/1F1(a,b,x)
- *
- * I think this suffers from typical "anomalous convergence".
- * I could not find a region where it was truly useful.
- */
- #if 0
- static
- int
- hyperg_1F1_CF1(const double a, const double b, const double x, double * result)
- {
- const double RECUR_BIG = GSL_SQRT_DBL_MAX;
- const int maxiter = 5000;
- int n = 1;
- double Anm2 = 1.0;
- double Bnm2 = 0.0;
- double Anm1 = 0.0;
- double Bnm1 = 1.0;
- double a1 = b - a - 1.0;
- double b1 = b - x - 2.0*(a+1.0);
- double An = b1*Anm1 + a1*Anm2;
- double Bn = b1*Bnm1 + a1*Bnm2;
- double an, bn;
- double fn = An/Bn;
- while(n < maxiter) {
- double old_fn;
- double del;
- n++;
- Anm2 = Anm1;
- Bnm2 = Bnm1;
- Anm1 = An;
- Bnm1 = Bn;
- an = (a + n - 1.0) * (b - a - n);
- bn = b - x - 2.0*(a+n);
- An = bn*Anm1 + an*Anm2;
- Bn = bn*Bnm1 + an*Bnm2;
- if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
- An /= RECUR_BIG;
- Bn /= RECUR_BIG;
- Anm1 /= RECUR_BIG;
- Bnm1 /= RECUR_BIG;
- Anm2 /= RECUR_BIG;
- Bnm2 /= RECUR_BIG;
- }
- old_fn = fn;
- fn = An/Bn;
- del = old_fn/fn;
-
- if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
- }
- *result = fn;
- if(n == maxiter)
- GSL_ERROR ("error", GSL_EMAXITER);
- else
- return GSL_SUCCESS;
- }
- #endif /* 0 */
- /* 1F1(a,b+1,x)/1F1(a,b,x)
- *
- * This seemed to suffer from "anomalous convergence".
- * However, I have no theory for this recurrence.
- */
- #if 0
- static
- int
- hyperg_1F1_CF1_b(const double a, const double b, const double x, double * result)
- {
- const double RECUR_BIG = GSL_SQRT_DBL_MAX;
- const int maxiter = 5000;
- int n = 1;
- double Anm2 = 1.0;
- double Bnm2 = 0.0;
- double Anm1 = 0.0;
- double Bnm1 = 1.0;
- double a1 = b + 1.0;
- double b1 = (b + 1.0) * (b - x);
- double An = b1*Anm1 + a1*Anm2;
- double Bn = b1*Bnm1 + a1*Bnm2;
- double an, bn;
- double fn = An/Bn;
- while(n < maxiter) {
- double old_fn;
- double del;
- n++;
- Anm2 = Anm1;
- Bnm2 = Bnm1;
- Anm1 = An;
- Bnm1 = Bn;
- an = (b + n) * (b + n - 1.0 - a) * x;
- bn = (b + n) * (b + n - 1.0 - x);
- An = bn*Anm1 + an*Anm2;
- Bn = bn*Bnm1 + an*Bnm2;
- if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
- An /= RECUR_BIG;
- Bn /= RECUR_BIG;
- Anm1 /= RECUR_BIG;
- Bnm1 /= RECUR_BIG;
- Anm2 /= RECUR_BIG;
- Bnm2 /= RECUR_BIG;
- }
- old_fn = fn;
- fn = An/Bn;
- del = old_fn/fn;
-
- if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
- }
- *result = fn;
- if(n == maxiter)
- GSL_ERROR ("error", GSL_EMAXITER);
- else
- return GSL_SUCCESS;
- }
- #endif /* 0 */
- /* 1F1(a,b,x)
- * |a| <= 1, b > 0
- */
- static
- int
- hyperg_1F1_small_a_bgt0(const double a, const double b, const double x, gsl_sf_result * result)
- {
- const double bma = b-a;
- const double oma = 1.0-a;
- const double ap1mb = 1.0+a-b;
- const double abs_bma = fabs(bma);
- const double abs_oma = fabs(oma);
- const double abs_ap1mb = fabs(ap1mb);
- const double ax = fabs(x);
- if(a == 0.0) {
- result->val = 1.0;
- result->err = 0.0;
- return GSL_SUCCESS;
- }
- else if(a == 1.0 && b >= 1.0) {
- return hyperg_1F1_1(b, x, result);
- }
- else if(a == -1.0) {
- result->val = 1.0 + a/b * x;
- result->err = GSL_DBL_EPSILON * (1.0 + fabs(a/b * x));
- result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
- return GSL_SUCCESS;
- }
- else if(b >= 1.4*ax) {
- return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
- }
- else if(x > 0.0) {
- if(x > 100.0 && abs_bma*abs_oma < 0.5*x) {
- return hyperg_1F1_asymp_posx(a, b, x, result);
- }
- else if(b < 5.0e+06) {
- /* Recurse backward on b from
- * a suitably high point.
- */
- const double b_del = ceil(1.4*x-b) + 1.0;
- double bp = b + b_del;
- gsl_sf_result r_Mbp1;
- gsl_sf_result r_Mb;
- double Mbp1;
- double Mb;
- double Mbm1;
- int stat_0 = gsl_sf_hyperg_1F1_series_e(a, bp+1.0, x, &r_Mbp1);
- int stat_1 = gsl_sf_hyperg_1F1_series_e(a, bp, x, &r_Mb);
- const double err_rat = fabs(r_Mbp1.err/r_Mbp1.val) + fabs(r_Mb.err/r_Mb.val);
- Mbp1 = r_Mbp1.val;
- Mb = r_Mb.val;
- while(bp > b+0.1) {
- /* Do backward recursion. */
- Mbm1 = ((x+bp-1.0)*Mb - x*(bp-a)/bp*Mbp1)/(bp-1.0);
- bp -= 1.0;
- Mbp1 = Mb;
- Mb = Mbm1;
- }
- result->val = Mb;
- result->err = err_rat * (fabs(b_del)+1.0) * fabs(Mb);
- result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mb);
- return GSL_ERROR_SELECT_2(stat_0, stat_1);
- }
- else if (fabs(x) < fabs(b) && fabs(a*x) < sqrt(fabs(b)) * fabs(b-x)) {
- return hyperg_1F1_largebx(a, b, x, result);
- } else {
- return hyperg_1F1_large2bm4a(a, b, x, result);
- }
- }
- else {
- /* x < 0 and b not large compared to |x|
- */
- if(ax < 10.0 && b < 10.0) {
- return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
- }
- else if(ax >= 100.0 && GSL_MAX(abs_ap1mb,1.0) < 0.99*ax) {
- return hyperg_1F1_asymp_negx(a, b, x, result);
- }
- else {
- return hyperg_1F1_luke(a, b, x, result);
- }
- }
- }
- /* 1F1(b+eps,b,x)
- * |eps|<=1, b > 0
- */
- static
- int
- hyperg_1F1_beps_bgt0(const double eps, const double b, const double x, gsl_sf_result * result)
- {
- if(b > fabs(x) && fabs(eps) < GSL_SQRT_DBL_EPSILON) {
- /* If b-a is very small and x/b is not too large we can
- * use this explicit approximation.
- *
- * 1F1(b+eps,b,x) = exp(ax/b) (1 - eps x^2 (v2 + v3 x + ...) + ...)
- *
- * v2 = a/(2b^2(b+1))
- * v3 = a(b-2a)/(3b^3(b+1)(b+2))
- * ...
- *
- * See [Luke, Mathematical Functions and Their Approximations, p.292]
- *
- * This cannot be used for b near a negative integer or zero.
- * Also, if x/b is large the deviation from exp(x) behaviour grows.
- */
- double a = b + eps;
- gsl_sf_result exab;
- int stat_e = gsl_sf_exp_e(a*x/b, &exab);
- double v2 = a/(2.0*b*b*(b+1.0));
- double v3 = a*(b-2.0*a)/(3.0*b*b*b*(b+1.0)*(b+2.0));
- double v = v2 + v3 * x;
- double f = (1.0 - eps*x*x*v);
- result->val = exab.val * f;
- result->err = exab.err * fabs(f);
- result->err += fabs(exab.val) * GSL_DBL_EPSILON * (1.0 + fabs(eps*x*x*v));
- result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
- return stat_e;
- }
- else {
- /* Otherwise use a Kummer transformation to reduce
- * it to the small a case.
- */
- gsl_sf_result Kummer_1F1;
- int stat_K = hyperg_1F1_small_a_bgt0(-eps, b, -x, &Kummer_1F1);
- if(Kummer_1F1.val != 0.0) {
- int stat_e = gsl_sf_exp_mult_err_e(x, 2.0*GSL_DBL_EPSILON*fabs(x),
- Kummer_1F1.val, Kummer_1F1.err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_K);
- }
- else {
- result->val = 0.0;
- result->err = 0.0;
- return stat_K;
- }
- }
- }
- /* 1F1(a,2a,x) = Gamma(a + 1/2) E(x) (|x|/4)^(-a+1/2) scaled_I(a-1/2,|x|/2)
- *
- * E(x) = exp(x) x > 0
- * = 1 x < 0
- *
- * a >= 1/2
- */
- static
- int
- hyperg_1F1_beq2a_pos(const double a, const double x, gsl_sf_result * result)
- {
- if(x == 0.0) {
- result->val = 1.0;
- result->err = 0.0;
- return GSL_SUCCESS;
- }
- else {
- gsl_sf_result I;
- int stat_I = gsl_sf_bessel_Inu_scaled_e(a-0.5, 0.5*fabs(x), &I);
- gsl_sf_result lg;
- int stat_g = gsl_sf_lngamma_e(a + 0.5, &lg);
- double ln_term = (0.5-a)*log(0.25*fabs(x));
- double lnpre_val = lg.val + GSL_MAX_DBL(x,0.0) + ln_term;
- double lnpre_err = lg.err + GSL_DBL_EPSILON * (fabs(ln_term) + fabs(x));
- int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
- I.val, I.err,
- result);
- return GSL_ERROR_SELECT_3(stat_e, stat_g, stat_I);
- }
- }
- /* Determine middle parts of diagonal recursion along b=2a
- * from two endpoints, i.e.
- *
- * given: M(a,b) and M(a+1,b+2)
- * get: M(a+1,b+1) and M(a,b+1)
- */
- #if 0
- inline
- static
- int
- hyperg_1F1_diag_step(const double a, const double b, const double x,
- const double Mab, const double Map1bp2,
- double * Map1bp1, double * Mabp1)
- {
- if(a == b) {
- *Map1bp1 = Mab;
- *Mabp1 = Mab - x/(b+1.0) * Map1bp2;
- }
- else {
- *Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2;
- *Mabp1 = (a * *Map1bp1 - b * Mab)/(a-b);
- }
- return GSL_SUCCESS;
- }
- #endif /* 0 */
- /* Determine endpoint of diagonal recursion.
- *
- * given: M(a,b) and M(a+1,b+2)
- * get: M(a+1,b) and M(a+1,b+1)
- */
- #if 0
- inline
- static
- int
- hyperg_1F1_diag_end_step(const double a, const double b, const double x,
- const double Mab, const double Map1bp2,
- double * Map1b, double * Map1bp1)
- {
- *Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2;
- *Map1b = Mab + x/b * *Map1bp1;
- return GSL_SUCCESS;
- }
- #endif /* 0 */
- /* Handle the case of a and b both positive integers.
- * Assumes a > 0 and b > 0.
- */
- static
- int
- hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * result)
- {
- double ax = fabs(x);
- if(a == b) {
- return gsl_sf_exp_e(x, result); /* 1F1(a,a,x) */
- }
- else if(a == 1) {
- return gsl_sf_exprel_n_e(b-1, x, result); /* 1F1(1,b,x) */
- }
- else if(b == a + 1) {
- gsl_sf_result K;
- int stat_K = gsl_sf_exprel_n_e(a, -x, &K); /* 1F1(1,1+a,-x) */
- int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x),
- K.val, K.err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_K);
- }
- else if(a == b + 1) {
- gsl_sf_result ex;
- int stat_e = gsl_sf_exp_e(x, &ex);
- result->val = ex.val * (1.0 + x/b);
- result->err = ex.err * (1.0 + x/b);
- result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b));
- result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
- return stat_e;
- }
- else if(a == b + 2) {
- gsl_sf_result ex;
- int stat_e = gsl_sf_exp_e(x, &ex);
- double poly = (1.0 + x/b*(2.0 + x/(b+1.0)));
- result->val = ex.val * poly;
- result->err = ex.err * fabs(poly);
- result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b) * (2.0 + fabs(x/(b+1.0))));
- result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
- return stat_e;
- }
- else if(b == 2*a) {
- return hyperg_1F1_beq2a_pos(a, x, result); /* 1F1(a,2a,x) */
- }
- else if( ( b < 10 && a < 10 && ax < 5.0 )
- || ( b > a*ax )
- || ( b > a && ax < 5.0 )
- ) {
- return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
- }
- else if(b > a && b >= 2*a + x) {
- /* Use the Gautschi CF series, then
- * recurse backward to a=0 for normalization.
- * This will work for either sign of x.
- */
- double rap;
- int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
- double ra = 1.0 + x/a * rap;
- double Ma = GSL_SQRT_DBL_MIN;
- double Map1 = ra * Ma;
- double Mnp1 = Map1;
- double Mn = Ma;
- double Mnm1;
- int n;
- for(n=a; n>0; n--) {
- Mnm1 = (n * Mnp1 - (2*n-b+x) * Mn) / (b-n);
- Mnp1 = Mn;
- Mn = Mnm1;
- }
- result->val = Ma/Mn;
- result->err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + 1.0) * fabs(Ma/Mn);
- return stat_CF1;
- }
- else if(b > a && b < 2*a + x && b > x) {
- /* Use the Gautschi series representation of
- * the continued fraction. Then recurse forward
- * to the a=b line for normalization. This will
- * work for either sign of x, although we do need
- * to check for b > x, for when x is positive.
- */
- double rap;
- int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
- double ra = 1.0 + x/a * rap;
- gsl_sf_result ex;
- int stat_ex;
- double Ma = GSL_SQRT_DBL_MIN;
- double Map1 = ra * Ma;
- double Mnm1 = Ma;
- double Mn = Map1;
- double Mnp1;
- int n;
- for(n=a+1; n<b; n++) {
- Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
- Mnm1 = Mn;
- Mn = Mnp1;
- }
- stat_ex = gsl_sf_exp_e(x, &ex); /* 1F1(b,b,x) */
- result->val = ex.val * Ma/Mn;
- result->err = ex.err * fabs(Ma/Mn);
- result->err += 4.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val);
- return GSL_ERROR_SELECT_2(stat_ex, stat_CF1);
- }
- else if(x >= 0.0) {
- if(b < a) {
- /* The point b,b is below the b=2a+x line.
- * Forward recursion on a from b,b+1 is possible.
- * Note that a > b + 1 as well, since we already tried a = b + 1.
- */
- if(x + log(fabs(x/b)) < GSL_LOG_DBL_MAX-2.0) {
- double ex = exp(x);
- int n;
- double Mnm1 = ex; /* 1F1(b,b,x) */
- double Mn = ex * (1.0 + x/b); /* 1F1(b+1,b,x) */
- double Mnp1;
- for(n=b+1; n<a; n++) {
- Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
- Mnm1 = Mn;
- Mn = Mnp1;
- }
- result->val = Mn;
- result->err = (x + 1.0) * GSL_DBL_EPSILON * fabs(Mn);
- result->err *= fabs(a-b)+1.0;
- return GSL_SUCCESS;
- }
- else {
- OVERFLOW_ERROR(result);
- }
- }
- else {
- /* b > a
- * b < 2a + x
- * b <= x (otherwise we would have finished above)
- *
- * Gautschi anomalous convergence region. However, we can
- * recurse forward all the way from a=0,1 because we are
- * always underneath the b=2a+x line.
- */
- gsl_sf_result r_Mn;
- double Mnm1 = 1.0; /* 1F1(0,b,x) */
- double Mn; /* 1F1(1,b,x) */
- double Mnp1;
- int n;
- gsl_sf_exprel_n_e(b-1, x, &r_Mn);
- Mn = r_Mn.val;
- for(n=1; n<a; n++) {
- Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
- Mnm1 = Mn;
- Mn = Mnp1;
- }
- result->val = Mn;
- result->err = fabs(Mn) * (1.0 + fabs(a)) * fabs(r_Mn.err / r_Mn.val);
- result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
- return GSL_SUCCESS;
- }
- }
- else {
- /* x < 0
- * b < a (otherwise we would have tripped one of the above)
- */
- if(a <= 0.5*(b-x) || a >= -x) {
- /* Gautschi continued fraction is in the anomalous region,
- * so we must find another way. We recurse down in b,
- * from the a=b line.
- */
- double ex = exp(x);
- double Manp1 = ex;
- double Man = ex * (1.0 + x/(a-1.0));
- double Manm1;
- int n;
- for(n=a-1; n>b; n--) {
- Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0));
- Manp1 = Man;
- Man = Manm1;
- }
- result->val = Man;
- result->err = (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(Man);
- result->err *= fabs(b-a)+1.0;
- return GSL_SUCCESS;
- }
- else {
- /* Pick a0 such that b ~= 2a0 + x, then
- * recurse down in b from a0,a0 to determine
- * the values near the line b=2a+x. Then recurse
- * forward on a from a0.
- */
- int a0 = ceil(0.5*(b-x));
- double Ma0b; /* M(a0,b) */
- double Ma0bp1; /* M(a0,b+1) */
- double Ma0p1b; /* M(a0+1,b) */
- double Mnm1;
- double Mn;
- double Mnp1;
- int n;
- {
- double ex = exp(x);
- double Ma0np1 = ex;
- double Ma0n = ex * (1.0 + x/(a0-1.0));
- double Ma0nm1;
- for(n=a0-1; n>b; n--) {
- Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0));
- Ma0np1 = Ma0n;
- Ma0n = Ma0nm1;
- }
- Ma0bp1 = Ma0np1;
- Ma0b = Ma0n;
- Ma0p1b = (b*(a0+x)*Ma0b + x*(a0-b)*Ma0bp1)/(a0*b);
- }
- /* Initialise the recurrence correctly BJG */
- if (a0 >= a)
- {
- Mn = Ma0b;
- }
- else if (a0 + 1>= a)
- {
- Mn = Ma0p1b;
- }
- else
- {
- Mnm1 = Ma0b;
- Mn = Ma0p1b;
- for(n=a0+1; n<a; n++) {
- Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
- Mnm1 = Mn;
- Mn = Mnp1;
- }
- }
- result->val = Mn;
- result->err = (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(Mn);
- result->err *= fabs(b-a)+1.0;
- return GSL_SUCCESS;
- }
- }
- }
- /* Evaluate a <= 0, a integer, cases directly. (Polynomial; Horner)
- * When the terms are all positive, this
- * must work. We will assume this here.
- */
- static
- int
- hyperg_1F1_a_negint_poly(const int a, const double b, const double x, gsl_sf_result * result)
- {
- if(a == 0) {
- result->val = 1.0;
- result->err = 1.0;
- return GSL_SUCCESS;
- }
- else {
- int N = -a;
- double poly = 1.0;
- int k;
- for(k=N-1; k>=0; k--) {
- double t = (a+k)/(b+k) * (x/(k+1));
- double r = t + 1.0/poly;
- if(r > 0.9*GSL_DBL_MAX/poly) {
- OVERFLOW_ERROR(result);
- }
- else {
- poly *= r; /* P_n = 1 + t_n P_{n-1} */
- }
- }
- result->val = poly;
- result->err = 2.0 * (sqrt(N) + 1.0) * GSL_DBL_EPSILON * fabs(poly);
- return GSL_SUCCESS;
- }
- }
- /* Evaluate negative integer a case by relation
- * to Laguerre polynomials. This is more general than
- * the direct polynomial evaluation, but is safe
- * for all values of x.
- *
- * 1F1(-n,b,x) = n!/(b)_n Laguerre[n,b-1,x]
- * = n B(b,n) Laguerre[n,b-1,x]
- *
- * assumes b is not a negative integer
- */
- static
- int
- hyperg_1F1_a_negint_lag(const int a, const double b, const double x, gsl_sf_result * result)
- {
- const int n = -a;
- gsl_sf_result lag;
- const int stat_l = gsl_sf_laguerre_n_e(n, b-1.0, x, &lag);
- if(b < 0.0) {
- gsl_sf_result lnfact;
- gsl_sf_result lng1;
- gsl_sf_result lng2;
- double s1, s2;
- const int stat_f = gsl_sf_lnfact_e(n, &lnfact);
- const int stat_g1 = gsl_sf_lngamma_sgn_e(b + n, &lng1, &s1);
- const int stat_g2 = gsl_sf_lngamma_sgn_e(b, &lng2, &s2);
- const double lnpre_val = lnfact.val - (lng1.val - lng2.val);
- const double lnpre_err = lnfact.err + lng1.err + lng2.err
- + 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val);
- const int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
- s1*s2*lag.val, lag.err,
- result);
- return GSL_ERROR_SELECT_5(stat_e, stat_l, stat_g1, stat_g2, stat_f);
- }
- else {
- gsl_sf_result lnbeta;
- gsl_sf_lnbeta_e(b, n, &lnbeta);
- if(fabs(lnbeta.val) < 0.1) {
- /* As we have noted, when B(x,y) is near 1,
- * evaluating log(B(x,y)) is not accurate.
- * Instead we evaluate B(x,y) directly.
- */
- const double ln_term_val = log(1.25*n);
- const double ln_term_err = 2.0 * GSL_DBL_EPSILON * ln_term_val;
- gsl_sf_result beta;
- int stat_b = gsl_sf_beta_e(b, n, &beta);
- int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err,
- lag.val, lag.err,
- result);
- result->val *= beta.val/1.25;
- result->err *= beta.val/1.25;
- return GSL_ERROR_SELECT_3(stat_e, stat_l, stat_b);
- }
- else {
- /* B(x,y) was not near 1, so it is safe to use
- * the logarithmic values.
- */
- const double ln_n = log(n);
- const double ln_term_val = lnbeta.val + ln_n;
- const double ln_term_err = lnbeta.err + 2.0 * GSL_DBL_EPSILON * fabs(ln_n);
- int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err,
- lag.val, lag.err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_l);
- }
- }
- }
- /* Handle negative integer a case for x > 0 and
- * generic b.
- *
- * Combine [Abramowitz+Stegun, 13.6.9 + 13.6.27]
- * M(-n,b,x) = (-1)^n / (b)_n U(-n,b,x) = n! / (b)_n Laguerre^(b-1)_n(x)
- */
- #if 0
- static
- int
- hyperg_1F1_a_negint_U(const int a, const double b, const double x, gsl_sf_result * result)
- {
- const int n = -a;
- const double sgn = ( GSL_IS_ODD(n) ? -1.0 : 1.0 );
- double sgpoch;
- gsl_sf_result lnpoch;
- gsl_sf_result U;
- const int stat_p = gsl_sf_lnpoch_sgn_e(b, n, &lnpoch, &sgpoch);
- const int stat_U = gsl_sf_hyperg_U_e(-n, b, x, &U);
- const int stat_e = gsl_sf_exp_mult_err_e(-lnpoch.val, lnpoch.err,
- sgn * sgpoch * U.val, U.err,
- result);
- return GSL_ERROR_SELECT_3(stat_e, stat_U, stat_p);
- }
- #endif
- /* Assumes a <= -1, b <= -1, and b <= a.
- */
- static
- int
- hyperg_1F1_ab_negint(const int a, const int b, const double x, gsl_sf_result * result)
- {
- if(x == 0.0) {
- result->val = 1.0;
- result->err = 0.0;
- return GSL_SUCCESS;
- }
- else if(x > 0.0) {
- return hyperg_1F1_a_negint_poly(a, b, x, result);
- }
- else {
- /* Apply a Kummer transformation to make x > 0 so
- * we can evaluate the polynomial safely. Of course,
- * this assumes b <= a, which must be true for
- * a<0 and b<0, since otherwise the thing is undefined.
- */
- gsl_sf_result K;
- int stat_K = hyperg_1F1_a_negint_poly(b-a, b, -x, &K);
- int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x),
- K.val, K.err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_K);
- }
- }
- /* [Abramowitz+Stegun, 13.1.3]
- *
- * M(a,b,x) = Gamma(1+a-b)/Gamma(2-b) x^(1-b) *
- * { Gamma(b)/Gamma(a) M(1+a-b,2-b,x) - (b-1) U(1+a-b,2-b,x) }
- *
- * b not an integer >= 2
- * a-b not a negative integer
- */
- static
- int
- hyperg_1F1_U(const double a, const double b, const double x, gsl_sf_result * result)
- {
- const double bp = 2.0 - b;
- const double ap = a - b + 1.0;
- gsl_sf_result lg_ap, lg_bp;
- double sg_ap;
- int stat_lg0 = gsl_sf_lngamma_sgn_e(ap, &lg_ap, &sg_ap);
- int stat_lg1 = gsl_sf_lngamma_e(bp, &lg_bp);
- int stat_lg2 = GSL_ERROR_SELECT_2(stat_lg0, stat_lg1);
- double t1 = (bp-1.0) * log(x);
- double lnpre_val = lg_ap.val - lg_bp.val + t1;
- double lnpre_err = lg_ap.err + lg_bp.err + 2.0 * GSL_DBL_EPSILON * fabs(t1);
- gsl_sf_result lg_2mbp, lg_1papmbp;
- double sg_2mbp, sg_1papmbp;
- int stat_lg3 = gsl_sf_lngamma_sgn_e(2.0-bp, &lg_2mbp, &sg_2mbp);
- int stat_lg4 = gsl_sf_lngamma_sgn_e(1.0+ap-bp, &lg_1papmbp, &sg_1papmbp);
- int stat_lg5 = GSL_ERROR_SELECT_2(stat_lg3, stat_lg4);
- double lnc1_val = lg_2mbp.val - lg_1papmbp.val;
- double lnc1_err = lg_2mbp.err + lg_1papmbp.err
- + GSL_DBL_EPSILON * (fabs(lg_2mbp.val) + fabs(lg_1papmbp.val));
- gsl_sf_result M;
- gsl_sf_result_e10 U;
- int stat_F = gsl_sf_hyperg_1F1_e(ap, bp, x, &M);
- int stat_U = gsl_sf_hyperg_U_e10_e(ap, bp, x, &U);
- int stat_FU = GSL_ERROR_SELECT_2(stat_F, stat_U);
- gsl_sf_result_e10 term_M;
- int stat_e0 = gsl_sf_exp_mult_err_e10_e(lnc1_val, lnc1_err,
- sg_2mbp*sg_1papmbp*M.val, M.err,
- &term_M);
- const double ombp = 1.0 - bp;
- const double Uee_val = U.e10*M_LN10;
- const double Uee_err = 2.0 * GSL_DBL_EPSILON * fabs(Uee_val);
- const double Mee_val = term_M.e10*M_LN10;
- const double Mee_err = 2.0 * GSL_DBL_EPSILON * fabs(Mee_val);
- int stat_e1;
- /* Do a little dance with the exponential prefactors
- * to avoid overflows in intermediate results.
- */
- if(Uee_val > Mee_val) {
- const double factorM_val = exp(Mee_val-Uee_val);
- const double factorM_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorM_val;
- const double inner_val = term_M.val*factorM_val - ombp*U.val;
- const double inner_err =
- term_M.err*factorM_val + fabs(ombp) * U.err
- + fabs(term_M.val) * factorM_err
- + GSL_DBL_EPSILON * (fabs(term_M.val*factorM_val) + fabs(ombp*U.val));
- stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Uee_val, lnpre_err+Uee_err,
- sg_ap*inner_val, inner_err,
- result);
- }
- else {
- const double factorU_val = exp(Uee_val - Mee_val);
- const double factorU_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorU_val;
- const double inner_val = term_M.val - ombp*factorU_val*U.val;
- const double inner_err =
- term_M.err + fabs(ombp*factorU_val*U.err)
- + fabs(ombp*factorU_err*U.val)
- + GSL_DBL_EPSILON * (fabs(term_M.val) + fabs(ombp*factorU_val*U.val));
- stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Mee_val, lnpre_err+Mee_err,
- sg_ap*inner_val, inner_err,
- result);
- }
- return GSL_ERROR_SELECT_5(stat_e1, stat_e0, stat_FU, stat_lg5, stat_lg2);
- }
- /* Handle case of generic positive a, b.
- * Assumes b-a is not a negative integer.
- */
- static
- int
- hyperg_1F1_ab_pos(const double a, const double b,
- const double x,
- gsl_sf_result * result)
- {
- const double ax = fabs(x);
- if( ( b < 10.0 && a < 10.0 && ax < 5.0 )
- || ( b > a*ax )
- || ( b > a && ax < 5.0 )
- ) {
- return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
- }
- else if( x < -100.0
- && GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.7*fabs(x)
- ) {
- /* Large negative x asymptotic.
- */
- return hyperg_1F1_asymp_negx(a, b, x, result);
- }
- else if( x > 100.0
- && GSL_MAX_DBL(fabs(b-a),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.7*fabs(x)
- ) {
- /* Large positive x asymptotic.
- */
- return hyperg_1F1_asymp_posx(a, b, x, result);
- }
- else if(fabs(b-a) <= 1.0) {
- /* Directly handle b near a.
- */
- return hyperg_1F1_beps_bgt0(a-b, b, x, result); /* a = b + eps */
- }
- else if(b > a && b >= 2*a + x) {
- /* Use the Gautschi CF series, then
- * recurse backward to a near 0 for normalization.
- * This will work for either sign of x.
- */
- double rap;
- int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
- double ra = 1.0 + x/a * rap;
- double Ma = GSL_SQRT_DBL_MIN;
- double Map1 = ra * Ma;
- double Mnp1 = Map1;
- double Mn = Ma;
- double Mnm1;
- gsl_sf_result Mn_true;
- int stat_Mt;
- double n;
- for(n=a; n>0.5; n -= 1.0) {
- Mnm1 = (n * Mnp1 - (2.0*n-b+x) * Mn) / (b-n);
- Mnp1 = Mn;
- Mn = Mnm1;
- }
- stat_Mt = hyperg_1F1_small_a_bgt0(n, b, x, &Mn_true);
- result->val = (Ma/Mn) * Mn_true.val;
- result->err = fabs(Ma/Mn) * Mn_true.err;
- result->err += 2.0 * GSL_DBL_EPSILON * (fabs(a)+1.0) * fabs(result->val);
- return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1);
- }
- else if(b > a && b < 2*a + x && b > x) {
- /* Use the Gautschi series representation of
- * the continued fraction. Then recurse forward
- * to near the a=b line for normalization. This will
- * work for either sign of x, although we do need
- * to check for b > x, which is relevant when x is positive.
- */
- gsl_sf_result Mn_true;
- int stat_Mt;
- double rap;
- int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
- double ra = 1.0 + x/a * rap;
- double Ma = GSL_SQRT_DBL_MIN;
- double Mnm1 = Ma;
- double Mn = ra * Mnm1;
- double Mnp1;
- double n;
- for(n=a+1.0; n<b-0.5; n += 1.0) {
- Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
- Mnm1 = Mn;
- Mn = Mnp1;
- }
- stat_Mt = hyperg_1F1_beps_bgt0(n-b, b, x, &Mn_true);
- result->val = Ma/Mn * Mn_true.val;
- result->err = fabs(Ma/Mn) * Mn_true.err;
- result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val);
- return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1);
- }
- else if(x >= 0.0) {
- if(b < a) {
- /* Forward recursion on a from a=b+eps-1,b+eps.
- */
- double N = floor(a-b);
- double eps = a - b - N;
- gsl_sf_result r_M0;
- gsl_sf_result r_M1;
- int stat_0 = hyperg_1F1_beps_bgt0(eps-1.0, b, x, &r_M0);
- int stat_1 = hyperg_1F1_beps_bgt0(eps, b, x, &r_M1);
- double M0 = r_M0.val;
- double M1 = r_M1.val;
- double Mam1 = M0;
- double Ma = M1;
- double Map1;
- double ap;
- double start_pair = fabs(M0) + fabs(M1);
- double minim_pair = GSL_DBL_MAX;
- double pair_ratio;
- double rat_0 = fabs(r_M0.err/r_M0.val);
- double rat_1 = fabs(r_M1.err/r_M1.val);
- for(ap=b+eps; ap<a-0.1; ap += 1.0) {
- Map1 = ((b-ap)*Mam1 + (2.0*ap-b+x)*Ma)/ap;
- Mam1 = Ma;
- Ma = Map1;
- minim_pair = GSL_MIN_DBL(fabs(Mam1) + fabs(Ma), minim_pair);
- }
- pair_ratio = start_pair/minim_pair;
- result->val = Ma;
- result->err = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Ma);
- result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Ma);
- result->err += 2.0 * GSL_DBL_EPSILON * fabs(Ma);
- return GSL_ERROR_SELECT_2(stat_0, stat_1);
- }
- else {
- /* b > a
- * b < 2a + x
- * b <= x
- *
- * Recurse forward on a from a=eps,eps+1.
- */
- double eps = a - floor(a);
- gsl_sf_result r_Mnm1;
- gsl_sf_result r_Mn;
- int stat_0 = hyperg_1F1_small_a_bgt0(eps, b, x, &r_Mnm1);
- int stat_1 = hyperg_1F1_small_a_bgt0(eps+1.0, b, x, &r_Mn);
- double Mnm1 = r_Mnm1.val;
- double Mn = r_Mn.val;
- double Mnp1;
- double n;
- double start_pair = fabs(Mn) + fabs(Mnm1);
- double minim_pair = GSL_DBL_MAX;
- double pair_ratio;
- double rat_0 = fabs(r_Mnm1.err/r_Mnm1.val);
- double rat_1 = fabs(r_Mn.err/r_Mn.val);
- for(n=eps+1.0; n<a-0.1; n++) {
- Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
- Mnm1 = Mn;
- Mn = Mnp1;
- minim_pair = GSL_MIN_DBL(fabs(Mn) + fabs(Mnm1), minim_pair);
- }
- pair_ratio = start_pair/minim_pair;
- result->val = Mn;
- result->err = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(a)+1.0) * fabs(Mn);
- result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Mn);
- result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
- return GSL_ERROR_SELECT_2(stat_0, stat_1);
- }
- }
- else {
- /* x < 0
- * b < a
- */
- if(a <= 0.5*(b-x) || a >= -x) {
- /* Recurse down in b, from near the a=b line, b=a+eps,a+eps-1.
- */
- double N = floor(a - b);
- double eps = 1.0 + N - a + b;
- gsl_sf_result r_Manp1;
- gsl_sf_result r_Man;
- int stat_0 = hyperg_1F1_beps_bgt0(-eps, a+eps, x, &r_Manp1);
- int stat_1 = hyperg_1F1_beps_bgt0(1.0-eps, a+eps-1.0, x, &r_Man);
- double Manp1 = r_Manp1.val;
- double Man = r_Man.val;
- double Manm1;
- double n;
- double start_pair = fabs(Manp1) + fabs(Man);
- double minim_pair = GSL_DBL_MAX;
- double pair_ratio;
- double rat_0 = fabs(r_Manp1.err/r_Manp1.val);
- double rat_1 = fabs(r_Man.err/r_Man.val);
- for(n=a+eps-1.0; n>b+0.1; n -= 1.0) {
- Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0));
- Manp1 = Man;
- Man = Manm1;
- minim_pair = GSL_MIN_DBL(fabs(Manp1) + fabs(Man), minim_pair);
- }
- /* FIXME: this is a nasty little hack; there is some
- (transient?) instability in this recurrence for some
- values. I can tell when it happens, which is when
- this pair_ratio is large. But I do not know how to
- measure the error in terms of it. I guessed quadratic
- below, but it is probably worse than that.
- */
- pair_ratio = start_pair/minim_pair;
- result->val = Man;
- result->err = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Man);
- result->err *= pair_ratio*pair_ratio + 1.0;
- return GSL_ERROR_SELECT_2(stat_0, stat_1);
- }
- else {
- /* Pick a0 such that b ~= 2a0 + x, then
- * recurse down in b from a0,a0 to determine
- * the values near the line b=2a+x. Then recurse
- * forward on a from a0.
- */
- double epsa = a - floor(a);
- double a0 = floor(0.5*(b-x)) + epsa;
- double N = floor(a0 - b);
- double epsb = 1.0 + N - a0 + b;
- double Ma0b;
- double Ma0bp1;
- double Ma0p1b;
- int stat_a0;
- double Mnm1;
- double Mn;
- double Mnp1;
- double n;
- double err_rat;
- {
- gsl_sf_result r_Ma0np1;
- gsl_sf_result r_Ma0n;
- int stat_0 = hyperg_1F1_beps_bgt0(-epsb, a0+epsb, x, &r_Ma0np1);
- int stat_1 = hyperg_1F1_beps_bgt0(1.0-epsb, a0+epsb-1.0, x, &r_Ma0n);
- double Ma0np1 = r_Ma0np1.val;
- double Ma0n = r_Ma0n.val;
- double Ma0nm1;
- err_rat = fabs(r_Ma0np1.err/r_Ma0np1.val) + fabs(r_Ma0n.err/r_Ma0n.val);
- for(n=a0+epsb-1.0; n>b+0.1; n -= 1.0) {
- Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0));
- Ma0np1 = Ma0n;
- Ma0n = Ma0nm1;
- }
- Ma0bp1 = Ma0np1;
- Ma0b = Ma0n;
- Ma0p1b = (b*(a0+x)*Ma0b+x*(a0-b)*Ma0bp1)/(a0*b); /* right-down hook */
- stat_a0 = GSL_ERROR_SELECT_2(stat_0, stat_1);
- }
-
- /* Initialise the recurrence correctly BJG */
- if (a0 >= a - 0.1)
- {
- Mn = Ma0b;
- }
- else if (a0 + 1>= a - 0.1)
- {
- Mn = Ma0p1b;
- }
- else
- {
- Mnm1 = Ma0b;
- Mn = Ma0p1b;
- for(n=a0+1.0; n<a-0.1; n += 1.0) {
- Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
- Mnm1 = Mn;
- Mn = Mnp1;
- }
- }
- result->val = Mn;
- result->err = (err_rat + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Mn);
- return stat_a0;
- }
- }
- }
- /* Assumes b != integer
- * Assumes a != integer when x > 0
- * Assumes b-a != neg integer when x < 0
- */
- static
- int
- hyperg_1F1_ab_neg(const double a, const double b, const double x,
- gsl_sf_result * result)
- {
- const double bma = b - a;
- const double abs_x = fabs(x);
- const double abs_a = fabs(a);
- const double abs_b = fabs(b);
- const double size_a = GSL_MAX(abs_a, 1.0);
- const double size_b = GSL_MAX(abs_b, 1.0);
- const int bma_integer = ( bma - floor(bma+0.5) < _1F1_INT_THRESHOLD );
- if( (abs_a < 10.0 && abs_b < 10.0 && abs_x < 5.0)
- || (b > 0.8*GSL_MAX(fabs(a),1.0)*fabs(x))
- ) {
- return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
- }
- else if( x > 0.0
- && size_b > size_a
- && size_a*log(M_E*x/size_b) < GSL_LOG_DBL_EPSILON+7.0
- ) {
- /* Series terms are positive definite up until
- * there is a sign change. But by then the
- * terms are small due to the last condition.
- */
- return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
- }
- else if( (abs_x < 5.0 && fabs(bma) < 10.0 && abs_b < 10.0)
- || (b > 0.8*GSL_MAX_DBL(fabs(bma),1.0)*abs_x)
- ) {
- /* Use Kummer transformation to render series safe.
- */
- gsl_sf_result Kummer_1F1;
- int stat_K = gsl_sf_hyperg_1F1_series_e(bma, b, -x, &Kummer_1F1);
- int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
- Kummer_1F1.val, Kummer_1F1.err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_K);
- }
- else if( x < -30.0
- && GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x)
- ) {
- /* Large negative x asymptotic.
- * Note that we do not check if b-a is a negative integer.
- */
- return hyperg_1F1_asymp_negx(a, b, x, result);
- }
- else if( x > 100.0
- && GSL_MAX_DBL(fabs(bma),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.99*fabs(x)
- ) {
- /* Large positive x asymptotic.
- * Note that we do not check if a is a negative integer.
- */
- return hyperg_1F1_asymp_posx(a, b, x, result);
- }
- else if(x > 0.0 && !(bma_integer && bma > 0.0)) {
- return hyperg_1F1_U(a, b, x, result);
- }
- else {
- /* FIXME: if all else fails, try the series... BJG */
- if (x < 0.0) {
- /* Apply Kummer Transformation */
- int status = gsl_sf_hyperg_1F1_series_e(b-a, b, -x, result);
- double K_factor = exp(x);
- result->val *= K_factor;
- result->err *= K_factor;
- return status;
- } else {
- int status = gsl_sf_hyperg_1F1_series_e(a, b, x, result);
- return status;
- }
- /* Sadness... */
- /* result->val = 0.0; */
- /* result->err = 0.0; */
- /* GSL_ERROR ("error", GSL_EUNIMPL); */
- }
- }
- /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
- int
- gsl_sf_hyperg_1F1_int_e(const int a, const int b, const double x, gsl_sf_result * result)
- {
- /* CHECK_POINTER(result) */
- if(x == 0.0) {
- result->val = 1.0;
- result->err = 0.0;
- return GSL_SUCCESS;
- }
- else if(a == b) {
- return gsl_sf_exp_e(x, result);
- }
- else if(b == 0) {
- DOMAIN_ERROR(result);
- }
- else if(a == 0) {
- result->val = 1.0;
- result->err = 0.0;
- return GSL_SUCCESS;
- }
- else if(b < 0 && (a < b || a > 0)) {
- /* Standard domain error due to singularity. */
- DOMAIN_ERROR(result);
- }
- else if(x > 100.0 && GSL_MAX_DBL(1.0,fabs(b-a))*GSL_MAX_DBL(1.0,fabs(1-a)) < 0.5 * x) {
- /* x -> +Inf asymptotic */
- return hyperg_1F1_asymp_posx(a, b, x, result);
- }
- else if(x < -100.0 && GSL_MAX_DBL(1.0,fabs(a))*GSL_MAX_DBL(1.0,fabs(1+a-b)) < 0.5 * fabs(x)) {
- /* x -> -Inf asymptotic */
- return hyperg_1F1_asymp_negx(a, b, x, result);
- }
- else if(a < 0 && b < 0) {
- return hyperg_1F1_ab_negint(a, b, x, result);
- }
- else if(a < 0 && b > 0) {
- /* Use Kummer to reduce it to the positive integer case.
- * Note that b > a, strictly, since we already trapped b = a.
- */
- gsl_sf_result Kummer_1F1;
- int stat_K = hyperg_1F1_ab_posint(b-a, b, -x, &Kummer_1F1);
- int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
- Kummer_1F1.val, Kummer_1F1.err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_K);
- }
- else {
- /* a > 0 and b > 0 */
- return hyperg_1F1_ab_posint(a, b, x, result);
- }
- }
- int
- gsl_sf_hyperg_1F1_e(const double a, const double b, const double x,
- gsl_sf_result * result
- )
- {
- const double bma = b - a;
- const double rinta = floor(a + 0.5);
- const double rintb = floor(b + 0.5);
- const double rintbma = floor(bma + 0.5);
- const int a_integer = ( fabs(a-rinta) < _1F1_INT_THRESHOLD && rinta > INT_MIN && rinta < INT_MAX );
- const int b_integer = ( fabs(b-rintb) < _1F1_INT_THRESHOLD && rintb > INT_MIN && rintb < INT_MAX );
- const int bma_integer = ( fabs(bma-rintbma) < _1F1_INT_THRESHOLD && rintbma > INT_MIN && rintbma < INT_MAX );
- const int b_neg_integer = ( b < -0.1 && b_integer );
- const int a_neg_integer = ( a < -0.1 && a_integer );
- const int bma_neg_integer = ( bma < -0.1 && bma_integer );
- /* CHECK_POINTER(result) */
- if(x == 0.0) {
- /* Testing for this before testing a and b
- * is somewhat arbitrary. The result is that
- * we have 1F1(a,0,0) = 1.
- */
- result->val = 1.0;
- result->err = 0.0;
- return GSL_SUCCESS;
- }
- else if(b == 0.0) {
- DOMAIN_ERROR(result);
- }
- else if(a == 0.0) {
- result->val = 1.0;
- result->err = 0.0;
- return GSL_SUCCESS;
- }
- else if(a == b) {
- /* case: a==b; exp(x)
- * It's good to test exact equality now.
- * We also test approximate equality later.
- */
- return gsl_sf_exp_e(x, result);
- } else if(fabs(b) < _1F1_INT_THRESHOLD && fabs(a) < _1F1_INT_THRESHOLD) {
- /* a and b near zero: 1 + a/b (exp(x)-1)
- */
- /* Note that neither a nor b is zero, since
- * we eliminated that with the above tests.
- */
-
- gsl_sf_result exm1;
- int stat_e = gsl_sf_expm1_e(x, &exm1);
- double sa = ( a > 0.0 ? 1.0 : -1.0 );
- double sb = ( b > 0.0 ? 1.0 : -1.0 );
- double lnab = log(fabs(a/b)); /* safe */
- gsl_sf_result hx;
- int stat_hx = gsl_sf_exp_mult_err_e(lnab, GSL_DBL_EPSILON * fabs(lnab),
- sa * sb * exm1.val, exm1.err,
- &hx);
- result->val = (hx.val == GSL_DBL_MAX ? hx.val : 1.0 + hx.val); /* FIXME: excessive paranoia ? what is DBL_MAX+1 ?*/
- result->err = hx.err;
- return GSL_ERROR_SELECT_2(stat_hx, stat_e);
- } else if (fabs(b) < _1F1_INT_THRESHOLD && fabs(x*a) < 1) {
- /* b near zero and a not near zero
- */
- const double m_arg = 1.0/(0.5*b);
- gsl_sf_result F_renorm;
- int stat_F = hyperg_1F1_renorm_b0(a, x, &F_renorm);
- int stat_m = gsl_sf_multiply_err_e(m_arg, 2.0 * GSL_DBL_EPSILON * m_arg,
- 0.5*F_renorm.val, 0.5*F_renorm.err,
- result);
- return GSL_ERROR_SELECT_2(stat_m, stat_F);
- }
- else if(a_integer && b_integer) {
- /* Check for reduction to the integer case.
- * Relies on the arbitrary "near an integer" test.
- */
- return gsl_sf_hyperg_1F1_int_e((int)rinta, (int)rintb, x, result);
- }
- else if(b_neg_integer && !(a_neg_integer && a > b)) {
- /* Standard domain error due to
- * uncancelled singularity.
- */
- DOMAIN_ERROR(result);
- }
- else if(a_neg_integer) {
- return hyperg_1F1_a_negint_lag((int)rinta, b, x, result);
- }
- else if(b > 0.0) {
- if(-1.0 <= a && a <= 1.0) {
- /* Handle small a explicitly.
- */
- return hyperg_1F1_small_a_bgt0(a, b, x, result);
- }
- else if(bma_neg_integer) {
- /* Catch this now, to avoid problems in the
- * generic evaluation code.
- */
- gsl_sf_result Kummer_1F1;
- int stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &Kummer_1F1);
- int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
- Kummer_1F1.val, Kummer_1F1.err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_K);
- }
- else if(a < 0.0 && fabs(x) < 100.0) {
- /* Use Kummer to reduce it to the generic positive case.
- * Note that b > a, strictly, since we already trapped b = a.
- * Also b-(b-a)=a, and a is not a negative integer here,
- * so the generic evaluation is safe.
- */
- gsl_sf_result Kummer_1F1;
- int stat_K = hyperg_1F1_ab_pos(b-a, b, -x, &Kummer_1F1);
- int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
- Kummer_1F1.val, Kummer_1F1.err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_K);
- }
- else if (a > 0) {
- /* a > 0.0 */
- return hyperg_1F1_ab_pos(a, b, x, result);
- } else {
- return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
- }
- }
- else {
- /* b < 0.0 */
- if(bma_neg_integer && x < 0.0) {
- /* Handle this now to prevent problems
- * in the generic evaluation.
- */
- gsl_sf_result K;
- int stat_K;
- int stat_e;
- if(a < 0.0) {
- /* Kummer transformed version of safe polynomial.
- * The condition a < 0 is equivalent to b < b-a,
- * which is the condition required for the series
- * to be positive definite here.
- */
- stat_K = hyperg_1F1_a_negint_poly((int)rintbma, b, -x, &K);
- }
- else {
- /* Generic eval for negative integer a. */
- stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &K);
- }
- stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
- K.val, K.err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_K);
- }
- else if(a > 0.0) {
- /* Use Kummer to reduce it to the generic negative case.
- */
- gsl_sf_result K;
- int stat_K = hyperg_1F1_ab_neg(b-a, b, -x, &K);
- int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
- K.val, K.err,
- result);
- return GSL_ERROR_SELECT_2(stat_e, stat_K);
- }
- else {
- return hyperg_1F1_ab_neg(a, b, x, result);
- }
- }
- }
-
- #if 0
- /* Luke in the canonical case.
- */
- if(x < 0.0 && !a_neg_integer && !bma_neg_integer) {
- double prec;
- return hyperg_1F1_luke(a, b, x, result, &prec);
- }
- /* Luke with Kummer transformation.
- */
- if(x > 0.0 && !a_neg_integer && !bma_neg_integer) {
- double prec;
- double Kummer_1F1;
- double ex;
- int stat_F = hyperg_1F1_luke(b-a, b, -x, &Kummer_1F1, &prec);
- int stat_e = gsl_sf_exp_e(x, &ex);
- if(stat_F == GSL_SUCCESS && stat_e == GSL_SUCCESS) {
- double lnr = log(fabs(Kummer_1F1)) + x;
- if(lnr < GSL_LOG_DBL_MAX) {
- *result = ex * Kummer_1F1;
- return GSL_SUCCESS;
- }
- else {
- *result = GSL_POSINF;
- GSL_ERROR ("overflow", GSL_EOVRFLW);
- }
- }
- else if(stat_F != GSL_SUCCESS) {
- *result = 0.0;
- return stat_F;
- }
- else {
- *result = 0.0;
- return stat_e;
- }
- }
- #endif
- /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
- #include "gsl_specfunc__eval.h"
- double gsl_sf_hyperg_1F1_int(const int m, const int n, double x)
- {
- EVAL_RESULT(gsl_sf_hyperg_1F1_int_e(m, n, x, &result));
- }
- double gsl_sf_hyperg_1F1(double a, double b, double x)
- {
- EVAL_RESULT(gsl_sf_hyperg_1F1_e(a, b, x, &result));
- }
|