Phase.hs

module Synthesizer.Basic.Phase
   (T,
    fromRepresentative,
    toRepresentative,
    increment,
    decrement,
    multiply,
   ) where

import qualified Algebra.RealRing              as RealRing
import qualified Algebra.Ring                  as Ring
import qualified Algebra.Additive              as Additive

import qualified Algebra.ToInteger             as ToInteger

import System.Random (Random(..))
import Test.QuickCheck (Arbitrary(arbitrary), choose)

import Foreign.Storable (Storable(..), )
import Foreign.Ptr (castPtr, )

-- import Data.Function.HT (powerAssociative, )
import Data.Tuple.HT (mapFst, )

import qualified NumericPrelude.Numeric as NP
import NumericPrelude.Numeric
import NumericPrelude.Base
import Prelude ()

import qualified GHC.Float as GHC


newtype T a = Cons {decons :: a}
   deriving Eq


instance Show a => Show (T a) where
   showsPrec p x =
      showParen (p >= 10)
         (showString "Phase.fromRepresentative " . showsPrec 11 (toRepresentative x))

instance Storable a => Storable (T a) where
   {-# INLINE sizeOf #-}
   sizeOf = sizeOf . toRepresentative
   {-# INLINE alignment #-}
   alignment = alignment . toRepresentative
   {-# INLINE peek #-}
   peek ptr = fmap Cons $ peek (castPtr ptr)
   {-# INLINE poke #-}
   poke ptr = poke (castPtr ptr) . toRepresentative


instance (Ring.C a, Random a) => Random (T a) where
   randomR = error "Phase.randomR makes no sense"
   random = mapFst Cons . randomR (zero, one)

instance (Ring.C a, Random a) => Arbitrary (T a) where
   arbitrary = fmap Cons $ choose (zero, one)



{-# INLINE fromRepresentative #-}
fromRepresentative :: RealRing.C a => a -> T a
fromRepresentative = Cons . RealRing.fraction

{-# INLINE toRepresentative #-}
toRepresentative :: T a -> a
toRepresentative = decons

{- test, how fast the function can be, if we assume that the increment is smaller than one
{-# INLINE increment #-}
increment :: RealRing.C a => a -> T a -> T a
increment d = (+ Cons d)

{-# INLINE decrement #-}
decrement :: RealRing.C a => a -> T a -> T a
decrement d = Additive.subtract (Cons d)
-}

{-# INLINE increment #-}
increment :: RealRing.C a => a -> T a -> T a
increment d = lift (d Additive.+)

{-# INLINE decrement #-}
decrement :: RealRing.C a => a -> T a -> T a
decrement d = lift (Additive.subtract d)


{-# INLINE add #-}
add :: (Ring.C a, Ord a) => T a -> T a -> T a
add (Cons x) (Cons y) =
   let z = x+y
   in  Cons $ if z>=one then z-one else z

{-# INLINE sub #-}
sub :: (Ring.C a, Ord a) => T a -> T a -> T a
sub (Cons x) (Cons y) =
   let z = x-y
   in  Cons $ if z<zero then z+one else z

{-# INLINE neg #-}
neg :: (Ring.C a, Ord a) => T a -> T a
neg (Cons x) =
   Cons $ if x==zero then zero else one-x

{-# INLINE multiply #-}
multiply :: (RealRing.C a, ToInteger.C b) => b -> T a -> T a
multiply n = lift (NP.fromIntegral n Ring.*)

{-
This implementation computes the fraction several times.
We hope that it can reduce cancellations,
but interim rounding errors seem to be equally bad.

It is certainly slower than 'multiply' and
it needs as many iterations as the number of bits of the multiplier.

> *Synthesizer.Basic.Phase> multiplyPrecise  (1000000::Integer) (fromRepresentative 2.3) :: T Double
> Phase.fromRepresentative 0.9999999998223643
> *Synthesizer.Basic.Phase> multiply  (1000000::Integer) (fromRepresentative 2.3) :: T Double
> Phase.fromRepresentative 0.999999999825377

{-# INLINE multiplyPrecise #-}
multiplyPrecise :: (RealRing.C a, ToInteger.C b) => b -> T a -> T a
multiplyPrecise n x =
   if n<zero
     then powerAssociative (+) zero (neg x) (fromIntegral (negate n))
     else powerAssociative (+) zero x (fromIntegral n)
-}


instance RealRing.C a => Additive.C (T a) where
   {-# INLINE zero #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   {-# INLINE negate #-}
   zero = Cons Additive.zero
   (+) = add
   (-) = sub
   negate = neg
{-
This implementation requires fromRepresentative,
that needs to do checks on the size of numbers
in order to choose between float2Int/int2Float and Prelude.properFraction
   (+) = lift2 (Additive.+)
   (-) = lift2 (Additive.-)
   negate = lift Additive.negate
-}

{-# INLINE lift #-}
lift :: (RealRing.C b) =>
   (a -> b) -> T a -> T b
lift f =
   fromRepresentative . f . toRepresentative

{-
{-# INLINE lift2 #-}
lift2 :: (RealRing.C c) =>
   (a -> b -> c) -> T a -> T b -> T c
lift2 f x y =
   fromRepresentative (f (toRepresentative x) (toRepresentative y))
-}


{-# INLINE customFromRepresentative #-}
customFromRepresentative ::
   (Additive.C a) =>
   (a -> i) -> (i -> a) -> a -> T a
customFromRepresentative toInt fromInt x =
   Cons (x Additive.- fromInt (toInt x))

{-# INLINE customLift #-}
customLift ::
   (Additive.C b) =>
   (b -> i) -> (i -> b) ->
   (a -> b) -> T a -> T b
customLift toInt fromInt f =
   customFromRepresentative toInt fromInt . f . toRepresentative

{-
{-# INLINE customLift2 #-}
customLift2 ::
   (Additive.C c) =>
   (c -> i) -> (i -> c) ->
   (a -> b -> c) -> T a -> T b -> T c
customLift2 toInt fromInt f x y =
   customFromRepresentative toInt fromInt $
   f (toRepresentative x) (toRepresentative y)
-}

{-# INLINE customMultiply #-}
customMultiply ::
   (Ring.C a, Ord a, ToInteger.C b) =>
   (a -> i) -> (i -> a) ->
   b -> T a -> T a
customMultiply toInt fromInt n (Cons x) =
   customFromRepresentative toInt fromInt $
   if n<zero && x>zero
     then (one-x) * NP.fromIntegral (NP.negate n)
     else x * NP.fromIntegral n


{- |
Optimization for the case,
that the integral part of the number is non-negative and fits in an Int.
This is the case for addition and integral scaling.

FIXME:
The increment and decrement routines are a bit dangerous,
because they fail if the increment value is large than maxBound::Int.
However, we will always use increments with absolute value below one.
-}
{-# RULES

     "Phase.multiply @ Float"  multiply = customMultiply GHC.float2Int  GHC.int2Float;
     "Phase.multiply @ Double" multiply = customMultiply GHC.double2Int GHC.int2Double;

     "Phase.increment @ Float"  increment = \d -> customLift GHC.float2Int  GHC.int2Float  (+d);
     "Phase.increment @ Double" increment = \d -> customLift GHC.double2Int GHC.int2Double (+d);

     "Phase.decrement @ Float"  decrement = \d -> customLift GHC.float2Int  GHC.int2Float  (subtract d);
     "Phase.decrement @ Double" decrement = \d -> customLift GHC.double2Int GHC.int2Double (subtract d);

  #-}

{-
     "Phase.+        @ Float"  (+) = customLift2 GHC.float2Int GHC.int2Float (+);
     "Phase.-        @ Float"  (-) = customLift2 GHC.float2Int GHC.int2Float (-);

     "Phase.+        @ Double" (+) = customLift2 GHC.double2Int GHC.int2Double (+);
     "Phase.-        @ Double" (-) = customLift2 GHC.double2Int GHC.int2Double (-);

-}