diff --git a/src/matrix.rs b/src/matrix.rs
index b45dca2..d2af6f0 100644
--- a/src/matrix.rs
+++ b/src/matrix.rs
@@ -16,6 +16,7 @@
 use rand::{Rand, Rng};
 use num_traits::{cast, NumCast};
 use std::fmt;
+use std::iter;
 use std::mem;
 use std::ops::*;
 use std::ptr;
@@ -967,6 +968,34 @@ macro_rules! impl_matrix {
             fn sub_assign(&mut self, other: $MatrixN<S>) { $(self.$field -= other.$field);+ }
         }
 
+        impl<S: BaseFloat> iter::Sum for $MatrixN<S> {
+            #[inline]
+            fn sum<I: Iterator<Item=Self>>(iter: I) -> Self {
+                iter.fold(Self::zero(), Add::add)
+            }
+        }
+
+        impl<'a, S: 'a + BaseFloat> iter::Sum<&'a Self> for $MatrixN<S> {
+            #[inline]
+            fn sum<I: Iterator<Item=&'a Self>>(iter: I) -> Self {
+                iter.fold(Self::zero(), Add::add)
+            }
+        }
+
+        impl<S: BaseFloat> iter::Product for $MatrixN<S> {
+            #[inline]
+            fn product<I: Iterator<Item=Self>>(iter: I) -> Self {
+                iter.fold(Self::identity(), Mul::mul)
+            }
+        }
+
+        impl<'a, S: 'a + BaseFloat> iter::Product<&'a Self> for $MatrixN<S> {
+            #[inline]
+            fn product<I: Iterator<Item=&'a Self>>(iter: I) -> Self {
+                iter.fold(Self::identity(), Mul::mul)
+            }
+        }
+
         impl_scalar_ops!($MatrixN<usize> { $($field),+ });
         impl_scalar_ops!($MatrixN<u8> { $($field),+ });
         impl_scalar_ops!($MatrixN<u16> { $($field),+ });
diff --git a/src/quaternion.rs b/src/quaternion.rs
index 2c6189f..92be09e 100644
--- a/src/quaternion.rs
+++ b/src/quaternion.rs
@@ -13,6 +13,7 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
+use std::iter;
 use std::mem;
 use std::ops::*;
 
@@ -187,6 +188,34 @@ impl<S: BaseFloat> One for Quaternion<S> {
     }
 }
 
+impl<S: BaseFloat> iter::Sum for Quaternion<S> {
+    #[inline]
+    fn sum<I: Iterator<Item=Self>>(iter: I) -> Self {
+        iter.fold(Self::zero(), Add::add)
+    }
+}
+
+impl<'a, S: 'a + BaseFloat> iter::Sum<&'a Self> for Quaternion<S> {
+    #[inline]
+    fn sum<I: Iterator<Item=&'a Self>>(iter: I) -> Self {
+        iter.fold(Self::zero(), Add::add)
+    }
+}
+
+impl<S: BaseFloat> iter::Product for Quaternion<S> {
+    #[inline]
+    fn product<I: Iterator<Item=Self>>(iter: I) -> Self {
+        iter.fold(Self::one(), Mul::mul)
+    }
+}
+
+impl<'a, S: 'a + BaseFloat> iter::Product<&'a Self> for Quaternion<S> {
+    #[inline]
+    fn product<I: Iterator<Item=&'a Self>>(iter: I) -> Self {
+        iter.fold(Self::one(), Mul::mul)
+    }
+}
+
 impl<S: BaseFloat> VectorSpace for Quaternion<S> {
     type Scalar = S;
 }
diff --git a/src/rotation.rs b/src/rotation.rs
index 4a25101..6769634 100644
--- a/src/rotation.rs
+++ b/src/rotation.rs
@@ -14,6 +14,7 @@
 // limitations under the License.
 
 use std::fmt;
+use std::iter;
 use std::ops::*;
 
 use structure::*;
@@ -33,6 +34,7 @@ pub trait Rotation<P: EuclideanSpace>: Sized + Copy + One where
     // FIXME: Ugly type signatures - blocked by rust-lang/rust#24092
     Self: ApproxEq<Epsilon = <P as EuclideanSpace>::Scalar>,
     <P as EuclideanSpace>::Scalar: BaseFloat,
