module vec;
import BigRational;
private import std.math;
private import std.c.math;
private import std.string;

alias std.math math;
template VecMath(VALUE) {
  VALUE sqrt(VALUE x) {
    return math.sqrt(x);
  }
}
template VecMath(VALUE:BigRational.Rational) {
  VALUE sqrt(VALUE x) {
    VALUE ret;
    assert(0);
    return ret;
  }
}
template VecMath(VALUE:int) {
  int sqrt(int x) {
    return cast(int)(.01+math.sqrt(cast(float)x));//cast(int)floor((double)math.sqrt(cast(double)x)+.01);
  }
}

struct vec(VALUE, int tsize){ 
  enum {nvalues=(tsize>0?tsize:1)}
  VALUE [nvalues]f;
  VALUE* getAt(uint ii) {
    return ii<tsize?&f[ii]:&f[tsize-1];
  }  
  VALUE x() {
    return f[0];
  }
  VALUE y() {
    return f[1<tsize?1:0];
  }
  VALUE z() {
    return f[2<tsize?2:0];
  }
  VALUE w() {
    return f[3<tsize?3:0];
  }
  VALUE x(VALUE i) {
    f[0]=i;
    return i;
  }
  VALUE y(VALUE i) {
    if (tsize>1)
      f[1]=i;
    return i;
  }
  VALUE z(VALUE i) {
    if (tsize>2)
      f[2]=i;
    return i;
  }
  VALUE w(VALUE i) {
    if (tsize>3)
      f[3]=i;
    return i;
  }
  VALUE i() {
    return f[0];
  }
  VALUE j() {
    return f[1<tsize?1:0];
  }
  VALUE k() {
    return f[2<tsize?2:0];
  }
  VALUE i(VALUE i) {
    f[0]=i;
    return i;
  }
  VALUE j(VALUE i) {
    if (tsize>1)
      f[1]=i;
    return i;
  }
  VALUE k(VALUE i) {
    if (tsize>2)
      f[2]=i;
    return i;
  }

    
  static VALUE dot (vec t, vec oth) {
    VALUE ret=t.f[0]*oth.f[0];
    for (int ii=1;ii<tsize;++ii) {
      ret+=t.f[ii]*oth.f[ii];
    }
    return ret;
  }
  vec opNeg() {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=-f[ii];
    }
    return ret;
  }
  vec opCmp() {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=(BigRational.Zero!(VALUE).zero()==f[ii])?BigRational.Zero!(VALUE).one():BigRational.Zero!(VALUE).zero();
    }
    return ret;
  }
  /*
  .vec!(int,tsize) opAnd(vec b) {
    .vec!(int,tsize) ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=(f[ii]&&b[ii])?1:0;
    }
    return ret;   
  }
  .vec!(int,tsize) opOr(vec b) {
    .vec!(int,tsize) ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=(f[ii]||b[ii])?1:0;
    }
    return ret;   
  }
  .vec!(int,tsize) opXor(vec b) {
    .vec!(int,tsize) ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=((f[ii]||b[ii])&&!(f[ii]&&b[ii]))?1:0;
    }
    return ret;
  }
  
  vec opAndAssign(vec b) {
    for (int ii=0;ii<nvalues;++ii) {
      f[ii]=(f[ii]&&b[ii])?1:0;
    }
    return *this;   
  }
  vec opOrAssign(vec b) {
    for (int ii=0;ii<nvalues;++ii) {
      f[ii]=(f[ii]||b[ii])?1:0;
    }
    return *this;   
  }
  vec opXorAssign(vec b) {
    for (int ii=0;ii<nvalues;++ii) {
      f[ii]=((f[ii]||b[ii])&&!(f[ii]&&b[ii]))?1:0;
    }
    return *this;
  }
  */
  VALUE length() {
    return VecMath!(VALUE).sqrt(dot(*this,*this));
  }
  VALUE opIndex(uint ii, VALUE v) {
    return (*getAt(ii))=v;
  }
  VALUE opIndex(uint ii) {
    return ii<tsize?f[ii]:f[tsize-1];
  }
  vec opMulAssign(vec b) {
    for (int ii=0;ii<nvalues;++ii) { 
       this.f[ii]*=b.f[ii];
    }
    return *this;
  }
  vec opMul(vec b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=this.f[ii]*b.f[ii];
    }
    return ret;
  }
  vec opMulAssign(VALUE b) {
    for (int ii=0;ii<nvalues;++ii) { 
       this.f[ii]*=b;
    }
    return *this;
  }
  vec opMul(VALUE b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=this.f[ii]*b;
    }
    return ret;
  }

  vec opAddAssign(vec b) {
    for (int ii=0;ii<nvalues;++ii) { 
       this.f[ii]+=b.f[ii];
    }
    return *this;
  }
  vec opAdd(vec b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=this.f[ii]+b.f[ii];
    }
    return ret;
  }
  vec opAddAssign(VALUE b) {
    for (int ii=0;ii<nvalues;++ii) { 
       this.f[ii]+=b;
    }
    return *this;
  }
  vec opAdd(VALUE b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=this.f[ii]+b;
    }
    return ret;
  }


  vec opSubAssign(vec b) {
    for (int ii=0;ii<nvalues;++ii) { 
       this.f[ii]-=b.f[ii];
    }
    return *this;
  }
  vec opSub(vec b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=this.f[ii]-b.f[ii];
    }
    return ret;
  }
  vec opSubAssign(VALUE b) {
    for (int ii=0;ii<nvalues;++ii) { 
       this.f[ii]-=b;
    }
    return *this;
  }
  vec opSub(VALUE b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=this.f[ii]-b;
    }
    return ret;
  }

  vec opSub_r(vec b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=b.f[ii]-this.f[ii];
    }
    return ret;
  }
  vec opSub_r(VALUE b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=b-this.f[ii];
    }
    return ret;
  }

  vec opDivAssign(vec b) {
    for (int ii=0;ii<nvalues;++ii) { 
       this.f[ii]/=b.f[ii];
    }
    return *this;
  }
  vec opDiv(vec b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=this.f[ii]/b.f[ii];
    }
    return ret;
  }
  vec opDivAssign(VALUE b) {
    for (int ii=0;ii<nvalues;++ii) { 
       this.f[ii]/=b;
    }
    return *this;
  }
  vec opDiv(VALUE b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=this.f[ii]/b;
    }
    return ret;
  }

  vec opDiv_r(vec b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=b.f[ii]/this.f[ii];
    }
    return ret;
  }
  vec opDiv_r(VALUE b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=b/this.f[ii];
    }
    return ret;
  }


  vec opModAssign(vec b) {
    for (int ii=0;ii<nvalues;++ii) { 
       this.f[ii]%=b.f[ii];
    }
    return *this;
  }
  vec min(vec b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) {
       ret[ii]=f[ii]<b[ii]?f[ii]:b[ii];
    }
    return ret;
  }
  vec max(vec b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) {
       ret[ii]=f[ii]<b[ii]?b[ii]:f[ii];
    }
    return ret;
  }
  vec opMod(vec b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=this.f[ii]%b.f[ii];
    }
    return ret;
  }
  vec opModAssign(VALUE b) {
    for (int ii=0;ii<nvalues;++ii) { 
       this.f[ii]%=b;
    }
    return *this;
  }
  vec opMod(VALUE b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=this.f[ii]%b;
    }
    return ret;
  }


  vec opMod_r(vec b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=b.f[ii]%this.f[ii];
    }
    return ret;
  }
  vec opMod_r(VALUE b) {
    vec ret;
    for (int ii=0;ii<nvalues;++ii) { 
       ret.f[ii]=b%this.f[ii];
    }
    return ret;
  }

  .vec!(int,tsize) opEq(vec inp) {
    .vec!(int,tsize) ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=(f[ii]==inp.f[ii]);
    }
    return ret;
  }

  .vec!(int,tsize) opLess(vec inp) {
    .vec!(int,tsize) ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=(f[ii]<inp.f[ii]);
    }
    return ret;
  }
  .vec!(int,tsize) opLess(VALUE inp) {
    .vec!(int,tsize) ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=(f[ii]<inp);
    }
    return ret;
  }
  .vec!(int,tsize) opLEqual(vec inp) {
    .vec!(int,tsize) ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=(f[ii]<inp.f[ii]);
    }
    return ret;
  }
  .vec!(int,tsize) opLEqual(VALUE inp) {
    .vec!(int,tsize) ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=(f[ii]<inp);
    }
    return ret;
  }
  .vec!(int,tsize) opGreater(vec inp) {
    .vec!(int,tsize) ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=(f[ii]>inp.f[ii]);
    }
    return ret;
  }
  .vec!(int,tsize) opGreater(VALUE inp) {
    .vec!(int,tsize) ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=(f[ii]>inp);
    }
    return ret;
  }

  .vec!(int,tsize) opGEqual(vec inp) {
    .vec!(int,tsize) ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=(f[ii]>=inp.f[ii]);
    }
    return ret;
  }
  .vec!(int,tsize) opGEqual(VALUE inp) {
    .vec!(int,tsize) ret;
    for (int ii=0;ii<nvalues;++ii) {
      ret.f[ii]=(f[ii]>=inp);
    }
    return ret;
  }

  //disallow!
  int opEqual(vec o){
    for (int ii=0;ii<nvalues;++ii) {
       if (o.f[ii]!=f[ii]) return 0;
    }
    return 1;
  }
  int opCmp(vec o){
    if (nvalues==1) {
      return this.opCmp(o.f[0]);
    }
    VALUE thus=dot(*this,*this);
    VALUE oth = dot(o,o);
    if (thus<oth) return -1;
    if (thus>oth) return 1;    
    return 0;
  }
  int opEqual(VALUE o){
    for (int ii=0;ii<nvalues;++ii) {
       if (o!=f[ii]) return 0;
    }
    return 1;
  }
  int opCmp(VALUE o){
    if (nvalues==1) {
      if (f[0]<o) return -1;
      if (f[0]>o) return 1;
      return 0;
    }
    VALUE thus=dot(*this,*this);
    VALUE oth = o>=BigRational.Zero!(VALUE).zero()?o:-o;
    if (thus<oth) return -1;
    if (thus>oth) return 1;    
    return 0;
  }
  int any() {
    for (int ii=1;ii<nvalues;++ii) {
      if (f[ii]!=BigRational.Zero!(VALUE).zero()) return 1;
    }
    return 0;
  }
  int all() {
    for (int ii=1;ii<nvalues;++ii) {
      if (f[ii]==BigRational.Zero!(VALUE).zero()) return 0;
    }
    return 1;
  }
  vec print () {
    printf ("[ ");
    for (int ii=0;ii<nvalues;++ii) { 
       printf("%lf ",cast(double)f[ii]);
       
    }
    printf ("]");
    return *this;
  }
  static vec opCall(VALUE a, VALUE b, VALUE c, VALUE d) {
     vec ret;
     ret[3]=d;
     ret[2]=c;
     ret[1]=b;
     ret[0]=a;
     return ret;
  }
  static vec opCall(VALUE a, VALUE b, VALUE c) {
     vec ret;
     ret[2]=c;
     ret[1]=b;
     ret[0]=a;
     return ret;
  }
  static vec opCall(VALUE a, VALUE b) {
     vec ret;
     ret[1]=b;
     ret[0]=a;
     return ret;
  }
  static vec opCall(vec v) {
     vec ret;
     for (int ii=nvalues-1;ii>=0;--ii) {
       ret[ii]=v[ii];
     }
     return ret;
  }
  static vec opCall(VALUE []v) {
     vec ret;
     for (int ii=v.length-1;ii>=0;--ii) {
       ret[ii]=v[ii];
     }
     return ret;
  }
  static vec opCall(VALUE a) {
     vec ret;
     for (int ii=0;ii<ret.nvalues;++ii) {
       ret[ii]=a;
     }
     return ret;
  }

  .vec!(VALUE,4) swizzle(int a, int b, int c, int d) {
     .vec!(VALUE,4) ret;
     ret.f[0]=*getAt(a);
     ret.f[1]=*getAt(b);
     ret.f[2]=*getAt(c);
     ret.f[3]=*getAt(d);
     return ret;
  }
  .vec!(VALUE,3) swizzle(int a, int b, int c) {
     .vec!(VALUE,3) ret;
     ret.f[0]=*getAt(a);
     ret.f[1]=*getAt(b);
     ret.f[2]=*getAt(c);
     return ret;
  }
  .vec!(VALUE,2) swizzle(int a, int b) {
     .vec!(VALUE,2) ret;
     ret.f[0]=*getAt(a);
     ret.f[1]=*getAt(b);
     return ret;
  }
  .vec!(VALUE,1) swizzle(int a) {
     .vec!(VALUE,1) ret;
     ret.f[0]=*getAt(a);
     return ret;
  }

  private:
  class Disallowed{};  
}

alias vec!(BigRational.Rational,4) rat4;
alias vec!(BigRational.Rational,3) rat3;
alias vec!(BigRational.Rational,2) rat2;

/*
alias vec!(int,4) rat4;
alias vec!(int,3) rat3;
alias vec!(int,2) rat2;
*/