/* 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;
}