|
=={{header|Haskell}}==
=== Exact rational function using Farey tree ===
In a lazy functional language Minkowski question mark function can be implemented using one of it's basic properties:
where p/q and r/s are fractions, such that |ps - rq| = 1.
This recursive definition can be implemented as lazy corecursion, i.e. by generating two infinite binary trees: '''mediant'''-based FareyStern-Brocot tree, containing all rationals, and '''mean'''-based tree with corresponding values of Minkowsky ?-function. There is one-to-one correspondence between these two trees so both {{math|?(x)}} and {{math|?<sup>-1</sup>(x)}} may be implemented as mapping between them. For details see the paper [[https://habr.com/ru/post/591949/]] (in Russian).
<lang haskell>import Data.Tree
First we define tools to handle trees.
<lang haskell>import Data.List (unfoldr)
import Data.Ratio
import Data.Tree (Tree (..), levels, unfoldTree)List
import Control.Monad.Zip (mzip)
intervalTree :: (a -> a -> a) -> (a, a) -> Tree a
intervalTree node = unfoldTree $
\(a, b) -> let m = node a b in (m, [(a,m), (m,b)])
Node a _ ==> Node b [] = const b
Node a [] ==> Node b _ = const b
Node a [l1, r1] ==> Node b [l2, r2] =
\x -> case x `compare` a of
LT -> (l1 ==> l2) x
EQ -> b
GT -> (r1 ==> r2) x
mkTreemirror :: (a ->Num a -=> a) -> a ->Tree a -> Tree a
mirror t = Node 0 [reflect (negate <$> t), t]
mkTree f a b = unfoldTree go (a, b)
where
goreflect (Node a [l,br]) = let m = fNode a b[reflect in (mr, [(a,m),reflect (m,b)l])
------------------------------------------------------------
pathBy :: Ord b => (a -> b) -> Tree a -> b -> [Either a a]
pathBy f (Node a [l,r]) x =
case x `compare` f a of
LT -> Left a : pathBy f l x
EQ -> [Right a]
GT -> Right a : pathBy f r x</lang>
sternBrocot :: Tree Rational
Now it is possible to define two trees:
sternBrocot = toRatio <$> intervalTree mediant ((0,1), (1,0))
where
mediant (p, q) (r, s) = (p + r, q + s)
<lang haskell>farey = toRatio <$> mkTree mediant (0p, 1q) (1,= 1)p % q
minkowski = toRatio <$> mkTree mean (0, 1) (1, 1)
minkowski :: Tree Rational
mediant (a,b) (c,d) = (a + c, b + d)
minkowski = toRatio <$> intervalTree mean ((0,1), (1,0))
mean (a,b) (c,d) = (a*d + c*b, 2*b*d)
toRatio (a, b) = a % b</lang>
mean (p, q) (1, 0) = (p+1, q)
mean (p, q) (r, s) = (p*s + q*r, 2*q*s)
questionMark, invQuestionMark :: Rational -> Rational
questionMark = mirror sternBrocot ==> mirror minkowski
invQuestionMark = mirror minkowski ==> mirror sternBrocot
------------------------------------------------------------
-- Floating point trees and functions
sternBrocotF :: Tree Double
sternBrocotF = mirror $ fromRational <$> sternBrocot
minkowskiF :: Tree Double
minkowskiF = mirror $ intervalTree mean (0, 1/0)
where
mean a b | isInfinite b = a + 1
| otherwise = (a + b) / 2
questionMarkF, invQuestionMarkF :: Double -> Double
questionMarkF = sternBrocotF ==> minkowskiF
invQuestionMarkF = minkowskiF ==> sternBrocotF</lang>
<pre>λ> mapM_ print $ take 4 $ levels farey
[1 % 16,3 % 16,5 % 16,7 % 16,9 % 16,11 % 16,13 % 16,15 % 16]</pre>
λ> questionMark (1/2)
Here is symmetric definitions of {{math|?(x)}} and {{math|?<sup>-1</sup>(x)}} for rational numbers:
<lang haskell>minkowskiQR :: Ratio Integer -> Ratio Integer
minkowskiQR = fromFraction . fmap transform . properFraction
where
transform = oddFunc $ lookupTree (mzip farey minkowski)
invMinkowskiQR :: Ratio Integer -> Ratio Integer
invMinkowskiQR = fromFraction . fmap transform . properFraction
where
transform = oddFunc $ lookupTree (mzip minkowski farey)
fromFraction (i, f) = fromIntegral i + f
lookupTree :: Ord a => Tree (a, c) -> a -> c
lookupTree t =
snd . either id id . last . pathBy fst t
oddFunc f 0 = 0
oddFunc f x = signum x * f (abs x)</lang>
<pre>λ> minkowskiQR (1/2)
1 % 2
λ> minkowskiQRquestionMark (2/7)
3 % 16
λ> minkowskiQRquestionMark (-22/7)
(-193) % 64
λ> invMinkowskiQRinvQuestionMark (3/16)
2 % 7
λ> invMinkowskiQRinvQuestionMark (13/256)
5 % 27</pre>
<pre>λ> questionMark $ (sqrt 5 + 1) / 2
=== Floating point function using Farey tree ===
Paths leading to numbers in Farey tree, give diadic representation of corresponding value of Minkowski ?-function and vice versa. So it is possible to use Farey tree to define Minkowski function and it's inverse for floating point numbers.
<lang haskell>minkowskiQF :: Double -> Double
minkowskiQF = oddFunc $ fromDiadic . fmap transform . properFraction
where
transform 0 = []
transform f = track (fromRational <$> farey) f
invMinkowskiQF :: Double -> Double
invMinkowskiQF = oddFunc $ fromFraction . fmap transform . toDiadic
where
transform [] = 0
transform f = follow (fromRational <$> farey) f
fromDiadic :: (Int, [Int]) -> Double
fromDiadic = fromFraction . fmap (foldr go 0 . take 55)
where
go x r = (r + fromIntegral x)/2
toDiadic :: Double -> (Int, [Int])
toDiadic = fmap (unfoldr go) . properFraction
where
go x = case properFraction (x * 2) of
(0, 0) -> Nothing
(i, f) -> Just (i `mod` 2, f)
track :: Ord a => Tree a -> a -> [Int]
track t = fmap (either (const 0) (const 1)) . pathBy id t
follow :: Tree a -> [Int] -> a
follow t lst = rootLabel $ foldl (\t -> (subForest t !!)) t $ init lst</lang>
<pre>λ> minkowskiQF (1/2)
0.5
λ> minkowskiQF (2/7)
0.1875
λ> minkowskiQF (-22/7)
-3.015625
λ> invMinkowskiQF (3/16)
0.2857142857142857
λ> invMinkowskiQF (13/256)
0.18518518518518517
λ> minkowskiQF (sqrt 2)
1.4000000000003183</pre>
The task and tests:
<lang haskell>-- sequence of all positive rationals
sternBrocot = toRatio <$> mkTree mediant (0, 1) (1, 0)
rationals = concat (levels sternBrocot)
testEq f g = all (\x -> f x == g x)
testEqF f g = all (\x -> abs (f x - g x) < 1e-11)
testIds :: [[Ratio Integer] -> Bool]
testIds =
[ testEq (invMinkowskiQR . minkowskiQR) id
, testEq (minkowskiQR . invMinkowskiQR) id . fmap minkowskiQR
, testEqF (invMinkowskiQF . minkowskiQF) id . fmap fromRational
, testEqF (minkowskiQF . invMinkowskiQF) id . fmap fromRational
, testEq (minkowskiQF . fromRational) (fromRational . minkowskiQR) ]</lang>
<pre>λ> minkowskiQF $ (sqrt 5 + 1) / 2
1.6666666666678793
λ> 5/3
1.6666666666666667
λ> invMinkowskiQFinvQuestionMark (-5/9)
-0.5657414540893351
λ> (sqrt 13 - 7)/6
-0.5657414540893352</pre>
λ> sequence testIds $ take 1000 rationals
[True,True,True,True,True]</pre>
=={{header|Julia}}==
|