Changeset 3988 for trunk/psLib/src/math/psSpline.c

Timestamp:

May 19, 2005, 12:53:47 PM (21 years ago)

Author:

gusciora

Message:

I fixed the bug where the Chebyshev poly had order 1.

File:

• trunk/psLib/src/math/psSpline.c (modified) (6 diffs)

trunk/psLib/src/math/psSpline.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.101 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-05-11 22:02:16 $
+ *  @version $Revision: 1.102 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-05-19 22:53:47 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
  *
- *  XXX: What happens if the polyEval functions are called with data of the wrong
- *       type?
  *  XXX: Should the "coeffErr[]" be used as well?  Bug ???.  Ignore coeffErr
  *
  *  XXX: In the various polyAlloc(n) functions, n is really the order of the
  *  polynomial plus 1.  To create a 2nd-order polynomial, n == 3.
- *
- *  XXX: potential bug: for a multi-dimensional polynomial with order (m, n)
- *  the functions in this file currently do not ignore many of the
- *  coefficients in the coeff matrix that ...
- *
  */
 /*****************************************************************************/
@@ -255 +248 @@
     PS_INT_CHECK_NON_NEGATIVE(maxChebyPoly, NULL);
 
+    printf("HERE: 00\n");
+
     psPolynomial1D **chebPolys = NULL;
 
@@ -269 +264 @@
     if (maxChebyPoly > 1) {
         chebPolys[1]->coeff[1] = 1;
-    }
-    for (psS32 i = 2; i < maxChebyPoly; i++) {
-        for (psS32 j = 0; j < chebPolys[i - 1]->n; j++) {
-            chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j];
-        }
-        for (psS32 j = 0; j < chebPolys[i - 2]->n; j++) {
-            chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
-        }
+
+        for (psS32 i = 2; i < maxChebyPoly; i++) {
+            for (psS32 j = 0; j < chebPolys[i - 1]->n; j++) {
+                chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j];
+            }
+            for (psS32 j = 0; j < chebPolys[i - 2]->n; j++) {
+                chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
+            }
+        }
+    } else {
+        // XXX: Code this.
+        printf("WARNING: %d-order chebyshev polynomials not correctly implemented.\n", maxChebyPoly);
     }
 
@@ -318 +317 @@
 {
     PS_FLOAT_CHECK_RANGE(x, -1.0, 1.0, 0.0);
+    // XXX: Create a macro for this in psConstants.h
+    if (myPoly->n < 1) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial as order %d.", myPoly->n);
+        return(NAN);
+    }
     psVector *d;
-    psS32 n;
+    psS32 n = myPoly->n;
     psS32 i;
-    psF32 tmp;
-
-    n = myPoly->n;
+    psF32 tmp = 0.0;
+
+    // Special case where the Chebyshev poly is constant.
+    if (n == 1) {
+        if (myPoly->mask[0] == 0) {
+            tmp += myPoly->coeff[0];
+        }
+        return(tmp);
+    }
+
+    // Special case where the Chebyshev poly is linear.
+    if (n == 2) {
+        if (myPoly->mask[0] == 0) {
+            tmp+= myPoly->coeff[0];
+        }
+        if (myPoly->mask[1] == 0) {
+            tmp+= myPoly->coeff[1] * x;
+        }
+        return(tmp);
+    }
+
+    // General case where the Chebyshev poly has 2 or more terms.
     d = psVectorAlloc(n, PS_TYPE_F32);
     if(myPoly->mask[n-1] == 0) {
@@ -330 +353 @@
         d->data.F32[n-1] = 0.0;
     }
+
     d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]);
     if(myPoly->mask[n-2] == 0) {
         d->data.F32[n-2] += myPoly->coeff[n-2];
     }
+
     for (i=n-3;i>=1;i--) {
         d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) -
@@ -350 +375 @@
     return(tmp);
 
-    /*
+    /* This is old code that does not use Clenshaw's formula.  Get rid of it.
 
     psS32 n;