Source code for sarkas.utilities.fdints
"""
Module of Numba'd functions for the calculation of the unnormalized Fermi-Dirac integrals of
half integer order p = -9/2,(1),21/2 and integer order p = 0,(1),10 via piecewise minimax rational approximation.\n
The Numba'd is not necessary for a one time calculation, however, we decided to add the Numba functionality so that
future users can use them in Numba'd functions.\n
The functions contained here are a Python translation of the code xfdh.txt :cite:`Fukushima2015b`.\n
All the credits go to Fukushima, T. email: Toshio.Fukushima at nao.ac.jp. \n
It was easier than using other packages and/or write our own.
"""
from numba import jit
from numba.core.types import float64
from numpy import exp, log, sqrt
from .exceptions import SarkasError
# __all__ = [
# "fdm1h",
# "fd1h",
# "fd3h",
# "fd5h",
# "fd7h",
# "fd9h",
# "fd11h",
# "fd13h",
# "fd15h",
# "fd17h",
# "fd19h",
# "fd21h",
# ]
[docs]def fd_doc_hparams(order):
frac = order.as_integer_ratio()
def fd_doc(func):
func.__doc__ = f"""
Double precision rational minimax approximation of Fermi-Dirac integral
of half integer order p={frac[0]}/{frac[1]}. Ref. :cite:`Fukushima2015b`.\n
Fukushima, T. email: Toshio.Fukushima at nao.ac.jp\n
Parameters
----------
x : float
Non-dimensional chemical potential, :math:`\\mu/(k_BT)`.
Returns
-------
fd : float
Value of the Fermi-Dirac integral.
"""
return func
return fd_doc
[docs]def fd_doc_iparams(order):
def fd_doc(func):
func.__doc__ = f"""
Double precision rational minimax approximation of Fermi-Dirac integral
of integer order p={order}. Ref. :cite:`Fukushima2015b`.\n
Fukushima, T. email: Toshio.Fukushima at nao.ac.jp\n
Parameters
----------
y : float
Non-dimensional chemical potential, :math:`\\mu/(k_BT)`.
Returns
-------
fd : float
Value of the Fermi-Dirac integral.
"""
return func
return fd_doc
[docs]def fermidirac_integral(p: float, eta: float) -> float:
"""Wrapper function to calculate the Fermi-Dirac integral using approximation found in Ref. :cite:`Fukushima2015b`.\n
Parameters
----------
p : float
Order of the unnormalized Fermi-Dirac integral.
eta: float
Value of the normalized chemical potential :math:`\\mu/(k_BT)`.
Returns
-------
fd : float
Value of the Fermi-Dirac integral.
Raises
------
`~sarkas.utilities.exceptions.SarkasError`
If the order of the integral is not supported.
Examples
--------
>>> from numpy import pi, sqrt
>>> # Electron density
>>> eta = -0.2860
>>> I = fermidirac_integral(p = 0.5, eta = eta)
>>> lambda_deB = 1.957093e-11 # [m]
>>> ne = 4.0/( sqrt(pi) * lambda_deB**3 ) * I
>>> f"{ne:.4e}"
'1.6201e+32'
"""
orders = [-9 / 2, -7 / 2, -5 / 2, -3 / 2, -1 / 2]
for ord in range(0, 22):
orders.append(ord / 2)
try:
if p in orders:
if p == -9 / 2:
fd = fdm9h(eta)
elif p == -7 / 2:
fd = fdm7h(eta)
elif p == -5 / 2:
fd = fdm5h(eta)
elif p == -3 / 2:
fd = fdm3h(eta)
elif p == -1 / 2:
fd = fdm1h(eta)
elif p == 0.0:
fd = fd0h(eta)
elif p == 1 / 2:
fd = fd1h(eta)
elif p == 1.0:
fd = fd2h(eta)
elif p == 3 / 2:
fd = fd3h(eta)
elif p == 4 / 2:
fd = fd4h(eta)
elif p == 5 / 2:
fd = fd5h(eta)
elif p == 6 / 2:
fd = fd6h(eta)
elif p == 7 / 2:
fd = fd7h(eta)
elif p == 8 / 2:
fd = fd8h(eta)
elif p == 9 / 2:
fd = fd9h(eta)
elif p == 10 / 2:
fd = fd10h(eta)
elif p == 11 / 2:
fd = fd11h(eta)
elif p == 12 / 2:
fd = fd12h(eta)
elif p == 13 / 2:
fd = fd13h(eta)
elif p == 14 / 2:
fd = fd14h(eta)
elif p == 15 / 2:
fd = fd15h(eta)
elif p == 16 / 2:
fd = fd16h(eta)
elif p == 17 / 2:
fd = fd17h(eta)
elif p == 18 / 2:
fd = fd18h(eta)
elif p == 19 / 2:
fd = fd19h(eta)
elif p == 20 / 2:
fd = fd20h(eta)
elif p == 21 / 2:
fd = fd21h(eta)
return fd
except:
raise SarkasError(f"Fermi-Dirac Integral of order {p} not yet implemented.")
@jit(float64(float64), nopython=True)
def invfd1h(u: float) -> float:
"""Approximate the inverse of the Fermi-Dirac integral :math:`I_{-1/2}(\\eta)` using the fits provided by Fukushima.
Function translated from Fukushima's code, see :cite:`Fukushima2015a`.
Parameters
----------
u : float
Normalized electron density :math:`=\\Lambda_{\\rm deB}^3 n_e \\sqrt{\\pi}/4`, \n
where :math:`\\Lambda_{\\rm deB}` is the de Broglie wavelength of the electron gas.
Examples
--------
>>> import numpy as np
>>> # Values taken from Tutorial Notebooks
>>> ne = 1.62e32 # [N/m^3]
>>> lambda_deB = 1.957093e-11 # [m]
>>> u = lambda_deB**3 * ne * np.sqrt(np.pi)/4.0
>>> eta = invfd1h(u)
>>> f"{eta:.4f}"
'-0.2860'
Returns
-------
difdih : float
Scaled chemical potential.
"""
u1 = 5.8936762325936050502
u2 = 20.292139527268839575
u3 = 69.680386701202787485
u4 = 246.24741525281437791
vc = 1543.4606939459401869
if u < u1:
v = u1 - u
r = (
u
* (
4.4225399543845577739e9
+ u
* (
+1.4318826531216930391e9
+ u
* (
+2.0024511162084252731e8
+ u
* (
+1.5771885953346837109e7
+ u
* (
+7.664073281017674960e5
+ u
* (
+2.3599362498847900809e4
+ u * (+4.4114601741712557348e2 + u * (+3.8844451277821727933e0))
)
)
)
)
)
)
/ (
+2.448084356710615572e9
+ v
* (
+2.063907695769060888e8
+ v
* (
+6.943821586626002503e6
+ v
* (+8.039397005856418743e4 + v * (-1.791261676399435220e3 + v * (-7.908051927048792349e1 - v)))
)
)
)
)
difd1h = log(r)
elif u < u2:
y = u - u1
difd1h = (
+5.849893914158469793e14
+ y
* (
+3.353389340896418967e14
+ y
* (
+7.300790845633384552e13
+ y
* (
+7.531271098292146768e12
+ y * (+3.726221594134586141e11 + y * (+7.827935737269045014e9 + y * (+5.021972425404123509e7)))
)
)
)
) / (
+1.4397298668421281743e14
+ y
* (
+6.440007889067504875e13
+ y
* (
+1.0498882082904393876e13
+ y
* (
+7.568424788316453035e11
+ y
* (
+2.3206228235771973103e10
+ y
* (
+2.4329782896354397638e8
+ y * (+3.6768664133860359837e5 + y * (-5.924317283823514482e2 + y))
)
)
)
)
)
)
elif u < u3:
y = u - u2
difd1h = (
+6.733834344762314874e18
+ y
* (
+1.1385291167086018856e18
+ y
* (
+7.441797125810403052e16
+ y
* (
+2.3556527595722738211e15
+ y
* (+3.6904107711114070061e13 + y * (+2.5927357055940595308e11 + y * (+5.989403440741097470e8)))
)
)
)
) / (
+6.968777783221497285e17
+ y
* (
+9.451599633557071205e16
+ y
* (
+4.7388759083089595117e15
+ y
* (
+1.0766510215928549449e14
+ y
* (
+1.0888539870400255904e12
+ y
* (
+4.0374047390260294467e9
+ y * (+2.3126814357531839818e6 + y * (-1.4788294703774470115e3 + y))
)
)
)
)
)
)
elif u < u4:
y = u - u3
difd1h = (
+7.884494095314249799e19
+ y
* (
+3.7486465573810023777e18
+ y
* (
+6.934193474730824900e16
+ y
* (
+6.302949477641708425e14
+ y
* (
+2.9299316609051704688e12
+ y
* (+6.591658047866512380e9 + y * (+6.082995857672390394e6 + y * (+1.5054843420905807932e3)))
)
)
)
)
) / (
+3.5593247304804720533e18
+ y
* (
+1.3505797700306451874e17
+ y
* (
+1.9160919212553016350e15
+ y
* (
+1.2652560651095328402e13
+ y
* (+3.9491055033213850687e10 + y * (+5.253083775042776560e7 + y * (+2.2252541165920236251e4 + y)))
)
)
)
)
else:
t = vc * (u ** -(4.0 / 3.0))
s = 1.0 - t
w = (+3.4330125059142833612e7 + s * (+8.713462091032493289e5 + s * (+2.4245560148256419080e3 + s))) / (
t * (+1.2961677595919532465e4 + s * (+3.2092883892793106635e2 + s * (+0.7192193760323717351e0)))
)
difd1h = sqrt(w)
return difd1h
@fd_doc_hparams(order=-9 / 2)
@jit(float64(float64), nopython=True)
def fdm9h(x: float) -> float:
factor = -2.0 / 7.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
0.270088205852269109
- ex
* (
98457.373700674016
+ t
* (
7895.1716828987114
+ t
* (
1170.39154599810719
+ t * (84.347808600628057 + t * (1.94779771326323869 + t * 0.00062804476066590963))
)
)
)
/ (
32220.8981607867840
+ t
* (
20608.4545396913128
+ t * (5234.8569367171782 + t * (658.75071261067250 + t * (40.9380138615149699 + t)))
)
)
)
elif x < 0.0:
s = -0.5 * x
fd = (
(
-1.8432713155273516
+ s
* (
-7.2420131297087391
+ s
* (
-2.2155513134917048
+ s
* (
0.04304761419364448
+ s
* (
-0.16668060134499780
+ s
* (
0.10403025949642735
+ s
* (
0.016505424797869711
+ s
* (
-0.0063406328104926592
+ s * (-0.0009078668865214013 + s * 0.00024692958093140298)
)
)
)
)
)
)
)
)
/ (
129.8291336259148
+ s
* (
-1.9134679779378
+ s
* (
249.6123960853833
+ s
* (
-3.3454980273582
+ s
* (
190.7566926415911
+ s
* (
-2.2190711173347
+ s
* (
72.29929273955658
+ s
* (-0.66382343556531 + s * (13.551723146833901 + s * (-0.0757657170023758 + s)))
)
)
)
)
)
)
)
* (x + 1.82715356570827096)
)
elif x < 2.0:
t = 0.5 * x
fd = (
(
27.9376225634323064
+ t
* (
17.6954702922794847
+ t
* (
-18.2160790482833106
+ t
* (
-2.57045608672749918
+ t
* (
2.20866272949862340
+ t
* (
-0.104202929779896185
+ t * (-0.177712759450918142 + t * (0.055743070525419128 - t * 0.0059642251814009821))
)
)
)
)
)
)
/ (
599.97601432112870
+ t
* (
-64.538219224509587
+ t
* (
1048.07497435809354
+ t
* (
-104.535376308139931
+ t
* (
694.31701140352634
+ t
* (
-61.402988847817226
+ t
* (207.271703307329581 + t * (-14.7241003912299059 + t * (23.5639226214135345 - t)))
)
)
)
)
)
)
* (x - 0.557104412737456300)
) # care for zero point
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
fd = (
-(
0.105736863491337914
+ t
* (
0.182686468618144972
+ t
* (
-0.0830333880807117100
+ t
* (
0.0123471104343183151
+ t
* (
0.0213438566082710205
+ t
* (
-0.0431206240183745280
+ t
* (
0.0311190249582737067
+ t
* (
-0.0115670162945002128
+ t * (0.00222720333535487902 - t * 0.000179769912955733772)
)
)
)
)
)
)
)
)
/ (
11.4369029896970731
+ t
* (
29.5790675947434036
+ t
* (
55.2958932613908644
+ t
* (
59.7192644653455683
+ t
* (
51.8646589384828875
+ t
* (
26.5061101661226671
+ t
* (
11.1515306816966566
+ t * (-0.140562900980389991 + t * (-0.566848363630388431 - t))
)
)
)
)
)
)
)
* (x - 3.79898975457237574)
) # care for zero point
elif x < 10.0:
t = 0.2 * x - 1.0
fd = -(
0.0131615741211877996
+ t
* (
0.0630819955070757629
+ t
* (
0.0347025335431948214
+ t
* (
0.00659613983600166143
+ t
* (
0.0239129779624141424
+ t
* (
0.000528952248906432005
+ t
* (
0.00411353218448862483
+ t
* (
-0.00119448333536781374
+ t
* (
0.000497581590580634773
+ t * (-0.0000934288845342437109 + t * 7.29489294899966855e-6)
)
)
)
)
)
)
)
)
) / (
12.2828253912160684
+ t
* (
63.2143680817101004
+ t
* (
155.880993448869246
+ t
* (
239.067686638304706
+ t
* (
254.656949674257383
+ t
* (
200.134181826630681
+ t
* (
120.784221710133443
+ t
* (56.5037249583191013 + t * (20.2196452284175417 + t * (5.22412387722958575 + t)))
)
)
)
)
)
)
)
elif x < 20.0:
t = 0.1 * x - 1.0
fd = -(
0.000243271785585373458
+ t
* (
0.000534699823019827761
+ t
* (
0.000824373965146324914
+ t
* (
0.00106321349769688454
+ t
* (
0.000488890808530810197
+ t
* (
0.000510404811590977338
+ t
* (
0.0000543711407291589868
+ t
* (
0.0000597066095966549458
+ t
* (
-9.52409541937957252e-6
+ t * (1.26499979908255436e-6 - t * 8.23879904655758674e-8)
)
)
)
)
)
)
)
)
) / (
1.89001039755418574
+ t
* (
12.0291520111708943
+ t
* (
36.3016148236874403
+ t
* (
69.0194326307265861
+ t
* (
92.5327210822279397
+ t
* (
92.3443444983825834
+ t
* (
70.3437819943334403
+ t
* (41.0546070181362401 + t * (17.9813587502434782 + t * (5.52467389311001165 + t)))
)
)
)
)
)
)
)
elif x < 40.0:
t = 0.05 * x - 1.0
fd = -(
5.76414987722574846e-6
+ t
* (
0.0000175838621123388022
+ t
* (
0.0000306497936833891724
+ t
* (
0.0000245465851262614617
+ t
* (
9.68986921084327441e-6
+ t
* (
-9.60429091051474115e-6
+ t
* (
-2.71258276002580026e-6
+ t
* (
2.14306002093744082e-7
+ t * (-2.02878733208210684e-8 + t * 1.11972866705178940e-9)
)
)
)
)
)
)
)
) / (
0.672666709110930353
+ t
* (
4.51114986457223474
+ t
* (
14.2130577502093159
+ t
* (
27.0924755033205922
+ t
* (
33.3161322251216687
+ t
* (
25.3582658290231325
+ t * (8.61719049364763904 + t * (-2.98613950390747585 + t * (-3.72974082078125757 - t)))
)
)
)
)
)
)
else:
w = 1.0 / (x * x)
s = 1.0 - 1600.0 * w
fd = (
factor
/ (sqrt(x) * x * x * x)
* (
1.0
+ w
* (
20360.4093941608469
+ s
* (
48690.0157659639679
+ s
* (
56148.8553220593957
+ s * (16506.5714937251887 + s * (1066.30280063583769 + s * 3.73715435019239458))
)
)
)
/ (
765.319028032818290
+ s
* (
1851.12566867527470
+ s * (2160.21729234460541 + s * (677.261952266814024 + s * (56.0179124057448609 + s)))
)
)
)
)
return fd
@fd_doc_hparams(order=-7 / 2)
@jit(float64(float64), nopython=True)
def fdm7h(x: float) -> float:
factor = -2.0 / 5.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
-0.945308720482941881
+ ex
* (
22521.6817696480830
+ t
* (
3006.50943115266211
+ t * (223.777970320362539 + t * (7.1096231598119406 - t * 0.00186000064211079826))
)
)
/ (
4211.64894546505314
+ t * (2132.92301652550982 + t * (400.935813755588393 + t * (33.0304914363496116 + t)))
)
)
elif x < 0.0:
s = -0.5 * x
fd = (
(
5.49736264373955718
+ s
* (
3.93868182790099018
+ s
* (
0.00680535603918657089
+ s
* (
-0.0610218211998865281
+ s
* (
-0.195616572333524455
+ s
* (
0.0172573301744899007
+ s
* (
0.0148839062617757320
+ s * (-0.00259857880173559174 + s * 0.0000170523346905109458)
)
)
)
)
)
)
)
/ (
48.4238197836838959
+ s
* (
8.18206322265668686
+ s
* (
73.9370055146725832
+ s
* (
10.7170918787834126
+ s
* (
42.0907078084017798
+ s * (4.72976153773643018 + s * (10.5965223842451046 + s * (0.704336166377684061 + s)))
)
)
)
)
)
* (x + 0.731435761340666539)
) # care for zero point
elif x < 2.0:
t = 0.5 * x
fd = (
143.837264704139569
+ t
* (
322.683997310438763
+ t
* (
-209.374575032955500
+ t
* (
-54.681932464048843
+ t
* (
-30.3617985514093593
+ t
* (
-5.7605448223972412
+ t
* (
3.65745825030136209
+ t * (0.145440718467921945 + t * (-0.340001125998911237 + t * 0.0446577756904410773))
)
)
)
)
)
)
) / (
1732.20737492411446
+ t
* (
97.952930820143076
+ t
* (
2660.50150777614871
+ t
* (
115.730278615026183
+ t
* (
1524.22026694958607
+ t
* (
40.4116469657132500
+ t * (385.907052068138736 + t * (1.45068964456455148 + t * (36.4544527882983039 - t)))
)
)
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
fd = (
-(
66.1873851488861175
+ t
* (
57.9197865216898287
+ t
* (
-10.6658833177171402
+ t
* (
10.2066946495410176
+ t
* (
1.19595124198914370
+ t
* (
-1.81286879962085966
+ t
* (
1.52800227451442796
+ t
* (
-0.546498161708704000
+ t * (0.0988077919447287808 - t * 0.00745279173673929794)
)
)
)
)
)
)
)
)
/ (
1525.13818882033488
+ t
* (
3966.07893650424194
+ t
* (
6837.27734004268945
+ t
* (
7086.45653543153192
+ t
* (
5612.71274301484627
+ t
* (
2928.56588821761611
+ t
* (1256.17847516062581 + t * (310.558929999645079 + t * (89.3091775022809646 - t)))
)
)
)
)
)
)
* (x - 2.59355543171509487)
) # care for zero point
elif x < 10.0:
t = 0.2 * x - 1.0
fd = -(
0.184179377827614791
+ t
* (
0.516542442597467614
+ t
* (
0.292799259677321999
+ t
* (
0.166303092691678597
+ t
* (
0.258615034442190300
+ t
* (
-0.0303157249245734253
+ t
* (
0.0373447714795897992
+ t
* (
-0.00332926057033243781
+ t * (0.000563408857595472728 - t * 0.0000334037183021934960)
)
)
)
)
)
)
)
) / (
18.2615126816087432
+ t
* (
85.1686926876657081
+ t
* (
181.378990080239764
+ t
* (
232.528454347497368
+ t
* (
200.564064511137621
+ t
* (
125.206874911765040
+ t * (58.6883968744415564 + t * (20.6768210567814309 + t * (5.02620698086580395 + t)))
)
)
)
)
)
)
elif x < 20.0:
t = 0.1 * x - 1.0
fd = -(
0.00521027824483942692
+ t
* (
0.0159223834333424812
+ t
* (
0.0290786393625134859
+ t
* (
0.0350554456604475612
+ t
* (
0.0276022017688728517
+ t
* (
0.0190056432579866130
+ t
* (
0.00725681379911524042
+ t
* (
0.00278845414163336780
+ t
* (
0.000285395695659601419
+ t * (-9.95541819385254380e-6 + t * 2.96666654254966665e-7)
)
)
)
)
)
)
)
)
) / (
3.40047012308625195
+ t
* (
20.3896764592442369
+ t
* (
58.0983301256082527
+ t
* (
104.466101980605120
+ t
* (
132.574327309248600
+ t
* (
125.279491711147983
+ t
* (
90.3637198347191251
+ t
* (49.8942487159252830 + t * (20.6309250519127594 + t * (5.96507377885482873 + t)))
)
)
)
)
)
)
)
elif x < 40.0:
t = 0.05 * x - 1.0
fd = -(
0.000222755850099679261
+ t
* (
0.00115329637370955141
+ t
* (
0.00290802444182764526
+ t
* (
0.00445062601336644735
+ t
* (
0.00436605816073874863
+ t
* (
0.00260703930450196953
+ t
* (
0.000708204675930321953
+ t
* (
0.0000387281377930932522
+ t * (-1.07635110532534644e-6 + t * 3.45876932703497152e-8)
)
)
)
)
)
)
)
) / (
0.958823449511661332
+ t
* (
7.43981964725756497
+ t
* (
27.2010614567095770
+ t
* (
61.4593095996041459
+ t
* (
94.2239246717949674
+ t
* (
100.959076669111257
+ t * (74.4305466149867624 + t * (35.3635458846636560 + t * (9.43423871492993288 + t)))
)
)
)
)
)
)
else:
w = 1.0 / (x * x)
s = 1.0 - 1600.0 * w
fd = (
factor
/ (sqrt(x) * x * x)
* (
1.0
+ w
* (
114389.966467649306
+ s
* (
109888.157408212290
+ s * (21209.9082486576965 + s * (896.625771647944984 + s * 2.27005555539139462))
)
)
/ (
7803.86522209355981
+ s * (7642.62657004470466 + s * (1584.72834559231955 + s * (86.0988387732258338 + s)))
)
)
)
return fd
@fd_doc_hparams(order=-5 / 2)
@jit(float64(float64), nopython=True)
def fdm5h(x: float) -> float:
factor = -2.0 / 3.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
2.36327180120735470
- ex
* (
55891.891151428516
+ t
* (
10900.4507673238614
+ t * (731.61871282829311 + t * (17.3576202287841595 + t * 0.00111269690244199535))
)
)
/ (8361.6144419623308 + t * (3709.66537451151168 + t * (598.60597628606205 + t * (41.1772384480908825 + t))))
)
elif x < 0.0:
s = -0.5 * x
fd = (
46.7273884990861909
+ s
* (
65.1460094720913271
+ s
* (
-23.3560330148569366
+ s
* (
16.3832748811699383
+ s
* (
-6.19159773777140396
+ s
* (
0.799041627235852088
+ s * (-0.136702158227027886 + s * (0.0592594484852497570 - s * 0.00813537888028213710))
)
)
)
)
)
) / (
166.600811683081505
+ s
* (
-14.3475443594245037
+ s
* (
207.623210767457661
+ s
* (
-16.0651771000195284
+ s
* (
91.9265825199149767
+ s * (-5.91543082309762067 + s * (16.6946542439476255 + s * (-0.714551880206870314 + s)))
)
)
)
)
)
elif x < 2.0:
t = 0.5 * x
fd = (
-(
55.096201117757232
+ t
* (
14.3294333176185961
+ t
* (
8.7057675255700556
+ t
* (
2.85717540287481278
+ t
* (
0.153440611811510950
+ t * (0.111339313377987928 + t * (0.052394928476702755 - t * 0.00308189450941126788))
)
)
)
)
)
/ (
217.840649625384406
+ t
* (
-13.7539547994891024
+ t
* (
264.475360271691825
+ t
* (
-16.3936568402444332
+ t
* (
112.575929318729487
+ t * (-6.5735021347402422 + t * (19.1324946771515682 + t * (-0.89572129073702419 + t)))
)
)
)
)
)
* (x - 1.10894923342229868)
) # care for zero point
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
fd = -(
4.91753613607304828
+ t
* (
20.4599278677959969
+ t
* (
14.7299603843029675
+ t
* (
6.07034924563989291
+ t
* (
2.84749597355562421
+ t
* (
1.42067667370006629
+ t * (0.194025010869198784 + t * (0.0296942990922311325 + t * 0.0000687311336894144097))
)
)
)
)
)
) / (
39.1942895728295075
+ t
* (
102.720679402185990
+ t
* (
163.814228574815100
+ t
* (
160.042363559709804
+ t
* (
112.736106158227044
+ t * (53.5881251269350506 + t * (19.3013585080264391 + t * (4.60099930165203599 + t)))
)
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
fd = -(
2.27928384600199307
+ t
* (
4.20924981693908136
+ t
* (
3.71772880855517432
+ t
* (
3.45602273475008183
+ t
* (
1.35886950179549082
+ t
* (
0.536408369715153184
+ t * (0.133773439887674841 + t * (0.0125358162552656152 - t * 0.000258659892241001720))
)
)
)
)
)
) / (
29.5538409791709387
+ t
* (
102.889070328418668
+ t
* (
171.483234327128671
+ t
* (
174.069526985568413
+ t
* (
121.638517061339171
+ t * (61.2471517251153565 + t * (22.6513730231817499 + t * (5.63585485331701445 + t)))
)
)
)
)
)
elif x < 20.0:
t = 0.1 * x - 1.0
fd = -(
0.138865047109769882
+ t
* (
0.495879942014454703
+ t
* (
0.975841030964797014
+ t
* (
1.19756546203286062
+ t
* (
1.04262591465861665
+ t
* (
0.636815947333856586
+ t
* (
0.267431185397286764
+ t
* (
0.0680541575349993884
+ t * (0.00339987875007404643 - t * 0.0000596473600392971721)
)
)
)
)
)
)
)
) / (
6.10307238045677639
+ t
* (
32.0683844393622288
+ t
* (
81.7709446404731603
+ t
* (
131.907213444156949
+ t
* (
148.807381678792024
+ t
* (
122.126271732385001
+ t * (73.4713378422568434 + t * (31.3211204388750577 + t * (8.63324324695581677 + t)))
)
)
)
)
)
)
elif x < 40.0:
t = 0.05 * x - 1.0
fd = -(
0.294313222065093082
+ t
* (
1.31885924825320721
+ t
* (
2.84488318326892294
+ t
* (
3.47529649317584998
+ t
* (
2.44599317338319296
+ t
* (
0.711594188260176831
+ t * (0.0568475228400688404 + t * (0.000310240841471278615 + t * 7.25941094961112814e-6))
)
)
)
)
)
) / (
38.8567648628867585
+ t
* (
233.714122704911173
+ t
* (
657.094399788013655
+ t
* (
1097.58958431940912
+ t
* (
1157.52830939122072
+ t * (742.054932624983009 + t * (251.716945386093114 + t * (34.4183681886106367 + t)))
)
)
)
)
)
else:
w = 1.0 / (x * x)
s = 1.0 - 1600.0 * w
fd = (
factor
/ (sqrt(x) * x)
* (
1.0
+ w
* (563249.531577994933 + s * (348717.634501234702 + s * (40765.2886083808789 + s * 865.059113888543852)))
/ (
90262.6676538587293
+ s * (56944.6907876652433 + s * (7175.67113380361324 + s * (207.376531486457753 + s)))
)
)
)
return fd
@fd_doc_hparams(order=-3 / 2)
@jit(float64(float64), nopython=True)
def fdm3h(x: float) -> float:
factor = -2.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
-3.54490770181103205
+ ex
* (
82737.595643818605
+ t
* (
18481.5553495836940
+ t * (1272.73919064487495 + t * (26.3420403338352574 - t * 0.00110648970639283347))
)
)
/ (16503.7625405383183 + t * (6422.0552658801394 + t * (890.85389683932154 + t * (51.251447078851450 + t))))
)
elif x < 0.0:
s = -0.5 * x
t = 1.0 - s
fd = -(
946.638483706348559
+ t
* (
76.3328330396778450
+ t
* (
62.7809183134124193
+ t
* (
83.8442376534073219
+ t
* (
23.2285755924515097
+ t
* (
3.21516808559640925
+ t * (1.58754232369392539 + t * (0.687397326417193593 + t * 0.111510355441975495))
)
)
)
)
)
) / (
889.4123665319664
+ s
* (
126.7054690302768
+ s
* (
881.4713137175090
+ s
* (
108.2557767973694
+ s
* (
289.38131234794585
+ s * (27.75902071820822 + s * (34.252606975067480 + s * (1.9592981990370705 + s)))
)
)
)
)
)
elif x < 2.0:
t = 0.5 * x
fd = -(
754.61690882095729
+ t
* (
565.56180911009650
+ t
* (
494.901267018948095
+ t
* (
267.922900418996927
+ t
* (
110.418683240337860
+ t
* (
39.4050164908951420
+ t * (10.8654460206463482 + t * (2.11194887477009033 + t * 0.246843599687496060))
)
)
)
)
)
) / (
560.03894899770103
+ t
* (
70.007586553114572
+ t
* (
582.42052644718871
+ t
* (
56.181678606544951
+ t
* (
205.248662395572799
+ t * (12.5169006932790528 + t * (27.2916998671096202 + t * (0.53299717876883183 + t)))
)
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
fd = -(
526.022770226139287
+ t
* (
631.116211478274904
+ t
* (
516.367876532501329
+ t
* (
267.894697896892166
+ t
* (
91.3331816844847913
+ t
* (
17.5723541971644845
+ t * (1.46434478819185576 + t * (1.29615441010250662 + t * 0.223495452221465265))
)
)
)
)
)
) / (
354.867400305615304
+ t
* (
560.931137013002977
+ t
* (
666.070260050472570
+ t
* (
363.745894096653220
+ t
* (
172.272943258816724
+ t * (23.7751062504377332 + t * (12.5916012142616255 + t * (-0.888604976123420661 + t)))
)
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
fd = -(
18.0110784494455205
+ t
* (
36.1225408181257913
+ t
* (
38.4464752521373310
+ t
* (
24.1477896166966673
+ t
* (
9.27772356782901602
+ t * (2.49074754470533706 + t * (0.163824586249464178 - t * 0.00329391807590771789))
)
)
)
)
) / (
18.8976860386360201
+ t
* (
49.3696375710309920
+ t
* (
60.9273314194720251
+ t * (43.6334649971575003 + t * (20.6568810936423065 + t * (6.11094689399482273 + t)))
)
)
)
elif x < 20.0:
t = 0.1 * x - 1.0
fd = -(
4.10698092142661427
+ t
* (
17.1412152818912658
+ t
* (
32.6347877674122945
+ t
* (
36.6653101837618939
+ t
* (
25.9424894559624544
+ t
* (
11.2179995003884922
+ t * (2.30099511642112478 + t * (0.0928307248942099967 - t * 0.00146397877054988411))
)
)
)
)
)
) / (
6.40341731836622598
+ t
* (
30.1333068545276116
+ t
* (
64.0494725642004179
+ t
* (
80.5635003792282196
+ t * (64.9297873014508805 + t * (33.3013900893183129 + t * (9.61549304470339929 + t)))
)
)
)
)
elif x < 40.0:
t = 0.05 * x - 1.0
fd = -(
95.2141371910496454
+ t
* (
420.050572604265456
+ t
* (
797.778374374075796
+ t
* (
750.378359146985564
+ t
* (
324.818150247463736
+ t
* (
50.3115388695905757
+ t * (0.372431961605507103 + t * (-0.103162211894757911 + t * 0.00191752611445211151))
)
)
)
)
)
) / (
212.232981736099697
+ t
* (
1043.79079070035083
+ t
* (
2224.50099218470684
+ t
* (
2464.84669868672670
+ t * (1392.55318009810070 + t * (346.597189642259199 + t * (22.7314613168652593 - t)))
)
)
)
)
else:
w = 1.0 / (x * x)
s = 1.0 - 1600.0 * w
fd = (
factor
/ sqrt(x)
* (
1.0
+ w
* (12264.3569103180524 + s * (3204.34872454052352 + s * (140.119604748253961 + s * 0.523918919699235590)))
/ (9877.87829948067200 + s * (2644.71979353906092 + s * (128.863768007644572 + s)))
)
)
return fd
@fd_doc_hparams(order=-1 / 2)
@jit(float64(float64), nopython=True)
def fdm1h(x: float) -> float:
"""!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=-1/2
!
! Reference: Fukushima, T. (2015) Appl. Math. Comp., 259, 708-729
! "Precise and fast computation of Fermi-Dirac integral of integer
! and half integer order by piecewise minimax rational approximation"
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
"""
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
1.77245385090551603
- ex
* (
40641.4537510284430
+ t
* (9395.7080940846442 + t * (649.96168315267301 + t * (12.7972295804758967 + t * 0.00153864350767585460)))
)
/ (32427.1884765292940 + t * (11079.9205661274782 + t * (1322.96627001478859 + t * (63.738361029333467 + t))))
)
elif x < 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
272.770092131932696
+ t
* (
30.8845653844682850
+ t
* (
-6.43537632380366113
+ t
* (
14.8747473098217879
+ t
* (
4.86928862842142635
+ t
* (
-1.53265834550673654
+ t * (-1.02698898315597491 + t * (-0.177686820928605932 - t * 0.00377141325509246441))
)
)
)
)
)
) / (
293.075378187667857
+ s
* (
305.818162686270816
+ s
* (
299.962395449297620
+ s
* (
207.640834087494249
+ s * (92.0384803181851755 + s * (37.0164914112791209 + s * (7.88500950271420583 + s)))
)
)
)
)
elif x < 2.0:
t = 0.5 * x
fd = (
3531.50360568243046
+ t
* (
6077.5339658420037
+ t
* (
6199.7700433981326
+ t
* (
4412.78701919567594
+ t
* (
2252.27343092810898
+ t * (811.84098649224085 + t * (191.836401053637121 + t * 23.2881838959183802))
)
)
)
)
) / (
3293.83702584796268
+ t
* (
1528.97474029789098
+ t
* (
2568.48562814986046
+ t
* (
925.64264653555825
+ t * (574.23248354035988 + t * (132.803859320667262 + t * (29.8447166552102115 + t)))
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
fd = (
4060.70753404118265
+ t
* (
10812.7291333052766
+ t
* (
13897.5649482242583
+ t
* (
10628.4749852740029
+ t
* (
5107.70670190679021
+ t * (1540.84330126003381 + t * (284.452720112970331 + t * 29.5214417358484151))
)
)
)
)
) / (
1564.58195612633534
+ t
* (
2825.75172277850406
+ t
* (
3189.16066169981562
+ t
* (
1955.03979069032571
+ t * (828.000333691814748 + t * (181.498111089518376 + t * (32.0352857794803750 + t)))
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
fd = (
1198.41719029557508
+ t
* (
3263.51454554908654
+ t
* (
3874.97588471376487
+ t
* (
2623.13060317199813
+ t
* (
1100.41355637121217
+ t * (267.469532490503605 + t * (25.4207671812718340 + t * 0.389887754234555773))
)
)
)
)
) / (
273.407957792556998
+ t
* (
595.918318952058643
+ t
* (
605.202452261660849
+ t * (343.183302735619981 + t * (122.187622015695729 + t * (20.9016359079855933 + t)))
)
)
)
elif x < 20.0:
t = 0.1 * x - 1.0
fd = (
9446.00169435237637
+ t
* (
36843.4448474028632
+ t
* (
63710.1115419926191
+ t
* (
62985.2197361074768
+ t
* (
37634.5231395700921
+ t * (12810.9898627807754 + t * (1981.56896138920963 + t * 81.4930171897667580))
)
)
)
)
) / (
1500.04697810133666
+ t
* (
5086.91381052794059
+ t
* (
7730.01593747621895
+ t
* (
6640.83376239360596
+ t * (3338.99590300826393 + t * (860.499043886802984 + t * (78.8565824186926692 + t)))
)
)
)
)
elif x < 40.0:
t = 0.05 * x - 1.0
fd = (
22977.9657855367223
+ t
* (
123416.616813887781
+ t
* (
261153.765172355107
+ t
* (
274618.894514095795
+ t
* (
149710.718389924860
+ t * (40129.3371700184546 + t * (4470.46495881415076 + t * 132.684346831002976))
)
)
)
)
) / (
2571.68842525335676
+ t
* (
12521.4982290775358
+ t
* (
23268.1574325055341
+ t
* (
20477.2320119758141
+ t * (8726.52577962268114 + t * (1647.42896896769909 + t * (106.475275142076623 + t)))
)
)
)
)
else:
w = 1.0 / (x * x)
t = 1600.0 * w
fd = (
sqrt(x)
* 2.0
* (
1.0
- w
* (
0.411233516712009968
+ t
* (
0.00110980410034088951
+ t
* (
0.0000113689298990173683
+ t * (2.56931790679436797e-7 + t * (9.97897786755446178e-9 + t * 8.67667698791108582e-10))
)
)
)
)
)
return fd
@fd_doc_iparams(order=0)
@jit(float64(float64), nopython=True)
def fd0h(y: float) -> float:
x = -abs(y)
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
1.0
- ex
* (
22696.2126132366633
+ t
* (5222.0667923565138 + t * (357.623326425354522 + t * (6.9167792879948140 + t * 0.00200096064827815813)))
)
/ (45392.4252264733267 + t * (14539.5980679273792 + t * (1611.36476693109675 + t * (71.072178562726798 + t))))
)
elif x <= 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
159.601717762460980
+ t
* (
23.7193942338278703
+ t
* (
0.377783268730614356
+ t
* (
10.5181677709577503
+ t
* (
3.78181326142271599
+ t
* (
-0.441998676614933572
+ t * (-0.450072959113928254 + t * (-0.0734798777056723512 + t * 0.000915454570009894267))
)
)
)
)
)
) / (
284.26032127745967
+ s
* (
315.2592651624449
+ s
* (
310.2713981221035
+ s
* (
206.21640678892182
+ s * (96.77898293084927 + s * (35.456591489081173 + s * (8.1762315442738975 + s)))
)
)
)
)
if y > 0.0:
fd += y
return fd
@fd_doc_hparams(order=1 / 2)
@jit(float64(float64), nopython=True)
def fd1h(x: float) -> float:
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
0.886226925452758014
- ex
* (
19894.4553386951666
+ t
* (
4509.64329955948557
+ t * (303.461789035142376 + t * (5.7574879114754736 + t * 0.00275088986849762610))
)
)
/ (63493.915041308052 + t * (19070.1178243603945 + t * (1962.19362141235102 + t * (79.250704958640158 + t))))
)
elif x < 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
149.462587768865243
+ t
* (
22.8125889885050154
+ t
* (
-0.629256395534285422
+ t
* (
9.08120441515995244
+ t
* (
3.35357478401835299
+ t
* (
-0.473677696915555805
+ t * (-0.467190913556185953 + t * (-0.0880610317272330793 - t * 0.00262208080491572673))
)
)
)
)
)
) / (
269.94660938022644
+ s
* (
343.6419926336247
+ s
* (
323.9049470901941
+ s
* (
218.89170769294024
+ s * (102.31331350098315 + s * (36.319337289702664 + s * (8.3317401231389461 + s)))
)
)
)
)
elif x < 2.0:
t = 0.5 * x
fd = (
71652.717119215557
+ t
* (
134954.734070223743
+ t
* (
153693.833350315645
+ t
* (
123247.280745703400
+ t
* (
72886.293647930726
+ t
* (
32081.2499422362952
+ t * (10210.9967337762918 + t * (2152.71110381320778 + t * 232.906588165205042))
)
)
)
)
)
) / (
105667.839854298798
+ t
* (
31946.0752989314444
+ t
* (
71158.788776422211
+ t
* (
15650.8990138187414
+ t
* (
13521.8033657783433
+ t * (1646.98258283527892 + t * (618.90691969249409 + t * (-3.36319591755394735 + t)))
)
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
fd = (
23744.8706993314289
+ t
* (
68257.8589855623002
+ t
* (
89327.4467683334597
+ t
* (
62766.3415600442563
+ t
* (
20093.6622609901994
+ t * (-2213.89084119777949 + t * (-3901.66057267577389 - t * 948.642895944858861))
)
)
)
)
) / (
9488.61972919565851
+ t
* (
12514.8125526953073
+ t
* (
9903.44088207450946
+ t
* (
2138.15420910334305
+ t * (-528.394863730838233 + t * (-661.033633995449691 + t * (-51.4481470250962337 + t)))
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
fd = (
(
311337.452661582536
+ t
* (
1.11267074416648198e6
+ t
* (
1.75638628895671735e6
+ t
* (
1.59630855803772449e6
+ t
* (
910818.935456183774
+ t * (326492.733550701245 + t * (65507.2624972852908 + t * 4809.45649527286889))
)
)
)
)
)
/ (
39721.6641625089685
+ t
* (
86424.7529107662431
+ t
* (
88163.7255252151780
+ t
* (
50615.7363511157353
+ t * (17334.9774805008209 + t * (2712.13170809042550 + t * (82.2205828354629102 - t)))
)
)
)
)
* 0.999999999999999877
)
elif x < 20.0:
t = 0.1 * x - 1.0
fd = (
7.26870063003059784e6
+ t
* (
2.79049734854776025e7
+ t
* (
4.42791767759742390e7
+ t
* (
3.63735017512363365e7
+ t * (1.55766342463679795e7 + t * (2.97469357085299505e6 + t * 154516.447031598403))
)
)
)
) / (
340542.544360209743
+ t
* (
805021.468647620047
+ t
* (
759088.235455002605
+ t
* (
304686.671371640343
+ t * (39289.4061400542309 + t * (582.426138126398363 + t * (11.2728194581586028 - t)))
)
)
)
)
elif x < 40.0:
t = 0.05 * x - 1.0
fd = (
4.81449797541963104e6
+ t
* (
1.85162850713127602e7
+ t
* (
2.77630967522574435e7
+ t
* (
2.03275937688070624e7
+ t * (7.41578871589369361e6 + t * (1.21193113596189034e6 + t * 63211.9545144644852))
)
)
)
) / (
80492.7765975237449
+ t
* (
189328.678152654840
+ t
* (
151155.890651482570
+ t * (48146.3242253837259 + t * (5407.08878394180588 + t * (112.195044410775577 - t)))
)
)
)
else:
w = 1.0 / (x * x)
s = 1.0 - 1600.0 * w
fd = (
x
* sqrt(x)
* 0.666666666666666667
* (
1.0
+ w
* (8109.79390744477921 + s * (342.069867454704106 + s * 1.07141702293504595))
/ (6569.98472532829094 + s * (280.706465851683809 + s))
)
)
return fd
@fd_doc_iparams(order=1)
@jit(float64(float64), nopython=True)
def fd2h(y: float) -> float:
x = -abs(y)
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
1.0
- ex
* (
22189.1070807945062
+ t
* (
4915.92700908746777
+ t * (322.901386168881348 + t * (5.9897442965804548 + t * 0.00397641173774375092))
)
)
/ (88756.428323178025 + t * (25002.3197546553836 + t * (2389.06277237306633 + t * (88.376214553692756 + t))))
)
elif x <= 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
145.488167182330098
+ t
* (
251.392824471576922
+ t
* (
56.6537141912783024
+ t
* (
17.9918985363509694
+ t
* (
20.1369115558099802
+ t
* (
7.09659390228556164
+ t * (0.199701180197912643 + t * (-0.403173132925886253 - t * 0.0792966701498222697))
)
)
)
)
)
) / (
606.0757707716040
+ s
* (
374.1806357435014
+ s
* (
252.1367051536344
+ s
* (
27.2746245830016
+ s
* (
-61.57766112137513
+ s * (-53.72117554363975 + s * (-25.678454878692950 + s * (-7.1995819520154718 - s)))
)
)
)
)
)
if y > 0.0:
fd = -fd + 1.64493406684822644 + 0.5 * y * y
return fd
@fd_doc_hparams(order=3 / 2)
@jit(float64(float64), nopython=True)
def fd3h(x: float) -> float:
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
1.32934038817913702
- ex
* (1346.14119046566636 + t * (199.946876779712712 + t * (6.5210149677288048 + t * 0.0108588591982722183)))
/ (5728.3481201778541 + t * (1132.17837281710987 + t * (64.805243148002602 + t)))
)
elif x < 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
631.667081787115831
+ t
* (
504.131655805666135
+ t
* (
113.449065431934917
+ t
* (
56.0939647772947784
+ t
* (
43.3374223200846752
+ t
* (
12.8047010597109577
+ t * (0.219164386586949410 + t * (-0.678659552658390139 - t * 0.126533769309899232))
)
)
)
)
)
) / (
1180.5112183558028
+ s
* (
1101.0159189871135
+ s
* (
864.4448234404281
+ s
* (
392.2018227840790
+ s
* (
89.58093202779063
+ s * (-9.95066218572899 + s * (-17.312068771997626 + s * (-6.4116162917822773 - s)))
)
)
)
)
)
elif x < 2.0:
t = 0.5 * x
fd = (
90122.488639370400
+ t
* (
157095.208147064037
+ t
* (
166879.962599668589
+ t
* (
125708.597728045460
+ t
* (
69968.278213181390
+ t
* (
29035.3292989055404
+ t * (8736.4439472398517 + t * (1747.16784760309227 + t * 180.132410666734053))
)
)
)
)
)
) / (
78176.777123671727
+ t
* (
-1681.44633240543085
+ t
* (
38665.7913035496031
+ t
* (
-2527.29685826087874
+ t
* (
5062.6683078100048
+ t * (-553.21165462054589 + t * (165.395637981775430 + t * (-18.0295465153725544 + t)))
)
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
fd = (
912944.432058014054
+ t
* (
3.28217091334054338e6
+ t
* (
5.59250227196369585e6
+ t
* (
5.76136129685687470e6
+ t
* (
3.84331519034749983e6
+ t * (1.65284168824947710e6 + t * (423452.676670436605 + t * 49835.4127241373113))
)
)
)
)
) / (
164873.145721762182
+ t
* (
257442.511191094986
+ t
* (
225604.160532840884
+ t
* (
99932.1955662320024
+ t * (24761.0878784286761 + t * (1398.26392212830777 + t * (-36.4450237523474167 + t)))
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
fd = (
(
1.88412548327216052e6
+ t
* (
8.08838896259910792e6
+ t
* (
1.56869793001790529e7
+ t
* (
1.79109792599373447e7
+ t
* (
1.31345142328147214e7
+ t
* (
6.29500412046744325e6
+ t * (1.89326213154091054e6 + t * (312372.643127575407 + t * 18814.7420442630170))
)
)
)
)
)
)
/ (
67768.3347951202583
+ t
* (
147635.914444221358
+ t
* (
151908.303165069423
+ t
* (
86671.1222110642970
+ t * (27855.9481608626219 + t * (3833.22697473114940 + t * (98.3384567064269554 - t)))
)
)
)
)
* 0.999999999999999876
)
elif x < 20.0:
t = 0.1 * x - 1.0
fd = (
(
1.59656593348660977e9
+ t
* (
7.32769737561517060e9
+ t
* (
1.42662658588280191e10
+ t
* (
1.51238422045169918e10
+ t
* (
9.27233604548095476e9
+ t * (3.18834513406577423e9 + t * (5.36061988605886123e8 + t * 3.03619219668246382e7))
)
)
)
)
)
/ (
1.18906980815759995e7
+ t
* (
2.62209219322122975e7
+ t
* (
2.28143701746618729e7
+ t
* (
8.57156701742181874e6
+ t * (1.13860063870524239e6 + t * (27091.7884208687379 + t * (-275.664733379090447 + t)))
)
)
)
)
* 0.999999999999999829
)
elif x < 40.0:
t = 0.05 * x - 1.0
fd = (
(
2.60437581212904589e8
+ t
* (
1.08771546307370080e9
+ t
* (
1.81531350939088943e9
+ t
* (
1.52833764636304939e9
+ t * (6.70684451492750149e8 + t * (1.40870639531414149e8 + t * 1.04957900377463854e7))
)
)
)
)
/ (
358448.871166784200
+ t
* (
611808.419702466190
+ t
* (
326307.561591723775
+ t * (58407.9904827573816 + t * (2049.50040323021794 + t * (-39.8767861209088081 + t)))
)
)
)
* 0.999999999999999828
)
else:
w = 1.0 / (x * x)
s = 1.0 - 1600.0 * w
fd = (
x
* x
* sqrt(x)
* 2.0
/ 5.0
* (
1.0
+ w
* (
6.16739021212286242
+ s
* (
0.00111530123694574981
+ s * (-2.79156524536560815e-6 + s * (2.95571462110856359e-8 - s * 6.70917556862133933e-10))
)
)
)
)
return fd
@fd_doc_iparams(order=2)
@jit(float64(float64), nopython=True)
def fd4h(y: float) -> float:
x = -abs(y)
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
2.0
- ex
* (1914.06748184935743 + t * (273.085756700981399 + t * (8.5861610217850095 + t * 0.0161890243763741414)))
/ (7656.2699273974454 + t * (1399.35442210906621 + t * (72.929152915475392 + t)))
)
elif x <= 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
2711.49678259128843
+ t
* (
1299.85460914884154
+ t
* (
222.606134197895041
+ t
* (
172.881855215582924
+ t
* (
112.951038040682055
+ t
* (
24.0376482128898634
+ t * (-2.68393549333878715 + t * (-2.14077421411719935 - t * 0.326188299771397236))
)
)
)
)
)
) / (
2517.1726659917047
+ s
* (
3038.7689794575778
+ s
* (
2541.7823512372631
+ s
* (
1428.0589853413436
+ s
* (
531.62378035996132
+ s * (122.54595216479181 + s * (8.395768655115050 + s * (-3.9142702096919080 - s)))
)
)
)
)
)
if y > 0.0:
fd = fd + y * (3.28986813369645287 + 0.333333333333333333 * y * y)
return fd
@fd_doc_hparams(order=5 / 2)
@jit(float64(float64), nopython=True)
def fd5h(x: float) -> float:
factor = 2 / 7.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
3.32335097044784255
- ex
* (3004.41138112148735 + t * (409.582975848055860 + t * (12.3422465543050559 + t * 0.0252600128832787650)))
/ (10227.9400759349389 + t * (1729.21912667445507 + t * (82.075224736463695 + t)))
)
elif x < 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
3273.17630052701392
+ t
* (
1698.67676090446982
+ t
* (
498.075924016264615
+ t
* (
286.399893838802205
+ t
* (
153.449411299292925
+ t * (46.3172989461474687 + t * (7.47269262560340991 + t * 0.512665223813025153))
)
)
)
)
) / (
1934.7654167995427
+ s
* (
2174.8582387970533
+ s
* (
1803.7716559637348
+ s
* (
970.4325125469826
+ s * (368.29988963760778 + s * (94.85059993048974 + s * (14.845914179579222 + s)))
)
)
)
)
elif x < 2.0:
t = 0.5 * x
s = 1.0 - t
fd = (
46746.0667679244140
+ t
* (
67556.152887188131
+ t
* (
58097.414724032516
+ t
* (
35613.1718933770642
+ t
* (
15859.9235610267359
+ t * (5040.9166211327297 + t * (1064.44778209849372 + t * 117.741363279815688))
)
)
)
)
) / (
13126.266942469915
+ s
* (
-1313.6945396119670
+ s
* (
2636.3767383046264
+ s
* (
453.2433336267639
+ s * (201.68980019364685 + s * (50.35944043583614 + s * (9.319958361334161 + s)))
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
fd = (
2.96446396463084245e6
+ t
* (
1.13985063437515144e7
+ t
* (
2.03551842990364919e7
+ t
* (
2.18469168036723457e7
+ t
* (
1.52194204967854489e7
+ t * (6.89564335059367232e6 + t * (1.88985082880325956e6 + t * 243639.263893338434))
)
)
)
)
) / (
169113.644493386066
+ t
* (
249598.556346488591
+ t
* (
162470.264413303189
+ t
* (
48830.6689634566512
+ t * (4080.41491690869679 + t * (-227.369223251313984 + t * (19.1547319985914212 - t)))
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
fd = (
(
5.98759698946959779e7
+ t
* (
2.87159410448355925e8
+ t
* (
5.99247230055005978e8
+ t
* (
7.05261579500244084e8
+ t
* (
5.04521860843522470e8
+ t * (2.18416589668256432e8 + t * (5.25509373690209805e7 + t * 5.33490422476777307e6))
)
)
)
)
)
/ (
469653.962099017877
+ t
* (
972160.483958338969
+ t
* (
696829.387839904810
+ t
* (
200021.997093042021
+ t * (12604.5784254924892 + t * (-464.256227099766735 + t * (26.4866799762414562 - t)))
)
)
)
)
* 0.999999999999999786
)
elif x < 20.0:
t = 0.1 * x - 1.0
fd = (
(
1.11471644837527339e9
+ t
* (
5.74225610798369390e9
+ t
* (
1.26551657157853213e10
+ t
* (
1.53960861052359422e10
+ t
* (
1.10946777927365321e10
+ t * (4.68865671930554474e9 + t * (1.05763202432846973e9 + t * 9.49124183370767767e7))
)
)
)
)
)
/ (
1.07734938834750844e6
+ t
* (
2.05459707311873616e6
+ t
* (
1.39824422108531691e6
+ t
* (
347716.997197363113
+ t * (17894.5194245484999 + t * (-553.162195184268593 + t * (28.1090136251865326 - t)))
)
)
)
)
* 0.999999999999999759
)
elif x < 40.0:
t = 0.05 * x - 1.0
fd = (
(
8.02670881104191218e9
+ t
* (
3.67306963017003546e10
+ t
* (
6.83085243661572356e10
+ t
* (
6.55107149365528043e10
+ t * (3.37406319128317264e10 + t * (8.67814222875818408e9 + t * 8.43844503352450216e8))
)
)
)
)
/ (
757905.984443885971
+ t
* (
868424.806294231233
+ t
* (
260670.917865642513
+ t * (14550.0472712579662 + t * (-476.120164041067762 + t * (25.8288614974100332 - t)))
)
)
)
* 0.999999999999999820
)
else:
w = 1.0 / (x * x)
t = 1600.0 * w
fd = (
x
* x
* x
* sqrt(x)
* factor
* (
1.0
+ w
* (
14.3931730849220041
+ t
* (
0.00776862867834285253
+ t * (3.78966458769690333e-6 + t * (2.09248095155530095e-8 + t * 3.52097438532254351e-10))
)
)
)
)
return fd
@fd_doc_iparams(order=3)
@jit(float64(float64), nopython=True)
def fd6h(y: float) -> float:
x = -abs(y)
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
6.0
- ex
* (5121.6401850302408 + t * (664.28706260743472 + t * (19.0856927562699544 + t * 0.0410982603688952131)))
/ (13657.7071600806539 + t * (2136.54222460571183 + t * (92.376788603062645 + t)))
)
elif x <= 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
7881.24597452900838
+ t
* (
4323.07526636309661
+ t
* (
1260.13125873282465
+ t
* (
653.359212389160499
+ t
* (
354.630774329461644
+ t * (113.373708671587772 + t * (19.9559488532742796 + t * 1.59407954898394322))
)
)
)
)
) / (
2570.7250703533430
+ s
* (
2972.7443644211129
+ s
* (
2393.9995533270879
+ s
* (
1259.0724833462608
+ s * (459.86413596901097 + s * (112.60906419590854 + s * (16.468882811659000 + s)))
)
)
)
)
if y > 0.0:
y2 = y * y
fd = -fd + 11.3643939539669510 + y2 * (4.93480220054467931 + y2 * 0.25)
return fd
@fd_doc_hparams(order=7 / 2)
@jit(float64(float64), nopython=True)
def fd7h(x: float) -> float:
factor = 2.0 / 9.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
11.6317283965674489
- ex
* (9371.8868848378853 + t * (1152.33581795824468 + t * (31.4225568934398918 + t * 0.069508514261022902)))
/ (18231.3053891121107 + t * (2639.60073595942887 + t * (103.985056337236794 + t)))
)
elif x < 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
12412.9547163940264
+ t
* (
6557.39898920194282
+ t
* (
1745.01251050433134
+ t
* (
883.903030170203178
+ t
* (
485.680209936338180
+ t
* (
150.395267497504100
+ t * (23.6841796284116177 + t * (1.07158463967401877 - t * 0.157435982722068683))
)
)
)
)
)
) / (
1990.3886647182248
+ s
* (
2450.5836860178091
+ s
* (
1959.3508463399357
+ s
* (
1039.6469110978007
+ s * (381.70374149314543 + s * (94.69566180593596 + s * (14.401931782580504 + s)))
)
)
)
)
elif x < 2.0:
t = 0.5 * x
s = 1.0 - t
fd = (
138724.961458493789
+ t
* (
169546.564556791848
+ t
* (
125336.240767414531
+ t
* (
68469.193938175092
+ t
* (
27692.3811560029092
+ t * (8054.8076265209314 + t * (1569.98144698357220 + t * 162.669853437143155))
)
)
)
)
) / (
7624.675802050973
+ s
* (
2096.5891309081106
+ s
* (
1735.7563869500611
+ s
* (
662.0787128606726
+ s * (218.91333929478294 + s * (55.27329667089387 + s * (9.904579892966869 + s)))
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
fd = (
2.20638302461716245e7
+ t
* (
7.19733683266519997e7
+ t
* (
1.11705692645922359e8
+ t
* (
1.06328515554623708e8
+ t
* (
6.68969329369869135e7
+ t * (2.78374130823093789e7 + t * (7.11915156250133512e6 + t * 870011.051674696596))
)
)
)
)
) / (
311792.108500755594
+ t
* (
206109.613235058227
+ t
* (
81814.2443130826819
+ t
* (
3513.24062956494342
+ t * (368.155702298381570 + t * (-106.091905490453452 + t * (15.0961925962965432 - t)))
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
fd = (
(
3.0709767203984838e9
+ t
* (
1.3145330402583295e10
+ t
* (
2.5004908399964440e10
+ t
* (
2.7299220224056906e10
+ t
* (
1.8396311263825673e10
+ t * (7.6098692716640943e9 + t * (1.7746200575457532e9 + t * 1.7740961200347283e8))
)
)
)
)
)
/ (
4.46811391083410494e6
+ t
* (
4.62189039302323844e6
+ t
* (
1.60932873023944955e6
+ t
* (
113458.965071032695
+ t * (-4751.36769214919833 + t * (351.424157693641750 + t * (-24.0392624870048634 + t)))
)
)
)
)
* 0.999999999999999779
)
elif x < 20.0:
t = 0.1 * x - 1.0
fd = (
(
4.91615192873221833e9
+ t
* (
2.30941143669878184e10
+ t
* (
4.52380364776557202e10
+ t
* (
4.68413409944225048e10
+ t
* (
2.65134543816397678e10
+ t * (7.35961621586797944e9 + t * (5.52546586238663078e8 - t * 7.67154234216335832e7))
)
)
)
)
)
/ (
550749.688306094769
+ t
* (
352810.964448809157
+ t
* (
12162.0407442978353
+ t * (-4414.91208383702609 + t * (288.821012224214490 + t * (-21.7191096434783647 + t)))
)
)
)
* 0.999999999999999647
)
elif x < 40.0:
t = 0.05 * x - 1.0
fd = (
(
1.23356074447405045e12
+ t
* (
6.8241092859415704e12
+ t
* (
1.5939237639946117e13
+ t
* (
2.0277830436278028e13
+ t
* (
1.5074920317497838e13
+ t * (6.4878414656043323e12 + t * (1.4761927839220043e12 + t * 1.3410778049868313e11))
)
)
)
)
)
/ (
7.2811112110966361e6
+ t
* (
8.4187326609915276e6
+ t
* (
2.59806064027658564e6
+ t
* (
153974.958146639672
+ t * (-5666.5977917092119 + t * (387.205694265311079 + t * (-25.1500172248070365 + t)))
)
)
)
)
* 0.999999999999999750
)
else:
w = 1.0 / (x * x)
t = 1600.0 * w
fd = (
x
* x
* x
* x
* sqrt(x)
* factor
* (
1.0
+ w
* (
25.9077115528595532
+ t
* (
0.0699176581158797143
+ t * (-0.0000113689722343055157 + t * (-2.69205155161558844e-8 - t * 2.81838487282327867e-10))
)
)
)
)
return fd
@fd_doc_iparams(order=4)
@jit(float64(float64), nopython=True)
def fd8h(y: float) -> float:
x = -abs(y)
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
24.0
- ex
* (18247.2542465629199 + t * (2120.56302902849207 + t * (54.659507299984584 + t * 0.121876197098273914)))
/ (24329.6723287505678 + t * (3261.01909656925521 + t * (117.071576489684022 + t)))
)
elif x <= 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
33093.9102025608137
+ t
* (
19031.4783357798975
+ t
* (
5431.08357152245897
+ t
* (
2393.94262931609398
+ t
* (
1268.40017257978070
+ t * (418.662172475927160 + t * (77.4108960876508190 + t * 6.67374311450268357))
)
)
)
)
) / (
2645.4885670047153
+ s
* (
3236.2237702166948
+ s
* (
2500.5977847175497
+ s
* (
1278.0109577275445
+ s * (448.99020896813485 + s * (105.86020755838874 + s * (15.216887271751039 + s)))
)
)
)
)
if y > 0.0:
y2 = y * y
fd = fd + y * (45.4575758158678040 + y2 * (6.57973626739290575 + y2 * 0.2))
return fd
@fd_doc_hparams(order=9 / 2)
@jit(float64(float64), nopython=True)
def fd9h(x: float) -> float:
factor = 2.0 / 11.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
52.3427777845535202
- ex
* (37544.8291647986266 + t * (4113.48125371067712 + t * (99.866490020209337 + t * 0.221035794126711913)))
/ (32460.7344732907336 + t * (4028.81223496037193 + t * (131.831032502212893 + t)))
)
elif x < 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
69576.0137991453185
+ t
* (
40945.8214269928341
+ t
* (
11636.1644381151941
+ t
* (
4627.86250549447612
+ t
* (
2350.13392862571985
+ t * (772.917515779084181 + t * (143.282768049011730 + t * 12.4161720203166775))
)
)
)
)
) / (
2535.8419026162019
+ s
* (
3176.5677917330904
+ s
* (
2403.1400307314004
+ s
* (
1205.9132705211921
+ s * (416.75320914741449 + s * (97.44802030226960 + s * (14.192502951043942 + s)))
)
)
)
)
elif x < 2.0:
t = 0.5 * x
s = 1.0 - t
fd = (
514822.53904173208
+ t
* (
522785.81416122782
+ t
* (
325144.232801546845
+ t
* (
158667.557186950475
+ t
* (
59144.454873935605
+ t * (15870.4677825005999 + t * (2843.37449801975925 + t * 272.600692278110437))
)
)
)
)
) / (
4614.934746888508
+ s
* (
2755.8568288489921
+ s
* (
1644.4943294235255
+ s
* (
712.2181577760027
+ s * (238.45583147179604 + s * (59.96344339125020 + s * (10.461476374827578 + s)))
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
fd = (
3.42730803487596668e7
+ t
* (
9.75806658972896811e7
+ t
* (
1.34864192616165924e8
+ t
* (
1.16686466196565430e8
+ t
* (
6.80744250678140765e7
+ t * (2.67765233037235380e7 + t * (6.59736845425507999e6 + t * 792820.676193592762))
)
)
)
)
) / (
98882.7648048876841
+ t
* (
8989.62788188007021
+ t
* (
9880.17356654342362
+ t
* (
-2706.34081341857522
+ t * (689.741654284866509 + t * (-130.873077409052406 + t * (16.2961910923971688 - t)))
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
fd = (
(
4.89775289337907328e9
+ t
* (
2.18308560190839931e10
+ t
* (
4.37929196886803725e10
+ t
* (
5.12376841913497741e10
+ t
* (
3.77548677711700977e10
+ t * (1.75249365659518690e10 + t * (4.74974437725451705e9 + t * 5.81264585063986499e8))
)
)
)
)
)
/ (
1.16603582458148748e6
+ t
* (
904387.524377664561
+ t
* (
128592.404626994618
+ t
* (
-10482.2282242040897
+ t * (1464.66129358886087 + t * (-193.730497094764856 + t * (19.0689920747890676 - t)))
)
)
)
)
* 0.999999999999999823
)
elif x < 20.0:
t = 0.1 * x - 1.0
fd = (
(
2.2676209606423816e13
+ t
* (
1.24113302112340516e14
+ t
* (
2.9428117821474437e14
+ t
* (
3.9045588584587502e14
+ t
* (
3.1161725148614816e14
+ t * (1.4855808500445768e14 + t * (3.8726601091729722e13 + t * 4.1673615587969372e12))
)
)
)
)
)
/ (
2.72286381952640596e8
+ t
* (
1.76997321474048891e8
+ t
* (
1.58695998922133285e7
+ t
* (
-757839.84310180434
+ t * (61065.315196529561 + t * (-4351.64276565798757 + t * (174.373773100542411 + t)))
)
)
)
)
* 0.999999999999999697
)
elif x < 40.0:
t = 0.05 * x - 1.0
fd = (
(
1.65627579413599954e13
+ t
* (
9.8868092265530109e13
+ t
* (
2.51632614382698224e14
+ t
* (
3.53003554017750602e14
+ t
* (
2.93707971643967604e14
+ t * (1.44188847117098389e14 + t * (3.83794461491062424e13 + t * 4.22247796386900843e12))
)
)
)
)
)
/ (
5.7640009494074862e6
+ t
* (
3.82118971474499980e6
+ t
* (
375185.320814949919
+ t
* (
-21189.7277758964890
+ t * (2225.13154291024219 + t * (-240.494104779578652 + t * (20.7062411414387229 - t)))
)
)
)
)
* 0.999999999999999697
)
else:
w = 1.0 / (x * x)
t = 1600.0 * w
fd = (
x
* x
* x
* x
* x
* sqrt(x)
* factor
* (
1.0
+ w
* (
40.7121181544936151
+ t
* (
0.256364746421960192
+ t * (0.000125058641885852090 + t * (5.92447825215879480e-8 + t * 3.38237202703194112e-10))
)
)
)
)
return fd
@fd_doc_iparams(order=5)
@jit(float64(float64), nopython=True)
def fd10h(y: float) -> float:
x = -abs(y)
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
120.0
- ex
* (81190.938912603315 + t * (8368.4990332049831 + t * (190.753830813273698 + t * 0.413800735538960261)))
/ (43301.8340867217726 + t * (4977.68099709243407 + t * (148.484432990224213 + t)))
)
elif x <= 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
159651.547840031111
+ t
* (
96307.2005742063042
+ t
* (
26923.5693307648389
+ t
* (
9274.54751848368696
+ t
* (
4445.76333033698006
+ t * (1461.45267097859337 + t * (272.164427980501432 + t * 23.6526046298891619))
)
)
)
)
) / (
2522.7839609396783
+ s
* (
3244.5527999477567
+ s
* (
2403.0532924519756
+ s
* (
1176.7202478443275
+ s * (397.7596246691212 + s * (91.84661231161838 + s * (13.491911254479298 + s)))
)
)
)
)
if y > 0.0:
y2 = y * y
fd = (
-fd
+ 236.532261911384425
+ y2 * (113.643939539669510 + y2 * (8.22467033424113218 + y2 * 0.166666666666666667))
)
return fd
@fd_doc_hparams(order=11 / 2)
@jit(float64(float64), nopython=True)
def fd11h(x: float) -> float:
factor = 2.0 / 13.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
287.885277815044361
- ex
* (183706.438585548927 + t * (17781.2198813711389 + t * (379.468701417596835 + t * 0.79823396490661890)))
/ (57756.370490854840 + t * (6150.6214356632731 + t * (167.282750039076672 + t)))
)
elif x < 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
399729.728975269086
+ t
* (
245048.290430557943
+ t
* (
64937.4257386015642
+ t
* (
17956.2886271385592
+ t
* (
8102.86054746075370
+ t * (2739.56619011305740 + t * (522.197272134368319 + t * 45.9656786040236010))
)
)
)
)
) / (
2594.1469579416305
+ s
* (
3458.1847294224010
+ s
* (
2522.2748339106920
+ s
* (
1201.9024264852580
+ s * (394.94834682669296 + s * (89.38046041124415 + s * (13.110503750759739 + s)))
)
)
)
)
elif x < 2.0:
t = 0.5 * x
s = 1.0 - t
fd = (
2.40644550206125935e6
+ t
* (
2.04903536564989627e6
+ t
* (
1.06804765545666824e6
+ t
* (
478219.980682203232
+ t
* (
171969.059741023943
+ t * (43927.0876916260199 + t * (7256.0914773211089 + t * 626.34457422853533))
)
)
)
)
) / (
3109.1225672233219
+ s
* (
2636.2757909353112
+ s
* (
1626.1798939585643
+ s
* (
743.1451899485297
+ s * (255.41213934754826 + s * (64.38943570079975 + s * (10.994511861529934 + s)))
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
fd = (
7.70852968624936390e8
+ t
* (
2.17991233625412353e9
+ t
* (
3.06662111607913460e9
+ t
* (
2.78373598102826860e9
+ t
* (
1.77557231061500545e9
+ t
* (
8.12988256861068453e8
+ t * (2.60886497846886753e8 + t * (5.37757509332685362e7 + t * 5.49029147083874521e6))
)
)
)
)
)
) / (
384975.653810726401
+ t
* (
-10859.8352547471390
+ t
* (
47234.6496298591730
+ t
* (
-15689.8880335197160
+ t
* (
4580.68161072037214
+ t * (-1041.20834072320464 + t * (172.979919201040627 + t * (-18.7593019718622089 + t)))
)
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
fd = (
(
1.71681959469725307e10
+ t
* (
7.11557910478070659e10
+ t
* (
1.35103963154256158e11
+ t
* (
1.52010553996923224e11
+ t
* (
1.09288906897533869e11
+ t * (5.01817283608162539e10 + t * (1.36404590384355901e10 + t * 1.69961340800732270e9))
)
)
)
)
)
/ (
600169.964027912926
+ t
* (
63989.8050690151317
+ t
* (
3306.49827047230707
+ t
* (
-2248.31637427985521
+ t * (649.075943625483938 + t * (-126.914632653258205 + t * (16.0289311039308149 - t)))
)
)
)
)
* 0.999999999999999735
)
elif x < 20.0:
t = 0.1 * x - 1.0
fd = (
(
1.05700178515928121e12
+ t
* (
6.12117465680051512e12
+ t
* (
1.55582880442596747e13
+ t
* (
2.24736526624247277e13
+ t
* (
1.99069164118763497e13
+ t * (1.08056362404911750e13 + t * (3.32622572257609455e12 + t * 4.47682895372084249e11))
)
)
)
)
)
/ (
1.27888301811042482e6
+ t
* (
318603.811469994946
+ t
* (
-33834.0645458632405
+ t
* (
6000.96391500182524
+ t * (-1070.68701259288375 + t * (161.612080255393227 + t * (-17.5331213409665917 + t)))
)
)
)
)
* 0.999999999999999653
)
elif x < 40.0:
t = 0.05 * x - 1.0
fd = (
(
1.06038260567473087e14
+ t
* (
6.8238792169646891e14
+ t
* (
1.88917697801412952e15
+ t
* (
2.91401892625408501e15
+ t
* (
2.70150682400302739e15
+ t * (1.50286942831009197e15 + t * (4.63446459437813772e14 + t * 6.0885766731916815e13))
)
)
)
)
)
/ (
2.08723140369481691e6
+ t
* (
445823.690247689681
+ t
* (
-42239.1125174643989
+ t
* (
6891.8137509933955
+ t * (-1158.79785029794968 + t * (168.017431473136672 + t * (-17.7743748531881933 + t)))
)
)
)
)
* 0.999999999999999637
)
else:
w = 1.0 / (x * x)
t = 1600.0 * w
fd = (
x
* x
* x
* x
* x
* x
* sqrt(x)
* factor
* (
1.0
+ w
* (
58.8063928898240827
+ t
* (
0.666548340699108862
+ t * (0.00162576228417182612 + t * (-2.56764926438640722e-7 - t * 6.18964665583859548e-10))
)
)
)
)
return fd
@fd_doc_iparams(order=6)
@jit(float64(float64), nopython=True)
def fd12h(y: float) -> float:
x = -abs(y)
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
720.0
- ex
* (433290.557356514403 + t * (39323.1678283287030 + t * (783.71684376947655 + t * 1.58412947146158337)))
/ (77029.432418935897 + t * (7600.9245809611507 + t * (188.511069473956679 + t)))
)
elif x <= 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
1.06809887479312876e6
+ t
* (
651074.246191348755
+ t
* (
152197.899924352192
+ t
* (
27269.0203707062592
+ t
* (
11937.9088600476726
+ t * (4659.76125900467198 + t * (964.100791156809939 + t * 88.5841245838029230))
)
)
)
)
) / (
2681.3731718905701
+ s
* (
3764.9048490469408
+ s
* (
2739.9504946219358
+ s
* (
1276.0933863294022
+ s * (406.93880737411049 + s * (89.70143752313466 + s * (12.997281214703279 + s)))
)
)
)
)
if y > 0.0:
y2 = y * y
fd = fd + y * (
1419.19357146830655 + y2 * (227.287879079339020 + y2 * (9.86960440108935862 + y2 * 0.142857142857142857))
)
return fd
@fd_doc_hparams(order=13 / 2)
@jit(float64(float64), nopython=True)
def fd13h(x: float) -> float:
factor = 2.0 / 15.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
s = 1.0 - t
fd = ex * (
1871.25430579778835
- ex
* (
10.2714824358192133
+ s
* (
0.0648393897767320738
+ s
* (
0.000971410002224865156
+ s * (0.0000232470542763426993 + s * (7.41058598210549775e-7 + s * 3.37424008840868627e-8))
)
)
)
)
elif x < 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
2.94504307654357499e6
+ t
* (
1.72142318914178000e6
+ t
* (
302664.670131641659
+ t
* (
-4392.91160671565775
+ t
* (
1886.14265293074511
+ t * (5292.33536016064607 + t * (1532.39748567462704 + t * 160.488909859013319))
)
)
)
)
) / (
2672.0262186872157
+ s
* (
4050.0534640335863
+ s
* (
2998.5673483066873
+ s
* (
1381.1929813352472
+ s * (429.69401904358381 + s * (92.16501861997658 + s * (13.091741877025778 + s)))
)
)
)
)
elif x < 2.0:
t = 0.5 * x
s = 1.0 - t
fd = (
1.45543624525939873e7
+ t
* (
1.09901378809591560e7
+ t
* (
5.0278309600981774e6
+ t
* (
2.17920053476107569e6
+ t
* (
801203.84933034017
+ t * (204535.562039585858 + t * (32251.6713113453623 + t * 2520.49966644549104))
)
)
)
)
) / (
2525.1970865635211
+ s
* (
2519.0633694107075
+ s
* (
1640.0689226392724
+ s
* (
779.5069098915619
+ s * (273.74080946976357 + s * (69.07956623994342 + s * (11.541634429054885 + s)))
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
s = 1.0 - t
fd = (
1.93509493140212933e9
+ t
* (
4.73031979570464235e9
+ t
* (
5.87850060979944364e9
+ t
* (
4.82135265081185025e9
+ t
* (
2.83819099819062720e9
+ t
* (
1.22325945046779759e9
+ t * (3.76487745749965788e8 + t * (7.58303164598034022e7 + t * 7.71973955404921694e6))
)
)
)
)
)
) / (
101717.15224358499
+ s
* (
25406.290690742596
+ s
* (
11612.352340539603
+ s
* (
4174.0221825716411
+ s
* (
1281.7558184421595
+ s * (331.85386698947176 + s * (69.634528227703809 + s * (10.868655326785347 + s)))
)
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
fd = (
(
5.2085062280477356e11
+ t
* (
2.2833424052547072e12
+ t
* (
4.6942573239060036e12
+ t
* (
5.8958834991587034e12
+ t
* (
4.9360061465395771e12
+ t
* (
2.8153080727145281e12
+ t
* (1.06730054566227826e12 + t * (2.4600629203104023e11 + t * 2.6466259504577480e10))
)
)
)
)
)
)
/ (
2.42061644583385800e6
+ t
* (
153096.878206723536
+ t
* (
38844.6431417054012
+ t
* (
-16885.0152281400560
+ t
* (
4865.05644899569952
+ t * (-1067.44742924509046 + t * (173.008767319757249 + t * (-18.6099093504480426 + t)))
)
)
)
)
)
* 0.999999999999999770
)
elif x < 20.0:
t = 0.1 * x - 1.0
fd = (
(
7.7494086475291441e13
+ t
* (
5.0061227151520991e14
+ t
* (
1.45146838868293917e15
+ t
* (
2.46379854838342831e15
+ t
* (
2.67523099025981612e15
+ t
* (
1.90126580903196655e15
+ t
* (8.6325468139568610e14 + t * (2.28901218799706877e14 + t * 2.71416555283189345e13))
)
)
)
)
)
)
/ (
8.9593891719102305e6
+ t
* (
2.23010123362691961e6
+ t
* (
-239561.734722929070
+ t
* (
43938.8978835860601
+ t
* (
-8408.3951285804076
+ t * (1448.44898008864648 + t * (-201.658017368512975 + t * (19.7129633131788616 - t)))
)
)
)
)
)
* 0.999999999999999586
)
elif x < 40.0:
t = 0.05 * x - 1.0
fd = (
(
1.22298919425291664e16
+ t
* (
8.9726625177108898e16
+ t
* (
2.89429751278928300e17
+ t
* (
5.3577901105140076e17
+ t
* (
6.2206348573306568e17
+ t
* (
4.63426239396617756e17
+ t
* (2.16065376561259960e17 + t * (5.7540345485220279e16 + t * 6.6835199866465357e15))
)
)
)
)
)
)
/ (
1.31871866498468464e7
+ t
* (
2.83885323472743888e6
+ t
* (
-275298.111415499841
+ t
* (
47078.0598194149218
+ t
* (
-8612.3789225919783
+ t * (1446.29041011108831 + t * (-199.315359842233309 + t * (19.5132779151984627 - t)))
)
)
)
)
)
* 0.999999999999999580
)
else:
w = 1.0 / (x * x)
t = 1600.0 * w
fd = (
x
* x
* x
* x
* x
* x
* x
* sqrt(x)
* factor
* (
1.0
+ w
* (
80.1905357588510589
+ t
* (
1.42831787292433216
+ t * (0.00812881145350198890 + t * (3.85165680520411590e-6 + t * 1.83598087001386478e-9))
)
)
)
)
return fd
@fd_doc_iparams(order=7)
@jit(float64(float64), nopython=True)
def fd14h(y: float) -> float:
x = -abs(y)
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
s = 1.0 - t
fd = ex * (
5040.0
- ex
* (
19.5849162780581217
+ s
* (
0.101236264600248646
+ s
* (
0.00131827373096852460
+ s * (0.0000283211236200358235 + s * (8.27819847443616991e-7 + s * 3.46664379401620231e-8))
)
)
)
)
elif x <= 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
8.05214238081846197e6
+ t
* (
4.23358827341165564e6
+ t
* (
262627.912342619479
+ t
* (
-323968.614686001584
+ t
* (
-98175.8116505823446
+ t * (-7562.84276896647246 + t * (1252.57372771721279 + t * 240.247242901336449))
)
)
)
)
) / (
2413.8835947192718
+ s
* (
4139.8252477104615
+ s
* (
3204.5076382872745
+ s
* (
1486.4720567970787
+ s * (456.50435752749115 + s * (95.80622002968487 + s * (13.317952574698873 + s)))
)
)
)
)
if y > 0.0:
y2 = y * y
fd = (
-fd
+ 10042.0286586746908
+ y2 * (4967.17750013907292 + y2 * (397.753788388843285 + y2 * (11.5145384679375851 + y2 * 0.125)))
)
return fd
@fd_doc_hparams(order=15 / 2)
@jit(float64(float64), nopython=True)
def fd15h(x: float) -> float:
factor = 2.0 / 17.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
14034.4072934834126
- ex
* (171704.258234536602 + t * (6049.5894628442878 + t * 18.1763510883857458))
/ (4429.36992795387362 + t * (175.155804618692908 + t))
)
elif x < 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
2.06790092766575034e7
+ t
* (
8.26618273142399596e6
+ t
* (
-2.08115753677480168e6
+ t
* (
-2.17301177649999531e6
+ t
* (
-641552.116750411401
+ t * (-91807.9577537648588 + t * (-5144.12042361256270 + t * 130.857943513312398))
)
)
)
)
) / (
1711.2966803109211
+ s
* (
3804.4823677048040
+ s
* (
3232.4717528599742
+ s
* (
1550.7278916841484
+ s * (478.50436506151174 + s * (99.39616883636118 + s * (13.585756243903013 + s)))
)
)
)
)
elif x < 2.0:
t = 0.5 * x
s = 1.0 - t
fd = (
1.14403769376860621e8
+ t
* (
8.4972567501965974e7
+ t
* (
3.81753609315202610e7
+ t
* (
1.65881819811148351e7
+ t
* (
6.1705051839981204e6
+ t * (1.58079386467916017e6 + t * (246507.752966115175 + t * 18602.2895319461278))
)
)
)
)
) / (
2572.2942565752562
+ s
* (
2659.5730458609857
+ s
* (
1732.1057976186518
+ s
* (
829.4655494427400
+ s * (293.34263153462293 + s * (73.73539965269387 + s * (12.058936149293461 + s)))
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
s = 1.0 - t
fd = (
6.26769725573910255e9
+ t
* (
1.32711699403202671e10
+ t
* (
1.45560142717062095e10
+ t
* (
1.07860795716164403e10
+ t
* (
5.86969120046064229e9
+ t
* (
2.38796511791514625e9
+ t * (7.07324240427174132e8 + t * (1.39791587840460999e8 + t * 1.42694865767192645e7))
)
)
)
)
)
) / (
30361.411562018280
+ s
* (
17707.813587539272
+ s
* (
8523.204548822773
+ s
* (
3390.1996842891790
+ s
* (
1127.1366883800605
+ s * (309.69425190845478 + s * (67.751604120864452 + s * (10.834010773474702 + s)))
)
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
s = 1.0 - t
fd = (
(
5.69341524462926732e11
+ t
* (
2.36489781301030086e12
+ t
* (
4.67723648136547218e12
+ t
* (
5.73677938363737470e12
+ t
* (
4.75966348844286492e12
+ t
* (
2.72977815388509051e12
+ t
* (1.05619383963810221e12 + t * (2.52484808268258635e11 + t * 2.87129101561538399e10))
)
)
)
)
)
)
/ (
233403.35536095456
+ s
* (
60565.447915936497
+ s
* (
18770.529358157892
+ s
* (
5494.0347508271098
+ s
* (
1458.2356321178129
+ s * (341.00432855822186 + s * (67.360611688552025 + s * (10.355894270482783 + s)))
)
)
)
)
)
* 0.999999999999999687
)
elif x < 20.0:
t = 0.1 * x - 1.0
s = 1.0 - t
fd = (
(
6.79824678667324436e13
+ t
* (
4.47889271605157194e14
+ t
* (
1.33539302658610566e15
+ t
* (
2.35081455707368313e15
+ t
* (
2.67105097393954020e15
+ t
* (
2.00614187037046236e15
+ t
* (9.73640078908216993e14 + t * (2.79804648854039015e14 + t * 3.66107518187916541e13))
)
)
)
)
)
)
/ (
585167.34458490542
+ s
* (
99101.77544747319
+ s
* (
23659.896891961921
+ s
* (
5847.0608231871316
+ s
* (
1394.1981365390388
+ s * (307.36846956466131 + s * (59.734384735038341 + s * (9.438593558054937 + s)))
)
)
)
)
)
* 0.999999999999999620
)
elif x < 40.0:
t = 0.05 * x - 1.0
s = 1.0 - t
fd = (
(
1.99221037305471970e16
+ t
* (
1.55301826780028291e17
+ t
* (
5.3562790972794281e17
+ t
* (
1.06787028331863552e18
+ t
* (
1.34677794609411928e18
+ t
* (
1.10118289362700600e18
+ t
* (5.7076444557408842e17 + t * (1.71793041083229002e17 + t * 2.30600690809359773e16))
)
)
)
)
)
)
/ (
959281.46439522639
+ s
* (
146009.54484661068
+ s
* (
31788.065669901676
+ s
* (
7257.6052635887384
+ s
* (
1618.0190118885021
+ s * (337.33224463190506 + s * (62.698857808593658 + s * (9.5893450572160768 + s)))
)
)
)
)
)
* 0.999999999999999649
)
else:
w = 1.0 / (x * x)
t = 1600.0 * w
fd = (
x
* x
* x
* x
* x
* x
* x
* x
* sqrt(x)
* factor
* (
1.0
+ w
* (
104.864546761574387
+ t
* (
2.69793375997248149
+ t * (0.0276379588953803004 + t * (0.0000654786006189425666 - t * 1.03209499898071117e-8))
)
)
)
)
return fd
@fd_doc_iparams(order=8)
@jit(float64(float64), nopython=True)
def fd16h(y: float) -> float:
x = -abs(y)
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
40320.0
- ex
* (438381.668209835602 + t * (14167.5566380822417 + t * 39.1239843501671193))
/ (5566.7513423478282 + t * (199.502568378336498 + t))
)
elif x <= 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
2.14460848478648315e8
+ t
* (
1.66513461191585055e8
+ t
* (
5.17491656951190442e7
+ t
* (
6.60393533940475286e6
+ t * (-301119.943009454383 + t * (-202842.867158865294 - t * 20669.3789063593406))
)
)
)
) / (
10903.784789014007
+ s
* (
14623.305826998018
+ s
* (
9179.893497142397
+ s
* (
3522.2940212081378
+ s * (905.1498999681050 + s * (159.07578981158691 + s * (18.015010064955669 + s)))
)
)
)
)
if y > 0.0:
y2 = y * y
fd = fd + y * (
80336.2292693975266
+ y2
* (13245.8066670375278 + y2 * (636.406061422149257 + y2 * (13.1594725347858115 + y2 * 0.111111111111111111)))
)
return fd
@fd_doc_hparams(order=17 / 2)
@jit(float64(float64), nopython=True)
def fd17h(x: float) -> float:
factor = 2.0 / 19.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
s = 1.0 - t
s2 = s * s
fd = ex * (
119292.461994609007
- ex
* (
(164.281538211493590 + s2 * (0.00397356161994318188 + s2 * 1.50191715167058393e-6))
+ s * (0.465418918711191171 + s2 * 0.0000615968461595874662)
)
)
elif x < 0.0:
s = -0.5 * x
s2 = s * s
t = 1.0 - s
t2 = t * t
fd = (
(
1.43233835112073902e8
+ t2 * (6.13766014358067700e7 + t2 * (2.43295219093432312e6 + t2 * 8461.46979813306275))
)
+ t * (1.39644553704144873e8 + t2 * (1.56246791599180488e7 + t2 * 218176.915876999296))
) / (
(3043.1981753546660 + s2 * (1761.0786885568746 + s2 * (112.08528557398394 + s2)))
+ s * (3390.9687009373892 + s2 * (550.68360665579008 + s2 * 14.640373761706403))
)
elif x < 2.0:
t = 0.5 * x
t2 = t * t
s = 1.0 - t
s2 = s * s
fd = (
(
1.21473360069923149e9
+ t2 * (4.95296056344158641e8 + t2 * (7.0192333337028334e7 + t2 * 2.53144559406161612e6))
)
+ t
* (
1.01176690309922349e9
+ t2 * (2.06442737728860908e8 + t2 * (1.68572449951275985e7 + t2 * 186718.775025914870))
)
) / (
(3455.142635040908 + s2 * (2023.9784836523565 + s2 * (308.93682588995655 + s2 * 12.338920623517048)))
+ s * (3415.073282258220 + s2 * (903.6597935624893 + s2 * (76.49008105409234 + s2)))
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
t2 = t * t
s = 1.0 - t
s2 = s * s
fd = (
(
2.67440099736518804e10
+ t2
* (
4.73732500104890339e10
+ t2 * (1.59135818804699047e10 + t2 * (1.74615616554871192e9 + t2 * 3.46053509281020592e7))
)
)
+ t
* (
4.91370104189578919e10
+ t2 * (3.16238224173646457e10 + t2 * (6.12061461496369111e9 + t2 * 3.38446352294444219e8))
)
) / (
(11121.047017722056 + s2 * (5748.503365098324 + s2 * (953.9183383047316 + s2 * (64.892044676666959 + s2))))
+ s * (9829.804681621007 + s2 * (2606.1477102230902 + s2 * (281.64461270648236 + s2 * 10.720667663673541)))
)
elif x < 10.0:
t = 0.2 * x - 1.0
t2 = t * t
s = 1.0 - t
s2 = s * s
fd = (
(
(
1.21689593180741239e12
+ t2
* (
9.12838410654446290e12
+ t2 * (8.99572172197880880e12 + t2 * (2.05751238857523311e12 + t2 * 6.19902241659841813e10))
)
)
+ t
* (
4.79662223812667505e12
+ t2 * (1.09342887556403655e13 + t2 * (5.19592600765981107e12 + t2 * 5.12310527144160372e11))
)
)
/ (
(
39255.392055313674
+ s2 * (9415.210433779872 + s2 * (1027.9861625160893 + s2 * (58.326091875891505 + s2)))
)
+ s
* (22199.471524635690 + s2 * (3354.6940213500329 + s2 * (269.15608801102034 + s2 * 9.690839898241749)))
)
* 0.999999999999999611
)
elif x < 20.0:
t = 0.1 * x - 1.0
t2 = t * t
s = 1.0 - t
s2 = s * s
fd = (
(
(
8.36136698862497927e13
+ t2
* (
1.73344511753855139e15
+ t2 * (3.79488820022826435e15 + t2 * (1.58129984938571763e15 + t2 * 7.24731492178964020e13))
)
)
+ t
* (
5.63471803594835697e14
+ t2 * (3.17713365223451161e15 + t2 * (3.02846836744419910e15 + t2 * 4.96500350365744464e14))
)
)
/ (
(
43633.188864890011
+ s2 * (8040.798889514530 + s2 * (794.3651890403060 + s2 * (46.655407061024955 + s2)))
)
+ s
* (21099.956516665469 + s2 * (2678.4837772160305 + s2 * (208.01106751201808 + s2 * 8.3567416627775934)))
)
* 0.999999999999999524
)
elif x < 40.0:
t = 0.05 * x - 1.0
t2 = t * t
s = 1.0 - t
s2 = s * s
fd = (
(
(
2.85846778936881832e16
+ t2
* (
8.67582358035872224e17
+ t2 * (2.53068942615802160e18 + t2 * (1.29359182644194722e18 + t2 * 6.72189933519645798e16))
)
)
+ t
* (
2.36061789976170730e17
+ t2 * (1.85563671344211925e18 + t2 * (2.25901769990104319e18 + t2 * 4.36848670280351452e17))
)
)
/ (
(
50522.192743560920
+ s2 * (8428.128974587599 + s2 * (782.3736655720094 + s2 * (44.899007032888684 + s2)))
)
+ s
* (23134.306793549278 + s2 * (2708.9200442813342 + s2 * (201.43688280680891 + s2 * 8.0956175513436493)))
)
* 0.999999999999999478
)
else:
w = 1.0 / (x * x)
t = 1600.0 * w
t2 = t * t
fd = (
x
* x
* x
* x
* x
* x
* x
* x
* x
* sqrt(x)
* factor
* (
1.0
+ w
* (
(132.828425897994463 + t2 * (0.0750173171258174240 + t2 * 1.95167395079296505e-7))
+ t * (4.66006740357572083 + t2 * 0.000414697262104018271)
)
)
)
return fd
@fd_doc_iparams(order=9)
@jit(float64(float64), nopython=True)
def fd18h(y: float) -> float:
x = -abs(y)
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
fd = ex * (
362880.0
- ex
* (3.11619728815890900e6 + t * (84459.263086739990 + t * 194.729051789880448))
/ (8793.5020477151129 + t * (258.970670923338654 + t))
)
elif x <= 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
5.34284039564434465e8
+ t
* (
5.33773366525970850e8
+ t
* (
2.42092992302655038e8
+ t
* (
6.42814979745526752e7
+ t * (1.06272445991848042e7 + t * (1.04623859092301352e6 + t * 48187.9003961411735))
)
)
)
) / (
3823.5401164205447
+ s
* (
4167.5033491759870
+ s
* (
2107.9832247353650
+ s * (639.46059745678970 + s * (125.53444672868253 + s * (15.642276104269218 + s)))
)
)
)
if y > 0.0:
y2 = y * y
fd = (
-fd
+ 725062.913034521571
+ y2
* (
361513.031712288870
+ y2 * (29803.0650008344375 + y2 * (954.609092133223885 + y2 * (14.8044066016340379 + y2 * 0.1)))
)
)
return fd
@fd_doc_hparams(order=19 / 2)
@jit(float64(float64), nopython=True)
def fd19h(x: float) -> float:
factor = 2.0 / 21.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
s = 1.0 - t
fd = ex * (
1.13327838894878557e6
- ex
* (
781.077395263087427
+ s
* (
1.48017751152694790
+ s * (0.00952643972760576450 + s * (0.000118887473010527875 + s * 2.40377778104657177e-6))
)
)
)
elif x < 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
1.98943681340803829e9
+ t
* (
2.02504701261695886e9
+ t
* (
9.40112997876683365e8
+ t
* (
2.57205491041485272e8
+ t * (4.42589158347859287e7 + t * (4.61202797242277802e6 + t * 232199.034580285755))
)
)
)
) / (
4645.3648503361484
+ s
* (
4976.8851201391641
+ s
* (
2464.6700288281232
+ s * (729.19053696862632 + s * (138.78343942961328 + s * (16.586263194372411 + s)))
)
)
)
elif x < 2.0:
t = 0.5 * x
s = 1.0 - t
fd = (
2.0904667585793652e10
+ t
* (
1.9423314417935927e10
+ t
* (
9.4466078397530459e9
+ t
* (
3.3380904172376415e9
+ t
* (
9.2799283366292994e8
+ t * (1.9279930343506925e8 + t * (2.6584305328117208e7 + t * 1.8721659680393183e6))
)
)
)
)
) / (
6510.453486834102
+ s
* (
6583.316420184059
+ s
* (
3538.873468101783
+ s
* (
1337.9108450581962
+ s * (388.7430265102922 + s * (85.60461680432659 + s * (12.861326654789701 + s)))
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
s = 1.0 - t
fd = (
1.46863114809614622e11
+ t
* (
2.34916998394298104e11
+ t
* (
1.98120796400802870e11
+ t
* (
1.18887551091630769e11
+ t
* (
5.55928633602927965e10
+ t
* (
2.03490255735994097e10
+ t * (5.60239395267172792e9 + t * (1.05980329750305191e9 + t * 1.07474670751003016e8))
)
)
)
)
)
) / (
4912.718774356158
+ s
* (
5658.595947000664
+ s
* (
3911.882548894156
+ s
* (
2003.2886387850459
+ s
* (
805.2407096838510
+ s * (255.52601531171211 + s * (62.070351460815245 + s * (10.603837914443282 + s)))
)
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
s = 1.0 - t
fd = (
(
4.12413857105017004e12
+ t
* (
1.53777683873017139e13
+ t
* (
2.80178767823520436e13
+ t
* (
3.25874492182871964e13
+ t
* (
2.64367160554342534e13
+ t
* (
1.53023378021225309e13
+ t
* (6.17571294140739606e12 + t * (1.59629111145863694e12 + t * 2.05057950138673222e11))
)
)
)
)
)
)
/ (
9870.066781117703
+ s
* (
8442.670237353289
+ s
* (
4641.969549915204
+ s
* (
1995.0747902033748
+ s
* (
707.3070308400987
+ s * (208.47240484486439 + s * (49.869236685972587 + s * (9.0113347734666514 + s)))
)
)
)
)
)
* 0.999999999999999624
)
elif x < 20.0:
t = 0.1 * x - 1.0
s = 1.0 - t
fd = (
(
4.83390017425457148e15
+ t
* (
3.49375860214157535e16
+ t
* (
1.17174577652359819e17
+ t
* (
2.38956570804562992e17
+ t
* (
3.26219845844040472e17
+ t
* (
3.09099843298003147e17
+ t
* (
2.03426110846076766e17
+ t
* (
8.98530419590583886e16
+ t * (2.42666818156258567e16 + t * 3.07539397786247273e15)
)
)
)
)
)
)
)
)
/ (
207639.00554015383
+ s
* (
101594.67784669785
+ s
* (
39318.472360728631
+ s
* (
13385.302295791795
+ s
* (
4099.2485074048775
+ s
* (
1128.0104124272086
+ s * (274.41287517597127 + s * (56.961169613859215 + s * (9.354997067031336 + s)))
)
)
)
)
)
)
* 0.999999999999999475
)
elif x < 40.0:
t = 0.05 * x - 1.0
s = 1.0 - t
fd = (
(
2.81147322924282096e18
+ t
* (
2.55872925207166842e19
+ t
* (
1.05126940075992923e20
+ t
* (
2.56086191531691937e20
+ t
* (
4.08006931722902169e20
+ t
* (
4.41541112498720876e20
+ t
* (
3.25219940350974164e20
+ t
* (
1.57672402355387744e20
+ t * (4.58547965512713670e19 + t * 6.13678661161381742e18)
)
)
)
)
)
)
)
)
/ (
253648.48647211156
+ s
* (
116226.83914995012
+ s
* (
42576.299330336023
+ s
* (
13857.868987762036
+ s
* (
4097.2448088308021
+ s
* (
1099.1091735270184
+ s * (263.31471897218326 + s * (54.436947056837050 + s * (9.0281927555123252 + s)))
)
)
)
)
)
)
* 0.999999999999999480
)
else:
w = 1.0 / (x * x)
t = 1600.0 * w
fd = (
x
* x
* x
* x
* x
* x
* x
* x
* x
* x
* sqrt(x)
* factor
* (
1.0
+ w
* (
164.082173168110588
+ t
* (
7.52780119040598426
+ t
* (
0.175040406459615661
+ t * (0.00174172934550196671 + t * (4.09368331231164635e-6 - t * 6.44677689509362889e-10))
)
)
)
)
)
return fd
@fd_doc_iparams(order=10)
@jit(float64(float64), nopython=True)
def fd20h(y: float) -> float:
x = -abs(y)
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
s = 1.0 - t
fd = ex * (
3628800.0
- ex
* (
1769.11836499240874
+ s
* (
2.74114748349296929
+ s * (0.0153129509121134183 + s * (0.000171415624076050339 + s * 3.15741345799313491e-6))
)
)
)
elif x <= 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
7.23923596166699396e9
+ t
* (
7.46011115766110324e9
+ t
* (
3.51560076972821513e9
+ t
* (
9.79963332403905796e8
+ t * (1.72722667850513422e8 + t * (1.85872779338803730e7 + t * 980119.027043354796))
)
)
)
) / (
5345.1738990257193
+ s
* (
5661.6190861926977
+ s
* (
2763.5410988949565
+ s * (803.48749805350521 + s * (149.54574471403198 + s * (17.324182954806273 + s)))
)
)
)
if y > 0.0:
y2 = y * y
fd = fd + y * (
7.25062913034521571e6
+ y2
* (
1.20504343904096290e6
+ y2
* (
59606.1300016688751
+ y2 * (1363.72727447603412 + y2 * (16.4493406684822644 + y2 * 0.0909090909090909091))
)
)
)
return fd
@fd_doc_hparams(order=21 / 2)
@jit(float64(float64), nopython=True)
def fd21h(x: float) -> float:
factor = 2.0 / 23.0
if x < -2.0:
ex = exp(x)
t = ex * 7.38905609893065023
s = 1.0 - t
fd = ex * (
1.18994230839622485e7
- ex
* (
4103.25503466860405
+ s
* (
5.19741906649655513
+ s * (0.0251967194122553158 + s * (0.000252954565705636841 + s * 4.24573915094571073e-6))
)
)
)
elif x < 0.0:
s = -0.5 * x
t = 1.0 - s
fd = (
2.57000236613575175e10
+ t
* (
2.66775768017824432e10
+ t
* (
1.26823522594150685e10
+ t
* (
3.57329482988205525e9
+ t * (6.38366360660921950e8 + t * (6.99165387034866263e7 + t * 3.77747237079641640e6))
)
)
)
) / (
5829.6137460125218
+ s
* (
6133.1261260721688
+ s
* (
2967.5333196500928
+ s * (853.66205862182849 + s * (156.68183527916256 + s * (17.793529538010685 + s)))
)
)
)
elif x < 2.0:
t = 0.5 * x
s = 1.0 - t
fd = (
3.9256463321926022e11
+ t
* (
3.5272008321450459e11
+ t
* (
1.4960338986194452e11
+ t
* (
4.0840964995951679e10
+ t
* (
8.3525752792093004e9
+ t * (1.3652900099806473e9 + t * (1.6653931710965140e8 + t * 1.13722223824443501e7))
)
)
)
)
) / (
10780.778916243366
+ s
* (
12211.902429240048
+ s
* (
6810.538241268853
+ s
* (
2444.9651837399508
+ s * (622.4594160475152 + s * (115.21708456233468 + s * (14.653186311439082 + s)))
)
)
)
)
elif x < 5.0:
t = 0.3333333333333333333 * (x - 2.0)
s = 1.0 - t
fd = (
1.02038438555518024e12
+ t
* (
1.44086845004445038e12
+ t
* (
1.06951376053758491e12
+ t
* (
5.83349119667798955e11
+ t
* (
2.59184513055987834e11
+ t
* (
9.22344989320630316e10
+ t * (2.46767885435893916e10 + t * (4.50298910194390990e9 + t * 4.40533249962065063e8))
)
)
)
)
)
) / (
2627.6481435132529
+ s
* (
3598.485403241250
+ s
* (
2816.3184982451496
+ s
* (
1591.2052134075520
+ s
* (
693.3250333689612
+ s * (234.74251456309597 + s * (59.835610870372227 + s * (10.529060212037859 + s)))
)
)
)
)
)
elif x < 10.0:
t = 0.2 * x - 1.0
s = 1.0 - t
fd = (
(
1.92539440493876721e13
+ t
* (
6.76427905453399646e13
+ t
* (
1.17204818714740525e14
+ t
* (
1.31312575057612394e14
+ t
* (
1.04191636638567130e14
+ t
* (
5.99784170385635521e13
+ t
* (2.45015915838649510e13 + t * (6.53409875052899427e12 + t * 8.86443190814042498e11))
)
)
)
)
)
)
/ (
3208.0283085434934
+ s
* (
3638.3347710577402
+ s
* (
2454.8591095140763
+ s
* (
1237.5705806010665
+ s
* (
499.5532233797705
+ s * (164.10632643205828 + s * (43.045792336177945 + s * (8.4174363747925276 + s)))
)
)
)
)
)
* 0.999999999999999649
)
elif x < 20.0:
t = 0.1 * x - 1.0
s = 1.0 - t
fd = (
(
1.35032839928658828e16
+ t
* (
9.84341732357224541e16
+ t
* (
3.35732954478237964e17
+ t
* (
7.02136022561204639e17
+ t
* (
9.91445873700748081e17
+ t
* (
9.80394164297009879e17
+ t
* (
6.79928848238746956e17
+ t
* (
3.19982319837818266e17
+ t * (9.33078432378309346e16 + t * 1.29960135147133318e16)
)
)
)
)
)
)
)
)
/ (
32616.320007579440
+ s
* (
25996.868562405952
+ s
* (
13806.569048984017
+ s
* (
5938.718527590630
+ s
* (
2189.8437434853080
+ s
* (
703.2974320942947
+ s * (195.48717008127781 + s * (45.691859712028360 + s * (8.3721801701244348 + s)))
)
)
)
)
)
)
* 0.999999999999999529
)
elif x < 40.0:
t = 0.05 * x - 1.0
s = 1.0 - t
fd = (
(
1.01458237235662575e19
+ t
* (
9.67883827639700444e19
+ t
* (
4.18744967361252695e20
+ t
* (
1.07961610867931522e21
+ t
* (
1.83106271883557796e21
+ t
* (
2.12354501642149972e21
+ t
* (
1.68956127221959607e21
+ t
* (
8.93514699913180709e20
+ t * (2.87036545075886579e20 + t * 4.31795348604572483e19)
)
)
)
)
)
)
)
)
/ (
32461.378904270757
+ s
* (
24848.502722479329
+ s
* (
12785.219096902698
+ s
* (
5373.990729362032
+ s
* (
1953.3664643688677
+ s
* (
624.0842567143388
+ s * (174.30734991710066 + s * (41.427358351772922 + s * (7.8380857354891532 + s)))
)
)
)
)
)
)
* 0.999999999999999445
)
else:
w = 1.0 / (x * x)
t = 1600.0 * w
fd = (
x**11
* sqrt(x)
* factor
* (
1.0
+ w
* (
198.625788571923339
+ t
* (
11.5426284919561671
+ t
* (
0.365993577138676996
+ t * (0.00572282501335319477 + t * (0.0000313848422097767888 + t * 1.47431097320479878e-8))
)
)
)
)
)
return fd