{-# 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,sort)
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 = absErr2FracI 0
{-# INLINE absErr2Frac #-}

errImpact :: (Integral a) => ErrorImpact -> Double -> a
errImpact n = 
  case compare n 0 of
    GT -> truncate
    LT -> ceiling
    _ -> round
{-# INLINE errImpact #-}

-- | A generalization of the 'absErr2Frac' with the possibility to check the sign.
absErr2FracI :: ErrorImpact -> Double -> Double -> Double
absErr2FracI n k y = 1.0 / y + 1.0 / fromIntegral ((errImpact n) (y/(k*y - 1.0))) - k
{-# INLINE absErr2FracI #-}

absErrUDecomp3 :: [Double] -> Double -> Double
absErrUDecomp3 = absErrUDecomp3I 0
{-# INLINE absErrUDecomp3 #-}

absErrUDecomp3I :: ErrorImpact -> [Double] -> Double -> Double
absErrUDecomp3I n [x,y,u] k = 1.0 / x + 1.0 / fromIntegral (ceiling y) + 1.0 / fromIntegral (errImpact n u) - k
absErrUDecomp3I n [x,y] k =  1.0 / x + 1.0 / fromIntegral (errImpact n y) - k
absErrUDecomp3I n [x] k = 1.0 / fromIntegral (errImpact n x) - k
absErrUDecomp3I _ _ _ = -1.0
{-# INLINE absErrUDecomp3I #-}

elemSolution2 :: (Integral a) => Int -> a -> Double -> Bool
elemSolution2 _ _ _ = True 
{-# INLINE elemSolution2 #-}

elemSolution2I :: (Integral a) => ErrorImpact -> a -> Double -> Bool
elemSolution2I n y k = n == 0 || (fromIntegral n * absErr2FracI n k (fromIntegral y)) >= 0
{-# INLINE elemSolution2I #-}


{- | 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 j0 = 1.0 / k
        j  
          | (fromIntegral . ceiling $ j0) == j0 = ceiling (1.0 / k) + 1
          | otherwise = ceiling (1.0 / k)
        p0 = 2.0/k
        p 
          | (fromIntegral . truncate $ p0) == p0 = truncate (2.0 / k) - 1
          | otherwise = truncate (2.0/k)
        j1 = fromIntegral j in
            if j == p 
                then if elemSolution2I 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 elemSolution2I n (errImpact n x) k && abs (absErr2FracI n k x) < abs (absErr2FracI n k y) 
                                                           then x 
                                                           else y) (-0.0001) $ [min j1 (fromIntegral p)..max j1 (fromIntegral p)]
 | otherwise = (0, -1.0)
{-# INLINE setOfSolutionsGmin #-}

{- | A variant of 'setOfSolutionsGmin' where possible 'Integral' values are not equal to the second argument.
-}
setOfSolutionsGmin3 :: (Integral a)  => ErrorImpact -> a -> Double -> (Double, Double)
setOfSolutionsGmin3 n noteq k 
 | isRangeN k = 
    let j0 = 1.0 / k
        j  
          | (fromIntegral . ceiling $ j0) == j0 = ceiling (1.0 / k) + 1
          | otherwise = ceiling (1.0 / k)
        p0 = 2.0/k
        p 
          | (fromIntegral . truncate $ p0) == p0 = truncate (2.0 / k) - 1
          | otherwise = truncate (2.0/k)
        j1 = fromIntegral j in
            if j == p 
                then if elemSolution2I 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 elemSolution2I n (errImpact n x) k && abs (absErr2FracI n k x) < abs (absErr2FracI n k y) 
                                                           then x 
                                                           else y) (-0.0001) . filter (\ch -> round ch /= noteq) $ [min j1 (fromIntegral p)..max j1 (fromIntegral p)]
 | otherwise = (0, -1.0)
{-# INLINE setOfSolutionsGmin3 #-}

-- | 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 (errImpact n y) - k) 
      where !(x,y) = setOfSolutionsGmin n k
{-# INLINE suitable21G #-}

{- | A variant of 'suitable21G' where possible 'Integral' values are not equal to the second argument.
 -}
suitable21G3 :: (Integral a) => ErrorImpact -> a -> Double -> Maybe ([Double],Double)
suitable21G3 n noteq k = 
   if abs x < 0.000001 
       then Nothing 
       else Just ([x, y], 1.0 / x + 1.0 / fromIntegral (errImpact n y) - k) 
      where !(x,y) = setOfSolutionsGmin3 n noteq k
{-# INLINE suitable21G3 #-}

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)
         ec = errImpact n c
         err = 1.0/fromIntegral ec - k
         cs = [fromIntegral ec]
           in Just (cs, absErrUDecomp3I n cs k) 
 | 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 -> Double -> Maybe ([Double],Double)
check3FracDecompPartialG n = check3FracDecompPartialPG n []
{-# 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 -> [Int] -> Double -> Maybe ([Double],Double)
check3FracDecompPartialPG n ns k 
 | k >= 0.005 && k <= 1 = 
     let !s2s = catMaybes . map (suitable21G n) $ lst2check in 
            if null s2s 
                then Nothing 
                else let ([a1,b1],err3) = minimumBy (\(_,er0) (_,er1) -> compare (abs er0) (abs er1)) s2s
                         u1 = fromIntegral . round $ 1.0/(k + err3 - (1.0/a1 + 1.0/fromIntegral (errImpact n b1))) 
                            in Just (sort [u1,a1,fromIntegral (errImpact n b1)],err3) 
 | otherwise = Nothing
      where !lst2check = (\ks -> let !upp = 2.0/3.0 in foldr (\t ys -> let !w = k - 1.0/t in if w >= 0.005 && w <= upp then w:ys else ys) [] ks) $ [fromIntegral (ceiling (1.0/k))..fromIntegral (truncate (3.0/k))] `mappend` map fromIntegral ns
{-# INLINE check3FracDecompPartialPG #-}

{- | A variant of 'check3FracDecompPartialPG' where there is no possible equal 'Integral' values as
 the result.
-}
check3FracDecompPartialPG3 :: ErrorImpact -> [Int] -> Double -> Maybe ([Double],Double)
check3FracDecompPartialPG3 n ns k 
 | k >= 0.005 && k <= 1 = 
    if null s2s 
        then Nothing 
        else let ([a1,b1],err3) = minimumBy (\(_,er0) (_,er1) -> compare (abs er0) (abs er1)) s2s
                 u1 = fromIntegral . round $ 1.0/(k + err3 - (1.0/a1 + 1.0/fromIntegral (errImpact n b1))) 
                    in Just (sort [u1,a1,fromIntegral (errImpact n b1)],err3) 
 | otherwise = Nothing
      where !s2s = (\ks -> let !upp = 2.0/3.0 in foldr (\t ys -> let !w = k - 1.0/fromIntegral t in if w >= 0.005 && w <= upp then let !tt = suitable21G3 n t w in if isJust tt then fromJust tt:ys else ys else ys) [] ks) $ [ceiling (1.0/k)..truncate (3.0/k)] `mappend` ns
{-# INLINE check3FracDecompPartialPG3 #-}


{- | 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, check3FracDecompPartialPG3 n 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 (errImpact n) . 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 :: 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 :: [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 = errImpact n 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 = errImpact n reverr
      !err5 = err4 + 1.0 / fromIntegral next in (ints `mappend` [next], err5)