mniip/DeRham.hs

## DeRham.hs
{-
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]    [ 1 0 0 ] [1]
                               [x] <- [ 0 a d ] [x]
                               [y]    [ 0 b c ] [y]

                               [1]    [ 1  0  0 ] [1]
                               [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]    [ 1 0 0 ] [1]
                               [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]    [ 1 0 0 ] [1]"
        , H.text "[x] <- [ 0 a d ] [x]"
        , H.text "[y]    [ 0 b c ] [y]"
        , H.empty
        , H.text "[1]    [ 1  0  0 ] [1]"
        , 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]    [ 1 0 0 ] [1]"
        , 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]
	{-
	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] [ 1 0 0 ] [1]
	[x] <- [ 0 a d ] [x]
	[y] [ 0 b c ] [y]

	[1] [ 1 0 0 ] [1]
	[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] [ 1 0 0 ] [1]
	[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] [ 1 0 0 ] [1]"
	, H.text "[x] <- [ 0 a d ] [x]"
	, H.text "[y] [ 0 b c ] [y]"
	, H.empty
	, H.text "[1] [ 1 0 0 ] [1]"
	, 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] [ 1 0 0 ] [1]"
	, 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]