[GRASS-dev] BLAS/LAPACK (Numerical Recipies in C)

[Glynn Clements]

Hamish wrote:

> Sören Gebbert wrote:
> > But the gmath lu solver is code from the numerical recipes,
> > we have to rewrite this method.
> 
> re N.R. code in GRASS, 
>   http://thread.gmane.org/gmane.comp.gis.grass.devel/18764/focus=18791
> 
> (this is just fyi, it is of course a good idea to remove the issue
> altogether by removing any code based on N.R.)
> 
> 
> On 2007-02-26, Stefano wrote:
> > Note that previous GRASS development team (CERL) had
> > an agreement with the authors of NR, to use freely NR:
> 
> [grass5/src.nonGPL/README]
>  /* Based on "Numerical Recipies in C; The Art of Scientific Computing"
>     (Cambridge University Press, 1988).  Copyright (C) 1986, 1988 by
>     Numerical Recipes Software.  Permission is granted for unlimited
>     use within GRASS only.  
>     [...]
>     Bill Press
>   */
> [sic]
> 
> If this agreement was with CERL I assume that means it predates GRASS
> moving to the GPL. The GPL audit alluded to in grass5/TODO.txt and
> grass5/src.nonGPL/ should have caught all known occurrences.
> 
> 
> grass-5.4.0$ grep -rI "Numerical Recipes" * | cut -f1 -d: | uniq
> TODO.txt
> html/grass_commandlist.html
> html/html/r.kineros.html
> html/html/r.surf.random.html
> html/html/r.tribs.html
> src/raster/r.surf.gauss/main.c
> src/raster/r.surf.gauss/README
> src/raster/r.surf.gauss/source.txt
> src/raster/r.surf.random/README
> src/raster/r.surf.random/SOURCE.TXT
> src.contrib/CERL/sites/s.surf.krig/NR.c
> src.contrib/CERL/sites/s.surf.krig/nrutil.c
> src.nonGPL/display/d.param.scale/CHANGES
> src.nonGPL/README
> 
> 
> A quick grep for "Numerical Recipes" in 6.3-cvs finds these references:
>  raster/r.surf.random
>  raster/r.surf.gauss
> 
> Although both modules' main.c ask:  /* still true ?*/

Personally, I think that NRS would have a hard time claiming copyright
over what is essentially a direct translation of the Box-Mueller
algorithm to C.

If you want to be able to cite a GPL source for the code, take the
version from GSL:

http://sources.redhat.com/cgi-bin/cvsweb.cgi/gsl/randist/gauss.c?rev=1.24&content-type=text/x-cvsweb-markup&cvsroot=gsl

	double
	gsl_ran_gaussian (const gsl_rng * r, const double sigma)
	{
	  double x, y, r2;
	
	  do
	    {
	      /* choose x,y in uniform square (-1,-1) to (+1,+1) */
	      x = -1 + 2 * gsl_rng_uniform_pos (r);
	      y = -1 + 2 * gsl_rng_uniform_pos (r);
	
	      /* see if it is in the unit circle */
	      r2 = x * x + y * y;
	    }
	  while (r2 > 1.0 || r2 == 0);
	
	  /* Box-Muller transform */
	  return sigma * y * sqrt (-2.0 * log (r2) / r2);
	}

Note: gsl_rng_uniform_pos() just returns a random double in the range
0-1, so you can use "(double) rand() / RAND_MAX" for that.

-- 
Glynn Clements <glynn at gclements.plus.com>