diff --git a/simgear/math/simd.hxx b/simgear/math/simd.hxx
index 26d7d84a..2a3f727d 100644
--- a/simgear/math/simd.hxx
+++ b/simgear/math/simd.hxx
@@ -20,306 +20,271 @@
 
 #include <cstring>
 
-# ifdef __SSE__
-#  include <xmmintrin.h>
-# endif
-# ifdef __SSE2__
-#  include <emmintrin.h>
-# endif
 
-
-template<typename T>
+template<typename T, int N>
 class simd4_t
 {
+private:
+    union {
+       T v4[4];
+       T vec[N];
+    };
+
 public:
+    simd4_t() {}
+    template<typename S>
+    simd4_t(S s) {
+        vec[0] = vec[1] = vec[2] = vec[3] = s;
+    }
+    explicit simd4_t(const T v[N]) {
+        std::memcpy(vec, v, sizeof(T[N]));
+    }
+    simd4_t(const simd4_t<T,N>& v) {
+        std::memcpy(vec, v.vec, sizeof(T[N]));
+    }
     ~simd4_t() {}
+
+    T (&ptr(void))[N] {
+        return vec;
+    }
+
+    const T (&ptr(void) const)[N] {
+        return vec;
+    }
+
+    inline operator const T*() const {
+        return vec;
+    }
+
+    inline operator T*() {
+        return vec;
+    }
+
+    template<typename S>
+    inline simd4_t<T,N>& operator=(S s) {
+        vec[0] = vec[1] = vec[2] = vec[3] = s;
+        return *this;
+    }
+    inline simd4_t<T,N>& operator=(const T v[N]) {
+        std::memcpy(vec, v, sizeof(T[N]));
+        return *this;
+    }
+    inline simd4_t<T,N>& operator=(const simd4_t<T,N>& v) {
+        std::memcpy(vec, v.vec, sizeof(T[N]));
+        return *this;
+    }
+
+    template<typename S>
+    inline simd4_t<T,N> operator+(S s)
+    {
+        simd4_t<T,N> r(*this);
+        r += s;
+        return r;
+    }
+    inline simd4_t<T,N> operator+(const T v[N])
+    {
+        simd4_t<T,N> r(v);
+        r += *this;
+        return r;
+    }
+    inline simd4_t<T,N> operator+(const simd4_t<T,N>& v)
+    {
+        simd4_t<T,N> r(v);
+        r += *this;
+        return r;
+    }
+
+    inline simd4_t<T,N> operator-()
+    {
+        simd4_t<T,N> r(0);
+        r -= *this;
+        return r;
+    }
+    template<typename S>
+    inline simd4_t<T,N> operator-(S s)
+    {
+        simd4_t<T,N> r(*this);
+        r -= s;
+        return r;
+    }
+    inline simd4_t<T,N> operator-(simd4_t<T,N>& v)
+    {
+        simd4_t<T,N> r(*this);
+        r -= v;
+        return r;
+    }
+
+    template<typename S>
+    inline simd4_t<T,N> operator*(S s)
+    {
+        simd4_t<T,N> r(s);
+        r *= *this;
+        return r;
+    }
+    inline simd4_t<T,N> operator*(const T v[N])
+    {
+        simd4_t<T,N> r(v);
+        r *= *this;
+        return r;
+    }
+    inline simd4_t<T,N> operator*(simd4_t<T,N>& v)
+    {
+        simd4_t<T,N> r(v);
+        r *= *this;
+        return r;
+    }
+
+    template<typename S>
+    inline simd4_t<T,N> operator/(S s)
+    {
+        simd4_t<T,N> r(1/T(s));
+        r *= this;
+        return r;
+    }
+    inline simd4_t<T,N> operator/(const T v[N])
+    {
+        simd4_t<T,N> r(*this);
+        r /= v;
+        return r;
+    }
+    inline simd4_t<T,N> operator/(simd4_t<T,N>& v)
+    {
+        simd4_t<T,N> r(*this);
+        r /= v;
+        return r;
+    }
+
+    template<typename S>
+    inline simd4_t<T,N>& operator+=(S s)
+    {
+        for (int i=0; i<N; i++) {
+            vec[i] += s;
+        }
+        return *this;
+    }
+    inline simd4_t<T,N>& operator+=(const T v[N])
+    {
+        simd4_t<T,N> r(v);
+        *this += r;
+        return *this;
+    }
+    inline simd4_t<T,N>& operator+=(const simd4_t<T,N>& v)
+    {
+        for (int i=0; i<N; i++) {
+           vec[i] += v[i];
+        }
+        return *this;
+    }
+
+    template<typename S>
+    inline simd4_t<T,N>& operator-=(S s)
+    {
+        for (int i=0; i<N; i++) {
+           vec[i] -= s;
+        }
+        return *this;
+    }
+ 
+    inline simd4_t<T,N>& operator-=(const T v[N])
+    {
+        simd4_t<T,N> r(v);
+        *this -= r;
+        return *this;
+    }
+    inline simd4_t<T,N>& operator-=(const simd4_t<T,N>& v)
+    {
+        for (int i=0; i<N; i++) {
+            vec[i] -= v[i];
+        }
+        return *this;
+    }
+
+    template<typename S>
+    inline simd4_t<T,N>& operator *=(S s)
+    {
+        for (int i=0; i<N; i++) {
+           vec[i] *= s;
+        }
+        return *this;
+    }
+    inline simd4_t<T,N>& operator*=(const T v[N])
+    {
+        simd4_t<T,N> r(v);
+        *this *= r;
+        return *this;
+    }
+    inline simd4_t<T,N>& operator*=(const simd4_t<T,N>& v)
+    {
+        for (int i=0; i<N; i++) {
+           vec[i] *= v[i];
+        }
+        return *this;
+    }
+
+    template<typename S>
+    inline simd4_t<T,N>& operator/=(S s)
+    {
+        s = (1/s);
+        for (int i=0; i<N; i++) {
+           vec[i] *= s;
+        }
+        return *this;
+    }
+    inline simd4_t<T,N>& operator/=(const T v[N])
+    {
+        simd4_t<T,N> r(v);
+        *this /= r;
+        return *this;
+    }
+    inline simd4_t<T,N>& operator/=(const simd4_t<T,N>& v)
+    {
+        for (int i=0; i<N; i++) {
+           vec[i] /= v[i];
+        }
+        return *this;
+    }
+
 };
 
