diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,7 @@
+# 0.1.5
+
+* Add numeric tests
+
 # 0.1.4
 
 * Add `ArrDW`
diff --git a/massiv-test.cabal b/massiv-test.cabal
--- a/massiv-test.cabal
+++ b/massiv-test.cabal
@@ -1,5 +1,5 @@
 name:                massiv-test
-version:             0.1.4
+version:             0.1.5
 synopsis:            Library that contains generators, properties and tests for Massiv Array Library.
 description:         This library is designed for users of massiv library that need random generators for writing custom property tests and reusing some of the predefined ones.
 homepage:            https://github.com/lehins/massiv
@@ -22,6 +22,7 @@
                     , Test.Massiv.Core.Mutable
                     , Test.Massiv.Array.Delayed
                     , Test.Massiv.Array.Mutable
+                    , Test.Massiv.Array.Numeric
                     , Test.Massiv.Utils
 
 
@@ -59,9 +60,11 @@
                     , Test.Massiv.Array.Delayed.WindowedSpec
                     , Test.Massiv.Array.DelayedSpec
                     , Test.Massiv.Array.MutableSpec
-                    , Test.Massiv.Array.Ops.TransformSpec
-                    , Test.Massiv.Array.Ops.SortSpec
+                    , Test.Massiv.Array.Ops.FoldSpec
                     , Test.Massiv.Array.Ops.MapSpec
+                    , Test.Massiv.Array.Ops.SortSpec
+                    , Test.Massiv.Array.Ops.TransformSpec
+                    , Test.Massiv.Array.NumericSpec
                     , Test.Massiv.Array.Numeric.IntegralSpec
                     , Test.Massiv.Array.StencilSpec
                     , Test.Massiv.VectorSpec
@@ -69,7 +72,6 @@
                     , Data.Massiv.Array.Manifest.VectorSpec
                     , Data.Massiv.Array.ManifestSpec
                     , Data.Massiv.Array.Ops.ConstructSpec
-                    , Data.Massiv.Array.Ops.FoldSpec
                     , Data.Massiv.Array.Ops.SliceSpec
                     , Data.Massiv.ArraySpec
   build-depends:      base
