star-wf

commit 09198e97681e4bf3ba62851a7f2e2252f6af643e
parent 3df263370df20f7b50d4645992cb41b8bd06f3af
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon, 11 Mar 2024 16:20:21 +0100

Improve the numerical precision of tab H

Do not normalize the function, as it is already normalized. The caller
can force renormalization if he tabulates over a specific range.
Avoiding normalization improves accuracy by two orders of magnitude when
H(x) is close to 1 by less than 1e-14.

Set x_max to 3.337. Above this, H(x) is 1 to an epsilon less than
1.e-14.

Diffstat:
M
src/swf.h
 | 
11
++++++++++-
M
src/swf_H.c
 | 
18
+++++++++---------
M
src/test_swf_H3d.c
 | 
6
+++---

3 files changed, 22 insertions(+), 13 deletions(-)
diff --git a/src/swf.h b/src/swf.h
@@ -51,9 +51,18 @@ struct swf_H_tabulate_args {
    * near 0 */
   double step; /* x_next = x + x*step: */
 
+  /* Force renormalization of the function. Note that H is already normalized to
+   * the correct interval, so normalization may be unnecessary. Worse still! It
+   * could introduce numerical inaccuracy. So don't normalize if you don't
+   * have to. */
+  int normalize;
+
   struct mem_allocator* allocator; /* NULL <=> use default allocator */
 };
-#define SWF_H_TABULATE_ARGS_DEFAULT__ {7.5e-3, 3.8, 1e-3, NULL}
+
+/* The H function's default tab values guarantee both numerical accuracy and a
+ * small memory footprint */
+#define SWF_H_TABULATE_ARGS_DEFAULT__ {7.5e-3, 3.337, 1e-3, 0, NULL}
 static const struct swf_H_tabulate_args SWF_H_TABULATE_ARGS_DEFAULT =
   SWF_H_TABULATE_ARGS_DEFAULT__;
 
diff --git a/src/swf_H.c b/src/swf_H.c
@@ -160,11 +160,8 @@ swf_H3d_tabulate
    struct swf_tabulation** out_tab)
 {
   struct swf_tabulation* tab = NULL;
-  struct item* items = NULL;
-  size_t nitems = 0;
   size_t i = 0;
   double x = 0;
-  double norm = 0;
   res_T res = RES_OK;
 
   if(!out_tab) { res = RES_BAD_ARG; goto error; }
@@ -201,14 +198,17 @@ swf_H3d_tabulate
     res = darray_item_push_back(&tab->items, &item);
     if(res != RES_OK) goto error;
   }
-
   ASSERT(darray_item_size_get(&tab->items) == i);
-  items = darray_item_data_get(&tab->items);
-  nitems = darray_item_size_get(&tab->items);
 
-  /* Normalize the tabulation */
-  norm = items[nitems - 1].fx;
-  FOR_EACH(i, 0, nitems) items[i].fx /= norm;
+
+  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;
diff --git a/src/test_swf_H3d.c b/src/test_swf_H3d.c
@@ -94,10 +94,10 @@ check_tabulation_inversion(void)
       printf("%e %.20e %e %e\n", x, Hx, errq, errl);
       if(1e-9 < Hx && Hx < (1.0 - 1e-9)) {
         CHK(errl < 1.e-5);
-        CHK(errq < 1.e-6);
+        CHK(errq < 1.e-7);
       } else {
-        CHK(errl < 1.e-1);
-        CHK(errq < 1.e-1);
+        CHK(errl < 1.e-3);
+        CHK(errq < 1.e-3);
       }
       pHx = Hx;
     }