vectormath_ng.rb

DEG2RAD=Math::PI/180
RAD2DEG=180/Math::PI

class Vector2
  attr_accessor(:x,:y)
  def initialize(x=0.0,y=0.0) @x=x;@y=y; end
  def serialize; {x:@x,y:@y} end
  def inspect; serialize.to_s end
  def to_a; serialize.values end
  def to_s; serialize.to_s end

  def add_f(r) new(@x+r,@y+r) end
  def add_f!(r) @x+=r;@y+=r;self end
  def sub_f(r) new(@x-r,@y-r) end
  def sub_f!(r) @x-=r;@y-=r;self end
  def mul_f(r) new(@x*r,@y*r) end
  def mul_f!(r) @x*=r;@y*=r;self end
  def add_v(r) new(@x+r.x,@y+r.y) end
  def add_v!(r) @x+=r.x;@y+=r.y;self end
  def sub_v(r) new(@x-r.x,@y-r.y) end
  def sub_v!(r) @x-=r.x;@y-=r.y;self end
  def mul_v(r) new(@x*r.x,@y*r.y) end
  def mul_v!(r) @x*=r.x;@y*=r.y;self end
  def mul_m(r) new(@x*r.xx+@y*r.yx+r.wx,@x*r.xy+@y*r.yy+r.wy) end
  def mul_m!(r) @x=@x*r.xx+@y*r.yx+r.wx;@y=@x*r.xy+@y*r.yy+r.wy;self end

  def dot(r) @x*r.x+@y*r.y end
  def cross(r) @x*r.y-@y*r.x end
  def length; Math.sqrt(@x*@x+@y*@y) end
  def length2; @x*@x+@y*@y end
  def negate; new(-@x,-@y) end
  def negate!; @x=-@x;@y=-@y;self end
  def normalize; l=1/Math.sqrt(@x*@x+@y*@y);new(@x*l,@y*l) end
  def normalize!; l=1/Math.sqrt(@x*@x+@y*@y);@x*=l;@y*=l;self end

  def angle_to(v2) Math.atan2(v2.y-@y,v2.x-@x) end
  def perpendicular; new(-@y, @x) end
  def perpendicular!; @y=-@y;@x=-@x;self end

  def self.lerp(l,r,t) new(l.x+(r.x-l.x)*t,l.y+(r.y-l.y)*t); end
  def self.catmull_rom(a,b,c,d,t) t2=t*t;t3=t2*t;new(0.5*((2.0*b.x)+(-a.x+c.x)*t+(2.0*a.x-5.0*b.x+4.0*c.x-d.x)*t2+(-a.x+3.0*b.x-3.0*c.x+d.x)*t3),0.5*((2.0*b.y)+(-a.y+c.y)*t+(2.0*a.y-5.0*b.y+4.0*c.y-d.y)*t2+(-a.y+3.0*b.y-3.0*c.y+d.y)*t3)) end
  def self.reflect(i,n) d=i.x*n.x+i.y*n.y;new(i.x-2.0*d*n.x,i.y-2.0*d*n.y) end
  def self.refract(i,n,r) d=i.x*n.x+i.y*n.y;k=1-r*r*(1-d*d);k<0.0?new:new(i.x*r-n.x*(r*d+Math.sqrt(k)),i.y*r-n.y*(r*d+Math.sqrt(k))) end
  def self.barycentric(a,b,c,f,g) new(a.x+(b.x-a.x)*f+(c.x-a.x)*g,a.y+(b.y-a.y)*f+(c.y-a.y)*g); end
  def self.min(a,b) new((a.x<b.x)?(a.x):(b.x),(a.y<b.y)?(a.y):(b.y)) end
  def self.max(a,b) new((a.x>b.x)?(a.x):(b.x),(a.y>b.y)?(a.y):(b.y)) end
end