-template<>
-class simd4_t<float>
+# ifdef __SSE__
+#  include <xmmintrin.h>
+
+template<int N>
+class simd4_t<float,N>
 {
 private:
-   typedef float  __vec4f_t[4];
+   typedef float  __vec4f_t[N];
 
-# ifdef __SSE__
     union {
         __m128 simd4;
         __vec4f_t vec;
     };
-# else
-    __vec4f_t vec;
-# endif
 
 public:
     simd4_t() {}
-    simd4_t(float f)
-    {
-# ifdef __SSE__
+    simd4_t(float f) {
         simd4 = _mm_set1_ps(f);
-# else
-        vec[0] = vec[1] = vec[2] = vec[3] = f;
-# endif
     }
-    explicit simd4_t(const __vec4f_t v)
-    {
-# ifdef __SSE__
+    explicit simd4_t(const __vec4f_t v) {
         simd4 = _mm_loadu_ps(v);
-# else
-        std::memcpy(vec, v, sizeof(float[4]));
-# endif
     }
-    simd4_t(const simd4_t<float>& v)
-    {
-# ifdef __SSE__
+    simd4_t(const simd4_t<float,N>& v) {
         simd4 = v.simd4;
-# else
-        std::memcpy(vec, v.vec, sizeof(float[4]));
-# endif
     }
+#if 0
+    simd4_t(const simd4_t<double,N>& v) {
+        simd4 = _mm_movelh_ps(_mm_cvtpd_ps(v.v4()[0]), _mm_cvtpd_ps(v.v4()[1]));
+    }
+#endif
 
-    inline float (&ptr(void))[4] {
+    float (&ptr(void))[N] {
         return vec;
     }
 
-    inline const float (&ptr(void) const)[4] {
+    const float (&ptr(void) const)[N] {
         return vec;
     }
 
-    inline simd4_t<float>& operator=(float f)
-    {
-# ifdef __SSE__
-        simd4 = _mm_set1_ps(f);
-# else
-        vec[0] = vec[1] = vec[2] = vec[3] = f;
-# endif
-        return *this;
-    }
-    inline simd4_t<float>& operator=(const __vec4f_t v)
-    {
-# ifdef __SSE__
-        simd4 = _mm_loadu_ps(v);
-# else
-        std::memcpy(vec, v, sizeof(float[4]));
-# endif
-        return *this;
-    }
-    inline simd4_t<float>& operator=(const simd4_t<float>& v)
-    {
-# ifdef __SSE__
-        simd4 = v.simd4;
-# else
-        std::memcpy(vec, v.vec, sizeof(float[4]));
-# endif
-        return *this;
-    }
-
-    inline simd4_t operator+(float f)
-    {
-        simd4_t r(*this);
-        r += f;
-        return r;
-    }
-    inline simd4_t operator+(const __vec4f_t v)
-    {
-        simd4_t r(*this);
-        r += v;
-        return r;
-    }
-    inline simd4_t operator+(const simd4_t<float>& v)
-    {
-        simd4_t r(*this);
-        r += v;
-        return r;
-    }
-
-    inline simd4_t operator-()
-    {
-        simd4_t r(0.0f);
-        r -= vec;
-        return r;
-    }
-    inline simd4_t operator-(float f)
-    {
-        simd4_t r(*this);
-        r -= f;
-        return r;
-    }
-    inline simd4_t operator-(const __vec4f_t v)
-    {
-        simd4_t r(*this);
-        r -= v;
-        return r;
-    }
-    inline simd4_t operator-(simd4_t<float>& v)
-    {
-        simd4_t r(*this);
-        r -= v;
-        return r;
-    }
-
-    inline simd4_t operator*(float f)
-    {
-        simd4_t r(*this);
-        r *= f;
-        return r;
-    }
-    inline simd4_t operator*(const __vec4f_t v)
-    {
-        simd4_t r(*this);
-        r *= v;
-        return r;
-    }
-    inline simd4_t operator*(simd4_t<float>& v)
-    {
-        simd4_t r(*this);
-        r *= v;
-        return r;
-    }
-
-    inline simd4_t operator/(float f)
-    {
-        simd4_t r(*this);
-        r /= f;
-        return r;
-    }
-    inline simd4_t operator/(const __vec4f_t v)
-    {
-        simd4_t r(*this);
-        r /= v;
-        return r;
-    }
-    inline simd4_t operator/(simd4_t<float>& v)
-    {
-        simd4_t r(*this);
-        r /= v;
-        return r;
-    }
-
-    inline simd4_t<float>& operator+=(float f)
-    {
-# ifdef __SSE__
-        simd4 += f;
-# else
-        vec[0] += f;
-        vec[1] += f;
-        vec[2] += f;
-        vec[3] += f;
-# endif
-        return *this;
-    }
-    inline simd4_t<float>& operator+=(const __vec4f_t v)
-    {
-        simd4_t r(v);
-        *this += r;
-        return *this;
-    }
-    inline simd4_t<float>& operator+=(const simd4_t<float>& v)
-    {
-# ifdef __SSE__
-        simd4 += v.simd4;
-# else
-        vec[0] += v[0];
-        vec[1] += v[1];
-        vec[2] += v[2];
-        vec[3] += v[3];
-# endif
-        return *this;
-    }
-
-    inline simd4_t<float>& operator-=(float f)
-    {
-# ifdef __SSE__
-        simd4 -= f;
-# else
-        vec[0] -= f;
-        vec[1] -= f;
-        vec[2] -= f;
-        vec[3] -= f;
-# endif
-        return *this;
-    }
-    inline simd4_t<float>& operator-=(const __vec4f_t v)
-    {
-        simd4_t r(v);
-        *this -= r;
-        return *this;
-    }
-    inline simd4_t<float>& operator-=(const simd4_t<float>& v)
-    {
-# ifdef __SSE__
-        simd4 -= v.simd4;
-# else
-        vec[0] -= v[0];
-        vec[1] -= v[1];
-        vec[2] -= v[2];
-        vec[3] -= v[3];
-# endif
-        return *this;
-    }
-
-    inline simd4_t<float>& operator *=(float f)
-    {
-# ifdef __SSE__
-        simd4 *= f;
-# else
-        vec[0] *= f;
-        vec[1] *= f;
-        vec[2] *= f;
-        vec[3] *= f;
-# endif
-        return *this;
-    }
-    inline simd4_t<float>& operator*=(const __vec4f_t v)
-    {
-        simd4_t r(v);
-        *this *= r;
-        return *this;
-    }
-    inline simd4_t<float>& operator*=(const simd4_t<float>& v)
-    {
-# ifdef __SSE__
-        simd4 *= v.simd4;
-# else
-        vec[0] *= v[0];
-        vec[1] *= v[1];
-        vec[2] *= v[2];
-        vec[3] *= v[3];
-# endif
-        return *this;
-    }
-
-    inline simd4_t<float>& operator/=(float f)
-    {
-        f = (1.0f/f);
-# ifdef __SSE__
-        simd4 *= f;
-# else
-        vec[0] *= f;
-        vec[1] *= f;
-        vec[2] *= f;
-        vec[3] *= f;
-# endif
-        return *this;
-    }
-    inline simd4_t<float>& operator/=(const __vec4f_t v)
-    {
-        simd4_t r(v);
-        *this /= r;
-        return *this;
-    }
-    inline simd4_t<float>& operator/=(const simd4_t<float>& v)
-    {
-# ifdef __SSE__
-        simd4 /= v.simd4;
-# else
-        vec[0] /= v[0];
-        vec[1] /= v[1];
-        vec[2] /= v[2];
-        vec[3] /= v[3];
-# endif
-        return *this;
-    }
-
     inline operator const float*() const {
         return vec;
     }
@@ -327,294 +292,102 @@ public:
     inline operator float*() {
         return vec;
     }
+
+    inline simd4_t<float,N>& operator=(float f) {
+        simd4 = _mm_set1_ps(f);
+        return *this;
+    }
+    inline simd4_t<float,N>& operator=(const __vec4f_t v) {
+        simd4 = _mm_loadu_ps(v);
+        return *this;
+    }
+    inline simd4_t<float,N>& operator=(const simd4_t<float,N>& v) {
+        simd4 = v.simd4;
+        return *this;
+    }
+
+    inline simd4_t<float,N>& operator+=(float f) {
+        simd4 += f;
+        return *this;
+    }
+    inline simd4_t<float,N>& operator+=(const simd4_t<float,N>& v) {
+        simd4 += v.simd4;
+        return *this;
+    }
+
+    inline simd4_t<float,N>& operator-=(float f) {
+        simd4 -= f;
+        return *this;
+    }
+    inline simd4_t<float,N>& operator-=(const simd4_t<float,N>& v) {
+        simd4 -= v.simd4;
+        return *this;
+    }
+
+    inline simd4_t<float,N>& operator *=(float f) {
+        simd4 *= f;
+        return *this;
+    }
+    inline simd4_t<float,N>& operator*=(const simd4_t<float,N>& v) {
+        simd4 *= v.simd4;
+        return *this;
+    }
+
+    inline simd4_t<float,N>& operator/=(float f) {
+        f = (1.0f/f);
+        simd4 *= f;
+        return *this;
+    }
+    inline simd4_t<float,N>& operator/=(const simd4_t<float,N>& v) {
+        simd4 /= v.simd4;
+        return *this;
+    }
 };
+# endif
 
 
-template<>
-class simd4_t<double>
+# ifdef __SSE2
+template<int N>
+class simd4_t<double,N>
 {
 private:
-   typedef double  __vec4d_t[4];
+   typedef double  __vec4d_t[N];
 
-# ifdef __SSE2__
     union {
-        __m128d simd4;
+        __m128d simd4[2];
         __vec4d_t vec;
     };
-# else
-    __vec4d_t vec;
-# endif
 
 public:
     simd4_t() {}
-    simd4_t(double d)
-    {
-# ifdef __SSE2__
-        simd4 = _mm_set1_pd(d);
-# else
-        vec[0] = vec[1] = vec[2] = vec[3] = d;
-# endif
+    simd4_t(double d) {
+        simd4[0] = simd4[1] = _mm_set1_pd(d);
     }
-    explicit simd4_t(const __vec4d_t v)
-    {
-# ifdef __SSE2__
-        simd4 = _mm_loadu_pd(v);
-# else
-        std::memcpy(vec, v, sizeof(double[4]));
-# endif
+    explicit simd4_t(const __vec4d_t v) {
+        simd4[0] = _mm_loadu_pd(v);
+        simd4[1] = _mm_loadu_pd(v+2);
     }
-    simd4_t(const simd4_t<double>& v)
-    {
-# ifdef __SSE2__
-        simd4 = v.simd4;
-# else
-        std::memcpy(vec, v.vec, sizeof(double[4]));
-# endif
+    simd4_t(const simd4_t<double,N>& v) {
+        simd4[0] = v.simd4[0];
+        simd4[1] = v.simd4[1];
     }
 
-    inline double (&ptr(void))[4] {
+    inline __m128d (&v4(void))[2] {
+        return simd4;
+    }
+
+    inline const __m128d (&v4(void) const)[2] {
+        return simd4;
+    }
+
+    inline double (&ptr(void))[N] {
         return vec;
     }
 
-    inline const double (&ptr(void) const)[4] {
+    inline const double (&ptr(void) const)[N] {
         return vec;
     }
 
-    inline simd4_t<double>& operator=(double d)
-    {
-# ifdef __SSE2__
-        simd4 = _mm_set1_pd(d);
-# else
-        vec[0] = vec[1] = vec[2] = vec[3] = d;
-# endif
-        return *this;
-    }
-    inline simd4_t<double>& operator=(const __vec4d_t v)
-    {
-# ifdef __SSE2__
-        simd4 = _mm_loadu_pd(v);
-# else
-        std::memcpy(vec, v, sizeof(double[4]));
-# endif
-        return *this;
-    }
-    inline simd4_t<double>& operator=(const simd4_t<double>& v)
-    {
-# ifdef __SSE2__
-        simd4 = v.simd4;
-# else
-        std::memcpy(vec, v.vec, sizeof(double[4]));
-# endif
-        return *this;
-    }
-
-    inline simd4_t operator+(double d)
-    {
-        simd4_t r(*this);
-        r += d;
-        return r;
-    }
-    inline simd4_t operator+(const __vec4d_t v)
-    {
-        simd4_t r(*this);
-        r += v;
-        return r;
-    }
-    inline simd4_t operator+(const simd4_t<double>& v)
-    {
-        simd4_t r(*this);
-        r += v;
-        return r;
-    }
-
-    inline simd4_t operator-()
-    {
-        simd4_t r(0.0f);
-        r -= vec;
-        return r;
-    }
-    inline simd4_t operator-(double d)
-    {
-        simd4_t r(*this);
-        r -= d;
-        return r;
-    }
-    inline simd4_t operator-(const __vec4d_t v)
-    {
-        simd4_t r(*this);
-        r -= v;
-        return r;
-    }
-    inline simd4_t operator-(simd4_t<double>& v)
-    {
-        simd4_t r(*this);
-        r -= v;
-        return r;
-    }
-
-    inline simd4_t operator*(double d)
-    {
-        simd4_t r(*this);
-        r *= d;
-        return r;
-    }
-    inline simd4_t operator*(const __vec4d_t v)
-    {
-        simd4_t r(*this);
-        r *= v;
-        return r;
-    }
-    inline simd4_t operator*(simd4_t<double>& v)
-    {
-        simd4_t r(*this);
-        r *= v;
-        return r;
-    }
-
-    inline simd4_t operator/(double d)
-    {
-        simd4_t r(*this);
-        r /= d;
-        return r;
-    }
-    inline simd4_t operator/(const __vec4d_t v)
-    {
-        simd4_t r(*this);
-        r /= v;
-        return r;
-    }
-    inline simd4_t operator/(simd4_t<double>& v)
-    {
-        simd4_t r(*this);
-        r /= v;
-        return r;
-    }
-
-    inline simd4_t<double>& operator+=(double d)
-    {
-# ifdef __SSE2__
-        simd4 += d;
-# else
-        vec[0] += d;
-        vec[1] += d;
-        vec[2] += d;
-        vec[3] += d;
-# endif
-        return *this;
-    }
-    inline simd4_t<double>& operator+=(const __vec4d_t v)
-    {
-        simd4_t r(v);
-        *this += r;
-        return *this;
-    }
-    inline simd4_t<double>& operator+=(const simd4_t<double>& v)
-    {
-# ifdef __SSE2__
-        simd4 += v.simd4;
-# else
-        vec[0] += v[0];
-        vec[1] += v[1];
-        vec[2] += v[2];
-        vec[3] += v[3];
-# endif
-        return *this;
-    }
-
-    inline simd4_t<double>& operator-=(double d)
-    {
-# ifdef __SSE2__
-        simd4 -= d;
-# else
-        vec[0] -= d;
-        vec[1] -= d;
-        vec[2] -= d;
-        vec[3] -= d;
-# endif
-        return *this;
-    }
-    inline simd4_t<double>& operator-=(const __vec4d_t v)
-    {
-        simd4_t r(v);
-        *this -= r;
-        return *this;
-    }
-    inline simd4_t<double>& operator-=(const simd4_t<double>& v)
-    {
-# ifdef __SSE2__
-        simd4 -= v.simd4;
-# else
-        vec[0] -= v[0];
-        vec[1] -= v[1];
-        vec[2] -= v[2];
-        vec[3] -= v[3];
-# endif
-        return *this;
-    }
-
-    inline simd4_t<double>& operator *=(double d)
-    {
-# ifdef __SSE2__
-        simd4 *= d;
-# else
-        vec[0] *= d;
-        vec[1] *= d;
-        vec[2] *= d;
-        vec[3] *= d;
-# endif
-        return *this;
-    }
-    inline simd4_t<double>& operator*=(const __vec4d_t v)
-    {
-        simd4_t r(v);
-        *this *= r;
-        return *this;
-    }
-    inline simd4_t<double>& operator*=(const simd4_t<double>& v)
-    {
-# ifdef __SSE2__
-        simd4 *= v.simd4;
-# else
-        vec[0] *= v[0];
-        vec[1] *= v[1];
-        vec[2] *= v[2];
-        vec[3] *= v[3];
-# endif
-        return *this;
-    }
-
-    inline simd4_t<double>& operator/=(double d)
-    {
-# ifdef __SSE2__
-        d = (1.0/d);
-        simd4 *= d;
-# else
-        vec[0] /= d;
-        vec[1] /= d;
-        vec[2] /= d;
-        vec[3] /= d;
-# endif
-        return *this;
-    }
-    inline simd4_t<double>& operator/=(const __vec4d_t v)
-    {
-        simd4_t r(v);
-        *this /= r;
-        return *this;
-    }
-    inline simd4_t<double>& operator/=(const simd4_t<double>& v)
-    {
-# ifdef __SSE2__
-        simd4 /= v.simd4;
-# else
-        vec[0] /= v[0];
-        vec[1] /= v[1];
-        vec[2] /= v[2];
-        vec[3] /= v[3];
-# endif
-        return *this;
-    }
-
     inline operator const double*() const {
         return vec;
     }
@@ -622,294 +395,104 @@ public:
     inline operator double*() {
         return vec;
     }
+
+    inline simd4_t<double,N>& operator=(double d) {
+        simd4[0] = simd[1] = _mm_set1_pd(d);
+        return *this;
+    }
+    inline simd4_t<double,N>& operator=(const __vec4d_t v) {
+        simd4[0] = _mm_loadu_pd(v);
+        simd4[1] = _mm_loadu_pd(v+2);
+        return *this;
+    }
+    inline simd4_t<double,N>& operator=(const simd4_t<double,N>& v) {
+        simd4[0] = v.simd4[0];
+        simd4[1] = v.simd4[1];
+        return *this;
+    }
+
+    inline simd4_t<double,N>& operator+=(double d) {
+        simd4[0] += d;
+        simd4[1] += d;
+        return *this;
+    }
+    inline simd4_t<double,N>& operator+=(const simd4_t<double,N>& v) {
+        simd4[0] += v.simd4[0];
+        simd4[1] += v.simd4[1];
+        return *this;
+    }
+
+    inline simd4_t<double,N>& operator-=(double d) {
+        simd4[0] -= d;
+        simd4[1] -= d;
+        return *this;
+    }
+    inline simd4_t<double,N>& operator-=(const simd4_t<double,N>& v) {
+        simd4[0] -= v.simd4[0];
+        simd4[1] -= v.simd4[1];
+        return *this;
+    }
+
+    inline simd4_t<double,N>& operator *=(double d) {
+        simd4[0] *= d;
+        simd4[1] *= d;
+        return *this;
+    }
+    inline simd4_t<double,N>& operator*=(const simd4_t<double,N>& v) {
+        simd4[0] *= v.simd4[0];
+        simd4[1] *= v.simd4[1];
+        return *this;
+    }
+
+    inline simd4_t<double,N>& operator/=(double d) {
+        d = (1.0/d);
+        simd4[0] *= d;
+        simd4[1] *= d;
+        return *this;
+    }
+    inline simd4_t<double,N>& operator/=(const simd4_t<double,N>& v) {
+        simd4[0] /= v.simd4[0];
+        simd4[1] /= v.simd4[1];
+        return *this;
+    }
 };
+# endif
 
 
-template<>
-class simd4_t<int>
-{
-private:
-   typedef int  __vec4i_t[4];
-
 # ifdef __SSE2__
+#  include <emmintrin.h>
+
+template<int N>
+class simd4_t<int,N>
+{
+private:
+   typedef int  __vec4i_t[N];
+
     union {
         __m128i simd4;
         __vec4i_t vec;
     };
-# else
-    __vec4i_t vec;
-# endif
 
 public:
     simd4_t() {}
-    simd4_t(int i)
-    {
-# ifdef __SSE2__
+    simd4_t(int i) {
         simd4 = _mm_set1_epi32(i);
-# else
-        vec[0] = vec[1] = vec[2] = vec[3] = i;
-# endif
     }
-    explicit simd4_t(const __vec4i_t v)
-    {
-# ifdef __SSE2__
+    explicit simd4_t(const __vec4i_t v) {
         simd4 = _mm_loadu_si128((__m128i*)v);
-# else
-        std::memcpy(vec, v, sizeof(int[4]));
-# endif
     }
-    simd4_t(const simd4_t<int>& v)
-    {
-# ifdef __SSE2__
+    simd4_t(const simd4_t<int,N>& v) {
         simd4 = v.simd4;
-# else
-        std::memcpy(vec, v.vec, sizeof(int[4]));
-# endif
     }
 
-    inline int (&ptr(void))[4] {
+    int (&ptr(void))[4] {
         return vec;
     }
 
-    inline const int (&ptr(void) const)[4] {
+    const int (&ptr(void) const)[4] {
         return vec;
     }
 
-    inline simd4_t<int>& operator=(int i)
-    {
-# ifdef __SSE2__
-        simd4 = _mm_set1_epi32(i);
-# else
-        vec[0] = vec[1] = vec[2] = vec[3] = i;
-# endif
-        return *this;
-    }
-    inline simd4_t<int>& operator=(const __vec4i_t v)
-    {
-# ifdef __SSE2__
-        simd4 = _mm_loadu_si128((__m128i*)v);
-# else
-        std::memcpy(vec, v, sizeof(int[4]));
-# endif
-        return *this;
-    }
-    inline simd4_t<int>& operator=(const simd4_t<int>& v)
-    {
-# ifdef __SSE2__
-        simd4 = v.simd4;
-# else
-        std::memcpy(vec, v.vec, sizeof(int[4]));
-# endif
-        return *this;
-    }
-
-    inline simd4_t operator+(int i)
-    {
-        simd4_t r(*this);
-        r += i;
-        return r;
-    }
-    inline simd4_t operator+(const __vec4i_t v)
-    {
-        simd4_t r(*this);
-        r += v;
-        return r;
-    }
-    inline simd4_t operator+(const simd4_t<int>& v)
-    {
-        simd4_t r(*this);
-        r += v;
-        return r;
-    }
-
-    inline simd4_t operator-()
-    {
-        simd4_t r(0.0f);
-        r -= vec;
-        return r;
-    }
-    inline simd4_t operator-(int i)
-    {
-        simd4_t r(*this);
-        r -= i;
-        return r;
-    }
-    inline simd4_t operator-(const __vec4i_t v)
-    {
-        simd4_t r(*this);
-        r -= v;
-        return r;
-    }
-    inline simd4_t operator-(simd4_t<int>& v)
-    {
-        simd4_t r(*this);
-        r -= v;
-        return r;
-    }
-
-    inline simd4_t operator*(int i)
-    {
-        simd4_t r(*this);
-        r *= i;
-        return r;
-    }
-    inline simd4_t operator*(const __vec4i_t v)
-    {
-        simd4_t r(*this);
-        r *= v;
-        return r;
-    }
-    inline simd4_t operator*(simd4_t<int>& v)
-    {
-        simd4_t r(*this);
-        r *= v;
-        return r;
-    }
-
-    inline simd4_t operator/(int i)
-    {
-        simd4_t r(*this);
-        r /= i;
-        return r;
-    }
-    inline simd4_t operator/(const __vec4i_t v)
-    {
-        simd4_t r(*this);
-        r /= v;
-        return r;
-    }
-    inline simd4_t operator/(simd4_t<int>& v)
-    {
-        simd4_t r(*this);
-        r /= v;
-        return r;
-    }
-
-    inline simd4_t<int>& operator+=(int i)
-    {
-# ifdef __SSE2__
-        simd4 += i;
-# else
-        vec[0] += i;
-        vec[1] += i;
-        vec[2] += i;
-        vec[3] += i;
-# endif
-        return *this;
-    }
-    inline simd4_t<int>& operator+=(const __vec4i_t v)
-    {
-        simd4_t r(v);
-        *this += r;
-        return *this;
-    }
-    inline simd4_t<int>& operator+=(const simd4_t<int>& v)
-    {
-# ifdef __SSE2__
-        simd4 += v.simd4;
-# else
-        vec[0] += v[0];
-        vec[1] += v[1];
-        vec[2] += v[2];
-        vec[3] += v[3];
-# endif
-        return *this;
-    }
-
-    inline simd4_t<int>& operator-=(int i)
-    {
-# ifdef __SSE2__
-        simd4 -= i;
-# else
-        vec[0] -= i;
-        vec[1] -= i;
-        vec[2] -= i;
-        vec[3] -= i;
-# endif
-        return *this;
-    }
-    inline simd4_t<int>& operator-=(const __vec4i_t v)
-    {
-        simd4_t r(v);
-        *this -= r;
-        return *this;
-    }
-    inline simd4_t<int>& operator-=(const simd4_t<int>& v)
-    {
-# ifdef __SSE2__
-        simd4 -= v.simd4;
-# else
-        vec[0] -= v[0];
-        vec[1] -= v[1];
-        vec[2] -= v[2];
-        vec[3] -= v[3];
-# endif
-        return *this;
-    }
-
-    inline simd4_t<int>& operator *=(int i)
-    {
-# ifdef __SSE2__
-        simd4 *= i;
-# else
-        vec[0] *= i;
-        vec[1] *= i;
-        vec[2] *= i;
-        vec[3] *= i;
-# endif
-        return *this;
-    }
-    inline simd4_t<int>& operator*=(const __vec4i_t v)
-    {
-        simd4_t r(v);
-        *this *= r;
-        return *this;
-    }
-    inline simd4_t<int>& operator*=(const simd4_t<int>& v)
-    {
-# ifdef __SSE2__
-        simd4 *= v.simd4;
-# else
-        vec[0] *= v[0];
-        vec[1] *= v[1];
-        vec[2] *= v[2];
-        vec[3] *= v[3];
-# endif
-        return *this;
-    }
-
-    inline simd4_t<int>& operator/=(int i)
-    {
-# ifdef __SSE2__
-        i = (1/i);
-        simd4 *= i;
-# else
-        vec[0] /= i;
-        vec[1] /= i;
-        vec[2] /= i;
-        vec[3] /= i;
-# endif
-        return *this;
-    }
-    inline simd4_t<int>& operator/=(const __vec4i_t v)
-    {
-        simd4_t r(v);
-        *this /= r;
-        return *this;
-    }
-    inline simd4_t<int>& operator/=(const simd4_t<int>& v)
-    {
-# ifdef __SSE2__
-        simd4 /= v.simd4;
-# else
-        vec[0] /= v[0];
-        vec[1] /= v[1];
-        vec[2] /= v[2];
-        vec[3] /= v[3];
-# endif
-        return *this;
-    }
-
     inline operator const int*() const {
         return vec;
     }
@@ -917,7 +500,58 @@ public:
     inline operator int*() {
         return vec;
     }
+
+    inline simd4_t<int,N>& operator=(int i) {
+        simd4 = _mm_set1_epi32(i);
+        return *this;
+    }
+    inline simd4_t<int,N>& operator=(const __vec4i_t v) {
+        simd4 = _mm_loadu_si128((__m128i*)v);
+        return *this;
+    }
+    inline simd4_t<int,N>& operator=(const simd4_t<int,N>& v) {
+        simd4 = v.simd4;
+        return *this;
+    }
+
+    inline simd4_t<int,N>& operator+=(int i) {
+        simd4 += i;
+        return *this;
+    }
+    inline simd4_t<int,N>& operator+=(const simd4_t<int,N>& v) {
+        simd4 += v.simd4;
+        return *this;
+    }
+
+    inline simd4_t<int,N>& operator-=(int i) {
+        simd4 -= i;
+        return *this;
+    }
+    inline simd4_t<int,N>& operator-=(const simd4_t<int,N>& v) {
+        simd4 -= v.simd4;
+        return *this;
+    }
+
+    inline simd4_t<int,N>& operator *=(int i) {
+        simd4 *= i;
+        return *this;
+    }
+    inline simd4_t<int,N>& operator*=(const simd4_t<int,N>& v) {
+        simd4 *= v.simd4;
+        return *this;
+    }
+
+    inline simd4_t<int,N>& operator/=(int i) {
+        i = (1/i);
+        simd4 *= i;
+        return *this;
+    }
+    inline simd4_t<int,N>& operator/=(const simd4_t<int,N>& v) {
+        simd4 /= v.simd4;
+        return *this;
+    }
 };
+# endif
 
 
 #endif /* __SIMD_H__ */