# 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
