{-# 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 :: Double -> Double
threeDigitsK Double
k = Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Double
kDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
1000)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
1000.0
{-# INLINE threeDigitsK #-}

setOfSolutions :: Double -> [(Double, Double)]
setOfSolutions :: Double -> [(Double, Double)]
setOfSolutions Double
k 
 | Double -> Bool
isRangeN Double
k = 
    let j :: Integer
j = Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double
1.0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
k)
        p :: Integer
p = Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
truncate (Double -> Double -> Double
forall a. Ord a => a -> a -> a
min (Double
2.0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
k) ((Double -> Double
forall a. Floating a => a -> a
sqrt (Double
kDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
16) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
k))) in
            if Integer
j Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
p then [(Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
j, let j1 :: Double
j1 = Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
j in Double
j1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
kDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
j1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1))] else [(Double
x, Double
xDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
kDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)) | Double
x <- [Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
j..Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
p]]
 | Bool
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 :: Double -> (Double, Double)
suitable2 Double
k 
 | [(Double, Double)] -> Bool
forall a. [a] -> Bool
null [(Double, Double)]
xs = (Double
forall a. HasCallStack => a
undefined, Double
forall a. HasCallStack => a
undefined)
 | Bool
otherwise = ((Double, Double) -> (Double, Double) -> Ordering)
-> [(Double, Double)] -> (Double, Double)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy (\(Double
_, Double
y1) (Double
_, Double
y2) -> Double -> Double -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Double -> Double
forall a. Num a => a -> a
abs (Double
y1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round Double
y1))) (Double -> Double
forall a. Num a => a -> a
abs (Double
y2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round Double
y2)))) [(Double, Double)]
xs
      where !xs :: [(Double, Double)]
xs = Double -> [(Double, Double)]
setOfSolutions Double
k
{-# INLINE suitable2 #-}

suitable21 :: Double -> Maybe ([Double],Double)
suitable21 :: Double -> Maybe ([Double], Double)
suitable21 Double
k  
 | [(Double, Double)] -> Bool
forall a. [a] -> Bool
null [(Double, Double)]
xs = Maybe ([Double], Double)
forall a. Maybe a
Nothing
 | Bool
otherwise = let (!Double
a,!Double
b) = ((Double, Double) -> (Double, Double) -> Ordering)
-> [(Double, Double)] -> (Double, Double)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy (\(Double
_, Double
y1) (Double
_, Double
y2) -> Double -> Double -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Double -> Double
forall a. Num a => a -> a
abs (Double
y1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round Double
y1))) (Double -> Double
forall a. Num a => a -> a
abs (Double
y2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round Double
y2)))) [(Double, Double)]
xs
    in ([Double], Double) -> Maybe ([Double], Double)
forall a. a -> Maybe a
Just ([Double
a,Double
b],Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
1.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round Double
b)))
      where !xs :: [(Double, Double)]
