Index: trunk/psLib/src/math/psFunctions.c =================================================================== --- trunk/psLib/src/math/psFunctions.c (revision 4580) +++ trunk/psLib/src/math/psFunctions.c (revision 4581) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * -* @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-07-19 02:55:54 $ +* @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-07-20 01:21:13 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -50,8 +50,4 @@ static void polynomial3DFree(psPolynomial3D* poly); static void polynomial4DFree(psPolynomial4D* poly); -static void dPolynomial1DFree(psDPolynomial1D* poly); -static void dPolynomial2DFree(psDPolynomial2D* poly); -static void dPolynomial3DFree(psDPolynomial3D* poly); -static void dPolynomial4DFree(psDPolynomial4D* poly); static void spline1DFree(psSpline1D *tmpSpline); static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x); @@ -166,71 +162,4 @@ } -static void dPolynomial1DFree(psDPolynomial1D* poly) -{ - psFree(poly->coeff); - psFree(poly->coeffErr); - psFree(poly->mask); -} - -static void dPolynomial2DFree(psDPolynomial2D* poly) -{ - for (unsigned int x = 0; x < poly->nX; x++) { - psFree(poly->coeff[x]); - psFree(poly->coeffErr[x]); - psFree(poly->mask[x]); - } - psFree(poly->coeff); - psFree(poly->coeffErr); - psFree(poly->mask); -} - -static void dPolynomial3DFree(psDPolynomial3D* poly) -{ - unsigned int x = 0; - unsigned int y = 0; - - for (x = 0; x < poly->nX; x++) { - for (y = 0; y < poly->nY; y++) { - psFree(poly->coeff[x][y]); - psFree(poly->coeffErr[x][y]); - psFree(poly->mask[x][y]); - } - psFree(poly->coeff[x]); - psFree(poly->coeffErr[x]); - psFree(poly->mask[x]); - } - - psFree(poly->coeff); - psFree(poly->coeffErr); - psFree(poly->mask); -} - -static void dPolynomial4DFree(psDPolynomial4D* poly) -{ - unsigned int x = 0; - unsigned int y = 0; - unsigned int z = 0; - - for (x = 0; x < poly->nX; x++) { - for (y = 0; y < poly->nY; y++) { - for (z = 0; z < poly->nZ; z++) { - psFree(poly->coeff[x][y][z]); - psFree(poly->coeffErr[x][y][z]); - psFree(poly->mask[x][y][z]); - } - psFree(poly->coeff[x][y]); - psFree(poly->coeffErr[x][y]); - psFree(poly->mask[x][y]); - } - psFree(poly->coeff[x]); - psFree(poly->coeffErr[x]); - psFree(poly->mask[x]); - } - - psFree(poly->coeff); - psFree(poly->coeffErr); - psFree(poly->mask); -} - /***************************************************************************** createChebyshevPolys(n): this routine takes as input the required order n, @@ -283,14 +212,14 @@ { psS32 loop_x = 0; - psF32 polySum = 0.0; - psF32 xSum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, - "---- Calling ordPolynomial1DEval(%f)\n", x); + "---- Calling ordPolynomial1DEval(%lf)\n", x); psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, "Polynomial order is %d\n", poly->n); for (loop_x = 0; loop_x < poly->n; loop_x++) { psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, - "Polynomial coeff[%d] is %f\n", loop_x, poly->coeff[loop_x]); + "Polynomial coeff[%d] is %lf\n", loop_x, poly->coeff[loop_x]); } @@ -298,5 +227,5 @@ if (poly->mask[loop_x] == 0) { psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 10, - "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, poly->coeff[loop_x]); + "polysum+= sum*coeff [%lf+= (%lf * %lf)\n", polySum, xSum, poly->coeff[loop_x]); polySum += xSum * poly->coeff[loop_x]; } @@ -312,5 +241,5 @@ static psF64 chebPolynomial1DEval(psF64 x, const psPolynomial1D* poly) { - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); // XXX: Create a macro for this in psConstants.h if (poly->n < 1) { @@ -321,5 +250,5 @@ psS32 n = poly->n; psS32 i; - psF32 tmp = 0.0; + psF64 tmp = 0.0; // Special case where the Chebyshev poly is constant. @@ -343,26 +272,26 @@ // General case where the Chebyshev poly has 2 or more terms. - d = psVectorAlloc(n, PS_TYPE_F32); + d = psVectorAlloc(n, PS_TYPE_F64); if(poly->mask[n-1] == 0) { - d->data.F32[n-1] = poly->coeff[n-1]; + d->data.F64[n-1] = poly->coeff[n-1]; } else { - d->data.F32[n-1] = 0.0; - } - - d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]); + 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.F32[n-2] += poly->coeff[n-2]; + d->data.F64[n-2] += poly->coeff[n-2]; } for (i=n-3;i>=1;i--) { - d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) - - (d->data.F32[i+2]); + d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) - + (d->data.F64[i+2]); if(poly->mask[i] == 0) { - d->data.F32[i] += poly->coeff[i]; - } - } - - tmp = (x * d->data.F32[1]) - - (d->data.F32[2]); + 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]); @@ -400,7 +329,7 @@ psS32 loop_x = 0; psS32 loop_y = 0; - psF32 polySum = 0.0; - psF32 xSum = 1.0; - psF32 ySum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; for (loop_x = 0; loop_x < poly->nX; loop_x++) { @@ -420,6 +349,6 @@ static psF64 chebPolynomial2DEval(psF64 x, psF64 y, const psPolynomial2D* poly) { - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0); PS_ASSERT_POLY_NON_NULL(poly, NAN); @@ -427,5 +356,5 @@ psS32 loop_y = 0; psS32 i = 0; - psF32 polySum = 0.0; + psF64 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -460,8 +389,8 @@ psS32 loop_y = 0; psS32 loop_z = 0; - psF32 polySum = 0.0; - psF32 xSum = 1.0; - psF32 ySum = 1.0; - psF32 zSum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; + psF64 zSum = 1.0; for (loop_x = 0; loop_x < poly->nX; loop_x++) { @@ -485,12 +414,12 @@ static psF64 chebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psPolynomial3D* poly) { - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(z, -1.0, 1.0, 0.0); psS32 loop_x = 0; psS32 loop_y = 0; psS32 loop_z = 0; psS32 i = 0; - psF32 polySum = 0.0; + psF64 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -533,9 +462,9 @@ psS32 loop_z = 0; psS32 loop_t = 0; - psF32 polySum = 0.0; - psF32 xSum = 1.0; - psF32 ySum = 1.0; - psF32 zSum = 1.0; - psF32 tSum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; + psF64 zSum = 1.0; + psF64 tSum = 1.0; for (loop_x = 0; loop_x < poly->nX; loop_x++) { @@ -563,8 +492,8 @@ static psF64 chebPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psPolynomial4D* poly) { - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(t, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(z, -1.0, 1.0, 0.0); + PS_ASSERT_DOUBLE_WITHIN_RANGE(t, -1.0, 1.0, 0.0); psS32 loop_x = 0; psS32 loop_y = 0; @@ -572,5 +501,5 @@ psS32 loop_t = 0; psS32 i = 0; - psF32 polySum = 0.0; + psF64 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -612,281 +541,4 @@ return(polySum); } - -/***************************************************************************** - Polynomial coefficients will be accessed in [w][x][y][z] fashion. - *****************************************************************************/ -static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* poly) -{ - psS32 loop_x = 0; - psF64 polySum = 0.0; - psF64 xSum = 1.0; - - for (loop_x = 0; loop_x < poly->n; loop_x++) { - if (poly->mask[loop_x] == 0) { - polySum += xSum * poly->coeff[loop_x]; - } - xSum *= x; - } - - return(polySum); -} - -// XXX: You can do this without having to psAlloc() vector d. -// XXX: How does the mask vector effect Crenshaw's formula? -static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - psVector *d; - psS32 n; - psS32 i; - psF64 tmp; - - n = poly->n; - d = psVectorAlloc(n, PS_TYPE_F64); - if(poly->mask[n-1] == 0) { - d->data.F64[n-1] = poly->coeff[n-1]; - } 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); - return(tmp); -} - -static psF64 dOrdPolynomial2DEval(psF64 x, - psF64 y, - const psDPolynomial2D* poly) -{ - psS32 loop_x = 0; - psS32 loop_y = 0; - psF64 polySum = 0.0; - psF64 xSum = 1.0; - psF64 ySum = 1.0; - - for (loop_x = 0; loop_x < poly->nX; loop_x++) { - ySum = xSum; - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - if (poly->mask[loop_x][loop_y] == 0) { - polySum += ySum * poly->coeff[loop_x][loop_y]; - } - ySum *= y; - } - xSum *= x; - } - - return(polySum); -} - -static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 i = 0; - psF64 polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - psS32 maxChebyPoly = 0; - - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = poly->nX; - if (poly->nY > maxChebyPoly) { - maxChebyPoly = poly->nY; - } - chebPolys = createChebyshevPolys(maxChebyPoly); - - for (loop_x = 0; loop_x < poly->nX; loop_x++) { - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - if (poly->mask[loop_x][loop_y] == 0) { - polySum += poly->coeff[loop_x][loop_y] * - psPolynomial1DEval(chebPolys[loop_x], x) * - psPolynomial1DEval(chebPolys[loop_y], y); - } - } - } - - for (i=0;inX; loop_x++) { - ySum = xSum; - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - zSum = ySum; - for (loop_z = 0; loop_z < poly->nZ; loop_z++) { - if (poly->mask[loop_x][loop_y][loop_z] == 0) { - polySum += zSum * poly->coeff[loop_x][loop_y][loop_z]; - } - zSum *= z; - } - ySum *= y; - } - xSum *= x; - } - - return(polySum); -} - -static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - psS32 i = 0; - psF64 polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - psS32 maxChebyPoly = 0; - - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = poly->nX; - if (poly->nY > maxChebyPoly) { - maxChebyPoly = poly->nY; - } - if (poly->nZ > maxChebyPoly) { - maxChebyPoly = poly->nZ; - } - chebPolys = createChebyshevPolys(maxChebyPoly); - - for (loop_x = 0; loop_x < poly->nX; loop_x++) { - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - for (loop_z = 0; loop_z < poly->nZ; loop_z++) { - if (poly->mask[loop_x][loop_y][loop_z] == 0) { - polySum += poly->coeff[loop_x][loop_y][loop_z] * - psPolynomial1DEval(chebPolys[loop_x], x) * - psPolynomial1DEval(chebPolys[loop_y], y) * - psPolynomial1DEval(chebPolys[loop_z], z); - } - } - } - } - - for (i=0;inX; loop_x++) { - ySum = xSum; - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - zSum = ySum; - for (loop_z = 0; loop_z < poly->nZ; loop_z++) { - tSum = zSum; - for (loop_t = 0; loop_t < poly->nT; loop_t++) { - if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) { - polySum += tSum * poly->coeff[loop_x][loop_y][loop_z][loop_t]; - } - tSum *= t; - } - zSum *= z; - } - ySum *= y; - } - xSum *= x; - } - - return(polySum); -} - -static psF64 dChebPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psDPolynomial4D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(t, -1.0, 1.0, 0.0); - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - psS32 loop_t = 0; - psS32 i = 0; - psF64 polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - psS32 maxChebyPoly = 0; - - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = poly->nX; - if (poly->nY > maxChebyPoly) { - maxChebyPoly = poly->nY; - } - if (poly->nZ > maxChebyPoly) { - maxChebyPoly = poly->nZ; - } - if (poly->nT > maxChebyPoly) { - maxChebyPoly = poly->nT; - } - chebPolys = createChebyshevPolys(maxChebyPoly); - - for (loop_x = 0; loop_x < poly->nX; loop_x++) { - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - for (loop_z = 0; loop_z < poly->nZ; loop_z++) { - for (loop_t = 0; loop_t < poly->nT; loop_t++) { - if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) { - polySum += poly->coeff[loop_x][loop_y][loop_z][loop_t] * - psPolynomial1DEval(chebPolys[loop_x], x) * - psPolynomial1DEval(chebPolys[loop_y], y) * - psPolynomial1DEval(chebPolys[loop_z], z) * - psPolynomial1DEval(chebPolys[loop_t], t); - } - } - } - } - } - - for (i=0;itype = type; newPoly->n = n; - newPoly->coeff = (psF32 *)psAlloc(n * sizeof(psF32)); - newPoly->coeffErr = (psF32 *)psAlloc(n * sizeof(psF32)); + newPoly->coeff = psAlloc(n * sizeof(psF64)); + newPoly->coeffErr = psAlloc(n * sizeof(psF64)); newPoly->mask = (char *)psAlloc(n * sizeof(char)); for (i = 0; i < n; i++) { @@ -1111,10 +763,10 @@ newPoly->nY = nY; - newPoly->coeff = (psF32 **)psAlloc(nX * sizeof(psF32 *)); - newPoly->coeffErr = (psF32 **)psAlloc(nX * sizeof(psF32 *)); + newPoly->coeff = psAlloc(nX * sizeof(psF64 *)); + newPoly->coeffErr = psAlloc(nX * sizeof(psF64 *)); newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); - newPoly->coeffErr[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); + newPoly->coeff[x] = psAlloc(nY * sizeof(psF64)); + newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64)); newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char)); } @@ -1150,14 +802,14 @@ newPoly->nZ = nZ; - newPoly->coeff = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); - newPoly->coeffErr = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); + newPoly->coeff = psAlloc(nX * sizeof(psF64 **)); + newPoly->coeffErr = psAlloc(nX * sizeof(psF64 **)); newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); - newPoly->coeffErr[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); + newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 *)); + newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 *)); newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); - newPoly->coeffErr[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); + newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64)); + newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64)); newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char)); } @@ -1199,18 +851,18 @@ newPoly->nT = nT; - newPoly->coeff = (psF32 ****)psAlloc(nX * sizeof(psF32 ***)); - newPoly->coeffErr = (psF32 ****)psAlloc(nX * sizeof(psF32 ***)); + newPoly->coeff = psAlloc(nX * sizeof(psF64 ***)); + newPoly->coeffErr = psAlloc(nX * sizeof(psF64 ***)); newPoly->mask = (char ****)psAlloc(nX * sizeof(char ***)); for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF32 ***)psAlloc(nY * sizeof(psF32 **)); - newPoly->coeffErr[x] = (psF32 ***)psAlloc(nY * sizeof(psF32 **)); + newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 **)); + newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 **)); newPoly->mask[x] = (char ***)psAlloc(nY * sizeof(char **)); for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = (psF32 **)psAlloc(nZ * sizeof(psF32 *)); - newPoly->coeffErr[x][y] = (psF32 **)psAlloc(nZ * sizeof(psF32 *)); + newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64 *)); + newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64 *)); newPoly->mask[x][y] = (char **)psAlloc(nZ * sizeof(char *)); for (z = 0; z < nZ; z++) { - newPoly->coeff[x][y][z] = (psF32 *)psAlloc(nT * sizeof(psF32)); - newPoly->coeffErr[x][y][z] = (psF32 *)psAlloc(nT * sizeof(psF32)); + newPoly->coeff[x][y][z] = psAlloc(nT * sizeof(psF64)); + newPoly->coeffErr[x][y][z] = psAlloc(nT * sizeof(psF64)); newPoly->mask[x][y][z] = (char *)psAlloc(nT * sizeof(char)); } @@ -1253,11 +905,11 @@ PS_ASSERT_POLY_NON_NULL(poly, NULL); PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); psVector *tmp; - tmp = psVectorAlloc(x->n, PS_TYPE_F32); + tmp = psVectorAlloc(x->n, PS_TYPE_F64); for (psS32 i=0;in;i++) { - tmp->data.F32[i] = psPolynomial1DEval(poly, x->data.F32[i]); + tmp->data.F64[i] = psPolynomial1DEval(poly, x->data.F64[i]); } @@ -1288,7 +940,7 @@ PS_ASSERT_POLY_NON_NULL(poly, NULL); PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); psVector *tmp; @@ -1301,9 +953,9 @@ // Create output vector to return - tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); // Evaluate the polynomial at the specified points for (psS32 i=0; idata.F32[i] = psPolynomial2DEval(poly,x->data.F32[i],y->data.F32[i]); + tmp->data.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]); } @@ -1336,9 +988,9 @@ PS_ASSERT_POLY_NON_NULL(poly, NULL); PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL); psVector *tmp; @@ -1354,12 +1006,12 @@ // Allocate output vector - tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); // Evaluate polynomial for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F32[i] = psPolynomial3DEval(poly, - x->data.F32[i], - y->data.F32[i], - z->data.F32[i]); + tmp->data.F64[i] = psPolynomial3DEval(poly, + x->data.F64[i], + y->data.F64[i], + z->data.F64[i]); } @@ -1389,378 +1041,4 @@ const psVector *z, const psVector *t) -{ - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(t, NULL); - PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, NULL); - - psVector *tmp; - psS32 vecLen=x->n; - - // Determine output vector size from min of input vectors - if (z->n < vecLen) { - vecLen = z->n; - } - if (y->n < vecLen) { - vecLen = y->n; - } - if (t->n < vecLen) { - vecLen = t->n; - } - - // Allocate output vector - tmp = psVectorAlloc(vecLen, PS_TYPE_F32); - - // Evaluate polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F32[i] = psPolynomial4DEval(poly, - x->data.F32[i], - y->data.F32[i], - z->data.F32[i], - t->data.F32[i]); - } - - // Return output vector - return(tmp); -} - - -psDPolynomial1D* psDPolynomial1DAlloc( int n, - psPolynomialType type) -{ - PS_ASSERT_INT_POSITIVE(n, NULL); - - unsigned int i = 0; - psDPolynomial1D* newPoly = NULL; - - newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D)); - psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial1DFree); - - newPoly->type = type; - newPoly->n = n; - newPoly->coeff = (psF64 *)psAlloc(n * sizeof(psF64)); - newPoly->coeffErr = (psF64 *)psAlloc(n * sizeof(psF64)); - newPoly->mask = (char *)psAlloc(n * sizeof(char)); - for (i = 0; i < n; i++) { - newPoly->coeff[i] = 0.0; - newPoly->coeffErr[i] = 0.0; - newPoly->mask[i] = 0; - } - - return(newPoly); -} - -psDPolynomial2D* psDPolynomial2DAlloc( int nX, int nY, - psPolynomialType type) -{ - PS_ASSERT_INT_POSITIVE(nX, NULL); - PS_ASSERT_INT_POSITIVE(nY, NULL); - - unsigned int x = 0; - unsigned int y = 0; - psDPolynomial2D* newPoly = NULL; - - newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D)); - psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial2DFree); - - newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - - newPoly->coeff = (psF64 **)psAlloc(nX * sizeof(psF64 *)); - newPoly->coeffErr = (psF64 **)psAlloc(nX * sizeof(psF64 *)); - newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); - newPoly->coeffErr[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); - newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char)); - } - for (x = 0; x < nX; x++) { - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = 0.0; - newPoly->coeffErr[x][y] = 0.0; - newPoly->mask[x][y] = 0; - } - } - - return(newPoly); -} - -psDPolynomial3D* psDPolynomial3DAlloc( int nX, int nY, int nZ, - psPolynomialType type) -{ - PS_ASSERT_INT_POSITIVE(nX, NULL); - PS_ASSERT_INT_POSITIVE(nY, NULL); - PS_ASSERT_INT_POSITIVE(nZ, NULL); - - unsigned int x = 0; - unsigned int y = 0; - unsigned int z = 0; - psDPolynomial3D* newPoly = NULL; - - newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D)); - psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial3DFree); - - newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - newPoly->nZ = nZ; - - newPoly->coeff = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); - newPoly->coeffErr = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); - newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); - newPoly->coeffErr[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); - newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); - newPoly->coeffErr[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); - newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char)); - } - } - for (x = 0; x < nX; x++) { - for (y = 0; y < nY; y++) { - for (z = 0; z < nZ; z++) { - newPoly->coeff[x][y][z] = 0.0; - newPoly->coeffErr[x][y][z] = 0.0; - newPoly->mask[x][y][z] = 0; - } - } - } - - return(newPoly); -} - -psDPolynomial4D* psDPolynomial4DAlloc( int nX, int nY, int nZ, int nT, - psPolynomialType type) -{ - PS_ASSERT_INT_POSITIVE(nX, NULL); - PS_ASSERT_INT_POSITIVE(nY, NULL); - PS_ASSERT_INT_POSITIVE(nZ, NULL); - PS_ASSERT_INT_POSITIVE(nT, NULL); - - unsigned int x = 0; - unsigned int y = 0; - unsigned int z = 0; - unsigned int t = 0; - psDPolynomial4D* newPoly = NULL; - - newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D)); - psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial4DFree); - - newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - newPoly->nZ = nZ; - newPoly->nT = nT; - - newPoly->coeff = (psF64 ****)psAlloc(nX * sizeof(psF64 ***)); - newPoly->coeffErr = (psF64 ****)psAlloc(nX * sizeof(psF64 ***)); - newPoly->mask = (char ****)psAlloc(nX * sizeof(char ***)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF64 ***)psAlloc(nY * sizeof(psF64 **)); - newPoly->coeffErr[x] = (psF64 ***)psAlloc(nY * sizeof(psF64 **)); - newPoly->mask[x] = (char ***)psAlloc(nY * sizeof(char **)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = (psF64 **)psAlloc(nZ * sizeof(psF64 *)); - newPoly->coeffErr[x][y] = (psF64 **)psAlloc(nZ * sizeof(psF64 *)); - newPoly->mask[x][y] = (char **)psAlloc(nZ * sizeof(char *)); - for (z = 0; z < nZ; z++) { - newPoly->coeff[x][y][z] = (psF64 *)psAlloc(nT * sizeof(psF64)); - newPoly->coeffErr[x][y][z] = (psF64 *)psAlloc(nT * sizeof(psF64)); - newPoly->mask[x][y][z] = (char *)psAlloc(nT * sizeof(char)); - } - } - } - for (x = 0; x < nX; x++) { - for (y = 0; y < nY; y++) { - for (z = 0; z < nZ; z++) { - for (t = 0; t < nT; t++) { - newPoly->coeff[x][y][z][t] = 0.0; - newPoly->coeffErr[x][y][z][t] = 0.0; - newPoly->mask[x][y][z][t] = 0; - } - } - } - } - - return(newPoly); -} - - -psF64 psDPolynomial1DEval(const psDPolynomial1D* poly, psF64 x) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial1DEval(x, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial1DEval(x, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psDPolynomial1DEvalVector(const psDPolynomial1D *poly, - const psVector *x) - -{ - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - - psVector *tmp; - - tmp = psVectorAlloc(x->n, PS_TYPE_F64); - for (psS32 i=0;in;i++) { - tmp->data.F64[i] = psDPolynomial1DEval(poly, - x->data.F64[i]); - } - - return(tmp); -} - - -psF64 psDPolynomial2DEval(const psDPolynomial2D* poly, - psF64 x, - psF64 y) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - if (poly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial2DEval(x, y, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial2DEval(x, y, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psDPolynomial2DEvalVector(const psDPolynomial2D *poly, - const psVector *x, - const psVector *y) -{ - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); - - psVector *tmp; - psS32 vecLen=x->n; - - // Determine the output vector length from minimum length of input vectors - if (y->n < vecLen) { - vecLen = y->n; - } - - // Allocate output vector - tmp = psVectorAlloc(vecLen, PS_TYPE_F64); - - // Evaluate the polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F64[i] = psDPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]); - } - - // Return output vector - return(tmp); -} - - -psF64 psDPolynomial3DEval(const psDPolynomial3D* poly, - psF64 x, - psF64 y, - psF64 z) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial3DEval(x, y, z, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial3DEval(x, y, z, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psDPolynomial3DEvalVector(const psDPolynomial3D *poly, - const psVector *x, - const psVector *y, - const psVector *z) - -{ - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL); - - psVector *tmp; - psS32 vecLen=x->n; - - // Determine the size of output vector from min of input vectors - if (y->n < vecLen) { - vecLen = y->n; - } - if (z->n < vecLen) { - vecLen = z->n; - } - - // Allocate output vector - tmp = psVectorAlloc(vecLen, PS_TYPE_F64); - - // Evaluate polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F64[i] = psDPolynomial3DEval(poly, - x->data.F64[i], - y->data.F64[i], - z->data.F64[i]); - } - - // Return output vector - return(tmp); -} - -psF64 psDPolynomial4DEval(const psDPolynomial4D* poly, - psF64 x, - psF64 y, - psF64 z, - psF64 t) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial4DEval(x,y,z,t, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial4DEval(x,y,z,t, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psDPolynomial4DEvalVector(const psDPolynomial4D *poly, - const psVector *x, - const psVector *y, - const psVector *z, - const psVector *t) { PS_ASSERT_POLY_NON_NULL(poly, NULL); @@ -1777,5 +1055,5 @@ psS32 vecLen=x->n; - // Determine the output vector size from min of input vectors + // Determine output vector size from min of input vectors if (z->n < vecLen) { vecLen = z->n; @@ -1791,11 +1069,11 @@ tmp = psVectorAlloc(vecLen, PS_TYPE_F64); - // Evaluate the polynomial + // Evaluate polynomial for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F64[i] = psDPolynomial4DEval(poly, - x->data.F64[i], - y->data.F64[i], - z->data.F64[i], - t->data.F64[i]); + tmp->data.F64[i] = psPolynomial4DEval(poly, + x->data.F64[i], + y->data.F64[i], + z->data.F64[i], + t->data.F64[i]); } @@ -1803,7 +1081,4 @@ return(tmp); } - - - //typedef struct {