Skip to content

Commit f349610

Browse files
oscardssmithoscardssmith
committed
update besselj1
1 parent c9c81d9 commit f349610

File tree

1 file changed

+84
-58
lines changed

1 file changed

+84
-58
lines changed

‎src/besselj.jl

Lines changed: 84 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@
2323
# 2009 19th IEEE Symposium on Computer Arithmetic. IEEE, 2009.
2424
#
2525

26-
const J0_ROOTS(::Type{Float64}) = (
26+
J0_ROOTS(::Type{Float64}) = (
2727
(2.4048255576957730, -1.176691651530894e-16),
2828
(3.8317059702075125, -1.5269184090088067e-16),
2929
(5.5200781102863110, 8.088597146146722e-17),
@@ -41,58 +41,58 @@ const J0_ROOTS(::Type{Float64}) = (
4141
(24.352471530749302, 9.169067133951066e-16),
4242
(25.903672087618382, 4.894530726419825e-16)
4343
)
44-
const J0_POLYS(::Type{Float64}) = (
45-
(0.0000000000000000000, -0.5191474972894669, 0.10793870175491979, 0.05660177443794914, -0.008657669593292222,
46-
-0.0021942003590739974, 0.0002643770365942964, 4.37291931443113e-5, -4.338825419759815e-6, -5.304927679598406e-7,
44+
@inline J0_POLYS(::Type{Float64}) = (
45+
(0.0000000000000000000, -0.5191474972894669, 0.10793870175491979, 0.05660177443794914, -0.008657669593292222,
46+
-0.0021942003590739974, 0.0002643770365942964, 4.37291931443113e-5, -4.338825419759815e-6, -5.304927679598406e-7,
4747
4.469819175606667e-8, 4.3284827345621115e-9, -3.135111000732148e-10, -2.628876489517534e-11),
48-
(-0.402759395702553, 2.476919088072758e-16, 0.20137969785127532, -0.017518715285670765, -0.013352611033152278,
48+
(-0.402759395702553, 2.476919088072758e-16, 0.20137969785127532, -0.017518715285670765, -0.013352611033152278,
4949
0.0010359438492839443, 0.00037218755624680334, -2.4952042421113972e-5, -5.776086353844158e-6, 3.374317129436161e-7,
5050
5.727482259215452e-8, -2.9561880489355444e-9, -3.905845672635605e-10, 1.971332566705736e-11),
51-
(0.000000000000000000, 0.34026480655836816, -0.030820651425593214, -0.05298855286760721, 0.004631042145890388,
52-
0.002257440229081131, -0.00017518572879518415, -4.6521091062878115e-5, 3.199785909661533e-6, 5.716500268232186e-7,
51+
(0.000000000000000000, 0.34026480655836816, -0.030820651425593214, -0.05298855286760721, 0.004631042145890388,
52+
0.002257440229081131, -0.00017518572879518415, -4.6521091062878115e-5, 3.199785909661533e-6, 5.716500268232186e-7,
5353
-3.5112898510841466e-8, -4.684643389757727e-9, 2.562685034682206e-10, 2.7958958795750104e-11),
54-
(0.30011575252613254, -1.6640272822046001e-16, -0.15005787626306408, 0.007129737603121546, 0.011742619737383848,
55-
-0.0006260583453094324, -0.00035093119008693753, 1.7929701912295164e-5, 5.6239324892485754e-6, -2.668437501970219e-7,
54+
(0.30011575252613254, -1.6640272822046001e-16, -0.15005787626306408, 0.007129737603121546, 0.011742619737383848,
55+
-0.0006260583453094324, -0.00035093119008693753, 1.7929701912295164e-5, 5.6239324892485754e-6, -2.668437501970219e-7,
5656
-5.6648488273749086e-8, 2.48117399780498e-9, 3.8876537586241154e-10, -1.6657136713437192e-11),
57-
(0.000000000000000000, -0.27145229992838193, 0.015684124960953488, 0.044033774963413, -0.0025093022271948434,
58-
-0.0020603351551475315, 0.00011243486771159352, 4.482303558813413e-5, -2.288390108003442e-6, -5.679383588459768e-7,
57+
(0.000000000000000000, -0.27145229992838193, 0.015684124960953488, 0.044033774963413, -0.0025093022271948434,
58+
-0.0020603351551475315, 0.00011243486771159352, 4.482303558813413e-5, -2.288390108003442e-6, -5.679383588459768e-7,
5959
2.693939234375692e-8, 4.737285529934781e-9, -2.0612709555352797e-10, -2.8163166483726606e-11),
60-
(-0.2497048770578432, 1.1807897766765572e-16, 0.12485243852891914, -0.0040907858517059345, -0.010102792347641438,
61-
0.00038536375952213334, 0.00031859711440332953, -1.2373899600646271e-5, -5.3013932979548665e-6, 2.001098153528186e-7,
60+
(-0.2497048770578432, 1.1807897766765572e-16, 0.12485243852891914, -0.0040907858517059345, -0.010102792347641438,
61+
0.00038536375952213334, 0.00031859711440332953, -1.2373899600646271e-5, -5.3013932979548665e-6, 2.001098153528186e-7,
6262
5.4711629662471434e-8, -1.9724572531751518e-9, -3.8121398193699247e-10, 1.3667679743782715e-11),
63-
(0.000000000000000000, 0.23245983136472478, -0.009857064513825458, -0.03818600911162367, 0.0016073972920762946,
64-
0.0018420433388794816, -7.581358465623415e-5, -4.159284549011371e-5, 1.650645590334915e-6, 5.425453494592871e-7,
63+
(0.000000000000000000, 0.23245983136472478, -0.009857064513825458, -0.03818600911162367, 0.0016073972920762946,
64+
0.0018420433388794816, -7.581358465623415e-5, -4.159284549011371e-5, 1.650645590334915e-6, 5.425453494592871e-7,
6565
-2.0556207467977526e-8, -4.620018928884712e-9, 1.642028058414746e-10, 2.7701605444102412e-11),
66-
(0.21835940724787295, -8.89726402965429e-17, -0.10917970362393398, 0.0027314677279632535, 0.008944552393700088,
67-
-0.00026391472261453583, -0.00028847875053074687, 8.858193371737123e-6, 4.9233776180403375e-6, -1.5077786827161215e-7,
66+
(0.21835940724787295, -8.89726402965429e-17, -0.10917970362393398, 0.0027314677279632535, 0.008944552393700088,
67+
-0.00026391472261453583, -0.00028847875053074687, 8.858193371737123e-6, 4.9233776180403375e-6, -1.5077786827161215e-7,
6868
-5.190218733666561e-8, 1.5539413886301204e-9, 3.674809363354973e-10, -1.1113645791216594e-11),
69-
(0.00000000000000000, -0.20654643307799603, 0.006916736034268416, 0.034115572697347704, -0.001137276252948717,
70-
-0.0016680057255530109, 5.4841792064182565e-5, 3.837965853474541e-5, -1.2335804050962046e-6, -5.106259295634553e-7,
69+
(0.00000000000000000, -0.20654643307799603, 0.006916736034268416, 0.034115572697347704, -0.001137276252948717,
70+
-0.0016680057255530109, 5.4841792064182565e-5, 3.837965853474541e-5, -1.2335804050962046e-6, -5.106259295634553e-7,
7171
1.592333632709497e-8, 4.423517565793139e-9, -1.3138837384184105e-10, -2.6809397212536384e-11),
72-
(-0.19646537146865717, 6.979167865106427e-17, 0.09823268573432613, -0.001988037402152532, -0.008095530671166083,
73-
0.00019440675128712672, 0.0002640383898036336, -6.666777683303928e-6, -4.5715696772304925e-6, 1.1666296153560847e-7,
72+
(-0.19646537146865717, 6.979167865106427e-17, 0.09823268573432613, -0.001988037402152532, -0.008095530671166083,
73+
0.00019440675128712672, 0.0002640383898036336, -6.666777683303928e-6, -4.5715696772304925e-6, 1.1666296153560847e-7,
7474
4.8913639696764225e-8, -1.2379867207945651e-9, -3.508930968415813e-10, 9.07632091591013e-12),
75-
(0.00000000000000000, 0.18772880304043943, -0.005194182350684612, -0.031096513233785917, 0.0008577442641273341,
76-
0.0015312251534677639, -4.184307585284775e-5, -3.5603170534217916e-5, 9.58002109234601e-7, 4.795250964600283e-7,
75+
(0.00000000000000000, 0.18772880304043943, -0.005194182350684612, -0.031096513233785917, 0.0008577442641273341,
76+
0.0015312251534677639, -4.184307585284775e-5, -3.5603170534217916e-5, 9.58002109234601e-7, 4.795250964600283e-7,
7777
-1.263434500625308e-8, -4.20550167685701e-9, 1.065791122326672e-10, 2.5727417675684577e-11),
78-
(0.18006337534431555, -5.638484737644332e-17, -0.09003168767215539, 0.0015299132863046024, 0.0074441453680234426,
79-
-0.00015060569680378768, -0.0002443398416353605, 5.227001519013193e-6, 4.267152607972633e-6, -9.295966007495808e-8,
78+
(0.18006337534431555, -5.638484737644332e-17, -0.09003168767215539, 0.0015299132863046024, 0.0074441453680234426,
79+
-0.00015060569680378768, -0.0002443398416353605, 5.227001519013193e-6, 4.267152607972633e-6, -9.295966007495808e-8,
8080
-4.610438011417262e-8, 1.0049092632275165e-9, 3.339442682325105e-10, -7.4991079301099e-12),
81-
(0.00000000000000000, -0.17326589422922986, 0.004084217951979124, 0.028749284970146657, -0.0006761643016121907,
82-
-0.0014215899173758441, 3.320978125391802e-5, 3.3264379323815026e-5, -7.684444100376941e-7, -4.515479818833484e-7,
81+
(0.00000000000000000, -0.17326589422922986, 0.004084217951979124, 0.028749284970146657, -0.0006761643016121907,
82+
-0.0014215899173758441, 3.320978125391802e-5, 3.3264379323815026e-5, -7.684444100376941e-7, -4.515479818833484e-7,
8383
1.0271599456634145e-8, 3.993903280527723e-9, -8.793975528824604e-11, -2.4610374225652004e-11),
84-
(-0.16718460047381803, 4.662597138876655e-17, 0.08359230023690671, -0.0012242529339116864, -0.006925682915280748,
85-
0.00012100729852042988, 0.00022821854128498174, -4.230799709959437e-6, -4.007618360337058e-6, 7.601217108010105e-8,
84+
(-0.16718460047381803, 4.662597138876655e-17, 0.08359230023690671, -0.0012242529339116864, -0.006925682915280748,
85+
0.00012100729852042988, 0.00022821854128498174, -4.230799709959437e-6, -4.007618360337058e-6, 7.601217108010105e-8,
8686
4.358483157250522e-8, -8.317621452517008e-10, -3.178805429572483e-10, 6.284538399593043e-12),
87-
(0.000000000000000000, 0.16170155068925002, -0.0033200234006036146, -0.02685937038656159, 0.0005505380905909195,
88-
0.0013316994659124243, -2.715683205388875e-5, -3.1289544794362e-5, 6.326900360777323e-7, 4.2697141781883964e-7,
87+
(0.000000000000000000, 0.16170155068925002, -0.0033200234006036146, -0.02685937038656159, 0.0005505380905909195,
88+
0.0013316994659124243, -2.715683205388875e-5, -3.1289544794362e-5, 6.326900360777323e-7, 4.2697141781883964e-7,
8989
-8.532447325590315e-9, -3.799009445641295e-9, 7.379552811590934e-11, 2.353654960834717e-11),
90-
(0.15672498625285222, 1.1464880342445208e-19, -0.07836249312642612, 0.0010083833270351796, 0.006501011610557426,
91-
-9.993664895655577e-5, -0.00021478298462967253, 3.511086890959515e-6, 3.7857022791122945e-6, -6.350944888510364e-8,
90+
(0.15672498625285222, 1.1464880342445208e-19, -0.07836249312642612, 0.0010083833270351796, 0.006501011610557426,
91+
-9.993664895655577e-5, -0.00021478298462967253, 3.511086890959515e-6, 3.7857022791122945e-6, -6.350944888510364e-8,
9292
-4.135588012575766e-8, 7.135988717747828e-10, 3.2075281131621564e-10, 2.208506585533542e-12)
9393
)
94-
const J0_POLY_PIO2(::Type{Float64}) = (
95-
1.000000000000000, -0.25, 0.01562499999999994, -0.00043402777777725544, 6.781684026082576e-6,
94+
J0_POLY_PIO2(::Type{Float64}) = (
95+
1.000000000000000, -0.25, 0.01562499999999994, -0.00043402777777725544, 6.781684026082576e-6,
9696
-6.781683757550061e-8, 4.709479394601058e-10, -2.4016837144506874e-12, 9.104258208703104e-15
9797
)
9898
"""
@@ -106,7 +106,7 @@ function besselj0(x::Float64)
106106

107107
if x < 26.0
108108
x < pi/2 && return evalpoly(x*x, J0_POLY_PIO2(T))
109-
n = round(Int, fma(TWOOPI(T), x, -1/2))
109+
n = unsafe_trunc(Int, TWOOPI(T)*x)
110110
root = @inbounds J0_ROOTS(T)[n]
111111
r = x - root[1] - root[2]
112112
return evalpoly(r, @inbounds J0_POLYS(T)[n])
@@ -126,7 +126,7 @@ function besselj0(x::Float64)
126126
p = evalpoly(x2, p2)
127127
q = evalpoly(x2, q2)
128128
end
129-
129+
130130
a = SQ2OPI(T) * sqrt(xinv) * p
131131
xn = muladd(xinv, q, - PIO4(T))
132132

@@ -160,6 +160,37 @@ function besselj0(x::Float32)
160160
end
161161
end
162162

163+
J1_ROOTS(::Type{Float64}) = (
164+
#(0.0,0.0),
165+
(1.8411837813406593, 4.7898393919093694e-18),
166+
(3.8317059702075125, -1.5269184090088067e-16),
167+
(5.3314427735250325, 1.5109105349471405e-16),
168+
(7.015586669815619, -9.414165653410389e-17),
169+
(8.536316366346286, -1.5433871213307537e-16),
170+
(10.173468135062722, 4.482162274768888e-16),
171+
(11.706004902592063, 7.121366942298246e-16),
172+
(13.323691936314223, 2.600408064718813e-16),
173+
(14.863588633909034, -6.265788988781879e-16),
174+
(16.470630050877634, -1.619019544798128e-15),
175+
(18.015527862681804, -1.1196999448424267e-16),
176+
(19.615858510468243, -1.004445634526616e-15),
177+
(21.16436985918879, 1.7024131380423588e-15),
178+
(22.760084380592772, -4.925749373614922e-16),
179+
(24.311326857210776, -2.614798558537172e-16),
180+
(25.903672087618382, 4.894530726419825e-16),
181+
#(26.1,0.0)
182+
)
183+
184+
@inline J1_POLYS(::Type{Float64}) = ((0.5818652242815964, 1.856296717590682e-18, -0.2051107121477717, 0.006058948324603371, 0.01380176980791771, -0.0003723170971546836, -0.0003949590731742279, 9.202949788197776e-6, 6.267287951499227e-6, -1.2678569694609744e-7, -6.325128726427676e-8, 1.1250631934718714e-9, 4.370887706315301e-10, -6.9196663897189545e-12), (0.0, -0.402759395702553, 0.05255614585697702, 0.053410444132724666, -0.0051797192456287736, -0.002233125339144248, 0.00017466429059297917, 4.6208701269187226e-5, -3.0368626730565984e-6, -5.727804542169931e-7, 3.248095002810744e-8, 4.735119508144983e-9, -2.337340608580871e-10, -2.7702936049865113e-11), (-0.3461262018537915, -3.8155536376435295e-18, 0.16697453550109173, -0.009678268542877979, -0.012099225779107702, 0.0006654009006392731, 0.0003541389004643254, -1.7427203115743037e-5, -5.6552920575085545e-6, 2.4842939771685583e-7, 5.70953426220719e-8, -2.2535737298086475e-9, -3.9378293189634726e-10, 1.4069960915236214e-11), (0.0, 0.30011575252613254, -0.021389212809341387, -0.04697047894974135, 0.0031302917260392776, 0.00210558714324504, -0.00012550790943706365, -4.499147527460248e-5, 2.4015801599445833e-6, 5.665268816114443e-7, -2.72717063313724e-8, -4.7202213932346135e-9, 2.042981112498402e-10, 2.771928613815164e-11), (0.2732999416331998, 3.561226959512617e-18, -0.13477468037992238, 0.005116340346487853, 0.010631861751951011, -0.0004487436837325132, -0.0003268000181954574, 1.338255595522175e-5, 5.363175640292698e-6, -2.064719381656635e-7, -5.499626200248995e-8, 1.9736740545225123e-9, 3.8273222037285807e-10, -1.278031314231906e-11), (0.0, -0.2497048770578432, 0.012272357555101357, 0.040411169390792624, -0.0019268187972532505, -0.0019115826893794257, 8.661729444387141e-5, 4.241116278258051e-5, -1.8009789135751702e-6, -5.471601795172287e-7, 2.168246459437719e-8, 4.63065563452336e-9, -1.693619439040066e-10, -2.7469138621341535e-11), (-0.23330441717143408, -3.1115030009112076e-18, 0.1158009224460766, -0.0032489977328225466, -0.009372527206018624, 0.0003036138211658651, 0.0002980455550043141, -9.813818566107743e-6, -5.02422851525659e-6, 1.6136260037331645e-7, 5.25161116324859e-8, -1.6179928529587045e-9, -3.703476148324225e-10, 1.0858871789346562e-11), (0.0, 0.21835940724787295, -0.008194403183877378, -0.035778209575030466, 0.0013195736128040353, 0.0017308725061717788, -6.200735153199724e-5, -3.938703739294411e-5, 1.3569937679279234e-6, 5.190655478726095e-7, -1.708250372429145e-8, -4.465083874723755e-9, 1.3830562434346436e-10, 2.6807910960740527e-11), (0.20701265272531905, 2.662491243039858e-18, -0.10303781563402645, 0.0022897295462871396, 0.008432435058832219, -0.0002196594196231479, -0.00027272154367319313, 7.369094613922164e-6, 4.682126912827974e-6, -1.2611225276326565e-7, -4.97328501372088e-8, 1.3106601105836534e-9, 3.5526726618793736e-10, -9.062586317137296e-12), (0.0, -0.19646537146865717, 0.005964112206447887, 0.032382122684890685, -0.0009720337562251067, -0.0015842303417534763, 4.666744216252213e-5, 3.657257351175072e-5, -1.0499611451317645e-6, -4.891789717262985e-7, 1.3609384856400248e-8, 4.264401444852941e-9, -1.1318488442382296e-10, -2.5888166054348253e-11), (-0.18801748852581776, -2.2677211286754067e-18, 0.09371909377243096, -0.0017233243363059443, -0.007714266760096894, 0.0001675586237678538, 0.0002518271743879865, -5.733517249250505e-6, -4.372424694665777e-6, 1.0045691177526798e-7, 4.6981261166072386e-8, -1.069605169831798e-9, -3.391413218535611e-10, 7.563513605572664e-12), (0.0, 0.18006337534431555, -0.0045897398589059635, -0.029776581472313393, 0.0007530284838489293, 0.0014660390526531409, -3.658900932512622e-5, -3.413723643387153e-5, 8.366324980452383e-7, 4.6108490576275497e-7, -1.1047178123633162e-8, -4.05903578695289e-9, 9.361926857267905e-11, 2.4863648392841034e-11), (0.1734590492857464, 1.9379085376817515e-18, -0.0865359019387601, 0.001356818985614292, 0.007147216520707795, -0.00013295776592214128, -0.0002346295928358077, 4.60314337029431e-6, 4.1031960356554604e-6, -8.183609100757509e-8, -4.443990880675131e-8, 8.854716308195176e-10, 3.233307974174374e-10, -6.363098067312356e-12), (0.0, -0.16718460047381806, 0.003672758801728573, 0.027702731661334842, -0.0006050364924615997, -0.0013693112504481982, 2.9615596888762443e-5, 3.20609618748376e-5, -6.841058698413295e-7, -4.3588784488057464e-7, 9.143749949363787e-9, 3.864215216864548e-9, -7.850681795998092e-11, -2.3833917223167445e-11), (-0.1618382095526585, -1.6679071380870433e-18, 0.08078219522574892, -0.0011038527651115032, -0.006686444726785339, 0.00010870535077447649, 0.00022030197342732733, -3.791805466469061e-6, -3.871147758354664e-6, 6.805972840337265e-8, 4.21598280598984e-8, -7.445362901530465e-10, -3.0852724534340563e-10, 5.412341302866906e-12), (0.0, 0.15672498625285222, -0.0030251499811055943, -0.026004046442226, 0.0004996832448036509, 0.0012886979076613464, -2.4577609435529096e-5, -3.028561997466637e-5, 5.715983737007265e-7, 4.1362763577190565e-7, -7.703837603662007e-9, -3.685491114666358e-9, 6.675799000603599e-11, 2.2850358357646402e-11), (0.030457219768533664, 0.15255483607689563, -0.018128761087154197, -0.025120198560798643, 0.0017425631959934952, 0.0012355509577850023, -6.532206878098611e-5, -2.8822686464838518e-5, 1.2845508487047282e-6, 3.908248079330429e-7, -1.5441649818886705e-8, -3.4581646719962977e-9, 1.2316553635263845e-10, 2.1284753070823124e-11))
185+
186+
J1_POLY_PIO2(::FLoat64) = (0.5, -0.0624999999999989, 0.002604166666657291, -5.42534721917933e-5, 6.781683542660179e-7, -5.651361336587487e-9, 3.36191211106159e-11, -1.4511302591871352e-13)
187+
for i in 2:2:18
188+
poly = Float64.(Tuple(ratfn_minimax(x->r2(x,pts[i]), (-pi/4, pi/4), 13, 0)[1]))
189+
push!(polys, poly)
190+
poly = Float64.(Tuple(ratfn_minimax(x->r(x,pts[i+1]), (-pi/4, pi/4), 12, 0)[1]))
191+
push!(polys, (0.0, poly...))
192+
end
193+
163194
"""
164195
besselj1(x::T) where T <: Union{Float32, Float64}
165196
@@ -171,33 +202,28 @@ function besselj1(x::Float64)
171202
isinf(x) && return zero(x)
172203

173204
if x <= 5.0
174-
z = x * x
175-
w = evalpoly(z, RP_j1(T)) / evalpoly(z, RQ_j1(T))
176-
w = w * x * (z - 1.46819706421238932572e1) * (z - 4.92184563216946036703e1)
177-
return w
178-
elseif x < 25.0
179-
w = 5.0 / x
180-
z = w * w
181-
p = evalpoly(z, PP_j1(T)) / evalpoly(z, PQ_j1(T))
182-
q = evalpoly(z, QP_j1(T)) / evalpoly(z, QQ_j1(T))
183-
xn = x - THPIO4(T)
184-
sc = sincos(xn)
185-
p = p * sc[2] - w * q * sc[1]
186-
return p * SQ2OPI(T) / sqrt(x)
205+
x <= pi/2 && return x*evalpoly(x*x, J1_POLY_PIO2(T))
206+
n = unsafe_trunc(Int, TWOOPI(T)*x)
207+
root = @inbounds J1_ROOTS(T)[n]
208+
r = x - root[1] - root[2]
209+
return evalpoly(r, @inbounds J1_POLYS(T)[n])
187210
else
211+
xinv = inv(x)
212+
x2 = xinv*xinv
188213
if x < 120.0
189-
p = (one(T), 3/16, -99/512, 6597/8192, -4057965/524288, 1113686901/8388608, -951148335159/268435456, 581513783771781/4294967296)
190-
q = (3/8, -21/128, 1899/5120, -543483/229376, 8027901/262144, -30413055339/46137344, 9228545313147/436207616)
214+
p1 = (one(T), 3/16, -99/512, 6597/8192, -4057965/524288, 1113686901/8388608, -951148335159/268435456, 581513783771781/4294967296)
215+
q1 = (3/8, -21/128, 1899/5120, -543483/229376, 8027901/262144, -30413055339/46137344, 9228545313147/436207616)
216+
p = evalpoly(x2, p1)
217+
q = evalpoly(x2, q1)
191218
else
192-
p = (one(T), 3/16, -99/512, 6597/8192)
193-
q = (3/8, -21/128, 1899/5120, -543483/229376)
219+
p2 = (one(T), 3/16, -99/512, 6597/8192)
220+
q2 = (3/8, -21/128, 1899/5120, -543483/229376)
221+
p = evalpoly(x2, p2)
222+
q = evalpoly(x2, q2)
194223
end
195-
xinv = inv(x)
196-
x2 = xinv*xinv
197224

198-
p = evalpoly(x2, p)
199225
a = SQ2OPI(T) * sqrt(xinv) * p
200-
xn = muladd(xinv, evalpoly(x2, q), - 3 * PIO4(T))
226+
xn = muladd(xinv, q, - 3 * PIO4(T))
201227

202228
# the following computes b = cos(x + xn) more accurately
203229
# see src/misc.jl

0 commit comments

Comments
 (0)