Hypergeometric functions

Specific Domains Numerics
,
1

Hello,
I need to evaluate the Gauss Hypergeometric functions _2F_1(a,b;c;z). Here are some specifications for the arguments of this function

  • a has compex values. The real part is \sim2 (and is strictly positive), while the complex part ranges from 0 to \sim500
  • b has compex values. The real part ranges from 10 to \sim5000, while the complex part ranges from 0 to \sim500
  • c is a real number, which ranges from 10 to 5000
  • z is a real number, which ranges from 0 to 1.

Which packages should I consider to evaluate _2F_1, given the previous requirements? I have found some existing packages, such as HypergeometricFunctions.jl or Nemo.jl, but I don’t know which one is the most stable/efficient given my requirements. I could also consider to call the mpmath module.
I also add that I value stability and precision, rather then speed. I could evaluate the hypergeometric functions for a fixed grid of a, b, c and z and then I could store the output.

If requested, I can give more details on the numerical integrals I am trying to solve.

Thank you for your help,
Marco

2 Likes
2

Specialfunctions.jl would be the first package I would look at, but it seems it does not support hypergeometric functions although the name would fit. :wink:

The next thing would probably be searching on Juliahub.
Hypergeometricfunctions.jl shows up first, but there is also StatsFuns and the Distributions package.

HypergeometricFunctions.jl looks tike the package to use, but opening an issue for SpecialFunctions.jl might also be a good idea.

2 Likes
3

Hi @marcobonici, HypergeometricFunctions.jl is a package based on

  • N. Michel and M. V. Stoitsov, Fast computation of the Gauss hypergeometric function with all its parameters complex with application to the Pöschl–Teller–Ginocchio potential wave functions, Comp. Phys. Commun., 178:535–551, 2008;
  • J. W. Pearson, S. Olver and M. A. Porter, Numerical methods for the computation of the confluent and Gauss hypergeometric functions, Numer. Algor., 74:821–866, 2017; and,
  • research on rational approximation algorithms that is not yet published.

The repository is not rigorous, but designed to be generic & flexible in the Julia-sense. (As a general comment, the parameter regime you’re interested in is problematic for stability and efficiency.)

If timing is no issue when creating your tables, you can always iterate over the precision until you are satisfied.

julia> using HypergeometricFunctions

julia> bits = 64
64

julia> setprecision(bits)
64

julia> ab = Complex{BigFloat}[2+500im, 2000+500im]
2-element Vector{Complex{BigFloat}}:
    2.0 + 500.0im
 2000.0 + 500.0im

julia> c = BigFloat[5000]
1-element Vector{BigFloat}:
 5000.0

julia> z = BigFloat(0.75)
0.75

julia> val = HypergeometricFunctions.pFqweniger(ab, c, z)
2.26285976373047298522e-18 - 3.87520480967281999701e-20im

julia> err = 1
1

julia> while err > 1e-100
           bits *= 2
           setprecision(bits)
           valnew = HypergeometricFunctions.pFqweniger(ab, c, z)
           err = abs(val-valnew)/abs(valnew)
           val = valnew
       end

julia> val
-1.260430203099602746599042376553865089118791199122916608865044929767940602836103359950826080051784754262000964653796917697258667174883250490398306584988606846792349409222203304725508986614660980506497947620292278913430773931273141234716825404266226618075265239797537643758196575342689058294488706792357746118743e-26 - 8.939487174343187315364167524433463458456915010715489477551635823318511013823554989260856845403228828796408064308487650539525096344830318147568045736429366146064820295315216835446407078108498963617764757578632658672273341724744581851370670269520147673054699790301299748194494370020751574979689007127930265104835e-27im

julia> bits
1024

julia> err
2.076078402922008073651649498895162272423799332685727180642421898948357431880875216390476692009157211189582065628002185295264430187924787259202316778393986184638156636159256782253864740290441853318534610756824276539837983656847940295587318289239282139268008182178475172235273054373920924841419749631097045267134e-127

julia> setprecision(2*bits)
2048

julia> valnew = HypergeometricFunctions.pFqweniger(ab, c, z)
-1.26043020309960274659904237655386508911879119912291660886504492976794060283610335995082608005178475426200096465379691769725866717488325049039830658498860684679234940922220330472550898661466098050649794762029227891343077393127314123471682540426622661807526523979753764375819657534261569928633223802354185464494766154356210819164851682735102308511168394254891629575747602082019751674097000612659852782541543187283051406119091268392627596536663849616173443388237855365033290101898906806866581073585456999741221237330897093677080632213302097866634063357324933021607334383575670291533241974365732598650727036823360449690715e-26 - 8.93948717434318731536416752443346345845691501071548947755163582331851101382355498926085684540322882879640806430848765053952509634483031814756804573642936614606482029531521683544640707810849896361776475757863265867227334172474458185137067026952014767305469979030129974819449437002083973023625932412489164085368327216895265713555927131363492307260309435071862316501683001952084907235694436090440504302825122722489400840488043243828999965455776325815456521543953804408993972162880630181260695521646554317594704584093357149321300308555709736707677357654506616735740843905578502021970625061659865515388282141676580517129449e-27im

julia> err = abs(val-valnew)/abs(valnew)
4.7815116374998721567868315521174153185485502894089632676089111656582693365399889037597171865806283497362103510757133134147391427406941262720724544482094324435197607776174553738763558509147494239172058700996953831291772660269475041149742146415317156429143626450415927245551344591547272099486309765129051745410803399821410144781863897956570143995873365635339243621007401773495519463044205561687377848034212284711504404377626436158196398780889673986572733672479233669247066863005421502031258354660773806141522987149577343472110951667313541558437584219864762958647100028190851057065469344294922605175527829315613973088368e-281

julia>
5 Likes
4

Nemo will be rigorous and do what you want. It calls into Arb, which is developed by @Fredrik_Johansson (also the author of mpmath). He might have more insight into the efficiency for your parameter range.

2 Likes
5

Hi, thank you all for your kind answers.

In the next few days I will try to test all the packages you have suggested and I’ll report here what I will find.
Cheers,
Marco

1 Like
6

As a wrapper of the GNU Scientific Library, GSL.jl also has Hypergeometric functions.

2 Likes
7

Nemo.jl (Arb) is strictly superior to mpmath here, so no need to compare mpmath.

Arb shouldn’t have any trouble for these inputs, though you will likely have to set the precision adaptively.

4 Likes