1/*-*-C-*-*/ 2 3#ifdef MEAN_FUNCTION 4static int MEAN_FUNCTION (VOID_STAR xp, unsigned int inc, unsigned int num, VOID_STAR yp) 5{ 6 GENERIC_TYPE *x, *xmax; 7 double c, m, err; 8 9 x = (GENERIC_TYPE *) xp; 10 xmax = x + num; 11 12 num = num/inc; 13 if (num == 0) 14 return 0; 15 16 c = (double) *x; 17 if (num == 1) 18 { 19 *(MEAN_RESULT_TYPE *)yp = (MEAN_RESULT_TYPE) c; 20 return 0; 21 } 22 23 err = 0.0; 24 m = c; 25 26 while (x < xmax) 27 { 28 double v = *x; 29 double dm = (v - c)/num; 30 double m1 = m + dm; 31 err += dm - (m1 - m); 32 m = m1; 33 x += inc; 34 } 35 *(MEAN_RESULT_TYPE *)yp = (MEAN_RESULT_TYPE) (m + err); 36 return 0; 37} 38# undef MEAN_FUNCTION 39# undef MEAN_RESULT_TYPE 40#endif 41 42#ifdef STDDEV_FUNCTION 43static int 44STDDEV_FUNCTION (VOID_STAR xp, unsigned int inc, unsigned int num, VOID_STAR s) 45{ 46 unsigned int i, n; 47 double mean_i, variance_i; 48 GENERIC_TYPE *x = (GENERIC_TYPE *) xp; 49 double err; 50 51 err = mean_i = variance_i = 0.0; 52 n = 1; 53 for (i = 0; i < num; i += inc) 54 { 55 double diff, x_i; 56 double dv, v1; 57 58 x_i = x[i]; 59 diff = x_i - mean_i; 60 mean_i += diff / n; 61 dv = diff * (x_i - mean_i); 62 v1 = variance_i + dv; 63 err += dv - (v1 - variance_i); 64 variance_i = v1; 65 n++; 66 } 67 variance_i += err; 68 n--; 69 if (n > 1) 70 variance_i = sqrt (variance_i / (n - 1)); 71 else 72 variance_i = 0.0; 73 74 *(STDDEV_RESULT_TYPE *)s = (STDDEV_RESULT_TYPE) variance_i; 75 76 return 0; 77} 78#undef STDDEV_RESULT_TYPE 79#undef STDDEV_FUNCTION 80#endif 81 82/* 83 * The following code is public domain. 84 * Algorithm by Torben Mogensen, implementation by N. Devillard. 85 * This code in public domain. 86 */ 87#ifdef NC_MEDIAN_FUNCTION 88static int NC_MEDIAN_FUNCTION (VOID_STAR ap, unsigned int inc, unsigned int num, VOID_STAR y) 89{ 90 unsigned int i, n, m; 91 GENERIC_TYPE *a; 92 GENERIC_TYPE min, max, guess, maxltguess, mingtguess; 93 unsigned int less, greater, equal; 94 95 n = num/inc; 96 if (n == 0) 97 { 98 SLang_set_error (SL_INVALID_PARM); 99 return -1; 100 } 101 a = (GENERIC_TYPE *)ap; 102 m = (n + 1)/2; 103 104 min = max = a[0]; 105 for (i = 0; i < num; i += inc) 106 { 107 GENERIC_TYPE ai = a[i]; 108 if (ai < min) min = ai; 109 if (ai > max) max = ai; 110 } 111 112 while (1) 113 { 114 less = 0; greater = 0; equal = 0; 115 /* We want guess=(min+max)/2. But min+max could overflow. So use 116 * (min+(max-min)+min)/2 = min + (max-min)/2 117 */ 118 guess = min + (max-min)/2; 119 maxltguess = min; 120 mingtguess = max; 121 122 for (i=0; i < num; i += inc) 123 { 124 GENERIC_TYPE ai = a[i]; 125 if (ai<guess) 126 { 127 less++; 128 if (ai>maxltguess) maxltguess = ai; 129 continue; 130 } 131 if (ai > guess) 132 { 133 greater++; 134 if (ai < mingtguess) mingtguess = ai; 135 continue; 136 } 137 equal++; 138 } 139 if ((less <= m) && (greater <= m)) 140 break; 141 142 if (less > greater) 143 max = maxltguess; 144 else 145 min = mingtguess; 146 } 147 if (less >= m) 148 guess = maxltguess; 149 else if (less + equal < m) 150 guess = mingtguess; 151 152 *(GENERIC_TYPE *)y = guess; 153 return 0; 154} 155#undef NC_MEDIAN_FUNCTION 156#endif 157 158#ifdef MEDIAN_FUNCTION 159static int 160MEDIAN_FUNCTION (VOID_STAR ap, unsigned int inc, unsigned int num, VOID_STAR y) 161{ 162 GENERIC_TYPE *aa = (GENERIC_TYPE *) ap; 163 GENERIC_TYPE *a; 164 unsigned int i, j, k, l, m, n; 165 166 n = num/inc; 167 if (n < 3) 168 { 169 if (n == 0) 170 { 171 SLang_set_error (SL_INVALID_PARM); 172 return -1; 173 } 174 if ((n == 1) 175 || (*aa < *(aa + inc))) 176 { 177 *(GENERIC_TYPE *)y = *aa; 178 return 0; 179 } 180 *(GENERIC_TYPE *)y = *(aa + inc); 181 return 0; 182 } 183 184 if (NULL == (a = (GENERIC_TYPE *) SLmalloc (sizeof (GENERIC_TYPE)*n))) 185 return -1; 186 187 for (i = 0; i < n; i++) 188 { 189 a[i] = *aa; 190 aa += inc; 191 } 192 193 if (n & 1) 194 k = n/2; 195 else 196 k = n/2 - 1; 197 198 l = 0; 199 m = n - 1; 200 201 while (l < m) 202 { 203 GENERIC_TYPE x = a[k]; 204 i = l; 205 j = m ; 206 do 207 { 208 while (a[i] < x) i++; 209 while (x < a[j]) j--; 210 if (i <= j) 211 { 212 GENERIC_TYPE tmp = a[i]; 213 a[i] = a[j]; 214 a[j] = tmp; 215 i++; 216 j--; 217 } 218 } 219 while (i<=j); 220 221 if (j < k) l=i; 222 if (k < i) m=j; 223 } 224 *(GENERIC_TYPE *) y = a[k]; 225 SLfree ((char *) a); 226 return 0; 227} 228#undef MEDIAN_FUNCTION 229#endif 230#undef GENERIC_TYPE 231