class Vector3
  attr_accessor(:x,:y,:z)
  def initialize(x=0.0,y=0.0,z=0.0) @x=x;@y=y;@z=z end
  def serialize; {x:@x,y:@y,z:@z} end
  def inspect; serialize.to_s end
  def to_a; serialize.values end
  def to_s; serialize.to_s end

  def add_f(r) new(@x+r,@y+r,@z+r) end
  def add_f!(r) @x+=r;@y+=r;@z+=r;self end
  def sub_f(r) new(@x-r,@y-r,@z-r) end
  def sub_f!(r) @x-=r;@y-=r;@z-=r;self end
  def mul_f(r) new(@x*r,@y*r,@z*r) end
  def mul_f!(r) @x*=r;@y*=r;@z*=r;self end
  def add_v(r) new(@x+r.x,@y+r.y,@z+r.z) end
  def add_v!(r) @x+=r.x;@y+=r.y;@z+=r.z;self end
  def sub_v(r) new(@x-r.x,@y-r.y,@z-r.z) end
  def sub_v!(r) @x-=r.x;@y-=r.y;@z-=r.z;self end
  def mul_v(r) new(@x*r.x,@y*r.y,@z*r.z) end
  def mul_v!(r) @x*=r.x;@y*=r.y;@z*=r.z;self end
  def mul_m(r) new(@x*r.xx+@y*r.yx+@z*r.zx+r.wx,@x*r.xy+@y*r.yy+@z*r.zy+r.wy,@x*r.xz+@y*r.yz+@z*r.zz+r.wz) end
  def mul_m!(r) @x=@x*r.xx+@y*r.yx+@z*r.zx+r.wx;@y=@x*r.xy+@y*r.yy+@z*r.zy+r.wy;@z=@x*r.xz+@y*r.yz+@z*r.zz+r.wz;self end

  def dot(r) @x*r.x+@y*r.y+@z*r.z end
  def cross(r) new(@y*r.z-@z*r.y,@z*r.x-@x*r.z,@x*r.y-@y*r.x) end
  def cross!(r) x=@y*r.z-@z*r.y;y=@z*r.x-@x*r.z;z=@x*r.y-@y*r.x;@x=x;@y=y;@z=z;self end
  def length; Math.sqrt(@x*@x+@y*@y+@z*@z) end
  def length2; @x*@x+@y*@y+@z*@z end
  def negate; new(-@x,-@y,-@z) end
  def negate!; @x=-@x;@y=-@y;@z=-@z;self end
  def normalize; l=1.0/Math.sqrt(@x*@x+@y*@y+@z*@z);new(@x*l,@y*l,@z*l) end
  def normalize!; l=1.0/Math.sqrt(@x*@x+@y*@y+@z*@z);@x*=l;@y*=l;@z*=l;self end

  def angle_to(v3) Math.acos((@x*v3.x+@y*v3.y+@z*v3.z)/(Math.sqrt(@x*@x+@y*@y+@z*@z)*Math.sqrt(v3.x*v3.x+v3.y*v3.y+v3.z*v3.z))) end
  def perpendicular; @x.abs>@y.abs ?new(-@z,0,@x):new(0,@z,-@y) end
  def perpendicular!; @x.abs>@y.abs ?(temp=@x;@x=-@z;@z=temp):temp=@y;@y=@z;@z=-temp; self end

  def self.lerp(l,r,t) new(l.x+(r.x-l.x)*t,l.y+(r.y-l.y)*t,l.z+(r.z-l.z)*t) end
  def self.catmull_rom(a,b,c,d,t) t2=t*t;t3=t2*t;new(0.5*((2.0*b.x)+(-a.x+c.x)*t+(2.0*a.x-5.0*b.x+4.0*c.x-d.x)*t2+(-a.x+3.0*b.x-3.0*c.x+d.x)*t3),0.5*((2.0*b.y)+(-a.y+c.y)*t+(2.0*a.y-5.0*b.y+4.0*c.y-d.y)*t2+(-a.y+3.0*b.y-3.0*c.y+d.y)*t3),0.5*((2.0*b.z)+(-a.z+c.z)*t+(2.0*a.z-5.0*b.z+4.0*c.z-d.z)*t2+(-a.z+3.0*b.z-3.0*c.z+d.z)*t3)) end
  def self.refract(i,n,r) d=n.x*i.x+n.y*i.y+n.z*i.z;k=1-r*r*(1-d*d);t=r*d+Math.sqrt(k);new(i.x*r-n.x*t,i.y*r-n.y*t,i.z*r-n.z*t) end
  def self.reflect(i,n) d=2*(i.x*n.x+i.y*n.y+i.z*n.z);new(i.x-d*n.x,i.y-d*n.y,i.z-d*n.z) end
  def self.barycentric(a,b,c,f,g) new(a.x+(b.x-a.x)*f+(c.x-a.x)*g,a.y+(b.y-a.y)*f+(c.y-a.y)*g,a.z+(b.z-a.z)*f+(c.z-a.z)*g) end
  def self.min(a,b) new((a.x<b.x)?(a.x):(b.x),(a.y<b.y)?(a.y):(b.y),(a.y<b.z)?(a.z):(b.z)) end
  def self.max(a,b) new((a.x>b.x)?(a.x):(b.x),(a.y>b.y)?(a.y):(b.y),(a.z>b.z)?(a.z):(b.z)) end
