GRASS Programmer's Manual
6.4.2(2012)
|
00001 #include <stdio.h> 00002 #include <stdlib.h> 00003 #include <math.h> 00004 #include "local_proto.h" 00005 00006 00007 double *dmax(double *x, int n) 00008 { 00009 static double y[2]; 00010 double *xcopy, sqrt2, sqrtn, mean = 0.0, sdx = 0.0, fx; 00011 double dp, dp_max, dm, dm_max; 00012 int i; 00013 00014 if ((xcopy = (double *)malloc(n * sizeof(double))) == NULL) { 00015 fprintf(stderr, "Memory error in dmax\n"); 00016 exit(EXIT_FAILURE); 00017 } 00018 00019 sqrt2 = sqrt((double)2.0); 00020 sqrtn = sqrt((double)n); 00021 00022 for (i = 0; i < n; ++i) { 00023 xcopy[i] = x[i]; 00024 mean += x[i]; 00025 sdx += x[i] * x[i]; 00026 } 00027 sdx = sqrt((n * sdx - mean * mean) / (n * (n - 1.0))); 00028 mean /= n; 00029 00030 qsort(xcopy, n, sizeof(double), dcmp); 00031 00032 for (i = 0; i < n; ++i) { 00033 xcopy[i] = (xcopy[i] - mean) / sdx; 00034 fx = 0.5 + normp(xcopy[i] / sqrt2) / 2.0; 00035 if (fx <= 1e-5) 00036 fx = 1e-5; 00037 00038 if (fx >= 0.99999) 00039 fx = 0.99999; 00040 00041 dp = (double)(i + 1) / (double)n - fx; 00042 dm = fx - i / (double)n; 00043 if (i == 0 || dp > dp_max) 00044 dp_max = dp; 00045 00046 if (i == 0 || dm > dm_max) 00047 dm_max = dm; 00048 } 00049 00050 y[0] = dp_max; 00051 y[1] = dm_max; 00052 00053 free(xcopy); 00054 00055 return y; 00056 }