[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[lmi-commits] [lmi] master 284efedf 08/11: For immediate reversion: bina
From: |
Greg Chicares |
Subject: |
[lmi-commits] [lmi] master 284efedf 08/11: For immediate reversion: binary exponentiation |
Date: |
Thu, 26 May 2022 18:14:25 -0400 (EDT) |
branch: master
commit 284efedf7ad41b03c2e40221a6410caf3565b42b
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>
For immediate reversion: binary exponentiation
Both SGI's power() and Knuth's Algorithm A in section 4.6.3 are
presumably well tested, and SGI's documentation implies that they
are equivalent, but that's hard to see at a glance.
Knuth calls this "right-to-left binary method for exponentiation".
The SGI documentation calls this the "Russian Peasant" algorithm,
but that's really something different (although related).
---
sandbox_test.cpp | 94 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 94 insertions(+)
diff --git a/sandbox_test.cpp b/sandbox_test.cpp
index 6baaf4c8..46dcacb2 100644
--- a/sandbox_test.cpp
+++ b/sandbox_test.cpp
@@ -23,7 +23,101 @@
#include "test_tools.hpp"
+#include <cstdio>
+
+// Binary method for exponentiation.
+//
+// See Knuth, TAOCP volume 2, section 4.6.3 (p. 442 in 2nd ed.);
+// and SGI's power(), present elsewhere in the lmi sources.
+//
+// The printf statements that aren't commented out print a table
+// like Knuth's example. Enable the ones that are commented out
+// to see the details of each multiplication. Or :g/printf/d
+// to remove all that clutter.
+//
+// Knuth's algorithm takes one more multiplication for 'Y *= Z'
+// when Y has its initial value of unity. SGI refactors it to
+// avoid a goto, but the result is harder to understand.
+
+#pragma GCC diagnostic ignored "-Wunused-label"
+
+// TAOCP, volume 2, section 4.6.3, page 442
+double Algorithm_A(double x, int n)
+{
+ if(n <= 0) throw "n must be positive";
+ int mult_count {0};
+ A1:
+ int N {n};
+ double Y {1.0};
+ double Z {x};
+ std::printf(" %3s %7s %7s\n" , "N", "Y", "Z");
+ std::printf("After step A1 %3i %7.f %7.f\n", N, Y, Z );
+ A2: // [Halve N.] (At this point, x^n = Y * Z^N .)
+ bool was_even = 0 == N % 2;
+ N /= 2; // integer division truncates
+ if(was_even) goto A5;
+ A3: // [Multiply Y by Z.]
+// std::printf("multiply #%i %7.f by %7.f yielding %7.f\n", mult_count, Y, Z,
Y*Z);
+ Y *= Z;
+ ++mult_count;
+ A4: // [N == 0?]
+ std::printf("After step A4 %3i %7.f %7.f\n", N, Y, Z );
+ if(0 == N)
+ {
+ std::printf("Algorithm A: %i multiplications\n", mult_count);
+ return Y;
+ }
+ A5: // [Square Z.]
+//std::printf("multiply #%i %7.f by %7.f yielding %7.f\n", mult_count, Z, Z,
Z*Z);
+ Z *= Z;
+ ++mult_count;
+ goto A2;
+}
+
+// SGI extension to STL, somewhat refactored for clarity
+double power(double x, int n)
+{
+ if(n < 0)
+ {
+ throw std::logic_error("power() called with negative exponent.");
+ }
+
+ int mult_count {0};
+ while(0 == n % 2) // while ((n & 1) == 0)
+ {
+ n /= 2; // n >>= 1;
+// std::printf("multiply #%i %7.f by %7.f yielding %7.f\n", mult_count, x, x,
x*x);
+ x *= x;
+ ++mult_count;
+std::printf("After step B1 %3i %7.f\n", n, x);
+ }
+ double result = x;
+ n /= 2; // n >>= 1;
+ while (n != 0)
+ {
+// std::printf("multiply #%i %7.f by %7.f yielding %7.f\n", mult_count, x, x,
x*x);
+ x *= x;
+ ++mult_count;
+std::printf("After step B2 %3i %7.f\n", n, x);
+ if(0 != n % 2) // if((n & 1) != 0)
+ {
+// std::printf("multiply #%i %7.f by %7.f yielding %7.f\n", mult_count,
result, x, result*x);
+ result *= x;
+ ++mult_count;
+ }
+ n /= 2; // n >>= 1;
+ }
+std::printf("power(): %i multiplications\n", mult_count);
+ return result;
+}
+
int test_main(int, char*[])
{
+ LMI_TEST(8388608 == Algorithm_A(2.0, 23));
+ LMI_TEST(8388608 == power (2.0, 23));
+
+ LMI_TEST(2 * 8388608 == Algorithm_A(2.0, 24));
+ LMI_TEST(2 * 8388608 == power (2.0, 24));
+
return 0;
}
- [lmi-commits] [lmi] master updated (47f54ead -> 77cc4d65), Greg Chicares, 2022/05/26
- [lmi-commits] [lmi] master 1b9395f2 03/11: Improve documentation, Greg Chicares, 2022/05/26
- [lmi-commits] [lmi] master 699308ab 01/11: Improve documentation, Greg Chicares, 2022/05/26
- [lmi-commits] [lmi] master 086ce6b3 06/11: Fix a probable defect introduced 20210314T1905Z, Greg Chicares, 2022/05/26
- [lmi-commits] [lmi] master 284efedf 08/11: For immediate reversion: binary exponentiation,
Greg Chicares <=
- [lmi-commits] [lmi] master f920e584 10/11: Improve documentation, Greg Chicares, 2022/05/26
- [lmi-commits] [lmi] master 43bbe463 05/11: Use lmi::expm1() and lmi::log1p() in a unit-test function, Greg Chicares, 2022/05/26
- [lmi-commits] [lmi] master 2c0ec871 07/11: Use lmi rather than std functions, Greg Chicares, 2022/05/26
- [lmi-commits] [lmi] master b1b2ed4f 09/11: Revert "For immediate reversion: binary exponentiation", Greg Chicares, 2022/05/26
- [lmi-commits] [lmi] master b6752650 02/11: Improve documentation, Greg Chicares, 2022/05/26
- [lmi-commits] [lmi] master a110f9a0 04/11: Write fdlibm forwarding functions out of line, Greg Chicares, 2022/05/26
- [lmi-commits] [lmi] master 77cc4d65 11/11: Prefer nonstd::power() to std::pow() in a particular case, Greg Chicares, 2022/05/26