Mercurial > forge

changeset 5797:d43bc9e56a0c octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
using internally long double for increased accuracy
author schloegl
date Tue, 30 Jun 2009 12:15:57 +0000
parents 5a8ad92a3e1c
children 4507bf3417b3
files extra/NaN/src/covm_mex.cpp extra/NaN/src/sumskipnan_mex.cpp
diffstat 2 files changed, 19 insertions(+), 18 deletions(-) [+]
line wrap: on
line diff
--- a/extra/NaN/src/covm_mex.cpp	Tue Jun 30 07:49:32 2009 +0000
+++ b/extra/NaN/src/covm_mex.cpp	Tue Jun 30 12:15:57 2009 +0000
@@ -43,6 +43,7 @@
 
 //#define NO_FLAG
 
+
 void mexFunction(int POutputCount,  mxArray* POutput[], int PInputCount, const mxArray *PInputs[]) 
 {
     	double 		*X0,*Y0=NULL,*X,*Y,*W=NULL;
@@ -200,8 +201,8 @@
 	    for (j=0; j<cY; j++) {
 		X = X0+i*rX;
 		Y = Y0+j*rY;
-		double cc=0.0;
-		double nn=0.0;
+		long double cc=0.0;
+		long double nn=0.0;
 		for (k=0; k<rX; k++) {
 			double z = X[k]*Y[k];
 			if (isnan(z)) {
@@ -222,7 +223,7 @@
 	    for (j=0; j<cY; j++) {
 		X = X0+i*rX;
 		Y = Y0+j*rY;
-		double cc=0.0;
+		long double cc=0.0;
 		size_t nn=0;
 		for (k=0; k<rX; k++) {
 			double z = X[k]*Y[k];
@@ -246,8 +247,8 @@
 	    for (j=i; j<cY; j++) {
 		X = X0+i*rX;
 		Y = Y0+j*rY;
-		double cc=0.0;
-		double nn=0.0;
+		long double cc=0.0;
+		long double nn=0.0;
 		for (k=0; k<rX; k++) {
 			double z = X[k]*Y[k];
 			if (isnan(z)) {
@@ -271,7 +272,7 @@
 	    for (j=i; j<cY; j++) {
 		X = X0+i*rX;
 		Y = Y0+j*rY;
-		double cc=0.0;
+		long double cc=0.0;
 		size_t nn=0;
 		for (k=0; k<rX; k++) {
 			double z = X[k]*Y[k];
--- a/extra/NaN/src/sumskipnan_mex.cpp	Tue Jun 30 07:49:32 2009 +0000
+++ b/extra/NaN/src/sumskipnan_mex.cpp	Tue Jun 30 12:15:57 2009 +0000
@@ -79,9 +79,9 @@
 
 	// check for proper number of input and output arguments
 	if ((PInputCount <= 0) || (PInputCount > 4))
-	        mexErrMsgTxt("SumSkipNan.MEX requires between 1 and 4 arguments.");
+	        mexErrMsgTxt("SUMSKIPNAN.MEX requires between 1 and 4 arguments.");
 	if (POutputCount > 4)
-	        mexErrMsgTxt("SumSkipNan.MEX has 1 to 3 output arguments.");
+	        mexErrMsgTxt("SUMSKIPNAN.MEX has 1 to 3 output arguments.");
 
 	// get 1st argument
 	if(mxIsDouble(PInputs[0]) && !mxIsComplex(PInputs[0]))
@@ -346,7 +346,7 @@
 #define stride 1 
 inline int __sumskipnan2__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN)
 {
-	register double sum=0; 
+	register long double sum=0; 
 	register size_t count=0; 
 	register char   flag=0; 
 	// LOOP  along dimension DIM
@@ -364,7 +364,7 @@
 			flag = 1; 
 #endif
 
-		data +=stride;	
+		data++;	// stride=1
 	}
 	while (data < end);
 	
@@ -378,8 +378,8 @@
 
 inline int __sumskipnan3__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN)
 {
-	register double sum=0; 
-	register double msq=0; 
+	register long double sum=0; 
+	register long double msq=0; 
 	register size_t count=0; 
 	register char   flag=0; 
 	// LOOP  along dimension DIM
@@ -412,8 +412,8 @@
 #define stride 1 
 inline int __sumskipnan2w__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W)
 {
-	register double sum=0; 
-	register double count=0; 
+	register long double sum=0; 
+	register long double count=0; 
 	register char   flag=0; 
 	// LOOP  along dimension DIM
 	
@@ -430,7 +430,7 @@
 			flag = 1; 
 #endif
 
-		data +=stride;	
+		data++;	// stride=1
 		W++;
 	}
 	while (data < end);
@@ -445,9 +445,9 @@
 
 inline int __sumskipnan3w__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W)
 {
-	register double sum=0; 
-	register double msq=0; 
-	register double count=0; 
+	register long double sum=0; 
+	register long double msq=0; 
+	register long double count=0; 
 	register char   flag=0; 
 	// LOOP  along dimension DIM