Końcowy krok w procesie tworzenia biblioteki numerycznej obliczającej relatywistyczne całki Fermiego-Diraca to testowanie. Ponieważ nie są znane żadne szczególne przypadki analityczne, to co można zrobić to porównanie jednej implementacji numerycznej z drugą. Jako ,,referencyjne'' zostało przyjęte bezpośrednie całkowanie adaptywnym algorytmem numerycznym działającym w dowolnej precyzji. Dla ujemnego eta (potencjału chemicznego) jest to kwadratura podwójnie eksponencjalna (Tanh-Sinh), dla dużego eta (silna degeneracja) adaptywny algorytm uruchamiany na trzech podprzedziałach: (0, eta-96), (eta-96,eta+96) oraz (eta+96, Infinity). Precyzja to 32 cyfry znaczące. Naiwne podejście polega na uruchomieniu obliczeń i czekaniu na wyniki.

Oprócz typowych zagadnień związanych z niereprezentowalnością liczb wymiernych w postaci zmiennoprzecinkowej typu double IEEE754 szybko zderzymy się z masą problemów. Większość ich wynika z zastosowania obliczeń dowolnej precyzji w masowej skali i związanego z tym przetwarzania symbolicznego. Licząc jedną całkę prawdopodobieństwo ich napotkania jest znikome. Licząc ich dziesiątki tysięcy staje się pewnością.

Problemy wynikają (prawdopodobnie) z rekursywnego zwiekszania precyzji numerycznej przez pewne elementarne algorytmy, np: obliczające funkcje typu Sin, Exp czy Log. Zdumiewająco łatwo algorytmy te potrafią zażądać kilkanaście GB pamięci, tylko po to aby ustalić w której ćwiartce należy obliczyć wartości f. trygonometrycznych. Jeżeli pamięci w systemie jest mniej, obliczenia de facto stają w miejscu. W najgorszym wariancie system operacyjny wiesza się lub przestaje odpowiadać. Sytuacja pogarsza się, gdy obliczenia chcemy wykonać na maszynie wielordzeniowej. Każdy z Kerneli, lub kilka naraz, mogą stanąć w miejscu albo ,,zjeść'' cały RAM. Spowolnione działanie systemu operacyjnego bywa uznane przez master Kernel jako zerwanie połączenia i przekierowanie tego samego rachunku na inny, wolny Kernel, co tylko pogarsza sprawę.

Dlatego kluczowe jest przewidzenie i odpowiednie przechwycenie powyższych sytuacji, czyli braku pamięci i braku czasu. W Mathematice służą do tego polecenia MemoryConstrained i TimeConstrained. Niestety, MemoryConstrained, nie blokuje użycia pamięci przez podprogramy wywoływane przez Kernel. Pierwsze testy na domowym komputerze z 16 GB RAM-u wywoływały błąd braku pamięci, który można było przechwycić poleceniem Catch[...., _SystemException] i w przypadku pojawienia się komunikatu SystemException["MemoryAllocationFailure"] zwrócić ten sam błąd co MemoryConstrained. Ten sam rachunek uruchomiony na maszynie z 376GB pamięci zawieszał system, bo Mathematica uznawała, że pobranie 25% tej pamięci jest nieszkodliwe. Ale takich Kerneli było kilkanaście...

Co w powyższej sytuacji można zrobić? Istnieje Linuksowe polecenie ulimit , które pozwala na ograniczenie pamięci dostępnej dla użytkownika. Standardowe ustawienia to ( ulimit -a ):
core file size (blocks, -c) 0
data seg size (kbytes, -d) unlimited
scheduling priority (-e) 0
file size (blocks, -f) unlimited
pending signals (-i) 31569
max locked memory (kbytes, -l) 65536
max memory size (kbytes, -m) unlimited
open files (-n) 1024
pipe size (512 bytes, -p) 8
POSIX message queues (bytes, -q) 819200
real-time priority (-r) 0
stack size (kbytes, -s) 8192
cpu time (seconds, -t) unlimited
max user processes (-u) 31569
virtual memory (kbytes, -v) unlimited
file locks (-x) unlimited
Mieszanie ustawieniami ulimit -m, ulimit -v oraz ulimit -S spowodowało uboczne skutki, w postaci ,,zabijania'' procesów tych Kerneli, które ośmieliły się zażądać zbyt dużo RAM-u. Nie o to chodziło.

Pomoc nadeszła z nieoczekiwanej strony. Przeszukując internet, natrafiłem na post niejakiego A. Odrzywołka sprzed prawie 3 lat, opisującego podobny problem. Otóż istnieją w Mathematice nieudokumentowane oficjalnie ustawienia Internal`$MaxExponent oraz Internal`$MinExponent określające jak głęboko może sięgać rekursywne zwiększanie zakresu obliczeń zmiennoprzecinkowych dowolnej precyzji. Są one ustawione, oczywiście, na Infinity. Wystarczy je zmienić na coś bardziej realistycznego, powiedzmy Internal`$MaxExponent = 4096 i problemy z pożeraniem pamięci znikają jak ręką odjął. Skutkiem ubocznym jest brak zbieżności algorytmu całkowania oraz niemożność wykonania obliczeń, gdy f. podcałkowa przyjmuje wartości poniżej/powyżej limitu wykładnika, czyli znaczące ograniczenie zakresu dostępnych warości referencyjnych. Prace nad zoptymalizowaniem kosztów obliczeniowych i zasobów względem oczekiwanych wyników trwają.

Na koniec omawiana funkcja w wersji naiwnej:
F[k_, Eta_, Theta_] :=
NIntegrate[x^k Sqrt[1 + Theta x/2]/(1 + Exp[x - Eta]), {x, 0, Infinity}, WorkingPrecision -> 32]
versus f. opakowana licznymi bezpiecznikami:
F[k_?NumericQ, Eta_?NumericQ, Theta_?NumericQ] :=
Block[{$MaxExtraPrecision = 8192, Internal`$MaxExponent = 4096, Internal`$MinExponent = -4096},
wynik = Catch[
MemoryConstrained[
TimeConstrained[
Check[
If[Eta > 4096,
NIntegrate[
x^k Sqrt[1 + Theta x/2]/(1 + Exp[x - Eta]), {x,
0, Eta - 3*32, Eta + 3*32, Infinity},
WorkingPrecision -> 32, MaxRecursion -> 12,
Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 8}],

NIntegrate[
x^k Sqrt[1 + Theta x/2]/(1 + Exp[x - Eta]), {x, 0,
Infinity}, WorkingPrecision -> 32, MaxRecursion -> 16,
Method -> {"DoubleExponential", "SymbolicProcessing" -> 8}]
],
-2]
, 8, -1],
1024*1024*1024, -4],
_SystemException];
If[wynik === SystemException["MemoryAllocationFailure"], -4, wynik]
];