#ifndef L1Algo_PSEUDO_F64vec2_H
#define L1Algo_PSEUDO_F64vec2_H

#include <iostream>
using namespace std; // iklm

#include <cmath>
#include "vec_arithmetic.h"

/**********************************
 *
 *   Vector of four doubles
 *
 **********************************/
class F64vec2 
{
 public:

  double v[2];

  double & operator[]( int i ){ return ((double*)&v)[i]; }
  double   operator[]( int i ) const { return ((double*)&v)[i]; }

  F64vec2( ){}
  F64vec2( const F64vec2 &a ) { 
    v[0] = a.v[0]; 
    v[1] = a.v[1]; 
  }
  F64vec2( const double &a ) { 
    v[0] = a; 
    v[1] = a;
  }

  F64vec2( const double &f0, const double &f1, const double &f2, const double &f3 ){ 
    v[0] = f0; 
    v[1] = f1;
  }

  friend double min( double x, double y ){ return x<y ?x :y; }
  friend double max( double x, double y ){ return x<y ?y :x; }
  friend double asgnb( double x, double y ){ return y>=0 ?fabs(x) :-fabs(x); }
  friend double rsqrt( double x ){ return 1./sqrt(x); }
  friend double rcp( double x ){ return 1./x; }
  friend double sgn( double x ){ return x>=0 ?1 :-1; }

#define _f2(A,B,F) F64vec2 z; z.v[0] = F(A.v[0],B.v[0]); z.v[1] = F(A.v[1],B.v[1]); \
  return z;
#define _f1(A,F)   F64vec2 z; z.v[0] = F(A.v[0]);        z.v[1] = F(A.v[1]); \
   return z;
#define _op(A,B,O) F64vec2 z; z.v[0] = A.v[0] O B.v[0];  z.v[1] = A.v[1] O B.v[1]; \
  return z; 
#define _op2(B,O) v[0] O B.v[0];  v[1] O B.v[1]; \
  return *this;
  #define _op3(A,B,O) F64vec2 z;                                          \
  z.v[0] = (A.v[0] O B.v[0]) ? 1. : 0.;                                 \
  z.v[1] = (A.v[1] O B.v[1]) ? 1. : 0.;                                 \
  z.v[2] = (A.v[2] O B.v[2]) ? 1. : 0.;                                 \
  z.v[3] = (A.v[3] O B.v[3]) ? 1. : 0.; return z;

//   bool operator bool(const F64vec2 &a){if }

  friend F64vec2 operator-( const F64vec2 &a ){ return F64vec2(0) - a; }
  friend F64vec2 operator+( const F64vec2 &a ){ return a; } 
  friend F64vec2 operator!( const F64vec2 &a ){ return F64vec2(1) - a; } // CHECKME
  
  /* Arithmetic Operators */
  friend F64vec2 operator +(const F64vec2 &a, const F64vec2 &b){ _op(a,b,+) }
  friend F64vec2 operator -(const F64vec2 &a, const F64vec2 &b){ _op(a,b,-) }
  friend F64vec2 operator *(const F64vec2 &a, const F64vec2 &b){ _op(a,b,*) }
  friend F64vec2 operator /(const F64vec2 &a, const F64vec2 &b){ _op(a,b,/) }
  friend F64vec2 operator <(const F64vec2 &a, const F64vec2 &b){ _op3(a,b,<) }
  friend F64vec2 operator >(const F64vec2 &a, const F64vec2 &b){ _op3(a,b,>) }
  friend F64vec2 operator>=(const F64vec2 &a, const F64vec2 &b){ _op3(a,b,>=) }
  friend F64vec2 operator &(const F64vec2 &a, const F64vec2 &b){ return a*b; }


//   friend bool operator ==(const F64vec2 &a, const F64vec2 &b){ _op(a,b,==) }
  F64vec2& operator +=(const F64vec2 &a) { _op2(a,+=) }
	F64vec2& operator -=(const F64vec2 &a) { _op2(a,-=) }   
	F64vec2& operator *=(const F64vec2 &a) { _op2(a,*=) } 
	F64vec2& operator /=(const F64vec2 &a) { _op2(a,/=) }
// 	F64vec2& operator &=(const F64vec2 &a) { _op2(a,&=) }
// 	F64vec2& operator |=(const F64vec2 &a) { _op2(a,|=) }
// 	F64vec2& operator ^=(const F64vec2 &a) { _op2(a,^=) }


  /* Functions */
  friend F64vec2 min( const F64vec2 &a, const F64vec2 &b ){ _f2(a,b,min) }
  friend F64vec2 max( const F64vec2 &a, const F64vec2 &b ){ _f2(a,b,max) }
//   friend F64vec2 asgnb( const F64vec2 &a, const F64vec2 &b ){ _f2(a,b,asgnb) }
  friend F64vec2 sqrt ( const F64vec2 &a ){ _f1(a,sqrt) }
//   friend F64vec2 rsqrt( const F64vec2 &a ){ _f1(a,rsqrt) }
//   friend F64vec2 rcp  ( const F64vec2 &a ){ _f1(a,rcp) }
  friend F64vec2 fabs (const F64vec2 &a) { _f1(a,fabs) }
//   friend F64vec2 sgn  ( const F64vec2 &a ){ _f1(a,sgn) }
  friend F64vec2 exp( const F64vec2 &a ){ _f1(a,exp) } 
  friend F64vec2 log( const F64vec2 &a ){ _f1(a,log) } 
  friend F64vec2 sin( const F64vec2 &a ){ _f1(a,sin) } 
  friend F64vec2 cos( const F64vec2 &a ){ _f1(a,cos) } 
#undef _f1
#undef _f2
#undef _op

  friend F64vec2 rsqrt(const F64vec2 &a) {
    return 1./sqrt(a);
  }
  friend F64vec2 rcp(const F64vec2 &a) {
    return 1/a;
  }

  friend F64vec2 sgn(const F64vec2 &a){
    return 2*(a>=0) - 1;
  }
  // friend F64vec2 asgnb(const F64vec2 &a, const F64vec2 &b ){ 
  //   return 0; 
  // }

  /* Define all operators for consistensy */
  
  vec_arithmetic(F64vec2,double);

  friend ostream & operator<<(ostream &strm, const F64vec2 &a ){
    strm<<a[0]<<" "<<a[1];
    return strm;
  }

  friend istream & operator>>(istream &strm, F64vec2 &a ){
    double tmp;
    strm>>tmp;
    a = tmp;
    return strm;
  }

} __attribute__ ((aligned(16)));;

typedef F64vec2 fvec;
typedef double fscal;
const int fvecLen = 2;
//#define fvec_true  _f32vec4_true
//#define fvec_false _f32vec4_false
#define _fvecalignment  


#endif