+    Self: iter::Product<Self>,
 {
     /// Create a rotation to a given direction with an 'up' vector.
     fn look_at(dir: P::Diff, up: P::Diff) -> Self;
@@ -157,6 +159,20 @@ impl<S: BaseFloat> From<Basis2<S>> for Matrix2<S> {
     fn from(b: Basis2<S>) -> Matrix2<S> { b.mat }
 }
 
+impl<S: BaseFloat> iter::Product for Basis2<S> {
+    #[inline]
+    fn product<I: Iterator<Item=Self>>(iter: I) -> Self {
+        iter.fold(Basis2 { mat: Matrix2::identity() }, Mul::mul)
+    }
+}
+
+impl<'a, S: 'a + BaseFloat> iter::Product<&'a Self> for Basis2<S> {
+    #[inline]
+    fn product<I: Iterator<Item=&'a Self>>(iter: I) -> Self {
+        iter.fold(Basis2 { mat: Matrix2::identity() }, Mul::mul)
+    }
+}
+
 impl<S: BaseFloat> Rotation<Point2<S>> for Basis2<S> {
     #[inline]
     fn look_at(dir: Vector2<S>, up: Vector2<S>) -> Basis2<S> {
@@ -263,6 +279,20 @@ impl<S: BaseFloat> From<Basis3<S>> for Quaternion<S> {
     fn from(b: Basis3<S>) -> Quaternion<S> { b.mat.into() }
 }
 
+impl<S: BaseFloat> iter::Product for Basis3<S> {
+    #[inline]
+    fn product<I: Iterator<Item=Self>>(iter: I) -> Self {
+        iter.fold(Basis3 { mat: Matrix3::identity() }, Mul::mul)
+    }
+}
+
+impl<'a, S: 'a + BaseFloat> iter::Product<&'a Self> for Basis3<S> {
+    #[inline]
+    fn product<I: Iterator<Item=&'a Self>>(iter: I) -> Self {
+        iter.fold(Basis3 { mat: Matrix3::identity() }, Mul::mul)
+    }
+}
+
 impl<S: BaseFloat> Rotation<Point3<S>> for Basis3<S> {
     #[inline]
     fn look_at(dir: Vector3<S>, up: Vector3<S>) -> Basis3<S> {
diff --git a/src/structure.rs b/src/structure.rs
index 5d7dc98..b29a1b3 100644
--- a/src/structure.rs
+++ b/src/structure.rs
@@ -17,6 +17,7 @@
 
 use num_traits::{cast, Float};
 use std::cmp;
+use std::iter;
 use std::ops::*;
 
 use approx::ApproxEq;
@@ -153,6 +154,7 @@ pub trait ElementWise<Rhs = Self> {
 /// ```
 pub trait VectorSpace: Copy + Clone where
     Self: Zero,
+    Self: iter::Sum<Self>,
 
     Self: Add<Self, Output = Self>,
     Self: Sub<Self, Output = Self>,
@@ -455,6 +457,7 @@ pub trait SquareMatrix where
     Self::Scalar: BaseFloat,
 
     Self: One,
+    Self: iter::Product,
 
     Self: Matrix<
         // FIXME: Can be cleaned up once equality constraints in where clauses are implemented
diff --git a/src/vector.rs b/src/vector.rs
index f645f3e..df3f7f8 100644
--- a/src/vector.rs
+++ b/src/vector.rs
@@ -16,6 +16,7 @@
 use rand::{Rand, Rng};
 use num_traits::NumCast;
 use std::fmt;
+use std::iter;
 use std::mem;
 use std::ops::*;
 
@@ -163,6 +164,20 @@ macro_rules! impl_vector {
             }
         }
 
+        impl<S: BaseNum> iter::Sum for $VectorN<S> {
+            #[inline]
+            fn sum<I: Iterator<Item=Self>>(iter: I) -> Self {
+                iter.fold(Self::zero(), Add::add)
+            }
+        }
+
+        impl<'a, S: 'a + BaseNum> iter::Sum<&'a Self> for $VectorN<S> {
+            #[inline]
+            fn sum<I: Iterator<Item=&'a Self>>(iter: I) -> Self {
+                iter.fold(Self::zero(), Add::add)
+            }
+        }
+
         impl<S: BaseNum> VectorSpace for $VectorN<S> {
             type Scalar = S;
         }
@@ -371,6 +386,20 @@ macro_rules! impl_vector_default {
             }
         }
 
+        impl<S: BaseNum> iter::Sum for $VectorN<S> {
+            #[inline]
+            fn sum<I: Iterator<Item=Self>>(iter: I) -> Self {
+                iter.fold(Self::zero(), Add::add)
+            }
+        }
+
+        impl<'a, S: 'a + BaseNum> iter::Sum<&'a Self> for $VectorN<S> {
+            #[inline]
+            fn sum<I: Iterator<Item=&'a Self>>(iter: I) -> Self {
+                iter.fold(Self::zero(), Add::add)
+            }
+        }
+
         impl<S: BaseNum> VectorSpace for $VectorN<S> {
             type Scalar = S;
         }
diff --git a/tests/matrix.rs b/tests/matrix.rs
index fc3cda3..43940d3 100644
--- a/tests/matrix.rs
+++ b/tests/matrix.rs
@@ -96,6 +96,18 @@ pub mod matrix2 {
         assert_eq!(A * B, &A * &B);
     }
 
