{- This program generates a poly-line approximation to a De Rham curve. The curve is specified using a collection of affine contractions. Under the assumption that all fixed points are in the unit circle, the program is able to guarantee an upper bound on the deviation of the actual curve from the poly-line approximation. The output is a series of lines with three columns each: the value of the parameter in [0, 1], the real part, and the imaginary part of the curve. Any real number can be specified as a fraction of two real numbers, e.g. 1/2, and any complex number can be specified in either of the forms: 1,2 1+2i 1+i2 Usage: DeRham [--q | --qout] [--cts] --tolerance TOL [--epsilon EPS] (--cesaro z | --koch-peano z | --generic [a:b:]c:d:e:f | --arrowhead a:b:c | --affine a:b:c:d:e:f | --triangle a:b:c:d:e:f) Available options: -h,--help Show this help text --q Do intermediate arithmetic in Q. This is slow but allows epsilon below 1e-14 --qout Also output the curve points as elements of Q, in a/b format --cts Don't insert empty lines in places where continuity could not be ensured --tolerance TOL Upper bound on how much the curve is allowed to deviate from the poly-line --epsilon EPS Precision with which points on the curve are calculated. Also continuity will not be ensured for parameterization steps smaller than EPS. (Default: 1e-14) --cesaro z Use Cesaro affine maps corresponding to z, a complex number, to generate Cesaro-Faber curves. Uses maps: w <- zw w <- z + (1-z)w This requires |z| < 1, |1-z| < 1. --koch-peano z Use Koch-Peano affine maps corresponding to z, a complex number, to generate Peano and Koch curves. Uses maps: w <- zw* w <- z + (1-z)w* This requires |z| < 1, |1-z| < 1. --generic [a:b:]c:d:e:f Use a generic pair of affine maps that have fixpoints at 0 and 1:  [ 1 0 0 ]  [x] <- [ 0 a d ] [x] [y] [ 0 b c ] [y]  [ 1 0 0 ]  [x] <- [ a 1-a e ] [x] [y] [ b -b f ] [y] where a, b, c, d, e, f are real numbers. This puts the midpoint of the curve at a,b. The values of a and b are 1/2,1 unless otherwise specified. --arrowhead a:b:c Use three maps that generates the arrowhead curve by mapping a triangle to three corners of itself. The complex numbers a, b, and c choose how the sides of the triangle are partitioned. Setting a=b=c=1/2 generates the standard Sierpinski triangle. Incompatible with --q. --affine a:b:c:d:e:f Specify an arbitrary affine map using real numbers a, b, c, d, e, f:  [ 1 0 0 ]  [x] <- [ e a b ] [x] [y] [ f c d ] [y] Repeat this option multiple times to specify multiple maps in this way. --triangle a:b:c:d:e:f Specify an arbitrary affine map that would map the triangle a,b,c to the triangle d,e,f; where a, b, c, d, e, f are complex numbers. Repeat this option multiple times to specify multiple maps in this way. -} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE DeriveFunctor #-} import Control.Monad import Data.List import Data.Ratio import System.IO import Options.Applicative import qualified Options.Applicative.Help as H import Text.ParserCombinators.ReadP hiding (option, many) import Numeric -- | ''Data.Complex'' requires 'RealFrac' for 'Num' data Complex a = a :+ a deriving (Eq, Show, Functor) infix 6 :+ instance Num a => Num (Complex a) where (a :+ b) + (c :+ d) = (a + c) :+ (b + d) (a :+ b) - (c :+ d) = (a - c) :+ (b - d) (a :+ b) * (c :+ d) = (a * c - b * d) :+ (a * d + b * c) negate (a :+ b) = (negate a :+ negate b) abs = error "abs" signum = error "signum" fromInteger n = (fromInteger n :+ 0) cis :: Floating a => a -> Complex a cis phi = cos phi :+ sin phi magnitudeSq :: Num a => Complex a -> a magnitudeSq (x :+ y) = x * x + y * y -- | Affine map of the plane data Affine a = Affine { dxdx :: !a , dxdy :: !a , dydx :: !a , dydy :: !a , x0 :: !a , y0 :: !a } deriving (Eq, Ord, Show, Functor) applyAffine :: Num a => Affine a -> Complex a -> Complex a applyAffine (Affine dxdp dxdq dydp dydq x0 y0) (dp :+ dq) = (x0 + dxdp * dp + dxdq * dq) :+ (y0 + dydp * dp + dydq * dq) -- | >>> affineCenter f = applyAffine f (0 :+ 0) affineCenter :: Affine a -> Complex a affineCenter (Affine _ _ _ _ x0 y0) = x0 :+ y0 -- | Frobenius norm squared affineNormSq :: Num a => Affine a -> a affineNormSq (Affine dxdx dxdy dydx dydy _ _) = dxdx^2 + dxdy^2 + dydx^2 + dydy^2 -- | >>> applyAffine (f <> g) v = applyAffine f (applyAffine g v) instance Num a => Semigroup (Affine a) where (Affine dxdp dxdq dydp dydq x0 y0) <> (Affine dpdx dpdy dqdx dqdy p0 q0) = Affine (dxdp * dpdx + dxdq * dqdx) (dxdp * dpdy + dxdq * dqdy) (dydp * dpdx + dydq * dqdx) (dydp * dpdy + dydq * dqdy) (x0 + dxdp * p0 + dxdq * q0) (y0 + dydp * p0 + dydq * q0) instance Num a => Monoid (Affine a) where mempty = Affine 1 0 0 1 0 0 -- | Assuming all fixpoints are in the unit circle around 0+0i, returns an affine map that turns the unit circle into a neighborhood of the curve whose radius is smaller than given tolerance deRhamCurveNbhd :: (Ord a, RealFrac a) => [Affine a] -> a -> a -> Affine a deRhamCurveNbhd maps !tolerance = go mempty where n = length maps go !a !t | affineNormSq a < tolerance^2 = a | i <- min (n - 1) $floor (fromIntegral n * t) = go (a <> maps !! i) (fromIntegral n * t - fromIntegral i) -- | Assuming all fixpoints are in the unit circle around 0+0i, returns an affine map that turns the unit circle into a neighborhood corresponding to the given uncertainty in the argument deRhamCurvePrecision :: (Ord a, RealFrac a) => [Affine a] -> a -> a -> Affine a deRhamCurvePrecision maps = go mempty where n = length maps go !a !dt !t | dt >= 1 = a | i <- min (n - 1)$ floor (fromIntegral n * t) = go (a <> maps !! i) (fromIntegral n * dt) (fromIntegral n * t - fromIntegral i) tabulateCurve :: (Ord a, RealFrac a) => [Affine a] -> a -> a -> (a, a) -> [Maybe (a, Complex a)] tabulateCurve maps !machineEps !tolerance (a, b) = let ea = eval a; eb = eval b in Just (a, ea) : go ea eb a b where eval !t = affineCenter $deRhamCurveNbhd maps machineEps t go !ea !eb !a !b | magnitudeSq (ea - eb) < tolerance^2 && affineNormSq precision < tolerance^2 = [Just (mid, em), Just (b, eb)] | abs (a - mid) < machineEps || abs (b - mid) < machineEps = if magnitudeSq (ea - eb) < tolerance^2 then [Just (b, eb)] else [Nothing, Just (b, eb)] | otherwise = go ea em a mid ++ go em eb mid b where dt = (b - a) / 2 mid = (a + b) / 2 precision = deRhamCurvePrecision maps dt mid em = eval mid -- | >>> applyAffine (affineFromImages (input1, output1) (input2, output2) (input3, output3)) input_i = output_i affineFromImages :: Fractional a => (Complex a, Complex a) -> (Complex a, Complex a) -> (Complex a, Complex a) -> Affine a affineFromImages (x1 :+ y1, z1 :+ w1) (x2 :+ y2, z2 :+ w2) (x3 :+ y3, z3 :+ w3) = Affine (a z1 z2 z3) (b z1 z2 z3) (a w1 w2 w3) (b w1 w2 w3) (e z1 z2 z3) (e w1 w2 w3) where det = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2) a p1 p2 p3 = (p1 * (y2 - y3) + p2 * (y3 - y1) + p3 * (y1 - y2)) / det b p1 p2 p3 = (p1 * (x3 - x2) + p2 * (x1 - x3) + p3 * (x2 - x1)) / det e p1 p2 p3 = (p1 * (x2 * y3 - x3 * y2) + p2 * (x3 * y1 - x1 * y3) + p3 * (x1 * y2 - x2 * y1)) / det cesaro :: Num a => Complex a -> [Affine a] cesaro (a :+ b) = [ Affine a (-b) b a 0 0 , Affine (1 - a) b (-b) (1 - a) a b ] kochPeano :: Num a => Complex a -> [Affine a] kochPeano (a :+ b) = [ Affine a b b (-a) 0 0 , Affine (1 - a) (-b) (-b) (a - 1) a b ] generic :: Num a => a -> a -> a -> a -> a -> a -> [Affine a] generic alpha beta delta epsilon zeta eta = [ Affine alpha delta beta epsilon 0 0 , Affine (1 - alpha) zeta (-beta) eta alpha beta ] arrowhead :: Floating a => Complex a -> Complex a -> Complex a -> [Affine a] arrowhead alpha beta gamma = [ affineFromImages (a, a) (b, q) (c, r) , affineFromImages (a, r) (b, b) (c, p) , affineFromImages (a, p) (b, q) (c, c) ] where a = 0 :+ 0 b = cis (pi / 3) c = 1 :+ 0 p = (1 - alpha) * b + alpha * c q = (1 - beta) * c + beta * a r = (1 - gamma) * a + gamma * b data OptCurves = Cesaro (Complex Rational) | KochPeano (Complex Rational) | Generic Rational Rational Rational Rational Rational Rational | Arrowhead (Complex Rational) (Complex Rational) (Complex Rational) | Maps [OptMap] deriving Show data OptMap = AffineMap Rational Rational Rational Rational Rational Rational | TriangleMap (Complex Rational, Complex Rational, Complex Rational) (Complex Rational, Complex Rational, Complex Rational) deriving Show parseCurve :: OptCurves -> Either [Affine Double] [Affine Rational] parseCurve (Cesaro z) = Right$ cesaro z parseCurve (KochPeano z) = Right $kochPeano z parseCurve (Generic a b c d e f) = Right$ generic a b c d e f parseCurve (Arrowhead a b c) = Left $arrowhead (fmap fromRational a) (fmap fromRational b) (fmap fromRational c) parseCurve (Maps maps) = Right$ map parseMap maps where parseMap (AffineMap a b c d e f) = Affine a b c d e f parseMap (TriangleMap (x1 :+ y1, x2 :+ y2, x3 :+ y3) (z1 :+ w1, z2 :+ w2, z3 :+ w3)) = affineFromImages (x1 :+ y1, z1 :+ w1) (x2 :+ y2, z2 :+ w2) (x3 :+ y3, z3 :+ w3) isContraction :: (Ord a, Fractional a) => a -> Affine a -> Bool isContraction q (Affine dxdx dxdy dydx dydy _ _) = if discr > 0 then abs (dxdx + dydy) < 2 * q && (2 * q - abs (dxdx + dydy))^2 > discr else 4 * q^2 > discr + (dxdx + dydy)^2 where discr = (dxdx - dydy)^2 + 4 * dxdy * dydx readP :: ReadP a -> ReadM a readP p = eitherReader $\str -> case readP_to_S p str of (x, ""):_ -> return x _ -> Left$ "no parse: " ++ show str fraction :: ReadP Rational fraction = liftA2 (/) (float <* skipSpaces <* char '/') (skipSpaces *> float) <|> float where float = readS_to_P readFloat signed :: Num a => ReadP a -> ReadP a signed c = (negate <$> (char '-' *> c)) <|> c point :: Num a => ReadP a -> ReadP (Complex a) point c = liftA2 (:+) (signed c <* skipSpaces <* char ',') (skipSpaces *> signed c) <|> expr where expr = liftA2 (+) (sterm <* skipSpaces <* char '+') (skipSpaces *> term) <|> liftA2 (-) (term <* skipSpaces <* char '-') (skipSpaces *> term) <|> sterm sterm = (negate <$> (char '-' *> skipSpaces *> term)) <|> term term = ((0 :+) <$> (char 'i' *> skipSpaces *> c)) <|> ((0 :+ 1) <$ char 'i') <|> ((0 :+) <$> (c <* skipSpaces <* char 'i')) <|> ((:+ 0) <$> c) (<:*:>) :: ReadP (a -> b) -> ReadP a -> ReadP b c <:*:> d = (c <* skipSpaces <* char ':') <*> (skipSpaces *> d) infixl 4 <:*:> data OptComp = DoubleComp | QComp | QOut deriving (Eq, Show) data Opts = Opts { optComp :: !OptComp , optContinuous :: !Bool , optTolerance :: !Rational , optEpsilon :: !Rational , optCurves :: OptCurves } deriving Show optsParser :: Parser Opts optsParser = pure Opts <*> (flag' QComp (long "q" <> help qHelp) <|> flag DoubleComp QOut (long "qout" <> help qOutHelp)) <*> switch (long "cts" <> help ctsHelp) <*> option fractionP (long "tolerance" <> metavar "TOL" <> help toleranceHelp) <*> (option fractionP (long "epsilon" <> metavar "EPS" <> help epsilonHelp) <|> pure 1e-14) <*> curvesParser where curvesParser = (Cesaro <$> option pointP (long "cesaro" <> metavar "z" <> helpDoc (Just cesaroHelp))) <|> (KochPeano <$> option pointP (long "koch-peano" <> metavar "z" <> helpDoc (Just kochPeanoHelp))) <|> ((\(a, b, c, d, e, f) -> Generic a b c d e f) <$> option genericP (long "generic" <> metavar "[a:b:]c:d:e:f" <> helpDoc (Just genericHelp))) <|> ((\(a, b, c) -> Arrowhead a b c) <$> option arrowheadP (long "arrowhead" <> metavar "a:b:c" <> helpDoc (Just arrowheadHelp))) <|> (Maps <$> some mapParser) fractionP = readP fraction pointP = readP$ point fraction genericP = readP $g4 <|> g6 where g4 = ((,,,,,) 0.5 1.0) <$> c <:*:> c <:*:> c <:*:> c g6 = (,,,,,) <$> c <:*:> c <:*:> c <:*:> c <:*:> c <:*:> c c = signed fraction arrowheadP = readP$ (,,) <$> c <:*:> c <:*:> c where c = point fraction mapParser = ((\(a, b, c, d, e, f) -> AffineMap a b c d e f) <$> option affineP (long "affine" <> metavar "a:b:c:d:e:f" <> helpDoc (Just affineHelp))) <|> ((\(a, b) -> TriangleMap a b) <$> option trianglesP (long "triangle" <> metavar "a:b:c:d:e:f" <> helpDoc (Just triangleHelp))) affineP = readP$ (,,,,,) <$> c <:*:> c <:*:> c <:*:> c <:*:> c <:*:> c where c = signed fraction trianglesP = readP$ (,) <$> ((,,) <$> c <:*:> c <:*:> c) <:*:> ((,,) <$> c <:*:> c <:*:> c) where c = point fraction qHelp = "Do intermediate arithmetic in Q. This is slow but allows epsilon below 1e-14" qOutHelp = "Also output the curve points as elements of Q, in a/b format" ctsHelp = "Don't insert empty lines in places where continuity could not be ensured" toleranceHelp = "Upper bound on how much the curve is allowed to deviate from the poly-line" epsilonHelp = "\ \Precision with which points on the curve are calculated. Also continuity \ \will not be ensured for parameterization steps smaller than EPS. \ \(Default: 1e-14)" cesaroHelp = H.vsep [ H.extractChunk$ H.paragraph "\ \Use Cesaro affine maps corresponding to z, a complex number, to \ \generate Cesaro-Faber curves. Uses maps:" , H.indent 4 (H.vsep [H.text "w <- zw", H.text "w <- z + (1-z)w"]) , H.extractChunk $H.paragraph "\ \This requires |z| < 1, |1-z| < 1." ] kochPeanoHelp = H.vsep [ H.extractChunk$ H.paragraph "\ \Use Koch-Peano affine maps corresponding to z, a complex number, to \ \generate Peano and Koch curves. Uses maps:" , H.indent 4 (H.vsep [H.text "w <- zw*", H.text "w <- z + (1-z)w*"]) , H.extractChunk $H.paragraph "\ \This requires |z| < 1, |1-z| < 1." ] matrix xss = H.hsep$ map (H.vsep . map H.text) ([replicate height "["] ++ xss ++ [replicate height "]"]) where height = case xss of [] -> 1 (xs:_) -> length xs genericHelp = H.vsep [ H.extractChunk $H.paragraph "\ \Use a generic pair of affine maps that have fixpoints at 0 and 1:" , H.indent 4$ H.vsep [ H.text " [ 1 0 0 ] " , H.text "[x] <- [ 0 a d ] [x]" , H.text "[y] [ 0 b c ] [y]" , H.empty , H.text " [ 1 0 0 ] " , H.text "[x] <- [ a 1-a e ] [x]" , H.text "[y] [ b -b f ] [y]" ] , H.extractChunk $H.paragraph "\ \where a, b, c, d, e, f are real numbers. This puts the midpoint of the \ \curve at a,b. The values of a and b are 1/2,1 unless otherwise \ \specified." ] arrowheadHelp = H.vsep [ H.extractChunk$ H.paragraph "\ \Use three maps that generates the arrowhead curve by mapping a triangle \ \to three corners of itself. The complex numbers a, b, and c choose how \ \the sides of the triangle are partitioned. Setting a=b=c=1/2 generates \ \the standard Sierpinski triangle. Incompatible with --q." ] affineHelp = H.vsep [ H.extractChunk $H.paragraph "\ \Specify an arbitrary affine map using real numbers a, b, c, d, e, f:" , H.indent 4$ H.vsep [ H.text " [ 1 0 0 ] " , H.text "[x] <- [ e a b ] [x]" , H.text "[y] [ f c d ] [y]" ] , H.extractChunk $H.paragraph "\ \Repeat this option multiple times to specify multiple maps in this way." ] triangleHelp = H.vsep [ H.extractChunk$ H.paragraph "\ \Specify an arbitrary affine map that would map the triangle a,b,c to \ \the triangle d,e,f; where a, b, c, d, e, f are complex numbers. Repeat \ \this option multiple times to specify multiple maps in this way." ] optsInfo :: ParserInfo Opts optsInfo = info (helper <*> optsParser) $headerDoc (Just$ H.vsep [ H.extractChunk $H.paragraph "\ \This program generates a poly-line approximation to a De Rham curve. The \ \curve is specified using a collection of affine contractions. Under the \ \assumption that all fixed points are in the unit circle, the program is \ \able to guarantee an upper bound on the deviation of the actual curve \ \from the poly-line approximation." , H.empty , H.extractChunk$ H.paragraph "\ \The output is a series of lines with three columns each: the value of the \ \parameter in [0, 1], the real part, and the imaginary part of the curve." , H.empty , H.extractChunk (H.paragraph "\ \Any real number can be specified as a fraction of two real numbers, e.g. \ \1/2, and any complex number can be specified in either of the forms:") , H.indent 4 (H.vsep [H.text "1,2", H.text "1+2i", H.text "1+i2"]) ]) main :: IO () main = do opts <- execParser optsInfo let maps = parseCurve (optCurves opts) case optComp opts of DoubleComp -> go opts id $either id (map$ fmap fromRational) maps _ -> case maps of Left maps -> do hPutStrLn stderr "Not rational maps: " hPutStrLn stderr $show maps hPutStrLn stderr "Try:" hPutStrLn stderr$ formatMaps maps Right maps -> go opts fromRational maps where go :: (RealFrac a, Show a) => Opts -> (a -> Double) -> [Affine a] -> IO () go opts fmt maps = do forM_ (filter (not . isContraction 0.999) maps) $\a -> do hPutStrLn stderr$ "Not a contraction (q > 0.999): " ++ show a forM_ (tabulateCurve maps (fromRational $optEpsilon opts) (fromRational$ optTolerance opts) (0, 1)) $\m -> case m of Just (t, x :+ y) -> do if optComp opts == QOut then putStrLn$ formatQ t ++ " " ++ formatQ x ++ " " ++ formatQ y else putStrLn $show (fmt t) ++ " " ++ show (fmt x) ++ " " ++ show (fmt y) Nothing -> unless (optContinuous opts)$ do putStrLn "" formatQ x = show (numerator q) ++ "/" ++ show (denominator q) where q = toRational x formatMaps maps = intercalate " " $map ("--affine " ++)$ map mkMap maps mkMap (Affine a b c d e f) = intercalate ":" \$ map show [a,b,c,d,e,f]