brush-strokes-0.1.0.0: src/lib/Math/Bezier/ArcLength.hs

{-# LANGUAGE AllowAmbiguousTypes  #-}
{-# LANGUAGE BangPatterns         #-}
{-# LANGUAGE DuplicateRecordFields #-}
{-# LANGUAGE ScopedTypeVariables  #-}

module Math.Bezier.ArcLength
  ( -- * Arc-length reparametrisation
    ArcLengthParametrisation(..)
  , curveArcLengthParametrisation

  -- ** Options
  , ArcLengthOptions(..), defaultArcLengthOptions
  , Integrator(..)
  , Regulariser(..), defaultRegulariser
  , RegularisationMethod(..)
  , NewtonRaphsonOptions(..), defaultNewtonRaphsonOptions

  -- * Spline arc-length parametrisation
  , SplineArcLengthParam(..)
  , splineArcLengthParam
  , splineTotalArcLength
  , splineParameterAtArcLength
  , splitSplineAtArcLengths
  )
 where

-- base
import Data.Complex
  ( realPart )
import Data.Maybe
  ( fromMaybe )
import GHC.Exts
  ( fromList, proxy# )
import GHC.TypeNats
  ( KnownNat, natVal', SomeNat (..), someNatVal )

-- acts
import Data.Act
  ( Torsor (..) )

-- containers
import Data.Sequence
  ( Seq )
import qualified Data.Sequence as Seq
  ( Seq(..) )

-- primitive
import Data.Primitive.PrimArray
  ( PrimArray )
import qualified Data.Primitive.PrimArray as Prim
  ( indexPrimArray )

-- brush-strokes
import qualified Math.Bezier.Cubic as Cubic
  ( Bezier(..), bezier', speedCriticalPoints, subdivide )
import qualified Math.Bezier.Quadratic as Quadratic
  ( Bezier(..), bezier', speedCriticalPoints, subdivide )
import Math.Bezier.Spline
  ( Curve(..), Curves(..), NextPoint(..), Spline(..), SplineType(..)
  , SplineTypeI
  , adjustSplineType, openCurveEnd, subdivideSplineAt
  )
import Math.Epsilon
  ( epsilon )
import Math.Module
  ( Module, Inner, norm )
import Math.Roots
  ( newtonRaphson )

--------------------------------------------------------------------------------
-- Curve arc-length

data ArcLengthOptions
  = ArcLengthOptions
  { arcLengthIntegrator :: !Integrator
      -- ^ Quadrature method used for arc-length integration.
  , splitAtCriticalPoints :: !Bool
      -- ^ Whether to split at (near) critical points.
      -- (A "near" critical point is the real part of a possibly complex critical point.)
      --
      -- (Ignored for the Gravesen method.)
  , newtonRaphsonOptions :: !NewtonRaphsonOptions
      -- ^ Options for the Newton–Raphson method used to compute the
      -- inverse of the arc-length function.
  }

-- | Method used to regularise the integrand near a low-speed endpoint.
data RegularisationMethod
  = SquareRoot
    -- ^ Substitute \( t = u^2 \) near the endpoint, so that the integrand
    --   \( \|B'(t)\| \sim \sqrt{t} \) becomes smooth in \( u \).
    --   (No other methods yet; more may be added in the future.)

-- | Cusp-endpoint regularisation settings.
data Regulariser
  = Regulariser
  { regularisationMethod    :: !RegularisationMethod
  , regularisationThreshold :: !Double
      -- ^ Speed threshold below which the substitution is applied at an endpoint.
  }

defaultRegulariser :: Regulariser
defaultRegulariser = Regulariser
  { regularisationMethod    = SquareRoot
  , regularisationThreshold = 0.5
  }

data Integrator
  = ClenshawCurtis
    { clenshawCurtisDegree :: !Int -- ^ degree (must be even)
    , regulariser          :: !( Maybe Regulariser )
        -- ^ Regularisation at cuspidal endpoints?
    }
  | TanhSinh
    { tanhSinhNodes    :: !Int
        -- ^ @N@: nodes per side (2N+1 nodes total).
        --
        -- Recommended value: such that @N*h ≥ 3.2@.
    , tanhSinhStepSize :: !Double
        -- ^ @h@: step size.
        --
        -- Recommended value: between @0.5@ (fast) and @0.1@ (precise).
    }
  | GaussLegendre
    { gaussLegendreDegree :: !Int
        -- ^ Number of quadrature nodes.
        --
        -- An @n@-point rule is exact for polynomials of degree up to @2n−1@.
    , regulariser         :: !( Maybe Regulariser )
        -- ^ Cusp-endpoint regularisation (@Nothing@ disables it).
    }
  | Gravesen
    { gravesenTol :: !Double
        -- ^ Stopping tolerance on the polygon-length minus chord-length gap
        -- for each piece.
        -- The tolerance is halved at each subdivision level, so the total
        -- arc-length error is bounded by @gravesenTol@.
    }

defaultArcLengthOptions :: ArcLengthOptions
defaultArcLengthOptions =
  ArcLengthOptions
    { arcLengthIntegrator   = ClenshawCurtis
        { clenshawCurtisDegree = 16
        , regulariser          = Just defaultRegulariser
        }
    , splitAtCriticalPoints = True
    , newtonRaphsonOptions  = defaultNewtonRaphsonOptions
    }

data NewtonRaphsonOptions
  = NewtonRaphsonOptions
  { newtonRaphsonDigits :: !Int  -- ^ desired bits of precision
  , newtonRaphsonMaxIterations :: !Word -- ^ maximum number of iterations (at least 1)
  }

defaultNewtonRaphsonOptions :: NewtonRaphsonOptions
defaultNewtonRaphsonOptions =
  NewtonRaphsonOptions
    { newtonRaphsonDigits = 50
    , newtonRaphsonMaxIterations = 64
    }

data ArcLengthParametrisation
  = ArcLengthParametrisation
  { totalArcLength :: !Double
  , parameterFromArcLength :: !( Double -> Double )
  }

-- | Compute the arc-length parametrisation of a 'Curve'.
curveArcLengthParametrisation
  :: forall diff ptData crvData
  .  ( Torsor diff ptData, Module Double diff, Inner Double diff
     , Show ptData )
  => ArcLengthOptions
  -> ptData
  -> Curve Open crvData ptData
  -> ArcLengthParametrisation
curveArcLengthParametrisation opts p0 = \case

  LineTo { curveEnd = NextPoint p1 } ->
    let !len = norm @diff ( p0 --> p1 )
    in ArcLengthParametrisation
        { totalArcLength = len
        , parameterFromArcLength =
            \ s ->
              if len > 0
              then max 0 $ min 1 $ s / len
              else 0
        }

  Bezier2To { controlPoint = p1, curveEnd = NextPoint p2 } ->
    let !bez    = Quadratic.Bezier { p0, p1, p2 }
        speed   = norm @diff . Quadratic.bezier' @diff bez
        !critPts =
          if splitAtCriticalPoints opts
          then Quadratic.speedCriticalPoints @diff bez
          else []
    in build speed $
         pickIntegrator ( gravesenQuadIntegrator @diff bez ) speed critPts

  Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 } ->
    let !bez    = Cubic.Bezier { p0, p1, p2, p3 }
        speed   = norm @diff . Cubic.bezier' @diff bez
        !critPts =
          if splitAtCriticalPoints opts
          then map realPart $ Cubic.speedCriticalPoints @diff bez
          else []
    in build speed $
         pickIntegrator ( gravesenCubicIntegrator @diff bez ) speed critPts

  where
    pickIntegrator
      :: ( Double -> Double -> Double ) -- ^ Gravesen cumulative arc-length integrator
      -> ( Double -> Double )   -- ^ speed function
      -> [ Double ]             -- ^ critical points
      -> ( Double -> Double )   -- ^ cumulative arc-length function
    pickIntegrator gravesenFn speed critPts =
      case arcLengthIntegrator opts of
        Gravesen { gravesenTol = eps }
          -> gravesenFn eps
        ClenshawCurtis { clenshawCurtisDegree = d, regulariser = reg }
          | SomeNat @n _ <- someNatVal ( fromIntegral d )
          -> makeCompositeCCIntegrator @n reg speed critPts
        TanhSinh { tanhSinhNodes = nodesPerSide, tanhSinhStepSize = tsH }
          -> makeCompositeTSIntegrator nodesPerSide tsH speed critPts
        GaussLegendre { gaussLegendreDegree = d, regulariser = reg }
          -> makeCompositeGLIntegrator d reg speed critPts

    -- | Build the final parametrisation from speed and cumulative arc-length.
    build speed integrate =
      let !lg = integrate 1
      in
        ArcLengthParametrisation
          { totalArcLength = lg
          , parameterFromArcLength =
              \ s ->
                if | s <= 0
                   -> 0
                   | s >= lg
                   -> 1
                   | otherwise
                   -> let t_0 = s / lg
                      in case newtonRaphson
                             ( newtonRaphsonMaxIterations $ newtonRaphsonOptions opts )
                             ( newtonRaphsonDigits $ newtonRaphsonOptions opts )
                             ( \ t -> (# integrate t - s, speed t #) )
                             t_0
                      of
                        Just t' -> t'
                        Nothing -> t_0
          }
{-# INLINEABLE curveArcLengthParametrisation #-}

--------------------------------------------------------------------------------
-- Spline arc-length

-- | Precomputed arc-length data for all curves of a spline.
data SplineArcLengthParam = SplineArcLengthParam
  { splineArcLengthCurveEnds   :: !( Seq Double )
      -- ^ Cumulative arc lengths at the end of each curve.
  , splineArcLengthCurveParams :: !( Seq ArcLengthParametrisation )
      -- ^ Per-curve arc-length parametrisations.
  }

-- | Compute the arc-length parametrisation of a spline.
splineArcLengthParam
  :: forall diff clo crvData ptData
  .  ( SplineTypeI clo
     , Torsor diff ptData, Module Double diff, Inner Double diff
     , Show ptData )
  => ArcLengthOptions
  -> Spline clo crvData ptData
  -> SplineArcLengthParam
splineArcLengthParam opts spl =
  let Spline { splineStart, splineCurves = OpenCurves curves } =
        adjustSplineType @Open spl
      pairs = go Seq.Empty 0 splineStart curves
  in SplineArcLengthParam
       { splineArcLengthCurveEnds   = fmap fst pairs
       , splineArcLengthCurveParams = fmap snd pairs
       }
  where
    go :: Seq (Double, ArcLengthParametrisation)
       -> Double
       -> ptData
       -> Seq ( Curve Open crvData ptData )
       -> Seq (Double, ArcLengthParametrisation)
    go acc _   _  Seq.Empty            = acc
    go acc cum pt ( crv Seq.:<| rest ) =
      let param = curveArcLengthParametrisation @diff opts pt crv
          cum'  = cum + totalArcLength param
      in  go ( acc Seq.:|> ( cum', param ) ) cum' ( openCurveEnd crv ) rest

-- | Total arc length of a spline.
splineTotalArcLength :: SplineArcLengthParam -> Double
splineTotalArcLength ( SplineArcLengthParam { splineArcLengthCurveEnds = ends } ) =
  case ends of
    Seq.Empty   -> 0
    _ Seq.:|> x -> x

-- | Find the @(curve index, Bézier parameter t ∈ [0,1])@ within the spline
-- at the given arc length.
splineParameterAtArcLength
  :: SplineArcLengthParam
  -> Double          -- ^ target arc length
  -> ( Int, Double ) -- ^ (curve index, Bézier parameter t ∈ [0,1])
splineParameterAtArcLength
  ( SplineArcLengthParam
    { splineArcLengthCurveEnds = ends
    , splineArcLengthCurveParams = crvParams
    }
  )
  s = go 0 0 ends crvParams
  where
    go i _cumStart Seq.Empty _
      = ( max 0 ( i - 1 ), 1 )
    go i  cumStart ( end Seq.:<| restEnds ) ( param Seq.:<| restParams )
      | s <= end || null restEnds
      = ( i, parameterFromArcLength param ( s - cumStart ) )
      | otherwise
      = go ( i + 1 ) end restEnds restParams
    go i _cumStart _ _
      = ( max 0 ( i - 1 ), 1 )

-- | Split an open spline at a sorted list of absolute arc-lengths.
splitSplineAtArcLengths
  :: forall diff crvData ptData
  .  ( Torsor diff ptData, Module Double diff, Inner Double diff
     , Show ptData -- debugging
     )
  => ArcLengthOptions
  -> [ Double ]   -- ^ sorted arc-length positions
  -> Spline Open crvData ptData
  -> [ Spline Open crvData ptData ]
splitSplineAtArcLengths opts positions spl = subdivideSplineAt @diff spl splitParams
  where
    !param      = splineArcLengthParam @diff opts spl
    splitParams = map ( splineParameterAtArcLength param ) positions

--------------------------------------------------------------------------------
-- Clenshaw–Curtis quadrature

-- | Chebyshev node of the second kind.
chebyshevNode
  :: forall n
  .  KnownNat n
  => Int -- ^ 0 <= k <= n
  -> Double
chebyshevNode k = cos ( pi * fromIntegral k / fromIntegral n )
  where
    n :: Int
    n = fromIntegral $ natVal' @n proxy#
{-# INLINEABLE chebyshevNode #-}

-- | Clenshaw–Curtis nodes.
ccNodes :: forall n. KnownNat n => PrimArray Double
ccNodes =
  fromList
    [ 0.5 * ( 1 + chebyshevNode @n k )
    | k <- [ 0 .. n ]
    ]
  where
    n :: Int
    n = fromIntegral $ natVal' @n proxy#
{-# INLINEABLE ccNodes #-}

-- | Pre computation of a Clenshaw–Curtis integrator.
makeCCIntegrator
  :: forall n
  .  KnownNat n
  => PrimArray Double
  -> ( Double -> Double )
makeCCIntegrator ys =
  -- Pre-computation
  let
    n :: Int
    n  = fromIntegral $ natVal' @n proxy#
    ep :: Int -> Double
    ep k = if k == 0 || k == n then 0.5 else 1

    -- Chebyshev coefficients (pre-computation)
    cs :: PrimArray Double
    cs =
      fromList
        [ ( 2 / fromIntegral n )
          * sum [ ep k
                * ( ys `Prim.indexPrimArray` k )
                * cos ( fromIntegral j * fromIntegral k * pi / fromIntegral n )
                | k <- [0..n]
                ]
        | j <- [0..n]
        ]

  -- The integrator proper
  in \t ->
  let y0 = 2*t - 1    -- map t ∈ [0,1] to y₀ ∈ [−1,1]

      -- Chebyshev recurrence (T_{j−1}, T_j) as plain
      go :: Double -> Double -> Double -> Int -> Double
      go !acc !tm1 !tj !j
        | j > n
        = acc
        | otherwise
        = let
            !tj1 = 2 * y0 * tj - tm1
            !cj   = cs `Prim.indexPrimArray` j
            !ej   = cj * ( if j == 0 || j == n then 0.5 else 1 )
            a | j == 0    = ej * ( y0 + 1)
              | j == 1    = ej * ( tj1 - 1 ) / 4
              | otherwise = ej * ( tj1 / fromIntegral ( 2 * ( j + 1 ))
                                 - tm1 / fromIntegral ( 2 * ( j - 1 ))
                                 - ( if even j then 1 else -1 ) / fromIntegral ( j * j - 1 ) )
          in go ( acc + a ) tj tj1 ( j + 1 )
  in 0.5 * go 0 y0 1 0
{-# INLINEABLE makeCCIntegrator #-}

--------------------------------------------------------------------------------
-- Endpoint regularisation

-- | Select the regularisation substitution for a single piece @[a,b]@.
-- Returns the transformed integrand on @[0,1]@ and the map @t ↦ u@.
applyRegulariser
  :: Maybe Regulariser
  -> Double               -- ^ @a@
  -> Double               -- ^ @b@
  -> ( Double -> Double ) -- ^ speed
  -> ( Double -> Double   -- ^ transformed integrand on @[0,1]@
     , Double -> Double ) -- ^ \( t \in [a,b] \mapsto u \in [0,1] \)
applyRegulariser Nothing a b speed =
  ( \u -> speed ( a + (b-a)*u ) * (b-a)
  , \x -> (x-a) / (b-a) )
applyRegulariser
  ( Just Regulariser { regularisationMethod = meth, regularisationThreshold = tol } )
  a b speed
  | sa <= tol && sa <= sb = case meth of
      SquareRoot ->
        ( \u -> speed ( a + (b-a)*u*u ) * 2*(b-a)*u
        , \x -> sqrt ( max 0 ( (x-a) / (b-a) ) ) )
  | sb <= tol && sb < sa  = case meth of
      SquareRoot ->
        ( \u -> speed ( b - (b-a)*(1-u)*(1-u) ) * 2*(b-a)*(1-u)
        , \x -> 1 - sqrt ( max 0 ( (b-x) / (b-a) ) ) )
  | otherwise =
      ( \u -> speed ( a + (b-a)*u ) * (b-a)
      , \x -> (x-a) / (b-a) )
  where
    sa = speed a
    sb = speed b
{-# INLINEABLE applyRegulariser #-}

--------------------------------------------------------------------------------
-- Clenshaw–Curtis quadrature (continued)

-- | Composite Clenshaw–Curtis integrator that splits @[0,1]@ at the critical
-- points of @|B′(t)|²@ and optionally applies an endpoint regularisation
-- substitution on each piece at the low-speed endpoint.
--
-- The splitting avoids bad convergence around critical points; the
-- regularisation restores exponential convergence when a boundary point is a
-- critical point (Bernstein ellipse expansion).
makeCompositeCCIntegrator
  :: forall n
  .  KnownNat n
  => Maybe Regulariser      -- ^ cusp-endpoint regularisation
  -> ( Double -> Double )   -- ^ speed function on @[0,1]@
  -> [ Double ]             -- ^ critical points of @|speed|²@ in @(0,1)@, sorted ascending
  -> ( Double -> Double )   -- ^ \( t \mapsto \int_0^t speed(x) dx \)
makeCompositeCCIntegrator reg speed critPts =
  let
    n = fromIntegral $ natVal' @n proxy#
    bps = 0 : critPts ++ [1]
    pieceData =
      [ ( a, b, integr, toU, integr 1 )
      | ( a, b ) <- zip bps ( drop 1 bps )
      , b - a > 1e-8
      , let ( f01, toU ) = applyRegulariser reg a b speed
            ys = fromList
                   [ f01 ( ccNodes @n `Prim.indexPrimArray` k )
                   | k <- [ 0 .. n ] ]
            integr = makeCCIntegrator @n ys
      ]
    !totalLen = sum [ fullLen | (_, _, _, _, fullLen) <- pieceData ]

  in \ t ->
    if | t <= 0
       -> 0
       | t >= 1
       -> totalLen
       | otherwise
       -> sum [ if b <= t then fullLen else integr ( toU t )
              | ( a, b, integr, toU, fullLen ) <- pieceData
              , a < t
              ]
{-# INLINEABLE makeCompositeCCIntegrator #-}

--------------------------------------------------------------------------------
-- Tanh–Sinh quadrature

-- | Composite Tanh–Sinh integrator that splits @[0,1]@ at the given subdivision
-- points (e.g. cusps), and applies Tanh–Sinh (double-exponential) quadrature
-- on each piece.
makeCompositeTSIntegrator
  :: Int                  -- ^ @N@: nodes per side (@2N+1@ nodes total)
  -> Double               -- ^ @h@: step size
  -> ( Double -> Double ) -- ^ speed function on @[0,1]@
  -> [ Double ]           -- ^ subdivision points
  -> ( Double -> Double ) -- ^ \( t \mapsto \int_0^t speed(x) dx \)
makeCompositeTSIntegrator nodesPerSide h speed subdivPts =
  let
    pi_2 :: Double
    pi_2 = pi / 2

    bps :: [ Double ]
    bps = 0 : subdivPts ++ [1]

    n :: Int
    n = 2 * nodesPerSide + 1  -- total number of nodes

    -- Pre-compute normalised Tanh–Sinh nodes and weights for [-1,1]:
    --   φ_k = tanh(π/2 · sinh(k·h))                      node, ∈ (−1,1)
    --   ψ_k = π/2 · cosh(k·h) / cosh²(π/2 · sinh(k·h))   weight kernel
    --
    -- Integral on [a,b] ≈ h · (b−a)/2 · Σ_k ψ_k · f( (a+b)/2 + (b−a)/2 · φ_k )
    φs :: PrimArray Double
    !φs = fromList
      [ tanh ( pi_2 * sinh ( fromIntegral k * h ) )
      | k <- [ -nodesPerSide .. nodesPerSide ]
      ]

    ψs :: PrimArray Double
    !ψs = fromList
      [ let skh = sinh ( fromIntegral k * h )
            ckh = cosh ( fromIntegral k * h )
            c   = cosh ( pi_2 * skh )
        in  pi_2 * ckh / ( c * c )
      | k <- [ -nodesPerSide .. nodesPerSide ]
      ]

    -- Tanh–Sinh quadrature of 'speed' on [a, b].
    tsIntegral :: Double -> Double -> Double
    tsIntegral a b =
      let mid  = ( a + b ) * 0.5
          half = ( b - a ) * 0.5
          go :: Double -> Int -> Double
          go !acc !i
            | i >= n    = acc
            | otherwise =
              let φ = φs `Prim.indexPrimArray` i
                  ψ = ψs `Prim.indexPrimArray` i
              in go ( acc + ψ * speed ( mid + half * φ ) ) ( i + 1 )
      in h * half * go 0 0

    pieceData :: [ ( Double, Double, Double ) ]
    pieceData =
      [ ( a, b, tsIntegral a b )
      | ( a, b ) <- zip bps ( drop 1 bps )
      , b - a > 1e-8
      ]

    totalLen :: Double
    !totalLen = sum [ fullLen | ( _, _, fullLen ) <- pieceData ]

  in \ t ->
    if
      | t <= 0
      -> 0
      | t >= 1
      -> totalLen
      | otherwise
      -> sum [ if t >= b then fullLen else tsIntegral a t
             | ( a, b, fullLen ) <- pieceData
             , a < t
             ]
{-# INLINEABLE makeCompositeTSIntegrator #-}

--------------------------------------------------------------------------------
-- Gauss–Legendre quadrature

-- | Compute @n@-point Gauss–Legendre nodes (on [−1,1]) and weights
-- using Newton iteration on the Legendre polynomial recurrence.
computeGLNodesWeights :: Int -> ( PrimArray Double, PrimArray Double )
computeGLNodesWeights n =
  let
    -- Evaluate (P_n(x), P_n′(x)) via the three-term recurrence:
    --   P_0 = 1,  P_0′ = 0
    --   P_1 = x,  P_1′ = 1
    --   (k+1) P_{k+1} = (2k+1) x P_k − k P_{k−1}
    legendre :: Double -> (# Double, Double #)
    legendre x
      | n == 0    = (# 1, 0 #)
      | otherwise = go 1 1 x 0 1
      where
        go :: Int -> Double -> Double -> Double -> Double -> (# Double, Double #)
        go !k !p_0 !p_1 !dp_0 !dp_1
          | k == n    = (# p_1, dp_1 #)
          | otherwise =
              let
                kd, p_2, dp_2 :: Double
                !kd  = fromIntegral k
                !p_2  = ( ( 2 * kd + 1 ) * x * p_1 - kd * p_0 ) / ( kd + 1 )
                !dp_2 = ( ( 2 * kd + 1 ) * ( p_1 + x * dp_1 ) - kd * dp_0 ) / ( kd + 1 )
              in go (k+1) p_1 p_2 dp_1 dp_2

    -- Newton iteration to refine a root of P_n.
    refineRoot :: Double -> Double
    refineRoot x0 = fromMaybe x0 $ newtonRaphson 50 50 legendre x0

    -- Compute the m largest positive roots (including 0 for odd n).
    m = ( n + 1 ) `div` 2

    halfData :: [ ( Double, Double ) ]
    halfData =
      [ let x_i  = refineRoot ( cos ( pi * ( fromIntegral k - 0.25 )
                                         / ( fromIntegral n + 0.5  ) ) )
            (# _, dp_i #) = legendre x_i
            w_i  = 2 / ( ( 1 - x_i * x_i ) * dp_i * dp_i )
        in ( x_i, w_i )
      | k <- [ 1 .. m ]
      ]

    -- Use symmetry to obtain all nodes
    ( xs, ws ) = unzip $
      concatMap
        ( \ ( x_i, w_i ) ->
            if abs x_i < epsilon
            then [ (  0  , w_i ) ]
            else [ (  x_i, w_i )
                 , ( -x_i, w_i ) ]
        ) halfData

  in ( fromList xs
     , fromList ws )
{-# INLINEABLE computeGLNodesWeights #-}

-- | Composite Gauss–Legendre integrator that splits @[0,1]@ at the critical
-- points of @|B′(t)|²@ and optionally applies an endpoint regularisation
-- substitution at the low-speed endpoint of each piece.
makeCompositeGLIntegrator
  :: Int                    -- ^ number of quadrature nodes
  -> Maybe Regulariser      -- ^ cusp-endpoint regularisation
  -> ( Double -> Double )   -- ^ speed function on @[0,1]@
  -> [ Double ]             -- ^ critical points of @|speed|²@ in @(0,1)@, sorted ascending
  -> ( Double -> Double )   -- ^ \( t \mapsto \int_0^t speed(x) dx \)
makeCompositeGLIntegrator nPts reg speed critPts =
  let
    bps = 0 : critPts ++ [ 1 ]

    ( glXs, glWs ) = computeGLNodesWeights nPts

    -- GL quadrature of f on [a, b].
    glQuad :: ( Double -> Double ) -> Double -> Double -> Double
    glQuad f a b =
      let mid  = ( a + b ) * 0.5
          half = ( b - a ) * 0.5
          go :: Double -> Int -> Double
          go !acc !i
            | i >= nPts = acc
            | otherwise =
                let xi = glXs `Prim.indexPrimArray` i
                    wi = glWs `Prim.indexPrimArray` i
                in go ( acc + wi * f ( mid + half * xi ) ) ( i + 1 )
      in half * go 0 0

    pieceData :: [ ( Double, Double, Double -> Double -> Double, Double -> Double, Double ) ]
    pieceData =
      [ ( a, b, glQuad f01, toU, fullLen )
      | ( a, b ) <- zip bps ( drop 1 bps )
      , b - a > 1e-8
      , let ( f01, toU ) = applyRegulariser reg a b speed
      , let fullLen = glQuad f01 0 1
      ]

    !totalLen = sum [ fullLen | ( _, _, _, _, fullLen ) <- pieceData ]

  in \ t ->
    if | t <= 0
       -> 0
       | t >= 1
       -> totalLen
       | otherwise
       -> sum [ if b <= t then fullLen else intgr 0 ( toU t )
              | ( a, b, intgr, toU, fullLen ) <- pieceData
              , a < t
              ]
{-# INLINEABLE makeCompositeGLIntegrator #-}

--------------------------------------------------------------------------------
-- Gravesen adaptive subdivision

-- | Cumulative arc-length integrator for a quadratic Bézier using the
-- Gravesen (1993) adaptive polygon/chord method.
gravesenQuadIntegrator
  :: forall diff ptData
  .  ( Torsor diff ptData, Module Double diff, Inner Double diff )
  => Quadratic.Bezier ptData
  -> Double                   -- ^ @eps@: stopping tolerance
  -> ( Double -> Double )
gravesenQuadIntegrator bez eps =
  let !total = go eps bez
  in \ t ->
    if t <= 0 then 0
    else if t >= 1 then total
    else go eps $ fst $ Quadratic.subdivide @diff bez t
  where
    go :: Double -> Quadratic.Bezier ptData -> Double
    go e b@( Quadratic.Bezier { p0, p1, p2 } ) =
      let l0 = norm @diff ( p0 --> p2 )
          l1 = norm @diff ( p0 --> p1 ) + norm @diff ( p1 --> p2 )
      in if l1 - l0 < e
         then ( 2 * l0 + l1 ) / 3
         else let ( lb, rb ) = Quadratic.subdivide @diff b 0.5
              in go ( e / 2 ) lb + go ( e / 2 ) rb
{-# INLINEABLE gravesenQuadIntegrator #-}

-- | Cumulative arc-length integrator for a cubic Bézier using the
-- Gravesen (1993) adaptive polygon/chord method.
gravesenCubicIntegrator
  :: forall diff ptData
  .  ( Torsor diff ptData, Module Double diff, Inner Double diff )
  => Cubic.Bezier ptData
  -> Double                   -- ^ @eps@: stopping tolerance
  -> ( Double -> Double )
gravesenCubicIntegrator bez eps =
  let !total = go eps bez
  in \ t ->
    if t <= 0 then 0
    else if t >= 1 then total
    else go eps $ fst $ Cubic.subdivide @diff bez t
  where
    go :: Double -> Cubic.Bezier ptData -> Double
    go e b@( Cubic.Bezier { p0, p1, p2, p3 } ) =
      let l0 = norm @diff ( p0 --> p3 )
          l1 = norm @diff ( p0 --> p1 ) + norm @diff ( p1 --> p2 ) + norm @diff ( p2 --> p3 )
      in if l1 - l0 < e
         then ( l0 + l1 ) / 2       -- = (2·l0 + 2·l1) / 4, degree-3 Gravesen estimate
         else let ( lb, rb ) = Cubic.subdivide @diff b 0.5
              in go ( e / 2 ) lb + go ( e / 2 ) rb
{-# INLINEABLE gravesenCubicIntegrator #-}