From 22f703cab05b7cd368f4de9e03991b7664dc5022 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Mon, 1 Sep 2014 13:56:46 +0200 Subject: Initial import of argyll version 1.5.1-8 --- numlib/sobol.c | 211 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 211 insertions(+) create mode 100644 numlib/sobol.c (limited to 'numlib/sobol.c') diff --git a/numlib/sobol.c b/numlib/sobol.c new file mode 100644 index 0000000..40a2c38 --- /dev/null +++ b/numlib/sobol.c @@ -0,0 +1,211 @@ + +/***************************************************/ +/* Sobol sub-random vector sequence generator */ +/***************************************************/ + +/* Code is an expression of the algorithm decsribed in */ +/* the SSOBOL.F fortran source file, with additional */ +/* guidance from "Numerical Recipes in C", by W.H.Press, B.P.Flannery, */ +/* S.A.Teukolsky & W.T.Vetterling. */ + +/* + * Copyright 2002 Graeme W. Gill + * All rights reserved. + * + * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :- + * see the License.txt file for licencing details. + */ + +#include "numsup.h" +#include "sobol.h" + +/* + * The array poly gives successive primitive + * polynomials coded in binary, e.g. + 45 = 100101 + * has bits 5, 2, and 0 set (counting from the + * right) and therefore represents + X**5 + X**2 + X**0 + + * These polynomials are in the order used by + * sobol in ussr comput. maths. math. phys. 16 (1977), + * 236-242. + */ + +static int sobol_poly[SOBOL_MAXDIM] = { + 1, 3, 7, 11, 13, 19, 25, 37, 59, 47, + 61, 55, 41, 67, 97, 91, 109, 103, 115, 131, + 193, 137, 145, 143, 241, 157, 185, 167, 229, 171, + 213, 191, 253, 203, 211, 239, 247, 285, 369, 299 +}; + +/* + * The initialization of the array vinit is from + * Sobol and Levitan, the production of points uniformly + * distributed in a multidimensional cube (in Russian), + * preprint ipm akad. nauk sssr, no. 40, moscow 1976. + * For a polynomial of degree m, m initial + * values are needed : these are the values given here. + * subsequent values are calculated during initialisation. + */ + +static int vinit[8][SOBOL_MAXDIM] = { + { + 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 + }, + { + 0, 0, 1, 3, 1, 3, 1, 3, 3, 1, + 3, 1, 3, 1, 3, 1, 1, 3, 1, 3, + 1, 3, 1, 3, 3, 1, 3, 1, 3, 1, + 3, 1, 1, 3, 1, 3, 1, 3, 1, 3 + }, + { + 0, 0, 0, 7, 5, 1, 3, 3, 7, 5, + 5, 7, 7, 1, 3, 3, 7, 5, 1, 1, + 5, 3, 3, 1, 7, 5, 1, 3, 3, 7, + 5, 1, 1, 5, 7, 7, 5, 1, 3, 3 + }, + { + 0, 0, 0, 0, 0, 1, 7, 9, 13, 11, + 1, 3, 7, 9, 5, 13, 13, 11, 3, 15, + 5, 3, 15, 7, 9, 13, 9, 1, 11, 7, + 5, 15, 1, 15, 11, 5, 3, 1, 7, 9 + }, + { + 0, 0, 0, 0, 0, 0, 0, 9, 3, 27, + 15, 29, 21, 23, 19, 11, 25, 7, 13, 17, + 1, 25, 29, 3, 31, 11, 5, 23, 27, 19, + 21, 5, 1, 17, 13, 7, 15, 9, 31, 9 + }, + { + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 37, 33, 7, 5, 11, 39, 63, + 27, 17, 15, 23, 29, 3, 21, 13, 31, 25, + 9, 49, 33, 19, 29, 11, 19, 27, 15, 25 + }, + { + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 13, + 33, 115, 41, 79, 17, 29, 119, 75, 73, 105, + 7, 59, 65, 21, 3, 113, 61, 89, 45, 107 + }, + { + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 7, 23, 39 + } +}; + +/* Get the next sobol vector */ +/* return nz if we've run out */ +static int next_sobol(sobol *s, double * v) +{ + int i, p; + unsigned int c; + + s->count++; + + /* Find the position of the right-hand zero in count */ + for (c = s->count, p = 0; (c & 1) == 0; p++, c >>= 1) + ; + + if(p > SOBOL_MAXBIT) + return 1; /* Run out */ + + for (i = 0; i < s->dim; i++) { + s->lastq[i] ^= s->dir[p][i]; + v[i] = s->lastq[i] * s->recipd; + } + + return 0; +} + +/* Free up the object */ +static void del_sobol(sobol *s) { + if (s != NULL) + free(s); +} + +/* reset the count */ +static void reset_sobol(sobol *s) { + int i; + + /* Set up first vector and values */ + s->count = 0; + for (i = 0; i < s->dim; i++) + s->lastq[i] = 0; +} + +/* Return NULL on error */ +sobol *new_sobol(int dim) { + sobol *s = NULL; + int i, j, p; + + if (dim < 1 || dim > SOBOL_MAXDIM) { + return NULL; + } + + if ((s = (sobol *)malloc(sizeof(sobol))) == NULL) { + return NULL; + } + + s->dim = dim; + s->next = next_sobol; + s->reset = reset_sobol; + s->del = del_sobol; + + /* Initialize the direction table */ + for (i = 0; i < dim; i++) { + + if (i == 0) { + for (j = 0; j < SOBOL_MAXBIT; j++) + s->dir[j][i] = 1; + } else { + int m; /* Degree */ + int pm; /* Polinomial mask */ + + /* Find degree of polynomial from binary encoding */ + for (m = 0, pm = sobol_poly[i] >> 1; pm != 0; m++, pm >>= 1) + ; + + /* The leading elements of row i come from vinit[][] */ + for (j = 0; j < m; j++) { + s->dir[j][i] = vinit[j][i]; + } + + /* Calculate remaining elements of row i as explained */ + /* in bratley and fox, section 2 */ + pm = sobol_poly[i]; + for (j = m; j < SOBOL_MAXBIT; j++) { + int k; + int newv = s->dir[j-m][i]; + for (k = 0; k < m; k++) { + if (pm & (1 << (m-k-1))) { + newv ^= s->dir[j-k-1][i] << (k+1); + } + } + s->dir[j][i] = newv; + } + } + } + /* Multiply columns of v by appropriate power of 2 */ + for (p = 2, j = SOBOL_MAXBIT-2; j >= 0; j--, p <<= 1) { + for (i = 0; i < dim; i++) + s->dir[j][i] *= p; + } + + /* recipd is 1/(common denominator of the elements in v) */ + s->recipd = 1.0/(1 << SOBOL_MAXBIT); + + /* Set up first vector and values */ + s->count = 0; + for (i = 0; i < dim; i++) + s->lastq[i] = 0; + + return s; +} + -- cgit v1.2.3