Changeset 5813 for trunk/psLib/src/math/psPolynomial.c

Timestamp:

Dec 19, 2005, 1:58:47 PM (20 years ago)

Author:

gusciora

Message:

Implemented the chi-squared algorithm for calculated Chebyshev poly
coefficients.

File:

• trunk/psLib/src/math/psPolynomial.c (modified) (7 diffs)

trunk/psLib/src/math/psPolynomial.c

@@ -7 +7 @@
 *  polynomials.  It also contains a Gaussian functions.
 *
-*  @version $Revision: 1.133 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-11-30 02:00:09 $
+*  @version $Revision: 1.134 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-12-19 23:58:47 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -165 +165 @@
 not nOrder.
  *****************************************************************************/
-static psPolynomial1D **createChebyshevPolys(psS32 maxChebyPoly)
+psPolynomial1D **createChebyshevPolys(psS32 numPolys)
+{
+    PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(numPolys, 1, NULL);
+
+    if (0) {
+        psPolynomial1D **chebPolys = (psPolynomial1D **) psAlloc(numPolys * sizeof(psPolynomial1D *));
+        for (psS32 i = 0; i < numPolys; i++) {
+            chebPolys[i] = psPolynomial1DAlloc(i, PS_POLYNOMIAL_ORD);
+            chebPolys[i]->coeff[i] = 1;
+        }
+        return (chebPolys);
+    } else {
+        psPolynomial1D **chebPolys = (psPolynomial1D **) psAlloc(numPolys * sizeof(psPolynomial1D *));
+        for (psS32 i = 0; i < numPolys; i++) {
+            chebPolys[i] = psPolynomial1DAlloc(i, PS_POLYNOMIAL_ORD);
+        }
+
+        // Create the Chebyshev polynomials.
+        // Polynomial i has i-th order.
+        chebPolys[0]->coeff[0] = 1.0;
+        if (numPolys >= 2) {
+            chebPolys[1]->coeff[1] = 1.0;
+
+            for (psS32 i = 2; i < numPolys; i++) {
+                for (psS32 j = 0; j < chebPolys[i - 1]->nX+1; j++) {
+                    chebPolys[i]->coeff[j + 1] = 2.0 * chebPolys[i - 1]->coeff[j];
+                }
+                for (psS32 j = 0; j < chebPolys[i - 2]->nX+1; j++) {
+                    chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
+                }
+            }
+        }
+
+        return (chebPolys);
+    }
+}
+
+/*
+static psPolynomial1D **createChebyshevPolysOld(psS32 maxChebyPoly)
 {
     PS_ASSERT_INT_NONNEGATIVE(maxChebyPoly, NULL);
-
+ 
     psPolynomial1D **chebPolys = NULL;
-
+ 
     chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
     for (psS32 i = 0; i < maxChebyPoly; i++) {
         chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
     }
-
+ 
     // Create the Chebyshev polynomials.
     // Polynomial i has i-th order.
     chebPolys[0]->coeff[0] = 1;
-
+ 
     // XXX: Bug 296
     if (maxChebyPoly > 1) {
         chebPolys[1]->coeff[1] = 1;
-
+ 
         for (psS32 i = 2; i < maxChebyPoly; i++) {
             for (psS32 j = 0; j < chebPolys[i - 1]->nX; j++) {
@@ -196 +234 @@
         printf("WARNING: %u-order chebyshev polynomials not correctly implemented.\n", maxChebyPoly);
     }
-
+ 
     return (chebPolys);
 }
+*/
 
 /*****************************************************************************
@@ -236 +275 @@
 // XXX: How does the mask vector effect Crenshaw's formula?
 // XXX: We assume that x is scaled between -1.0 and 1.0;
-static psF64 chebPolynomial1DEval(psF64 x,
-                                  const psPolynomial1D* poly)
-{
-    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
-    // XXX: Create a macro for this in psConstants.h
-    if (poly->nX < 1) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial is order %u.", poly->nX);
-        return(NAN);
-    }
+// XXX: Create a faster version for low-order Chebyshevs.
+static psF64 chebPolynomial1DEval(
+    psF64 x,
+    const psPolynomial1D* poly)
+{
+    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, NAN);
+    PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(poly->nX, 0, NAN);
     psVector *d;
+
     // XXX: n should be nTerms here (for clarity).  Or get rid of the variable.
-    unsigned int n = 1 + poly->nX;
+    unsigned int nTerms = 1 + poly->nX;
     unsigned int i;
     psF64 tmp = 0.0;
 
     // Special case where the Chebyshev poly is constant.
-    if (n == 1) {
+    if (nTerms == 1) {
         if (poly->mask[0] == 0) {
             tmp += poly->coeff[0];
@@ -260 +298 @@
 
     // Special case where the Chebyshev poly is linear.
-    if (n == 2) {
+    if (nTerms == 2) {
         if (poly->mask[0] == 0) {
             tmp+= poly->coeff[0];
@@ -270 +308 @@
     }
 
-    // General case where the Chebyshev poly has 2 or more terms.
-    d = psVectorAlloc(n, PS_TYPE_F64);
-    if(poly->mask[n-1] == 0) {
-        d->data.F64[n-1] = poly->coeff[n-1];
+    if (1) {
+        // General case where the Chebyshev poly has 2 or more terms.
+        d = psVectorAlloc(nTerms, PS_TYPE_F64);
+        if(poly->mask[nTerms-1] == 0) {
+            d->data.F64[nTerms-1] = poly->coeff[nTerms-1];
+        } else {
+            d->data.F64[nTerms-1] = 0.0;
+        }
+
+        d->data.F64[nTerms-2] = (2.0 * x * d->data.F64[nTerms-1]);
+        if(poly->mask[nTerms-2] == 0) {
+            d->data.F64[nTerms-2] += poly->coeff[nTerms-2];
+        }
+
+        for (i=nTerms-3;i>=1;i--) {
+            d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) -
+                             (d->data.F64[i+2]);
+            if(poly->mask[i] == 0) {
+                d->data.F64[i] += poly->coeff[i];
+            }
+        }
+
+        tmp = (x * d->data.F64[1]) -
+              (d->data.F64[2]);
+        if(poly->mask[0] == 0) {
+            tmp += (0.5 * poly->coeff[0]);
+        }
+        psFree(d);
     } else {
-        d->data.F64[n-1] = 0.0;
-    }
-
-    d->data.F64[n-2] = (2.0 * x * d->data.F64[n-1]);
-    if(poly->mask[n-2] == 0) {
-        d->data.F64[n-2] += poly->coeff[n-2];
-    }
-
-    for (i=n-3;i>=1;i--) {
-        d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) -
-                         (d->data.F64[i+2]);
-        if(poly->mask[i] == 0) {
-            d->data.F64[i] += poly->coeff[i];
-        }
-    }
-
-    tmp = (x * d->data.F64[1]) -
-          (d->data.F64[2]);
-    if(poly->mask[0] == 0) {
-        tmp += (0.5 * poly->coeff[0]);
-    }
-    psFree(d);
+        // This is old code that does not use Clenshaw's formula.  Get rid of it.
+        psPolynomial1D **chebPolys = createChebyshevPolys(1 + poly->nX);
+
+        tmp = 0.0;
+        for (psS32 i=0;i<(1 + poly->nX);i++) {
+            tmp+= (poly->coeff[i] * psPolynomial1DEval(chebPolys[i], x));
+        }
+        tmp-= (poly->coeff[0]/2.0);
+
+        for (psS32 i=0;i<(1 + poly->nX);i++) {
+            psFree(chebPolys[i]);
+        }
+        psFree(chebPolys);
+    }
+
     return(tmp);
-
-    /* This is old code that does not use Clenshaw's formula.  Get rid of it.
-
-    psS32 i;
-    psF32 tmp;
-    psPolynomial1D **chebPolys = NULL;
-
-    chebPolys = createChebyshevPolys(1 + poly->nX);
-
-    tmp = 0.0;
-    for (i=0;i<(1 + poly->nX);i++) {
-        tmp+= (poly->coeff[i] * psPolynomial1DEval(x, chebPolys[i]));
-    }
-    tmp-= (poly->coeff[0]/2.0);
-
-
-    return(tmp);
-    */
 }
 
@@ -628 +665 @@
     newPoly->type = type;
     newPoly->nX = n;
+    newPoly->p_min = NAN;
+    newPoly->p_max = NAN;
     newPoly->coeff = psAlloc((n + 1) * sizeof(psF64));
     newPoly->coeffErr = psAlloc((n + 1) * sizeof(psF64));