fixes

pow 0.5 to sqrt

diff --git a/src/spirals.cpp b/src/spirals.cpp
index d051facd8aee0e4b6f48f3dcc3057bb1cec29d78..f74a8d6f1748d1fbeb2bdf7e4bc807ccc22b3abf 100644 (file)
--- a/src/spirals.cpp
+++ b/src/spirals.cpp
@@ -54,10 +54,10 @@ long double Fx (long double x0, long double y0, long double z0, long double p1,
     long double ret = 0;
     for (int i = 0; i < x.size(); i++) {
         ret +=
-           std::sqrt (  pow (p2 * (x[i][2] - z0) - p3 * (x[i][1] - y0), 2) +
+           sqrt (   pow (p2 * (x[i][2] - z0) - p3 * (x[i][1] - y0), 2) +
                         pow (p3 * (x[i][0] - x0) - p1 * (x[i][2] - z0), 2) +
                         pow (p1 * (x[i][1] - y0) - p2 * (x[i][0] - x0), 2)) /
-           std::sqrt (p1 * p1 + p2 * p2 + p3 * p3);
+           sqrt (p1 * p1 + p2 * p2 + p3 * p3);
     }
     return ret;
 }
@@ -67,10 +67,10 @@ long double fx0 (long double x0, long double y0, long double z0, long double p1,
     for (int i = 0; i < x.size(); i++) {
         ret +=
            (2 * p2 * (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1])) + 2 * p3 * (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]))) /
-           (2 *     pow (   pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
+           (2 *     sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
-                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2), 0.5) *
-                    pow (p1 * p1 + p2 * p2 + p3 * p3, 0.5));
+                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
+                    sqrt (p1 * p1 + p2 * p2 + p3 * p3));
     }
     return ret;
 }
@@ -80,10 +80,10 @@ long double fy0 (long double x0, long double y0, long double z0, long double p1,
     for (int i = 0; i < x.size(); i++) {
         ret +=
            -(2 * p1 * (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1])) - 2 * p3 * (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]))) /
-            (2 *    pow (   pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
+            (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
-                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2), 0.5) *
-                    pow (p1 * p1 + p2 * p2 + p3 * p3, 0.5));
+                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
+                    sqrt (p1 * p1 + p2 * p2 + p3 * p3));
     }
     return ret;
 }
@@ -93,10 +93,10 @@ long double fz0 (long double x0, long double y0, long double z0, long double p1,
     for (int i = 0; i < x.size(); i++) {
         ret +=
            -(2 * p1 * (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2])) + 2 * p2 * (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]))) /
-            (2 *    pow (   pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
+            (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
-                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2), 0.5) *
-                    pow (p1 * p1 + p2 * p2 + p3 * p3, 0.5));
+                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
+                    sqrt (p1 * p1 + p2 * p2 + p3 * p3));
     }
     return ret;
 }
@@ -106,13 +106,13 @@ long double fp1 (long double x0, long double y0, long double z0, long double p1,
     for (int i = 0; i < x.size(); i++) {
         ret +=
            -(2 * (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1])) * (y0 - x[i][1]) + 2 * (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2])) * (z0 - x[i][2])) /
-            (2 *    pow (   pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
+            (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
-                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2), 0.5) *
-                    pow (p1 * p1 + p2 * p2 + p3 * p3, 0.5)) -
-            (p1 *   pow (   pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
+                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
+                    sqrt (p1 * p1 + p2 * p2 + p3 * p3)) -
+            (p1 *   sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
                             pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
-                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2), 0.5)) /
+                            pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2))) /
                     pow (p1 * p1 + p2 * p2 + p3 * p3, 1.5);
     }
     return ret;
@@ -123,13 +123,13 @@ long double fp2 (long double x0, long double y0, long double z0, long double p1,
     for (int i = 0; i < x.size(); i++) {
         ret +=
            (2 * (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1])) * (x0 - x[i][0]) - 2 * (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2])) * (z0 - x[i][2])) /
-           (2 *    pow (   pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
+           (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
                            pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
-                           pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2), 0.5) *
-                   pow (p1 * p1 + p2 * p2 + p3 * p3, 0.5)) -
-           (p2 *   pow (   pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
+                           pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
+                   sqrt (p1 * p1 + p2 * p2 + p3 * p3)) -
+           (p2 *   sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
                            pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
-                           pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2), 0.5)) /
+                           pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2))) /
                    pow (p1 * p1 + p2 * p2 + p3 * p3, 1.5);
     }
     return ret;
@@ -140,13 +140,13 @@ long double fp3 (long double x0, long double y0, long double z0, long double p1,
     for (int i = 0; i < x.size(); i++) {
         ret +=
            (2 * (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2])) * (x0 - x[i][0]) + 2 * (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2])) * (y0 - x[i][1])) /
-           (2 *    pow (   pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
+           (2 *    sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
                            pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
-                           pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2), 0.5) *
-                   pow (p1 * p1 + p2 * p2 + p3 * p3, 0.5)) -
-           (p3 *   pow (   pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
+                           pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2)) *
+                   sqrt (p1 * p1 + p2 * p2 + p3 * p3)) -
+           (p3 *   sqrt (  pow (p2 * (x0 - x[i][0]) - p1 * (y0 - x[i][1]), 2) +
                            pow (p3 * (x0 - x[i][0]) - p1 * (z0 - x[i][2]), 2) +
-                           pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2), 0.5)) /
+                           pow (p3 * (y0 - x[i][1]) - p2 * (z0 - x[i][2]), 2))) /
                    pow (p1 * p1 + p2 * p2 + p3 * p3, 1.5);
     }
     return ret;
@@ -171,42 +171,42 @@ void linear_kernel_search (long double &x0, long double &y0, long double &z0, lo
         while (Fx(x0 - L1 * FX0, y0, z0, p1, p2, p3, x) - FX >= 0) {
             L1 /= 2;
             if (x0 - L1 * FX0 < epsi) {
-                break;
+                L1 = 0;
             }
         }
         L2 = 1;
         while (Fx(x0, y0 - L2 * FY0, z0, p1, p2, p3, x) - FX >= 0) {
             L2 /= 2;
             if (y0 - L2 * FY0 < epsi) {
-                break;
+                L2 = 0;
             }
         }
         L3 = 1;
         while (Fx(x0, y0, z0 - L3 * FZ0, p1, p2, p3, x) - FX >= 0) {
             L3 /= 2;
             if (z0 - L3 * FZ0 < epsi) {
-                break;
+                L3 = 0;
             }
         }
         L4 = 1;
         while (Fx(x0, y0, z0, p1 - L4 * FP1, p2, p3, x) - FX >= 0) {
             L4 /= 2;
             if (p1 - L4 * FP1 < epsi) {
-                break;
+                L4 = 0;
             }
         }
         L5 = 1;
         while (Fx(x0, y0, z0, p1, p2 - L5 * FP2, p3, x) - FX >= 0) {
             L5 /= 2;
             if (p2 - L5 * FP2 < epsi) {
-                break;
+                L5 = 0;
             }
         }
         L6 = 1;
         while (Fx(x0, y0, z0, p1, p2, p3 - L6 * FP3, x) - FX >= 0) {
             L6 /= 2;
             if (p3 - L6 * FP3 < epsi) {
-                break;
+                L6 = 0;
             }
         }
         x0 -= L1 * FX0;