1/*
2 *
3 *   Copyright (c) International Business Machines  Corp., 2002
4 *
5 *   This program is free software;  you can redistribute it and/or modify
6 *   it under the terms of the GNU General Public License as published by
7 *   the Free Software Foundation; either version 2 of the License, or
8 *   (at your option) any later version.
9 *
10 *   This program is distributed in the hope that it will be useful,
11 *   but WITHOUT ANY WARRANTY;  without even the implied warranty of
12 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
13 *   the GNU General Public License for more details.
14 *
15 *   You should have received a copy of the GNU General Public License
16 *   along with this program;  if not, write to the Free Software
17 *   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20/* Group Bull & IBM Corporation */
21/* 11/20/2002	Port to LTP	robbiew@us.ibm.com */
22/*                                               jacky.malcles@bull.net */
23/* IBM Corporation */
24/* 06/30/2001	Port to Linux	nsharoff@us.ibm.com */
25
26/*
27 * fptest01.c -- Floating point test.
28 *
29 * It is taken from a benchmark called "barsim".
30 *
31 * If the computation arrives at the expected values this routine
32 * prints a "passed" message and exits 0.  If an incorrect value is
33 * computed a "failed" message is printed and the routine exits 1.
34 */
35
36#include <stdio.h>
37#include <errno.h>
38#include <math.h>
39#include <stdlib.h>
40#include <unistd.h>
41#include <sys/types.h>
42#include <sys/stat.h>
43#include <fcntl.h>
44
45#define MAGIC1	1632.796126
46#define DIFF1	0.001
47#define MAGIC2	0.777807
48#define DIFF2	0.001
49#define EVENTMX	256
50#define BIG 1.e50
51#define FALSE 0
52#define TRUE  1
53#define TRYCRIT   1
54#define ENTERCRIT 2
55#define LEAVECRIT 3
56#define ATBARRIER 4
57#define ENTERWORK 5
58#define LEAVEWORK 6
59#define NULLEVENT 999
60
61/** LTP Port **/
62#include "test.h"
63
64char *TCID = "fptest01";	/* Test program identifier.    */
65int TST_TOTAL = 1;		/* Total number of test cases. */
66/**************/
67
68struct event {
69	int proc;
70	int type;
71	double time;
72};
73
74static int init(void);
75static int doevent(struct event *);
76static int term(void);
77static int addevent(int, int, double);
78
79static void gaussinit(double, double);
80static double gauss(void);
81
82struct event eventtab[EVENTMX];
83struct event rtrevent;
84int waiting[EVENTMX];		/* array of waiting processors */
85int nwaiting;			/* number of waiting processors */
86double global_time;		/* global clock */
87double lsttime;			/* time used for editing */
88double dtc, dts, alpha;		/* timing parameters */
89int nproc;			/* number of processors */
90int barcnt;			/* number of processors ATBARRIER */
91int ncycle;			/* number of cycles completed */
92int ncycmax;			/* number of cycles to run */
93int critfree;			/* TRUE if critical section not occupied */
94int gcount;			/* # calls to gauss */
95
96static struct event *nextevent(void);
97
98int main(int argc, char **argv)
99{
100	struct event *ev;
101
102	nproc = 128;
103	ncycmax = 10;
104	dtc = 0.01;
105	dts = 0.0;
106	alpha = 0.1;
107
108	init();
109
110	while ((ev = nextevent()) != NULL) {
111		doevent(ev);
112	}
113
114	term();
115	tst_resm(TPASS, "PASS");
116	tst_exit();
117}
118
119/*
120	initialize all processes to "entering work section"
121*/
122static int init(void)
123{
124	int p;
125	double dtw, dtwsig;
126
127	ncycle = 0;
128	global_time = 0;
129	lsttime = 0;
130	barcnt = 0;
131	nwaiting = 0;
132	critfree = TRUE;
133
134	dtw = 1. / nproc;	/* mean process work time */
135	dtwsig = dtw * alpha;	/* std deviation of work time */
136	gaussinit(dtw, dtwsig);
137
138	for (p = 1; p <= nproc; p++) {
139		eventtab[p].type = NULLEVENT;
140	}
141
142	for (p = 1; p <= nproc; p++) {
143		addevent(ENTERWORK, p, global_time);
144	}
145
146	return (0);
147}
148
149/*
150	print edit quantities
151*/
152static int term(void)
153{
154	double avgspd;
155	double t_total = 0.0;
156	double v;
157	int i;
158
159	for (i = 0; i < nproc; i++)
160		t_total += eventtab[i].time;
161
162	avgspd = ncycle / global_time;
163
164	v = t_total - MAGIC1;
165	if (v < 0.0)
166		v *= -1.0;
167
168	if (v > DIFF1) {
169		tst_resm(TFAIL, "FAIL");
170		v = t_total - MAGIC1;
171		tst_resm(TINFO, "t_total = %.15f\n", t_total);
172		tst_resm(TINFO, "expected  %.15f\n", MAGIC1);
173		tst_resm(TINFO, "diff = %.15f\n", v);
174		tst_exit();
175	}
176
177	v = avgspd - MAGIC2;
178	if (v < 0.0)
179		v *= -1.0;
180
181	if (v > DIFF2) {
182		tst_resm(TFAIL, "FAIL");
183		v = avgspd - MAGIC2;
184		tst_resm(TINFO, "avgspd  = %.15f\n", avgspd);
185		tst_resm(TINFO, "expected  %.15f\n", MAGIC2);
186		tst_resm(TINFO, "diff = %.15f\n", v);
187		tst_exit();
188	}
189	return (0);
190}
191
192/*
193	add an event to the event queue
194*/
195static int addevent(int type, int proc, double t)
196{
197	int i;
198	int ok = FALSE;
199
200	for (i = 1; i <= nproc; i++) {
201		if (eventtab[i].type == NULLEVENT) {
202			eventtab[i].type = type;
203			eventtab[i].proc = proc;
204			eventtab[i].time = t;
205			ok = TRUE;
206			break;
207		}
208	}
209	if (ok)
210		return (0);
211	else
212		tst_brkm(TBROK, NULL, "No room for event");
213	return (0);
214}
215
216/*
217	get earliest event in event queue
218*/
219static struct event *nextevent(void)
220{
221	double mintime = BIG;
222	int imin = 0;
223	int i;
224
225	for (i = 1; i <= nproc; i++) {
226		if (eventtab[i].type != NULLEVENT && eventtab[i].time < mintime) {
227			imin = i;
228			mintime = eventtab[i].time;
229		}
230	}
231
232	if (imin) {
233		rtrevent.type = eventtab[imin].type;
234		rtrevent.proc = eventtab[imin].proc;
235		rtrevent.time = eventtab[imin].time;
236		eventtab[imin].type = NULLEVENT;
237		return (&rtrevent);
238	} else
239		return (NULL);
240}
241
242/*
243	add a processor to the waiting queue
244*/
245static int addwaiting(int p)
246{
247	waiting[++nwaiting] = p;
248	return (0);
249}
250
251/*
252	remove the next processor from the waiting queue
253*/
254static int getwaiting(void)
255{
256	if (nwaiting)
257		return (waiting[nwaiting--]);
258	else
259		return (0);
260}
261
262static double dtcrit(void)
263{
264	return (dtc);
265}
266
267static double dtspinoff(void)
268{
269	return (dts);
270}
271
272static double dtwork(void)
273{
274	return (gauss());
275}
276
277/*
278	take the action prescribed by 'ev', update the clock, and
279	generate any subsequent events
280*/
281static int doevent(struct event *ev)
282{
283	double nxttime;
284	int i, p, proc;
285
286	global_time = ev->time;
287	proc = ev->proc;
288
289	switch (ev->type) {
290	case TRYCRIT:
291		if (critfree == TRUE)
292			addevent(ENTERCRIT, proc, global_time);
293		else
294			addwaiting(proc);
295		break;
296	case ENTERCRIT:
297		critfree = FALSE;
298		nxttime = global_time + dtcrit();
299		addevent(LEAVECRIT, proc, nxttime);
300		break;
301	case LEAVECRIT:
302		critfree = TRUE;
303		addevent(ATBARRIER, proc, global_time);
304		if ((p = getwaiting()) != 0) {
305			nxttime = global_time;
306			addevent(ENTERCRIT, p, nxttime);
307		}
308		break;
309	case ATBARRIER:
310		barcnt++;
311		if (barcnt == nproc) {
312			nxttime = global_time;
313			for (i = 1; i <= nproc; i++) {
314				nxttime += dtspinoff();
315				addevent(ENTERWORK, i, nxttime);
316			}
317			barcnt = 0;
318			ncycle++;
319		}
320		break;
321	case ENTERWORK:
322		nxttime = global_time + dtwork();
323		if (ncycle < ncycmax)
324			addevent(LEAVEWORK, proc, nxttime);
325		break;
326	case LEAVEWORK:
327		addevent(TRYCRIT, proc, global_time);
328		break;
329	default:
330		tst_brkm(TBROK, NULL, "Illegal event");
331		break;
332	}
333	return (0);
334}
335
336static int alternator = 1;
337static double mean;
338static double stdev;
339static double u1, u2;
340static double twopi;
341
342static void gaussinit(double m, double s)
343{
344	mean = m;
345	stdev = s;
346	twopi = 2. * acos((double)-1.0);
347	u1 = twopi / 400.0;
348	u2 = twopi / 500.0;
349}
350
351static double gauss(void)
352{
353	double x1, x2;
354
355	gcount++;
356
357	u1 += u2;
358	if (u1 > 0.99)
359		u1 = twopi / 500.0;
360	u2 += u1;
361	if (u2 > 0.99)
362		u2 = twopi / 400.0;
363
364	if (alternator == 1) {
365		alternator = -1;
366		x1 = sqrt(-2.0 * log(u1)) * cos(twopi * u2);
367		return (mean + stdev * x1);
368	} else {
369		alternator = 1;
370		x2 = sqrt(-2.0 * log(u1)) * sin(twopi * u2);
371		return (mean + stdev * x2);
372	}
373}
374