--------------------------------------------------------------------------------------------
--   Posit Numbers
--   Copyright   :  (C) 2022 Nathan Waivio
--   License     :  BSD3
--   Maintainer  :  Nathan Waivio <nathan.waivio@gmail.com>
--   Stability   :  Stable
--   Portability :  Portable
--
-- | Library implementing standard Posit Numbers (Posit Standard version
--   3.2.0.0, with some improvements) a fixed width word size of
--   2^es bytes.
-- 
---------------------------------------------------------------------------------------------


{-# LANGUAGE GADTs #-} --   For our main type Posit (es :: ES)
{-# LANGUAGE DataKinds #-}  --   For our ES kind and the constructors Z, I, II, III, IV, V for exponent size type
{-# LANGUAGE KindSignatures #-}  --   For defining the type of kind ES that indexes the GADT
{-# LANGUAGE ViewPatterns #-}  --   To decode the posit in the pattern
{-# LANGUAGE BangPatterns #-}  --   Added Strictness for some fixed point algorithms
{-# LANGUAGE PatternSynonyms #-}  --   for a nice NaR interface
{-# LANGUAGE FlexibleInstances #-} --   To make instances for each specific type [Posit8 .. Posit256]
{-# LANGUAGE FlexibleContexts #-}  --   Allow non-type variables in the constraints
{-# 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 MultiParamTypeClasses #-}  --   To convert between Posit Types
{-# LANGUAGE ScopedTypeVariables #-} --   To reduce some code duplication
{-# LANGUAGE UndecidableInstances #-}  --   To reduce some code duplication, I think the code is decidable but GHC is not smart enough ;), like there being only 1 instance that is polymorphic and works for all of my types.
{-# LANGUAGE CPP #-} --   To remove Storable instances to remove noise when performing analysis of Core
{-# OPTIONS_GHC -Wno-unticked-promoted-constructors #-}  --   Turn off noise
{-# OPTIONS_GHC -Wno-type-defaults #-}  --   Turn off noise
{-# OPTIONS_GHC -Wno-unused-top-binds #-}  --   Turn off noise

-- ----
--  Posit numbers implementing:
--
--    * Show
--    * Eq
--    * Ord  -- compare as an integer representation
--    * Num  -- Addition, subtraction, multiplication, and other operations
--    * Enum  -- Successor and Predecessor
--    * Fractional  -- division, divide by zero is Not a Real (NaR) number
--    * Real
--    * Bounded
--    * FusedOps  -- dot product and others
--    * Convertible  -- Conversions between different posit formats
--    * AltShow
--    * Read
--    * Storable  -- Formats for binary data, for computation and data interchange
--    * RealFrac
--    * RealFloat
--    * Floating  -- Mathematical functions such as logarithm, exponential, trigonometric, and hyperbolic functions. Warning! May induce trance.
--
-- ----

module Posit
(-- * Main Exported Types
 Posit8, -- |An 8-bit Posit number with 'es' ~ 'Z'
 Posit16, -- |An 16-bit Posit number with 'es' ~ 'I'
 Posit32, -- |An 32-bit Posit number with 'es' ~ 'II'
 Posit64, -- |An 64-bit Posit number with 'es' ~ 'III'
 Posit128, -- |An 128-bit Posit number with 'es' ~ 'IV'
 Posit256, -- |An 256-bit Posit number with 'es' ~ 'V'
 
 -- * Patterns for Matching Exported Types
 pattern NaR,  -- |A pattern for Exception handling when a value is Not a Real number (NaR).
 pattern R,  -- |A pattern for the non-Exceptional case, yielding a Rational, will make a total function when paired with NaR, if the Rational implementation is total.
 
 -- * Fused Operation Interface defined by the Posit Standard
 FusedOps(..),
 
 -- * Posits are Convertable between different Posit representations
 Convertible(..),
 
 -- * Additional functions to show the Posit in different formats
 AltShow(..),
 
 -- * Additional Special Functions
 AltFloating(..),
 
 -- * Functions to lift functions of Integers or Rationals to operate on Posit Types
 viaIntegral,
 viaRational,
 viaRational2,
 viaRational3,
 viaRational4,
 viaRational6,
 viaRational8,
#ifdef O_TEST
 -- * Alternative algorithms for test purposes
 funExp,
 funExp2,
 funExpTaylor,
 funLogTaylor,
 funExpTuma,
 funLogTuma,
 funLogDomainReduction,
 funPi1,
 funPi2,
 funPi3,
 funPi4,
 funPsiSha1,
 funPsiSha2,
 funPsiSha3
#endif
 ) where


import Prelude hiding (rem)

-- Imports for Show and Read Instances
import Data.Scientific (Scientific
                       ,scientificP
                       ,fromRationalRepetendUnlimited
                       ,formatScientific
                       ,FPFormat(Generic)) -- Used to print/show and read the rational value

import Text.Read (Lexeme(Ident)
                 ,readPrec
                 ,readListPrec
                 ,(+++)
                 ,pfail
                 ,readListPrecDefault
                 ,lexP
                 ,lift
                 ,parens) -- Used to read a Posit value

-- Imports for Vectorization Class Instances
import Data.Foldable (toList)  -- Used for fused operations on foldable/lists

-- Imports for Storable Instance
import Foreign.Storable (Storable, sizeOf, alignment, peek, poke)  -- Used for Storable Instances of Posit
import Foreign.Ptr (Ptr, plusPtr, castPtr)  -- Used for dealing with Pointers for the Posit Storable Instance


-- 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 (Rational, (%))  -- 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


-- =====================================================================
-- ===                  Posit Implementation                         ===
-- =====================================================================

-- The machine implementation of the Posit encoding/decoding
import Posit.Internal.PositC (ES(..), PositC(..))  -- The main internal implementation details


-- |Base GADT rapper type, that uses the Exponent Size kind to index the various implementations
data Posit (es :: ES) where
     Posit :: PositC es => !(IntN es) -> Posit es

-- |Not a Real Number, the Posit is like a Maybe type, it's either a real number or not
pattern NaR :: (PositC es) => Posit es
pattern NaR <- (Posit (decode -> Nothing)) where
  NaR = Posit unReal
--

--
-- |A Real or at least Rational Number, rounded to the nearest Posit Rational representation
pattern R :: PositC es => Rational -> Posit es
pattern R r <- (Posit (decode -> Just r)) where
  R r = Posit (encode $ Just r)
--

-- Posit functions are complete if the following two patterns are completely defined.
{-# COMPLETE NaR, R #-}

-- Concrete types exported for use.
type Posit8 = Posit Z
type Posit16 = Posit I
type Posit32 = Posit II
type Posit64 = Posit III
type Posit128 = Posit IV
type Posit256 = Posit V

#ifndef O_NO_SHOW
-- Show
--
instance forall es. (PositC es) => Show (Posit es) where
  show NaR = "NaR"
  show (R r) = formatScientific Generic (Just $ decimalPrec @es) (fst.fromRationalRepetendUnlimited $ r)
--
#endif



-- Two Posit Numbers are Equal if their Finite Precision Integer representation is Equal
--
-- All things equal I would rather write it like this:
instance forall es. (Eq (IntN es)) => Eq (Posit es) where
  (Posit int1) == (Posit int2) = int1 == int2
--



-- Two Posit Numbers are ordered by their Finite Precision Integer representation
--
-- Ordinarily I would only like one instance to cover them all
instance forall es. (Ord (IntN es), PositC es) => Ord (Posit es) where
  compare (Posit int1) (Posit int2) = compare int1 int2
--



-- Num
--
-- I'm num trying to get this definition:
instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => Num (Posit es) where
  -- Addition
  (+) = viaRational2 (+)
  -- Multiplication
  (*) = viaRational2 (*)
  -- 'abs', Absolute Value, it's like a magnitude of sorts, abs of a posit is the same as abs of the integer representation
  abs = viaIntegral abs
  -- 'signum' it is a kind of an representation of directionality, the sign of a number for instance
  signum = viaRational signum
  -- 'fromInteger' rounds the integer into the closest posit number
  fromInteger int = Posit $ encode (Just $ fromInteger int)
  -- 'negate', Negates the sign of the directionality. negate of a posit is the same as negate of the integer representation
  negate = viaIntegral negate
--

-- deriving via Integral Class, for the Integral representation of the posit
viaIntegral :: forall es. PositC es => (IntN es -> IntN es) -> Posit es -> Posit es
viaIntegral f (Posit int) = Posit $ f int
--



-- Enum-ish, A Posit has a Successor and Predecessor so its an ordinal number, as per Posit standard next, prior
-- The Posit Standard requires 2's complement integer overflow to be ignored
instance forall es. (Num (IntN es), Enum (IntN es), Ord (IntN es), PositC es) => Enum (Posit es) where
  -- succ (Posit int) = Posit (int + 1)
  succ = viaIntegral (+1)
  -- succ = viaIntegral succ  -- Non-compliant, runtime error pred NaR, and worse it is Int64 for types of greater precision, probably because of Preludes gross abomination of toEnum/fromEnum
  -- pred (Posit int) = Posit (int - 1)
  pred = viaIntegral (subtract 1)
  -- pred = viaIntegral pred  -- Non-compliant, runtime error pred NaR, and worse it is Int64 for types of greater precision, probably because of Preludes gross abomination of toEnum/fromEnum
  -- enumFrom :: Posit es -> [Posit es]
  enumFrom n = enumFromTo n maxBound
  enumFromTo n m
    | n == m = [n]
    | n < m = n : enumFromTo (succ n) m
    | otherwise = []
  -- enumFromThen n m :: Posit es -> Posit es -> [Posit es]
  enumFromThen NaR _ = [NaR]
  enumFromThen _ NaR = [NaR]
  enumFromThen n m = n : go n
    where
      step = m - n
      go :: Posit es -> [Posit es]
      go NaR = [NaR]
      go !l = case compare step 0 of
                LT -> let !n' = l + step  -- rounding occurs here, because the next comparison needs it, it wouldn't make sense otherwise...
                      in if n' - l > step
                         then []
                         else n' : go n'
                EQ -> [n, m]
                GT -> let !n' = l + step
                      in if n' - l < step
                         then []  -- with tapered resolution this algorithm can reach a fixed point where the next value is equal to the previous value
                         else n' : go n'
  enumFromThenTo NaR  _   _  = [NaR]
  enumFromThenTo  _  NaR  _  = [NaR]
  enumFromThenTo  _   _  NaR = [NaR]
  enumFromThenTo  e1  e2  e3 = takeWhile predicate (enumFromThen e1 e2)
    where
      mid = (e2 - e1) / 2
      predicate | e2 >= e1  = (<= e3 + mid)
                | otherwise = (>= e3 + mid)
--



-- Fractional Instances; (Num => Fractional)
--
-- How the Frac do I get this definition:
instance forall es. (Num (IntN es), Ord (IntN es), Eq (IntN es), PositC es) => Fractional (Posit es) where
  fromRational = R
 
  recip 0 = NaR
  recip p = viaRational recip p
--

-- Rational Instances; Num & Ord Instanced => Real
--
-- I for real want this definition:
instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => Real (Posit es) where
  toRational NaR = error "Your input is Not a Real or Rational (NaR) number, please try again!"
  toRational (R r) = r
--

-- Implementing instances via Rational Data Type's instance,
-- The function checks for NaR, to protect against the runtime error 'toRational' would generate if called with a NaR value
-- Unary::Arity NaR guarded pass through with wrapping and unwrapping use of a Rational function
viaRational :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational) -> Posit es -> Posit es
viaRational _ NaR = NaR
viaRational f (R r) = fromRational $ f r

-- Binary NaR guarded pass through with wrapping and unwrapping use of a Rational function
viaRational2 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es
viaRational2 _ NaR  _  = NaR
viaRational2 _  _  NaR = NaR
viaRational2 f (R r1) (R r2) = R $ r1 `f` r2

-- Ternary NaR guarded pass through with wrapping and unwrapping use of a Rational function
viaRational3 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es
viaRational3 _ NaR  _   _  = NaR
viaRational3 _  _  NaR  _  = NaR
viaRational3 _  _   _  NaR = NaR
viaRational3 f (R r1) (R r2) (R r3) = R $ f r1 r2 r3

-- Quaternary NaR guarded pass through with wrapping and unwrapping use of a Rational function
viaRational4 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es
viaRational4 _ NaR  _   _   _  = NaR
viaRational4 _  _  NaR  _   _  = NaR
viaRational4 _  _   _  NaR  _  = NaR
viaRational4 _  _   _   _  NaR = NaR
viaRational4 f (R r0) (R r1) (R r2) (R r3) = R $ f r0 r1 r2 r3

-- Senary NaR guarded pass through with wrapping and unwrapping use of a Rational function
viaRational6 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es
viaRational6 _ NaR  _   _   _   _   _  = NaR
viaRational6 _  _  NaR  _   _   _   _  = NaR
viaRational6 _  _   _  NaR  _   _   _  = NaR
viaRational6 _  _   _   _  NaR  _   _  = NaR
viaRational6 _  _   _   _   _  NaR  _  = NaR
viaRational6 _  _   _   _   _   _  NaR = NaR
viaRational6 f (R a1) (R a2) (R a3) (R b1) (R b2) (R b3) = R $ f a1 a2 a3 b1 b2 b3

-- Octonary NaR guarded pass through with wrapping and unwrapping use of a Rational function
viaRational8 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es
viaRational8 _ NaR  _   _   _   _   _   _   _  = NaR
viaRational8 _  _  NaR  _   _   _   _   _   _  = NaR
viaRational8 _  _   _  NaR  _   _   _   _   _  = NaR
viaRational8 _  _   _   _  NaR  _   _   _   _  = NaR
viaRational8 _  _   _   _   _  NaR  _   _   _  = NaR
viaRational8 _  _   _   _   _   _  NaR  _   _  = NaR
viaRational8 _  _   _   _   _   _   _  NaR  _  = NaR
viaRational8 _  _   _   _   _   _   _   _  NaR = NaR
viaRational8 f (R a0) (R a1) (R a2) (R a3) (R b0) (R b1) (R b2) (R b3) = R $ f a0 a1 a2 a3 b0 b1 b2 b3



-- Bounded, bounded to what?!? To the ℝ! NaR is out of bounds!!!
--
-- I'm bound to want this definition:
instance forall es. PositC es => Bounded (Posit es) where
  -- 'minBound' the most negative number represented
  minBound = Posit mostNegVal
  -- 'maxBound' the most positive number represented
  maxBound = Posit mostPosVal
--


-- =====================================================================
-- ===                    Fused Operations                           ===
-- =====================================================================

-- |A class that delays the rounding operation until the end for some operations
class Num a => FusedOps a where
  -- |Fused Multiply Add: (a * b) + c
  fma :: a -> a -> a -> a
  -- |Fused Add Multiply: (a + b) * c
  fam :: a -> a -> a -> a
  -- |Fused Multiply Multiply Subtract: (a * b) - (c * d)
  fmms :: a -> a -> a -> a -> a
  -- |Fused Sum of 3 values: a + b + c
  fsum3 :: a -> a -> a -> a
  -- |Fused Sum of 4 values: a + b + c + d
  fsum4 :: a -> a -> a -> a -> a
  -- |Fused Sum of a List of Posits
  fsumL :: Foldable t => t a -> a
  -- |Fused Dot Product of 3 element vector: (a1 * b1) + (a2 * b2) + (a3 * b3)
  fdot3 :: a -> a -> a -> a -> a -> a -> a
  -- |Fused Dot Product of 4 element vector: (a0 * b0) + (a1 * b1) + (a2 * b2) + (a3 * b3)
  fdot4 :: a -> a -> a -> a -> a -> a -> a -> a -> a
  -- |Fused Dot Product of Two Lists
  fdotL :: Foldable t => t a -> t a -> a
  -- |Fused Subtract Multiply: a - (b * c)
  fsm :: a -> a -> a -> a
 


-- Rational Instance
instance FusedOps Rational where
  fsm a b c = a - (b * c)
  fma a b c = (a * b) + c
  fam a b c = (a + b) * c
  fmms a b c d = (a * b) - (c * d)
  fsum3 a b c = a + b + c
  fsum4 a b c d = a + b + c + d
  fsumL (toList -> l) = go l 0
    where
      go [] acc = acc
      go (x : xs) acc = go xs (acc + x)
  fdot3 a1 a2 a3 b1 b2 b3 = (a1 * b1) + (a2 * b2) + (a3 * b3)
  fdot4 a0 a1 a2 a3 b0 b1 b2 b3 = (a0 * b0) + (a1 * b1) + (a2 * b2) + (a3 * b3)
  fdotL (toList -> l1) (toList -> l2) = go l1 l2 0
    where
      go [] [] acc = acc
      go []  _  _  = error "Lists not the same length"
      go _  []  _  = error "Lists not the same length"
      go (b : bs) (c : cs) acc = go bs cs (fma b c acc)
--

--
instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => FusedOps (Posit es) where
  -- Fused Subtract Multiply
  fsm = viaRational3 fsm
  -- Fuse Multiply Add
  fma = viaRational3 fma
  -- Fuse Add Multiply
  fam = viaRational3 fam
  -- Fuse Multiply Multiply Subtract
  fmms = viaRational4 fmms
  -- Fuse Sum of 3 Posits
  fsum3 = viaRational3 fsum3
  -- Fuse Sum of 4 Posits
  fsum4 = viaRational4 fsum4
  -- Fuse Sum of a List
  fsumL (toList -> l) = Posit $ encode (Just $ go l 0)
    where
      go :: [Posit es] -> Rational -> Rational
      go [] !acc = acc
      go ((Posit int) : xs) !acc = case decode int of
                                     Nothing -> error "Posit List contains NaR"
                                     Just r -> go xs (acc + r)
  -- Fuse Dot Product of a 3-Vector
  fdot3 = viaRational6 fdot3
  -- Fuse Dot Product of a 4-Vector
  fdot4 = viaRational8 fdot4
  -- Fuse Dot Product of two Lists
  fdotL (toList -> l1) (toList -> l2) = Posit $ encode (Just $ go l1 l2 0)
    where
      go [] [] !acc = acc
      go []  _   _  = error "Lists not the same length"
      go _  []   _  = error "Lists not the same length"
      go ((Posit int1) : bs) ((Posit int2) : cs) !acc = case decode int1 of
                                                          Nothing -> error "First Posit List contains NaR"
                                                          Just r1 -> case decode int2 of
                                                                       Nothing -> error "Second Posit List contains NaR"
                                                                       Just r2 -> go bs cs (acc + (r1 * r2))
--




-- =====================================================================
-- ===                  Conversion Between Posits Types              ===
-- =====================================================================

-- |A Convertible class that will cast or 'convert' between two different Posit es types
class Convertible a b where
  convert :: a -> b

instance forall es1 es2. (PositC es1, PositC es2, Ord (IntN es1), Ord (IntN es2), Num (IntN es1), Num (IntN es2)) => Convertible (Posit es1) (Posit es2) where
  convert NaR = NaR
  convert (R r) = R r
--


#ifndef O_NO_SHOW
-- =====================================================================
-- ===                Alternative Show Formats                       ===
-- =====================================================================

-- |A Alternative to the typical 'Show' class to assist in displaying the Posit es type in different formats
class AltShow a where
  -- |Display the Posit in its Binary Representation
  displayBinary :: a -> String
  -- |Display the Posit in its Integral Representation
  displayIntegral :: a -> String
  -- |Display the Posit as a Rational
  displayRational :: a -> String
  -- |Display the Posit as a Decimal until the Repetend occurs
  displayDecimal :: a -> String
--

--
instance forall es. (Show (IntN es), Ord (IntN es), Num (IntN es), PositC es) => AltShow (Posit es) where
  displayBinary (Posit int) = displayBin int
 
  displayIntegral (Posit int) = show int
 
  displayRational = viaShowable id
 
  displayDecimal = viaShowable (fst.fromRationalRepetendUnlimited)
--

viaShowable :: forall es a. (Show a, Ord (IntN es), Num (IntN es), PositC es) => (Rational -> a) -> Posit es -> String
viaShowable _ NaR = "NaR"
viaShowable f (R r) = show $ f r
#endif

#ifndef O_NO_READ
-- =====================================================================
-- ===                         Read Posit                            ===
-- =====================================================================

--
instance forall es. (PositC es) => Read (Posit es) where
  readPrec =
    parens $ do
      x <- lexP
      case x of
        Ident "NaR" -> return NaR
        _ -> pfail
      +++
      do
        s <- lift scientificP
        return $ R (toRational s)
 
  readListPrec = readListPrecDefault
--
#endif


-- =====================================================================
-- ===                  Storable Instances                           ===
-- =====================================================================
--
#ifndef O_NO_STORABLE
--
instance forall es. (Storable (IntN es), PositC es) => Storable (Posit es) where
  sizeOf _ = fromIntegral $ nBytes @es
  alignment _ = fromIntegral $ nBytes @es
  peek ptr = do
    int <- peek (castPtr ptr :: Ptr (IntN es))
    return $ Posit int
  poke ptr (Posit int) = do
    poke (castPtr ptr :: Ptr (IntN es)) int
--
#endif


-- =====================================================================
-- ===                        Real Frac                              ===
-- =====================================================================

--
instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => RealFrac (Posit es) where
  -- properFraction :: Integral b => a -> (b, a)
  properFraction = viaRationalErrTrunkation "NaR value is not a RealFrac" properFraction
--

viaRationalErrTrunkation :: forall es a. (Num (IntN es), (Ord (IntN es)), PositC es, Integral a) => String -> (Rational -> (a, Rational)) -> Posit es -> (a, Posit es)
viaRationalErrTrunkation err _ NaR = error err
viaRationalErrTrunkation _ f (R r) =
  let (int, r') = f r
  in (int, R r')

-- =====================================================================
-- ===                         Real Float                            ===
-- =====================================================================
--
instance forall es. (Eq (IntN es), Ord (IntN es), Num (IntN es), Floating (Posit es), PositC es) => RealFloat (Posit es) where
  isIEEE _ = False
  isDenormalized _ = False
  isNegativeZero _ = False
 
  isNaN NaR = True
  isNaN  _  = False
 
  isInfinite NaR = True
  isInfinite _ = False
 
  -- 'atan2' of y x is the argument "arg function" (also called phase or angle) of the complex number x + i y.
  -- angle from an x basis vector to some other vector
  --
  --     Y
  --     ^
  --     |    ^ (x,y)
  --     |   /
  --     |  / <-  alpha (radians)
  --     | /                      \
  --      /                        |
  --      -----------------------------------> X
  --
  --
  atan2 NaR  _  = NaR
  atan2  _  NaR = NaR
  atan2 y x
    | x == 0 && y == 0 = NaR
    | x > 0             = atan (y/x)
    | x < 0  && y >= 0  = atan (y/x) + pi
    | x < 0  && y  < 0  = atan (y/x) - pi
    | x == 0 && y  > 0  = pi / 2
    | x == 0 && y  < 0  = negate $ pi / 2
    | otherwise = error "What!?!?!" -- The case where x == 0 && y == 0
 
  floatRadix _ = 2
  floatDigits _ = undefined
  floatRange _ = (negate maxExponent, maxExponent)
    where
      maxExponent = fromIntegral $ (nBytes @es) * ((nBits @es) - 2)
  decodeFloat = undefined
  encodeFloat = undefined
--



-- =====================================================================
-- ===                         Floating                              ===
-- =====================================================================


instance Floating Posit8 where
  pi = convert (pi :: Posit256) :: Posit8
  exp x = convert (exp (convert x) :: Posit256) :: Posit8
  log x = convert (log (convert x) :: Posit256) :: Posit8
  x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit8
  sin x = convert (sin (convert x) :: Posit256) :: Posit8
  cos x = convert (cos (convert x) :: Posit256) :: Posit8
  asin x = convert (asin (convert x) :: Posit256) :: Posit8
  acos x = convert (acos (convert x) :: Posit256) :: Posit8
  atan x = convert (atan (convert x) :: Posit256) :: Posit8
  sinh x = convert (sinh (convert x) :: Posit256) :: Posit8
  cosh x = convert (cosh (convert x) :: Posit256) :: Posit8
  asinh x = convert (asinh (convert x) :: Posit256) :: Posit8
  acosh x = convert (acosh (convert x) :: Posit256) :: Posit8
  atanh x = convert (atanh (convert x) :: Posit256) :: Posit8

instance Floating Posit16 where
  pi = convert (pi :: Posit256) :: Posit16
  exp x = convert (exp (convert x) :: Posit256) :: Posit16
  log x = convert (log (convert x) :: Posit256) :: Posit16
  x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit16
  sin x = convert (sin (convert x) :: Posit256) :: Posit16
  cos x = convert (cos (convert x) :: Posit256) :: Posit16
  asin x = convert (asin (convert x) :: Posit256) :: Posit16
  acos x = convert (acos (convert x) :: Posit256) :: Posit16
  atan x = convert (atan (convert x) :: Posit256) :: Posit16
  sinh x = convert (sinh (convert x) :: Posit256) :: Posit16
  cosh x = convert (cosh (convert x) :: Posit256) :: Posit16
  asinh x = convert (asinh (convert x) :: Posit256) :: Posit16
  acosh x = convert (acosh (convert x) :: Posit256) :: Posit16
  atanh x = convert (atanh (convert x) :: Posit256) :: Posit16

instance Floating Posit32 where
  pi = convert (pi :: Posit256) :: Posit32
  exp x = convert (exp (convert x) :: Posit256) :: Posit32
  log x = convert (log (convert x) :: Posit256) :: Posit32
  x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit32
  sin x = convert (sin (convert x) :: Posit256) :: Posit32
  cos x = convert (cos (convert x) :: Posit256) :: Posit32
  asin x = convert (asin (convert x) :: Posit256) :: Posit32
  acos x = convert (acos (convert x) :: Posit256) :: Posit32
  atan x = convert (atan (convert x) :: Posit256) :: Posit32
  sinh x = convert (sinh (convert x) :: Posit256) :: Posit32
  cosh x = convert (cosh (convert x) :: Posit256) :: Posit32
  asinh x = convert (asinh (convert x) :: Posit256) :: Posit32
  acosh x = convert (acosh (convert x) :: Posit256) :: Posit32
  atanh x = convert (atanh (convert x) :: Posit256) :: Posit32

instance Floating Posit64 where
  pi = convert (pi :: Posit256) :: Posit64
  exp x = convert (exp (convert x) :: Posit256) :: Posit64
  log x = convert (log (convert x) :: Posit256) :: Posit64
  x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit64
  sin x = convert (sin (convert x) :: Posit256) :: Posit64
  cos x = convert (cos (convert x) :: Posit256) :: Posit64
  asin x = convert (asin (convert x) :: Posit256) :: Posit64
  acos x = convert (acos (convert x) :: Posit256) :: Posit64
  atan x = convert (atan (convert x) :: Posit256) :: Posit64
  sinh x = convert (sinh (convert x) :: Posit256) :: Posit64
  cosh x = convert (cosh (convert x) :: Posit256) :: Posit64
  asinh x = convert (asinh (convert x) :: Posit256) :: Posit64
  acosh x = convert (acosh (convert x) :: Posit256) :: Posit64
  atanh x = convert (atanh (convert x) :: Posit256) :: Posit64

instance Floating Posit128 where
  pi = convert (pi :: Posit256) :: Posit128
  exp x = convert (exp (convert x) :: Posit256) :: Posit128
  log x = convert (log (convert x) :: Posit256) :: Posit128
  x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit128
  sin x = convert (sin (convert x) :: Posit256) :: Posit128
  cos x = convert (cos (convert x) :: Posit256) :: Posit128
  asin x = convert (asin (convert x) :: Posit256) :: Posit128
  acos x = convert (acos (convert x) :: Posit256) :: Posit128
  atan x = convert (atan (convert x) :: Posit256) :: Posit128
  sinh x = convert (sinh (convert x) :: Posit256) :: Posit128
  cosh x = convert (cosh (convert x) :: Posit256) :: Posit128
  asinh x = convert (asinh (convert x) :: Posit256) :: Posit128
  acosh x = convert (acosh (convert x) :: Posit256) :: Posit128
  atanh x = convert (atanh (convert x) :: Posit256) :: Posit128

instance Floating Posit256 where
  pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286 :: Posit256
  exp = funExp
  log = funLogDomainReduction funLogTaylor
  (**) = funPow
  sin = funSin
  cos = funCos
  asin = funAsin
  acos = funAcos
  atan = funAtan
  sinh = funSinh
  cosh = funCosh
  asinh = funAsinh
  acosh = funAcosh
  atanh = funAtanh





class AltFloating p where
  phi :: p
  gamma :: p -> p
  sinc :: p -> p
  expm1 :: p -> p

instance AltFloating Posit8 where
  phi = convert (phi :: Posit256) :: Posit8
  gamma x = convert (gamma (convert x) :: Posit256) :: Posit8
  sinc x = convert (sinc (convert x) :: Posit256) :: Posit8
  expm1 x =
    let b = atanh $ x / 2
    in (2 * b) / (1 - b)

instance AltFloating Posit16 where
  phi = convert (phi :: Posit256) :: Posit16
  gamma x = convert (gamma (convert x) :: Posit256) :: Posit16
  sinc x = convert (sinc (convert x) :: Posit256) :: Posit16
  expm1 x =
    let b = atanh $ x / 2
    in (2 * b) / (1 - b)

instance AltFloating Posit32 where
  phi = convert (phi :: Posit256) :: Posit32
  gamma x = convert (gamma (convert x) :: Posit256) :: Posit32
  sinc x = convert (sinc (convert x) :: Posit256) :: Posit32
  expm1 x =
    let b = atanh $ x / 2
    in (2 * b) / (1 - b)

instance AltFloating Posit64 where
  phi = convert (phi :: Posit256) :: Posit64
  gamma x = convert (gamma (convert x) :: Posit256) :: Posit64
  sinc x = convert (sinc (convert x) :: Posit256) :: Posit64
  expm1 x =
    let b = atanh $ x / 2
    in (2 * b) / (1 - b)

instance AltFloating Posit128 where
  phi = convert (phi :: Posit256) :: Posit128
  gamma x = convert (gamma (convert x) :: Posit256) :: Posit128
  sinc x = convert (sinc (convert x) :: Posit256) :: Posit128
  expm1 x =
    let b = atanh $ x / 2
    in (2 * b) / (1 - b)

instance AltFloating Posit256 where
  phi = funPhi 1.6
  gamma = funGammaSeries
  sinc = funSinc
  expm1 x =
    let b = atanh $ x / 2
    in (2 * b) / (1 - b)


-- | 'phi' fixed point recursive algorithm,
funPhi :: Posit256 -> Posit256
funPhi  px@(Posit x)
    | x == x' = Posit x
    | otherwise = funPhi (Posit x')
      where
        (Posit x') = (px^2 + 2*px) / (px^2 + 1)


-- 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

-- seems pretty close to 1 ULP with the input of 0.7813
funAtan :: Posit256 -> Posit256
funAtan NaR = NaR
funAtan x
  | abs x < 1/2^122 = x  -- small angle approximaiton, found emperically
  | x < 0 = negate.funAtan $ negate x  -- if negative turn it positive, it reduces the other domain reductions by half, found from Universal CORDIC
  | x > 1 = pi/2 - funAtan (recip x)  -- if larger than one use the complementary angle, found from Universal CORDIC
  | x > twoMsqrt3 = pi/6 + funAtan ((sqrt 3 * x - 1)/(sqrt 3 + x))  -- another domain reduction, using an identity, found from https://mathonweb.com/help_ebook/html/algorithms.htm
  | otherwise = funArcTanTaylor x
--

twoMsqrt3 :: Posit256
twoMsqrt3 = 2 - sqrt 3

--
funArcTanTaylor :: Posit256 -> Posit256
funArcTanTaylor x = go 0 0
  where
    go !k !acc
      | acc == (acc + term k) = acc
      | otherwise = go (k+1) (acc + term k)
    term :: Integer -> Posit256
    term k = ((-1)^k * x^(2 * k + 1)) / fromIntegral (2 * k + 1)
--

--
funAsin :: Posit256 -> Posit256
funAsin NaR = NaR
funAsin x
  | abs x > 1 = NaR
  | x == 1 = pi/2
  | x == -1 = -pi/2
  | otherwise = funAtan w
    where
      w = x / sqrt (1 - x^2)
--

--
funAcos :: Posit256 -> Posit256
funAcos NaR = NaR
funAcos x
  | abs x > 1 = NaR
  | x < 0 = pi + funAtan invw
  | x == 0 = pi/2
  | x > 0 = funAtan invw
  | otherwise = error "Prove it covers for Rational Numbers."
    where
      invw = sqrt (1 - x^2) / x
--

-- 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


fac :: Natural -> Natural
fac 0 = 1
fac n = n * fac (n - 1)

--
funAsinh :: Posit256 -> Posit256
funAsinh NaR = NaR
funAsinh x = log $ x + sqrt (x^2 + 1)
--

--
funAcosh :: Posit256 -> Posit256
funAcosh NaR = NaR
funAcosh x
  | x < 1 = NaR
  | otherwise = log $ x + sqrt (x^2 - 1)
--

--
funAtanh :: Posit256 -> Posit256
funAtanh NaR = NaR
funAtanh x
  | abs x >= 1 = NaR
  | x < 0 = negate.funAtanh.negate $ x  -- make use of odd parity to only calculate the positive part
  | otherwise = 0.5 * log ((1+t) / (1-t)) - (fromIntegral ex / 2) * lnOf2
    where
      (ex, sig) = (int * fromIntegral (nBytes @V) + fromIntegral nat + 1, fromRational rat / 2)
      (_,int,nat,rat) = (posit2TupPosit @V).toRational $ x' -- sign should always be positive
      x' = 1 - x
      t = (2 - sig - x') / (2 + sig - x')
--

--
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)
--

--
funSin :: Posit256 -> Posit256
funSin NaR = NaR
funSin 0 = 0
funSin x = funSin' $ x / (2*pi)
--
-- funSin' is sine normalized by 2*pi
funSin' :: Posit256 -> Posit256
funSin' x
  | x == 0 = 0
  | x == 0.25 = 1
  | x == 0.5 = 0
  | x == 0.75 = -1
  | x == 1 = 0
  | x < 0 = negate.funSin'.negate $ x
  | x > 1 =
    let (_,rem) = properFraction x
    in funSin' rem
  | x > 0.75 && x < 1 = negate.funSin' $ 1 - x -- reduce domain by quadrant symmetry
  | x > 0.5 && x < 0.75 = negate.funSin' $ x - 0.5
  | x > 0.25 && x < 0.5 = funSin' $ 0.5 - x
  | x > 0.125 && x < 0.25 = funCosTuma $ 2*pi * (0.25 - x) -- reduce domain and use cofunction
  | otherwise = funSinTuma $ 2*pi * x
--

-- 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)
--

--
funSinTuma :: Posit256 -> Posit256
funSinTuma NaR = NaR
funSinTuma z = go 19 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)
--

--
funCos :: Posit256 -> Posit256
funCos NaR = NaR
funCos 0 = 1
funCos x = funCos' $ x / (2*pi)
--
-- funCos' is cosine normalized for 2*pi
funCos' :: Posit256 -> Posit256
funCos' NaR = NaR
funCos' x
  | x == 0 = 1
  | x == 0.25 = 0
  | x == 0.5 = -1
  | x == 0.75 = 0
  | x == 1 = 1
  | x < 0 = funCos'.negate $ x  -- reduce domain by symmetry across 0 to turn x positive
  | x > 1 = -- reduce domain by using perodicity
    let (_,rem) = properFraction x
    in funCos' rem
  | x > 0.75 && x < 1 = funCos' $ 1 - x  -- reduce domain by quadrant symmetry
  | x > 0.5 && x < 0.75 = negate.funCos' $ x - 0.5
  | x > 0.25 && x < 0.5 = negate.funCos' $ 0.5 - x
  | x > 0.125 && x < 0.25 = funSinTuma $ 2*pi * (0.25 - x) -- reduce domain and use cofunction
  | otherwise = funCosTuma $ 2*pi * x --
--

-- 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)
--

--
funCosTuma :: Posit256 -> Posit256
funCosTuma NaR = NaR
funCosTuma z = go 19 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)
--

-- ~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)
--


--
funLog :: Posit256 -> Posit256
funLog x = funLog2 x * lnOf2
--

--
-- Use the constant, for performance
lnOf2 :: Posit256
lnOf2 = Posit 28670435363615573179632300308403400109260626501925370561166468529302554498548
--

--
-- 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)
--

--
funLog2 :: Posit256 -> Posit256
funLog2 NaR = NaR
funLog2 z
  | z <= 0 = NaR -- includes the NaR case
  | otherwise = go (fromInteger ex) 1 sig  -- domain reduction
    where
      go :: Posit256 -> Posit256 -> Posit256 -> Posit256
      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 (nBytes @V) + fromIntegral nat, fromRational rat)
      (_,int,nat,rat) = (posit2TupPosit @V).toRational $ z -- sign should always be positive
      findSquaring m s
        | s >= 2 && s < 4 = (m, s/2)
        | otherwise = findSquaring (m+1) (s^2)
--


--  Gauss–Legendre algorithm, Seems only accurate to 2-3 ULP, but really slow
funPi1 :: Posit256
funPi1 = go 0 3 1 (recip.sqrt $ 2) (recip 4) 1
  where
    go :: Posit256 -> Posit256 -> Posit256 -> Posit256 -> Posit256 -> Posit256 -> Posit256
    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 :: Posit256
funPi2 = recip $ go 0 0 0.5 (5 / phi^3)
  where
    go :: Posit256 -> Natural -> Posit256 -> Posit256 -> Posit256
    go !prev !n !a !s
      | prev == a = a
      | 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 (n+1) (trace (show a') a') s'
--


-- Bailey–Borwein–Plouffe (BBP) formula, to 1-2 ULP, and blazing fast, converges in 60 iterations
funPi3 :: Posit256
funPi3 = go 0 0
  where
    go :: Integer -> Posit256 -> Posit256
    go !k !acc
      | acc == acc + term k = acc
      | otherwise = go (k+1) (acc + term k)
    term :: Integer -> Posit256
    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 :: Posit256
funPi4 = (1/2^6) * go 0 0
  where
    go :: Integer -> Posit256 -> Posit256
    go !k !acc
      | acc == acc + term k = acc
      | otherwise = go (k+1) (acc + term k)
    term :: Integer -> Posit256
    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)))
--



--
-- looks to be about 4 ULP accurate at -100, right on the money at -1000
funExp :: Posit256 -> Posit256
funExp x = funExp2 funExpTaylor (x / lnOf2)
--

--
--
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)



--
-- 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)
--

--
-- 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)
--

--
--
funPow :: Posit256 -> Posit256 -> Posit256
NaR `funPow` _ = NaR
_ `funPow` NaR = NaR
funPow 0 y
  | y < 0 = NaR -- NaR: Divide by Zero
  | y == 0 = NaR -- NaR: Indeterminate
  | y > 0 = 0
funPow x y
  | y < 0 = recip $ funPow x (negate y)
  | x < 0 = -- NaR if y is not an integer
    let (int,rem) = properFraction y
    in if rem == 0
       then x^^int
       else NaR -- NaR: Imaginary Number
  | otherwise = exp $ y * log x
--

-- Looks like 1 ULP for 0.7813
funSinc :: Posit256 -> Posit256
funSinc NaR = NaR
funSinc 0 = 1  -- Why the hell not!
funSinc theta = sin theta / theta
--

-- 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
-- Shannon wavelet
funPsiSha1 :: Posit256 -> Posit256
funPsiSha1 NaR = NaR
funPsiSha1 t = 2 * funSinc (2 * t) - funSinc t
--

-- Shannon wavelet
funPsiSha2 :: Posit256 -> Posit256
funPsiSha2 NaR = NaR
funPsiSha2 t = funSinc (t/2) * cos (3*pi*t/2)
--

-- Shannon wavelet
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)
--

-- funPsiSha1.(pi*) === funPsiSha3

-- 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 (nBytes @V) + fromIntegral nat + 1, fromRational rat / 2) -- move significand range from 1,2 to 0.5,1
      (_,int,nat,rat) = (posit2TupPosit @V).toRational $ x -- sign should always be positive
     
 

-- 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 :: Posit256 -> Posit256
funGammaRamanujan z = sqrt pi * (x / exp 1)**x * (8*x^3 + 4*x^2 + x + (1/30))**(1/6)
  where
    x = z - 1

--
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]

funGammaSeries :: Posit256 -> Posit256
funGammaSeries z = sqrt(2 * pi) * (z**(z - 0.5)) * exp (negate z) * (1 + series)
  where
    series :: Posit256
    series = sum $ 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."
--

funGammaCalc :: Posit256 -> Posit256
funGammaCalc z = sqrt (2*pi / z) * ((z / exp 1) * sqrt (z * sinh (recip z) + recip (810 * z^6)))**z


funGammaNemes :: Posit256 -> Posit256
funGammaNemes z = sqrt (2*pi / z) * (recip (exp 1) * (z + recip (12 * z - recip (10 * z))))**z

funGammaYang :: Posit256 -> Posit256
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 :: Posit256 -> Posit256
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