diff --git a/src/Test/Massiv/Array/Numeric.hs b/src/Test/Massiv/Array/Numeric.hs
new file mode 100644
--- /dev/null
+++ b/src/Test/Massiv/Array/Numeric.hs
@@ -0,0 +1,307 @@
+{-# LANGUAGE AllowAmbiguousTypes #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE MonoLocalBinds #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
+module Test.Massiv.Array.Numeric
+  ( -- * Spec for safe Mutable instance
+    prop_MatrixMatrixMultiply
+  , mutableNumericSpec
+  , mutableNumericFloatSpec
+  ) where
+
+import Data.Massiv.Array as A
+import Test.Massiv.Utils as T
+import Test.Massiv.Core.Common ()
+
+
+naiveMatrixMatrixMultiply ::
+     (Num e, Source (R r1) Ix1 e, Source (R r2) Ix1 e, OuterSlice r1 Ix2 e, InnerSlice r2 Ix2 e)
+  => Array r1 Ix2 e
+  -> Array r2 Ix2 e
+  -> Array D Ix2 e
+naiveMatrixMatrixMultiply arr1 arr2
+  | n1 /= m2 =
+    error $
+    "(|*|): Inner array dimensions must agree, but received: " ++
+    show (size arr1) ++ " and " ++ show (size arr2)
+  | otherwise =
+    makeArrayR D Seq (Sz (m1 :. n2)) $ \(i :. j) ->
+      A.foldlS (+) 0 (A.zipWith (*) (arr1 !> i) (arr2 <! j))
+  where
+    Sz2 m1 n1 = size arr1
+    Sz2 m2 n2 = size arr2
+{-# INLINE naiveMatrixMatrixMultiply #-}
+
+
+prop_MatrixMatrixMultiply ::
+     forall r e. (Numeric r e, Mutable r Ix2 e, Eq e, Show e)
+  => Fun e e
+  -> Matrix r e
+  -> Property
+prop_MatrixMatrixMultiply f arr = expectProp $ do
+  let arr' = A.transpose (A.map (applyFun f) arr)
+  arr !><! arr' `shouldBe` naiveMatrixMatrixMultiply (delay arr) arr'
+  arr !><! transpose arr `shouldBe` naiveMatrixMatrixMultiply (delay arr) (transpose arr)
+  let Sz2 m n = size arr
+  when (m /= n) $
+    arr .><. arr `shouldThrow` (== SizeMismatchException (size arr) (Sz2 m n))
+
+prop_MatrixVectorMultiply ::
+     forall r e.
+     ( Numeric r e
+     , InnerSlice r Ix2 e
+     , Mutable r Ix2 e
+     , Source (R r) Ix1 e
+     , Source r Ix1 e
+     , Construct r Ix1 e
+     , Eq e
+     , Show e
+     )
+  => Fun Int e
+  -> Matrix r e
+  -> Property
+prop_MatrixVectorMultiply f arr =
+  expectProp $ do
+    let Sz2 _ n = size arr
+        v = makeArray Seq (Sz n) (applyFun f)
+    arr !>< v `shouldBe` flatten (naiveMatrixMatrixMultiply (delay arr) (resize' (Sz2 n 1) v))
+    arr .>< makeArray Seq (Sz (n + 1)) (applyFun f) `shouldThrow`
+      (== SizeMismatchException (size arr) (Sz2 (n + 1) 1))
+
+prop_VectorMatrixMultiply ::
+     forall r e.
+     ( Numeric r e
+     , OuterSlice r Ix2 e
+     , Mutable r Ix2 e
+     , Source (R r) Ix1 e
+     , Mutable r Ix1 e
+     , Eq e
+     , Show e
+     )
+  => Fun Int e
+  -> Matrix r e
+  -> Property
+prop_VectorMatrixMultiply f arr =
+  expectProp $ do
+    let Sz2 m _ = size arr
+        v = makeArray Seq (Sz m) (applyFun f)
+    v ><! arr `shouldBe` flatten (naiveMatrixMatrixMultiply (resize' (Sz2 1 m) v) (delay arr))
+    makeArray Seq (Sz (m + 1)) (applyFun f) ><. arr `shouldThrow`
+      (== SizeMismatchException (Sz2 1 (m + 1)) (size arr))
+
+prop_DotProduct ::
+     forall r e. (Numeric r e, Mutable r Ix1 e, Eq e, Show e)
+  => Fun e e
+  -> Vector r e
+  -> Property
+prop_DotProduct f v =
+  expectProp $ do
+    let v' = A.map (applyFun f) v
+    v !.! compute v' `shouldBe` A.sum (A.zipWith (*) v v')
+    dotM v (makeArray Seq (size v + 1) (const 0)) `shouldThrow`
+      (== SizeMismatchException (size v) (size v + 1))
+
+prop_Norm ::
+     forall r e. (NumericFloat r e, Mutable r Ix1 e, RealFloat e, Show e)
+  => e
+  -> Vector r e
+  -> Property
+prop_Norm eps v = epsilonEq eps (sqrt (v !.! v)) (normL2 v)
+
+
+
+prop_Plus ::
+     forall r e.
+     (Numeric r e, Mutable r Ix2 e, Show (Array r Ix2 e), Eq (Array r Ix2 e))
+  => Fun e e
+  -> Matrix r e
+  -> e
+  -> Property
+prop_Plus f arr e = expectProp $ do
+  arr .+ e `shouldBe` compute (A.map (+ e) arr)
+  e +. arr `shouldBe` arr .+ e
+  let arr' = compute (A.map (applyFun f) arr)
+  arr !+! arr' `shouldBe` compute (A.zipWith (+) arr arr')
+  let Sz2 m n = size arr
+  when (m /= n) $
+    arr .+. compute (transpose arr) `shouldThrow` (== SizeMismatchException (size arr) (Sz2 n m))
+
+prop_Minus ::
+     forall r e.
+     (Numeric r e, Mutable r Ix2 e, Show (Array r Ix2 e), Eq (Array r Ix2 e))
+  => Fun e e
+  -> Matrix r e
+  -> e
+  -> Property
+prop_Minus f arr e = expectProp $ do
+  arr .- e `shouldBe` compute (A.map (subtract e) arr)
+  e -. arr `shouldBe` negateA (arr .- e)
+  let arr' = compute (A.map (applyFun f) arr)
+  arr !-! arr' `shouldBe` compute (A.zipWith (-) arr arr')
+  let Sz2 m n = size arr
+  when (m /= n) $
+    arr .-. compute (transpose arr) `shouldThrow` (== SizeMismatchException (size arr) (Sz2 n m))
+
+prop_Times ::
+     forall r e.
+     (Numeric r e, Mutable r Ix2 e, Show (Array r Ix2 e), Eq (Array r Ix2 e))
+  => Fun e e
+  -> Matrix r e
+  -> e
+  -> Property
+prop_Times f arr e = expectProp $ do
+  arr .* e `shouldBe` compute (A.map (* e) arr)
+  e *. arr `shouldBe` arr .* e
+  let arr' = compute (A.map (applyFun f) arr)
+  arr !*! arr' `shouldBe` compute (A.zipWith (*) arr arr')
+  let Sz2 m n = size arr
+  when (m /= n) $
+    arr .*. compute (transpose arr) `shouldThrow` (== SizeMismatchException (size arr) (Sz2 n m))
+
+prop_Divide ::
+     forall r e.
+     ( NumericFloat r e
+     , Mutable r Ix2 e
+     , Show e
+     , RealFloat e
+     , Show (Array r Ix2 e)
+     , Eq (Array r Ix2 e)
+     )
+  => e -- ^ Epsilon
+  -> Fun e e
+  -> Matrix r e
+  -> e
+  -> Property
+prop_Divide eps f arr e = e /= 0 ==> expectProp $ do
+  arr ./ e `shouldBe` compute (A.map (/ e) arr)
+  epsilonFoldableExpect eps (delay (e /. arr)) (delay (e *. recipA arr))
+  let arr' = compute (A.map (applyFun f) arr)
+  unless (A.or (A.zipWith (\x y -> x == 0 && y == 0) arr arr')) $
+    arr !/! arr' `shouldBe` compute (A.zipWith (/) arr arr')
+  let Sz2 m n = size arr
+  when (m /= n) $
+    arr ./. compute (transpose arr) `shouldThrow` (== SizeMismatchException (size arr) (Sz2 n m))
+
+prop_Floating ::
+     forall r e. (RealFloat e, Source r Ix2 e, NumericFloat r e, Show e)
+  => e
+  -> Matrix r e
+  -> Property
+prop_Floating eps arr = expectProp $ do
+  epsilonFoldableExpect eps (delay (absA arr)) (A.map abs arr)
+  epsilonFoldableExpect eps (delay (signumA arr)) (A.map signum arr)
+  epsilonFoldableExpect eps (delay (recipA arr)) (A.map recip arr)
+  epsilonFoldableExpect eps (delay (expA arr)) (A.map exp arr)
+  epsilonFoldableExpect eps (delay (sqrtA arr)) (A.map sqrt arr)
+  epsilonFoldableExpect eps (delay (logA arr)) (A.map log arr)
+  epsilonFoldableExpect eps (delay (sinA arr)) (A.map sin arr)
+  epsilonFoldableExpect eps (delay (cosA arr)) (A.map cos arr)
+  epsilonFoldableExpect eps (delay (tanA arr)) (A.map tan arr)
+  epsilonFoldableExpect eps (delay (asinA arr)) (A.map asin arr)
+  epsilonFoldableExpect eps (delay (acosA arr)) (A.map acos arr)
+  epsilonFoldableExpect eps (delay (atanA arr)) (A.map atan arr)
+  epsilonFoldableExpect eps (delay (sinhA arr)) (A.map sinh arr)
+  epsilonFoldableExpect eps (delay (coshA arr)) (A.map cosh arr)
+  epsilonFoldableExpect eps (delay (tanhA arr)) (A.map tanh arr)
+  epsilonFoldableExpect eps (delay (asinhA arr)) (A.map asinh arr)
+  epsilonFoldableExpect eps (delay (acoshA arr)) (A.map acosh arr)
+  epsilonFoldableExpect eps (delay (atanhA arr)) (A.map atanh arr)
+
+prop_Floating2 ::
+     forall r e. (RealFloat e, Mutable r Ix2 e, NumericFloat r e, Show e)
+  => e
+  -> Matrix r e
+  -> Fun e e
+  -> Property
+prop_Floating2 eps arr1 f = expectProp $ do
+  let arr2 = compute (A.map (applyFun f) arr1)
+  epsilonFoldableExpect eps (delay (logBaseA arr1 arr2)) (A.zipWith logBase arr1 arr2)
+  epsilonFoldableExpect eps (delay (arr1 .** arr2)) (A.zipWith (**) arr1 arr2)
+  res <- atan2A arr1 arr2
+  epsilonFoldableExpect eps (delay res) (A.zipWith atan2 arr1 arr2)
+
+
+mutableNumericSpec ::
+     forall r e.
+     ( Numeric r e
+     , Mutable r Ix2 e
+     , InnerSlice r Ix2 e
+     , OuterSlice r Ix2 e
+     , Source (R r) Ix1 e
+     , Mutable r Ix1 e
+     , Eq e
+     , Show e
+     , Function e
+     , CoArbitrary e
+     , Arbitrary e
+     , Arbitrary (Array r Ix1 e)
+     , Arbitrary (Array r Ix2 e)
+     , Show (Array r Ix2 e)
+     , Eq (Array r Ix2 e)
+     , Show (Array r Ix1 e)
+     )
+  => Spec
+mutableNumericSpec =
+  describe "Numerc Operations" $ do
+    prop "Plus" $ prop_Plus @r @e
+    prop "Minus" $ prop_Minus @r @e
+    prop "Times" $ prop_Times @r @e
+    prop "DotProduct" $ prop_DotProduct @r @e
+    prop "Power" $ \(arr :: Array r Ix2 e) (NonNegative p) -> expectProp $
+      arr .^ p `shouldBe` compute (A.map (^ p) arr)
+    prop "MatrixMatrixMultiply" $ prop_MatrixMatrixMultiply @r @e
+    prop "MatrixVectorMultiply" $ prop_MatrixVectorMultiply @r @e
+    prop "VectorMatrixMultiply" $ prop_VectorMatrixMultiply @r @e
+    prop "Identity" $ \ n -> expectProp $ do
+      computeIO (identityMatrix (Sz n)) `shouldReturn`
+        makeArray @r Seq (Sz2 n n) (\ (i :. j) -> if i == j then 1 else 0 :: e)
+    prop "LowerTriangular" $ \ comp n f -> expectProp $ do
+      computeIO (lowerTriangular comp (Sz n) (applyFun f . fromIx2)) `shouldReturn`
+        makeArray @r Seq (Sz2 n n) (\ (i :. j) -> if i >= j then applyFun f (i, j) else 0 :: e)
+    prop "UpperTriangular" $ \ comp n f -> expectProp $ do
+      computeIO (upperTriangular comp (Sz n) (applyFun f . fromIx2)) `shouldReturn`
+        makeArray @r Seq (Sz2 n n) (\ (i :. j) -> if i <= j then applyFun f (i, j) else 0 :: e)
+
+mutableNumericFloatSpec ::
+     forall r.
+     ( NumericFloat r Float
+     , Mutable r Ix1 Float
+     , Mutable r Ix2 Float
+     , Arbitrary (Array r Ix1 Float)
+     , Arbitrary (Array r Ix2 Float)
+     , Show (Array r Ix1 Float)
+     , Show (Array r Ix2 Float)
+     , Eq (Array r Ix2 Float)
+     , NumericFloat r Double
+     , Mutable r Ix1 Double
+     , Mutable r Ix2 Double
+     , Arbitrary (Array r Ix1 Double)
+     , Arbitrary (Array r Ix2 Double)
+     , Show (Array r Ix1 Double)
+     , Show (Array r Ix2 Double)
+     , Eq (Array r Ix2 Double)
+     )
+  => Spec
+mutableNumericFloatSpec = do
+  let ef = 1e-6 :: Float
+      ed = 1e-12 :: Double
+  describe "NumericFloat Operations" $ do
+    describe "Float" $ do
+      prop "Divide" $ prop_Divide @r ef
+      prop "Floating" $ prop_Floating @r ef
+      prop "Floating2" $ prop_Floating2 @r ef
+      prop "Norm" $ prop_Norm @r ef
+      prop "Power" $ prop_Power @r ef
+    describe "Double" $ do
+      prop "Divide" $ prop_Divide @r ed
+      prop "Floating" $ prop_Floating @r ed
+      prop "Floating2" $ prop_Floating2 @r ed
+      prop "Norm" $ prop_Norm @r ed
+      prop "Power" $ prop_Power @r ed
+
+prop_Power ::
+     (Numeric r e, Source r Ix2 e, RealFloat e, Show e) => e -> Matrix r e -> Int -> Property
+prop_Power eps arr p = expectProp $
+  epsilonFoldableExpect eps (delay (arr .^^ p)) (A.map (^^ p) arr)
diff --git a/src/Test/Massiv/Utils.hs b/src/Test/Massiv/Utils.hs
--- a/src/Test/Massiv/Utils.hs
+++ b/src/Test/Massiv/Utils.hs
@@ -12,9 +12,18 @@
   , toStringException
   , ExpectedException(..)
   , applyFun2Compat
+  , expectProp
+  -- * Epsilon comparison
+  , epsilonExpect
+  , epsilonFoldableExpect
+  , epsilonMaybeEq
+  , epsilonEq
+  , epsilonEqDouble
+  , epsilonEqFloat
   , module X
   ) where
 
+import qualified Data.Foldable as F
 import Control.Monad as X
 import Control.Monad.ST as X
 import Data.Maybe as X (fromMaybe, isJust, isNothing)
@@ -93,3 +102,72 @@
 instance Function Word where
   function = functionMap fromIntegral fromInteger
 #endif
+
+-- | Convert an hspec Expectation to a quickcheck Property.
+--
+-- @since 1.5.0
+expectProp :: Expectation -> Property
+expectProp = monadicIO . run
+
+
+epsilonExpect ::
+     (HasCallStack, Show a, RealFloat a)
+  => a -- ^ Epsilon, a maximum tolerated error. Sign is ignored.
+  -> a -- ^ Expected result.
+  -> a -- ^ Tested value.
+  -> Expectation
+epsilonExpect epsilon x y =
+  X.forM_ (epsilonMaybeEq epsilon x y) $ \errMsg ->
+    expectationFailure $ "Expected: " ++ show x ++ " but got: " ++ show y ++ "\n   " ++ errMsg
+
+
+epsilonFoldableExpect ::
+     (HasCallStack, Foldable f, Show (f e), Show e, RealFloat e) => e -> f e -> f e -> Expectation
+epsilonFoldableExpect epsilon x y = do
+  F.length x `shouldBe` F.length y
+  unless (F.null x) $
+    X.forM_ (zipWithM (epsilonMaybeEq epsilon) (F.toList x) (F.toList y)) $ \errMsgs ->
+      expectationFailure $
+      "Expected: " ++ show x ++ " but got: " ++ show y ++ "\n" ++ unlines (fmap ("    " ++) errMsgs)
+
+
+epsilonMaybeEq ::
+     (Show a, RealFloat a)
+  => a -- ^ Epsilon, a maximum tolerated error. Sign is ignored.
+  -> a -- ^ Expected result.
+  -> a -- ^ Tested value.
+  -> Maybe String
+epsilonMaybeEq epsilon x y
+  | isNaN x && not (isNaN y) = Just $ "Expected NaN, but got: " ++ show y
+  | x == y = Nothing
+  | diff > n = Just $ concat [show x, " /= ", show y, " (Tolerance: ", show diff, " > ", show n, ")"]
+  | otherwise = Nothing
+  where
+    (absx, absy) = (abs x, abs y)
+    n = epsilon * (1 + max absx absy)
+    diff = abs (y - x)
+
+
+epsilonEq ::
+     (Show a, RealFloat a)
+  => a -- ^ Epsilon, a maximum tolerated error. Sign is ignored.
+  -> a -- ^ Expected result.
+  -> a -- ^ Tested value.
+  -> Property
+epsilonEq epsilon x y = property $ epsilonExpect epsilon x y
+
+epsilonEqDouble ::
+     Double -- ^ Expected result.
+  -> Double -- ^ Tested value.
+  -> Property
+epsilonEqDouble = epsilonEq epsilon
+  where
+    epsilon = 1e-12
+
+epsilonEqFloat ::
+     Float -- ^ Expected result.
+  -> Float -- ^ Tested value.
+  -> Property
+epsilonEqFloat = epsilonEq epsilon
+  where
+    epsilon = 1e-6
diff --git a/tests/Data/Massiv/Array/Ops/FoldSpec.hs b/tests/Data/Massiv/Array/Ops/FoldSpec.hs
deleted file mode 100644
--- a/tests/Data/Massiv/Array/Ops/FoldSpec.hs
+++ /dev/null
@@ -1,85 +0,0 @@
-{-# LANGUAGE AllowAmbiguousTypes #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE MonoLocalBinds #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE TypeApplications #-}
-{-# LANGUAGE TypeFamilies #-}
-module Data.Massiv.Array.Ops.FoldSpec (spec) where
-
-import qualified Data.Foldable as F
-import Data.Massiv.Array as A
-import Data.Semigroup
-import Prelude hiding (map, product, sum)
-import Test.Massiv.Core
-
-
-
-prop_SumSEqSumP :: Index ix => Array D ix Int -> Bool
-prop_SumSEqSumP arr = sum arr == sum (setComp Par arr)
-
-
-prop_ProdSEqProdP :: Index ix => Array D ix Int -> Bool
-prop_ProdSEqProdP arr = product arr == product (setComp Par arr)
-
-
-foldOpsProp ::
-     (Source P ix Int)
-  => Fun Int Bool
-  -> ArrTinyNE P ix Int
-  -> Property
-foldOpsProp f (ArrTinyNE arr) =
-  (A.maximum' arr === getMax (foldMono Max arr)) .&&.
-  (A.minimum' arr === getMin (foldSemi Min maxBound arr)) .&&.
-  (A.sum arr === F.sum ls) .&&.
-  (A.product (A.map ((+ 0.1) . (fromIntegral :: Int -> Double)) arr) ===
-   getProduct (foldMono (Product . (+ 0.1) . fromIntegral) arr)) .&&.
-  (A.all (apply f) arr === F.all (apply f) ls) .&&.
-  (A.any (apply f) arr === F.any (apply f) ls) .&&.
-  (A.or (A.map (apply f) arr) === F.or (fmap (apply f) ls)) .&&.
-  (A.and (A.map (apply f) arr) === F.and (fmap (apply f) ls))
-  where
-    ls = toList arr
-
-
-prop_NestedFoldP :: Array D Ix1 (Array D Ix1 Int) -> Bool
-prop_NestedFoldP arr = sum (setComp Par (map sum $ setComp Par arr)) == sum (map sum arr)
-
-
-specFold ::
-     forall ix. (Arbitrary ix, Index ix, Show (Array D ix Int), Show (Array P ix Int))
-  => String
-  -> Spec
-specFold dimStr =
-  describe dimStr $ do
-    it "sumS Eq sumP" $ property $ prop_SumSEqSumP @ix
-    it "prodS Eq prodP" $ property $ prop_ProdSEqProdP @ix
-    it "foldOps" $ property $ foldOpsProp @ix
-
-
-prop_foldOuterSliceToList ::
-     (Elt P ix Int ~ Array M (Lower ix) Int, OuterSlice P ix Int, Index (Lower ix))
-  => ArrTiny P ix Int
-  -> Property
-prop_foldOuterSliceToList (ArrTiny arr) =
-  foldOuterSlice A.toList arr === A.fold (A.map pure arr)
-
-
-spec :: Spec
-spec = do
-  specFold @Ix1 "Ix1"
-  specFold @Ix2 "Ix2"
-  specFold @Ix3 "Ix3"
-  it "Nested Parallel Fold" $ property prop_NestedFoldP
-  describe "Foldable Props" $ do
-    prop "Ix2" $ prop_foldOuterSliceToList @Ix2
-    prop "Ix3" $ prop_foldOuterSliceToList @Ix3
-    prop "Ix4" $ prop_foldOuterSliceToList @Ix4
-  describe "Exceptions" $ do
-    let emptySelector :: forall ix . Index ix => SizeException -> Bool
-        emptySelector = (== SizeEmptyException (Sz (zeroIndex :: ix)))
-    it "maximumM" $ maximumM (A.empty :: Array D Ix1 Int) `shouldThrow` emptySelector @Ix1
-    it "minimumM" $ minimumM (A.empty :: Array D Ix2 Int) `shouldThrow` emptySelector @Ix2
-    it "maximum'" $ (pure $! maximum' (A.empty :: Array D Ix3 Int)) `shouldThrow` emptySelector @Ix3
-    it "minimum'" $ (pure $! minimum' (A.empty :: Array D Ix4 Int)) `shouldThrow` emptySelector @Ix4
diff --git a/tests/Test/Massiv/Array/NumericSpec.hs b/tests/Test/Massiv/Array/NumericSpec.hs
new file mode 100644
--- /dev/null
+++ b/tests/Test/Massiv/Array/NumericSpec.hs
@@ -0,0 +1,18 @@
+{-# LANGUAGE TypeApplications #-}
+
+module Test.Massiv.Array.NumericSpec
+  ( spec
+  ) where
+
+import Data.Massiv.Array as A
+import Test.Massiv.Array.Numeric
+import Test.Massiv.Core
+
+spec :: Spec
+spec = do
+  mutableNumericSpec @P @Int
+  mutableNumericSpec @P @Float
+  mutableNumericFloatSpec @P
+  mutableNumericSpec @S @Int
+  mutableNumericSpec @S @Float
+  mutableNumericFloatSpec @S
diff --git a/tests/Test/Massiv/Array/Ops/FoldSpec.hs b/tests/Test/Massiv/Array/Ops/FoldSpec.hs
new file mode 100644
--- /dev/null
+++ b/tests/Test/Massiv/Array/Ops/FoldSpec.hs
@@ -0,0 +1,81 @@
+{-# LANGUAGE AllowAmbiguousTypes #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MonoLocalBinds #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
+{-# LANGUAGE TypeFamilies #-}
+module Test.Massiv.Array.Ops.FoldSpec (spec) where
+
+import qualified Data.Foldable as F
+import Data.Massiv.Array as A
+import Data.Semigroup
+import Prelude hiding (map, product, sum)
+import Test.Massiv.Core
+
+
+
+prop_SumSEqSumP :: Index ix => Array D ix Int -> Bool
+prop_SumSEqSumP arr = sum arr == sum (setComp Par arr)
+
+
+prop_ProdSEqProdP :: Index ix => Array D ix Int -> Bool
+prop_ProdSEqProdP arr = product arr == product (setComp Par arr)
+
+
+foldOpsProp :: Source P ix Int => Fun Int Bool -> ArrTinyNE P ix Int -> Expectation
+foldOpsProp f (ArrTinyNE arr) = do
+  A.maximum' arr `shouldBe` getMax (foldMono Max arr)
+  A.minimum' arr `shouldBe` getMin (foldSemi Min maxBound arr)
+  A.sum arr `shouldBe` F.sum ls
+  A.product (A.map ((+ 0.1) . (fromIntegral :: Int -> Double)) arr) `shouldBe`
+    getProduct (foldMono (Product . (+ 0.1) . fromIntegral) arr)
+  A.all (apply f) arr `shouldBe` F.all (apply f) ls
+  A.and (A.map (apply f) arr) `shouldBe` F.and (fmap (apply f) ls)
+  A.any (apply f) arr `shouldBe` F.any (apply f) ls
+  A.or (A.map (apply f) arr) `shouldBe` F.or (fmap (apply f) ls)
+  where
+    ls = toList arr
+
+
+prop_NestedFoldP :: Array D Ix1 (Array D Ix1 Int) -> Bool
+prop_NestedFoldP arr = sum (setComp Par (map sum $ setComp Par arr)) == sum (map sum arr)
+
+
+specFold ::
+     forall ix. (Arbitrary ix, Index ix, Show (Array D ix Int), Show (Array P ix Int))
+  => String
+  -> Spec
+specFold dimStr =
+  describe dimStr $ do
+    prop "sumS Eq sumP" $ prop_SumSEqSumP @ix
+    prop "prodS Eq prodP" $ prop_ProdSEqProdP @ix
+    prop "foldOps" $ foldOpsProp @ix
+
+
+prop_foldOuterSliceToList ::
+     (Elt P ix Int ~ Array M (Lower ix) Int, OuterSlice P ix Int, Index (Lower ix))
+  => ArrTiny P ix Int
+  -> Property
+prop_foldOuterSliceToList (ArrTiny arr) =
+  foldOuterSlice A.toList arr === A.fold (A.map pure arr)
+
+
+spec :: Spec
+spec = do
+  specFold @Ix1 "Ix1"
+  specFold @Ix2 "Ix2"
+  specFold @Ix3 "Ix3"
+  it "Nested Parallel Fold" $ property prop_NestedFoldP
+  describe "Foldable Props" $ do
+    prop "Ix2" $ prop_foldOuterSliceToList @Ix2
+    prop "Ix3" $ prop_foldOuterSliceToList @Ix3
+    prop "Ix4" $ prop_foldOuterSliceToList @Ix4
+  describe "Exceptions" $ do
+    let emptySelector :: forall ix . Index ix => SizeException -> Bool
+        emptySelector = (== SizeEmptyException (Sz (zeroIndex :: ix)))
+    it "maximumM" $ maximumM (A.empty :: Array D Ix1 Int) `shouldThrow` emptySelector @Ix1
+    it "minimumM" $ minimumM (A.empty :: Array D Ix2 Int) `shouldThrow` emptySelector @Ix2
+    it "maximum'" $ (pure $! maximum' (A.empty :: Array D Ix3 Int)) `shouldThrow` emptySelector @Ix3
+    it "minimum'" $ (pure $! minimum' (A.empty :: Array D Ix4 Int)) `shouldThrow` emptySelector @Ix4
diff --git a/tests/Test/Massiv/VectorSpec.hs b/tests/Test/Massiv/VectorSpec.hs
--- a/tests/Test/Massiv/VectorSpec.hs
+++ b/tests/Test/Massiv/VectorSpec.hs
@@ -828,12 +828,16 @@
             V.tail' arr !!==!! VP.tail (toPrimitiveVector arr)
           prop "take" $ \n (arr :: Array P Ix1 Word) ->
             V.take (Sz n) arr !==! VP.take n (toPrimitiveVector arr)
+          prop "takeWhile" $ \f (arr :: Array P Ix1 Word) ->
+            V.takeWhile (applyFun f) arr !==! VP.takeWhile (applyFun f) (toPrimitiveVector arr)
           prop "take'" $ \sz@(Sz n) (arr :: Array P Ix1 Word) ->
             V.take' sz arr !!==!! VP.slice 0 n (toPrimitiveVector arr)
           prop "stake" $ \n (arr :: Array P Ix1 Word) ->
             V.stake (Sz n) arr !==! VP.take n (toPrimitiveVector arr)
           prop "drop" $ \n (arr :: Array P Ix1 Word) ->
             V.drop (Sz n) arr !==! VP.drop n (toPrimitiveVector arr)
+          prop "dropWhile" $ \f (arr :: Array P Ix1 Word) ->
+            V.dropWhile (applyFun f) arr !==! VP.dropWhile (applyFun f) (toPrimitiveVector arr)
           prop "drop'" $ \sz@(Sz n) (arr :: Array P Ix1 Word) ->
             V.drop' sz arr !!==!! VP.slice n (unSz (size arr) - n) (toPrimitiveVector arr)
           prop "sdrop" $ \n (arr :: Array P Ix1 Word) ->
@@ -905,6 +909,9 @@
           prop "sconcat" $ \(vs :: [Vector P Int]) ->
             V.sconcat vs !==! VP.concat (fmap toPrimitiveVector vs)
       describe "Predicates" $ do
+        describe "Searching" $ do
+          prop "sfilter" $ \(v :: Vector P Word) (f :: Fun Word Bool) ->
+            V.findIndex (apply f) v === VP.findIndex (apply f) (toPrimitiveVector v)
         describe "Filtering" $ do
           prop "sfilter" $ \(v :: Vector P Word) (f :: Fun Word Bool) ->
             V.sfilter (apply f) v !==! VP.filter (apply f) (toPrimitiveVector v)