/* magstep.c 2.9.0 92/07/06 -- simple magstep calculations -- avoids math.h */ #include "stdc.h" #include "magstep.h" /* Copyright (c) 1991 Damian Cugley. */ /* * A faster routine (bit shifts, square if odd etc.) would * be a waste of time & space as this'll never be used with i > 10: */ double magstep(i) int i; { double f = 1.0; if (i < 0) return 1/magstep(-i); while (i-- > 0) f *= MAGSTEPONE; return f; } /* * Use e(x) for the difference between x and f. Want to minimize this. * assume a < f < Ma , * e(a) < e(Ma) <=> f - a < Ma - f <=> 2f < (M + 1)a */ double stepmag(f) double f; { double twof = 2*f, result; if (f <= 0.0) return 0.0; if (f < 1.0) return 1/stepmag(1/f); /* special cases for magstephalf (here M = m(1/2) ): */ if (twof < (1.0 + MAGSTEPHALF)) return 1.0; if (f < MAGSTEPHALF || twof < (1.0 + MAGSTEPHALF)*MAGSTEPHALF) return MAGSTEPHALF; /* now M = m(1) */ for (result = MAGSTEPONE; result < f; result *= MAGSTEPONE) if (twof < (1.0 + MAGSTEPONE)*result) return result; return result; }