posit-2022.0.0.0: test/Test/Algorithms.hs
{-# LANGUAGE TypeApplications #-} -- To apply types: @Type, it seems to select the specific class instance, when GHC is not able to reason about things, commenting this out shows an interesting interface
{-# LANGUAGE ScopedTypeVariables #-} -- To reduce some code duplication, this is important
{-# LANGUAGE FlexibleContexts #-} -- to talk about class constraints like: (PositC es, PositC (Next es)) =>
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DataKinds #-} -- For our ES kind and the constructors Z, I, II, III, IV, V for exponent size type, post-pended with the version.
{-# OPTIONS_GHC -Wno-type-defaults #-} -- Turn off noise
{-# OPTIONS_GHC -Wno-unused-top-binds #-} -- Turn off noise
module Test.Algorithms
( funLogDomainReduction
, funLogTaylor
, funExp2
, funExpTaylor
, funLogTuma
, funExpTuma
, funGammaSeriesFused
, funGammaRamanujan
, funGammaCalc
, funGammaNemes
, funGammaYang
, funGammaChen
, funGammaXminus1
-- , funGammaViaLngamma
, funPi1
, funPi2
, funPi3
, funPi4
, funPi5
, funPi6
) where
import Posit -- run with -O_TEST CPP directive
import Prelude hiding (rem)
-- would like to:
-- import Posit.Internal.ElementaryFunctions
-- Perhaps on the chopping block if we are moving to ElementaryFunctions
-- Imports for implementing the Transcendental Functions
import GHC.Natural (Natural) -- Import the Natural Numbers ℕ (u+2115) for some of the Transcendental Functions
import Data.Ratio ((%)) -- Import the Rational Numbers ℚ (u+211A), ℚ can get arbitrarily close to Real numbers ℝ (u+211D), used for some of the Transcendental Functions
import Debug.Trace (trace) -- temporary for debug purposes
-- The machine implementation of the Posit encoding/decoding
import Posit.Internal.PositC -- The main internal implementation details
-- Algorithms in Type: `Posit es`
-- ==============================================================
-- Other functions:
-- ==============================================================
-- Approximation of log2 "Log Base 2"
approx_log2 :: forall es. PositC es => Posit es -> Posit es
approx_log2 NaR = NaR
approx_log2 z
| z <= 0 = NaR -- includes the NaR case
| otherwise = go (fromInteger ex) 1 sig -- domain reduction
where
go :: Posit es -> Posit es -> Posit es -> Posit es
go !acc !mak !sig' -- fixed point iteration, y is [1,2) :: Posit256
| sig == 1 = acc
| acc == (acc + mak * 2^^(negate.fst.term $ sig')) = acc -- stop when fixed point is reached
| otherwise = go (acc + mak * 2^^(negate.fst.term $ sig')) (mak * 2^^(negate.fst.term $ sig')) (snd.term $ sig')
term = findSquaring 0 -- returns (m,s') m the number of times to square, and the new significand
(ex, sig) = (int * fromIntegral (2^(exponentSize @es)) + fromIntegral nat, fromRational rat)
(_,int,nat,rat) = (posit2TupPosit @es).toRational $ z -- sign should always be positive
findSquaring m s
| s >= 2 && s < 4 = (m, s/2)
| otherwise = findSquaring (m+1) (s^2)
-- calculate atan(1/2^n)
-- sum k=0 to k=inf of the terms, iterate until a fixed point is reached
funArcTan :: Natural -> Posit256
funArcTan 0 = pi / 4
funArcTan n
| n <= 122 = go 0 0
| otherwise = z -- at small z... (atan z) == z "small angle approximation"
where
go !k !acc
| acc == (acc + term k) = acc
| otherwise = go (k+1) (acc + term k)
term :: Integer -> Posit256
term k = ((-1)^k * z^(2 * k + 1)) / fromIntegral (2 * k + 1)
z = 1 / 2^n -- recip $ 2^n :: Posit256 -- inv2PowN
-- fI2PN = (1 /) . (2 ^)
funInv2PowN :: Natural -> Posit256
funInv2PowN n = 1 / 2^n
-- calculate atanh(1/2^n)
-- sum k=0 to k=inf of the terms, iterate until a fixed point is reached
funArcHypTan :: Natural -> Posit256
funArcHypTan 0 = NaR
funArcHypTan n
| n <= 122 = go 0 0
| otherwise = z -- at small z... (atan z) == z "small angle approximation"
where
go !k !acc
| acc == (acc + term k) = acc
| otherwise = go (k+1) (acc + term k)
term :: Integer -> Posit256
term k = (z^(2 * k + 1)) / fromIntegral (2 * k + 1)
z = 1 / 2^n
--
funAtanhTaylor :: Posit256 -> Posit256
funAtanhTaylor NaR = NaR
funAtanhTaylor x
| abs x >= 1 = NaR
| abs x < 1/2^122 = x -- small angle approximaiton, found emperically
| x < 0 = negate.funAtanhTaylor.negate $ x
| otherwise = go 0 0
where
go !k !acc
| acc == (acc + term k) = acc
| otherwise = go (k+1) (acc + term k)
term :: Integer -> Posit256
term k = (x^(2 * k + 1)) / fromIntegral (2 * k + 1)
--
-- Taylor series expansion and fixed point algorithm, most accurate near zero
funSinTaylor :: Posit256 -> Posit256
funSinTaylor NaR = NaR
funSinTaylor z = go 0 0
where
go :: Natural -> Posit256 -> Posit256
go !k !acc
| acc == (acc + term k) = acc
| otherwise = go (k+1) (acc + term k)
term :: Natural -> Posit256
term k = (-1)^k * z^(2*k+1) / (fromIntegral.fac $ 2*k+1)
--
-- Taylor series expansion and fixed point algorithm, most accurate near zero
funCosTaylor :: Posit256 -> Posit256
funCosTaylor NaR = NaR
funCosTaylor z = go 0 0
where
go :: Natural -> Posit256 -> Posit256
go !k !acc
| acc == (acc + term k) = acc
| otherwise = go (k+1) (acc + term k)
term :: Natural -> Posit256
term k = (-1)^k * z^(2*k) / (fromIntegral.fac $ 2*k)
--
-- ~16 ULP for 42
funSinh :: Posit256 -> Posit256
funSinh NaR = NaR
funSinh x = (exp x - exp (negate x))/2
--
-- ~2 ULP for 42
funSinhTaylor :: Posit256 -> Posit256
funSinhTaylor NaR = NaR
funSinhTaylor z = go 0 0
where
go :: Natural -> Posit256 -> Posit256
go !k !acc
| acc == (acc + term k) = acc
| otherwise = go (k+1) (acc + term k)
term :: Natural -> Posit256
term k = z^(2*k+1) / (fromIntegral.fac $ 2*k+1)
--
--
funSinhTuma :: Posit256 -> Posit256
funSinhTuma NaR = NaR
funSinhTuma 0 = 0
funSinhTuma z | z < 0 = negate.funSinhTuma.negate $ z
funSinhTuma z | z > 80 = 0.5 * funExpTuma z
funSinhTuma z = go 256 1
where
go :: Natural -> Posit256 -> Posit256
go 1 !acc = z * acc
go !k !acc = go (k-1) (1 + (z^2 / fromIntegral ((2*k-2) * (2*k-1))) * acc)
--
-- ~17 ULP for 42
funCosh :: Posit256 -> Posit256
funCosh NaR = NaR
funCosh x = (exp x + exp (negate x))/2
--
-- ~3 ULP for 42
funCoshTaylor :: Posit256 -> Posit256
funCoshTaylor NaR = NaR
funCoshTaylor z = go 0 0
where
go :: Natural -> Posit256 -> Posit256
go !k !acc
| acc == (acc + term k) = acc
| otherwise = go (k+1) (acc + term k)
term :: Natural -> Posit256
term k = z^(2*k) / (fromIntegral.fac $ 2*k)
--
--
funCoshTuma :: Posit256 -> Posit256
funCoshTuma NaR = NaR
funCoshTuma 0 = 1
funCoshTuma z | z < 0 = funCoshTuma.negate $ z
funCoshTuma z | z > 3 = 0.5 * (funExpTuma z + funExpTuma (negate z))
funCoshTuma z = go 20 1
where
go :: Natural -> Posit256 -> Posit256
go 1 !acc = acc
go !k !acc = go (k-1) (1 + (z^2 / fromIntegral ((2*k-3)*(2*k-2)))*acc)
--
{-
-- | 'phi' fixed point recursive algorithm,
approx_phi :: (PositC es) => Posit es -> Posit es
approx_phi px@(Posit x)
| x == x' = Posit x
| otherwise = approx_phi (Posit x')
where
(Posit x') = (px^2 + 2*px) / (px^2 + 1)
-- LiquidHaskell is telling me this is unsafe if px is imaginary
-- lucky for us Posit256 is not imaginary
-}
--
-- Some series don't converge reliably, this one does
funLnOf2 :: Posit256
funLnOf2 = go 1 0
where
go :: Natural -> Posit256 -> Posit256
go !k !acc
| acc == (acc + term k) = acc
| otherwise = go (k+1) (acc + term k)
term :: Natural -> Posit256
term k = 1 / fromIntegral (2^k * k)
--
--
-- Gauss–Legendre algorithm, Seems only accurate to 2-3 ULP, but really slow
funPi1 :: forall es. (PositC es, PositC (Next es)) => Posit es
funPi1 = go 0 3 1 (recip.sqrt $ 2) (recip 4) 1
where
go :: Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es
go !prev !next !a !b !t !p
| prev == next = next
| otherwise =
let a' = (a + b) / 2
b' = sqrt $ a * b
t' = t - p * (a - ((a + b) / 2))^2
p' = 2 * p
in go next ((a' + b')^2 / (4 * t')) a' b' t' p'
--
-- Borwein's algorithm, with quintic convergence,
-- gets to 7 ULP in 4 iterations, but really slow due to expensive function evaluations
-- quite unstable and will not converge if sqrt is not accurate, which means log must be accurate
funPi2 :: forall es. (PositC es, PositC (Next es)) => Posit es
funPi2 = recip $ go 0 0 0 0.5 (5 / phi^3)
where
go :: Posit es -> Posit es -> Natural -> Posit es -> Posit es -> Posit es
go !prevA !prevS !n !a !s
| prevA == a = a
| prevS == s = a
| abs (prevA - a) < eps = a -- P256 will not reach a fixed point where `prevA == a` it sort of oscelates until divergence occurs, if we test for less than eps it can stop early
| otherwise =
let x = 5 / s - 1
y = (x - 1)^2 + 7
z = (0.5 * x * (y + sqrt (y^2 - 4 * x^3)))**(1/5)
a' = s^2 * a - (5^n * ((s^2 - 5)/2 + sqrt (s * (s^2 - 2*s + 5))))
s' = 25 / ((z + x/z + 1)^2 * s)
in go a s (n+1) (trace ("ΔA: " ++ show (a' - a)) a') (trace ("ΔS: " ++ show (s' - s)) s')
--
-- Bailey–Borwein–Plouffe (BBP) formula, to 1-2 ULP, and blazing fast, converges in 60 iterations
funPi3 :: forall es. (PositC es) => Posit es
funPi3 = go 0 0
where
go :: Integer -> Posit es -> Posit es
go !k !acc
| acc == acc + term k = acc
| otherwise = go (k+1) (acc + term k)
term :: Integer -> Posit es
term k = fromRational $ (1 % 16^k) * ((120 * k^2 + 151 * k + 47) % (512 * k^4 + 1024 * k^3 + 712 * k^2 + 194 * k + 15))
--
-- Fabrice Bellard improvement on the BBP, 2-3 ULP, even faster, converges in 25 iterations, really fast
funPi4 :: forall es. (PositC es) => Posit es
funPi4 = (1/2^6) * go 0 0
where
go :: Integer -> Posit es -> Posit es
go !k !acc
| acc == acc + term k = acc
| otherwise = go (k+1) (acc + term k)
term :: Integer -> Posit es
term k = fromRational $ ((-1)^k % (2^(10*k))) * ((1 % (10 * k + 9)) - (2^2 % (10 * k + 7)) - (2^2 % (10 * k + 5)) - (2^6 % (10 * k + 3)) + (2^8 % (10 * k + 1)) - (1 % (4 * k + 3)) - (2^5 % (4 * k + 1)))
--
-- Borwin's Quadradic Alogrithm 1985
funPi5 :: forall es. (PositC es, PositC (Next es)) => Posit es
funPi5 = recip $ go 0 0 1 (6 - 4 * sqrt 2) (sqrt 2 - 1)
where
go :: Posit es -> Posit es -> Natural -> Posit es -> Posit es -> Posit es
go !prevA !prevY !n a y
| prevA == a = a
| prevY == y = a
| otherwise =
let f = (1 - y^4)**(1/4)
y' = (1 - f) / (1 + f)
a' = a * (1 + y')^4 - 2^(2 * n + 1) * y' * (1 + y' + y'^2)
in if n == 3
then a'
else go a y (n+1) (trace ("A: " ++ show a') a') (trace ("Y: " ++ show y') y')
--
-- 3.14159265358979323846264338327950288419716939937510582097494459231
-- ULP: -97
-- Borwin's Cubic Algirthm
funPi6 :: forall es. (PositC es, PositC (Next es)) => Posit es
funPi6 = recip $ go 0 0 1 (1/3) ((sqrt 3 - 1) / 2)
where
go :: Posit es -> Posit es -> Natural -> Posit es -> Posit es -> Posit es
go !prevA !prevS !n !a !s
| prevA == a = a
| prevS == s = a
| otherwise =
let r = 3 / (1 + 2 * (1 - s^3)**(1/3))
s'= (r - 1) / 2
a'= r^2 * a - 3^(n-1) * (r^2 - 1)
in if n == 4
then a'
else go a s (n+1) a' s'
-- 3.14159265358979323846264338327950288419716939937510582097494459231
-- ULP: 216
--
--
--
-- calculate exp, its most accurate near zero
-- use the Nested Series of Jan J Tuma
funExpTuma :: Posit256 -> Posit256
funExpTuma NaR = NaR
funExpTuma 0 = 1
funExpTuma z = go 57 1 -- was 66
where
go :: Natural -> Posit256 -> Posit256
go !k !acc
| k == 0 = acc
| otherwise = go (k-1) (1 + (z / fromIntegral k) * acc)
--
--
-- Interestingly enough, wikipedia defines two alternative solutions
-- for the Shannon Wavelet, eventhough there are infinite solutions
-- where the functions are equal, they are not equal. It a class of
-- functions with the charicteristic of being a band pass filter in the
-- frequency space.
-- Shannon wavelet
funPsiSha1 :: Posit256 -> Posit256
funPsiSha1 NaR = NaR
funPsiSha1 t = 2 * sinc (2 * t) - sinc t
--
-- Shannon wavelet
funPsiSha2 :: Posit256 -> Posit256
funPsiSha2 NaR = NaR
funPsiSha2 t = sinc (t/2) * cos (3*pi*t/2)
--
-- Shannon wavelet, same as funPsiSha1 but with a factor of pi, with the
-- Law: funPsiSha1.(pi*) === funPsiSha3
-- or : funPsiSha1 === funpsiSha3.(/pi)
-- Posit256 seems to hold to a few ULP
funPsiSha3 :: Posit256 -> Posit256
funPsiSha3 NaR = NaR
funPsiSha3 0 = 1 -- Why the hell not!
funPsiSha3 t =
let pit = pi * t
invpit = recip pit
in invpit * (sin (2 * pit) - sin pit)
--
--
-- Using the CORDIC domain reduction and some approximation function
funLogDomainReduction :: (Posit256 -> Posit256) -> Posit256 -> Posit256
funLogDomainReduction _ NaR = NaR
funLogDomainReduction _ 1 = 0
funLogDomainReduction f x
| x <= 0 = NaR
| otherwise = f sig + (fromIntegral ex * lnOf2)
where
(ex, sig) = (int * fromIntegral (2^(exponentSize @V_3_2)) + fromIntegral nat + 1, fromRational rat / 2) -- move significand range from 1,2 to 0.5,1
(_,int,nat,rat) = (posit2TupPosit @V_3_2).toRational $ x -- sign should always be positive
--
-- Use the constant, for performance
lnOf2 :: PositC es => Posit es
lnOf2 = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202
--
--
-- calculate exp, its most accurate near zero
-- sum k=0 to k=inf of the terms, iterate until a fixed point is reached
funExpTaylor :: Posit256 -> Posit256
funExpTaylor NaR = NaR
funExpTaylor 0 = 1
funExpTaylor z = go 0 0
where
go :: Natural -> Posit256 -> Posit256
go !k !acc
| acc == (acc + term k) = acc -- if x == x + dx then terminate and return x
| otherwise = go (k+1) (acc + term k)
term :: Natural -> Posit256
term k = (z^k) / (fromIntegral.fac $ k)
--
--
--
funExp2 :: (Posit256 -> Posit256) -> Posit256 -> Posit256
funExp2 _ NaR = NaR
funExp2 _ 0 = 1
funExp2 f x
| x < 0 = recip.funExp2 f.negate $ x -- always calculate the positive method
| otherwise = case properFraction x of
(int,rem) -> fromIntegral (2^int) * f (lnOf2 * rem)
funGammaSeriesFused :: forall es. (PositC es, PositC (Next es)) => Posit es -> Posit es
funGammaSeriesFused z = sqrt(2 * pi) * (z**(z - 0.5)) * exp (negate z) * (1 + series)
where
series :: Posit es
series = fsumL $ zipWith (*) [fromRational (a % b) | (a,b) <- zip a001163 a001164] [recip $ z^n | n <- [1..len]] -- zipWith (\x y -> ) a001163 a001164
lenA = length a001163
lenB = length a001164
len = if lenA == lenB
then lenA
else error "Seiries Numerator and Denominator do not have the same length."
--
--
a001163 :: [Integer] -- Numerator
a001163 = [1, 1, -139, -571, 163879, 5246819, -534703531, -4483131259, 432261921612371, 6232523202521089, -25834629665134204969, -1579029138854919086429, 746590869962651602203151, 1511513601028097903631961, -8849272268392873147705987190261, -142801712490607530608130701097701]
a001164 :: [Integer] -- Denominator
a001164 = [12, 288, 51840, 2488320, 209018880, 75246796800, 902961561600, 86684309913600, 514904800886784000, 86504006548979712000, 13494625021640835072000, 9716130015581401251840000, 116593560186976815022080000, 2798245444487443560529920000, 299692087104605205332754432000000, 57540880724084199423888850944000000]
--
-- natural log with log phi acurate to 9 ULP
funLogTaylor :: Posit256 -> Posit256
funLogTaylor NaR = NaR
funLogTaylor 1 = 0
funLogTaylor x | x <= 0 = NaR
funLogTaylor x
| x <= 2 = go 1 0
| otherwise = error "The funLogTaylor algorithm is being used improperly"
where
go :: Natural -> Posit256 -> Posit256
go !k !acc
| acc == (acc + term k) = acc
| otherwise = go (k + 1) (acc + term k)
term :: Natural -> Posit256
term k = (-1)^(k+1) * (x - 1)^k / fromIntegral k
-- natural log the Jan J Tuma way
funLogTuma :: Posit256 -> Posit256
funLogTuma NaR = NaR
funLogTuma 1 = 0 -- domain reduced input is [0.5,1) and/or , where funLogTuma 1 = 0
funLogTuma x | x <= 0 = NaR -- zero and less than zero is NaR
funLogTuma x
= go 242 1
where
xM1 = x - 1 -- now [-0.5, 0)
go :: Natural -> Posit256 -> Posit256
go !k !acc
| k == 0 = xM1 * acc
| otherwise = go (k-1) (recip (fromIntegral k) - xM1 * acc)
--
--
funGammaRamanujan :: (PositC es, PositC (Next es)) => Posit es -> Posit es
funGammaRamanujan z = sqrt pi * (x / exp 1)**x * (8*x^3 + 4*x^2 + x + (1/30))**(1/6)
where
x = z - 1
--
--
funGammaCalc :: (PositC es, PositC (Next es)) => Posit es -> Posit es
funGammaCalc z = sqrt (2*pi / z) * ((z / exp 1) * sqrt (z * sinh (recip z) + recip (810 * z^6)))**z
funGammaNemes :: (PositC es, PositC (Next es)) => Posit es -> Posit es
funGammaNemes z = sqrt (2*pi / z) * (recip (exp 1) * (z + recip (12 * z - recip (10 * z))))**z
funGammaYang :: (PositC es, PositC (Next es)) => Posit es -> Posit es
funGammaYang z = sqrt (2 * pi * x) * (x / exp 1)**x * (x * sinh (recip x))**(x/2) * exp (fromRational (7 % 324) * recip (x^3 * (35 * x^2 + 33)))
where
x = z - 1
funGammaChen :: (PositC es, PositC (Next es)) => Posit es -> Posit es
funGammaChen z = sqrt (2 * pi * x) * (x / exp 1)**x * (1 + recip (12*x^3 + (24/7)*x - 0.5))**(x^2 + fromRational (53 % 210))
where
x = z - 1
funGammaXminus1 :: (PositC es, PositC (Next es)) => Posit es -> Posit es
funGammaXminus1 x = go (x - 1)
where
go z = sqrt (2 * pi) * exp z ** (negate z) * z ** (z + 0.5)
{-
funGammaInfProd :: Posit es -> Posit es
funGammaInfProd
funGammaViaInv :: Posit es -> Posit es
funGammaViaInv
-}
{-
funGammaViaLngamma :: forall es. (PositC es, PositC (Next es)) => Posit es -> Posit es
funGammaViaLngamma z = exp $ lngamma
where
lngamma :: Posit es
lngamma = negate eulersConstant * z - log z + go 0 1
go :: Posit es -> Integer -> Posit es
go NaR _ = NaR
go prev k | prev == prev + next k = prev
| otherwise = go (trace ("Next: " ++ show (prev + next k)) (prev + next k)) (k+1)
next :: Integer -> Posit es
next k = z / fromIntegral k - (log $ 1 + z / fromIntegral k)
eulersConstant :: PositC es => Posit es
eulersConstant = 0.57721566490153286060651209008240243104215933593992
-}
fac :: Natural -> Natural
fac 0 = 1
fac n = n * fac (n - 1)
--