{-# LANGUAGE BangPatterns, NoImplicitPrelude #-}

module UnitFractionsDecomposition2 where

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

threeDigitsK :: Double -> Double
threeDigitsK k = fromIntegral (round (k*1000)) / 1000.0
{-# INLINE threeDigitsK #-}

setOfSolutions :: Double -> [(Double, Double)]
setOfSolutions k0 
 | isRangeN k0 = 
    let k = threeDigitsK k0
        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 k0 
 | null xs = (undefined, undefined)
 | otherwise = minimumBy (\(_, y1) (_, y2) -> compare (abs (y1 - fromIntegral (round y1))) (abs (y2 - fromIntegral (round y2)))) xs
      where !xs = setOfSolutions k0
{-# INLINE suitable2 #-}

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

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

check1FracDecomp :: Double -> Maybe ([Double], Double)
check1FracDecomp k0 
 | k0 >= 0.005 && k0 <= 0.501 = let k = threeDigitsK k0
                                    c = (1.0/k) in Just ([c], k - fromIntegral (round c))
 | otherwise = Nothing
{-# INLINE check1FracDecomp #-}

check3FracDecompPartial :: Bool -> Double -> Maybe ([Double],Double)
check3FracDecompPartial direction k0 
 | k0 >= 0.005 && k0 < 1 = let k = threeDigitsK k0
                               u = (\us -> if null us then Nothing else Just (if direction then last us else head us)) . catMaybes . map (\t -> let w = k - 1.0/t in if w >= 0.005 && w <= (2.0/3.0) then Just t else Nothing) $ [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],err3) = fromJust s2 in
        Just ([fromJust u,a1,b1],k0 - 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]

lessErrDenoms :: Double -> [Integer]
lessErrDenoms = (\(_,ks,_) -> if isNothing ks then [] else map round . fst . fromJust $ ks) . lessErrSimpleDecomp
{-# INLINE lessErrDenoms #-}