binaries/src/disembl/Tisean_3.0.1/source_c/routines/rand_arb_dist.c

   1 /*
   2  *   This file is part of TISEAN
   3  *
   4  *   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   5  *
   6  *   TISEAN is free software; you can redistribute it and/or modify
   7  *   it under the terms of the GNU General Public License as published by
   8  *   the Free Software Foundation; either version 2 of the License, or
   9  *   (at your option) any later version.
  10  *
  11  *   TISEAN is distributed in the hope that it will be useful,
  12  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  *   GNU General Public License for more details.
  15  *
  16  *   You should have received a copy of the GNU General Public License
  17  *   along with TISEAN; if not, write to the Free Software
  18  *   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  19  */
  20 /*Author: Rainer Hegger Last modified: Feb 11, 2006 */
  21 /*Changes:
  22   Feb 11, 2006: First version
  23 */
  24 /*Comment:
  25   Creates a sequence of random numbers with arbitrary distribution
  26   Input Paramters: x=original data defining the distribution
  27                    nx= length of x
  28                    nc= number of random numbers to create
  29                    nb= number of bins for the distribution
  30                    iseed = seed for the random number generator
  31   Return: rand = array containing the nc random numbers
  32 */
  33
  34 #ifndef _STDLIB_H
  35 #include <stdlib.h>
  36 #endif
  37
  38 #ifndef _LIMITS_H
  39 #include <limits.h>
  40 #endif
  41
  42 #ifndef _TIME_H
  43 #include <time.h>
  44 #endif
  45
  46 extern void rescale_data(double*,unsigned long,double*,double*);
  47 extern void check_alloc(void*);
  48 extern unsigned long rnd_long(void);
  49 extern void  rnd_init(unsigned long);
  50
  51 double *rand_arb_dist(double *x,unsigned long nx,unsigned long nc,
  52                       unsigned int nb,unsigned long iseed)
  53 {
  54   double h,min,inter,*randarb,drnd,epsinv=1.0/(double)nb;
  55   unsigned long i,j,*box,hrnd,nall=nx+nb;
  56
  57   rescale_data(x,nx,&min,&inter);
  58
  59   check_alloc(box=(unsigned long*)malloc(sizeof(unsigned long)*nb));
  60   for (i=0;i<nb;i++)
  61     box[i]=1;
  62
  63   for (i=0;i<nx;i++) {
  64     h=x[i];
  65     if (h >= 1.0)
  66       h -= epsinv/2.0;
  67     j=(unsigned int)(h*nb);
  68     box[j]++;
  69   }
  70   for (i=1;i<nb;i++)
  71     box[i] += box[i-1];
  72
  73   check_alloc(randarb=(double*)malloc(sizeof(double)*nc));
  74
  75   if (iseed == 0)
  76     iseed=(unsigned long)time((time_t*)&iseed);
  77
  78   rnd_init(iseed);
  79   for (i=0;i<1000;i++)
  80     rnd_long();
  81
  82   for (i=0;i<nc;i++) {
  83     hrnd=rnd_long()%nall;
  84     for (j=0;j<nb;j++)
  85       if (box[j] >= hrnd)
  86         break;
  87     drnd=(double)rnd_long()/(double)ULONG_MAX*epsinv;
  88     randarb[i]=min+((double)j*epsinv+drnd)*inter;
  89   }
  90
  91   free(box);
  92
  93   return randarb;
  94 }