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
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
//========================================================================
//
// FixedPoint.cc
//
// Fixed point type, with C++ operators.
//
// Copyright 2004 Glyph & Cog, LLC
//
//========================================================================
#include <config.h>
#if USE_FIXEDPOINT
#ifdef USE_GCC_PRAGMAS
#pragma implementation
#endif
#include "FixedPoint.h"
FixedPoint FixedPoint::sqrt(FixedPoint x) {
FixedPoint y0, y1, z;
if (x.val <= 0) {
y1.val = 0;
} else {
y1.val = x.val >> 1;
do {
y0.val = y1.val;
z = x / y0;
y1.val = (y0.val + z.val) >> 1;
} while (::abs(y0.val - y1.val) > 1);
}
return y1;
}
//~ this is not very accurate
FixedPoint FixedPoint::pow(FixedPoint x, FixedPoint y) {
FixedPoint t, t2, lnx0, lnx, z0, z;
int d, i;
if (y.val <= 0) {
z.val = 0;
} else {
// y * ln(x)
t = (x - 1) / (x + 1);
t2 = t * t;
d = 1;
lnx = 0;
do {
lnx0 = lnx;
lnx += t / d;
t *= t2;
d += 2;
} while (::abs(lnx.val - lnx0.val) > 2);
lnx.val <<= 1;
t = y * lnx;
// exp(y * ln(x))
t2 = t;
d = 1;
i = 1;
z = 1;
do {
z0 = z;
z += t2 / d;
t2 *= t;
++i;
d *= i;
} while (::abs(z.val - z0.val) > 2 && d < (1 << fixptShift));
}
return z;
}
int FixedPoint::mul(int x, int y) {
#if 1 //~tmp
return ((FixPtInt64)x * y) >> fixptShift;
#else
int ah0, ah, bh, al, bl;
ah0 = x & fixptMaskH;
ah = x >> fixptShift;
al = x - ah0;
bh = y >> fixptShift;
bl = y - (bh << fixptShift);
return ah0 * bh + ah * bl + al * bh + ((al * bl) >> fixptShift);
#endif
}
int FixedPoint::div(int x, int y) {
#if 1 //~tmp
return ((FixPtInt64)x << fixptShift) / y;
#else
#endif
}
#endif // USE_FIXEDPOINT
|
