{-# LANGUAGE BangPatterns, NoImplicitPrelude #-}

module UnitFractionsDecomposition2 where

import GHC.Base
import GHC.Num ((+),(-),(*),abs,Integer)
import GHC.List (null,last,head,length,filter)
import Data.List (minimumBy)
import GHC.Real (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

{-| 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,
if it is equal to -1 - the absolute error can be only negative, otherwise - only positive
 
 - -}
setOfSolutionsG :: ErrorImpact -> Double -> [(Double, Double)]
setOfSolutionsG n k 
 | isRangeN k = 
    let j = ceiling (1.0 / k)
        p = truncate (min (2.0 / k) ((sqrt (k*k + 16) - k + 4)/(4*k))) in
            if j == p then [(fromIntegral j, let j1 = fromIntegral j in j1/(k*j1 - 1.0)) | f j ] else [(x, x/(k*x - 1.0)) | x <- [fromIntegral j..fromIntegral p], f (round x)]
 | otherwise = []
              where g y = 1.0 / y + 1.0 / fromIntegral (round (y/(k*y - 1.0))) - k
                    f y = case n of 
                            0 -> True
                            -1 -> g (fromIntegral y) < 0
                            _ -> g (fromIntegral y) > 0
{-# INLINE setOfSolutionsG #-}

setOfSolutions :: Double -> [(Double, Double)]
setOfSolutions = setOfSolutionsG 0
{-# INLINE setOfSolutions #-}

-- | Partially defined function, if there is no solutions then returns a tuple of 'undefined'. Beter
-- to use 'suitable21'
suitable2 :: Double -> (Double,Double)
suitable2 k 
 | null xs = (undefined, undefined)
 | otherwise = minimumBy (\(_, y1) (_, y2) -> compare (abs (y1 - fromIntegral (round y1))) (abs (y2 - fromIntegral (round y2)))) xs
      where !xs = setOfSolutions k
{-# 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,
if it is equal to -1 - the absolute error can be only negative, otherwise - only positive
-}
suitable21G :: ErrorImpact -> Double -> Maybe ([Double],Double)
suitable21G n k  
 | null xs = Nothing
 | otherwise = let (!a,!b) = minimumBy (\(_, y1) (_, y2) -> compare (abs (y1 - fromIntegral (round y1))) (abs (y2 - fromIntegral (round y2)))) xs
    in Just ([a,b],1.0/a + 1.0/fromIntegral (round b) - k)
      where !xs = setOfSolutionsG 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,
if it is equal to -1 - the absolute error can be only negative, otherwise - only positive
-}
check1FracDecompG :: ErrorImpact -> Double -> Maybe ([Double], Double)
check1FracDecompG n k 
 | k >= 0.005 && k <= 0.501 = 
     let c = (1.0/k)
         err = k - fromIntegral (round c)
           in if err == 0 || n == 0 || (n == -1 && err < 0) || (n /= -1 && err >= 0) then  Just ([c], fromIntegral (round c) - 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,
if it is equal to -1 - the absolute error can be only negative, otherwise - only positive
-}
check3FracDecompPartialG :: ErrorImpact -> Bool -> Double -> Maybe ([Double],Double)
check3FracDecompPartialG n direction 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] in if isNothing u then Nothing else let s2 = suitable21 (k - 1.0/fromJust u) in 
                                     if isNothing s2 then Nothing else let ([a1,b1],_) = fromJust s2 in let err = 1.0/a1 + 1.0/fromIntegral (round b1) + 1.0/fromJust u - k in
                                         if err == 0 || n == 0 || (n == -1 && err < 0) || (n /= -1 && err >= 0) then Just ([fromJust u,a1,b1],
                                              1.0/a1 + 1.0/fromIntegral (round b1) + 1.0/fromJust u - k) else Nothing
 | otherwise = Nothing
{-# 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,
if it is equal to -1 - the absolute error can be only negative, otherwise - only positive
-}
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 s2 = suitable21 (k - 1.0/fromJust u) in 
                                       if isNothing s2 then Nothing else let ([a1,b1],_) = fromJust s2 in let err = 1.0/a1 + 1.0/fromIntegral (round b1) + 1.0/fromJust u - k in
                                         if err == 0 || n == 0 || (n == -1 && err < 0) || (n /= -1 && err >= 0) then 
                                          Just ([fromJust u,a1,b1], 1.0/a1 + 1.0/fromIntegral (round b1) + 1.0/fromJust u - 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,
if it is equal to -1 - the absolute error can be only negative, otherwise - only positive
-}
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,
if it is equal to -1 - the absolute error can be only negative, otherwise - only positive
-}
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 #-}