{-# LANGUAGE BangPatterns, NoImplicitPrelude #-}

module UnitFractionsDecomposition2 where

import GHC.Base
import GHC.Num ((+),(-),(*),abs,Integer)
import GHC.List (null,last,head,length,filter,sum,cycle,zip)
import Data.List (minimumBy)
import GHC.Real (Integral,round,fromIntegral,(/),truncate,ceiling)
import GHC.Float (sqrt)
import Data.Maybe (isNothing,isJust,fromJust,catMaybes)
import Data.Ord (comparing)
import Data.Tuple (fst,snd)

-- | Rounding to thousandth.
threeDigitsK :: Double -> Double
threeDigitsK k = fromIntegral (round (k*1000)) / 1000.0
{-# INLINE threeDigitsK #-}

-- | Characterizes the impact of the absolute error sign on the approximation. 
type ErrorImpact = Int

-- | Absolute error with sign for the two unit fractions approximations and the first argument
--  (a in the related paper) being taken as the second parameter for the function. 
--  The second argument here  is expected to be 'fromIntegral' @a0@ where 'Data.List.elem' @a0@ @[2..] == @ 'True'.
absErr2Frac :: Double -> Double -> Double
absErr2Frac k y = 1.0 / y + 1.0 / fromIntegral (round (y/(k*y - 1.0))) - k
{-# INLINE absErr2Frac #-}

absErrUDecomp3 :: [Double] -> Double -> Double
absErrUDecomp3 [x,y,u] k = 1.0 / x + 1.0 / fromIntegral (round y) + 1.0 / u - k
absErrUDecomp3 [x,y] k =  1.0 / x + 1.0 / fromIntegral (round y) - k
absErrUDecomp3 [x] k = 1.0 / fromIntegral (round x) - k
absErrUDecomp3 _ _ = -1.0
{-# INLINE absErrUDecomp3 #-}

elemSolution2 :: (Integral a) => Int -> a -> Double -> Bool
elemSolution2 n y k = n == 0 || compare n 0 == compare (absErr2Frac k (fromIntegral y)) 0
{-# INLINE elemSolution2 #-}


{- | Searches for the minimum absolute error solution to two unit fractions decomposition 
 (approximation) for the fraction in the 'isRangeN' 'True' values with taking into account
the sign of the absolute error. 

If the 'ErrorImpact' parameter is equal to 0 , then absolute error can be of any sign,
otherwise the sign of it is the same as the sign of the 'ErrorImpact' argument.
-}
setOfSolutionsGmin :: ErrorImpact -> Double -> (Double, Double)
setOfSolutionsGmin n k 
 | isRangeN k = 
    let j = ceiling (1.0 / k)
        p = truncate (min (2.0 / k) ((sqrt (k*k + 16) - k + 4)/(4*k)))
        j1 = fromIntegral j in
            if j == p 
                then if elemSolution2 n j k
                         then (j1,j1/(k*j1 - 1.0)) 
                         else (0, -1.0) 
                else (\t -> if abs (t + 0.0001) < 0.000001
                                then (0, 0) 
                                else (t, t/(k*t - 1.0))) . 
                                        foldr (\x y -> if elemSolution2 n (round x) k && abs (absErr2Frac k x) < abs (absErr2Frac k y) 
                                                           then x 
                                                           else y) (-0.0001) $ [j1..fromIntegral p]
 | otherwise = (0, -1.0)
{-# INLINE setOfSolutionsGmin #-}


-- | Searches for the minimum absolute error solution to two unit fractions decomposition 
-- (approximation) for the fraction in @k@ the 'isRangeN' @k = @ 'True'.
suitable2 :: Double -> (Double,Double)
suitable2 = setOfSolutionsGmin 0
{-# INLINE suitable2 #-}

{- | Allows to take  into account the sign of the absolute error of the aproximation. If the 
'ErrorImpact' parameter is equal to 0 , then absolute error can be of any sign,
otherwise the sign of it is the same as the sign of the 'ErrorImpact' argument.
-}
suitable21G :: ErrorImpact -> Double -> Maybe ([Double],Double)
suitable21G n k = 
   if abs x < 0.000001 
       then Nothing 
       else Just ([x, y], 1.0 / x + 1.0 / fromIntegral (round y) - k) 
      where !(x,y) = setOfSolutionsGmin n k
{-# INLINE suitable21G #-}

suitable21 :: Double -> Maybe ([Double],Double)
suitable21 = suitable21G 0
{-# INLINE suitable21 #-}

isRangeN :: Double -> Bool
isRangeN k = k >= 0.005 && k <= 0.9
{-# INLINE isRangeN #-}

-- | The preferable range of the argument for 'suitable2' and 'suitable21' functions. For arguments
-- in this range the functions always have informative results.
isRangeNPref :: Double -> Bool
isRangeNPref k = k >= 0.005 && k < (2.0/3.0)
{-# INLINE isRangeNPref #-}

{- | Allows to take  into account the sign of the absolute error of the aproximation. If the 
'ErrorImpact' parameter is equal to 0 , then absolute error can be of any sign,
otherwise the sign of it is the same as the sign of the 'ErrorImpact' argument.
-}
check1FracDecompG :: ErrorImpact -> Double -> Maybe ([Double], Double)
check1FracDecompG n k 
 | k >= 0.005 && k <= 0.501 = 
     let c = (1.0/k)
         err = fromIntegral (round c) - k
           in if err * fromIntegral n >= 0  
                  then let cs = [c] in Just (cs, absErrUDecomp3 cs k) 
                  else Nothing
 | otherwise = Nothing
{-# INLINE check1FracDecompG #-}

{- | Allows to take  into account the sign of the absolute error of the aproximation. If the 
'ErrorImpact' parameter is equal to 0 , then absolute error can be of any sign,
otherwise the sign of it is the same as the sign of the 'ErrorImpact' argument.
-}
check3FracDecompPartialG :: ErrorImpact -> Bool -> Double -> Maybe ([Double],Double)
check3FracDecompPartialG n direction = check3FracDecompPartialPG n direction []
{-# INLINE check3FracDecompPartialG #-}

{- | Allows to take  into account the sign of the absolute error of the aproximation. If the 
'ErrorImpact' parameter is equal to 0 , then absolute error can be of any sign, 
otherwise the sign of it is the same as the sign of the 'ErrorImpact' argument.
-}
check3FracDecompPartialPG :: ErrorImpact -> Bool -> [Int] -> Double -> Maybe ([Double],Double)
check3FracDecompPartialPG n direction ns k 
 | k >= 0.005 && k <= 1 = 
     let u = (\us -> if null us 
                         then Nothing 
                         else Just (if direction then last us else head us)) . filter (\t -> let w = k - 1.0/t in 
                            w >= 0.005 && w <= (2.0/3.0)) $ [2.0..10.0] `mappend` map fromIntegral ns in 
                               if isNothing u 
                                   then Nothing 
                                   else let u1 = fromJust u in let s2 = suitable21 (k - 1.0 / u1) in 
                                           if isNothing s2 
                                               then Nothing 
                                               else let ([a1,b1],_) = fromJust s2 in let err = absErrUDecomp3 [a1, b1, u1] k in
                                                       if err * fromIntegral n >= 0  
                                                           then Just ([u1,a1,b1], absErrUDecomp3 [a1, b1, u1] k) 
                                                           else Nothing
 | otherwise = Nothing
{-# INLINE check3FracDecompPartialPG #-}

{- | Allows to take  into account the sign of the absolute error of the aproximation. If the 
'ErrorImpact' parameter is equal to 0 , then absolute error can be of any sign,
otherwise the sign of it is the same as the sign of the 'ErrorImpact' argument.
-}
lessErrSimpleDecompPG :: ErrorImpact -> [Int] -> Double -> (Int,Maybe ([Double],Double),Double)
lessErrSimpleDecompPG n ns k = 
  (\ts -> if null ts 
              then (0,Nothing,-1.0) 
              else let p = minimumBy (comparing (abs . snd)) ts in (length (fst p), Just p, snd p))  . 
                              catMaybes $ [check1FracDecompG n k,suitable21G n k, check3FracDecompPartialPG n True ns k, check3FracDecompPartialPG n False ns k]

{-| Allows to take  into account the sign of the absolute error of the aproximation. If the 
'ErrorImpact' parameter is equal to 0 , then absolute error can be of any sign,
otherwise the sign of it is the same as the sign of the 'ErrorImpact' argument.
-}
lessErrDenomsPG :: ErrorImpact -> [Int] -> Double -> [Integer]
lessErrDenomsPG n ns = 
  (\(_,ks,_) -> if isNothing ks 
                    then [] 
                    else map round . fst . fromJust $ ks) . lessErrSimpleDecompPG n ns
{-# INLINE lessErrDenomsPG #-}

--------------------------------------------

-- | Tries to approximate the fraction by just one unit fraction. Can be used for the numbers
-- between 0.005 and 0.501.
check1FracDecomp :: Double -> Maybe ([Double], Double)
check1FracDecomp = check1FracDecompG 0
{-# INLINE check1FracDecomp #-}

{- | Function to find the less by absolute value error approximation. One of the denominators is
taken from the range [2..10]. The two others are taken from the appropriate 'suitable21' 
applicattion.
-}
check3FracDecompPartial :: Bool -> Double -> Maybe ([Double],Double)
check3FracDecompPartial = check3FracDecompPartialG 0
{-# INLINE check3FracDecompPartial #-}

{- | Extended version of the 'check3FracDecompPartial' with the first denominator being taken not
 - only from the [2..10], but also from the custom user provided list. 
 - -}
check3FracDecompPartialP :: Bool -> [Int] -> Double -> Maybe ([Double],Double)
check3FracDecompPartialP = check3FracDecompPartialPG 0
{-# INLINE check3FracDecompPartialP #-}


lessErrSimpleDecompP :: [Int] -> Double -> (Int,Maybe ([Double],Double),Double)
lessErrSimpleDecompP = lessErrSimpleDecompPG 0

-- | Returns a list of denominators for fraction decomposition using also those ones from the the list provided as the first argument. Searches for the least error from the checked ones.
lessErrDenomsP :: [Int] -> Double -> [Integer]
lessErrDenomsP = lessErrDenomsPG 0
{-# INLINE lessErrDenomsP #-}

-- | Returns a unit fraction decomposition into 4 unit fractions. For the cases where the fraction
-- can be represented as a sum of three or less unit fractions, the last element in the resulting
-- list is really related to the floating-point arithmetic rounding errors and it depends on the
-- machine architecture and the realization of the floating-point arithmetic operations.
-- Almost any 'Double' value inside the [0.005, 2/3] can be rather well approximated by the sum of 4
-- unit fractions from the function here. There are also some intervals in the (2/3, 1] that have
-- also good approximations, though not every fraction is approximated well here.
lessErrSimpleDecomp4PG :: ErrorImpact -> [Int] -> Double -> ([Integer], Double)
lessErrSimpleDecomp4PG n ns k = 
  let !ints = lessErrDenomsPG (- 1) ns k 
      !revs = map (\t -> 1.0/ fromIntegral t) ints
      !s = sum revs
      !err = s - k
      !reverr = 1.0 / abs err
      !next 
        | n > 0 = truncate reverr
        | n < 0 = ceiling reverr
        | otherwise = round reverr
      !err4 = err + 1.0 / fromIntegral next in (ints `mappend` [next], err4)

-- | Returns a unit fraction decomposition into 5 unit fractions. For the cases where the fraction
-- can be represented as a sum of three or less unit fractions, the last element(s) in the resulting
-- list is (are) really related to the floating-point arithmetic rounding errors and it depends on the
-- machine architecture and the realization of the floating-point arithmetic operations.
-- Almost any 'Double' value inside the [0.005, 2/3] can be rather well approximated by the sum of 5
-- unit fractions from the function here. There are also some intervals in the (2/3, 1] that have
-- also good approximations, though not every fraction is approximated well here.
lessErrSimpleDecomp5PG :: ErrorImpact -> [Int] -> Double -> ([Integer], Double)
lessErrSimpleDecomp5PG n ns k = 
  let (!ints, !err4) = lessErrSimpleDecomp4PG (- 1) ns k 
      !reverr = 1.0 / abs err4
      !next 
        | n > 0 = truncate reverr
        | n < 0 = ceiling reverr
        | otherwise = round reverr
      !err5 = err4 + 1.0 / fromIntegral next in (ints `mappend` [next], err5)

----------------------------------------

-- | The argument should be greater or equal than 0.005 (1/200) though it is not checked. Returns the
-- representation of the fraction using canonical ancient Egyptian representation and its error as
-- 'Double' value in the resulting tuple.
egyptianFractionDecomposition :: Double -> ([(Integer,Integer)], Double)
egyptianFractionDecomposition k 
  | k > 2.0/3.0 && k <= 1.0 = let (ks, err) = lessErrSimpleDecomp5PG 0 [] (k - 2.0/3.0) in ((2,3) : zip (cycle [1]) ks, err)
  | otherwise = let (ks, err) = lessErrSimpleDecomp5PG 0 [] k in (zip (cycle [1]) ks, err)