htswanalysis/MACS/lib/gsl/gsl-1.11/qrng/halton.c

   1 /* Authors: O. Teytaud
   2  * Copyright (C) 2007  O. Teytaud
   3  * (all comments welcome at olivier.teytaud@inria.fr)
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Implementation for Halton generator.  See [J.H. Halton, On the
  21  * efficiency of certain quasi-random sequences of points in
  22  * evaluating multi-dimensional integrals Numerische Mathematik, 1960]
  23  */
  24
  25 #include <config.h>
  26 #include <gsl/gsl_qrng.h>
  27
  28 /* maximum allowed space dimension */
  29 #define HALTON_MAX_DIMENSION 1229
  30
  31 /* prototypes for generator type functions */
  32 static size_t halton_state_size (unsigned int dimension);
  33 static int halton_init (void *state, unsigned int dimension);
  34 static int halton_get (void *state, unsigned int dimension, double *v);
  35
  36 /* global Halton generator type object */
  37 static const gsl_qrng_type halton_type = {
  38   "halton",
  39   HALTON_MAX_DIMENSION,
  40   halton_state_size,
  41   halton_init,
  42   halton_get
  43 };
  44
  45 const gsl_qrng_type *gsl_qrng_halton = &halton_type;
  46
  47 /* prime numbers (thanks to trolltech http://doc.trolltech.com/3.2/primes.html)
  48 */
  49 static const int prime_numbers[HALTON_MAX_DIMENSION] = {
  50   2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
  51   31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
  52   73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
  53   127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
  54   179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
  55   233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
  56   283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
  57   353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
  58   419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
  59   467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
  60   547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
  61   607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
  62   661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
  63   739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
  64   811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
  65   877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
  66   947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013,
  67   1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069,
  68   1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151,
  69   1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,
  70   1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291,
  71   1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373,
  72   1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451,
  73   1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
  74   1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583,
  75   1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657,
  76   1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733,
  77   1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,
  78   1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889,
  79   1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987,
  80   1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053,
  81   2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,
  82   2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213,
  83   2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287,
  84   2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
  85   2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
  86   2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531,
  87   2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617,
  88   2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687,
  89   2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,
  90   2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819,
  91   2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903,
  92   2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999,
  93   3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
  94   3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181,
  95   3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257,
  96   3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331,
  97   3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
  98   3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511,
  99   3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571,
 100   3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643,
 101   3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,
 102   3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821,
 103   3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907,
 104   3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989,
 105   4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057,
 106   4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139,
 107   4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231,
 108   4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
 109   4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409,
 110   4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493,
 111   4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583,
 112   4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657,
 113   4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751,
 114   4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831,
 115   4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937,
 116   4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003,
 117   5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087,
 118   5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179,
 119   5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279,
 120   5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387,
 121   5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443,
 122   5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521,
 123   5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639,
 124   5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693,
 125   5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791,
 126   5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857,
 127   5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939,
 128   5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053,
 129   6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133,
 130   6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221,
 131   6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301,
 132   6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367,
 133   6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473,
 134   6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571,
 135   6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673,
 136   6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761,
 137   6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833,
 138   6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917,
 139   6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997,
 140   7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103,
 141   7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207,
 142   7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297,
 143   7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411,
 144   7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499,
 145   7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561,
 146   7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643,
 147   7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723,
 148   7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829,
 149   7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919,
 150   7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017,
 151   8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111,
 152   8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219,
 153   8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291,
 154   8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387,
 155   8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501,
 156   8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597,
 157   8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677,
 158   8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741,
 159   8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831,
 160   8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929,
 161   8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011,
 162   9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109,
 163   9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199,
 164   9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283,
 165   9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377,
 166   9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439,
 167   9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533,
 168   9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631,
 169   9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733,
 170   9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811,
 171   9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887,
 172   9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973
 173 };
 174
 175 /* Halton generator state.
 176  *   sequence_count       = number of calls with this generator
 177  */
 178 typedef struct
 179 {
 180   unsigned int sequence_count;
 181 }
 182 halton_state_t;
 183
 184 static size_t
 185 halton_state_size (unsigned int dimension)
 186 {
 187   return sizeof (halton_state_t);
 188 }
 189
 190 static int
 191 halton_init (void *state, unsigned int dimension)
 192 {
 193   halton_state_t *h_state = (halton_state_t *) state;
 194
 195   h_state->sequence_count = 0;
 196
 197   if (dimension < 1 || dimension > HALTON_MAX_DIMENSION)
 198     {
 199       return GSL_EINVAL;
 200     }
 201
 202   return GSL_SUCCESS;
 203 }
 204
 205 static double
 206 vdcorput (int x, int b)
 207 {
 208   double r = 0.;
 209   double v = 1.;
 210   double binv = 1. / (double) b;
 211
 212   while (x > 0)
 213     {
 214       v *= binv;
 215       r += v * (double) (x % b);
 216       x /= b;
 217     }
 218   return r;
 219 }
 220
 221 static int
 222 halton_get (void *state, unsigned int dimension, double *v)
 223 {
 224   halton_state_t *h_state = (halton_state_t *) state;
 225   unsigned int i;
 226
 227   if (dimension < 1 || dimension > HALTON_MAX_DIMENSION)
 228     {
 229       return GSL_EINVAL;
 230     }
 231   h_state->sequence_count++;
 232
 233   for (i = 0; i < dimension; i++)
 234     {
 235       v[i] = vdcorput (h_state->sequence_count, prime_numbers[i]);
 236     }
 237
 238   return GSL_SUCCESS;
 239 }