end

class Vector4
  attr_accessor(:x,:y,:z,:w)
  def initialize(x=0,y=0,z=0,w=1) @x=x;@y=y;@z=z;@w=w end
  def serialize; {x:@x,y:@y,z:@z,w:@w} end
  def inspect; serialize.to_s end
  def to_a; serialize.values end
  def to_s; serialize.to_s end

  def add_f(r) new(@x+r,@y+r,@z+r,@w+r) end
  def add_f!(r) @x+=r;@y+=r;@z+=r;@w+=r;self end
  def sub_f(r) new(@x-r,@y-r,@z-r,@w-r) end
  def sub_f!(r) @x-=r;@y-=r;@z-=r;@w-=r;self end
  def mul_f(r) new(@x*r,@y*r,@z*r,@w*r) end
  def mul_f!(r) @x*=r;@y*=r;@z*=r;@w*=r;self end
  def add_v(r) new(@x+r.x,@y+r.y,@z+r.z,@w+r.w) end
  def add_v!(r) @x+=r.x;@y+=r.y;@z+=r.z;@w+=r.w;self end
  def sub_v(r) new(@x-r.x,@y-r.y,@z-r.z,@w-r.w) end
  def sub_v!(r) @x-=r.x;@y-=r.y;@z-=r.z;@w-=r.w;self end
  def mul_v(r) new(@x*r.x,@y*r.y,@z*r.z,@w*r.w) end
  def mul_v!(r) @x*=r.x;@y*=r.y;@z*=r.z;@w*=r.w;self end
  def mul_m(r) new(@x*r.xx+@y*r.yx+@z*r.zx+@w*r.wx,@x*r.xy+@y*r.yy+ @z*r.zy+@w*r.wy,@x*r.xz+@y*r.yz+@z*r.zz+@w*r.wz,@x*r.xw+@y*r.yw+@z*r.zw+@w*r.ww) end
  def mul_m!(r) @x=@x*r.xx+@y*r.yx+@z*r.zx+@w*r.wx;@y=@x*r.xy+@y*r.yy+ @z*r.zy+@w*r.wy;@z=@x*r.xz+@y*r.yz+@z*r.zz+@w*r.wz;@w=@x*r.xw+@y*r.yw+@z*r.zw+@w*r.ww;self end

  def dot(r) @x*r.x+@y*r.y+@z*r.z+@w*r.w end
  def cross(r) new(@y*r.z-@z*r.y,@z*r.x-@x*r.z,@x*r.y-@y*r.x,0) end
  def cross!(r) @x=@y*r.z-@z*r.y;@y=@z*r.x-@x*r.z;@z=@x*r.y-@y*r.x;@w=0;self end
  def length; Math.sqrt(@x*@x+@y*@y+@z*@z+@w*@w) end
  def length2; @x*@x+@y*@y+@z*@z+@w*@w end
  def negate; new(-@x,-@y,-@z,-@w) end
  def negate!; @x=-@x;@y=-@y;@z=-@z;@w=-@w;self end
  def normalize; l=1/Math.sqrt(@x*@x+@y*@y+@z*@z+@w*@w);Vector4.new(@x*l,@y*l,@z*l,@w*l) end
  def normalize!; l=1/Math.sqrt(@x*@x+@y*@y+@z*@z+@w*@w);@x*=l;@y*=l;@z*=l;@w*=l;self end

  def angle_to(v4) Math.acos((@x*v4.x+@y*v4.y+@z*v4.z+@w*v4.w)/(Math.sqrt(@x*@x+@y*@y+@z*@z+@w*@w)*Math.sqrt(v4.x*v4.x+v4.y*v4.y+v4.z*v4.z+v4.w*v4.w))) end
  def perpendicular; @x.abs>@y.abs ?new(-@z,0,@x,0):new(0,@z,-@y,0) end
  def perpendicular!; @x.abs>@y.abs ?(t=@x;@x=-@z;@z=t;@y=0;@w=0):(t=@y;@y=@z;@z=-t;@x=0;@w=0); self end

  def self.lerp(l,r,t) new(l.x+(r.x-l.x)*t,l.y+(r.y-l.y)*t,l.z+(r.z-l.z)*t,l.w+(r.w-l.w)*t) end
  def self.catmull_rom(a,b,c,d,t) t2=t*t;t3=t2*t;new(0.5*((2.0*b.x)+(-a.x+c.x)*t+(2.0*a.x-5.0*b.x+4.0*c.x-d.x)*t2+(-a.x+3.0*b.x-3.0*c.x+d.x)*t3),0.5*((2.0*b.y)+(-a.y+c.y)*t+(2.0*a.y-5.0*b.y+4.0*c.y-d.y)*t2+(-a.y+3.0*b.y-3.0*c.y+d.y)*t3),0.5*((2.0*b.z)+(-a.z+c.z)*t+(2.0*a.z-5.0*b.z+4.0*c.z-d.z)*t2+(-a.z+3.0*b.z-3.0*c.z+d.z)*t3),0.5*((2.0*b.w)+(-a.w+c.w)*t+(2.0*a.w-5.0*b.w+4.0*c.w-d.w)*t2+(-a.w+3.0*b.w-3.0*c.w+d.w)*t3)) end
  def self.refract(i,n,r); c=i.dot(n);c=c<0 ? -c:c;c=c>1 ?1:c;c=-c;e=1.0;f=r;f=f/e;k=1.0-(f*f)*(1.0-c*c);k<0.0 ?new:new(i.x*f+n.x*(f*c-Math.sqrt(k)),i.y*f+n.y*(f*c-Math.sqrt(k)),i.z*f+n.z*(f*c-Math.sqrt(k)),i.w*f+n.w*(f*c-Math.sqrt(k))) end
  def self.reflect(i,n) new(i.x-2.0*(i.x*n.x),i.y-2.0*(i.y*n.y),i.z-2.0*(i.z*n.z),i.w-2.0*(i.w*n.w)) end
  def self.barycentric(a,b,c,f,g) new(a.x+(b.x-a.x)*f+(c.x-a.x)*g,a.y+(b.y-a.y)*f+(c.y-a.y)*g,a.z+(b.z-a.z)*f+(c.z-a.z)*g,a.w+(b.w-a.w)*f+(c.w-a.w)*g) end
  def self.min(a,b) new((a.x<b.x)?(a.x):(b.x),(a.y<b.y)?(a.y):(b.y),(a.x<b.z)?(a.z):(b.z),(a.w<b.w)?(a.w):(b.w)) end
  def self.max(a,b) new((a.x>b.x)?(a.x):(b.x),(a.y>b.y)?(a.y):(b.y),(a.z>b.z)?(a.z):(b.z),(a.w>b.w)?(a.w):(b.w)) end
