#|
Compute the polynomial part of Stirling's approximation
|#

#|
(fllog-gamma x) = (+ (flstirling x)
                     (* x (fllog x))
                     (- x)
                     (* 0.5 (fllog (/ (* 2.0 pi) x))))

(flgamma x) = (* (flexp-stirling x)
                 (flexpt x x)
                 (flexp (- x))
                 (flsqrt (/ (* 2.0 pi) x)))
|#

(module math.functions.stirling-error (fpstirling
                                       fpexp-stirling)
  (import scheme
          chicken.base
          chicken.type
          chicken.flonum
          math.base
          math.flonum.functions
          math.flonum.polyfun
          math.polynomial.chebyshev
          math.racket-shim)


(: fpstirling-15-27 (float -> float))
(define (fpstirling-15-27 x)
  (fp/ x ((inline-chebyshev-fppoly-fun
           15.0 27.0
           (11016.79948099981773627209876968885362204
            3024.00014800260474631678762744430635965
            215.9999679230867472929981033536003535621
            6.197575655013797795140267332559371048751e-6
            -1.123482005324992955884916740741637685148e-6
            1.955324929224623311339400035192381837496e-7
            -3.307831282170910797335684021141647395959e-8
            5.479661724561398559219886067154634410417e-9
            -8.931610760513212080760137735821652297713e-10
            1.437126105618637970121702240104527823475e-10
            -2.288053816318008707333061098567897054395e-11
            3.610695995893100197394046453166382200035e-12
            -5.652126896275849502824249809697970828599e-13
            8.590674950213391909741768525418882339807e-14))
          x)))

(: fpstirling-8-15 (float -> float))
(define (fpstirling-8-15 x)
  (fp/ x ((inline-chebyshev-fppoly-fun
           8.0 15.0
           (3321.798247797932459144974149656409495232
            966.0005298901628397215169376532205709987
            73.49987816788887093424892777593993146966
            2.494733140666379193517223900773536294113e-5
            -4.787074309492002066556664123607488277298e-6
            8.807687744072022638350005765442426267097e-7
            -1.573092394772476574574240114018318793612e-7
            2.747656271318921269727872382117004258265e-8
            -4.715970675936475871998141501178029715698e-9
            7.980147231866560869562868608196647541232e-10
            -1.334466185130166386044851438572428687415e-10
            2.209136885186544813721858028010037067248e-11
            -3.625194612470337046498524409054506251461e-12
            5.89961828339772115897107672853833190289e-13
            -9.302513801056728846652906602460974539922e-14))
          x)))

(: fpstirling-4.5-8 (float -> float))
(define (fpstirling-4.5-8 x)
  (fp/ x ((inline-chebyshev-fppoly-fun
           4.5 8.0
           (975.0443031907733133389158279812611499587
            262.5015627132377516143395761208780936736
            18.3746764141027136751383456417283990405
            5.938082786318118011724703482868275638379e-5
            -1.016730602988921204230885968680626948872e-5
            1.662515635381875817794816512099954740469e-6
            -2.628928792442607777063166011284861608561e-7
            4.050822127038092335113827822346792254026e-8
            -6.112487132453923578260625403517465865551e-9
            9.06364081362638907112852405510257653843e-10
            -1.323998679460683118169932740761715566757e-10
            1.908950450244905822866094061486038536357e-11
            -2.720586588977906872800608377555774608691e-12
            3.837044294086912170963236288848111628104e-13
            -5.358780479116698179655468776097177855545e-14
            7.28528236204856563656658475438575529055e-15))
          x)))

(: fpstirling-2.75-4.5 (float -> float))
(define (fpstirling-2.75-4.5 x)
  (fp/ x ((inline-chebyshev-fppoly-fun
           2.75 4.5
           (325.3467595768007572161779845085023817787
            76.12859544638192139459089603623245542144
            4.593138374223773581286859767098963154362
            9.131333736750069011356870865537489178965e-5
            -1.261380743745780228407758034482080830216e-5
            1.651593135371361224051430190039696407515e-6
            -2.077106380119501929519954989431521108645e-7
            2.529653473668130040283553911283657922475e-8
            -2.999709265632592433769951001453104388127e-9
            3.476863330588570351919415335411696337828e-10
            -3.950193609609407291091876605408976762817e-11
            4.408617684109024142013268467557430326734e-12
            -4.841203680484443364188480509836354990709e-13
            5.237422339007086632193845068503408696153e-14
            -5.586740617646839083565522161226600650702e-15
            5.819151022309982950144818176644406281669e-16))
          x)))

(: fpstirling-1.75-2.75 (float -> float))
(define (fpstirling-1.75-2.75 x)
  (fp/ x ((inline-chebyshev-fppoly-fun
           1.75 2.75
           (125.2629230772482028509095781215440101143
            27.0072859608297381233699496757487406466
            1.498952409460891689508422629474544515631
            1.303994200168118662091579592936614198868e-4
            -1.484656128624609178802098025636180332568e-5
            1.586277471840025118149864944019821162874e-6
            -1.613450944444201653180199914617430303249e-7
            1.576196350202965718756195765505628839837e-8
            -1.487529043170189643984145466700307624064e-9
            1.361396867785753372873967352158707789021e-10
            -1.211139487543343017696681192325640474208e-11
            1.048488661663920779172817513037101780685e-12
            -8.830149928864984808376749130424857464651e-14
            7.17803116949183817976881871715031957752e-15))
          x)))

(: fpstirling-1-1.75 (float -> float))
(define (fpstirling-1-1.75 x)
  (fp/ x ((inline-chebyshev-fppoly-fun
           1.0 1.75
           (47.777314237104471778486946812101904408
            12.39297816793955917929120829286561089801
            0.8410330212103364337410690481702160644618
            3.491490532470952098652713628408706066953e-4
            -4.041038030805226839115326915532433949952e-5
            4.326596765075402520137649667194718278897e-6
            -4.345154866576609604930256428231288658095e-7
            4.119728067409649315207082217960099828835e-8
            -3.688214630251487303562262057578976086001e-9
            3.093381138554054580872113820898196059744e-10
            -2.375075896907996495765945952278720383373e-11
            1.564882039047195989899222646804037705211e-12
            -6.828937818750811464962353946511915270149e-14
            -2.634943298133669306444329836343107130006e-15
            1.246448725672418721361338401869402106648e-15))
          x)))

(: fpstirling-0.5-1 (float -> float))
(define (fpstirling-0.5-1 x)
  (fp/ (fp/ ((inline-chebyshev-fppoly-fun
              0.5 1.0
              (0.1194805870082255169722538097612297787364
               0.02135330110301585139295937835983603110036
               -3.196025725983946277932392443714902585582e-5
               -7.124517413424091772083413675848741107899e-7
               6.740176044636219727019993384779194346055e-7
               -1.521613091258975783466558010817018177872e-7
               2.670498821139674650221603532664569395499e-8
               -4.213924812826837370071299439461672434822e-9
               6.29854364241318961398679285759803936539e-10
               -9.147744510967378381527882109301861271705e-11
               1.30962582774325190596052614795108510579e-11
               -1.864252676220728470927650730310462587513e-12
               2.652845480456961816308338061503354760458e-13
               -3.786115147736501189751766017357422211683e-14
               5.429913212512962097859788722473828366346e-15
               -7.833776841686429592284675227535973006279e-16
               1.137470065654636255090782592536038800769e-16
               -1.661731102056490017849945642459553975973e-17
               2.391847876948248930795972764225994129252e-18))
             x) x) x))

(: fplog-gamma1p-taylor-0.25 (float -> float))
(define (fplog-gamma1p-taylor-0.25 x)
  (fp* x ((make-fppolyfun
           (-3.930873456872526458554152107865433609025e-1
            6.625352492439489510635402587594391341412e-1
            -2.554826879615743257115227963673084935423e-1
            1.367707934883504139541130662399073826264e-1
            -8.33925057553356520179116267871200690717e-2
            5.439040175688482279740808974984708405172e-2
            -3.692956435369459194525902820018967491214e-2
            2.57438576595354260979378543783296811343e-2
            -1.828207612943360135796036054648999203789e-2
            1.316309381436802721935242956408783247211e-2
            -9.579129355901410837410770822828309502828e-3
            7.030888316926896871662033116131894701662e-3
            -5.19700944806749099503387802706683780591e-3
            3.864280363537579178918868416954592735695e-3
            -2.88790175803055025196965940208247852597e-3
            2.167711417892748993423507278259085698205e-3
            -1.63339016229319926066901111641248008815e-3
            1.234965287649728727056306104178162094552e-3
            -9.365596024752160424737600167029560750915e-4
            7.121915716550684289534550746747099228174e-4
            -5.429051782355057412755621566583833039375e-4
            4.147802845104665920833274078635131825501e-4
            -3.175365018203068467021525584288060261718e-4))
          (fp- x 0.25))))

(: fpstirling (float -> float))
(define (fpstirling x)
  (cond [(fp<= x 0.0)  (if (fp= x 0.0) +inf.0 +nan.0)]
        [(fp< x 1.0)
         (cond [(fp< x 1e-300)
                (fp- (fp* -0.5 (fplog (fp* 2.0 pi)))
                     (fp* 0.5 (fplog x)))]
               [(fp< x 0.5)
                (+ (fplog-gamma1p-taylor-0.25 x)
                   (fp* -0.5 (fplog (* 2.0 pi x)))
                   (- (fp* x (fplog x)))
                   x)]
               [else  (fpstirling-0.5-1 x)])]
        [(fp< x 8.0)
         (cond [(fp< x 2.75)
                (cond [(fp< x 1.75)  (fpstirling-1-1.75 x)]
                      [else  (fpstirling-1.75-2.75 x)])]
               [(fp< x 4.5)   (fpstirling-2.75-4.5 x)]
               [else  (fpstirling-4.5-8 x)])]
        [(fp< x 27.0)
         (cond [(fp< x 15.0)  (fpstirling-8-15 x)]
               [else  (fpstirling-15-27 x)])]
        [(fp< x 4000.0)
         (fp/ ((make-fppolyfun (#i1/12 #i-1/360 #i1/1260 #i-1/1680 #i1/1188))
               (fp/ (fp/ 1.0 x) x))
              x)]
        [(fp< x 2e7)
         (fp/ ((make-fppolyfun (#i1/12 #i-1/360)) (fp/ (fp/ 1.0 x) x)) x)]
        [else
         (fp/ #i1/12 x)]))

(: fpexp-stirling (float -> float))
(define (fpexp-stirling x)
  (cond [(fp<= x 0.0)  (if (fp= x 0.0) +inf.0 +nan.0)]
        [(fp< x 1e-300)
         (fp/ (fp/ 1.0 (fpsqrt (fp* 2.0 pi)))
              (fpsqrt x))]
        [(fp< x 1e-17)
         ;; Using Gamma(x) ~ 1/x near 0
         (fp/ (fp* (fpexp x) (fpexpt x (- x)))
              (fpsqrt (* 2.0 pi x)))]
        [(fp< x 0.03)
         (fp/ (* (fpexp (fplog-gamma1p-taylor-0.25 x))
                 (fpexpt x (- x))
                 (fpexp x))
              (fpsqrt (* 2.0 pi x)))]
        [else
         ;; Error is <= 1 ulp
         (fpexp (fpstirling x))])))