commit 596ae52a870c494f3a2d9971abea4a48c0502564
parent b20d98f6b06ea0e385280a38cbe515f7ce7461fd
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 13 Mar 2024 15:43:17 +0100
Add and test the 2D H function
This is the 2D version of the H-function, i.e. the cumulative
probability density function used to sample a passage time for a circle.
Its API is the same as that of its 3D version: the function can be
evaluated inverted and tabulated to improve inversion efficiency.
The tests are the same as for the 3D H function. We check its API and
ensure its numerical accuracy.
Diffstat:
| M | Makefile | | | 4 | ++-- |
| M | config.mk | | | 2 | +- |
| M | src/swf.h | | | 36 | +++++++++++++++++++++++++++++------- |
| M | src/swf_H.c | | | 232 | ++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------- |
| A | src/swf_HXd.h | | | 127 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | src/test_swf_H2d.c | | | 35 | +++++++++++++++++++++++++++++++++++ |
| M | src/test_swf_H3d.c | | | 97 | +++++-------------------------------------------------------------------------- |
| A | src/test_swf_HXd.h | | | 136 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
8 files changed, 489 insertions(+), 180 deletions(-)
diff --git a/Makefile b/Makefile
@@ -110,7 +110,7 @@ lint:
################################################################################
# Tests
################################################################################
-TEST_SRC = src/test_swf_H3d.c
+TEST_SRC = src/test_swf_H2d.c src/test_swf_H3d.c
TEST_OBJ = $(TEST_SRC:.c=.o)
TEST_DEP = $(TEST_SRC:.c=.d)
@@ -137,5 +137,5 @@ $(TEST_DEP): config.mk swf-local.pc
$(TEST_OBJ): config.mk swf-local.pc
$(CC) $(CFLAGS_EXE) $(SWF_CFLAGS) $(RSYS_CFLAGS) -c $(@:.o=.c) -o $@
-test_swf_H3d: config.mk swf-local.pc $(LIBNAME)
+test_swf_H2d test_swf_H3d: config.mk swf-local.pc $(LIBNAME)
$(CC) $(CFLAGS_EXE) -o $@ src/$@.o $(LDFLAGS_EXE) $(SWF_LIBS) $(RSYS_LIBS) -lm
diff --git a/config.mk b/config.mk
@@ -47,7 +47,7 @@ CFLAGS_HARDENED =\
-fstack-protector-strong
CFLAGS_COMMON =\
- -std=c89\
+ -std=c99\
-pedantic\
-fvisibility=hidden\
-fstrict-aliasing\
diff --git a/src/swf.h b/src/swf.h
@@ -62,6 +62,9 @@ struct swf_H_tabulate_args {
/* The H function's default tab values guarantee both numerical accuracy and a
* small memory footprint */
+#define SWF_H2D_TABULATE_ARGS_DEFAULT__ {8.0e-3, 6.05, 1e-3, 0, NULL}
+static const struct swf_H_tabulate_args SWF_H2D_TABULATE_ARGS_DEFAULT =
+ SWF_H2D_TABULATE_ARGS_DEFAULT__;
#define SWF_H3D_TABULATE_ARGS_DEFAULT__ {7.5e-3, 3.337, 1e-3, 0, NULL}
static const struct swf_H_tabulate_args SWF_H3D_TABULATE_ARGS_DEFAULT =
SWF_H3D_TABULATE_ARGS_DEFAULT__;
@@ -71,19 +74,25 @@ struct swf_tabulation;
BEGIN_DECLS
-/* This function evaluates the 3D H function, as specified in:
+/* These functions evaluate the H function, as specified in:
*
* "The floating random walk and its application to Monte Carlo solutions of
* heat equations" - Haji-Sheikh A. and Sparrow E.M., SIAM Journal on Applied
* Mathematics, 14(2): 370-389, 1966
*
- * This function is expressed as an infinite series, which is numerically
- * evaluated by summing all terms, until convergence is reached :
+ * These function are expressed as an infinite series, which is numerically
+ * evaluated by summing all terms, until convergence is reached.
+ * Two functions are defined here, H2d and H3d, which are used in 2D and 3D
+ * respectively:
+ * +INF
+ * __
+ * H2d(x) = 1 - 2 >_ exp(-x*ak^2) / (ak*J1(ak)); ak is the k^th root of J0
+ * k=1
*
- * +INF
- * __
- * H(x) = 1 + 2 >_ (-1)^k * exp(-(PI*k)^2 * x)
- * k=1
+ * +INF
+ * __
+ * H3d(x) = 1 + 2 >_ (-1)^k * exp(-(PI*k)^2*x)
+ * k=1
*
* The H function is the cumulated function of the probability density used to
* sample passage times for a given sphere radius. Its parameter x is positive
@@ -92,6 +101,10 @@ BEGIN_DECLS
* - tau > 0 the time interval [s]
* - r > 0 is the radius of the sphere [m] */
SWF_API double
+swf_H2d_eval
+ (const double x);
+
+SWF_API double
swf_H3d_eval
(const double x);
@@ -99,6 +112,15 @@ SWF_API double
swf_H3d_inverse
(const double y); /* in [0, 1[ */
+SWF_API double
+swf_H2d_inverse
+ (const double y); /* in [0, 1[ */
+
+SWF_API res_T
+swf_H2d_tabulate
+ (const struct swf_H_tabulate_args* args,
+ struct swf_tabulation** tab);
+
SWF_API res_T
swf_H3d_tabulate
(const struct swf_H_tabulate_args* args,
diff --git a/src/swf_H.c b/src/swf_H.c
@@ -13,6 +13,8 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
+#define _XOPEN_SOURCE /* j1 */
+
#include "swf.h"
#include "swf_tabulation.h"
@@ -21,6 +23,60 @@
#include <limits.h>
#include <math.h>
+/* Location of the first positive zero of the Bessel function J0.
+ *
+ * They are precalculated using the gsl_sf_bessel_zero_J0 function and defined
+ * in their hexadecimal floating-point representation, to ensure that they match
+ * the value calculated by the GSL to the nearest bit. The C program used to
+ * calculate these values can be summarized as follows:
+ *
+ * #include <stdio.h>
+ * #include <gsl/gsl_sf_bessel.h>
+ * int main(void)
+ * {
+ * for(unsigned i=1; i<=100; printf("%a\n", gsl_sf_bessel_zero_J0(i++)));
+ * return 0;
+ * }
+ */
+static const double j0_roots[] = {
+ NaN,
+ 0x1.33d152e971b3bp+1, 0x1.6148f5b2c2e3ep+2, 0x1.14eb56cccded1p+3,
+ 0x1.79544008272abp+3, 0x1.ddca13ef271e1p+3, 0x1.212313f8a19edp+4,
+ 0x1.5362dd173f795p+4, 0x1.85a3b930156e8p+4, 0x1.b7e54a5fd5f1cp+4,
+ 0x1.ea27591cbbed8p+4, 0x1.0e34e13a66fe6p+5, 0x1.275637a9619eap+5,
+ 0x1.4077a7ed62935p+5, 0x1.59992c65d0d86p+5, 0x1.72bac0f810807p+5,
+ 0x1.8bdc6293f064ep+5, 0x1.a4fe0ee444c7p+5, 0x1.be1fc41a4c5fcp+5,
+ 0x1.d74180c9e41eap+5, 0x1.f06343d0971c9p+5, 0x1.04c28621f11ep+6,
+ 0x1.11536cb22d724p+6, 0x1.1de4554a1c2d6p+6, 0x1.2a753fa82047ap+6,
+ 0x1.37062b9535d1p+6, 0x1.439718e2e3795p+6, 0x1.50280769a218fp+6,
+ 0x1.5cb8f7079c7aep+6, 0x1.6949e79fb1f06p+6, 0x1.75dad918abf93p+6,
+ 0x1.826bcb5c9b61dp+6, 0x1.8efcbe585425p+6, 0x1.9b8db1fb017fbp+6,
+ 0x1.a81ea635cd31dp+6, 0x1.b4af9afb9610bp+6, 0x1.c1409040b2ea7p+6,
+ 0x1.cdd185fabf635p+6, 0x1.da627c2070f25p+6, 0x1.e6f372a97287p+6,
+ 0x1.f384698e45aa8p+6, 0x1.000ab06414167p+7, 0x1.06532c287ecd1p+7,
+ 0x1.0c9ba8119dec5p+7, 0x1.12e4241ced385p+7, 0x1.192ca04822089p+7,
+ 0x1.1f751c9124fd1p+7, 0x1.25bd98f60c82fp+7, 0x1.2c06157518087p+7,
+ 0x1.324e920cabc8fp+7, 0x1.38970ebb4d19ep+7, 0x1.3edf8b7f9f282p+7,
+ 0x1.4528085860164p+7, 0x1.4b708544666ep+7, 0x1.51b902429edb7p+7,
+ 0x1.58017f520a27dp+7, 0x1.5e49fc71bb6b2p+7, 0x1.649279a0d6701p+7,
+ 0x1.6adaf6de8e417p+7, 0x1.7123742a23dep+7, 0x1.776bf182e50dcp+7,
+ 0x1.7db46ee82b546p+7, 0x1.83fcec595afe4p+7, 0x1.8a4569d5e2445p+7,
+ 0x1.908de75d3884dp+7, 0x1.96d664eedd8edp+7, 0x1.9d1ee28a58fdap+7,
+ 0x1.a367602f39a33p+7, 0x1.a9afdddd15001p+7, 0x1.aff85b9386c73p+7,
+ 0x1.b640d952306b7p+7, 0x1.bc895718b8b88p+7, 0x1.c2d1d4e6cb72dp+7,
+ 0x1.c91a52bc19007p+7, 0x1.cf62d0985618dp+7, 0x1.d5ab4e7b3b7a4p+7,
+ 0x1.dbf3cc6485a69p+7, 0x1.e23c4a53f4a4p+7, 0x1.e884c8494bc31p+7,
+ 0x1.eecd46445169ap+7, 0x1.f515c444cee1p+7, 0x1.fb5e424a90284p+7,
+ 0x1.00d3602ab1e5p+8, 0x1.03f79f328d5a5p+8, 0x1.071bde3cc40b1p+8,
+ 0x1.0a401d49409dp+8, 0x1.0d645c57eeb4ep+8, 0x1.10889b68bae7ap+8,
+ 0x1.13acda7b92acep+8, 0x1.16d119906452p+8, 0x1.19f558a71eee5p+8,
+ 0x1.1d1997bfb2581p+8, 0x1.203dd6da0f199p+8, 0x1.236215f626682p+8,
+ 0x1.26865513ea1a6p+8, 0x1.29aa94334ca02p+8, 0x1.2cced35440fa3p+8,
+ 0x1.2ff31276bab31p+8, 0x1.3317519aadd77p+8, 0x1.363b90c00ef04p+8,
+ 0x1.395fcfe6d2fbfp+8
+};
+static const size_t j0_nroots = sizeof(j0_roots)/sizeof(*j0_roots);
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -44,6 +100,70 @@ check_swf_H_tabulate_args(const struct swf_H_tabulate_args* args)
return RES_OK;
}
+/* H(x) = 1 - 2 Sum(k=1..INF)[exp(-x*ak^2) / (ak*J1(ak))]
+ * H'(x) = 2 Sum(k=1..INF)[exp(-x*ak^2) * ak/J1(ak)]
+ * H"(x) = - 2 Sum(k=1..INF)[exp(-x*ak^2) * ak^3/J1(ak)] */
+static INLINE double
+H2d
+ (const double x,
+ double* dHx, /* First derivative. NULL <=> do not compute it */
+ double* d2Hx) /* Second derivative. NULL <=> do not compute it */
+{
+ double Hx = 0; /* Value of the function */
+
+ /* Sum of the terms in the series */
+ double sum = 0; /* Sum of the terms */
+ double sumd = 0; /* Sum of the derivative terms */
+ double sumd2 = 0; /* Sum of the second derivative terms */
+
+ unsigned k = 1; /* Index of the serie term */
+
+ /* Have the calculations converged or are there errors? */
+ int error = 1;
+ int errord = dHx != NULL;
+ int errord2 = d2Hx != NULL;
+
+ ASSERT(x >= 0);
+
+ if(x == 0) return 0;
+
+ do {
+ double ak, j, t, term, sum_next;
+ CHK(k < j0_nroots);
+
+ ak = j0_roots[k];
+ j = j1(ak);
+ t = exp(-x*ak*ak);
+ term = t / (ak*j);
+ sum_next = sum + term;
+
+ error = sum_next != sum;
+ sum = sum_next;
+
+ if(dHx) { /* Derivative */
+ const double termd = t/j * ak;
+ const double sumd_next = sumd + termd;
+ errord = sumd_next != sumd;
+ sumd = sumd_next;
+ }
+
+ if(d2Hx) { /* Second derivative */
+ const double termd2 = t/j * ak*ak*ak;
+ const double sumd2_next = sumd2 + termd2;
+ errord2 = sumd2_next != sumd2;
+ sumd2 = sumd2_next;
+ }
+
+ } while((error || errord || errord2) && ++k < UINT_MAX);
+
+ Hx = 1.0 - 2.0 * sum;
+
+ if(dHx) *dHx = 2.0 * sumd; /* Derivative */
+ if(d2Hx) *d2Hx = -2.0 * sumd2; /* Second Derivative */
+
+ return Hx;
+}
+
/* H(x) = 1 + 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x)]
* H'(x) = - 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x) * (PI*k)^2]
* H"(x) = + 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x) * (PI*k)^4] */
@@ -107,36 +227,51 @@ H3d
return Hx;
}
+/* Define generic dimension functions */
+#define SWF_DIMENSION 2
+#include "swf_HXd.h"
+#define SWF_DIMENSION 3
+#include "swf_HXd.h"
+
/*******************************************************************************
* Exported symbols
******************************************************************************/
double
+swf_H2d_eval(const double x)
+{
+ return H2d(x, NULL, NULL);
+}
+
+res_T
+swf_H2d_tabulate
+ (const struct swf_H_tabulate_args* args,
+ struct swf_tabulation** out_tab)
+{
+ return H_tabulate2d(args, out_tab);
+}
+
+double
swf_H3d_eval(const double x)
{
return H3d(x, NULL, NULL);
}
double
-swf_H3d_inverse(const double y)
+swf_H2d_inverse(const double y)
{
- const double epsilon = 1.0e-12;
- double x_min = SWF_H_TABULATE_ARGS_DEFAULT.x_min;
- double x_max = SWF_H_TABULATE_ARGS_DEFAULT.x_max;
- ASSERT(y >= 0 && y < 1);
-
- if(y == 0) return 0;
-
- while(x_max - x_min > epsilon) {
- const double x = (x_min + x_max) * 0.5;
- const double f_x = swf_H3d_eval(x);
+ return H_inverse2d(y,
+ SWF_H2D_TABULATE_ARGS_DEFAULT.x_min,
+ SWF_H2D_TABULATE_ARGS_DEFAULT.x_max,
+ 1.0e-12); /* Epsilon */
+}
- if(f_x < y) {
- x_min = x;
- } else {
- x_max = x;
- }
- }
- return (x_min + x_max) * 0.5;
+double
+swf_H3d_inverse(const double y)
+{
+ return H_inverse3d(y,
+ SWF_H3D_TABULATE_ARGS_DEFAULT.x_min,
+ SWF_H3D_TABULATE_ARGS_DEFAULT.x_max,
+ 1.0e-12); /* Epsilon */
}
res_T
@@ -144,64 +279,5 @@ swf_H3d_tabulate
(const struct swf_H_tabulate_args* args,
struct swf_tabulation** out_tab)
{
- struct swf_tabulation* tab = NULL;
- size_t i = 0;
- double x = 0;
- res_T res = RES_OK;
-
- if(!out_tab) { res = RES_BAD_ARG; goto error; }
-
- if((res = check_swf_H_tabulate_args(args)) != RES_OK) goto error;
- if((res = tabulation_create(args->allocator, &tab)) != RES_OK) goto error;
-
- /* Setup the tabulation */
- i = 0;
- for(x = args->x_min; x <= args->x_max; x += x*args->step) {
- const struct item* pitem = darray_item_cdata_get(&tab->items) + i - 1;
- struct item item = ITEM_NULL;
-
- /* Evaluate H(x) and its derivatives up to order 2 */
- item.x = MMIN(x, args->x_max);
- item.fx = H3d(item.x, &item.dfx, &item.d2fx);
-
- if(i > 0) {
- /* Detect oscillations due to numerical problems.
- * These oscillations occur close to 0 */
- if(item.fx < pitem->fx) {
- res = RES_BAD_OP;
- goto error;
- }
- /* Do not tabulate argument whose value is (numerically) equal to the
- * previous one: it would artificially add staircase steps */
- if(item.fx == pitem->fx) {
- continue;
- }
- }
-
- ++i;
-
- res = darray_item_push_back(&tab->items, &item);
- if(res != RES_OK) goto error;
- }
- ASSERT(darray_item_size_get(&tab->items) == i);
-
-
- if(args->normalize) {
- struct item* items = darray_item_data_get(&tab->items);
- const size_t nitems = darray_item_size_get(&tab->items);
- const double norm = items[nitems - 1].fx;
-
- /* Normalize the tabulation */
- FOR_EACH(i, 0, nitems) items[i].fx /= norm;
- }
-
-exit:
- if(out_tab) *out_tab = tab;
- return res;
-error:
- if(tab) {
- SWF(tabulation_ref_put(tab));
- tab = NULL;
- }
- goto exit;
+ return H_tabulate3d(args, out_tab);
}
diff --git a/src/swf_HXd.h b/src/swf_HXd.h
@@ -0,0 +1,127 @@
+/* Copyright (C) 2024 |Méso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "swf.h"
+#include "swf_tabulation.h"
+
+#include <rsys/rsys.h>
+
+#ifndef SWF_DIMENSION
+ #error "The SWF_DIMENSION macro must be defined"
+#endif
+
+#if SWF_DIMENSION != 2 && SWF_DIMENSION != 3
+ #error "Invalid SWF_DIMEMSION"
+#endif
+
+#define Xd(Name) CONCAT(CONCAT(Name, SWF_DIMENSION), d)
+
+/*******************************************************************************
+ * Helper generic dimension functions
+ ******************************************************************************/
+static INLINE double
+Xd(H_inverse)
+ (const double y,
+ const double x_min,
+ const double x_max,
+ const double epsilon)
+{
+ double x0 = x_min;
+ double x1 = x_max;
+ ASSERT(y >= 0 && y < 1 && x_min < x_max && epsilon > 0);
+
+ if(y == 0) return 0;
+
+ while(x1 - x0 > epsilon) {
+ const double x = (x0 + x1) * 0.5;
+ const double f_x = Xd(H)(x, NULL, NULL);
+
+ if(f_x < y) {
+ x0 = x;
+ } else {
+ x1 = x;
+ }
+ }
+ return (x0 + x1) * 0.5;
+}
+
+static INLINE res_T
+Xd(H_tabulate)
+ (const struct swf_H_tabulate_args* args,
+ struct swf_tabulation** out_tab)
+{
+ struct swf_tabulation* tab = NULL;
+ size_t i = 0;
+ double x = 0;
+ res_T res = RES_OK;
+
+ if(!out_tab) { res = RES_BAD_ARG; goto error; }
+
+ if((res = check_swf_H_tabulate_args(args)) != RES_OK) goto error;
+ if((res = tabulation_create(args->allocator, &tab)) != RES_OK) goto error;
+
+ /* Setup the tabulation */
+ i = 0;
+ for(x = args->x_min; x <= args->x_max; x += x*args->step) {
+ const struct item* pitem = darray_item_cdata_get(&tab->items) + i - 1;
+ struct item item = ITEM_NULL;
+
+ /* Evaluate H(x) and its derivatives up to order 2 */
+ item.x = x;
+ item.fx = Xd(H)(item.x, &item.dfx, &item.d2fx);
+
+ if(i > 0) {
+ /* Detect oscillations due to numerical problems.
+ * These oscillations occur close to 0 */
+ if(item.fx < pitem->fx) {
+ res = RES_BAD_OP;
+ goto error;
+ }
+ /* Do not tabulate argument whose value is (numerically) equal to the
+ * previous one: it would artificially add staircase steps */
+ if(item.fx == pitem->fx) {
+ continue;
+ }
+ }
+
+ ++i;
+
+ res = darray_item_push_back(&tab->items, &item);
+ if(res != RES_OK) goto error;
+ }
+ ASSERT(darray_item_size_get(&tab->items) == i);
+
+ if(args->normalize) {
+ struct item* items = darray_item_data_get(&tab->items);
+ const size_t nitems = darray_item_size_get(&tab->items);
+ const double norm = items[nitems - 1].fx;
+
+ /* Normalize the tabulation */
+ FOR_EACH(i, 0, nitems) items[i].fx /= norm;
+ }
+
+exit:
+ if(out_tab) *out_tab = tab;
+ return res;
+error:
+ if(tab) {
+ SWF(tabulation_ref_put(tab));
+ tab = NULL;
+ }
+ goto exit;
+}
+
+#undef SWF_DIMENSION
+#undef Xd
diff --git a/src/test_swf_H2d.c b/src/test_swf_H2d.c
@@ -0,0 +1,35 @@
+/* Copyright (C) 2024 |Méso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "swf.h"
+#include <rsys/mem_allocator.h>
+
+/* Define generic dimension helper functions */
+#define TEST_SWF_DIMENSION 2
+#include "test_swf_HXd.h"
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ (void)argc, (void)argv;
+
+ check_tabulation_creation_2d();
+ check_inversion_2d();
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_swf_H3d.c b/src/test_swf_H3d.c
@@ -14,98 +14,11 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include "swf.h"
-
-#include <rsys/rsys.h>
-#include <rsys/math.h>
#include <rsys/mem_allocator.h>
-#include <stdlib.h>
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static INLINE void
-check_tabulation_creation(void)
-{
- struct swf_H_tabulate_args args = SWF_H3D_TABULATE_ARGS_DEFAULT;
- struct swf_tabulation* tab = NULL;
-
- CHK(swf_H3d_tabulate(NULL, &tab) == RES_BAD_ARG);
- CHK(swf_H3d_tabulate(&args, NULL) == RES_BAD_ARG);
- CHK(swf_H3d_tabulate(&args, &tab) == RES_OK);
-
- CHK(swf_tabulation_ref_get(NULL) == RES_BAD_ARG);
- CHK(swf_tabulation_ref_get(tab) == RES_OK);
- CHK(swf_tabulation_ref_put(NULL) == RES_BAD_ARG);
- CHK(swf_tabulation_ref_put(tab) == RES_OK);
- CHK(swf_tabulation_ref_put(tab) == RES_OK);
-
- args.x_min = SWF_H3D_TABULATE_ARGS_DEFAULT.x_max;
- args.x_max = SWF_H3D_TABULATE_ARGS_DEFAULT.x_min;
- CHK(swf_H3d_tabulate(&args, &tab) == RES_BAD_ARG);
-
- args.x_max = args.x_min;
- CHK(swf_H3d_tabulate(&args, &tab) == RES_BAD_ARG);
-
- args.x_min = SWF_H3D_TABULATE_ARGS_DEFAULT.x_min;
- args.x_max = SWF_H3D_TABULATE_ARGS_DEFAULT.x_max;
- args.step = 0;
- CHK(swf_H3d_tabulate(&args, &tab) == RES_BAD_ARG);
-
- args.step = 1.0e-2;
- CHK(swf_H3d_tabulate(&args, &tab) == RES_OK);
- CHK(swf_tabulation_ref_put(tab) == RES_OK);
-}
-
-static void
-check_inversion(void)
-{
- const double nsteps = 10;
-
- struct swf_H_tabulate_args args = SWF_H3D_TABULATE_ARGS_DEFAULT;
- struct swf_tabulation* tab = NULL;
- double pHx = 0;
- double x0 = 0;
- size_t i = 0;
-
- CHK(swf_H3d_tabulate(&args, &tab) == RES_OK);
-
- for(x0 = args.x_min; x0 <= args.x_max; x0 += x0*args.step) {
- const double x1 = MMIN(x0 + x0*args.step, args.x_max);
- const double delta_x = (x1 - x0) / (double)nsteps;
-
- for(i = 0; i < nsteps; ++i) {
- const double x = x0 + (double)i*delta_x;
- const double Hx = swf_H3d_eval(x);
- const double xl = swf_tabulation_inverse(tab, SWF_LINEAR, Hx);
- const double xq = swf_tabulation_inverse(tab, SWF_QUADRATIC, Hx);
- const double xi = swf_H3d_inverse(Hx);
- const double errl = fabs((x - xl) / x);
- const double errq = fabs((x - xq) / x);
- const double erri = fabs((x - xi) / x);
-
- /* Do not check the value of x if the corresponding H is indistinguishable
- * from an H calculated from a previous x. In other words, the H's are
- * numerically equal and any of the corresponding x-values is a valid
- * inversion result, numerically speaking. */
- if(Hx == pHx) continue;
-
- printf("%e %e %e %e %e\n", x, Hx, errq, errl, erri);
- if(1e-9 < Hx && Hx < (1.0 - 1e-9)) {
- CHK(errl < 1.e-5);
- CHK(errq < 1.e-7);
- CHK(erri < 1.e-7);
- } else {
- CHK(errl < 1.e-3);
- CHK(errq < 1.e-3);
- CHK(erri < 1.e-3);
- }
- pHx = Hx;
- }
- }
-
- CHK(swf_tabulation_ref_put(tab) == RES_OK);
-}
+/* Define generic dimension helper functions */
+#define TEST_SWF_DIMENSION 3
+#include "test_swf_HXd.h"
/*******************************************************************************
* The test
@@ -115,8 +28,8 @@ main(int argc, char** argv)
{
(void)argc, (void)argv;
- check_tabulation_creation();
- check_inversion();
+ check_tabulation_creation_3d();
+ check_inversion_3d();
CHK(mem_allocated_size() == 0);
return 0;
}
diff --git a/src/test_swf_HXd.h b/src/test_swf_HXd.h
@@ -0,0 +1,136 @@
+/* Copyright (C) 2024 |Méso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "swf.h"
+
+#include <rsys/rsys.h>
+#include <rsys/math.h>
+
+#include <stdlib.h>
+
+#ifndef TEST_SWF_DIMENSION
+ #error "The TEST_SWF_DIMENSION macro must be defined"
+#endif
+
+/* 2D H function */
+#if TEST_SWF_DIMENSION == 2
+ #define SWF_HXD_TABULATE_ARGS_DEFAULT SWF_H2D_TABULATE_ARGS_DEFAULT
+ #define swf_HXd_eval swf_H2d_eval
+ #define swf_HXd_inverse swf_H2d_inverse
+ #define swf_HXd_tabulate swf_H2d_tabulate
+
+/* 3D H function */
+#elif TEST_SWF_DIMENSION == 3
+ #define SWF_HXD_TABULATE_ARGS_DEFAULT SWF_H3D_TABULATE_ARGS_DEFAULT
+ #define swf_HXd_eval swf_H3d_eval
+ #define swf_HXd_inverse swf_H3d_inverse
+ #define swf_HXd_tabulate swf_H3d_tabulate
+#else
+ #error "Invalid TEST_SWF_DIMENSION"
+#endif
+
+/* Macro defining a generic function name for the dimension*/
+#define Xd(Name) CONCAT(CONCAT(CONCAT(Name, _), TEST_SWF_DIMENSION), d)
+
+static INLINE void
+Xd(check_tabulation_creation)(void)
+{
+ struct swf_H_tabulate_args args = SWF_HXD_TABULATE_ARGS_DEFAULT;
+ struct swf_tabulation* tab = NULL;
+
+ CHK(swf_HXd_tabulate(NULL, &tab) == RES_BAD_ARG);
+ CHK(swf_HXd_tabulate(&args, NULL) == RES_BAD_ARG);
+ CHK(swf_HXd_tabulate(&args, &tab) == RES_OK);
+
+ CHK(swf_tabulation_ref_get(NULL) == RES_BAD_ARG);
+ CHK(swf_tabulation_ref_get(tab) == RES_OK);
+ CHK(swf_tabulation_ref_put(NULL) == RES_BAD_ARG);
+ CHK(swf_tabulation_ref_put(tab) == RES_OK);
+ CHK(swf_tabulation_ref_put(tab) == RES_OK);
+
+ args.x_min = SWF_HXD_TABULATE_ARGS_DEFAULT.x_max;
+ args.x_max = SWF_HXD_TABULATE_ARGS_DEFAULT.x_min;
+ CHK(swf_HXd_tabulate(&args, &tab) == RES_BAD_ARG);
+
+ args.x_max = args.x_min;
+ CHK(swf_HXd_tabulate(&args, &tab) == RES_BAD_ARG);
+
+ args.x_min = SWF_HXD_TABULATE_ARGS_DEFAULT.x_min;
+ args.x_max = SWF_HXD_TABULATE_ARGS_DEFAULT.x_max;
+ args.step = 0;
+ CHK(swf_HXd_tabulate(&args, &tab) == RES_BAD_ARG);
+
+ args.step = 1.0e-2;
+ CHK(swf_HXd_tabulate(&args, &tab) == RES_OK);
+ CHK(swf_tabulation_ref_put(tab) == RES_OK);
+}
+
+static void
+Xd(check_inversion)(void)
+{
+ const double nsteps = 10;
+
+ struct swf_H_tabulate_args args = SWF_HXD_TABULATE_ARGS_DEFAULT;
+ struct swf_tabulation* tab = NULL;
+ double pHx = 0;
+ double x0 = 0;
+ size_t i = 0;
+
+ CHK(swf_HXd_tabulate(&args, &tab) == RES_OK);
+
+ for(x0 = args.x_min; x0 <= args.x_max; x0 += x0*args.step) {
+ const double x1 = MMIN(x0 + x0*args.step, args.x_max);
+ const double delta_x = (x1 - x0) / (double)nsteps;
+
+ for(i = 0; i < nsteps; ++i) {
+ const double x = x0 + (double)i*delta_x;
+ const double Hx = swf_HXd_eval(x);
+ const double xl = swf_tabulation_inverse(tab, SWF_LINEAR, Hx);
+ const double xq = swf_tabulation_inverse(tab, SWF_QUADRATIC, Hx);
+ const double xi = swf_HXd_inverse(Hx);
+ const double errl = fabs((x - xl) / x);
+ const double errq = fabs((x - xq) / x);
+ const double erri = fabs((x - xi) / x);
+
+ /* Do not check the value of x if the corresponding H is indistinguishable
+ * from an H calculated from a previous x. In other words, the H's are
+ * numerically equal and any of the corresponding x-values is a valid
+ * inversion result, numerically speaking. */
+ if(Hx == pHx) continue;
+
+ printf("%e %e %e %e %e\n", x, Hx, errq, errl, erri);
+ if(1e-9 < Hx && Hx < (1.0 - 1e-9)) {
+ CHK(errl < 1.e-5);
+ CHK(errq < 1.e-7);
+ CHK(erri < 1.e-7);
+ } else {
+ CHK(errl < 1.e-3);
+ CHK(errq < 1.e-3);
+ CHK(erri < 1.e-3);
+ }
+ pHx = Hx;
+ }
+ }
+
+ CHK(swf_tabulation_ref_put(tab) == RES_OK);
+}
+
+/* Undefine macros used for genericity */
+#undef TEST_SWF_DIMENSION
+#undef SWF_HXD_TABULATE_ARGS_DEFAULT
+#undef swf_HXd_eval
+#undef swf_HXd_inverse
+#undef swf_HXd_tabulate
+#undef Xd
