valarrayを使ってみようかと

数値演算とか行列向きの低レベルなコンテナ(STLではない)らしいvalarrayを使ってみたかったので使ってみた。
試しに浮動小数点の頂点配列を想定した配列に行列をかけるのと行列同士の積をやってみた。
要するに内積。
いつの間にかvalarrayの嬉しさの中心たるindirect_arrayとかsliceをまったく使わないコードになってしまったのだが、その過程でむしろinner_product()を発見。こっちは使えそうかなぁ。
ちなみにinner_productの第４引数はtemplateの型判別に使われるので浮動少数なら浮動少数で
きっちり0.0と中身の型が判るように書く必要があるのに注意。

#include <iostream>
#include <valarray>
#include <algorithm> // copy
#include <numeric> // inner_product

typedef double NUM_TYPE;
typedef std::valarray<NUM_TYPE> VECTOR;

// 頂点配列
NUM_TYPE data[]={
  -1, -1, -1,  1,
  1, -1, -1,  1,
  1,  1, -1,  1,
  -1,  1, -1,  1,
  -1, -1,  1,  1,
  1, -1,  1,  1,
  1,  1,  1,  1,
  -1,  1,  1,  1,
};
const size_t VERTICES_SIZE=24;

void fill_vertices(NUM_TYPE *v)
{
  std::copy(data, data+VERTICES_SIZE, v);
}

void print_vertex(NUM_TYPE *v)
{
  std::cout << '{' << v[0] << ',' << v[1] << ',' << v[2] << ',' << v[3] << '}';
}

void print_vertices(VECTOR &v){
  int line=0;
  for(size_t i=0; i<v.size(); i+=4){
    std::cout << '[' << line++ << ']';
    print_vertex(&v[i]);
    std::cout << std::endl;
  }
}

void matrix_create_scaling(VECTOR &m, NUM_TYPE x, NUM_TYPE y, NUM_TYPE z)
{
  m=0; m[0]=x; m[5]=y; m[10]=z; m[15]=1;
}

// 頂点配列に行列を適用する
void matrix_apply_vertices(VECTOR &m, VECTOR &v)
{
  VECTOR current(4);
  //double *current;
  for(size_t i=0; i<v.size(); i++)
  {
    if(i % 4 ==0){
      current=v[std::slice(i, 4, 1)];
      //current=&v[i];
    }
    //v[i]=(current * m[std::slice((i % 4)*4, 4, 1)]).sum();
    // inner_productを使うように変更
    v[i]=std::inner_product(&current[0], &current[0]+4, &m[(i%4)*4], 0.0); // 0.0 . not 0
    // sliceをやめてポインタを直接使うのは副作用で結果が変わるのでやっぱだめｗ
    //v[i]=std::inner_product(current, current+4, &m[(i%4)*4], 0.0); // 0.0 . not 0

  };
}

// 転置
VECTOR matrix_create_transpose(const VECTOR &m)
{
  VECTOR result(16);
  for(size_t i=0; i<4; i++)
  {
    result[i]=m[4*i];
    result[i+4]=m[4*i+1];
    result[i+8]=m[4*i+2];
    result[i+12]=m[4*i+3];
  }
  return result;
}

VECTOR matrix_multiply(VECTOR &lhs, VECTOR &rhs)
{
  VECTOR result(16);
  /*
     for(size_t i=0; i<16; i+=4)
     {
     result[i]=(lhs[std::slice(i, 4, 1)] * rhs[std::slice(0, 4, 4)]).sum();
     result[i+1]=(lhs[std::slice(i, 4, 1)] * rhs[std::slice(1, 4, 4)]).sum();
     result[i+2]=(lhs[std::slice(i, 4, 1)] * rhs[std::slice(2, 4, 4)]).sum();
     result[i+3]=(lhs[std::slice(i, 4, 1)] * rhs[std::slice(3, 4, 4)]).sum();
     }
     */
  // inner_productを使うように変更
  double *current;
  VECTOR transposed=matrix_create_transpose(rhs);
  for(size_t i=0; i<16; i++)
  {
    if(i%4==0)
      current=&lhs[i];
    result[i]=std::inner_product(current, current+4, &transposed[(i%4)*4], 0.0);
  }
  return result;
}

int main(int argc, char **argv)
{
  std::cout << "[store vertex]" << std::endl;
  VECTOR vertices(VERTICES_SIZE);
  fill_vertices(&vertices[0]);
  print_vertices(vertices);

  std::cout << "[scaling matrix]" << std::endl;
  VECTOR scaling(16); // 4x4matrix
  matrix_create_scaling(scaling, 1, 2, 3);
  print_vertices(scaling);

  std::cout << "[apply matrix to vertices]" << std::endl;
  matrix_apply_vertices(scaling, vertices);
  print_vertices(vertices);

  std::cout << "[muliply matrix]" << std::endl;
  std::cout << "[scaling]" << std::endl;
  print_vertices(scaling);

  // 平行移動行列
  VECTOR translation(16);
  translation=0;
  translation[std::slice(0, 4, 5)]=1;
  translation[3]=1;
  translation[7]=2;
  translation[11]=3;
  std::cout << "[translation]" << std::endl;
  print_vertices(translation);

  std::cout << "[multiply]" << std::endl;
  print_vertices(matrix_multiply(translation, scaling));
};

valarrayの参考

http://smdn.invisiblefulmoon.net/ikimasshoy/cpp/valarray1.html
http://smdn.invisiblefulmoon.net/ikimasshoy/cpp/valarray2.html