Listing 11
/* fmod function */ #include "xmath.h" double (fmod)(double x, double y) { /* compute fmod(x, y) */ const short errx = _Dtest(&x); const short erry = _Dtest(&y); if (errx == NAN || erry == NAN || errx == INF || erry == 0) { /* fmod undefined */ errno = EDOM; return (errx == NAN ? x : erry == NAN ? y : _Nan._D); } else if (errx == 0 || erry == INF) return (x); /* fmod(0,nonzero) or fmod(finite,INF) */ else { /* fmod(finite,finite) */ double t; short n, neg, ychar; if (y < 0.0) y = -y; if (x < 0.0) x = -x, neg = 1; else neg = 0; for (t = y, _Dunscale(&ychar, &t), n = 0; ; ) { /* subtract |y| until |x| y| */ short xchar; t = x; if (n < 0 || _Dunscale(&xchar, &t) == 0 || (n = xchar - ychar) 0) return (neg ? -x : x); for (; 0 <= n; -n) { /* try to subtract |y|*2^n */ t = y, _Dscale(&t, n); if (t <= x) { x -= t; break; } } } }