more ~/poisson_1d_c.txt
#include <stdlib.h> #include <stdio.h> #include <math.h> ////////////////////////////////////////////////// // algorithm poisson random number (Knuth): // // init: // // Let L ← exp(-lambda), k ← 0 and p ← 1. // // do: // // k ← k + 1. // // Generate uniform random number u in [0,1] // // and let p ← p × u. // // while p ≥ L. // // return k − 1. // ////////////////////////////////////////////////// int poisson ( double lmbd ) { double L = exp( -lmbd ); int k = 0; double p = 1; do { double u = drand48(); p *= u; k += 1; } while ( p >= L ); return( k-1 ); } int main ( int argc, char** argv ) { int i; float s = atof( argv[1] ); // initial seed // int p = atoi( argv[2] ); // sampling points // int z = atoi( argv[3] ); // total size // int *v; // vector of sampling points // float ld = (float) z / (float) p; float w = 2.0; // inital weight // int k = 0; // currently found sampling points // v = ( int* ) malloc( z*sizeof( z ) ); srand48( s ); // initialize seed // do { i = 0; k = 0; // N.B. first point always acqured // // we use the nomenclature of first point == 0 // while ( i < z ) { v[k] = i; // store point to be acquired // i += 1; // put pointer at next point // k += 1; // next index of acqured point // // Now make a gap : // i += poisson( (ld-1.0)*w*sin((float)i/(float)(z+1)*M_PI) ); } // if more points T.B. acquired found than than wanted: try again with weight 2% larger // // if fewer points T.B. acquired found than than wanted: try again with weight 2% smaller // if ( k > p ) w *= 1.02; if ( k < p ) w /= 1.02; } while ( k != p ); // try until correct number points to be acquired found // // print the data on standard output // for ( k = 0 ; k < p ; k++ ) printf( "%d\n", v[k] ); }