blob: ed3ad8a447eac9bcada5caa8691a35e68897a4cf (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
|
/* 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;
}
|