/*-------------------------------------------------------------

Copyright (C) 2000 Peter Clote. 
All Rights Reserved.

Permission to use, copy, modify, and distribute this
software and its documentation for NON-COMMERCIAL purposes
and without fee is hereby granted provided that this
copyright notice appears in all copies.


THE AUTHOR AND PUBLISHER MAKE NO REPRESENTATIONS OR
WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS
AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
THIS SOFTWARE OR ITS DERIVATIVES.

-------------------------------------------------------------*/

/*
   It is an open problem to characterize the distribution of the root
   mean square deviation of normed distance matrix of a random
   self-avoiding walk from normed distance matrix of a fixed protein as
   given by PDB data. (Self-avoiding walks on face centered cubic lattice
   in 3 dimensions.) This distribution resembles the easier problem of
   generating 2 random walks in 2 dimensions on a square lattice, where
   each coordinate increment DeltaX and DeltaY are chosen uniformly from
   0,1,-1. This program generates such walks. 
 
   output for each pair of walks: 
   mean square distance and root mean square distance
*/



#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int randint(int);

main() {
  int i,j;
  int M = 1000; // number of runs
  int n = 50; // 20 points in walk
  int x,y,u,v; // (x,y) and (u,v) are the points of the two walks
  double distSq; // sum of squares of distance

  // generate M times 2 walks
  for (i=0;i<M;i++){
    x=y=u=v=0; //start from origin
    distSq=0;
    // each walk has length n
    for (j=0;j<n;j++){
       switch( randint(3) ) {
         case 0: x += 1; break;
         case 1: x -= 1; break;
         default: break;  // x = x
       } 
       switch( randint(3) ) {
         case 0: y += 1; break;
         case 1: y -= 1; break;
         default: break;  // y = y
       } 
       switch( randint(3) ) {
         case 0: u += 1; break;
         case 1: u -= 1; break;
         default: break;  // u = u
       } 
       switch( randint(3) ) {
         case 0: v += 1; break;
         case 1: v -= 1; break;
         default: break;  // v = v
       } 
       distSq += (x-u)*(x-u) + (y-v)*(y-v);
       }
    printf("%lf\t%lf\n",distSq/n,sqrt(distSq/n));
  }
}

// returns uniform random in [0,1]
double rreal() {
  return ((double)rand())  / (RAND_MAX+1.0);
}


// returns uniform random int in {0,...,m-1}
int randint(int m) {
  return (int) (rreal()*m);
}