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

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