/*************************************************************************** Hyperspherical distribution                          main.c  -  description                             -------------------    begin                : 2003-07-04    last update      : 2003-07-04    copyright            : (C) 2003 by     email                : Maurice.Clerc@WriteMe.com ***************************************************************************//*************************************************************************** *                                                                         * *   This program is free software; you can redistribute it and/or modify  * *   it under the terms of the GNU General Public License as published by  * *   the Free Software Foundation; either version 2 of the License, or     * *   (at your option) any later version.                                   * *                                                                         * ***************************************************************************/ /* It generates an uniform distribution inside a D-sphere (radius 1, center: origin). The method itself is proved. However, there may be some bugs. */                                                 #include <stdio.h>#include <stdlib.h>#include <math.h>#include <string.h>#include <syscall.h>#define Max_D 			100 // Max number of dimensions// Structuresstruct vector		{int size;double x[Max_D];};// Subroutinesdouble 			alea_float(double min, double max);double  alea_normal(double mean,double std_dev);void display_vector(struct vector x);void  save_vector(FILE *f_run,struct vector x);// Global variablesdouble   pi;//=============================================================================== int main(int argc, char *argv[]){ FILE *f_run; int                 D; int                 i; int                 j; double           length; int                 N; double           pw; double           r; struct  vector   x; pi=2*acos(0);  f_run =fopen("run.txt","w");  // File to save result     //---------------------------------------   Interactive input                            input:  printf("\n Dimension?  (0=stop):  ");  scanf("%i",&D);  if(D==0) goto end;  if (D>Max_D )  {     printf("\n No more than %i dimensions, please.",Max_D);     goto input;  }   printf("\n Number of points to generate? (0=stop): "); scanf("%i",&N); if (N==0) goto end; x.size=D; pw=1/(double)D; //------------------------------------------------------------------------------------------- LOOP on points to generate // This just to see how "looks like" the distribution.  Normally, in real applications, // you generate just _one_  pointfor (i=0;i<N;i++){// ----------------------------------- Step 1.  Direction    length=0;   for (j=0;j<D;j++)   {          x.x[j]=alea_normal(0,1);          length=length+  x.x[j]*x.x[j];   }   length=sqrt(length); //----------------------------------- Step 2.   Random radius  r=alea_float(0,1);  r=pow(r,pw);       for (j=0;j<D;j++)   {          x.x[j]=r*x.x[j]/length;   }  //-------------------------------------------------------- display_vector(x); save_vector(f_run,x); } //------------------------------------------------------------------------ End of loop on points  goto input;   // Ask for a new set of points end: printf("\n %i generated points saved on run.txt" ,N); return EXIT_SUCCESS;}// ========================================================================= ALEA_FLOATdouble alea_float(double min, double max) /* Random  number between min and max */{double 	r;// Normally, RAND_MAX = 32767 = 2^31-1	r=(double)rand()/RAND_MAX; // r in [0,1]	r=min+r*(max-min);	return r;}//==========================================================================  ALEA_NORMALdouble alea_normal(double mean, double std_dev) { /*     Use the polar form of the Box-Muller transformation to obtain     a pseudo random number from a Gaussian distribution */        double x1, x2, w, y1;        //double  y2;         do {                 x1 = 2.0 * alea_float(0,1) - 1.0;                 x2 = 2.0 * alea_float(0,1) - 1.0;                 w = x1 * x1 + x2 * x2;         } while ( w >= 1.0);         w = sqrt( (-2.0 * log( w ) ) / w );         y1 = x1 * w;        // y2 = x2 * w;          y1=y1*std_dev+mean;         return y1;  }//========================================================================== DISPLAY_VECTORvoid display_vector(struct vector x){int				d;int				ncar=0;//printf("\n Vector size %i\n ",x.size);printf("\n");for (d=0;d<x.size;d++){	printf("%6.6f ",x.x[d]);	ncar=ncar+10;	if (ncar>75)	{		ncar=0;		printf("\n");	}}} //========================================================================== SAVE_VECTORvoid  save_vector(FILE *f_run,struct vector x){int				d;int				ncar=0;fprintf(f_run,"\n");for (d=0;d<x.size;d++){	fprintf(f_run,"%6.6f ",x.x[d]);	ncar=ncar+10;	if (ncar>75)	{		ncar=0;		fprintf(f_run,"\n");	}}}