end

class Quaternion
  attr_accessor(:x,:y,:z,:w)
  def initialize(x=0,y=0,z=0,w=1) @x=x;@y=y;@z=z;@w=w end
  def serialize; {x:@x,y:@y,z:@z,w:@w} end
  def inspect; serialize.to_s end
  def to_a; serialize.values end
  def to_s; serialize.to_s end

  def mul_q(r)x=r.x;y=r.y;z=r.z;w=r.w;Quaternion.new(w*@x+x*@w+y*@z-z*@y,w*@y-x*@z+y*@w+z*@x,w*@z+x*@y-y*@x+z*@w,w*@w-x*@x-y*@y-z*@z) end
  def mul_q!(r) x=r.x;y=r.y;z=r.z;w=r.w;_x=w*@x+x*@w+y*@z-z*@y;_y=w*@y-x*@z+y*@w+z*@x;_z=w*@z+x*@y-y*@x+z*@w;_w=w*@w-x*@x-y*@y-z*@z;@x=_x;@y=_y;@z=_z;@w=_w;self end

  def dot(r) @x*r.x+@y*r.y+@z*r.z+@w*r.w end

  def length; Math.sqrt(@x*@x+@y*@y+@z*@z+@w*@w) end
  def length2; @x*@x+@y*@y+@z*@z+@w*@w end
  def negate; Quaternion.new(-@x,-@y,-@z,-@w) end
  def negate!; @x=-@x;@y=-@y;@z=-@z;@w=-@w;self end
  def normalize; l=1/Math.sqrt(@x*@x+@y*@y+@z*@z+@w*@w);Quaternion.new(@x*l,@y*l,@z*l,@w*l) end
  def normalize!; l=1/Math.sqrt(@x*@x+@y*@y+@z*@z+@w*@w);@x*=l;@y*=l;@z*=l;@w*=l;self end

  def rotate_x!(r) s=Math.sin(r*0.5);c=Math.cos(r*0.5);_x=c*@x+s*@w;_y=c*@y-s*@z;_z=c*@z+s*@y;_w=c*@w-s*@x;@x=_x;@y=_y;@z=_z;@w=_w;self end
  def rotate_y!(r) s=Math.sin(r*0.5);c=Math.cos(r*0.5);_x=c*@x+s*@z;_y=c*@y-s*@w;_z=c*@z-s*@x;_w=c*@w+s*@y;@x=_x;@y=_y;@z=_z;@w=_w;self end
  def rotate_z!(r) s=Math.sin(r*0.5);c=Math.cos(r*0.5);_x=c*@x-s*@y;_y=c*@y+s*@x;_z=c*@z-s*@w;_w=c*@w+s*@z;@x=_x;@y=_y;@z=_z;@w=_w;self end

  def self.slerp(a,b,t) d=a.dot(b);d=d>1 ?1:d;d=d<-1 ?-1:d;theta=Math.acos(d);sin_theta=Math.sin(theta);s0=Math.sin((1-t)*theta)/sin_theta;s1=Math.sin(t*theta)/sin_theta;new(s0*a.x+s1*b.x,s0*a.y+s1*b.y,s0*a.z+s1*b.z,s0*a.w+s1*b.w) end
  def self.rotation_m3(m) (m.zz<=0)?(m.xx>m.yy)?(t=1+m.xx-m.yy-m.zz;u=0.5/Math.sqrt(t);@x=t*u;@y=(m.xy+m.yx)*u;@z=(m.zx+m.xz)*u;@w=(m.yz-m.zy)*u):(t=1-m.xx-m.yy+m.zz;u=0.5/Math.sqrt(t);@x=(m.xy+m.yx)*u;@y=t*u;@z=(m.yz+m.zy)*u;@w=(m.zx-m.xz)*u):(-m.yy>m.xx)?(t=1-m.xx-m.yy+m.zz;u=0.5/Math.sqrt(t);@x=(m.zx+m.xz)*u;@y=(m.yz+m.zy)*u;@z=t*u;@w=(m.xy-m.yx)*u):(t=1+m.xx+m.yy+m.zz;u=0.5/Math.sqrt(t));new((m.yz-m.zy)*u,(m.zx-m.xz)*u,(m.xy-m.yx)*u,t*u) end
  def self.rotation_v3(r) hp=r.x*0.5;cp=Math.cos(hp);sp=Math.sin(hp);hy=r.y*0.5;cy=Math.cos(hy);sy=Math.sin(hy);hr=r.z*0.5;cr=Math.cos(hr);sr=Math.sin(hr);new(cr*sp*cy+sr*cp*sy,cr*cp*sy-sr*sp*cy,sr*cp*cy-cr*sp*sy,cr*cp*cy+sr*sp*sy) end # yaw pitch roll (radians)
