
  '================
  'MATHUTIL.INC
  'matrix maths etc
  '================

  'Example of matrix layout mapped to an array of 4x4 cells
  '
  '  ordinals: (cell location offsets) in hexadecimal
  '
  '  0 4 8 C
  '  1 5 9 D
  '  2 6 A E
  '  3 7 B F
  '


  macro ShowBuffer(n)
  ===================
  string bu=space n
  end macro

  macro ShowTitle(v)
  ==================
  mid bu,c,v cr : c+=2+len v
  end macro

  macro ShowElement(v)
  ====================
  mid bu,c,"  "+str(v,4) : c+=14
  end macro

  macro ShowNextLine()
  ====================
  mid bu,c,cr : c+=2
  end macro


  macro ShowReturn()
  ==================
  return left bu,c-1
  end macro


  function ShowVec(float *v, sys n=3) as string
  =============================================
  indexbase 0
  ShowBuffer 256
  int c=1,i=0
  n--
  ShowTitle "Vector:"
  for i=0 to n
    ShowElement v[i]
  next
  ShowNextLine
  ShowReturn
  end function


  function ShowMat(float*m,sys r=4) as string
  ===========================================
  indexbase 0
  ShowBuffer 512
  string s
  sys c,e,i,j,k,n
  static double min=1.0e-7
  c=1
  n=r-1
  e=r*r-1
  ShowTitle "Matrix:"
  for i=0 to n 'y
    for j=0 to e step r 'x
      k=i+j
      ShowElement m[k]
    next
    ShowNextLine
  next
  ShowReturn
  end function


  'function Transpose(float *r, *m, int n=4)
  '=========================================
  'indexbase 0
  'int i,j,k1,k2
  'int e=n-1
  'for i=0 to e
  '  k2=i
  '  for j=0 to e
  '    r[k1]=m[k2] : k1++ : k2+=n
  '  next
  'next  
  'end function


  function Transpose(float *r, *m, int nx=4,ny=0) as float*
  =========================================================
  ' r  resulting transposed matrix ptr
  ' m  source matrix ptr
  ' nx source matrix columns
  ' ny source matrix rows
  '
  'using striding pointers instead of indexes
  '
  int i  ' major iterator
  int j  ' minor iterator
  int mx ' column stride
  int my ' row stride
  sys*mm ' striding source ptr
  sys*rr ' striding transpose ptr
  '
  if ny=0
    ny=nx
  endif
  mx=nx*sizeof float
  my=sizeof float
  @rr=@r
  @mm=@m
  for i=0 to <nx   'transposed rows
    for j=0 to <ny 'tranposed column
      rr=mm
      @rr+=my 
      @mm+=mx
    next
    @mm=@mm+my-(ny*mx)
    '@mm=@m
  next
  return @r
  end function


  function SetIdentity(float *m, int n=4)
  =======================================
  indexbase 0
  int i=0, h=n+1, e=n*n-1
  for i=0 to e
    m[i]=0.0
  next
  for i=0 to e step h
    m[i]=1.0
  next
  end function


  'generic with striding pointers
  'def typf float
  'typedef float typf
  '
  function MatMul(float *r,*a,*b, int n=4) 'NxN MATRIX : N=1..
  ============================================================
  int ystep=sizeof float
  int xstep=n*sizeof float
  int i,j,k
  sys px
  for i=1 to n
    px=@a
    for j=1 to n
      r=0
      for k=1 to n
        r+=(a*b)
        @a+=xstep
        @b+=ystep
      next
      @r+=ystep
      px+=ystep
      @a=px
      @b-=xstep
    next
    @a-=xstep
    @b+=xstep
  next
  end function


  function MatCol float *r,*m, int c=0,n=4)
  =========================================
  int rstep=4
  int cstep=n*4
  int i
  @m=@m+c*cstep
  for i=1 to n
    r=m
    @r+=rstep  'NEXT CELL
    @m+=rstep  'NEXT ROW Y
  next
  end function

  
  function MatRow float *r,*m, int rw=0,n=4)
  ==========================================
  int rstep=4
  int cstep=n*4
  int i
  @m=@m+rw*rstep
  for i=1 to n
    r=m
    @r+=rstep  'NEXT CELL
    @m+=cstep  'NEXT COLUMN X
  next
  end function


  function InverseMat(float *m, sys n=4) 'NxN MATRIX : N=2..16
  ============================================================

  'Uses the Gauss-Jordan Method.
  'Inverts matrix, m[], "in-place".

  'limit to max 16x16 matrix

  int     pivot,row,col,maxrow
  double  test,colmax,divisor,factor,temp
  int     rowswitch,switchcount
  int     i,j
  int     rs[128] 'for limit array n=64

  '
  'ARRAY DEFINITIONS (COL,ROW) == (X,Y)
  '====================================
  '
  indexbase 0
  #define mat(row,col) m[ col*n+row ]
  #define rowswitches(row,col) rs[ col*64+row ]
  int first=0
  int last=n-1
  '
  'TRAVERSE COLUMNS


  for pivot=first to last 'COLUMN
  '
  '
  'GET THE MAXIMUM ELEMENT ON EACH PIVOT COLUMN
  '
  colmax=0 'ABSOLUTE MAX
  maxrow=pivot
  '
  for row=pivot to last
    test=abs (mat(row,pivot))
    if test>colmax then
      colmax=test
      maxrow=row
    end if
  next

  if maxrow<>pivot then
    '
    '
    'SWAP ROWS TO GET THE MAX VALUE AT THE PIVOT
    '
    for col=first to last
      temp=mat(pivot,col)
      mat(pivot,col)=mat(maxrow,col)
      mat(maxrow,col)=temp
    next
    '
    '
    'REMEMBER THE ROW SWAP
    '
    rowswitches(switchcount,0)=pivot
    rowswitches(switchcount,1)=maxrow
    switchcount++
  end if
  '
  '
  'DIVIDE THE PIVOT ROW BY THE PIVOT DIAGONAL
  '
  'also setting the pivot diagonal to 1
  '
  divisor=1.0/mat(pivot,pivot)
  mat(pivot,pivot)=1
  '
  for col=first to last
    mat(pivot,col)=mat(pivot,col) * divisor
  next


  'APPLY PIVOT FACTORS TO ALL ELEMENTS (OTHER THAN PIVOT ROW)
  '
  'also setting each pivot row element to 0
  '
  for row=first to last
    if row<>pivot then
      factor=-mat(row,pivot)
      mat(row,pivot)=0
      for col=first to last
        mat(row,col)=mat(row,col)+mat(pivot,col) * factor
      next 'COL
    end if
  next 'ROW
  '
  next 'PIVOT
  '
  'SWAP COLUMN VALUES ACCORDING TO ORIGINAL ROW POSITIONS
  '
  for rowswitch=switchcount-1 to 0 step-1
    i=rowswitches(rowswitch,0)
    j=rowswitches(rowswitch,1)
    for row=first to last
      temp = mat(row,i)
      mat(row,i)=mat(row,j)
      mat(row,j)=temp
    next 
  next
  '
  end function


  function det2x2sub(float *m, int i0,i1,i2,i3) as float
  ======================================================
  indexbase 0
  return m[i0] * m[i3] - m[i2] * m[i1]
  end function


  function InverseMat3x3(float *r, *m)
  ====================================
  indexbase 0
  float det
  '
  det += m[0] * det2x2sub(m, 4, 5, 7, 8)
  det -= m[3] * det2x2sub(m, 1, 2, 7, 8)
  det += m[6] * det2x2sub(m, 1, 2, 4, 5)
  '
  r[0] =  det2x2sub(m, 4, 5, 7, 8) / det
  r[1] = -det2x2sub(m, 1, 2, 7, 8) / det
  r[2] =  det2x2sub(m, 1, 2, 4, 5) / det
  r[3] = -det2x2sub(m, 3, 5, 6, 8) / det
  r[4] =  det2x2sub(m, 0, 2, 6, 8) / det
  r[5] = -det2x2sub(m, 0, 2, 3, 5) / det
  r[6] =  det2x2sub(m, 3, 4, 6, 7) / det
  r[7] = -det2x2sub(m, 0, 1, 6, 7) / det
  r[8] =  det2x2sub(m, 0, 1, 3, 4) / det
  end function


  function Mat4x4ToMat3x3(float *r,*m)
  ====================================
  indexbase 0
  'copy top left section
  r[0] = m[0]
  r[3] = m[4]
  r[6] = m[8]
  r[1] = m[1]
  r[4] = m[5]
  r[7] = m[9]
  r[2] = m[2]
  r[5] = m[6]
  r[8] = m[10]
  end function


  function mat4x4Mul(float *r,*m1,*m2)
  ====================================
  indexbase 0
  'inline version
  ' r[0]  = m1[0] * m2[0]  + m1[4] * m2[1]  + m1[8]  * m2[2]  + m1[12] * m2[3]
  ' r[1]  = m1[1] * m2[0]  + m1[5] * m2[1]  + m1[9]  * m2[2]  + m1[13] * m2[3]
  ' r[2]  = m1[2] * m2[0]  + m1[6] * m2[1]  + m1[10] * m2[2]  + m1[14] * m2[3]
  ' r[3]  = m1[3] * m2[0]  + m1[7] * m2[1]  + m1[11] * m2[2]  + m1[15] * m2[3]
  ' r[4]  = m1[0] * m2[4]  + m1[4] * m2[5]  + m1[8]  * m2[6]  + m1[12] * m2[7]
  ' r[5]  = m1[1] * m2[4]  + m1[5] * m2[5]  + m1[9]  * m2[6]  + m1[13] * m2[7]
  ' r[6]  = m1[2] * m2[4]  + m1[6] * m2[5]  + m1[10] * m2[6]  + m1[14] * m2[7]
  ' r[7]  = m1[3] * m2[4]  + m1[7] * m2[5]  + m1[11] * m2[6]  + m1[15] * m2[7]
  ' r[8]  = m1[0] * m2[8]  + m1[4] * m2[9]  + m1[8]  * m2[10] + m1[12] * m2[11]
  ' r[9]  = m1[1] * m2[8]  + m1[5] * m2[9]  + m1[9]  * m2[10] + m1[13] * m2[11]
  ' r[10] = m1[2] * m2[8]  + m1[6] * m2[9]  + m1[10] * m2[10] + m1[14] * m2[11]
  ' r[11] = m1[3] * m2[8]  + m1[7] * m2[9]  + m1[11] * m2[10] + m1[15] * m2[11]
  ' r[12] = m1[0] * m2[12] + m1[4] * m2[13] + m1[8]  * m2[14] + m1[12] * m2[15]
  ' r[13] = m1[1] * m2[12] + m1[5] * m2[13] + m1[9]  * m2[14] + m1[13] * m2[15]
  ' r[14] = m1[2] * m2[12] + m1[6] * m2[13] + m1[10] * m2[14] + m1[14] * m2[15]
  ' r[15] = m1[3] * m2[12] + m1[7] * m2[13] + m1[11] * m2[14] + m1[15] * m2[15]
  'exit function
  '
  'generic version
  'matmul(r,m1,m2,4)
  'exit function
  '
  'assembler version 
  addr rcx,r
  addr rsi,m1
  addr rdi,m2
  jmp fwd begin
  ----
  row:
  ----
  fld  dword [esi   ]
  fmul dword [edi   ]
  '
  fld  dword [esi+16]
  fmul dword [edi+ 4]
  faddp
  fld  dword [esi+32]
  fmul dword [edi+ 8]
  faddp
  fld  dword [esi+48]
  fmul dword [edi+12]
  faddp
  fstp dword [ecx]
  add ecx,4
  add esi,4
  ret
  ------
  block:
  ------
  call row
  call row
  call row
  call row
  sub esi,16
  add edi,16
  ret
  ------
  begin:
  ------
  call block
  call block
  call block
  call block
  end function


  function transformVec4(float *r,*m,*v, sys n=0)
  ===============================================
  'using mat4x4
  indexbase 0
  sys i
  for i=0 to n
    r[0]  = m[0] * v[0]  + m[4] * v[1]  + m[8]  * v[2]  + m[12] * v[3]
    r[1]  = m[1] * v[0]  + m[5] * v[1]  + m[9]  * v[2]  + m[13] * v[3]
    r[2]  = m[2] * v[0]  + m[6] * v[1]  + m[10] * v[2]  + m[14] * v[3]
    r[3]  = m[3] * v[0]  + m[7] * v[1]  + m[11] * v[2]  + m[15] * v[3]
    @v+=16
    @r+=16
  next
  /*
  'assembler version
  addr rcx,r
  addr rsi,m
  addr rdi,v
  jmp fwd begin
  ----
  row:
  ----
  fld  dword [esi   ]
  fmul dword [edi   ]
  fld  dword [esi+16]
  fmul dword [edi+ 4]
  fld  dword [esi+32]
  fmul dword [edi+ 8]
  fld  dword [esi+48]
  fmul dword [edi+12]
  faddp st1,st0
  faddp st1,st0
  faddp st1,st0
  fstp dword [ecx]
  add ecx,4
  add esi,4
  ret
  ------
  block:
  ------
  call row
  call row
  call row
  call row
  sub esi,16 'm
  add edi,16 'v
  ret
  ------
  begin:
  ------
  (
  call block
  dec i
  jge repeat
  )
  */
 end function


  function transformVec3(float *r,*m,*v, sys n=1)
  ===============================================
  'using mat4x4
  indexbase 0
  sys i
  n--
  for i=0 to n
    r[0]  = m[0] * v[0]  + m[4] * v[1]  + m[8]  * v[2]  + m[12]
    r[1]  = m[1] * v[0]  + m[5] * v[1]  + m[9]  * v[2]  + m[13]
    r[2]  = m[2] * v[0]  + m[6] * v[1]  + m[10] * v[2]  + m[14]
    @v+=12
    @r+=12
  next
  end function