+    #[test]
+    fn test_sum_matrix() {
+        let res: Matrix2<f64> = [A, B, C].iter().sum();
+        assert_eq!(res, A + B + C);
+    }
+
+    #[test]
+    fn test_product_matrix() {
+        let res: Matrix2<f64> = [A, B, C].iter().product();
+        assert_eq!(res, A * B * C);
+    }
+
     #[test]
     fn test_determinant() {
         assert_eq!(A.determinant(), -2.0f64)
@@ -258,6 +270,18 @@ pub mod matrix3 {
         assert_eq!(A * B, &A * &B);
     }
 
+    #[test]
+    fn test_sum_matrix() {
+        let res: Matrix3<f64> = [A, B, C, D].iter().sum();
+        assert_eq!(res, A + B + C + D);
+    }
+
+    #[test]
+    fn test_product_matrix() {
+        let res: Matrix3<f64> = [A, B, C, D].iter().product();
+        assert_eq!(res, A * B * C * D);
+    }
+
     #[test]
     fn test_determinant() {;
         assert_eq!(A.determinant(), 0.0f64);
@@ -615,6 +639,18 @@ pub mod matrix4 {
         assert_eq!(A * B, &A * &B);
     }
 
+    #[test]
+    fn test_sum_matrix() {
+        let res: Matrix4<f64> = [A, B, C, D].iter().sum();
+        assert_eq!(res, A + B + C + D);
+    }
+
+    #[test]
+    fn test_product_matrix() {
+        let res: Matrix4<f64> = [A, B, C, D].iter().product();
+        assert_eq!(res, A * B * C * D);
+    }
+
     #[test]
     fn test_determinant() {
         assert_eq!(A.determinant(), 0.0f64);
diff --git a/tests/quaternion.rs b/tests/quaternion.rs
index d2e9a8c..d2fbff5 100644
--- a/tests/quaternion.rs
+++ b/tests/quaternion.rs
@@ -52,6 +52,16 @@ mod operators {
     fn test_div() {
         impl_test_div!(2.0f32, Quaternion::from(Euler { x: Rad(1f32), y: Rad(1f32), z: Rad(1f32) }));
     }
+
+    #[test]
+    fn test_iter_product() {
+        let q1 = Quaternion::from(Euler { x: Rad(2f32), y: Rad(1f32), z: Rad(1f32) });
+        let q2 = Quaternion::from(Euler { x: Rad(1f32), y: Rad(2f32), z: Rad(1f32) });
+        let q3 = Quaternion::from(Euler { x: Rad(1f32), y: Rad(1f32), z: Rad(2f32) });
+
+        let res: Quaternion<f32> = [q1, q2, q3].iter().product();
+        assert_eq!(res, q1 * q2 * q3);
+    }
 }
 
 mod to_from_euler {
diff --git a/tests/vector.rs b/tests/vector.rs
index 307a5de..1bee749 100644
--- a/tests/vector.rs
+++ b/tests/vector.rs
@@ -20,6 +20,7 @@ extern crate cgmath;
 
 use cgmath::*;
 use std::f64;
+use std::iter;
 
 #[test]
 fn test_constructor() {
@@ -87,6 +88,14 @@ macro_rules! impl_test_rem {
     )
 }
 
+macro_rules! impl_test_iter_sum {
+    ($VectorN:ident { $($field:ident),+ }, $ty:ty, $s:expr, $v:expr) => (
+        let res: $VectorN<$ty> = iter::repeat($v).take($s as usize).sum();
+        assert_eq!(res,
+                   $VectorN::new($($v.$field * $s),+));
+    )
+}
+
 #[test]
 fn test_add() {
     impl_test_add!(Vector4 { x, y, z, w }, 2.0f32, vec4(2.0f32, 4.0, 6.0, 8.0));
@@ -140,6 +149,17 @@ fn test_sum() {
     assert_eq!(Vector4::new(5.0f64, 6.0f64, 7.0f64, 8.0f64).sum(), 26.0f64);
 }
 
+#[test]
+fn test_iter_sum() {
+    impl_test_iter_sum!(Vector4 { x, y, z, w }, f32, 2.0f32, vec4(2.0f32, 4.0, 6.0, 8.0));
+    impl_test_iter_sum!(Vector3 { x, y, z }, f32, 2.0f32, vec3(2.0f32, 4.0, 6.0));
+    impl_test_iter_sum!(Vector2 { x, y }, f32, 2.0f32, vec2(2.0f32, 4.0));
+
+    impl_test_iter_sum!(Vector4 { x, y, z, w }, usize, 2usize, vec4(2usize, 4, 6, 8));
+    impl_test_iter_sum!(Vector3 { x, y, z }, usize, 2usize, vec3(2usize, 4, 6));
+    impl_test_iter_sum!(Vector2 { x, y }, usize, 2usize, vec2(2usize, 4));
+}
+
 #[test]
 fn test_product() {
     assert_eq!(Vector2::new(1isize, 2isize).product(), 2isize);