#include <math.h>
#include <stdlib.h>
typedef double R;

// Gram-Schmidt
void gs(int m, int n, R a[m][n]){
int j; for(j=0; j<m; ++j) {int k; for(k=0; k<j; ++k) {
   R f = 0; {int l=n; while(l--) f += a[j][l]*a[k][l];}
   {int i=n; while(i--) a[j][i] -= f*a[k][i];}}
   R s = 0; {int i=n; while(i--) s += a[j][i]*a[j][i];}
   {R m = 1/sqrt(s); {int k=n; while(k--) a[j][k] *= m;}}}}

// Random Normal Deviate
R rnd(){R const scale = 0x1.p-30;
  static int f = 0;
  static R m;
  if(f) {f = 0; return m;}
  while (1){R x = scale*random()-1, y = scale*random()-1,
    s = x*x +y*y;
    if(s<1) {R b = sqrt(- 2*log(s)/s);
    f = 1; m = x*b; return y*b;}}}

// Create Random Orthogonal Matrix
void grom(int m, int n, R a[m][n]){
 {int i, j; for(i=0; i<m; ++i) for(j=0; j<n; ++j) a[i][j] = rnd();} gs(m, n, a);}

// Matrix Multiply
void mm(int i, int j, int k, R a[i][j], R b[j][k], R c[i][k]){
  {int I, J, K; for(I=0; I<i; ++I) for(K=0; K<j; ++K) {
    R x = 0; for(J=0; J<j; ++J) x += a[I][J]*b[J][K];
    c[I][K] = x;}}}


/*(define (ip a b) (let p ((a a) (b b) (s 0)) (if (null? a) s
   (p (cdr a) (cdr b) (+ s (* (car a) (car b)))))))*/