Re: Subnormal cancellation correctness

[Stefan Westerfeld <stefan space twc de>]

   Hi!

On Mon, Apr 03, 2006 at 06:31:48PM +0200, Tim Janik wrote:
> - i plan on writing a set of short standard macros for tests, to be used 
> instead
>   of your g_print/fflush/g_assert. however, i'll adapt your code accordingly
>   once i have the header in a properly includable location.

Agreed. I wrote something similar for aRts, so if you're using C++, you
might get some inspiration if you look at

http://websvn.kde.org/branches/arts/1.5/arts/tests/

Interesting files might be README.test, test.h and as an example on how
to write test code using the framework, testbuffer.cc.

> - there're soem coding style issues, commented in patch.
I fixed those now.

> - if subnormals.cc is really going to be used for long-standing performance
>   tests, it can still be seperated out. so adding correctness tests to it 
>   now
>   makes sense. thanks for tackling.
Ok.

> >+inline double test2d (double v) { return bse_double_zap_denormal (v); }
> >+inline double test3d (double v) { BSE_DOUBLE_FLUSH_with_cond (v); return 
> >v; }
> >+inline double test4d (double v) { BSE_DOUBLE_FLUSH_with_if (v); return v; 
> >}
> >+inline double test5d (double v) { BSE_DOUBLE_FLUSH_with_threshold (v); 
> >return v; }
> 
> i think with these added, you should rename the above to test1f, test2f,...

Ok. This also affected subnormals-aux.cc, but I made sure everything
still works either way.

