diff --git a/exact-real.cabal b/exact-real.cabal
--- a/exact-real.cabal
+++ b/exact-real.cabal
@@ -1,5 +1,5 @@
 name:         exact-real
-version:      0.8.0.2
+version:      0.9.0
 synopsis:     Exact real arithmetic
 description:
   A type to represent exact real number using a fast binary Cauchy sequence
diff --git a/src/Data/CReal/Internal.hs b/src/Data/CReal/Internal.hs
--- a/src/Data/CReal/Internal.hs
+++ b/src/Data/CReal/Internal.hs
@@ -28,9 +28,11 @@
   , recipBounded
   , shiftL
   , shiftR
+  , square
 
     -- ** Exponential
   , expBounded
+  , expPosNeg
   , logBounded
 
     -- ** Trigonometric
@@ -221,10 +223,12 @@
                  | otherwise -> atanBounded x
 
   -- TODO: benchmark replacing these with their series expansion
-  sinh x = (exp x - exp (-x)) / 2
-  cosh x = (exp x + exp (-x)) / 2
+  sinh x = let (expX, expNegX) = expPosNeg x
+           in (expX - expNegX) / 2
+  cosh x = let (expX, expNegX) = expPosNeg x
+           in (expX + expNegX) / 2
   tanh x = let e2x = exp (2 * x)
-           in (e2x - 1) / (e2x + 1)
+           in (e2x - 1) *. recipBounded (e2x + 1)
 
   asinh x = log (x + sqrt (x * x + 1))
   acosh x = log (x + sqrt (x + 1) * sqrt (x - 1))
@@ -348,14 +352,20 @@
                                     n = x (p + 2 * s + 2)
                                 in 2^(2 * p + 2 * s + 2) /. n)
 
+-- | Return the square of the input, more efficient than @('*')@
+square :: CReal n -> CReal n
+square (CR x) = CR (\p -> let s = log2 (abs (x 0) + 2) + 3
+                              n = x (p + s)
+                          in (n * n) /. 2^(p + 2 * s))
+
 --
--- Bounded exponential functions
+-- Bounded exponential functions and expPosNeg
 --
 
 -- | A more efficient 'exp' with the restriction that the input must be in the
 -- closed range [-1..1]
 expBounded :: CReal n -> CReal n
-expBounded x = let q = [1 % (n!) | n <- [0..]]
+expBounded x = let q = (1%) <$> scanl' (*) 1 [1..]
                in powerSeries q (max 5) x
 
 -- | A more efficient 'log' with the restriction that the input must be in the
@@ -365,6 +375,16 @@
                    y = (x - 1) .* recip x
                in y .* powerSeries q id y
 
+-- | @expPosNeg x@ returns @(exp x, exp (-x))#
+expPosNeg :: CReal n -> (CReal n, CReal n)
+expPosNeg x = let CR o = x / ln2
+                  l = o 0
+                  y = x - fromInteger l * ln2
+              in if l == 0
+                   then (expBounded x, expBounded (-x))
+                   else (expBounded y `shiftL` fromInteger l,
+                         expBounded (negate y) `shiftR` fromInteger l)
+
 --
 -- Bounded trigonometric functions
 --
@@ -393,6 +413,8 @@
 -- Multiplication with powers of two
 --
 
+infixl 8 `shiftL`, `shiftR`
+
 -- | @x \`shiftR\` n@ is equal to @x@ divided by 2^@n@
 --
 -- @n@ can be negative or zero
@@ -476,10 +498,6 @@
         satisfied r  = sq r <= x && sq (r + 1) > x
         initialGuess = 2 ^ (log2 x `div` 2)
         sq r         = r * r
-
--- | Factorial function
-(!) :: Integer -> Integer
-(!) x = product [2..x]
 
 --
 -- Searching
diff --git a/test/Test.hs b/test/Test.hs
--- a/test/Test.hs
+++ b/test/Test.hs
@@ -71,6 +71,12 @@
 test_boundedFunctions :: [TestTree]
 test_boundedFunctions = [ boundedFunctions (undefined :: CReal Precision) ]
 
+prop_expPosNeg :: CReal Precision -> Property
+prop_expPosNeg x = expPosNeg x === (exp x, exp (-x))
+
+prop_square :: CReal Precision -> Property
+prop_square x = square x === x * x
+
 main :: IO ()
 main = $(defaultMainGenerator)