From e67e7ecdfc2cc3eeba896aff365a11d6c7641aac Mon Sep 17 00:00:00 2001
From: Marc Glisse <Marc.Glisse@inria.fr>
Date: Mon, 23 Apr 2012 13:09:15 +0000
Subject: [PATCH] Untested SSE Vector_2

---
 NewKernel_d/include/CGAL/Vector/avx4.h |   2 +
 NewKernel_d/include/CGAL/Vector/sse2.h | 109 +++++++++++++++++++++++++
 2 files changed, 111 insertions(+)
 create mode 100644 NewKernel_d/include/CGAL/Vector/sse2.h

diff --git a/NewKernel_d/include/CGAL/Vector/avx4.h b/NewKernel_d/include/CGAL/Vector/avx4.h
index a11bbbe8538..4dbc6122605 100644
--- a/NewKernel_d/include/CGAL/Vector/avx4.h
+++ b/NewKernel_d/include/CGAL/Vector/avx4.h
@@ -78,6 +78,7 @@ namespace CGAL {
 
       struct Values_divide {
 	Vector operator()(double h,double a,double b,double c,double d) const {
+	  // {a,b,c,d}/{h,h,h,h} should be roughly the same
 	  Vector r = { a/h, b/h, c/h, d/h };
 	  return r;
 	}
@@ -168,6 +169,7 @@ namespace CGAL {
       return z+p[2];
     }
     public:
+    // Note: without AVX2, is it faster than the scalar computation?
     static double
       determinant_of_vectors_omit_last(Vector a, Vector b, Vector c) {
       __m256d x=a*avx3_right(b)-avx3_right(a)*b;
diff --git a/NewKernel_d/include/CGAL/Vector/sse2.h b/NewKernel_d/include/CGAL/Vector/sse2.h
new file mode 100644
index 00000000000..1b7b515157c
--- /dev/null
+++ b/NewKernel_d/include/CGAL/Vector/sse2.h
@@ -0,0 +1,109 @@
+#ifndef CGAL_VECTOR_SSE2_H
+#define CGAL_VECTOR_SSE2_H
+
+// Requires at least my patches for [], so maybe 4.8. Check intel, clang and MS.
+#ifndef __SSE2__ || (__GNUC_MAJOR__*100+__GNUC_MINOR__ < 408)
+#error Only supported on gcc 4.8+
+#endif
+#include <x86intrin.h>
+
+#include <CGAL/Dimension.h>
+#include <CGAL/enum.h> // CGAL::Sign
+#include <CGAL/number_utils.h> // CGAL::sign
+
+
+
+namespace CGAL {
+
+  struct Sse_vector_2 {
+    typedef double NT;
+    typedef Dimension_tag<2> Dimension;
+    typedef Dimension_tag<2> Max_dimension;
+    // No Rebind_dimension, this is a building block
+    template<class,bool=true> struct Property : boost::false_type {};
+    template<bool b> struct Property<Has_vector_plus_minus_tag,b>
+      : boost::true_type {};
+    /* MAYBE?
+       template<bool b> struct Property<Has_vector_scalar_ops_tag,b>
+       : boost::true_type {};
+       */
+    template<bool b> struct Property<Has_determinant_of_vectors_tag,b>
+      : boost::true_type {};
+
+    typedef __m128d Vector;
+    struct Construct_vector {
+      struct Dimension {
+	// Initialize with NaN?
+	Vector operator()(unsigned d) const {
+	  CGAL_assertion(d==2);
+	  return Vector();
+	}
+      };
+
+      struct Iterator {
+	template<typename Iter>
+	  Vector operator()(unsigned d,Iter const& f,Iter const& e) const {
+	    CGAL_assertion(d==2);
+	    double x0 = *f;
+	    double x1 = *++f;
+	    CGAL_assertion(++f==e);
+	    Vector a = { x0, x1 };
+	    return a;
+	  }
+      };
+
+      struct Iterator_and_last {
+	template<typename Iter,typename T>
+	  Vector operator()(unsigned d,Iter const& f,Iter const& e,double t) const {
+	    CGAL_assertion(d==2);
+	    Vector a = { *f, t };
+	    CGAL_assertion(++f==e);
+	    return a;
+	  }
+      };
+
+      struct Values {
+	  Vector operator()(double a,double b) const {
+	    Vector r = { a, b };
+	    return r;
+	  }
+      };
+
+      struct Values_divide {
+	Vector operator()(double h,double a,double b) const {
+	  // {a,b}/{h,h} is probably slower
+	  Vector r = { a/h, b/h };
+	  return r;
+	}
+      };
+    };
+
+    typedef double const* Vector_const_iterator;
+    static inline Vector_const_iterator vector_begin(Vector const&a){
+      return static_cast<Vector_const_iterator>&a;
+    }
+    static inline Vector_const_iterator vector_end(Vector const&a){
+      return static_cast<Vector_const_iterator>&a+2;
+    }
+    static inline unsigned size_of_vector(Vector){
+      return 2;
+    }
+
+    static double determinant_of_vectors(Vector a, Vector b) {
+      __m128d c = _mm_shuffle_pd (b, b, 1); // b1, b0
+      __m128d d = a * c; // a0*b1, a1*b0
+#ifdef __SSE3__
+      __m128d e = _mm_hsub_pd (d, d);
+      return e[0];
+#else
+      return d[0]-d[1];
+#endif
+    }
+    static CGAL::Sign sign_of_determinant_of_vectors(Vector a, Vector b) {
+      return CGAL::sign(determinant_of_vectors(a,b));
+    }
+
+  };
+
+}
+#endif