Idea stworzenia biblioteki obliczającej relatywistyczne całki Fermiego-Diraca w całym zakresie parametrów dopuszczalnych przez liczby zmiennoprzecinkowe podwójnej precyzji okazała się bardziej pracochłonna niż to zakładałem. Oczywiście! Kwadratura tanh-sinh nie wystarcza. Dla ekstremalnych wartości degeneracji konieczne jest użycie przybliżenia Sommerfelda . Niestety, wyraz wiodący zawiera funkcje hipergeometryczne. Tak samo jest w przypadku niezdegenerowanym, ale tym razem współczynniki szeregu zawierają konfluentną funkcję hipergeometryczną Kummera U(1+k, 5/2+k, x), dla k>-1 i x>0.

Zasadniczo U(a,b,z) można obliczyć dla małego |z|<16 jako różnicę (tzw. hiperkombinacja) dwóch funkcji Hypergeometric1F1, dla których z kolei działa standardowy szereg hipergeometryczny. Za wyjątkiem przypadku połówkowego k, np: k=-1/2, 1/2, 3/2 itd. Zbieżność szeregu jest dziwaczna, jest on szybkozbieżny dla wszystkich innych k, a rozbieżny dla np: k=1/2. Dokładność (względna) przybliżenia przez szereg ma więc dziwaczną postać gładkiej funkcji poprzerywanej ,,pochodnymi delty Diraca''. Matematycy łatwo przechodzą nad tym do porządku, argumentując, że ponieważ U(1+k,5/2+k,x) jest ciągła w k, to wystaczy obliczyć ją dla k+-epsilon i wyciągnąć średnią. W praktyce oblicza się jeszcze pochodną. Ale w ten sposób nie udało mi się osiągnąć wymaganej dokładności na poziomie kilku ULP.

Inny sposób wynika z obserwacji, że dla k=-1/2,1/2,3/2,... U(1+k,5/2+k,x) wyraża się dosyć skomplikowaną kombinacją funkcji BesselK, rzędu zero i jeden. Na pierwszy rzut oka wydaje się, że nie ma w tym żadnej reguły. Ale dzięki NIST DLMF można dowiedzieć się, że istnieje dwustopniowa rekurencja ,,w górę'':
U(k) = ( U(k-2) + (0.5+k-x) U(k-1) )/k/x.
Teraz wystarczy podać wzory dla U(-1/2) oraz U(1/2) i gotowe.

Niestety, K0(x) i K1(x) nie są w standardowej bibliotece matematycznej . Istnieje kilka implementacji, opartych jak mi się wydaje, na przeskalowaniu przez exp(x) oraz aproksymacji Czybyszewa. Kody GNU GSL wyglądają koszmarnie i trudno się w nich połapać, zwłaszcza kto od kogo je skopiował. Na szczęście od czasu do czasu na śmietniku można znaleźć radzieckie podręczniki, np: Tichonow, Samarski, Równania fizyki matematycznej, PWN 1963. Otóż w dodatku, na str. 538, znajdujemy wzór (22), który jest ewidentnie słabo znany:

Jest to oczywiście wzór w postaci gotowej dla metody podwójnie eksponencjalnej, czyli metody trapezów bez końców przedziału. Użyłem tu dokładnie tego samego kodu co dla całki Fermiego-Diraca, podniemiając jedynie funkcję podcałkową. Dokładność jest imponująca. Dla 3 poziomów rekursji, przy kroku 2^-3 = 1/8 = 0.125 prawie nigdy nie przekracza sie błędu 128ULP. Prawie, czyli wszędzie poda obszarem denormalnym . Tu być może asymptotyka byłaby lepsza, ale na razie dosyć!

Dysponując funkcjami K0(x) i K1(x) możemy wrócić do pierwotnego podproblemu obliczenia U(1+k,5/2+k,x). Testowanie wykazało, że metoda działa, ale tylko dla małych x<16. Dla x>>16, dostajemy wyrażenia typu 0*infinity. Same wartości funkcji BesselK(n,x) są zbyt małe dla x>700 aby dało się je przedstawić w typie double. Dlatego istniejące implementacje występują w dwóch wersjach (GNU GSL) albo mają flagi (Fortran) włączające skalowanie przez exp(x/2). Na szczęście okazuje się, że szereg asymptotyczny U(a,b,x) jest niewrażliwy na problemy z połówkowymi wartościami i problem mamy z głowy.

I to byłoby na tyle w kwestii funkcji hipergeometrycznej U(1+k,5/2+k,x), gdyby nie jedna zadziwiająca obserwacja. Otóż błąd powyższej implementacji w przypadku gdy x jest liczbą naturalną jest 10000 większy niż dla pozostałych liczb zmiennoprzecinkowych. Nie mam dla tego żadnego sensownego wyjaśnienia.