/*
  macro rotsetup
  ==============
  indexbase 0
  float m[16]
  a = rad(a)
  float ca = cos(a), sa = sin(a)
  end macro

  function mat4x4rotateVecX(float *r, *v, a)
  ==========================================
  '
  'x  1  0  0  0
  'y  0  c -s  0
  'z  0  s  c  0
  '   0  0  0  0
  '
  rotsetup
  m[0]  = 1
 'm[1]  = 0
 'm[2]  = 0
 'm[3]  = 0
  ----------
 'm[4]  = 0
  m[5]  = ca
  m[6]  = sa
 'm[7]  = 0
  ----------
 'm[8]  = 0
  m[9]  = -sa
  m[10] = ca
 'm[11] = 0
  transformVec3 r,m,v
  end function


  function mat4x4rotateVecY(float *r, *v, a)
  ==========================================
  '
  'x  c  0  s  0
  'y  0  1  0  0
  'z -s  0  c  0
  '   0  0  0  0
  '
  rotsetup
  m[0]  = ca
 'm[1]  = 0
  m[2]  = -sa
 'm[3]  = 0
  ----------
 'm[4]  = 0
  m[5]  = 1
 'm[6]  = 0
 'm[7]  = 0
  ----------
  m[8]  = sa
 'm[9]  = 0
  m[10] = ca
 'm[11] = 0
  transformVec3 r,m,v
  end function


  function mat4x4rotateVecZ(float *r, *v, a)
  ==========================================
  '
  'x  c -s  0  0
  'y  s  c  0  0
  'z  0  0  1  0
  '   0  0  0  0
  '
  rotsetup
  m[0]  = ca
  m[1]  = sa
 'm[2]  = 0
 'm[3]  = 0
  ----------
  m[4]  = -sa
  m[5]  = ca
 'm[6]  = 0
 'm[7]  = 0
  ----------
 'm[8]  = 0
 'm[9]  = 0
  m[10] = 1
 'm[11] = 0
  transformVec3 r,m,v
  end function
*/


  function det3x3sub(float *m, int i0,i1,i2,i3,i4,i5,i6,i7,i8) as float
  =====================================================================
  indexbase 0
  float det
  det += m[i0] * det2x2sub(m, i4, i5, i7, i8)
  det -= m[i3] * det2x2sub(m, i1, i2, i7, i8)
  det += m[i6] * det2x2sub(m, i1, i2, i4, i5)
  return det;
  end function


  function mat4x4Inverse(float *r,*m)
  ===================================
  indexbase 0
  float det
  det += m[0]  * det3x3sub(m, 5, 6, 7, 9, 10, 11, 13, 14, 15)
  det -= m[4]  * det3x3sub(m, 1, 2, 3, 9, 10, 11, 13, 14, 15)
  det += m[8]  * det3x3sub(m, 1, 2, 3, 5,  6,  7, 13, 14, 15)
  det -= m[12] * det3x3sub(m, 1, 2, 3, 5,  6,  7,  9, 10, 11)

  r[0]  =  det3x3sub(m, 5, 6, 7, 9, 10, 11, 13, 14, 15) / det;
  r[1]  = -det3x3sub(m, 1, 2, 3, 9, 10, 11, 13, 14, 15) / det;
  r[2]  =  det3x3sub(m, 1, 2, 3, 5,  6,  7, 13, 14, 15) / det;
  r[3]  = -det3x3sub(m, 1, 2, 3, 5,  6,  7,  9, 10, 11) / det;
  r[4]  = -det3x3sub(m, 4, 6, 7, 8, 10, 11, 12, 14, 15) / det;
  r[5]  =  det3x3sub(m, 0, 2, 3, 8, 10, 11, 12, 14, 15) / det;
  r[6]  = -det3x3sub(m, 0, 2, 3, 4,  6,  7, 12, 14, 15) / det;
  r[7]  =  det3x3sub(m, 0, 2, 3, 4,  6,  7,  8, 10, 11) / det;
  r[8]  =  det3x3sub(m, 4, 5, 7, 8,  9, 11, 12, 13, 15) / det;
  r[9]  = -det3x3sub(m, 0, 1, 3, 8,  9, 11, 12, 13, 15) / det;
  r[10] =  det3x3sub(m, 0, 1, 3, 4,  5,  7, 12, 13, 15) / det;
  r[11] = -det3x3sub(m, 0, 1, 3, 4,  5,  7,  8,  9, 11) / det;
  r[12] = -det3x3sub(m, 4, 5, 6, 8,  9, 10, 12, 13, 14) / det;
  r[13] =  det3x3sub(m, 0, 1, 2, 8,  9, 10, 12, 13, 14) / det;
  r[14] = -det3x3sub(m, 0, 1, 2, 4,  5,  6, 12, 13, 14) / det;
  r[15] =  det3x3sub(m, 0, 1, 2, 4,  5,  6,  8,  9, 10) / det;
  end function


  function mat4x4scale(float *r, x, y, z)
  =======================================
  indexbase 0
  r[0]  = x
  r[5]  = y
  r[10] = z
  r[15] = 1.0
  end function


  function mat4x4Translate(float *r, x, y, z)
  ===========================================
  indexbase 0
  r[ 0] = 1.0
  r[ 5] = 1.0
  r[10] = 1.0
  r[15] = 1.0
  r[12] = x
  r[13] = y
  r[14] = z
  end function


  function mat4x4rotate(float *r, a, x, y, z)
  ===========================================
  indexbase 0
  a = rad(a)
  'NORMALISE
  float f ' =sqr(x*x+y*y+z*z)
    fld dword x
    fmul st0
    fld dword y
    fmul st0
    fld dword z
    fmul st0
    faddp
    faddp
    fsqrt
    fstp dword f
  if f then f=1/f
  x*=f : y*=f : z*=f
  '
  'BASE CALCS
  float sa ' = sin(a)
  float ca ' = cos(a)
  float ka ' = 1-cos(a)
    fld dword  a
    fsincos
    fst dword ca
    fld1
    fsubrp st1
    fstp dword ka
    fstp dword sa
  '
  '
  float sx=sa*x,    sy=sa*y,    sz=sa*z
  float kxy=ka*x*y, kxz=ka*x*z, kyz=ka*y*z
  float kxx=ka*x*x, kyy=ka*y*y, kzz=ka*z*z
  '
  https://docs.gl/gl3/glRotate
  '---------------------------------
  'MATRIX LAYOUT 3D
  '=================================
  ' kxx+ca , kxy-sz , kxz+sy  , 0.0
  '---------------------------------
  ' kxy+sz , kyy+ca , kyz-sx  , 0.0
  '---------------------------------
  ' kxz-sy , kyz+sx , kzz+ca  , 0.0
  '---------------------------------
  ' 0.0     , 0.0     , 0.0   , 1.0
  '=================================
  '
  r[0]  = kxx + ca
  r[1]  = kxy + sz
  r[2]  = kxz - sy
  r[3]  = 0.0
  '
  r[4]  = kxy - sz
  r[5]  = kyy + ca
  r[6]  = kyz + sx
  r[7]  = 0.0
  '
  r[8]  = kxz + sy
  r[9]  = kyz - sx
  r[10] = kzz + ca
  r[11] = 0.0
  '
  r[12] = 0.0
  r[13] = 0.0
  r[14] = 0.0
  r[15] = 1.0
  '
  end function


  type vec3 float x,y,z
  =====================

  function Normalize(vec3 *r,*v)
  ==============================
  double f=sqr(v.x*v.x+v.y*v.y+v.z*v.z)
  if f then f=1/f
  r.x=v.x*f
  r.y=v.y*f
  r.z=v.z*f
  end function

  function Normalizef(float *r,*v)
  ================================
  indexbase 0
  double f=sqr(v[0]*v[0]+v[1]*v[1]+v[2]*v[2])
  if f then f=1/f
  r[0]=v[0]*f
  r[1]=v[1]*f
  r[2]=v[2]*f
  end function


  function dot(vec3 *u,*v) as float
  =================================
  return u.x * v.x + u.y * v.y + u.z * v.z
  end function


  function cross(vec3 *r,*u,*v)
  =============================
  r.x = u.y * v.z - u.z * v.y
  r.y = u.z * v.x - u.x * v.z
  r.z = u.x * v.y - u.y * v.x
  end function


  function diff(vec3 *r,*u,*v)
  ============================
  r.x=u.x - v.x
  r.y=u.y - v.y
  r.z=u.z - v.z
  end function

  function mat4x4Look(float *r, vec3 *eye,*center,*up)
  ====================================================
  indexbase 0
  vec3 vx,vy,vz,vt
  diff(vt,eye,center)
  normalize(vz,vt)
  '
  cross(vt,vz,up)
  normalize(vx,vt)
  '
  cross(vy,vz,vx)
  '
  r[0]  = vx.x
  r[1]  = vy.x
  r[2]  = vz.x
  r[4]  = vx.y
  r[5]  = vy.y
  r[6]  = vz.y
  r[8]  = vx.z
  r[9]  = vy.z
  r[10] = vz.z
  r[12] = -dot(vx, eye)
  r[13] = -dot(vy, eye)
  r[14] = -dot(vz, eye)
  end function


  function CopyFloat(float *r,*s, int n=1)
  ========================================
  copy @r,@s,n*4
  end function


  function CopyVec(float *r,*s, int n=1)
  ======================================
  copy @r,@s,n*12
  end function


  function CopyMat(float *r,*s, int n=1)
  ======================================
  'copy @r,@s,n*64
  mov ecx,n
  addr rsi,s
  addr rdi,r
  (
    dec ecx
    jl exit
    movups xmm0,[rsi   ]
    movups [rdi   ],xmm0
    movups xmm0,[rsi+16]
    movups [rdi+16],xmm0
    movups xmm0,[rsi+32]
    movups [rdi+32],xmm0
    movups xmm0,[rsi+48]
    movups [rdi+48],xmm0
    add rsi,64
    add rdi,64
    repeat
  )
  end function



  '
  'CAMERA SPACE
  'http://ogldev.atspace.co.uk/www/tutorial13/tutorial13.html
  '
  function SetViewMatrix(float *r, vec3 *vp,*vu,*vv,*vn)
  ======================================================
  indexbase 0
  vec3  u,v,n
  float p[16], m[16]
  '
  SetIdentity p,4
  p[12]={-vp.x, -vp.y, -vp.z}
  '
  normalize u,vu
  normalize v,vv
  normalize n,vn
  'itr cross product u from n and v
  '
  '      U    V      N
  '      UP   RIGHT  LOOKAT
  --------------------------------------
  m[0] = {u.x, v.x,  n.x, 0.0} 'column 0
  m[4] = {u.y, v.y,  n.y, 0.0} 'column 1
  m[8] = {u.z, v.z,  n.z, 0.0} 'column 2
  m[12]= {0.0, 0.0,  0.0, 1.0} 'column 3
  --------------------------------------
  '
  Mat4x4Mul r,m,p
  '
  end function



  'http://stackoverflow.com/questions/18404890/how-to-build-perspective-projection-matrix-no-api
  'http://gamedev.stackexchange.com/questions/56609/how-to-create-a-projection-matrix-in-opengl-es-2-0
  '**
  'http://www.songho.ca/opengl/gl_projectionmatrix.html
  'http://www.terathon.com/gdc07_lengyel.pdf
  '
  'http://www.swiftless.com/tutorials/opengl4/3-opengl-4-matrices.html
  '
  '
  function SetProjectionMatrix(float *r, fov, aspect, near, far, int orth=0, leftHanded=1)
  =========================================================================================
  'FOV RECEIVED IN DEGREES
  'ASPECT == VIEWPORT WIDTH/HEIGHT
  '
  indexbase 0
  SetIdentity r,4
  float h = 1 : if not lefthanded then h=-1
  float i = 1 / (far-near)
  float f = 1 / tan(rad(fov/2))
  '
  r[5] = aspect * f 'uh
  r[0] = h * f      'uw
  '
  if orth then
  '
  'ORTHOGRAPHIC PROJECTION
  -------------------------------------------------
  '  uw        0         0            0
  '  0         uh        0            0
  '  0         0         -2/(f-n)     -(f+n)/(f-n)
  '  0         0         0            1
  -------------------------------------------------
  '
  r[10] = -2.0 * i
 'r[11] = 0
  r[14] = -(far+near) * i
  r[15] = 1
  exit function
  end if
  '
  'PERSPECTIVE PROJECTION
  -------------------------------------------------
  '  uw        0         0            0
  '  0         uh        0            0
  '  0         0         -(f+n)/(f-n) -2fn/(f-n)
  '  0         0         -1           0
  -------------------------------------------------
  r[10] = -(far+near) * i
  r[11] = -1
  r[14] = -2.0 * far * near * i
 'r[15] = 0
  end function


  ==============================================
  'FUNCTIONS EMULATING CLASSIC OPENGL TRANSFORMS
  ==============================================


  float MatStack[800], sys MatIdx 'stack depth: 50

  SetIdentity MatStack,4 : MatIdx=0

  def mats MatStack[MatIdx]

  function matf() as float*
  =========================
  indexbase 0
  return @matstack[matidx]
  end function


  function scalef(float x,y,z)
  ============================
  indexbase 0
  float m[32]
  mat4x4scale m[0],x,y,z
  mat4x4mul m[16],MatStack[MatIdx],m[0]
  copymat MatStack[MatIdx],m[16]
  end function


  function translatef(float x,y,z)
  ================================
  indexbase 0
  float m[32]
  mat4x4translate m[0],x,y,z
  mat4x4mul m[16], MatStack[MatIdx],m[0]
  copymat MatStack[MatIdx],m[16]
  end function


  function rotatef(float a,x,y,z)
  ===============================
  indexbase 0
  float m[32]
  mat4x4rotate m[0],a,x,y,z
  mat4x4mul m[16],MatStack[MatIdx],m[0]
  copymat MatStack[MatIdx],m[16]
  end function


  function PushMatrix()
  =====================
  indexbase 0
  sys msd=MatIdx+16
  if msd<800 then
    copymat matstack[msd], matstack[MatIdx]
    MatIdx=msd
  end if
  end function


  function PopMatrix()
  ====================
  if MatIdx>0 then MatIdx-=16
  end function


  function LoadMatrix(float*m)
  ============================
  indexbase 0
  copymat matstack[MatIdx], m
  end function


  function SaveMatrix(float*m)
  ============================
  indexbase 0
  copymat m,matstack[MatIdx]
  end function


  function LoadIdentity()
  =======================
  indexbase 0
  SetIdentity matstack[MatIdx],4
  end function


  'SUPPORTING FUNCTIONS


  function MulMatrix(float *m)
  ============================
  indexbase 0
  float r[16]
  mat4x4mul r,MatStack[MatIdx],m
  copymat MatStack[MatIdx],r
  end function


  function inversef()
  ===================
  indexbase 0
  float m[16]
  mat4x4inverse m[0],MatStack[MatIdx]
  copymat MatStack[MatIdx],m[0]
  end function



  function transformf(float *r,*v, int n=1,st=12,oo=0)
  ====================================================
  '
  indexbase 0
  float m at (@MatStack[MatIdx])
  sys i
  @v+=oo 'OFFSET BYTES
  n-- 'FINAL VERTEX INDEX
  for i=0 to n
    r[0]  = m[0] * v[0]  + m[4] * v[1]  + m[8]  * v[2]  + m[12]
    r[1]  = m[1] * v[0]  + m[5] * v[1]  + m[9]  * v[2]  + m[13]
    r[2]  = m[2] * v[0]  + m[6] * v[1]  + m[10] * v[2]  + m[14]
    @v+=st 'STRIDE BYTES
    @r+=st
  next
  end function

  function NormalMatf(float *r,*s)
  ===============================
  'FROM GLSL 
  'mat4 normalMatrix = transpose(inverse(modelView));
  'vertex.normal   = normalize ( transpose ( inverse ( mat3 ( m ) ) ) * v_normal );
  '
  indexbase 0
  float k[9],m[9]
  m={s[0],s[1],s[2],s[4],s[5],s[6],s[8],s[9],s[10]}
  InverseMat3x3(k,m)
  Transpose(r,k,3,3)
  end function


  function TranslationMatf(float *r,*s)
  =====================================
  'BillBoarding
  indexbase 0
  float o=0.0, l=1.0
  r={l,o,o,o,  o,l,o,o,  o,o,l,o,  s[12],s[13],s[14],l}
  end function


  function transnormf(float *r,*v, int n=1,st=12,oo=0)
  ====================================================
  '
  indexbase 0
  float m at (@MatStack[MatIdx])
  sys i
  @v+=oo 'OFFSET BYTES
  n-- 'FINAL VERTEX INDEX
  for i=0 to n
    r[0]  = m[0] * v[0]  + m[4] * v[1]  + m[8]  * v[2]
    r[1]  = m[1] * v[0]  + m[5] * v[1]  + m[9]  * v[2]
    r[2]  = m[2] * v[0]  + m[6] * v[1]  + m[10] * v[2]
    normalizef r,r
    @v+=st 'STRIDE BYTES
    @r+=st
  next
  end function


  function transnormfA(float *r,*v, int n=1,st=12,oo=0)
  =====================================================
  '
  indexbase 0
  'float s at (@MatStack[MatIdx])
  float m[9]
  sys i
  '
  NormalMatf(m,MatStack[Matidx])
  '
  @v+=oo 'OFFSET BYTES
  n-- 'FINAL VERTEX INDEX
  '
 
  for i=0 to n
    r[0]  = m[0] * v[0]  + m[3] * v[1]  + m[6]  * v[2]
    r[1]  = m[1] * v[0]  + m[4] * v[1]  + m[7]  * v[2]
    r[2]  = m[2] * v[0]  + m[5] * v[1]  + m[8]  * v[2]
    normalizef r,r
    @v+=st 'STRIDE BYTES
    @r+=st
  next
  end function


  function translormf(float *r,*v, int n=1,st=12,oo=0)
  =================================================
  '
  indexbase 0
  float m at (@MatStack[MatIdx])
  sys i
  @v+=oo 'OFFSET BYTES
  n-- 'FINAL VERTEX INDEX
  for i=0 to n
    r[0] = v[0] + m[12]
    r[1] = v[1] + m[13]
    r[2] = v[2] + m[14]
    @v+=st 'STRIDE BYTES
    @r+=st
  next
  end function


