From 6f83519a1d8593297899d5f456d73df18c5b72e9 Mon Sep 17 00:00:00 2001
From: Markus Hitter <mah@jump-ing.de>
Date: Sun, 6 Jul 2014 17:04:24 +0200
Subject: [PATCH] Add preprocessor math.

For now this is a square root function which should solve entirely
in the preprocessor. Test results described in the file.

Test code for runtime results, inserted right before the main loop
in mendel.c:

  for (uint32_t i = 0; i < 10000000; i++) {
    uint32_t mathlib = (uint32_t)(sqrt(i) + .5);
    uint32_t preprocessor = (uint32_t)(SQRT(i) + .5);

    if (mathlib != preprocessor) {
      sersendf_P(PSTR("%lu: %lu %lu\n"), i, mathlib, preprocessor);
      break;
    }

    if ((i & 0x00001fff) == 0)
      sersendf_P(PSTR("%lu\n"), i);
  }
  sersendf_P(PSTR("Square root check done.\n"));

Test code for compile time results:

  sersendf_P(PSTR("10000000: %lu\n"), (uint32_t)SQRT(10000000));
  sersendf_P(PSTR("10000000: %lu\n"), (uint32_t)sqrt(10000000));
  sersendf_P(PSTR("20000000: %lu\n"), (uint32_t)SQRT(20000000));
  sersendf_P(PSTR("20000000: %lu\n"), (uint32_t)sqrt(20000000));
  sersendf_P(PSTR("30000000: %lu\n"), (uint32_t)SQRT(30000000));
  sersendf_P(PSTR("30000000: %lu\n"), (uint32_t)sqrt(30000000));
  sersendf_P(PSTR("40000000: %lu\n"), (uint32_t)SQRT(40000000));
  sersendf_P(PSTR("40000000: %lu\n"), (uint32_t)sqrt(40000000));
  sersendf_P(PSTR("50000000: %lu\n"), (uint32_t)SQRT(50000000));
  sersendf_P(PSTR("50000000: %lu\n"), (uint32_t)sqrt(50000000));
  sersendf_P(PSTR("60000000: %lu\n"), (uint32_t)SQRT(60000000));
  sersendf_P(PSTR("60000000: %lu\n"), (uint32_t)sqrt(60000000));
  sersendf_P(PSTR("70000000: %lu\n"), (uint32_t)SQRT(70000000));
  sersendf_P(PSTR("70000000: %lu\n"), (uint32_t)sqrt(70000000));
  sersendf_P(PSTR("80000000: %lu\n"), (uint32_t)SQRT(80000000));
  sersendf_P(PSTR("80000000: %lu\n"), (uint32_t)sqrt(80000000));
  sersendf_P(PSTR("90000000: %lu\n"), (uint32_t)SQRT(90000000));
  sersendf_P(PSTR("90000000: %lu\n"), (uint32_t)sqrt(90000000));
---
 preprocessor_math.h | 51 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 51 insertions(+)
 create mode 100644 preprocessor_math.h

diff --git a/preprocessor_math.h b/preprocessor_math.h
new file mode 100644
index 0000000..d4fe60a
--- /dev/null
+++ b/preprocessor_math.h
@@ -0,0 +1,51 @@
+/** \file
+
+  \brief Math functions meant to be calculated in the C preprocessor.
+
+  \details Math functions presented here avoid library calls, which means they
+           can be solved at compile time and as such used to initialise
+           constants. When used for intitialising, their cost at runtime is
+           zero, as they resolve into a single number.
+*/
+
+#ifndef _PREPROCESSOR_MATH_H
+#define _PREPROCESSOR_MATH_H
+
+
+/*! Preprocessor square root.
+
+  (uint32_t)(SQRT(i) + .5)
+  equals
+  (uint32_t)(sqrt(i) + .5)
+
+  These two provide identical results up to 1'861'860 (tested at runtime) and
+  up to 10'000'000 at compile time. At 90'000'000, deviation is about 5%,
+  increasing further at even higher numbers. This "+ .5" is for rounding and
+  not crucial. Casting to other sizes is also possible.
+
+  Can be used for calculations at runtime, too, where it costs 944(!) bytes
+  binary size and takes roughly 10'000(!) clock cycles.
+
+  Initial version found on pl.comp.lang.c, posted by Jean-Louis PATANE.
+*/
+#define SQR00(x) (((double)(x)) / 2)
+#define SQR01(x) ((SQR00(x) + ((x) / SQR00(x))) / 2)
+#define SQR02(x) ((SQR01(x) + ((x) / SQR01(x))) / 2)
+#define SQR03(x) ((SQR02(x) + ((x) / SQR02(x))) / 2)
+#define SQR04(x) ((SQR03(x) + ((x) / SQR03(x))) / 2)
+#define SQR05(x) ((SQR04(x) + ((x) / SQR04(x))) / 2)
+#define SQR06(x) ((SQR05(x) + ((x) / SQR05(x))) / 2)
+#define SQR07(x) ((SQR06(x) + ((x) / SQR06(x))) / 2)
+#define SQR08(x) ((SQR07(x) + ((x) / SQR07(x))) / 2)
+#define SQR09(x) ((SQR08(x) + ((x) / SQR08(x))) / 2)
+#define SQR10(x) ((SQR09(x) + ((x) / SQR09(x))) / 2)
+#define SQR11(x) ((SQR10(x) + ((x) / SQR10(x))) / 2)
+#define SQR12(x) ((SQR11(x) + ((x) / SQR11(x))) / 2)
+// You can add more of these lines here, each additional line increases
+// accurate range by about factor 2 and costs additional 40 bytes binary size
+// in the non-constant case. But beware, the length of the preprocessed term
+// explodes, leading to several seconds compile time above about SQR13.
+#define SQRT(x) ((SQR12(x) + ((x) / SQR12(x))) / 2)
+
+
+#endif /* _PREPROCESSOR_MATH_H */