streamlining Horner's rule functionality across Universal

include/universal/functions/binomial.hpp

@@ -1,7 +1,8 @@
 #pragma once
 // binomial.hpp: definition of recursive binomial coefficient function
 //
-// Copyright (C) 2017-2020 Stillwater Supercomputing, Inc.
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+// SPDX-License-Identifier: MIT
 //
 // This file is part of the universal numbers project, which is released under an MIT Open Source license.
 

include/universal/functions/ddpoly.hpp

@@ -1,7 +1,8 @@
 #pragma once
 // ddpoly.hpp: Evaluate a polynomial of degree N at point x as well as its ND derivatives
 //
-// Copyright (C) 2017-2020 Stillwater Supercomputing, Inc.
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+// SPDX-License-Identifier: MIT
 //
 // This file is part of the universal numbers project, which is released under an MIT Open Source license.
 namespace sw { namespace universal {

include/universal/functions/factorial.hpp

@@ -1,14 +1,15 @@
 #pragma once
 // factorial.hpp: definition of recursive and iterative factorial functions
 //
-// Copyright (C) 2017-2021 Stillwater Supercomputing, Inc.
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+// SPDX-License-Identifier: MIT
 //
 // This file is part of the universal numbers project, which is released under an MIT Open Source license.
 
 namespace sw::function {
 
 	// these factorials can take a Real type and thus could have a very funky behavior
-	// TODO: do we ceil that incoming argument or test on integer properties?
+	// TODO: create a Universal wide is_integral<> type trait and guard the template with it
 
 // factorial implemented using recursion. Should yield reasonable results even for Real types
 // as left-to-right evaluation starts with the smallest values first.

include/universal/functions/functions.hpp

@@ -1,15 +1,18 @@
 #pragma once
 // functions.hpp: aggregation of all Universal number library functions
 //
-// Copyright (C) 2017-2020 Stillwater Supercomputing, Inc.
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+// SPDX-License-Identifier: MIT
 //
 // This file is part of the universal numbers project, which is released under an MIT Open Source license.
 
 // properties of a number
-#include "isrepresentable.hpp"
-#include "twosum.hpp"
+#include <universal/functions/isrepresentable.hpp>
 
 // special functions
-#include "factorial.hpp"
-#include "binomial.hpp"
-#include "loss.hpp"
+#include <universal/functions/factorial.hpp>
+#include <universal/functions/binomial.hpp>
+#include <universal/functions/loss.hpp>
+
+// generic mathematical functions
+#include <universal/functions/horners.hpp>

include/universal/functions/horners.hpp

@@ -0,0 +1,78 @@
+#pragma once
+// horners.hpp: generic Horner's polynomial evaluation and root finding functions 
+//
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+// SPDX-License-Identifier: MIT
+//
+// This file is part of the universal numbers project, which is released under an MIT Open Source license.
+
+namespace sw { namespace universal {
+
+    /// <summary>
+    /// polyeval evaluates a given n-th degree polynomial at x using Horner's rule.
+    /// The polynomial is given by the array of (n+1) coefficients.
+    /// </summary>
+    /// <param name="c">polynomial coefficients</param>
+    /// <param name="n">portion of the polynomial to evaluate</param>
+    /// <param name="x">value to evaluate</param>
+    /// <returns>polynomial at x</returns>
+    template<typename Vector, typename Scalar>
+    inline Scalar polyeval(const Vector& coefficients, int n, const Scalar& x) {
+        // Horner's method of polynomial evaluation
+        assert(coefficients.size() > n);
+        Scalar r{ coefficients[n] };
+
+        for (int i = n - 1; i >= 0; --i) {
+            r *= x;
+            r += coefficients[i];
+        }
+
+        return r;
+    }
+
+    /// <summary>
+    /// polyroot finds a root close to the initial guess x0.
+    /// Will only find a single root as it is using a Newton iteration
+    /// </summary>
+    /// <param name="c">n-th degree polynomial</param>
+    /// <param name="x0">initial guess of the root of interest</param>
+    /// <param name="max_iter">stopping criterium</param>
+    /// <param name="threshold">stopping cirterium</param>
+    /// <returns>root closest to initial guess x0</returns>
+    template<typename Vector, typename Scalar>
+    inline Scalar polyroot(const Vector& c, const Scalar& x0, int max_iter, double threshold = 1e-16) {
+        if (threshold == 0.0) threshold = std::numeric_limits<Scalar>::epsilon();
+
+        int n = c.size() - 1;
+        double max_c = std::abs(double(c[0]));
+        // Compute the coefficients of the derivatives
+        Vector derivatives{ c };
+        for (int i = 1; i <= n; ++i) {
+            double v = std::abs(double(c[i]));
+            if (v > max_c) max_c = v;
+            derivatives[i - 1] = c[i] * static_cast<double>(i);
+        }
+        threshold *= max_c;
+
+        // Newton iteration
+        bool converged{ false };
+        Scalar x = x0;
+        for (int i = 0; i < max_iter; ++i) {
+            Scalar f = polyeval(c, n, x);
+
+            if (abs(f) < threshold) {
+                converged = true;
+                break;
+            }
+            x -= (f / polyeval(derivatives, n - 1, x));
+        }
+
+        if (!converged) {
+           std::cerr << "polyroot: failed to converge\n";
+            return Scalar(SpecificValue::snan);
+        }
+
+        return x;
+    }
+
+}} // namespace sw::universal

include/universal/functions/isrepresentable.hpp

@@ -1,7 +1,8 @@
 #pragma once
 // isrepresentable.hpp: test to see if a ratio of Integer type can be represented by a binary real
 //
-// Copyright (C) 2017-2021 Stillwater Supercomputing, Inc.
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+// SPDX-License-Identifier: MIT
 //
 // This file is part of the universal numbers project, which is released under an MIT Open Source license.
 #include <string>

include/universal/functions/lerp.hpp

@@ -1,7 +1,8 @@
 #pragma once
 // lerp.hpp: definition of a linear interpolation function
 //
-// Copyright (C) 2017-2022 Stillwater Supercomputing, Inc.
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+// SPDX-License-Identifier: MIT
 //
 // This file is part of the universal numbers project, which is released under an MIT Open Source license.
 

include/universal/functions/loss.hpp

@@ -1,7 +1,8 @@
 #pragma once
 // loss.hpp: definition of different loss functions
 //
-// Copyright (C) 2017-2022 Stillwater Supercomputing, Inc.
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+// SPDX-License-Identifier: MIT
 //
 // This file is part of the universal numbers project, which is released under an MIT Open Source license.
 

include/universal/number/cfloat/cfloat.hpp

@@ -72,6 +72,7 @@
 ///////////////////////////////////////////////////////////////////////////////////////
 /// elementary math functions library
 #include <universal/number/cfloat/mathlib.hpp>
+#include <universal/number/cfloat/mathext.hpp>
 
 ///////////////////////////////////////////////////////////////////////////////////////
 /// aliases for industry standard floating point configurations

include/universal/number/cfloat/mathext.hpp

@@ -0,0 +1,9 @@
+#pragma once
+// mathext.hpp: definition of mathematical extension functions for the classic floats
+//
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+// SPDX-License-Identifier: MIT
+//
+// This file is part of the universal numbers project, which is released under an MIT Open Source license.
+
+#include <universal/number/cfloat/mathext/horners.hpp>

include/universal/number/cfloat/mathext/horners.hpp

@@ -0,0 +1,77 @@
+#pragma once
+// horners.hpp: Horner's polynomial evaluation and root finding functions for double-double (dd) floating-point
+//
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+// SPDX-License-Identifier: MIT
+//
+// This file is part of the universal numbers project, which is released under an MIT Open Source license.
+#include <vector>
+#include <universal/number/cfloat/cfloat_fwd.hpp>
+
+namespace sw { namespace universal {
+
+    /// <summary>
+    /// polyeval evaluates a given n-th degree polynomial at x using Horner's rule.
+    /// The polynomial is given by the array of (n+1) coefficients.
+    /// </summary>
+    /// <param name="c">polynomial coefficients</param>
+    /// <param name="n">portion of the polynomial to evaluate</param>
+    /// <param name="x">value to evaluate</param>
+    /// <returns>polynomial at x</returns>
+    template<unsigned nbits, unsigned es, typename BlockType, bool hasSubnormals, bool hasSupernormals, bool isSaturating>
+    inline cfloat<nbits, es, BlockType, hasSubnormals, hasSupernormals, isSaturating> polyeval(const std::vector<cfloat<nbits, es, BlockType, hasSubnormals, hasSupernormals, isSaturating>>& coefficients, int n, const cfloat<nbits, es, BlockType, hasSubnormals, hasSupernormals, isSaturating>& x) {
+        // Horner's method of polynomial evaluation
+        assert(coefficients.size() > n);
+        cfloat<nbits, es, BlockType, hasSubnormals, hasSupernormals, isSaturating> r = coefficients[n];
+
+        for (int i = n - 1; i >= 0; --i) {
+            r *= x;
+            r += coefficients[i];
+        }
+
+        return r;
+    }
+
+    /* polyroot(c, n, x0)
+       Given an n-th degree polynomial, finds a root close to
+       the given guess x0.  Note that this uses simple Newton
+       iteration scheme, and does not work for multiple roots.  */
+
+    template<unsigned nbits, unsigned es, typename BlockType, bool hasSubnormals, bool hasSupernormals, bool isSaturating>
+    inline cfloat<nbits, es, BlockType, hasSubnormals, hasSupernormals, isSaturating> polyroot(const std::vector<cfloat<nbits, es, BlockType, hasSubnormals, hasSupernormals, isSaturating>>& c, const cfloat<nbits, es, BlockType, hasSubnormals, hasSupernormals, isSaturating>& x0, int n, int max_iter, double threshold) {
+	    using Cfloat = cfloat<nbits, es, BlockType, hasSubnormals, hasSupernormals, isSaturating>;
+        if (threshold == 0.0) threshold = std::numeric_limits<Cfloat>::epsilon();
+
+        int n = c.size() - 1;
+        double max_c = std::abs(double(c[0]));
+        // Compute the coefficients of the derivatives
+        std::vector<Cfloat> derivatives{ c };
+        for (int i = 1; i <= n; i++) {
+            double v = std::abs(double(c[i]));
+            if (v > max_c) max_c = v;
+            derivatives[i - 1] = c[i] * static_cast<double>(i);
+        }
+        threshold *= max_c;
+
+        // Newton iteration
+        bool converged{ false };
+        Cfloat x = x0;
+        for (int i = 0; i < max_iter; i++) {
+            Cfloat f = polyeval(c, n, x);
+
+            if (abs(f) < threshold) {
+                converged = true;
+                break;
+            }
+            x -= (f / polyeval(derivatives, n - 1, x));
+        }
+
+        if (!converged) {
+           std::cerr << "polyroot: failed to converge\n";
+            return Cfloat(SpecificValue::snan);
+        }
+
+        return x;
+    }
+
+}} // namespace sw::universal

include/universal/number/dd/dd_fwd.hpp

@@ -6,6 +6,7 @@
 //
 // This file is part of the universal numbers project, which is released under an MIT Open Source license.
 #include <cstdint>
+#include <string>
 
 namespace sw { namespace universal {
 

include/universal/number/dd/mathext/horners.hpp

@@ -25,47 +25,50 @@ namespace sw { namespace universal {
     inline dd polyeval(const std::vector<dd>& coefficients, int n, const dd& x) {
         // Horner's method of polynomial evaluation
         assert(coefficients.size() > n);
-        dd r = coefficients[n];
+        dd r{ coefficients[n] };
 
-        for (int i = n - 1; i >= 0; i--) {
+        for (int i = n - 1; i >= 0; --i) {
             r *= x;
             r += coefficients[i];
         }
 
         return r;
     }
 
-    /* polyroot(c, n, x0)
-       Given an n-th degree polynomial, finds a root close to
-       the given guess x0.  Note that this uses simple Newton
-       iteration scheme, and does not work for multiple roots.  */
-    inline dd polyroot(const std::vector<dd>& c, const dd& x0, int n, int max_iter, double threshold) {
-        dd x = x0;
-
-        std::vector<dd> d{ c };
-        bool converged{ false };
-        double max_c = std::abs(double(c[0]));
-        double v;
-
+    /// <summary>
+    /// polyroot finds a root close to the initial guess x0.
+    /// Will only find a single root as it is using a Newton iteration
+    /// </summary>
+    /// <param name="c">n-th degree polynomial</param>
+    /// <param name="x0">initial guess of the root of interest</param>
+    /// <param name="max_iter">stopping criterium</param>
+    /// <param name="threshold">stopping cirterium</param>
+    /// <returns></returns>
+    inline dd polyroot(const std::vector<dd>& c, const dd& x0, int max_iter, double threshold = 1e-16) {
         if (threshold == 0.0) threshold = dd_eps;
 
+        int n = c.size() - 1;
+        double max_c = std::abs(double(c[0]));
         // Compute the coefficients of the derivatives
-        for (int i = 1; i <= n; i++) {
-            v = std::abs(double(c[i]));
+        std::vector<dd> derivatives{ c };
+        for (int i = 1; i <= n; ++i) {
+            double v = std::abs(double(c[i]));
             if (v > max_c) max_c = v;
-            d[i - 1] = c[i] * static_cast<double>(i);
+            derivatives[i - 1] = c[i] * static_cast<double>(i);
         }
         threshold *= max_c;
 
         // Newton iteration
-        for (int i = 0; i < max_iter; i++) {
+        bool converged{ false };
+        dd x = x0;
+        for (int i = 0; i < max_iter; ++i) {
             dd f = polyeval(c, n, x);
 
             if (abs(f) < threshold) {
                 converged = true;
                 break;
             }
-            x -= (f / polyeval(d, n - 1, x));
+            x -= (f / polyeval(derivatives, n - 1, x));
         }
 
         if (!converged) {