Re: [Help-gsl] Problems with #include options in GSL

[Joseph Wakeling]

Apologies to Martin who will be getting this mail twice.  The first time 
I replied, I accidentally sent it to him only...


Martin Jansche wrote:
> > Note that I deliberately *don't* want to just compile to object form the
> > whole GSL.
>
> Why not?

Well, we could consider that to be merely part of my inexperience at 
working with large-scale code projects. ;-)

More seriously, since my projects are relatively small-scale, it's nice 
to be able to take only a few functions from GSL since most of them I 
don't right now need.  Then if I was going to pass on that code to a 
friend, I needn't require him to install the whole GSL, I can just pass 
on appropriate directories with the functions.

As it is I decided that I might take the time to rewrite a few of the 
modules as standalone pieces.  I've done this with the ran1 random 
number generator so as to be able to easily integrate it into my own code.

Here's what I wrote.  Does the notice put at the beginning look OK?  I 
was a bit unsure about the copyright thing but I guess I should take 
some responsibility for the changes.  I'm kind of new to software 
licensing issues so I want to make sure I'm not treading on any toes.

Any comments welcome.  Note that I *deliberately* removed the passing of 
the random seed structure because in my own code it's faster and the 
code uses one random number generator "universally" throughout the 
program.  Though I would definitely appreciate thoughts on the pros and 
cons of this.

N.B. is it really necessary to have the various static const variables 
declared at the start of the module, or would it be better to just use 
#define for them?

        -- Joe



/* gsl-_ran1.c
*
* This progam is a modified version of the rng/ran1.c code in the GNU
* Scientific Library (GSL), version 1.7.  It will stand alone as a
* module and does not need any code other than is contained herein!
*                                      -- Joseph Rushton Wakeling, 2005
*
* Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
* Copyright (C) 2005 Joseph Rushton Wakeling
*
* 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.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 
02110-1301, USA.
*/


#define N_SHUFFLE 32
#define N_DIV (1 + 2147483646/N_SHUFFLE)

static const long int m = 2147483647, a = 16807, q = 127773, r = 2836;
static const double z_max = 1 - 1.2e-7f;

static unsigned long int x;
static unsigned long int n;
static unsigned long int shuffle[N_SHUFFLE];


static unsigned long int ran1_get (void) {
  const long int h = x / q;
  const long int t = a * (x - h * q) - h * r;

  if(t<0)
      x = t + m;
  else
      x = t;

  {
      unsigned long int j = n / N_DIV;
      n = shuffle[j];
      shuffle[j] = x;
  }

  return(n);
}


double ran1_get_double (void) {
  double z = ran1_get() / 2147483647.0f ;
  if(z > z_max)
      return(z_max);
  return(z);
}


void ran1_set (unsigned long int s) {
  int i;
  long int h, t;

  if(s==0)
      s = 1;
     for(i=0;i<8;++i) {
      h = s / q;
      t = a * (s - h * q) - h * r;
      if(t < 0)
          t += m;
      s = t;
  }

  for(i=N_SHUFFLE-1;i>=0;--i) {
      h = s/q;
      t = a * (s - h * q) - h * r;
      if(t < 0)
          t += m;
      s = t;
      shuffle[i] = s;
  }

  x = s;
  n = s;

  return;
}