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

Timestamp:

Sep 22, 2005, 10:45:43 AM (21 years ago)

Author:

gusciora

Message:

Implementing the polynomial nTerm -> nOrder switch,

File:

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

trunk/psLib/src/math/psPolynomial.c

@@ -7 +7 @@
 *  polynomials.  It also contains a Gaussian functions.
 *
-*  @version $Revision: 1.124 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-22 03:00:31 $
+*  @version $Revision: 1.125 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-22 20:45:43 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -298 +298 @@
     /* This is old code that does not use Clenshaw's formula.  Get rid of it.
 
-    psS32 n;
     psS32 i;
     psF32 tmp;
     psPolynomial1D **chebPolys = NULL;
 
-    n = 1 + poly->COOL_1D_n;
-    chebPolys = createChebyshevPolys(n);
+    chebPolys = createChebyshevPolys(1 + poly->COOL_1D_n);
 
     tmp = 0.0;
@@ -360 +358 @@
     // Determine how many Chebyshev polynomials
     // are needed, then create them.
-    maxChebyPoly = 1 + poly->COOL_2D_nX;
-    if ((1 + poly->COOL_2D_nY) > maxChebyPoly) {
-        maxChebyPoly = (1 + poly->COOL_2D_nY);
-    }
-    chebPolys = createChebyshevPolys(maxChebyPoly);
+    maxChebyPoly = poly->COOL_2D_nX;
+    if (poly->COOL_2D_nY > maxChebyPoly) {
+        maxChebyPoly = poly->COOL_2D_nY;
+    }
+    chebPolys = createChebyshevPolys(maxChebyPoly + 1);
 
     for (loop_x = 0; loop_x < (1 + poly->COOL_2D_nX); loop_x++) {
@@ -375 +373 @@
         }
     }
-    for (i=0;i<maxChebyPoly;i++) {
+    for (i=0;i<maxChebyPoly+1;i++) {
         psFree(chebPolys[i]);
     }
@@ -453 +451 @@
     }
 
-    for (i=0;i<maxChebyPoly;i++) {
+    for (i=0;i<maxChebyPoly+1;i++) {
         psFree(chebPolys[i]);
     }
@@ -529 +527 @@
         maxChebyPoly = poly->COOL_4D_nT;
     }
-    chebPolys = createChebyshevPolys(maxChebyPoly);
+    chebPolys = createChebyshevPolys(maxChebyPoly + 1);
 
     for (loop_x = 0; loop_x < (1 + poly->COOL_4D_nX); loop_x++) {
@@ -536 +534 @@
                 for (loop_t = 0; loop_t < (1 + poly->COOL_4D_nZ); loop_t++) {
                     if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) {
+                        //printf("HMMM (%d, %d, %d, %d)\n", loop_x, loop_y, loop_z, loop_t);
+                        //printf("HMMM (%d, %d, %d, %d)\n", chebPolys[loop_x]->COOL_1D_n, chebPolys[loop_y]->COOL_1D_n, chebPolys[loop_z]->COOL_1D_n, chebPolys[loop_t]->COOL_1D_n);
+
                         polySum += poly->coeff[loop_x][loop_y][loop_z][loop_t] *
                                    psPolynomial1DEval(chebPolys[loop_x], x) *
@@ -547 +548 @@
     }
 
-    for (i=0;i<maxChebyPoly;i++) {
+    for (i=0;i<maxChebyPoly+1;i++) {
         psFree(chebPolys[i]);
     }
@@ -618 +619 @@
                                     psPolynomialType type)
 {
-    PS_ASSERT_INT_POSITIVE(n, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(n, NULL);
 
     unsigned int i = 0;
@@ -628 +629 @@
     newPoly->type = type;
     newPoly->COOL_1D_n = n;
-    newPoly->coeff = psAlloc(n * sizeof(psF64));
-    newPoly->coeffErr = psAlloc(n * sizeof(psF64));
+    newPoly->coeff = psAlloc((n + 1) * sizeof(psF64));
+    newPoly->coeffErr = psAlloc((n + 1) * sizeof(psF64));
     newPoly->mask = (psMaskType *)psAlloc((n + 1) * sizeof(psMaskType));
     for (i = 0; i < (n + 1); i++) {