{-# LANGUAGE ForeignFunctionInterface #-}
{-# OPTIONS_GHC -fno-warn-type-defaults #-}

module Math.Reference
    ( err
    , tgamma
    , lgamma
    ) where

import Data.Number.BigFloat (BigFloat, Prec50)
import Math.Gamma.Lanczos   (gammaLanczos, reflect)

-- type Ref = BigFloat (PrecPlus20 (PrecPlus20 (PrecPlus20 Prec50)))
type Ref = BigFloat Prec50

refGamma' :: Ref -> Ref
refGamma' x
    | x > 0 = gammaLanczos g cs x
    | otherwise = case properFraction x of
        (n,0) -> if n > 0 then fromInteger (product [1..n-1]) else 0/0
        _     -> reflect (gammaLanczos g cs) x

g::Fractional a => a
cs::Fractional a => [a]

-- g = 4*pi
-- cs = [9.99999999999999999999999996859805857094649110255858638007521499626667910350917912841651754219621023562696911332e-1,2.41021935915031575671085869856051212251086781069965411063137445014791649207976600743220302513537180075741255296e5,-9.39683442214927629464202549615391497731809604507151681446997378587954930487740163463037083316646661765306020767e5,1.50949462278130416707216346338262717752378737207648778592042488727849492084941928909939218772224986523524827271e6,-1.29196480388800562971760352435290649542587719650708873592882554634240911710117984734995323885290793478347646330e6,6.36694128920477129400076563202799705143678072113371399208674087118282528558129436201334367030862320422658246991e5,-1.82471418180570786650926611167028035492537813853558673501139223570434908719610942315092027198546110180213204434e5,2.93229025802547026821102620826974735041884874040373833440904730504747656548457832570174326117298554369079127589e4,-2.42104856309484391962703141393650395582185463486998318684467564163688570420245166352198654798010724012815552809e3,8.70263976705814825513391156852897437425645148258979915072545297700880110517285556763650793819190847901635114815e1,-9.90350651051541084034863017450573967496724797720902585856207732922619895739215627256300090899964209233833558059e-1,1.77111462470160503721207454939884167682509869539871034455375282246563600388989404816104050917245493955854296122e-3,-6.07979338168558073081122926060936198873479110416640952420808116982632249326161497364670524879511535821535057852e-8,2.41177003165992726532912316132624759434987463273639173113847102284810676871579183303359325047268888367893750488e-11,-3.29368781819147402947052101243060083111222984364433304208981956498545487091815683441661638260978104468249346673e-11,2.24453115729169763203131076060244943071934132012991691257619864969664661641236008681074860160604017290650304689e-11,-9.42100302390236927123870056289010276539140677886848718235340783874188089528640301130431598050200664465171227938e-12,2.61294791407497062121812006700436259418110616015207527385647574707939568053359405099396866356159720169068766212e-12,-4.32715100268426644428894578150555672315616608786268711314025106092325719775657429204891211631229932892886729036e-13,2.66596243373922858727920730746965698819332806551704556636516706915515953154082281969632365538824477332141826858e-14,1.77321356734488300746927561707710031581163520443393326611478105950839047661826135657152610805311480453596020540e-15]

-- g = 5*pi
-- cs = [9.99999999999999999999999999993886292491620692316160598686599867606051722036264024489964769505713926399764566018e-1,6.26152503787335403652715374769447739961891146555890980978381175134572365593726917595153339118464456879028200299e6,-3.16335440512758530015177280043484225515804122481774619669315038580009080307294118319420059535265919967069772033e7,6.88835095268754493780533335352036097825929449727064875766761598755949365626134351595870710713568770641108217825e7,-8.46956580333456325307500093500724512731750456929320232801425203134598715275876674318720648707327950233819645139e7,6.47329754881842199772339353831484565977006032033099687741144470787747326331686889977791028752629813425598481151e7,-3.19268433191141469791874286821381470267094549847221319695694408120309390209523512867822254931389981528998092501e7,1.02237086392200672155858881540024824965037923176630968557758516768451829832835513494013895382528680690315693015e7,-2.08890272303961042498315400960754400018306410875319695466462560199939031674854106506443962138469812371142769643e6,2.61321424050856457123433663461415238822132698699274898698444195581028556256292204118246298184152640660180308848e5,-1.86315094230899060297422858598853271269991701222133662269812151190073738494135161781746541881517543444088517450e4,6.72898907328811639584664793087658734023965461622305086223114830302184472501053920325357705077659100329375061164e2,-1.00943157793388998698906942778955472820868316720216759533786330365527456664479196396376802036991852763615117162e1,4.38168530207585179518852579327326813565547855423827923455725302731661798688344684949123570662576923314956409652e-2,-2.56700022740805144017893109688146755394089853324102182680986568703519793504112567874052887016745671988158860509e-5,2.36957641603015216306634255574107168474578275256678299025812704092568833214425267053473328220567647312858972842e-10,1.15694083038124833573384461491621092697227999368580657280434988874851263898593841144273580910524072051004097355e-13,-5.65316256384788292348472745512356223934162359491454541731051455366797335011238891664213587003781560857089287039e-14,8.49474213006453212163254972956476450659399535047675685153282680706022067567488560825213520889992328369290274917e-15,2.28070446151983292921304376186430296779950069002580275483997708509382375273278926321059842672372010645544369814e-15,-1.10200839918603550819509828894654530385656333624024070792512173994075523430080371272123208796268388120390605999e-15,1.38549328695392343672202630014049296547598488298392633928436570519822436737509703285623296516943877681573658856e-16]

-- g = 16
-- cs = [1.00000000000000000000000000000197282505854632603079830840695425875357874385280180091442294078665440768378249098e0,8.46526161663305542149634069187867594722046855412304900854881988813185784117115334278649283694595520258911681345e6,-4.36749138271579683480136186837344589623877333072263981316925882855709622493593700542050638104932540688970449938e7,9.74293930563290544819101625509993544487332726502764966578661696058106791782596806612874880755827781636515874517e7,-1.23204021755002536116556671349951637042719838766566833809791162301292663164161640991555755772951701480576547648e8,9.73252768070985015486610374238589008544120710425762737932789956744946587088886115434931192302762714053106117041e7,-4.99300527228931690175325509502458107530242214409813242610195665492694332478957451089840946485109879221508006660e7,1.67717577881504247993599407438579716255280180251024655266088962746483984030648719508692194681618566416623305775e7,-3.63600580914162233698153144687999173036255278366266794893601829961188765108175483037322789095070584870495716554e6,4.90477280895675364498279606024677805104335312501788362116122862426761656206254084389332939810112070827212068184e5,-3.86134008421465254942234421187656226690536946890041785481282535098425957705374976942206419371390945617575261325e4,1.59799318597802486242335864923532052679307039395883719223156368291075399965515365715900128098360767487863914335e3,-2.92431010362234704149212047766530662168493540855402650510149777639810317548261822822398843227566686564791963281e1,1.74363004063413808172424300916410794128403362491029945418010113517139134404174529869116082145355711405933435336e-1,-1.83889135553342355036509760427653033829857278652806446898154652952072080208456364635472214925432683895177618165e-4,7.33270898586375274415891722892287635561670631889884007144730226480410116396830758987236380519411854625863157502e-9,6.45278000390671865545527420895219728252383858385036218254715732532367519947191407848574538847659484635436688368e-14,-6.94479620536106946656438481553865602897856177541063529072614183743481416853611141761279751208940833601449239078e-14,4.05844190171647299618790214616008990498764985510143852349397107297416322977567534872136697380793297187751912405e-14,-1.31733313666599555356105188113267668954033474011729096786502800127065887481914541186166120448932913511034917752e-14,2.41974657267717455470266193816400294930708183176527991681763060922805701926472008462458376871397921380772169745e-15,-2.02601149216618151544716348082831639958354755920147356248364247315675789908515062799667823982813959353893675328e-16]

g = 16
cs = [1.00000000000000000000000000028760575603001265508442130506014379513956597321465491452116725029663713017318281155302254025482646249861597225426258553986702198310282111660814e0,8.46526161663305542149634069187868312963789113165756291473121826624097304045877402442031648610441688439038240651491328941139813852264197305778456591548908224120094894167246e6,-4.36749138271579683480136186837373556118742078558640362066413127517940638355237243621032816244371668677072024055467961309289136667985789351821077916283681043429310492245048e7,9.74293930563290544819101625512895178907295253257275119766847678373728139529380141034950496249773443820391664483106071950636988772945064231618026141173313776653501642702854e7,-1.23204021755002536116556671362722479624891929606934016862583690470559737742517159607918386688863050322288324639470143518965122467234214549961896872832625644312760192739883e8,9.73252768070985015486610377348101550835843155972616405114822066847313715073349332685175985019283585094322001442327102841230844094173203436533130957487267133091342982903057e7,-4.99300527228931690175325556901821725603229843211023737041708258158778148291264448338140176495149974138608744659105384002613928150272910583070070974922893028621993200088839e7,1.67717577881504247993599895442617443938121576155949723063296052767997809323676003179411440697229440626328003480936049898390416473418451439991990895198019408921657042045389e7,-3.63600580914162233698188824414224911677838005543078337270028990910418595747251549012507094099291095605470928507660671913700304145005691137775901293400972065312587564986433e6,4.90477280895675364500197057167325367117745996005366347188604951009746418243949483785855713959732995794028978290947548233313916565466830346275333882723371675976425249785277e5,-3.86134008421465255019809114269109448642105844599649825606338160546076754535370306095846501745502899782090098798450780912185670412952137486962203827925941120665671522549534e4,1.59799318597802488643745988279123446913929676490559975824444279070965368412618612119375002635982295898117052061965727441889307408959072186532547684742798077234919121659107e3,-2.92431010362235278765051929750170216644701333755426298013832714882892025277018457867985013916936895806752284909319872665280990487229040647950671740315052814101554630827301e1,1.74363004063520607047510940102445356013586745700563922221376454404331547889244236047131697225586656833395848432904781434760496518195945993246362481962863428301094957069334e-1,-1.83889135707484897639922506331631327352704695134176610912010000427967001903568852783866323303740089259923460414793348580157969596224404574122968841894479671002508273424642e-4,7.33288074820458553070627189736492168841920567015366487065477349916549899486012064868307264427889502768168730777757015382371906002557890047275383537043961908550989186406985e-9,-8.14003545649332564829862837358721864946644372854188383148183234349174404580929791883572722401611593795796287442225059328665383671935403293667079588128976652734012111434247e-14,2.30444691867487566543942418193243246136284662851351267620236991760404161776920402164496692134690158656181354762481193169147273885408030319356754409112626454347965677288216e-14,-1.60281647995499836731602960261827820681271025273058699260173363290101149359291286878982917703312074881795295123893833653667997796285463313440572145476822619163620650773297e-15,-1.51961151263436286650314929591024470402675431137010398063941973564941174782047280699515350734586316425940330788036418213975493171611105468630277933898275614760821849251065e-16]

foreign import ccall "math.h lgamma" lgamma :: Double -> Double
foreign import ccall "math.h tgamma" tgamma :: Double -> Double

err :: RealFrac a => (a -> a) -> a -> (a,a)
err approxGamma x = (realToFrac absErr, realToFrac relErr)
    where
        x' = realToFrac x
        
        y = realToFrac (approxGamma x)
        y' = refGamma' x'
        absErr = y - y'
        relErr = absErr / max (abs y) (abs y')