swf_H.c (8959B)
1 /* Copyright (C) 2024 |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 #define _XOPEN_SOURCE /* j1 */ 17 18 #include "swf.h" 19 #include "swf_tabulation.h" 20 21 #include <rsys/math.h> /* PI */ 22 23 #include <limits.h> 24 #include <math.h> 25 26 /* Location of the first positive zeros of the Bessel function J0. 27 * 28 * They are precalculated using the gsl_sf_bessel_zero_J0 function and defined 29 * in their hexadecimal floating-point representation, to ensure that they match 30 * the value calculated by the GSL to the nearest bit. The C program used to 31 * calculate these values can be summarized as follows: 32 * 33 * #include <stdio.h> 34 * #include <gsl/gsl_sf_bessel.h> 35 * int main(void) 36 * { 37 * for(unsigned i=1; i<=100; printf("%a\n", gsl_sf_bessel_zero_J0(i++))); 38 * return 0; 39 * } 40 */ 41 static const double j0_roots[] = { 42 NaN, 43 0x1.33d152e971b3bp+1, 0x1.6148f5b2c2e3ep+2, 0x1.14eb56cccded1p+3, 44 0x1.79544008272abp+3, 0x1.ddca13ef271e1p+3, 0x1.212313f8a19edp+4, 45 0x1.5362dd173f795p+4, 0x1.85a3b930156e8p+4, 0x1.b7e54a5fd5f1cp+4, 46 0x1.ea27591cbbed8p+4, 0x1.0e34e13a66fe6p+5, 0x1.275637a9619eap+5, 47 0x1.4077a7ed62935p+5, 0x1.59992c65d0d86p+5, 0x1.72bac0f810807p+5, 48 0x1.8bdc6293f064ep+5, 0x1.a4fe0ee444c7p+5, 0x1.be1fc41a4c5fcp+5, 49 0x1.d74180c9e41eap+5, 0x1.f06343d0971c9p+5, 0x1.04c28621f11ep+6, 50 0x1.11536cb22d724p+6, 0x1.1de4554a1c2d6p+6, 0x1.2a753fa82047ap+6, 51 0x1.37062b9535d1p+6, 0x1.439718e2e3795p+6, 0x1.50280769a218fp+6, 52 0x1.5cb8f7079c7aep+6, 0x1.6949e79fb1f06p+6, 0x1.75dad918abf93p+6, 53 0x1.826bcb5c9b61dp+6, 0x1.8efcbe585425p+6, 0x1.9b8db1fb017fbp+6, 54 0x1.a81ea635cd31dp+6, 0x1.b4af9afb9610bp+6, 0x1.c1409040b2ea7p+6, 55 0x1.cdd185fabf635p+6, 0x1.da627c2070f25p+6, 0x1.e6f372a97287p+6, 56 0x1.f384698e45aa8p+6, 0x1.000ab06414167p+7, 0x1.06532c287ecd1p+7, 57 0x1.0c9ba8119dec5p+7, 0x1.12e4241ced385p+7, 0x1.192ca04822089p+7, 58 0x1.1f751c9124fd1p+7, 0x1.25bd98f60c82fp+7, 0x1.2c06157518087p+7, 59 0x1.324e920cabc8fp+7, 0x1.38970ebb4d19ep+7, 0x1.3edf8b7f9f282p+7, 60 0x1.4528085860164p+7, 0x1.4b708544666ep+7, 0x1.51b902429edb7p+7, 61 0x1.58017f520a27dp+7, 0x1.5e49fc71bb6b2p+7, 0x1.649279a0d6701p+7, 62 0x1.6adaf6de8e417p+7, 0x1.7123742a23dep+7, 0x1.776bf182e50dcp+7, 63 0x1.7db46ee82b546p+7, 0x1.83fcec595afe4p+7, 0x1.8a4569d5e2445p+7, 64 0x1.908de75d3884dp+7, 0x1.96d664eedd8edp+7, 0x1.9d1ee28a58fdap+7, 65 0x1.a367602f39a33p+7, 0x1.a9afdddd15001p+7, 0x1.aff85b9386c73p+7, 66 0x1.b640d952306b7p+7, 0x1.bc895718b8b88p+7, 0x1.c2d1d4e6cb72dp+7, 67 0x1.c91a52bc19007p+7, 0x1.cf62d0985618dp+7, 0x1.d5ab4e7b3b7a4p+7, 68 0x1.dbf3cc6485a69p+7, 0x1.e23c4a53f4a4p+7, 0x1.e884c8494bc31p+7, 69 0x1.eecd46445169ap+7, 0x1.f515c444cee1p+7, 0x1.fb5e424a90284p+7, 70 0x1.00d3602ab1e5p+8, 0x1.03f79f328d5a5p+8, 0x1.071bde3cc40b1p+8, 71 0x1.0a401d49409dp+8, 0x1.0d645c57eeb4ep+8, 0x1.10889b68bae7ap+8, 72 0x1.13acda7b92acep+8, 0x1.16d119906452p+8, 0x1.19f558a71eee5p+8, 73 0x1.1d1997bfb2581p+8, 0x1.203dd6da0f199p+8, 0x1.236215f626682p+8, 74 0x1.26865513ea1a6p+8, 0x1.29aa94334ca02p+8, 0x1.2cced35440fa3p+8, 75 0x1.2ff31276bab31p+8, 0x1.3317519aadd77p+8, 0x1.363b90c00ef04p+8, 76 0x1.395fcfe6d2fbfp+8 77 }; 78 static const size_t j0_nroots = sizeof(j0_roots)/sizeof(*j0_roots); 79 80 /******************************************************************************* 81 * Helper functions 82 ******************************************************************************/ 83 static INLINE res_T 84 check_swf_H_tabulate_args(const struct swf_H_tabulate_args* args) 85 { 86 if(!args) return RES_BAD_ARG; 87 88 /* X cannot be negative */ 89 if(args->x_min < 0) 90 return RES_BAD_ARG; 91 92 /* X range cannot be degenerated */ 93 if(args->x_min >= args->x_max) 94 return RES_BAD_ARG; 95 96 /* Delta X cannot be null */ 97 if(args->step <= 0) 98 return RES_BAD_ARG; 99 100 return RES_OK; 101 } 102 103 /* H(x) = 1 - 2 Sum(k=1..INF)[exp(-x*ak^2) / (ak*J1(ak))] 104 * H'(x) = 2 Sum(k=1..INF)[exp(-x*ak^2) * ak/J1(ak)] 105 * H"(x) = - 2 Sum(k=1..INF)[exp(-x*ak^2) * ak^3/J1(ak)] */ 106 static INLINE double 107 H2d 108 (const double x, 109 double* dHx, /* First derivative. NULL <=> do not compute it */ 110 double* d2Hx) /* Second derivative. NULL <=> do not compute it */ 111 { 112 double Hx = 0; /* Value of the function */ 113 114 /* Sum of the terms in the series */ 115 double sum = 0; /* Sum of the terms */ 116 double sumd = 0; /* Sum of the derivative terms */ 117 double sumd2 = 0; /* Sum of the second derivative terms */ 118 119 unsigned k = 1; /* Index of the serie term */ 120 121 /* Have the calculations converged or are there errors? */ 122 int error = 1; 123 int errord = dHx != NULL; 124 int errord2 = d2Hx != NULL; 125 126 ASSERT(x >= 0); 127 128 if(x == 0) return 0; 129 130 do { 131 double ak, j, t, term, sum_next; 132 CHK(k < j0_nroots); 133 134 ak = j0_roots[k]; 135 j = j1(ak); 136 t = exp(-x*ak*ak); 137 term = t / (ak*j); 138 sum_next = sum + term; 139 140 error = sum_next != sum; 141 sum = sum_next; 142 143 if(dHx) { /* Derivative */ 144 const double termd = t/j * ak; 145 const double sumd_next = sumd + termd; 146 errord = sumd_next != sumd; 147 sumd = sumd_next; 148 } 149 150 if(d2Hx) { /* Second derivative */ 151 const double termd2 = t/j * ak*ak*ak; 152 const double sumd2_next = sumd2 + termd2; 153 errord2 = sumd2_next != sumd2; 154 sumd2 = sumd2_next; 155 } 156 157 } while((error || errord || errord2) && ++k < UINT_MAX); 158 159 Hx = 1.0 - 2.0 * sum; 160 161 if(dHx) *dHx = 2.0 * sumd; /* Derivative */ 162 if(d2Hx) *d2Hx = -2.0 * sumd2; /* Second Derivative */ 163 164 return Hx; 165 } 166 167 /* H(x) = 1 + 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x)] 168 * H'(x) = - 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x) * (PI*k)^2] 169 * H"(x) = + 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x) * (PI*k)^4] */ 170 static INLINE double 171 H3d 172 (const double x, 173 double* dHx, /* Derivative. NULL <=> do not compute it */ 174 double* d2Hx) /* Second derivative. NULL <=> do not compute it */ 175 { 176 double Hx = 0; /* Value of the function */ 177 178 /* Sum of the terms in the series */ 179 double sum = 0; /* Sum of the terms */ 180 double sumd = 0; /* Sum of the derivative terms */ 181 double sumd2 = 0; /* Sum of the second derivative terms */ 182 183 double sign = -1; /* Sign of the serie term */ 184 unsigned k = 1; /* Index of the serie term */ 185 186 /* Have the calculations converged or are there errors? */ 187 int error = 1; 188 int errord = dHx != NULL; 189 int errord2 = d2Hx != NULL; 190 191 ASSERT(x >= 0); /* Check pre-condition */ 192 193 /* Do not attempt to calculate 0; it would be numerically catastrophic. 194 * Simply return its value */ 195 if(x == 0) return 0; 196 197 do { 198 const double t = PI*(double)k; 199 const double t2 = t*t; 200 const double term = sign * exp(-t2*x); 201 const double sum_next = sum + term; 202 203 error = sum_next != sum; 204 sum = sum_next; 205 sign = -sign; 206 207 if(dHx) { /* Derivative */ 208 const double termd = term * t2; 209 const double sumd_next = sumd + termd; 210 errord = sumd_next != sumd; 211 sumd = sumd_next; 212 } 213 214 if(d2Hx) { /* Second derivative */ 215 const double termd2 = term * t2*t2; 216 const double sumd2_next = sumd2 + termd2; 217 errord2 = sumd2_next != sumd2; 218 sumd2 = sumd2_next; 219 } 220 221 } while((error || errord || errord2) && ++k < UINT_MAX); 222 223 Hx = 1.0 + 2.0 * sum; 224 225 if(dHx) *dHx = -2.0 * sumd; /* Derivative */ 226 if(d2Hx) *d2Hx = 2.0 * sumd2; /* Second derivative */ 227 return Hx; 228 } 229 230 /* Define generic dimension functions */ 231 #define SWF_DIMENSION 2 232 #include "swf_HXd.h" 233 #define SWF_DIMENSION 3 234 #include "swf_HXd.h" 235 236 /******************************************************************************* 237 * Exported symbols 238 ******************************************************************************/ 239 double 240 swf_H2d_eval(const double x) 241 { 242 return H2d(x, NULL, NULL); 243 } 244 245 double 246 swf_H3d_eval(const double x) 247 { 248 return H3d(x, NULL, NULL); 249 } 250 251 double 252 swf_H2d_inverse(const double y) 253 { 254 return H_inverse2d(y, 255 SWF_H2D_TABULATE_ARGS_DEFAULT.x_min, 256 SWF_H2D_TABULATE_ARGS_DEFAULT.x_max, 257 1.0e-12); /* Epsilon */ 258 } 259 260 double 261 swf_H3d_inverse(const double y) 262 { 263 return H_inverse3d(y, 264 SWF_H3D_TABULATE_ARGS_DEFAULT.x_min, 265 SWF_H3D_TABULATE_ARGS_DEFAULT.x_max, 266 1.0e-12); /* Epsilon */ 267 } 268 269 res_T 270 swf_H2d_tabulate 271 (const struct swf_H_tabulate_args* args, 272 struct swf_tabulation** out_tab) 273 { 274 return H_tabulate2d(args, out_tab); 275 } 276 277 res_T 278 swf_H3d_tabulate 279 (const struct swf_H_tabulate_args* args, 280 struct swf_tabulation** out_tab) 281 { 282 return H_tabulate3d(args, out_tab); 283 }