#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-1; ++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 scale = 0x1.p-30;
r rnd(){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);}