[Boost-commit] svn:boost r75054 - in sandbox/gtl: boost/polygon/detail libs/polygon/test

Subject: [Boost-commit] svn:boost r75054 - in sandbox/gtl: boost/polygon/detail libs/polygon/test
From: sydorchuk.andriy_at_[hidden]
Date: 2011-10-18 16:55:26

Author: asydorchuk
Date: 2011-10-18 16:55:24 EDT (Tue, 18 Oct 2011)
New Revision: 75054
URL: http://svn.boost.org/trac/boost/changeset/75054

Log:
Minor robust_sqrt_expr updates (reduces structure size, optimized performance).
Text files modified: 
   sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp   |    92 ++++++++++++++++++++++++--------------- 
   sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp |     4                                         
   2 files changed, 58 insertions(+), 38 deletions(-)

Modified: sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp	(original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp	2011-10-18 16:55:24 EDT (Tue, 18 Oct 2011)
@@ -514,70 +514,90 @@
     public:
         // Evaluates expression:
         // A[0] * sqrt(B[0]).
-        mpf eval1(mpt *A, mpt *B) {
-            lhs[0] = A[0];
-            rhs[0] = B[0];
-            numer[0] = lhs[0] * sqrt(rhs[0]);
-            return numer[0];
+        mpf& eval1(mpt *A, mpt *B) {
+            a[0] = A[0];
+            b[0] = B[0];
+            b[0] = sqrt(b[0]);
+            return b[0] *= a[0];
         }
 
         // Evaluates expression:
         // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) with
         // 7 * EPS relative error in the worst case.
-        mpf eval2(mpt *A, mpt *B) {
-            lhs[1] = eval1(A, B);
-            rhs[1] = eval1(A + 1, B + 1);
-            if ((lhs[1] >= 0 && rhs[1] >= 0) || (lhs[1] <= 0 && rhs[1] <= 0))
-                return lhs[1] + rhs[1];
-            numer[1] = A[0] * A[0] * B[0] - A[1] * A[1] * B[1];
-            return numer[1] / (lhs[1] - rhs[1]);
+        mpf& eval2(mpt *A, mpt *B) {
+            a[1] = eval1(A, B);
+            b[1] = eval1(A + 1, B + 1);
+            if ((a[1] >= 0 && b[1] >= 0) || (a[1] <= 0 && b[1] <= 0))
+                return b[1] += a[1];
+            for (int i = 0; i < 2; ++i) {
+                temp[i] = A[i] * A[i];
+                temp[i] *= B[i];
+            }
+            a[1] -= b[1];
+            b[1] = temp[0] - temp[1];
+            return b[1] /= a[1];
         }
 
         // Evaluates expression:
         // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2])
         // with 16 * EPS relative error in the worst case.
-        mpf eval3(mpt *A, mpt *B) {
-            lhs[2] = eval2(A, B);
-            rhs[2] = eval1(A + 2, B + 2);
-            if ((lhs[2] >= 0 && rhs[2] >= 0) || (lhs[2] <= 0 && rhs[2] <= 0))
-                return lhs[2] + rhs[2];
-            cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
-            cA[0] -= A[2] * A[2] * B[2];
+        mpf& eval3(mpt *A, mpt *B) {
+            a[2] = eval2(A, B);
+            b[2] = eval1(A + 2, B + 2);
+            if ((a[2] >= 0 && b[2] >= 0) || (a[2] <= 0 && b[2] <= 0))
+                return b[2] += a[2];
+            for (int i = 0; i < 3; ++i) {
+                temp[i] = A[i] * A[i];
+                temp[i] *= B[i];
+            }
+            cA[0] = temp[0] + temp[1];
+            cA[0] -= temp[2];
             cB[0] = 1;
-            cA[1] = A[0] * A[1] * 2;
+            cA[1] = A[0] * A[1];
+            cA[1] += cA[1];
             cB[1] = B[0] * B[1];
-            numer[2] = eval2(cA, cB);
-            return numer[2] / (lhs[2] - rhs[2]);
+            a[2] -= b[2];
+            b[2] = eval2(cA, cB);
+            return b[2] /= a[2];
         }
 
         
         // Evaluates expression:
         // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]) + A[3] * sqrt(B[3])
         // with 25 * EPS relative error in the worst case.
-        mpf eval4(mpt *A, mpt *B) {
-            lhs[3] = eval2(A, B);
-            rhs[3] = eval2(A + 2, B + 2);
-            if ((lhs[3] >= 0 && rhs[3] >= 0) || (lhs[3] <= 0 && rhs[3] <= 0))
-                return lhs[3] + rhs[3];
-            dA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
-            dA[0] -= A[2] * A[2] * B[2] + A[3] * A[3] * B[3];
+        mpf& eval4(mpt *A, mpt *B) {
+            a[3] = eval2(A, B);
+            b[3] = eval2(A + 2, B + 2);
+            if ((a[3] >= 0 && b[3] >= 0) || (a[3] <= 0 && b[3] <= 0))
+                return b[3] += a[3];
+            for (int i = 0; i < 4; ++i) {
+                temp[i] = A[i] * A[i];
+                temp[i] *= B[i];
+            }
+            dA[0] = temp[0] + temp[1];
+            dA[0] -= temp[2];
+            dA[0] -= temp[3];
             dB[0] = 1;
-            dA[1] = A[0] * A[1] * 2;
+            dA[1] = A[0] * A[1];
+            dA[1] += dA[1];
             dB[1] = B[0] * B[1];
-            dA[2] = A[2] * A[3] * -2;
+            dA[2] = A[2] * A[3];
+            dA[2] = -dA[2];
+            dA[2] += dA[2];
             dB[2] = B[2] * B[3];
-            numer[3] = eval3(dA, dB);
-            return numer[3] / (lhs[3] - rhs[3]);
+            a[3] -= b[3];
+            b[3] = eval3(dA, dB);
+            return b[3] /= a[3];
         }
 
     private:
-        mpf lhs[4];
-        mpf rhs[4];
-        mpf numer[4];
+        mpf a[4];
+        mpf b[4];
         mpt cA[2];
         mpt cB[2];
         mpt dA[3];
         mpt dB[3];
+        mpt temp[4];
     };
 
 } // detail

Modified: sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp	(original)
+++ sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp	2011-10-18 16:55:24 EDT (Tue, 18 Oct 2011)
@@ -260,8 +260,8 @@
                     expected_val -= a[j] * sqrt(b[j]);
                 }
             }
-            mpf received_val = (sqrt_expr_.eval4(A, B));
-            ret_val &= almost_equal(expected_val, received_val.get_d(), 1);
+            double received_val = get_d(sqrt_expr_.eval4(A, B));
+            ret_val &= almost_equal(expected_val, received_val, 1);
         }
         return ret_val;
     }