test_ssp_ran_uniform_disk.h (3902B)
1 /* Copyright (C) 2015-2025 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #ifndef TEST_SSP_RAN_UNIFORM_DISK_H 17 #define TEST_SSP_RAN_UNIFORM_DISK_H 18 19 #include "ssp.h" 20 #include "test_ssp_utils.h" 21 #include <rsys/math.h> 22 23 #include <string.h> 24 25 #define NBS 1000000 26 27 #endif /* TEST_SSP_RAN_UNIFORM_DISK_H */ 28 29 #if TYPE_FLOAT==0 30 #define REAL double 31 #define TEST test_double 32 #define RNG_UNIFORM_R ssp_rng_uniform_double 33 #define RNG_UNIFORM_DISK_LOCAL ssp_ran_uniform_disk_local 34 #define RNG_UNIFORM_DISK ssp_ran_uniform_disk 35 #define R33_BASIS d33_basis 36 #define R33_MULR3 d33_muld3 37 #define R3_EQ_EPS d3_eq_eps 38 #define R3_NORMALIZE d3_normalize 39 #define SQRT sqrt 40 41 #elif TYPE_FLOAT==1 42 #define REAL float 43 #define TEST test_float 44 #define RNG_UNIFORM_R ssp_rng_uniform_float 45 #define RNG_UNIFORM_DISK_LOCAL ssp_ran_uniform_disk_float_local 46 #define RNG_UNIFORM_DISK ssp_ran_uniform_disk_float 47 #define R33_BASIS f33_basis 48 #define R33_MULR3 f33_mulf3 49 #define R3_EQ_EPS f3_eq_eps 50 #define R3_NORMALIZE f3_normalize 51 #define SQRT (float)sqrt 52 #else 53 #error "TYPE_FLOAT must be defined either 0 or 1" 54 #endif 55 56 static void 57 TEST() 58 { 59 struct ssp_rng* rng0, *rng1; 60 struct mem_allocator allocator; 61 int i; 62 int r, c, nb = 0; 63 REAL pt[3], pt2[3], up[3]; 64 REAL frame[9]; 65 REAL x_sum = 0, x2_sum = 0; 66 REAL mean, std; 67 REAL exp_mean = NBS * (20 * 20 / ((REAL)PI * 100 * 100)), exp_std = 0; 68 unsigned counts[10][10]; 69 uint64_t seed = 1234; 70 71 mem_init_proxy_allocator(&allocator, &mem_default_allocator); 72 73 CHK(ssp_rng_create(&allocator, SSP_RNG_THREEFRY, &rng0) == RES_OK); 74 CHK(ssp_rng_create(&allocator, SSP_RNG_THREEFRY, &rng1) == RES_OK); 75 76 up[0] = RNG_UNIFORM_R(rng1, -1, 1); 77 up[1] = RNG_UNIFORM_R(rng1, -1, 1); 78 up[2] = RNG_UNIFORM_R(rng1, -1, 1); 79 R3_NORMALIZE(up, up); 80 R33_BASIS(frame, up); 81 82 ssp_rng_set(rng0, seed); 83 ssp_rng_set(rng1, seed); 84 85 memset(counts, 0, sizeof(counts)); 86 FOR_EACH(i, 0, NBS) { 87 REAL tmp[3]; 88 RNG_UNIFORM_DISK_LOCAL(rng0, 100, pt, NULL); 89 RNG_UNIFORM_DISK(rng1, 100, up, pt2, NULL); 90 R33_MULR3(tmp, frame, pt); 91 CHK(R3_EQ_EPS(tmp, pt2, (REAL)1.e-6) == 1); 92 ASSERT(pt[2] == 0); 93 /* locate pt in a 10x10 grid */ 94 r = (int)((100 + pt[0]) / 20); 95 c = (int)((100 + pt[1]) / 20); 96 ++counts[r][c]; 97 } 98 99 FOR_EACH(r, 0, 10) { 100 int x = (r >= 5 ? r - 4 : r - 5) * 20; 101 FOR_EACH(c, 0, 10) { 102 int y = (c >= 5 ? c - 4 : c - 5) * 20; 103 int r2 = x * x + y * y; 104 if(r2 > 100 * 100) 105 /* this square is not (fully) in the disk */ 106 continue; 107 ++nb; 108 x_sum += (REAL)counts[r][c]; 109 x2_sum += (REAL)(counts[r][c] * counts[r][c]); 110 } 111 } 112 mean = x_sum / (REAL)nb; 113 std = x2_sum / (REAL)nb - mean*mean; 114 std = std > 0 ? SQRT(std) : 0; 115 printf("%g %g\n", mean, std); 116 CHK(fabs(mean - exp_mean) < 10); 117 CHK(fabs(std - exp_std) < 200); 118 119 CHK(ssp_rng_ref_put(rng0) == RES_OK); 120 CHK(ssp_rng_ref_put(rng1) == RES_OK); 121 122 check_memory_allocator(&allocator); 123 mem_shutdown_proxy_allocator(&allocator); 124 CHK(mem_allocated_size() == 0); 125 } 126 127 #undef REAL 128 #undef TEST 129 #undef RNG_UNIFORM_R 130 #undef RNG_UNIFORM_DISK_LOCAL 131 #undef RNG_UNIFORM_DISK 132 #undef R33_BASIS 133 #undef R33_MULR3 134 #undef R3_EQ_EPS 135 #undef R3_NORMALIZE 136 #undef SQRT 137 #undef TYPE_FLOAT