Changeset 2394Tweet

Timestamp:

Feb 11, 2013, 4:27:49 PM (10 years ago)

Author:

sam

Message:

math: minor improvements to the Remez exchange algorithm.

Location:

trunk/src

Files:

• lol/math/remez.h (modified) (11 diffs)
• math/real.cpp (modified) (4 diffs)

trunk/src/lol/math/remez.h

@@ -57 +57 @@
             printf("Step %i error: ", n);
             newerror.print(m_decimals);
+            printf("\n");
 
             Step();
@@ -77 +78 @@
     }
 
-    T ChebyEval(T const &x)
+    T EvalCheby(T const &x)
     {
         T ret = 0.0, xn = 1.0;
@@ -100 +101 @@
         {
             zeroes[i] = (T)(2 * i - ORDER) / (T)(ORDER + 1);
-            fxn[i] = Value(zeroes[i]);
+            fxn[i] = EvalFunc(zeroes[i]);
         }
 
@@ -146 +147 @@
 
             left.value = control[i];
-            left.error = ChebyEval(left.value) - Value(left.value);
+            left.error = EvalCheby(left.value) - EvalFunc(left.value);
             right.value = control[i + 1];
-            right.error = ChebyEval(right.value) - Value(right.value);
+            right.error = EvalCheby(right.value) - EvalFunc(right.value);
 
             static T limit = ldexp((T)1, -500);
@@ -155 +156 @@
             {
                 mid.value = (left.value + right.value) / 2;
-                mid.error = ChebyEval(mid.value) - Value(mid.value);
-
-                if ((left.error < zero && mid.error < zero)
-                     || (left.error > zero && mid.error > zero))
+                mid.error = EvalCheby(mid.value) - EvalFunc(mid.value);
+
+                if ((left.error <= zero && mid.error <= zero)
+                     || (left.error >= zero && mid.error >= zero))
                     left = mid;
                 else
@@ -185 +186 @@
                 b = zeroes[i];
 
-            T maxerror = 0, maxweight = 0;
-            int best = -1;
-
             for (int round = 0; ; round++)
             {
+                T maxerror = 0, maxweight = 0;
+                int best = -1;
+
                 T c = a, delta = (b - a) / 4;
                 for (int k = 0; k <= 4; k++)
@@ -195 +196 @@
                     if (round == 0 || (k & 1))
                     {
-                        T error = ChebyEval(c) - Value(c);
-                        T weight = Weight(c);
-                        if (fabs(error * maxweight) >= fabs(maxerror * weight))
+                        T error = fabs(EvalCheby(c) - EvalFunc(c));
+                        T weight = fabs(Weight(c));
+                        /* if error/weight >= maxerror/maxweight */
+                        if (error * maxweight >= maxerror * weight)
                         {
                             maxerror = error;
@@ -218 +220 @@
                     b = a + delta * (best + 1);
                     a = a + delta * (best - 1);
-                    best = 2;
                     break;
                 }
@@ -224 +225 @@
                 if (delta < m_epsilon)
                 {
-                    T e = maxerror / maxweight;
+                    T e = fabs(maxerror / maxweight);
                     if (e > final)
                         final = e;
@@ -241 +242 @@
         T fxn[ORDER + 2];
         for (int i = 0; i < ORDER + 2; i++)
-            fxn[i] = Value(control[i]);
+            fxn[i] = EvalFunc(control[i]);
 
         /* We build a matrix of Chebishev evaluations: row i contains the
@@ -335 +336 @@
         }
 
-        printf("Polynomial estimate:\n");
+        printf("Polynomial estimate: ");
         for (int j = 0; j < ORDER + 1; j++)
         {
             if (j)
-                printf("+");
-            printf("x**%i*", j);
+                printf(" + x**%i * ", j);
             an[j].print(m_decimals);
         }
-        printf("\n");
-    }
-
-    T Value(T const &x)
+        printf("\n\n");
+    }
+
+    T EvalFunc(T const &x)
     {
         return m_func(x * m_k2 + m_k1);

trunk/src/math/real.cpp

@@ -1362 +1362 @@
     for (int i = 0; i < BIGITS; i++)
         std::printf(" %08x", m_mantissa[i]);
-    std::printf("\n");
 }
 
@@ -1371 +1370 @@
     char *buf = new char[ndigits + 32 + 10];
     real::sprintf(buf, ndigits);
-    std::printf("%s\n", buf);
+    std::printf("%s", buf);
     delete[] buf;
 }
@@ -1387 +1386 @@
     if (!x)
     {
-        std::strcpy(str, "0.0\n");
+        std::strcpy(str, "0.0");
         return;
     }
@@ -1415 +1414 @@
     /* Print exponent information */
     if (exponent)
-        str += std::sprintf(str, "e%c%i", exponent > 0 ? '+' : '-', -exponent);
+        str += std::sprintf(str, "e%c%i",
+                            exponent >= 0 ? '+' : '-', lol::abs(exponent));
 
     *str++ = '\0';