+ typedef enum {ADD, SUBTRACT, MULTIPLY, DIVIDE} Optype;
+
+ /** Performs an operation, returning if the result was exact **/
+ // NOTE: DIFFERENT to ParanoidOp (although that wraps to this...)
+ template <class T> bool TrustingOp(T & a, const T & b, Optype op);
+
+ /** Performs an operation _only_ if the result would be exact **/
+ template <class T> bool ParanoidOp(T & a, const T & b, Optype op)
+ {
+ T cpy(a);
+ if (TrustingOp<T>(cpy, b, op))
+ {
+ a = cpy;
+ return true;
+ }
+ return false;
+ }
+
+
+ template <> bool TrustingOp<float>(float & a, const float & b, Optype op);
+ template <> bool TrustingOp<double>(double & a, const double & b, Optype op);
+ template <> bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op);
+
+ // Attempt to comine two terms: a*b + c*d or a/b + c/d
+ template <class T> bool CombineTerms(T & aa, Optype aop, T & bb, T & cc, Optype cop, T & dd)
+ {
+ T a(aa); T b(bb); T c(cc); T d(dd);
+ if (aop == MULTIPLY && cop == MULTIPLY) // a*b + c*d
+ {
+ if ((ParanoidOp<T>(c, b, DIVIDE) || ParanoidOp(d, b, DIVIDE))
+ && TrustingOp<T>(c, d, MULTIPLY) && TrustingOp<T>(a,c,ADD)
+ && TrustingOp<T>(a, b, MULTIPLY)) // (a + (cd)/b) * b
+ {
+ aa = a;
+ bb = 1;
+ cc = 1;
+ dd = 1;
+ return true;
+ }
+ if ((ParanoidOp<T>(a, d, DIVIDE) || ParanoidOp(b, d, DIVIDE))
+ && TrustingOp<T>(a, b, MULTIPLY) && TrustingOp<T>(a,c,ADD)
+ && TrustingOp<T>(a, d, MULTIPLY)) // ((ab)/d + c)*d
+ {
+ aa = a;
+ bb = 1;
+ cc = 1;
+ dd = 1;
+ return true;
+ }
+ return false;
+ }
+ else if (aop == DIVIDE && cop == DIVIDE)
+ {
+ if (TrustingOp<T>(a, d, MULTIPLY) && TrustingOp<T>(c, b, MULTIPLY)
+ && TrustingOp<T>(a, c, ADD) && TrustingOp<T>(b, d, MULTIPLY))
+ {
+ cc = 1;
+ dd = 1;
+ if (ParanoidOp<T>(a, b, DIVIDE))
+ {
+ aa = a;
+ bb = 1;
+ return true;
+ }
+ aa = a;
+ bb = b;
+ return true;
+ }
+ return false;
+ }
+ return false;
+ }
+