euler.hs

{-# LANGUAGE Rank2Types #-}

class Numbery a where
  floating    :: a -> Float
  add         :: a -> a -> a
  mul         :: a -> a -> a

sigma :: Numbery a => [a] -> a
sigma = foldr1 add

prod :: Numbery a => [a] -> a
prod = foldr1 mul

data N = One | S N deriving Eq

instance Numbery N where
  floating One   = 1
  floating (S n) = 1 + floating n

  add One   One   = S One
  add (S n) One   = S . S $ n
  add One   (S n) = S . S $ n
  add (S n) (S m) = S . S $ add n m

  mul One   n     = n
  mul (S n) One   = S n
  mul (S n) (S m) = add (S m) $ mul n (S m)

instance Show N where
  show = show . floating

data I = Pos N | Z | Neg N deriving Eq

subN :: N -> N -> I
subN One   One   = Z
subN (S n) One   = Pos n
subN One   (S n) = Neg n
subN (S n) (S m) = subN n m

divide :: N -> N -> Maybe N
divide n m = try One where
  try x = case subN n (mul m x) of
    (Pos _) -> try (S x)
    Z       -> Just x
    (Neg _) -> Nothing

factor :: N -> [N]
factor One = []
factor x = try (S One) where
  try y = if x == y then [x] else case divide x y of
    Nothing  -> try (S y)
    (Just r) -> y : factor r

instance Numbery I where
  floating (Pos n) = floating n
  floating Z       = 0
  floating (Neg n) = -1 * floating n

  add (Pos n) (Pos m) = Pos $ add n m
  add (Pos n) Z       = Pos n
  add (Pos n) (Neg m) = subN n m
  add Z       i       = i
  add (Neg n) (Pos m) = subN m n
  add (Neg n) Z       = Neg n
  add (Neg n) (Neg m) = Neg $ add n m

  mul (Pos n) (Pos m) = Pos $ mul n m
  mul (Pos _) Z       = Z
  mul (Pos n) (Neg m) = Neg $ mul n m
  mul Z       _       = Z
  mul (Neg n) (Pos m) = Neg $ mul n m
  mul (Neg _) Z       = Z
  mul (Neg n) (Neg m) = Pos $ mul n m

instance Show I where
  show = show . floating

data Q = Div I N deriving Eq

remove :: Eq a => a -> [a] -> [a]
remove _ [] = []
remove y (x:xs) = if x == y then xs else x : remove y xs

difference :: Eq a => [a] -> [a] -> [a]
difference [] _      = []
difference l  []     = l
difference l  (y:ys) = difference (remove y l) ys

reduce :: Q -> Q
reduce (Div (Pos n) m) = Div (Pos $ r n m) $ r m n where
  r x y = let l = difference (factor x) (factor y) in if null l then One else prod l
reduce (Div Z _)       = Div Z One
reduce (Div (Neg n) m) = Div (Neg $ r n m) $ r m n where
  r x y = let l = difference (factor x) (factor y) in if null l then One else prod l

instance Numbery Q where
  floating (Div i n) = floating i / floating n

  add (Div i0 n0) (Div i1 n1) = reduce $ Div (add (mul i0 $ Pos n1) (mul i1 $ Pos n0)) (mul n0 n1)

  mul (Div i0 n0) (Div i1 n1) = reduce $ Div (mul i0 i1) (mul n0 n1)

instance Show Q where
  show = show . floating

class Profunctor p where
  dimap :: (s -> a) -> (b -> t) -> p a b -> p s t

type R = forall p. Profunctor p => p Q N -> p N N

data Forget r a b = Forget (a -> r)

instance Profunctor (Forget r) where
  dimap f _ (Forget g) = Forget $ g . f

precision :: R -> N -> Q
precision r = (\(Forget f) -> f) . r $ Forget id

upTo :: N -> [N]
upTo One   = [One]
upTo (S n) = S n : upTo n

euler :: R
euler = flip dimap id $ sigma . map (Div (Pos One) . prod . upTo) . (One :) . upTo

main :: IO ()
main = print . floating $ precision euler (S . S . S . S . S . S $ One)