/* math.c -- runtime support of arithmetic and math builtins */ #include #include "rts.h" /* mpd_rtor returns the exponent x**y, where x and y are Real. * Error if (x=0 and y<=0) or (x<0 and y not int). */ Real mpd_rtor (locn, x, y) char *locn; Real x, y; { if ((x == 0) && (y <= 0)) mpd_loc_abort (locn, "0.0**y with y<=0"); else if ((x < 0) && (ceil (y) != y)) mpd_loc_abort (locn, "x**y with x<0 and y not an integer"); return (Real) pow (x, y); } /* returns the exponent i**j, where i and j are integers. j must * be positive. Error if (x=0 and y<=0). */ Real mpd_rtoi (locn, x, j) char *locn; Real x; int j; { Real term = x; Real result = 1; int bits; if ((x == 0) && (j <= 0)) mpd_loc_abort (locn, "0.0**i with i<=0"); for (bits = (j < 0 ? -j : j); bits; bits >>= 1, term *= term) if (bits & 1) result *= term; return j >= 0 ? result : 1.0 / result; } /* returns the exponent i**j, where i and j are integers. * Error if (x=0 and y<=0) or if j is negative. */ int mpd_itoi (locn, i, j) char *locn; int i, j; { int term = i; int result = 1; if (j < 0) mpd_loc_abort (locn, "i**j with j<0"); else if ((i == 0) && (j <= 0)) mpd_loc_abort (locn, "0**i with i<=0"); for (; j; j>>=1, term *= term) if (j & 1) result *= term; return result; } /* * return the nearest integer to x (as a Real), and if it is exactly * between two Reals then return the even ``neighbor.'' */ Real mpd_round (locn, x) char *locn; Real x; { Real xx, y; xx = x + 0.5; y = floor (xx); /* if halfway, and y ``odd'', then choose the even one */ if ((xx == y) && (mpd_rmod (locn, y, 2.0) == 1.0)) y -= 1.0; return y; } /* * x mod y == x - y * floor (x / y) * (from Knuth vol1 p38). * thus, * a) if y>0, 0 <= x mod y < y * b) if y<0, 0 >= x mod y > y * c) if y == 0 then error * d) note that x mod y is always of the same sign as y */ Real mpd_rmod (locn, x, y) char *locn; Real x, y; { if (y == 0.0) mpd_loc_abort (locn, "x mod y with y=0.0"); return x - y * ((Real) floor (x / y)); } /* * see the documentation for mpd_rmod for mpd_imod. */ int mpd_imod (locn, x, y) char *locn; int x, y; { if (y == 0) mpd_loc_abort (locn, "i mod j with j=0"); if (x >= 0) { if (y > 0) return x % y; else { int res; res = x % (-y); if (res == 0) return 0; else return res + y; } } else { if (y > 0) { int res; res = (-x) % y; if (res == 0) return 0; else return -res + y; } else return - ((-x) % (-y)); } } /* * mpd_imax (nargs, v, ...) -- return maximum of n integer values. */ int mpd_imax (int n, ...) { va_list ap; int r, v; va_start (ap, n); if (n <= 0) mpd_malf ("no args to mpd_imax"); r = va_arg (ap, int); /* pick off the first one */ n--; while (n--) { v = va_arg (ap, int); if (v > r) r = v; } va_end (ap); return r; } /* * mpd_imin (nargs, v, ...) -- return minimum of n integer values. */ int mpd_imin (int n, ...) { va_list ap; int r, v; va_start (ap, n); if (n <= 0) mpd_malf ("no args to mpd_imin"); r = va_arg (ap, int); /* pick off the first one */ n--; while (n--) { v = va_arg (ap, int); if (v < r) r = v; } va_end (ap); return r; } /* * mpd_rmax (nargs, v, ...) -- return maximum of n real values. */ Real mpd_rmax (int n, ...) { va_list ap; Real r, v; va_start (ap, n); if (n <= 0) mpd_malf ("no args to mpd_rmax"); r = va_arg (ap, Real); /* pick off the first one */ n--; while (n--) { v = va_arg (ap, Real); if (v > r) r = v; } va_end (ap); return r; } /* * mpd_rmin (nargs, v, ...) -- return minimum of n real values. */ Real mpd_rmin (int n, ...) { va_list ap; Real r, v; va_start (ap, n); if (n <= 0) mpd_malf ("no args to mpd_rmin"); r = va_arg (ap, Real); /* pick off the first one */ n--; while (n--) { v = va_arg (ap, Real); if (v < r) r = v; } va_end (ap); return r; } /* * Random number generator based on Knuth. * See Seminumerical Algorithms, 2nd edn, section 3.6. * Uses the subtractive series X[n] = (X[n-55] - X[n-24]) mod M, n >= 55. */ #define RANDMAX 1000000000 #define RAND_L 24 #define RAND_K 55 #define RAND_D 21 static int X[RAND_K]; /* history array */ static int cur = 0; /* Current index in array. */ static int curL = RAND_K - RAND_L; /* always (cur - RAND_L) mod RAND_K */ static Mutex random_mutex; /* protection for above */ /* * mpd_init_random () -- initialize random number generation. */ void mpd_init_random () { INIT_LOCK (random_mutex, "random_mutex"); mpd_seed (0.0); } /* * mpd_seed (x) -- seed the random number generator based on the value of x. * * If x is zero, an irreproducible value is used. */ void mpd_seed (x) Real x; { /* these are also protected by random mutex */ static int num_seed_call; static int seed_pid; static int seed_time; int i, j, k, seed, index; time_t time (); LOCK (random_mutex, "mpd_seed"); num_seed_call++; if (num_seed_call == 1) { seed_pid = getpid (); seed_time = time ((time_t *) 0); } if (x == 0.0) /* irreproducible */ seed = (997 * seed_pid + seed_time + 23347 * num_seed_call); else if (x > -1 && x < 1) seed = 1 / x; else seed = x; if (seed < 0) seed = -seed; seed = seed % RANDMAX; /* * According to Knuth, "This computes a Fibonacci-like sequence; * multiplication of indices by 21 [the constant RAND_D] spreads the * values around so as to alleviate initial nonrandomness problems..." */ X[0] = j = seed; k = 1; for (i = 1; i < RAND_K; i++) { index = (RAND_D * i) % RAND_K; X[index] = k; k = j - k; if (k < 0) k += RANDMAX; j = X[index]; } cur = 0; curL = RAND_K - RAND_L; UNLOCK (random_mutex, "mpd_seed"); /* * `warm up' the generator */ for (i = 0; i < 5 * RAND_K; i += 1) mpd_random (0.0, 1.0); } /* mpd_random (x, y) returns a random number r such that x <= r < y */ Real mpd_random (x, y) double x, y; { int iresult; LOCK (random_mutex, "mpd_random"); iresult = X[cur] - X[curL]; if (iresult < 0) iresult += RANDMAX; X[cur] = iresult; cur += 1; if (cur == RAND_K) cur = 0; curL += 1; if (curL == RAND_K) curL = 0; UNLOCK (random_mutex, "mpd_random"); return x + (y - x) * iresult * (1.0 / RANDMAX); }