{-# 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 #-}

setOfSolutions :: Double -> [(Double, Double)]
setOfSolutions 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))] else [(x, x/(k*x - 1)) | x <- [fromIntegral j..fromIntegral p]]
 | otherwise = []
{-# 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 #-}

suitable21 :: Double -> Maybe ([Double],Double)
suitable21 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],k - (1.0/a + 1.0/fromIntegral (round b)))
      where !xs = setOfSolutions k
{-# 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 #-}

-- | 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 k 
 | k >= 0.005 && k <= 0.501 = let c = (1.0/k) in Just ([c], k - fromIntegral (round c))
 | otherwise = Nothing
{-# 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 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 
                                        Just ([fromJust u,a1,b1],k - 1.0/a1 - 1.0/fromIntegral (round b1) - 1.0/fromJust u)
 | otherwise = Nothing
{-# INLINE check3FracDecompPartial #-}

lessErrSimpleDecomp :: Double -> (Int,Maybe ([Double],Double),Double)
lessErrSimpleDecomp 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 $ [check1FracDecomp k,suitable21 k, check3FracDecompPartial True k, check3FracDecompPartial False k]

-- | A list of denominators for fraction decomposition. If the list has 3 elements (likely in most
-- cases) then the furst one is from the range [2..10].
lessErrDenoms :: Double -> [Integer]
lessErrDenoms = (\(_,ks,_) -> if isNothing ks then [] else map round . fst . fromJust $ ks) . lessErrSimpleDecomp
{-# INLINE lessErrDenoms #-}