content: {-# LANGUAGE TypeApplications     #-}
  
  -- Check that the standard analytic functions are correctly
  -- inverted by the corresponding inverse functions.
  
  main :: IO ()
  main = mapM_ print
   [ -- @recip@ is self-inverse on @ℝ\\{0}@.
     invDeviation @Double recip recip <$> [-1e20, -1e3, -1, -1e-40, 1e-40, 1e90]
   , invDeviation @Float  recip recip <$> [-1e10, -10, -1, -1e-20, 1e-20, 1e30]
   , -- @exp@ is invertible on @ℝ <-> [0…∞[@, but grows very fast.
     invDeviation @Double exp log <$> [-10, -5 .. 300]
   , invDeviation @Float  exp log <$> [-10 .. 60]
     -- @sin@ is only invertible on @[-π/2…π/2] <-> [-1…1]@.
   , invDeviation @Double sin asin <$> [-1.5, -1.4 .. 1.5]
   , invDeviation @Float  sin asin <$> [-1.5, -1.4 .. 1.5]
     -- @cos@ is invertible on @[0…π] <-> [-1…1]@.
   , invDeviation @Double cos acos <$> [0, 0.1 .. 3]
   , invDeviation @Float  cos acos <$> [0, 0.1 .. 3]
     -- @tan@ is invertible on @]-π/4…π/4[ <-> ]-∞…∞[@.
   , invDeviation @Double tan atan <$> [-0.7, -0.6 .. 0.7]
   , invDeviation @Float  tan atan <$> [-0.7, -0.6 .. 0.7]
     -- @sinh@ is invertible on @ℝ <-> ℝ@, but grows very fast.
   , invDeviation @Double sinh asinh <$> [-700, -672 .. 700]
   , invDeviation @Float  sinh asinh <$> [-80, -71 .. 80]
     -- @cosh@ is invertible on @[0…∞[ <-> [1…∞[@, but grows fast
   , invDeviation @Double cosh acosh <$> [0, 15 .. 700]
   , invDeviation @Float  cosh acosh <$> [0, 15 .. 80]
     -- @tanh@ is invertible on @ℝ <-> ]-1…1[@.
   , invDeviation @Double atanh tanh <$> [-0.99, -0.87 .. 0.9]
   , invDeviation @Float  atanh tanh <$> [-0.99, -0.87 .. 0.9]
   ]
  
  invDeviation :: KnownNumDeviation a
            => (a -> a) -- ^ Some numerical function @f@.
            -> (a -> a) -- ^ Inverse @g = f⁻¹@ of that function.
            -> a        -- ^ Value @x@ which to compare with @g (f x)@.
            -> Double   -- ^ Relative discrepancy between original/expected
                        --   value and actual function result.
  invDeviation f g 0 = rmNumericDeviation (g (f 0) + 1) - 1
  invDeviation f g x = rmNumericDeviation (g (f x) / x) - 1
  
  -- | We need to round results to some sensible precision,
  --   because floating-point arithmetic generally makes
  --   it impossible to /exactly/ invert functions.
  --   What precision this is depends on the type. The bounds
  --   here are rather generous; the functions should usually
  --   perform substantially better than that.
  class (Floating a, Eq a) => KnownNumDeviation a where
    rmNumericDeviation :: a -> Double
  
  instance KnownNumDeviation Double where
    rmNumericDeviation x = fromIntegral (round $ x * 2^36) / 2^36
  
  instance KnownNumDeviation Float where
    rmNumericDeviation x = fromIntegral (round $ x * 2^16) / 2^16
  