'https://learnopengl.com/#!Advanced-Lighting/Normal-Mapping
'http://www.opengl-tutorial.org/intermediate-tutorials/tutorial-13-normal-mapping/
'https://gamedev.stackexchange.com/questions/68612/how-to-compute-tangent-and-bitangent-vectors


  sub diff2(float *r,*v,*w)
    r[0]=v[0]=w[0]
    r[1]=v[1]=w[1]
  end sub

  sub diff3(float *r,*v,*w)
    r[0]=v[0]-w[0]
    r[1]=v[1]-w[1]
    r[2]=v[2]-w[2]
  end sub

  sub unify(float *v)
    float f = sqr( v[0]*v[0] + v[1]*v[1] + v[2]*v[2] )
    if f then f=1/f
    v[0]*=f : v[1]*=f : v[2]*=f
  end sub


  function TanBitanNorm(float *r, *p, *t)
  =======================================
  indexbase 0
  float e[9],u[6]
  '
  diff3 e[0],p[3],p[0]
  diff3 e[3],p[6],p[0]
  diff2 u[0],t[2],t[0]
  diff2 u[2],t[4],t[0]
  float f= (u[0]*u[3]) - (u[2]*u[1])
  if f then f=1/f
  'TAN
  'tangent1.x = f * (deltaUV2.y * edge1.x - deltaUV1.y * edge2.x);
  'tangent1.y = f * (deltaUV2.y * edge1.y - deltaUV1.y * edge2.y);
  'tangent1.z = f * (deltaUV2.y * edge1.z - deltaUV1.y * edge2.z);
  r[0] = f * ( ( u[3] * e[0]) - (u[1] * e[3]) )  
  r[1] = f * ( ( u[3] * e[1]) - (u[1] * e[4]) )  
  r[2] = f * ( ( u[3] * e[2]) - (u[1] * e[5]) )  
  unify r[0]
  'BITAN
  'bitangent1.x = f * (-deltaUV2.x * edge1.x + deltaUV1.x * edge2.x);
  'bitangent1.y = f * (-deltaUV2.x * edge1.y + deltaUV1.x * edge2.y);
  'bitangent1.z = f * (-deltaUV2.x * edge1.z + deltaUV1.x * edge2.z);
  r[3] = f * ( (-u[2] * e[0]) + (u[0] * e[3]) )  
  r[4] = f * ( (-u[2] * e[1]) + (u[0] * e[4]) )  
  r[5] = f * ( (-u[2] * e[2]) + (u[0] * e[5]) )  
  unify r[3]
  'NORMAL
  diff3(e[0],p[3],p[0])
  diff3(e[3],p[6],p[0])
  r[6]=(e[1]*e[5])-(e[2]*e[4]) ' py*qz - pz*qy
  r[7]=(e[2]*e[3])-(e[0]*e[5]) ' pz*qx - px*qz
  r[8]=(e[0]*e[4])-(e[1]*e[3]) ' px*qy - py*qx
  unify r[6]
  end function  


  function sortindex(float *v, sys *idx, sys mx)
  ==============================================
  'MERGE SORT
  'idx must have 2*mx capacity
  '
  indexbase 0
  sys i,j,k
  sys src,dest,pile,group,sw,kk
  sys bi,bj,ei,ej,ni,nj
  '
  'INITIALISE INDEX
  for i=0 to <mx
    idx[i]=i
  next
  '
  src=0
  dest=mx
  pile=1
  group=pile+pile
  do
    if pile>=mx then
      if src then
        for i=0 to <mx
          j=src+i
          idx[i]=idx[j]
        next
      end if
      exit do
    end if
    bi=src
    bj=bi+pile
    ei=bi+pile
    ej=bj+pile
    ni=bi
    nj=bj
    k=src+mx
    if ei>k then ei=k 'CLAMP
    if ej>k then ej=k 'CLAMP
    kk=0
    do
      sw=0
      if kk>mx then exit do
      if ni>=ei then
        if nj>=ej then
          'MOVE TO NEXT GROUP
          k=src+mx 'BOUNDARY
          bi+=group : ni=bi : ei+=group
          bj+=group : nj=bj : ej+=group
          if ni>=k then exit do
          if ei>k then ei=k 'CLAMP
          if ej>k then ej=k 'CLAMP
          continue do
        else
          sw=1  'TAKE FROM J
        end if
      elseif nj>=ej then
        'TAKE FROM I
      else
        'COMPARE
        i=idx[ni]
        j=idx[nj]
        if v[i]>v[j] then sw=1 'ASCENDING
      end if
      if sw then
        j=idx[nj] : nj++
      else
        j=idx[ni] : ni++
      end if
      k=dest+kk : idx[k]=j : kk++
    end do
    j=src : src=dest : dest=j 'SWAP DEST AND SOURCE INDEX BUFFER POSITION
    pile=group : group+=group 'DOUBLE STRIDE
  end do
  end function



  #ifdef MathUtilExtra



  '===========
  'QUATERNIONS
  '===========
  '
  'http://ogldev.atspace.co.uk/www/tutorial15/tutorial15.html
  'http://www.cprogramming.com/tutorial/3d/quaternions.html
  '
  ' Hamilton
  '
  'SLERP, or Spherical Linear intERPolation
  '
  ' Q = xi + yj + zk + w
  '
  ' i*i == j*j == k*k = ijk == -1
  '
  ' 1/Q == -xi - yi - zk +w
  '
  'ROTATION QUATERNION
  '
  ' W = QVQ^-1
  '
  ' Q = {
  ' v.x * sin(a/2),
  ' v.y * sin(a/2),
  ' v.z * sin(a/2),
  ' cos(a/2)
  ' }

  function NormQuaternion(float *r)
  =================================
  indexbase 0
  double f = sqr(r[0]*r[0] + r[1]*r[1] + r[2]*r[2] + r[3]*r[3])
  if f then f=1/f
  r[0]*=f
  r[1]*=f
  r[2]*=f
  r[3]*=f
  end function


  function SetQuaternion(float *r,a,x,y,z)
  ========================================
  indexbase 0
  float ha=rad(a/2) 'degrees to radians of half angle
  float sa=sin(ha)
  double f=sqr(x*x+y*y+z*z)
  if f then f=1/f
  x*=f : y*=f : z*=z
  r[0]=x*sa
  r[1]=y*sa
  r[2]=z*sa  
  r[3]=cos(ha)
  NormQuaternion r
  end function


  function ResetQuaternion(float *r)
  ==================================
  r={0,0,0,1}
  end function


  function MulQuaternion(float *r,*a,*b)
  ======================================
  '0 x
  '1 y
  '2 z
  '3 w
  indexbase 0
  'itr test first for deviation from unit
  NormQuaternion a
  NormQuaternion b
  '
  r[3] = a[3]*b[3] - a[0]*b[0] - a[1]*b[1] - a[2]*[b[2]
  r[0] = a[3]*b[0] + a[0]*b[3] + a[1]*b[2] - a[2]*[b[1]
  r[1] = a[3]*b[1] - a[0]*b[2] + a[1]*b[3] + a[2]*[b[0]
  r[2] = a[3]*b[2] + a[0]*b[1] - a[1]*b[0] + a[2]*[b[3]
  end function

  'generate local_rotation (details later)
  'Qtotal = Qlocal * Qtotal  //multiplication order matters on this line
  '
  '
  function MulQuaternion(float *qTotal,*qLocal)
  =============================================
  float q[4]
  MulQuaternion q,qLocal,qTotal
  copyfloat qTotal,q,4
  end function


'http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToMatrix/

  function MatFromQuat(float *r,*q)
  =================================
  indexbase 0
  float qx=q[0], qy=q[1], qz=q[2], qw=q[3] 'assume already normalised
  '
  float qxx=2*qx*qx, qyy=2*qy*qy, qzz=2*qz*qz
  float qxw=2*qx*qw, qyw=2*qy*qw, qzw=2*qz*qw
  float qxy=2*qx*qy, qxz=2*qx*qz, qyz=2*qy*qz
  '
  '----------------------------------------
  ' 1-qyy-qzz ,   qxy-qzw ,   qxz+qyw  , 0
  '----------------------------------------
  '   qxy+qzw , 1-qxx-qzz ,   qyz-qxw  , 0
  '----------------------------------------
  '   qxz-qyw ,   qyz+qxw , 1-qxx-qyy  , 0
  '----------------------------------------
  '   0       ,   0       ,   0        , 1
  '----------------------------------------
  '
  '
  r[0]  = 1 - qyy - qzz
  r[1]  = qxy + qzw
  r[2]  = qxz - qyw
  r[3]  = 0
  r[4]  = qxy - qzw
  r[5]  = 1 - qxx - qzz
  r[6]  = qyz + qxw
  r[7]  = 0
  r[8]  = qxz + qyw
  r[9]  = qyz - qxw
  r[10] = 1 - qxx - qyy
  r[11] = 0
  r[12] = 0
  r[13] = 0
  r[14] = 0
  r[15] = 1
  '
  'r[0]  = 1 - 2*qy*qy - 2*qz*qz
  'r[1]  = 2*qx*qy + 2*qz*qw
  'r[2]  = 2*qx*qz - 2*qy*qw
  'r[3]  = 0
  'r[4]  = 2*qx*qy - 2*qz*qw
  'r[5]  = 1 - 2*qx*qx - 2*qz*qz
  'r[6]  = 2*qy*qz + 2*qx*qw
  'r[7]  = 0
  'r[8]  = 2*qx*qz + 2*qy*qw
  'r[9]  = 2*qy*qz - 2*qx*qw
  'r[10] = 1 - 2*qx*qx - 2*qy*qy
  'r[11] = 0
  'r[12] = 0
  'r[13] = 0
  'r[14] = 0
  'r[15] = 1
  end function

  function RotateAxis(float a,x,y,z)
  ==================================
  float q[4], m[16]
  SetQuaternion q,a, x,y,z
  MatFromQuat   m,q
  MulMatrix     m
  end function

========================
'HALF-FLOAT / HALF FLOAT
========================
'
'see IEEE 754
'https://en.wikipedia.org/wiki/IEEE_754-1985
'
'TYPE     SIZE   SIGN  EXPONENT  FRACTION
'========================================
'HALF     16     1     5         10
'SINGLE   32     1     8         23
'DOUBLE   64     1     11        52
'EXTENDED 80     1     15        63


function FloatToHalfFloat(sys pf,ph,n) as sys
=============================================
'
'USE OF REGISTERS:
'edx sign transform
'eax exponent transform
'ecx fraction transform
'rsi source pointer
'rdi dest pointer
'
mov rsi,pf
mov rdi,ph
(
  dec dword n 'iterator down count
  jl exit
  mov eax,[rsi]
  mov edx,eax
  and edx,0x80000000 'hold sign bit only
  shr edx,0x10       'shift sign bit down 16
  and eax,0x7fffffff 'remove sign bit
  mov ecx,eax        'for significand
  '
  '
  'TEST FOR ZERO
  (
   cmp eax,0
   jz fwd nzero
  )
  'TEST FOR NAN
  (
   cmp eax,0x7f800000
   jle exit
   mov eax, 0x7fff 'set NAN
   jmp fwd nzero
  )
  'TEST FOR INFINITY
  (
   jl exit
   mov eax, 0x7c00 'set infinity
   jmp fwd nzero
  )
  shr eax,0x17       'shift exponent down 23
  shr ecx,0x0d       'reduce fraction 23 bits to 10 bits
  sub eax,0x70       'adjust exponent bias -112 == (15-127)
  '
  'TEST NEGATIVE EXPONENT BIAS (LOSING PRECISION)
  (
    jge exit     'exclude zero or positive bias
    cmp eax,-10
    (
      jg exit
      mov eax,0 'SET ZERO
      jmp fwd nzero
    )
    xchg ecx,eax
    neg ecx      'make positive
    shr eax,cl   'downshift fraction
    mov ecx,eax
    mov eax,0    'zero exponent
  )   
  '
  'TEST EXPONENT FOR OVERFLOW
  (
    cmp eax,0x1f
    jle exit
    mov eax,0x7c00 'CLAMP INFINITY
    jmp fwd nzero
  )  
  shl eax,0x0a       'place exponent 10 bits up
  and ecx,0x3ff      'mask significand 10 bits
  or  eax,ecx        'combine exponent and significand
  nzero:
  or  eax,edx        'combine sign
  mov [rdi],ax       'store
  add rsi,4          'stride next float
  add rdi,2          'stride next half-float
  repeat
)
return ph
end function


function HalfFloatToFloat(sys ph,pf,n) as sys
=============================================
'
'USE OF REGISTERS:
'edx sign transform
'eax exponent transform
'ecx significand transform
'rsi source pointer
'rdi dest pointer
'
mov rsi,ph
mov rdi,pf
(
  dec dword n 'iterator down count
  jl exit
  xor eax,eax
  mov ax,[rsi]
  mov edx,eax
  and edx,0x8000     'hold sign bit only
  shl edx,0x10       'shift sign bit up 16
  and eax,0x7fff     'remove sign bit
  mov ecx,eax        'for significand
  '
  '
  'TEST FOR ZERO
  (
   cmp eax,0
   jz fwd nzero
  )
  'TEST FOR NAN
  (
   cmp eax,0x7c00
   jle exit
   mov eax, 0x7fffffff 'set NAN
   jmp fwd nzero
  )
  'TEST FOR INFINITY
  (
   jl exit
   mov eax, 0x7f800000 'set infinity
   jmp fwd nzero
  )
  shr eax,0x0A       'shift exponent down 10
  add eax,0x70       'adjust exponent bias +112 == (127-15)
  shl eax,23         'shift exponent into final position
  and ecx,0x3ff      'mask significand 10 bits
  shl ecx,0x0d       'shift significand from 10 bits to 23 bits (13)
  '
  or  eax,ecx        'combine exponent and significand
  nzero:
  or  eax,edx        'combine sign
  mov [rdi],eax      'store
  add rsi,2          'stride next half-float
  add rdi,4          'stride next float
  repeat
)
return ph
end function



=================
class Transform4D
=================
'5x5 matrix
indexbase 0
float f[500] '5x5x20
int idx
method constructor()
  idx=0
end method
method destructor()
end method
method clear()
  int i=idx
  for i=idx to idx+24
    f[i]=0.0
  next
end method
method identity()
  clear()
  float g at @f+(idx*4)
  g[ 0]=1.0
  g[ 6]=1.0
  g[12]=1.0
  g[18]=1.0
  g[24]=1.0
end method
method pushmatrix()
  if idx<=475
    idx+=25
    copy @f+(idx*4),@f+(idx*4)-100,100 '25 floats
  endif
end method
'
method popmatrix()
  if idx>=25
    idx-=25
  endif
end method
'
method mulmatrix()
  'assembler version
  sys p=@f+(idx*4)-200 '4x25x2
  mov rcx,p
  mov rsi,rcx
  mov rdi,rcx 'matrix R
  add rsi,100 'matrix A
  add rdi,200 'matrix B
  jmp fwd begin
  ----
  row:
  ----
  fld  dword [esi   ] 'A col0
  fmul dword [edi   ] 'B row0
  '
  fld  dword [esi+20] 'A col1
  fmul dword [edi+ 4] 'B row1
  faddp
  fld  dword [esi+40] 'A col2
  fmul dword [edi+ 8] 'B row2
  faddp
  fld  dword [esi+60] 'A col3
  fmul dword [edi+12] 'B row3
  faddp
  fld  dword [esi+80] 'A col4
  fmul dword [edi+16] 'B row4
  faddp
  fstp dword [ecx]
  add ecx,4 'R next row
  add esi,4 'A next row
  ret
  ------
  block:
  ------
  call row
  call row
  call row
  call row
  call row
  sub esi,20 'A prev col
  add edi,20 'B next col
  ret
  ------
  begin:
  ------
  call block
  call block
  call block
  call block
  call block
  idx-=50 '-25x2 'drop matrix A & B
  end method
  '
  method translate(float x,y,z,w)
  pushmatrix()
  idx+=25
  identity()
  float g at @f+(idx*4)
  g[20]=x
  g[21]=y
  g[22]=z
  g[23]=w
  mulmatrix()
end method
'
method scale(float x,y,z,w)
  pushmatrix()
  idx+=25
  clear()
  float g at @f+(idx*4)
  g[0]=x
  g[6]=y
  g[12]=z
  g[18]=w
  g[24]=1.0
  mulmatrix()
end method
'
method rotate(float a,int m)
  'm rotation planes
  pushmatrix()
  idx+=25
  identity()
  a=rad(a)
  float sa,ca
  fld dword  a
  fsincos
  fstp dword ca
  fstp dword sa
  '
  float g at @f+(idx*4)
  '
  select m
  case 1 'xy
    g[0]=ca   ' cz00
    g[6]=ca   ' sc00
    g[1]=sa   ' 0010
    g[5]=-sa  ' 0001
  case 2 'xz
    g[0]=ca   ' c0z0
    g[12]=ca  ' 0100
    g[2]=sa   ' s0c0
    g[10]=-sa ' 0001
  case 3 'yz
    g[6]=ca   ' 1000
    g[12]=ca  ' 0cz0
    g[7]=sa   ' 0sc0
    g[11]=-sa ' 0001
  case 4 'xw
    g[0]=ca   ' c00z
    g[18]=ca  ' 0100
    g[3]=sa   ' 0010
    g[15]=-sa ' s00c
  case 5 'yw
    g[6]=ca   ' 1000
    g[18]=ca  ' 0c0z
    g[8]=sa   ' 0010
    g[16]=-sa ' 0s0c
  case 6 'zw
    g[12]=ca  ' 1000
    g[18]=ca  ' 0100
    g[13]=sa  ' 00cz
    g[17]=-sa ' 00sc
  end select
  mulmatrix()
end method
'
method TransformVec3(float *r,*v, int n=1,st=12,oo=0)
  '
  float m at @f+(idx*100)
  int i,j
  @v+=oo 'OFFSET BYTES
  for i=0 to n-1
    for j=0 to 2
      r = m[0]*v[0] + m[5]*v[1] +
          m[10]*v[2] +
          m[20]
      @r+=4
      @m+=4
    next
    @v+=st-12 'st stride bytes
    @r+=st-12
    @m-=12
  next
end method
'
method TransformVec4(float *r,*v, int n=1,st=16,oo=0)
  '
  float m at @f+(idx*100)
  int i,j
  @v+=oo 'OFFSET BYTES
  for i=0 to n-1
    for j=0 to 3
      r = m[0]*v[0] + m[5]*v[1] +
          m[10]*v[2] + m[15]*v[3] +
          m[20]
      @r+=4
      @m+=4
    next
    @v+=st-16 'st stride bytes
    @r+=st-16
    @m-=16
  next
end method
'
method show(string he)
  print he cr
  int x,y,i
  float g at @f+idx*4
  for y=0 to 4
    for x=0 to 4
      i=y+x*5
      print tab str(g[i],4)
    next
    print cr
  next
  print cr
end method
'
method ShowVec3(float *v, string he="", int n=1,st=12,oo=0)
  if he
    print he cr
  endif
  int i
  @v+=oo 'OFFSET BYTES
  for i=0 to n-1
    int x
    float g at @v
    for x=0 to 2
      print tab str(g[x],4)
    next
    print cr
    @v+=st-12 'st stride bytes
  next
end method
'
method ShowVec4(float *v, string he="", int n=1,st=16,oo=0)
  if he
    print he cr
  endif
  int i
  @v+=oo 'OFFSET BYTES
  for i=0 to n-1
    int x
    float g at @v
    for x=0 to 3
      print tab str(g[x],4)
    next
    print cr
    @v+=st-16 'st stride bytes
  next
end method
'
end class 'transform4d

/*
'
'TESTS
======
'
print "matrix 5x5 tests" cr cr
new transform4d m
m.clear                  : m.show "empty"
m.identity               : m.show "identity"
m.pushmatrix             : m.show "pushmatrix"
m.popmatrix              : m.show "popmatrix"
m.pushmatrix
m.pushmatrix
m.mulmatrix()            : m.show "mulmatrix"
m.scale(1,2,3,4)         : m.show "scale"
m.translate(10,20,30,40) : m.show "translate"
m.identity
m.rotate(30,1)           : m.show "rotate xy"
m.identity
m.rotate(30,2)           : m.show "rotate xz"
m.identity
m.rotate(30,3)           : m.show "rotate yz"
m.identity
m.rotate(30,4)           : m.show "rotate xw"
m.identity
m.rotate(30,5)           : m.show "rotate yw"
m.identity
m.rotate(30,6)           : m.show "rotate zw"
m.identity
m.rotate(-30,1)           : m.show "rotate xy"
'm.scale(1,2,3,4)         : m.show "scale"
'm.translate(10,20,30,40) : m.show "translate"
'
type vec4 float x,y,z,w
vec4 vi={1,0,0,0}
vec4 vo
m.transformVec4(vo,vi)    : m.showvec4(vi,"vector4 trans")
                          : m.showvec4(vo,"vector4 trans")
                          : m.showvec3(vo,"vector3 trans")
*/

'#ifdef MathUtilExtra

