lib/distribution/math_extension/chebyshev_series.rb
# Added by John O. Woods, SciRuby project.
# Derived from GSL-1.9 source files in the specfunc/ and cheb/ dirs.
module Distribution
module MathExtension
# Adapted from GSL-1.9.
class ChebyshevSeries
DATA = { lopx: [2.16647910664395270521272590407,
-0.28565398551049742084877469679,
0.01517767255690553732382488171,
-0.00200215904941415466274422081,
0.00019211375164056698287947962,
-0.00002553258886105542567601400,
2.9004512660400621301999384544e-06,
-3.8873813517057343800270917900e-07,
4.7743678729400456026672697926e-08,
-6.4501969776090319441714445454e-09,
8.2751976628812389601561347296e-10,
-1.1260499376492049411710290413e-10,
1.4844576692270934446023686322e-11,
-2.0328515972462118942821556033e-12,
2.7291231220549214896095654769e-13,
-3.7581977830387938294437434651e-14,
5.1107345870861673561462339876e-15,
-7.0722150011433276578323272272e-16,
9.7089758328248469219003866867e-17,
-1.3492637457521938883731579510e-17,
1.8657327910677296608121390705e-18],
lopxmx: [-1.12100231323744103373737274541,
0.19553462773379386241549597019,
-0.01467470453808083971825344956,
0.00166678250474365477643629067,
-0.00018543356147700369785746902,
0.00002280154021771635036301071,
-2.8031253116633521699214134172e-06,
3.5936568872522162983669541401e-07,
-4.6241857041062060284381167925e-08,
6.0822637459403991012451054971e-09,
-8.0339824424815790302621320732e-10,
1.0751718277499375044851551587e-10,
-1.4445310914224613448759230882e-11,
1.9573912180610336168921438426e-12,
-2.6614436796793061741564104510e-13,
3.6402634315269586532158344584e-14,
-4.9937495922755006545809120531e-15,
6.8802890218846809524646902703e-16,
-9.5034129794804273611403251480e-17,
1.3170135013050997157326965813e-17],
gstar_a: [2.16786447866463034423060819465,
-0.05533249018745584258035832802,
0.01800392431460719960888319748,
-0.00580919269468937714480019814,
0.00186523689488400339978881560,
-0.00059746524113955531852595159,
0.00019125169907783353925426722,
-0.00006124996546944685735909697,
0.00001963889633130842586440945,
-6.3067741254637180272515795142e-06,
2.0288698405861392526872789863e-06,
-6.5384896660838465981983750582e-07,
2.1108698058908865476480734911e-07,
-6.8260714912274941677892994580e-08,
2.2108560875880560555583978510e-08,
-7.1710331930255456643627187187e-09,
2.3290892983985406754602564745e-09,
-7.5740371598505586754890405359e-10,
2.4658267222594334398525312084e-10,
-8.0362243171659883803428749516e-11,
2.6215616826341594653521346229e-11,
-8.5596155025948750540420068109e-12,
2.7970831499487963614315315444e-12,
-9.1471771211886202805502562414e-13,
2.9934720198063397094916415927e-13,
-9.8026575909753445931073620469e-14,
3.2116773667767153777571410671e-14,
-1.0518035333878147029650507254e-14,
3.4144405720185253938994854173e-15,
-1.0115153943081187052322643819e-15],
gstar_b: [0.0057502277273114339831606096782,
0.0004496689534965685038254147807,
-0.0001672763153188717308905047405,
0.0000615137014913154794776670946,
-0.0000223726551711525016380862195,
8.0507405356647954540694800545e-06,
-2.8671077107583395569766746448e-06,
1.0106727053742747568362254106e-06,
-3.5265558477595061262310873482e-07,
1.2179216046419401193247254591e-07,
-4.1619640180795366971160162267e-08,
1.4066283500795206892487241294e-08,
-4.6982570380537099016106141654e-09,
1.5491248664620612686423108936e-09,
-5.0340936319394885789686867772e-10,
1.6084448673736032249959475006e-10,
-5.0349733196835456497619787559e-11,
1.5357154939762136997591808461e-11,
-4.5233809655775649997667176224e-12,
1.2664429179254447281068538964e-12,
-3.2648287937449326771785041692e-13,
7.1528272726086133795579071407e-14,
-9.4831735252566034505739531258e-15,
-2.3124001991413207293120906691e-15,
2.8406613277170391482590129474e-15,
-1.7245370321618816421281770927e-15,
8.6507923128671112154695006592e-16,
-3.9506563665427555895391869919e-16,
1.6779342132074761078792361165e-16,
-6.0483153034414765129837716260e-17],
e11: [-16.11346165557149402600,
7.79407277874268027690,
-1.95540581886314195070,
0.37337293866277945612,
-0.05692503191092901938,
0.00721107776966009185,
-0.00078104901449841593,
0.00007388093356262168,
-0.00000620286187580820,
0.00000046816002303176,
-0.00000003209288853329,
0.00000000201519974874,
-0.00000000011673686816,
0.00000000000627627066,
-0.00000000000031481541,
0.00000000000001479904,
-0.00000000000000065457,
0.00000000000000002733,
-0.00000000000000000108],
e12: [-0.03739021479220279500,
0.04272398606220957700,
-0.13031820798497005440,
0.01441912402469889073,
-0.00134617078051068022,
0.00010731029253063780,
-0.00000742999951611943,
0.00000045377325690753,
-0.00000002476417211390,
0.00000000122076581374,
-0.00000000005485141480,
0.00000000000226362142,
-0.00000000000008635897,
0.00000000000000306291,
-0.00000000000000010148,
0.00000000000000000315],
ae11: [0.121503239716065790,
-0.065088778513550150,
0.004897651357459670,
-0.000649237843027216,
0.000093840434587471,
0.000000420236380882,
-0.000008113374735904,
0.000002804247688663,
0.000000056487164441,
-0.000000344809174450,
0.000000058209273578,
0.000000038711426349,
-0.000000012453235014,
-0.000000005118504888,
0.000000002148771527,
0.000000000868459898,
-0.000000000343650105,
-0.000000000179796603,
0.000000000047442060,
0.000000000040423282,
-0.000000000003543928,
-0.000000000008853444,
-0.000000000000960151,
0.000000000001692921,
0.000000000000607990,
-0.000000000000224338,
-0.000000000000200327,
-0.000000000000006246,
0.000000000000045571,
0.000000000000016383,
-0.000000000000005561,
-0.000000000000006074,
-0.000000000000000862,
0.000000000000001223,
0.000000000000000716,
-0.000000000000000024,
-0.000000000000000201,
-0.000000000000000082,
0.000000000000000017],
ae12: [0.582417495134726740,
-0.158348850905782750,
-0.006764275590323141,
0.005125843950185725,
0.000435232492169391,
-0.000143613366305483,
-0.000041801320556301,
-0.000002713395758640,
0.000001151381913647,
0.000000420650022012,
0.000000066581901391,
0.000000000662143777,
-0.000000002844104870,
-0.000000000940724197,
-0.000000000177476602,
-0.000000000015830222,
0.000000000002905732,
0.000000000001769356,
0.000000000000492735,
0.000000000000093709,
0.000000000000010707,
-0.000000000000000537,
-0.000000000000000716,
-0.000000000000000244,
-0.000000000000000058],
ae13: [-0.605773246640603460,
-0.112535243483660900,
0.013432266247902779,
-0.001926845187381145,
0.000309118337720603,
-0.000053564132129618,
0.000009827812880247,
-0.000001885368984916,
0.000000374943193568,
-0.000000076823455870,
0.000000016143270567,
-0.000000003466802211,
0.000000000758754209,
-0.000000000168864333,
0.000000000038145706,
-0.000000000008733026,
0.000000000002023672,
-0.000000000000474132,
0.000000000000112211,
-0.000000000000026804,
0.000000000000006457,
-0.000000000000001568,
0.000000000000000383,
-0.000000000000000094,
0.000000000000000023],
ae14: [-0.18929180007530170,
-0.08648117855259871,
0.00722410154374659,
-0.00080975594575573,
0.00010999134432661,
-0.00001717332998937,
0.00000298562751447,
-0.00000056596491457,
0.00000011526808397,
-0.00000002495030440,
0.00000000569232420,
-0.00000000135995766,
0.00000000033846628,
-0.00000000008737853,
0.00000000002331588,
-0.00000000000641148,
0.00000000000181224,
-0.00000000000052538,
0.00000000000015592,
-0.00000000000004729,
0.00000000000001463,
-0.00000000000000461,
0.00000000000000148,
-0.00000000000000048,
0.00000000000000016,
-0.00000000000000005],
erfc_xlt1: [1.06073416421769980345174155056,
-0.42582445804381043569204735291,
0.04955262679620434040357683080,
0.00449293488768382749558001242,
-0.00129194104658496953494224761,
-0.00001836389292149396270416979,
0.00002211114704099526291538556,
-5.23337485234257134673693179020e-7,
-2.78184788833537885382530989578e-7,
1.41158092748813114560316684249e-8,
2.72571296330561699984539141865e-9,
-2.06343904872070629406401492476e-10,
-2.14273991996785367924201401812e-11,
2.22990255539358204580285098119e-12,
1.36250074650698280575807934155e-13,
-1.95144010922293091898995913038e-14,
-6.85627169231704599442806370690e-16,
1.44506492869699938239521607493e-16,
2.45935306460536488037576200030e-18,
-9.29599561220523396007359328540e-19],
erfc_x15: [0.44045832024338111077637466616,
-0.143958836762168335790826895326,
0.044786499817939267247056666937,
-0.013343124200271211203618353102,
0.003824682739750469767692372556,
-0.001058699227195126547306482530,
0.000283859419210073742736310108,
-0.000073906170662206760483959432,
0.000018725312521489179015872934,
-4.62530981164919445131297264430e-6,
1.11558657244432857487884006422e-6,
-2.63098662650834130067808832725e-7,
6.07462122724551777372119408710e-8,
-1.37460865539865444777251011793e-8,
3.05157051905475145520096717210e-9,
-6.65174789720310713757307724790e-10,
1.42483346273207784489792999706e-10,
-3.00141127395323902092018744545e-11,
6.22171792645348091472914001250e-12,
-1.26994639225668496876152836555e-12,
2.55385883033257575402681845385e-13,
-5.06258237507038698392265499770e-14,
9.89705409478327321641264227110e-15,
-1.90685978789192181051961024995e-15,
3.50826648032737849245113757340e-16],
erfc_x510: [1.11684990123545698684297865808,
0.003736240359381998520654927536,
-0.000916623948045470238763619870,
0.000199094325044940833965078819,
-0.000040276384918650072591781859,
7.76515264697061049477127605790e-6,
-1.44464794206689070402099225301e-6,
2.61311930343463958393485241947e-7,
-4.61833026634844152345304095560e-8,
8.00253111512943601598732144340e-9,
-1.36291114862793031395712122089e-9,
2.28570483090160869607683087722e-10,
-3.78022521563251805044056974560e-11,
6.17253683874528285729910462130e-12,
-9.96019290955316888445830597430e-13,
1.58953143706980770269506726000e-13,
-2.51045971047162509999527428316e-14,
3.92607828989125810013581287560e-15,
-6.07970619384160374392535453420e-16,
9.12600607264794717315507477670e-17],
sin: [-0.3295190160663511504173,
0.0025374284671667991990,
0.0006261928782647355874,
-4.6495547521854042157541e-06,
-5.6917531549379706526677e-07,
3.7283335140973803627866e-09,
3.0267376484747473727186e-10,
-1.7400875016436622322022e-12,
-1.0554678305790849834462e-13,
5.3701981409132410797062e-16,
2.5984137983099020336115e-17,
-1.1821555255364833468288e-19]
}
PARAMS = {
lopx: [20, -1, 1, 10],
lopxmx: [19, -1, 1, 9],
gstar_a: [29, -1, 1, 17],
gstar_b: [29, -1, 1, 18],
e11: [18, -1, 1, 13],
e12: [15, -1, 1, 10],
ae11: [38, -1, 1, 20],
ae12: [24, -1, 1, 15],
ae13: [24, -1, 1, 15],
ae14: [25, -1, 1, 13],
erfc_xlt1: [19, -1, 1, 12],
erfc_x15: [24, -1, 1, 16],
erfc_x510: [19, -1, 1, 12],
sin: [11, -1, 1, 11]
}
def initialize(coefficients, expansion_order, lower_interval_point, upper_interval_point, single_precision_order)
@c = coefficients.is_a?(Symbol) ? DATA[coefficients] : coefficients
@order = expansion_order
@lower_interval_point = lower_interval_point
@upper_interval_point = upper_interval_point
@single_precision_order = single_precision_order
end
# double * c; /* coefficients */
# int order; /* order of expansion */
# double a; /* lower interval point */
# double b; /* upper interval point */
# int order_sp; /* effective single precision order */
attr_reader :lower_interval_point, :upper_interval_point, :single_precision_order, :order, :c
def a
@lower_interval_point
end
def b
@upper_interval_point
end
def order_sp
@single_precision_order
end
def coefficients(idx)
@c[idx]
end
class << self
# cheb_eval_e in specfunc/cheb_eval.c (gsl-1.9)
def evaluate(series, x, with_error = false)
cs = Math.const_get "#{series.to_s.upcase}_CS"
fail(ArgumentError, "Unrecognized series #{series}") if cs.nil?
d = 0.0
dd = 0.0
y = (2 * x - cs.a - cs.b).quo(cs.b - cs.a)
y2 = 2 * y
e = 0.0
cs.order.downto(1) do |j|
temp = d
d = y2 * d - dd + cs.c[j]
e += (y2 * temp).abs + dd.abs + cs.c[j].abs
dd = temp
end
begin
temp = d
d = y * d - dd + 0.5 * cs.c[0]
e += (y * temp).abs + dd.abs + 0.5 * cs.c[0].abs
end
with_error ? [d, Float::EPSILON * e + (cs.c[cs.order])] : d
end
end
end
ChebyshevSeries::DATA.keys.each do |series|
Math.const_set "#{series.to_s.upcase}_CS", ChebyshevSeries.new(series, *(ChebyshevSeries::PARAMS[series]))
end
end
end