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