From 7efb895ee30b771f487c11fd8cbd50473aef4a6c Mon Sep 17 00:00:00 2001
From: Markus Hitter <mah@jump-ing.de>
Date: Tue, 24 Apr 2012 03:15:36 +0200
Subject: [PATCH] Introduce an integer multiply-divide algorithm.

We have multiplies followed by divides all over the place and
most of them are difficult to handle regarding overflows. This
new algorithm handles this fine in all cases, as long as all
three operators and the overall result fits into 32 bits.
---
 Makefile    |  2 +-
 dda_maths.c | 60 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 dda_maths.h | 11 ++++++++++
 3 files changed, 72 insertions(+), 1 deletion(-)
 create mode 100644 dda_maths.c
 create mode 100644 dda_maths.h

diff --git a/Makefile b/Makefile
index 5391157..d33ccf2 100644
--- a/Makefile
+++ b/Makefile
@@ -80,7 +80,7 @@ PROGID = arduino
 
 PROGRAM = mendel
 
-SOURCES = $(PROGRAM).c dda.c gcode_parse.c gcode_process.c timer.c temp.c sermsg.c dda_queue.c watchdog.c debug.c sersendf.c heater.c analog.c intercom.c pinio.c clock.c home.c crc.c delay.c
+SOURCES = $(PROGRAM).c gcode_parse.c gcode_process.c dda.c dda_maths.c dda_queue.c timer.c temp.c sermsg.c watchdog.c debug.c sersendf.c heater.c analog.c intercom.c pinio.c clock.c home.c crc.c delay.c
 
 ARCH = avr-
 CC = $(ARCH)gcc
diff --git a/dda_maths.c b/dda_maths.c
new file mode 100644
index 0000000..1c6d8db
--- /dev/null
+++ b/dda_maths.c
@@ -0,0 +1,60 @@
+
+/** \file
+	\brief Mathematic algorithms for the digital differential analyser (DDA).
+*/
+
+#include "dda_maths.h"
+
+#include <stdlib.h>
+#include <stdint.h>
+
+/*!
+  Integer multiply-divide algorithm.
+
+  \param multiplicand
+  \param multiplier
+  \param divisor
+  \return rounded result of multiplicand * multiplier / divisor
+
+  Calculate a * b / c, without overflowing and without using 64-bit integers.
+  Doing this the standard way, a * b could easily overflow, even if the correct
+  overall result fits into 32 bits. This algorithm avoids this intermediate
+  overflow and delivers valid results for all cases where each of the three
+  operators as well as the result fits into 32 bits.
+
+  Found on  http://stackoverflow.com/questions/4144232/
+  how-to-calculate-a-times-b-divided-by-c-only-using-32-bit-integer-types-even-i
+*/
+const uint32_t muldiv(uint32_t multiplicand, uint32_t multiplier,
+                      uint32_t divisor) {
+  uint32_t quotient = 0;
+  uint32_t remainder = 0;
+  uint32_t qn = multiplier / divisor;
+  uint32_t rn = multiplier % divisor;
+
+  while(multiplicand) {
+    if (multiplicand & 1) {
+      quotient += qn;
+      remainder += rn;
+      if (remainder >= divisor) {
+        quotient++;
+        remainder -= divisor;
+      }
+    }
+    multiplicand  >>= 1;
+    qn <<= 1;
+    rn <<= 1;
+    if (rn >= divisor) {
+      qn++; 
+      rn -= divisor;
+    }
+  }
+
+  // rounding
+  if (remainder > divisor / 2)
+    quotient++;
+
+  // remainder is valid here, but not returned
+  return quotient;
+}
+
diff --git a/dda_maths.h b/dda_maths.h
new file mode 100644
index 0000000..9efcf7a
--- /dev/null
+++ b/dda_maths.h
@@ -0,0 +1,11 @@
+#ifndef	_DDA_MATHS_H
+#define	_DDA_MATHS_H
+
+#include	<stdint.h>
+
+
+// return rounded result of multiplicand * multiplier / divisor
+const uint32_t muldiv(uint32_t multiplicand, uint32_t multiplier,
+                      uint32_t divisor);
+
+#endif	/* _DDA_MATHS_H */