end

class Matrix4
  attr_accessor(:xx,:xy,:xz,:xw,:yx,:yy,:yz,:yw,:zx,:zy,:zz,:zw,:wx,:wy,:wz,:ww)
  def initialize(xx=1,xy=0,xz=0,xw=0,yx=0,yy=1,yz=0,yw=0,zx=0,zy=0,zz=1,zw=0,wx=0,wy=0,wz=0,ww=1) @xx=xx;@xy=xy;@xz=xz;@xw=xw;@yx=yx;@yy=yy;@yz=yz;@yw=yw;@zx=zx;@zy=zy;@zz=zz;@zw=zw;@wx=wx;@wy=wy;@wz=wz;@ww=ww end
  def serialize; {xx:@xx,xy:@xy,xz:@xz,xw:@xw,yx:@yx,yy:@yy,yz:@yz,yw:@yw,zx:@zx,zy:@zy,zz:@zz,zw:@zw,wx:@wx,wy:@wy,wz:@wz,ww:@ww} end
  def inspect; serialize.to_s end
  def to_a; serialize.values end
  def to_s; serialize.to_s end

  def mul_m(m) xx=m.xx;xy=m.xy;xz=m.xz;xw=m.xw;yx=m.yx;yy=m.yy;yz=m.yz;yw=m.yw;zx=m.zx;zy=m.zy;zz=m.zz;zw=m.zw;wx=m.wx;wy=m.wy;wz=m.wz;ww=m.ww;Matrix4.new(@xx*xx+@xy*yx+@xz*zx+@xw*wx,@xx*xy+@xy*yy+@xz*zy+@xw*wy,@xx*xz+@xy*yz+@xz*zz+@xw*wz,@xx*xw+@xy*yw+@xz*zw+@xw*ww,@yx*xx+@yy*yx+@yz*zx+@yw*wx,@yx*xy+@yy*yy+@yz*zy+@yw*wy,@yx*xz+@yy*yz+@yz*zz+@yw*wz,@yx*xw+@yy*yw+@yz*zw+@yw*ww,@zx*xx+@zy*yx+@zz*zx+@zw*wx,@zx*xy+@zy*yy+@zz*zy+@zw*wy,@zx*xz+@zy*yz+@zz*zz+@zw*wz,@zx*xw+@zy*yw+@zz*zw+@zw*ww,@wx*xx+@wy*yx+@wz*zx+@ww*wx,@wx*xy+@wy*yy+@wz*zy+@ww*wy,@wx*xz+@wy*yz+@wz*zz+@ww*wz,@wx*xw+@wy*yw+@wz*zw+@ww*ww) end
  def mul_m!(m) xx=m.xx;xy=m.xy;xz=m.xz;xw=m.xw;yx=m.yx;yy=m.yy;yz=m.yz;yw=m.yw;zx=m.zx;zy=m.zy;zz=m.zz;zw=m.zw;wx=m.wx;wy=m.wy;wz=m.wz;ww=m.ww; _xx=@xx*xx+@xy*yx+@xz*zx+@xw*wx;_xy=@xx*xy+@xy*yy+@xz*zy+@xw*wy;_xz=@xx*xz+@xy*yz+@xz*zz+@xw*wz;_xw=@xx*xw+@xy*yw+@xz*zw+@xw*ww;_yx=@yx*xx+@yy*yx+@yz*zx+@yw*wx;_yy=@yx*xy+@yy*yy+@yz*zy+@yw*wy;_yz=@yx*xz+@yy*yz+@yz*zz+@yw*wz;_yw=@yx*xw+@yy*yw+@yz*zw+@yw*ww;_zx=@zx*xx+@zy*yx+@zz*zx+@zw*wx;_zy=@zx*xy+@zy*yy+@zz*zy+@zw*wy;_zz=@zx*xz+@zy*yz+@zz*zz+@zw*wz;_zw=@zx*xw+@zy*yw+@zz*zw+@zw*ww;_wx=@wx*xx+@wy*yx+@wz*zx+@ww*wx;_wy=@wx*xy+@wy*yy+@wz*zy+@ww*wy;_wz=@wx*xz+@wy*yz+@wz*zz+@ww*wz;_ww=@wx*xw+@wy*yw+@wz*zw+@ww*ww;@xx=_xx;@xy=_xy;@xz=_xz;@xw=_xw;@yx=_yx;@yy=_yy;@yz=_yz;@yw=_yw;@zx=_zx;@zy=_zy;@zz=_zz;@zw=_zw;@wx=_wx;@wy=_wy;@wz=_wz;@ww=_ww;self end

  def inverse; d=@xw*@yz*@zy*@wx-@xz*@yw*@zy*@wx-@xw*@yy*@zz*@wx+@xy*@yw*@zz*@wx+@xz*@yy*@zw*@wx-@xy*@yz*@zw*@wx-@xw*@yz*@zx*@wy+@xz*@yw*@zx*@wy+@xw*@yx*@zz*@wy-@xx*@yw*@zz*@wy-@xz*@yx*@zw*@wy+@xx*@yz*@zw*@wy+@xw*@yy*@zx*@wz-@xy*@yw*@zx*@wz-@xw*@yx*@zy*@wz+@xx*@yw*@zy*@wz+@xy*@yx*@zw*@wz-@xx*@yy*@zw*@wz-@xz*@yy*@zx*@ww+@xy*@yz*@zx*@ww+@xz*@yx*@zy*@ww-@xx*@yz*@zy*@ww-@xy*@yx*@zz*@ww+@xx*@yy*@zz*@ww;Matrix4.new(d*(@yz*@zw*@wy-@yw*@zz*@wy+@yw*@zy*@wz-@yy*@zw*@wz-@yz*@zy*@ww+@yy*@zz*@ww),d*(@xw*@zz*@wy-@xz*@zw*@wy-@xw*@zy*@wz+@xy*@zw*@wz+@xz*@zy*@ww-@xy*@zz*@ww),d*(@xz*@yw*@wy-@xw*@yz*@wy+@xw*@yy*@wz-@xy*@yw*@wz-@xz*@yy*@ww+@xy*@yz*@ww),d*(@xw*@yz*@zy-@xz*@yw*@zy-@xw*@yy*@zz+@xy*@yw*@zz+@xz*@yy*@zw-@xy*@yz*@zw),d*(@yw*@zz*@wx-@yz*@zw*@wx-@yw*@zx*@wz+@yx*@zw*@wz+@yz*@zx*@ww-@yx*@zz*@ww),d*(@xz*@zw*@wx-@xw*@zz*@wx+@xw*@zx*@wz-@xx*@zw*@wz-@xz*@zx*@ww+@xx*@zz*@ww),d*(@xw*@yz*@wx-@xz*@yw*@wx-@xw*@yx*@wz+@xx*@yw*@wz+@xz*@yx*@ww-@xx*@yz*@ww),d*(@xz*@yw*@zx-@xw*@yz*@zx+@xw*@yx*@zz-@xx*@yw*@zz-@xz*@yx*@zw+@xx*@yz*@zw),d*(@yy*@zw*@wx-@yw*@zy*@wx+@yw*@zx*@wy-@yx*@zw*@wy-@yy*@zx*@ww+@yx*@zy*@ww),d*(@xw*@zy*@wx-@xy*@zw*@wx-@xw*@zx*@wy+@xx*@zw*@wy+@xy*@zx*@ww-@xx*@zy*@ww),d*(@xy*@yw*@wx-@xw*@yy*@wx+@xw*@yx*@wy-@xx*@yw*@wy-@xy*@yx*@ww+@xx*@yy*@ww),d*(@xw*@yy*@zx-@xy*@yw*@zx-@xw*@yx*@zy+@xx*@yw*@zy+@xy*@yx*@zw-@xx*@yy*@zw),d*(@yz*@zy*@wx-@yy*@zz*@wx-@yz*@zx*@wy+@yx*@zz*@wy+@yy*@zx*@wz-@yx*@zy*@wz),d*(@xy*@zz*@wx-@xz*@zy*@wx+@xz*@zx*@wy-@xx*@zz*@wy-@xy*@zx*@wz+@xx*@zy*@wz),d*(@xz*@yy*@wx-@xy*@yz*@wx-@xz*@yx*@wy+@xx*@yz*@wy+@xy*@yx*@wz-@xx*@yy*@wz),d*(@xy*@yz*@zx-@xz*@yy*@zx+@xz*@yx*@zy-@xx*@yz*@zy-@xy*@yx*@zz+@xx*@yy*@zz)) end
  def transpose; xy=@xy;xz=@xz;xw=@xw;yx=@yx;yz=@yz;yw=@yw;zx=@zx;zy=@yz,zw=@zw;wx=@wx;wy=@wy;wz=@wz;Matrix4.new(@xx,xy,xz,xw,yx,@yy,yz,yw,zx,zy,@zz,zw,wx,wy,wz,@ww) end
  def transpose!; xy=@xy;xz=@xz;xw=@xw;yx=@yx;yz=@yz;yw=@yw;zx=@zx;zy=@zy;zw=@zw;wx=@wx;wy=@wy;wz=@wz;@xy=yx;@xz=zx;@xw=wx;@yx=xy;@yz=zy;@yw=wy;@zx=xz;@zy=yz;@zw=wz;@wx=xw;@wy=yw;@wz=zw end

  def self.look_at(p,t,u) px=p.x;py=p.y;pz=p.z;tx=t.x;ty=t.y;tz=t.z;ux=u.x;uy=u.y;uz=u.z;zx=tx-px;zy=ty-py;zz=tz-pz;l=1/Math.sqrt(zx*zx+zy*zy+zz*zz);zx*=l;zy*=l;zz*=l;xx=uy*zz-uz*zy;xy=uz*zx-ux*zz;xz=ux*zy-uy*zx;l=1/Math.sqrt(xx*xx+xy*xy+xz*xz);xx*=l;xy*=l;xz*=l;yx=zy*xz-zz*xy;yy=zz*xx-zx*xz;yz=zx*xy-zy*xx;new(xx,yx,zx,0,xy,yy,zy,0,xz,yz,zz,0,-(xx*px+xy*py+xz*pz),-(yx*px+yy*py+yz*pz),-(zx*px+zy*py+zz*pz),1) end
  def self.frustum(l,r,b,t,n,f) new(2*n/(r-l),0,0,0,0,2*n/(t-b),0,0,(r+l)/(r-l),(t+b)/(t-b),-(f+n)/(f-n),-1,0,0,-2*f*n/(f-n),1) end

  def self.perspective(y,a,n,f) c=Math.cos(y*0.5);s=Math.sin(y*0.5);h=c/s;r=f/(f-n);new(h/a,0,0,0,0,h,0,0,0,0,r,1,0,0,-r*n,0) end
  def self.orthographic(w,h,n,f) r=1/(f-n);new(2/w,0,0,0,0,2/h,0,0,0,0,r,0,0,0,-r*n,1) end
  def self.viewport(x,y,w,h) new(w>>1,0,0,x+(w>>1),0,h>>1,0,y+(h>>1),0,0,1,0,0,0,0,1) end

  def self.transform(t,r,s) r=rotation_q(r);sx=s.x;sy=s.y;sz=s.z;new(sx*r.xx,sx*r.xy,sx*r.xz,0,sy*r.yx,sy*r.yy,sy*r.yz,0,sz*r.zx,sz*r.zy,sz*r.zz,0,t.x,t.y,t.z,1) end
  def self.rotation_q(r) x=r.x;y=r.y;z=r.z;w=r.w;d=x*x+y*y+z*z+w*w;s=(d!=0)?2/d:1;xs=x*s;ys=y*s;zs=z*s;wx=w*xs;wy=w*ys;wz=w*zs;xx=x*xs;xy=x*ys;xz=x*zs;yy=y*ys;yz=y*zs;zz=z*zs;new(1-yy-zz,xy+wz,xz-wy,0,xy-wz,1-xx-zz,yz+wx,0,xz+wy,yz-wx,1-xx-yy,0,0,0,0,1) end
  def self.rotation_x(r) s=Math.sin(r);c=Math.cos(r);new(1,0,0,0,0,c,s,0,0,-s,c,0,0,0,0,1) end
  def self.rotation_y(r) s=Math.sin(r);c=Math.cos(r);new(c,0,-s,0,0,1,0,0,s,0,c,0,0,0,0,1) end
  def self.rotation_z(r) s=Math.sin(r);c=Math.cos(r);new(c,s,0,0,-s,c,0,0,0,0,1,0,0,0,0,1) end
  def self.translation(p) new(1,0,0,0,0,1,0,0,0,0,1,0,p.x,p.y,p.z,1) end
  def self.scale(s) new(s.x,0,0,0,0,s.y,0,0,0,0,s.z,0,0,0,0,1) end
end