Index: trunk/psLib/src/math/psPolynomial.c =================================================================== --- trunk/psLib/src/math/psPolynomial.c (revision 6101) +++ trunk/psLib/src/math/psPolynomial.c (revision 6186) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * -* @version $Revision: 1.135 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-01-21 02:43:31 $ +* @version $Revision: 1.136 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-01-23 22:25:31 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -14,6 +14,4 @@ * 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. */ /*****************************************************************************/ @@ -154,5 +152,5 @@ /***************************************************************************** -createChebyshevPolys(n): this routine takes as input the required order n, +p_psCreateChebyshevPolys(n): this routine takes as input the required order n, and returns as output as a pointer to an array of n psPolynomial1D structures, corresponding to the first n Chebyshev polynomials. @@ -162,8 +160,6 @@ outer coefficients of the Chebyshev polynomials. -XXX: From the poly nOrder/nTerm changem the maxChebyPoly here is still nTerms, -not nOrder. *****************************************************************************/ -psPolynomial1D **createChebyshevPolys(psS32 numPolys) +psPolynomial1D **p_psCreateChebyshevPolys(psS32 numPolys) { PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(numPolys, 1, NULL); @@ -202,46 +198,11 @@ } -/* -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++) { - chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j]; - } - for (psS32 j = 0; j < chebPolys[i - 2]->nX; j++) { - chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j]; - } - } - } else { - // XXX: Code this. - printf("WARNING: %u-order chebyshev polynomials not correctly implemented.\n", maxChebyPoly); - } - - return (chebPolys); -} -*/ /***************************************************************************** Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ -static psF64 ordPolynomial1DEval(psF64 x, - const psPolynomial1D* poly) +static psF64 ordPolynomial1DEval( + psF64 x, + const psPolynomial1D* poly) { unsigned int loop_x = 0; @@ -249,11 +210,8 @@ psF64 xSum = 1.0; - psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4, - "---- Calling ordPolynomial1DEval(%lf)\n", x); - psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4, - "Polynomial order is %u\n", poly->nX); + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + psTrace(__func__, 4, "Polynomial order is %u\n", poly->nX); for (loop_x = 0; loop_x < poly->nX+1; loop_x++) { - psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4, - "Polynomial coeff[%u] is %lf\n", loop_x, poly->coeff[loop_x]); + psTrace(__func__, 4, "Polynomial coeff[%u] is %lf\n", loop_x, poly->coeff[loop_x]); } @@ -262,5 +220,5 @@ // XXX: If you set the tracelevel to 10 here, and later set the tracelevel to // 2 or higher in the test code, you get seg faults. - psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 8, + psTrace(__func__, 8, "polysum+= sum*coeff [%lf+= (%lf * %lf)\n", polySum, xSum, poly->coeff[loop_x]); polySum += xSum * poly->coeff[loop_x]; @@ -269,4 +227,5 @@ } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); return(polySum); } @@ -284,5 +243,4 @@ psVector *d; - // XXX: n should be nTerms here (for clarity). Or get rid of the variable. unsigned int nTerms = 1 + poly->nX; unsigned int i; @@ -338,5 +296,5 @@ } else { // This is old code that does not use Clenshaw's formula. Get rid of it. - psPolynomial1D **chebPolys = createChebyshevPolys(1 + poly->nX); + psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(1 + poly->nX); tmp = 0.0; @@ -402,5 +360,5 @@ maxChebyPoly = poly->nY; } - chebPolys = createChebyshevPolys(maxChebyPoly + 1); + chebPolys = p_psCreateChebyshevPolys(maxChebyPoly + 1); for (loop_x = 0; loop_x < (1 + poly->nX); loop_x++) { @@ -476,5 +434,5 @@ maxChebyPoly = poly->nZ; } - chebPolys = createChebyshevPolys(maxChebyPoly + 1); + chebPolys = p_psCreateChebyshevPolys(maxChebyPoly + 1); for (loop_x = 0; loop_x < (1 + poly->nX); loop_x++) { @@ -567,6 +525,6 @@ maxChebyPoly = poly->nT; } - // XXX: Add 1 since createChebyshevPolys() takes nTerms, not nOrder. - chebPolys = createChebyshevPolys(maxChebyPoly + 1); + // Add 1 since p_psCreateChebyshevPolys() takes nTerms, not nOrder. + chebPolys = p_psCreateChebyshevPolys(maxChebyPoly + 1); for (loop_x = 0; loop_x < (1 + poly->nX); loop_x++) { @@ -648,27 +606,22 @@ /***************************************************************************** This routine must allocate memory for the polynomial structures. - XXX: replaces nterms variables with nOrder. Lots of potential for bugs - here. Probably should create separately named private variables - in these functions. *****************************************************************************/ -psPolynomial1D* psPolynomial1DAlloc(unsigned int n, - psPolynomialType type) +psPolynomial1D* psPolynomial1DAlloc( + unsigned int n, + psPolynomialType type) { PS_ASSERT_INT_NONNEGATIVE(n, NULL); - - unsigned int i = 0; - psPolynomial1D* newPoly = NULL; - - newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); + psU32 nOrder = n; + psPolynomial1D *newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); psMemSetDeallocator(newPoly, (psFreeFunc) polynomial1DFree); newPoly->type = type; - newPoly->nX = n; + newPoly->nX = nOrder; newPoly->p_min = NAN; newPoly->p_max = NAN; - 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++) { + newPoly->coeff = psAlloc((nOrder + 1) * sizeof(psF64)); + newPoly->coeffErr = psAlloc((nOrder + 1) * sizeof(psF64)); + newPoly->mask = (psMaskType *)psAlloc((nOrder + 1) * sizeof(psMaskType)); + for (psU32 i = 0; i < (nOrder + 1); i++) { newPoly->coeff[i] = 0.0; newPoly->coeffErr[i] = 0.0; @@ -967,9 +920,9 @@ // XXX: The output of this routine is always psF64 while 1D and 2D are // dependent on the input vectors. -psVector *psPolynomial3DEvalVector(const psPolynomial3D *poly, - const psVector *x, - const psVector *y, - const psVector *z) - +psVector *psPolynomial3DEvalVector( + const psPolynomial3D *poly, + const psVector *x, + const psVector *y, + const psVector *z) { PS_ASSERT_POLY_NON_NULL(poly, NULL); @@ -1036,9 +989,10 @@ } -psF64 psPolynomial4DEval(const psPolynomial4D* poly, - psF64 x, - psF64 y, - psF64 z, - psF64 t) +psF64 psPolynomial4DEval( + const psPolynomial4D* poly, + psF64 x, + psF64 y, + psF64 z, + psF64 t) { PS_ASSERT_POLY_NON_NULL(poly, NAN); @@ -1056,9 +1010,10 @@ } -psVector *psPolynomial4DEvalVector(const psPolynomial4D *poly, - const psVector *x, - const psVector *y, - const psVector *z, - const psVector *t) +psVector *psPolynomial4DEvalVector( + const psPolynomial4D *poly, + const psVector *x, + const psVector *y, + const psVector *z, + const psVector *t) { PS_ASSERT_POLY_NON_NULL(poly, NULL);