> >+{
> >+  g_print ("testing algorithm %s for correctness... ", algo_name);
> >+  fflush (stdout);
> >+  const int n = 1000000;
> 
> hm, this is the third time we have n = 1000000 here now, right?
> we should probably make that a global const then, what do you think?

For the correctness tests I used 10^6 iterations (so that they are not
too slow), for the speed tests we use between 10^7 and 10^8 iterations,
so actually its not the same constant everywhere. I have a small
preference for keeping it as it is (it easier to see what the code does
without reading the whole file), but if you really want to change it,
feel free to do so.

I am attaching the patch as I committed it now, so that you can see if
I missed something.

   Cu... Stefan

Index: ChangeLog
===================================================================
RCS file: /cvs/gnome/beast/bse/ChangeLog,v
retrieving revision 1.605
diff -u -p -r1.605 ChangeLog
--- ChangeLog	2 Apr 2006 16:20:29 -0000	1.605
+++ ChangeLog	3 Apr 2006 17:40:15 -0000
@@ -1,3 +1,8 @@
+Mon Apr  3 15:50:14 2006  Stefan Westerfeld <stefan space twc de>
+
+	* tests/subnormals.cc tests/subnormals-aux.cc: Added code which
+	verifies that subnormal elimination functions are correct.
+
 Sun Apr  2 18:19:42 2006  Tim Janik  <timj gtk org>
 
 	* bse/Makefile.am: some more dependency and minor build fixes.
Index: tests/subnormals-aux.cc
===================================================================
RCS file: /cvs/gnome/beast/bse/tests/subnormals-aux.cc,v
retrieving revision 1.3
diff -u -p -r1.3 subnormals-aux.cc
--- tests/subnormals-aux.cc	1 Apr 2006 14:38:02 -0000	1.3
+++ tests/subnormals-aux.cc	3 Apr 2006 17:40:15 -0000
@@ -20,33 +20,33 @@
 #include <bse/bseieee754.h>
 
 float
-test1 (float v)
+test1f (float v)
 {
   return v;
 }
 
 float
-test2 (float v)
+test2f (float v)
 {
   return bse_float_zap_denormal (v);
 }
 
 float
-test3 (float v)
+test3f (float v)
 {
   BSE_FLOAT_FLUSH_with_cond (v);
   return v;
 }
 
 float
-test4 (float v)
+test4f (float v)
 {
   BSE_FLOAT_FLUSH_with_if (v);
   return v;
 }
 
 float
-test5 (float v)
+test5f (float v)
 {
   BSE_FLOAT_FLUSH_with_threshold (v);
   return v;
Index: tests/subnormals.cc
===================================================================
RCS file: /cvs/gnome/beast/bse/tests/subnormals.cc,v
retrieving revision 1.5
diff -u -p -r1.5 subnormals.cc
--- tests/subnormals.cc	1 Apr 2006 14:42:43 -0000	1.5
+++ tests/subnormals.cc	3 Apr 2006 17:40:15 -0000
@@ -1,5 +1,6 @@
 /* BSE - Bedevilled Sound Engine
  * Copyright (C) 2006 Tim Janik
+ * Copyright (C) 2006 Stefan Westerfeld
  *
  * This library is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public
@@ -18,21 +19,67 @@
  */
 #include <bse/bse.h>
 #include <bse/bseieee754.h>
+#include <stdio.h>
 
 #if 1
-inline float    test1 (float v) { return v;     }
-inline float    test2 (float v) { return bse_float_zap_denormal (v); }
-inline float    test3 (float v) { BSE_FLOAT_FLUSH_with_cond (v); return v; }
-inline float    test4 (float v) { BSE_FLOAT_FLUSH_with_if (v); return v; }
-inline float    test5 (float v) { BSE_FLOAT_FLUSH_with_threshold (v); return v; }
+inline float  test1f (float v) { return v;     }
+inline float  test2f (float v) { return bse_float_zap_denormal (v); }
+inline float  test3f (float v) { BSE_FLOAT_FLUSH_with_cond (v); return v; }
+inline float  test4f (float v) { BSE_FLOAT_FLUSH_with_if (v); return v; }
+inline float  test5f (float v) { BSE_FLOAT_FLUSH_with_threshold (v); return v; }
 #else
-extern float    test1 (float v);
-extern float    test2 (float v);
-extern float    test3 (float v);
-extern float    test4 (float v);
-extern float    test5 (float v);
+extern float  test1f (float v);
+extern float  test2f (float v);
+extern float  test3f (float v);
+extern float  test4f (float v);
+extern float  test5f (float v);
 #endif
 
+inline double test2d (double v) { return bse_double_zap_denormal (v); }
+inline double test3d (double v) { BSE_DOUBLE_FLUSH_with_cond (v); return v; }
+inline double test4d (double v) { BSE_DOUBLE_FLUSH_with_if (v); return v; }
+inline double test5d (double v) { BSE_DOUBLE_FLUSH_with_threshold (v); return v; }
+
+template<float Func (float)> void
+test_correct_subnormal_elimination (const char* algo_name)
+{
+  g_print ("testing algorithm %s for correctness... ", algo_name);
+  fflush (stdout);
+  const int n = 1000000;
+  for (int i = 1; i < n; i++)
+    {
+      float value = BSE_FLOAT_MAX_SUBNORMAL * i / n;
+      g_assert (BSE_FLOAT_IS_SUBNORMAL (value));
+
+      float normalized_positive_value = Func (value);
+      g_assert (!BSE_FLOAT_IS_SUBNORMAL (normalized_positive_value));
+
+      float normalized_negative_value = Func (-value);
+      g_assert (!BSE_FLOAT_IS_SUBNORMAL (normalized_negative_value));
+    }
+  g_print ("PASSED\n");
+}
+
+template<double Func (double)> void
+test_correct_subnormal_elimination (const char* algo_name)
+{
+  g_print ("testing algorithm %s for correctness... ", algo_name);
+  fflush (stdout);
+  const int n = 1000000;
+  for (int i = 1; i < n; i++)
+    {
+      double value = BSE_DOUBLE_MAX_SUBNORMAL * i / n;
+      g_assert (BSE_DOUBLE_IS_SUBNORMAL (value));
+
+      double normalized_positive_value = Func (value);
+      g_assert (!BSE_DOUBLE_IS_SUBNORMAL (normalized_positive_value));
+
+      double normalized_negative_value = Func (-value);
+      g_assert (!BSE_DOUBLE_IS_SUBNORMAL (normalized_negative_value));
+    }
+  g_print ("PASSED\n");
+}
+
 int
 main (int   argc,
       char *argv[])
@@ -56,8 +103,8 @@ main (int   argc,
   for (float i = 0; i <= n; i += 1)
     {
       float v = max_sub * i / n;
-      buffer[j++] = test1 (v);
-      buffer[j++] = test1 (-v);
+      buffer[j++] = test1f (v);
+      buffer[j++] = test1f (-v);
       j %= blen;
     }
   volatile_accu += sum;
@@ -70,8 +117,8 @@ main (int   argc,
   for (float i = 0; i <= n; i += 1)
     {
       float v = max_sub * i / n;
-      buffer[j++] = test2 (v);
-      buffer[j++] = test2 (-v);
+      buffer[j++] = test2f (v);
+      buffer[j++] = test2f (-v);
       j %= blen;
     }
   volatile_accu += sum;
@@ -84,8 +131,8 @@ main (int   argc,
   for (float i = 0; i <= n; i += 1)
     {
       float v = max_sub * i / n;
-      buffer[j++] = test3 (v);
-      buffer[j++] = test3 (-v);
+      buffer[j++] = test3f (v);
+      buffer[j++] = test3f (-v);
       j %= blen;
     }
   volatile_accu += sum;
@@ -98,8 +145,8 @@ main (int   argc,
   for (float i = 0; i <= n; i += 1)
     {
       float v = max_sub * i / n;
-      buffer[j++] = test4 (v);
-      buffer[j++] = test4 (-v);
+      buffer[j++] = test4f (v);
+      buffer[j++] = test4f (-v);
       j %= blen;
     }
   volatile_accu += sum;
@@ -112,8 +159,8 @@ main (int   argc,
   for (float i = 0; i <= n; i += 1)
     {
       float v = max_sub * i / n;
-      buffer[j++] = test5 (v);
-      buffer[j++] = test5 (-v);
+      buffer[j++] = test5f (v);
+      buffer[j++] = test5f (-v);
       j %= blen;
     }
   volatile_accu += sum;
@@ -122,6 +169,16 @@ main (int   argc,
 
   g_print ("subnormal cancellation times: keep=%fs zap=%fs inlined-cond=%fs if-cond=%fs arithmetic=%f\n",
            test1_time, test2_time, test3_time, test4_time, test5_time);
+
+  test_correct_subnormal_elimination<test2f> ("zap");
+  test_correct_subnormal_elimination<test3f> ("inlined-cond");
+  test_correct_subnormal_elimination<test4f> ("if-cond");
+  test_correct_subnormal_elimination<test5f> ("arithmetic");
+
+  test_correct_subnormal_elimination<test2d> ("zap-double");
+  test_correct_subnormal_elimination<test3d> ("inlined-cond-double");
+  test_correct_subnormal_elimination<test4d> ("if-cond-double");
+  test_correct_subnormal_elimination<test5d> ("arithmetic-double");
 
   return 0;
 }
-- 
Stefan Westerfeld, Hamburg/Germany, http://space.twc.de/~stefan