xs = Double -> [(Double, Double)]
setOfSolutions Double
k
{-# INLINE suitable21 #-}

isRangeN :: Double -> Bool
isRangeN :: Double -> Bool
isRangeN Double
k = Double
k Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0.005 Bool -> Bool -> Bool
&& Double
k Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
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 :: Double -> Bool
isRangeNPref Double
k = Double
k Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0.005 Bool -> Bool -> Bool
&& Double
k Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< (Double
2.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
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 :: Double -> Maybe ([Double], Double)
check1FracDecomp Double
k 
 | Double
k Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0.005 Bool -> Bool -> Bool
&& Double
k Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0.501 = let c :: Double
c = (Double
1.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
k) in ([Double], Double) -> Maybe ([Double], Double)
forall a. a -> Maybe a
Just ([Double
c], Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round Double
c))
 | Bool
otherwise = Maybe ([Double], Double)
forall a. Maybe a
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 :: Bool -> Double -> Maybe ([Double], Double)
check3FracDecompPartial Bool
direction Double
k 
 | Double
k Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0.005 Bool -> Bool -> Bool
&& Double
k Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1 = let u :: Maybe Double
u = (\[Double]
us -> if [Double] -> Bool
forall a. [a] -> Bool
null [Double]
us then Maybe Double
forall a. Maybe a
Nothing else Double -> Maybe Double
forall a. a -> Maybe a
Just (if Bool
direction then [Double] -> Double
forall a. [a] -> a
last [Double]
us else [Double] -> Double
forall a. [a] -> a
head [Double]
us)) ([Double] -> Maybe Double)
-> ([Double] -> [Double]) -> [Double] -> Maybe Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
filter (\Double
t -> let w :: Double
w = Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
t in 
                                  Double
w Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0.005 Bool -> Bool -> Bool
&& Double
w Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= (Double
2.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
3.0)) ([Double] -> Maybe Double) -> [Double] -> Maybe Double
forall a b. (a -> b) -> a -> b
$ [Double
2.0..Double
10.0] in if Maybe Double -> Bool
forall a. Maybe a -> Bool
isNothing Maybe Double
u then Maybe ([Double], Double)
forall a. Maybe a
Nothing else let s2 :: Maybe ([Double], Double)
s2 = Double -> Maybe ([Double], Double)
suitable21 (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Maybe Double -> Double
forall a. HasCallStack => Maybe a -> a
fromJust Maybe Double
u) in 
                                     if Maybe ([Double], Double) -> Bool
forall a. Maybe a -> Bool
isNothing Maybe ([Double], Double)
s2 then Maybe ([Double], Double)
forall a. Maybe a
Nothing else let ([Double
a1,Double
b1],Double
_) = Maybe ([Double], Double) -> ([Double], Double)
forall a. HasCallStack => Maybe a -> a
fromJust Maybe ([Double], Double)
s2 in 
                                        ([Double], Double) -> Maybe ([Double], Double)
forall a. a -> Maybe a
Just ([Maybe Double -> Double
forall a. HasCallStack => Maybe a -> a
fromJust Maybe Double
u,Double
a1,Double
b1],Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round Double
b1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Maybe Double -> Double
forall a. HasCallStack => Maybe a -> a
fromJust Maybe Double
u)
 | Bool
otherwise = Maybe ([Double], Double)
forall a. Maybe a
Nothing
{-# INLINE check3FracDecompPartial #-}

lessErrSimpleDecomp :: Double -> (Int,Maybe ([Double],Double),Double)
lessErrSimpleDecomp :: Double -> (Int, Maybe ([Double], Double), Double)
lessErrSimpleDecomp Double
k = (\[([Double], Double)]
ts -> if [([Double], Double)] -> Bool
forall a. [a] -> Bool
null [([Double], Double)]
ts then (Int
0,Maybe ([Double], Double)
forall a. Maybe a
Nothing,-Double
1.0) else let p :: ([Double], Double)
p = (([Double], Double) -> ([Double], Double) -> Ordering)
-> [([Double], Double)] -> ([Double], Double)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy ((([Double], Double) -> Double)
-> ([Double], Double) -> ([Double], Double) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (Double -> Double
forall a. Num a => a -> a
abs (Double -> Double)
-> (([Double], Double) -> Double) -> ([Double], Double) -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Double], Double) -> Double
forall a b. (a, b) -> b
snd)) [([Double], Double)]
ts in ([Double] -> Int
forall a. [a] -> Int
length (([Double], Double) -> [Double]
forall a b. (a, b) -> a
fst ([Double], Double)
p), ([Double], Double) -> Maybe ([Double], Double)
forall a. a -> Maybe a
Just ([Double], Double)
p, ([Double], Double) -> Double
forall a b. (a, b) -> b
snd ([Double], Double)
p))  ([([Double], Double)] -> (Int, Maybe ([Double], Double), Double))
-> ([Maybe ([Double], Double)] -> [([Double], Double)])
-> [Maybe ([Double], Double)]
-> (Int, Maybe ([Double], Double), Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. 
                            [Maybe ([Double], Double)] -> [([Double], Double)]
forall a. [Maybe a] -> [a]
catMaybes ([Maybe ([Double], Double)]
 -> (Int, Maybe ([Double], Double), Double))
-> [Maybe ([Double], Double)]
-> (Int, Maybe ([Double], Double), Double)
forall a b. (a -> b) -> a -> b
$ [Double -> Maybe ([Double], Double)
check1FracDecomp Double
k,Double -> Maybe ([Double], Double)
suitable21 Double
k, Bool -> Double -> Maybe ([Double], Double)
check3FracDecompPartial Bool
True Double
k, Bool -> Double -> Maybe ([Double], Double)
check3FracDecompPartial Bool
False Double
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 :: Double -> [Integer]
lessErrDenoms = (\(Int
_,Maybe ([Double], Double)
ks,Double
_) -> if Maybe ([Double], Double) -> Bool
forall a. Maybe a -> Bool
isNothing Maybe ([Double], Double)
ks then [] else (Double -> Integer) -> [Double] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round ([Double] -> [Integer])
-> (Maybe ([Double], Double) -> [Double])
-> Maybe ([Double], Double)
-> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Double], Double) -> [Double]
forall a b. (a, b) -> a
fst (([Double], Double) -> [Double])
-> (Maybe ([Double], Double) -> ([Double], Double))
-> Maybe ([Double], Double)
-> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Maybe ([Double], Double) -> ([Double], Double)
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe ([Double], Double) -> [Integer])
-> Maybe ([Double], Double) -> [Integer]
forall a b. (a -> b) -> a -> b
$ Maybe ([Double], Double)
ks) ((Int, Maybe ([Double], Double), Double) -> [Integer])
-> (Double -> (Int, Maybe ([Double], Double), Double))
-> Double
-> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> (Int, Maybe ([Double], Double), Double)
lessErrSimpleDecomp
{-# INLINE lessErrDenoms #-}