1 //========================================================================
2 //
3 // FixedPoint.cc
4 //
5 // Fixed point type, with C++ operators.
6 //
7 // Copyright 2004 Glyph & Cog, LLC
8 //
9 //========================================================================
10
11 #include <aconf.h>
12
13 #if USE_FIXEDPOINT
14
15 #ifdef USE_GCC_PRAGMAS
16 #pragma implementation
17 #endif
18
19 #include "FixedPoint.h"
20
sqrt(FixedPoint x)21 FixedPoint FixedPoint::sqrt(FixedPoint x) {
22 FixedPoint y0, y1, z;
23
24 if (x.val <= 0) {
25 y1.val = 0;
26 } else {
27 y1.val = x.val >> 1;
28 do {
29 y0.val = y1.val;
30 z = x / y0;
31 y1.val = (y0.val + z.val) >> 1;
32 } while (::abs(y0.val - y1.val) > 1);
33 }
34 return y1;
35 }
36
37 //~ this is not very accurate
pow(FixedPoint x,FixedPoint y)38 FixedPoint FixedPoint::pow(FixedPoint x, FixedPoint y) {
39 FixedPoint t, t2, lnx0, lnx, z0, z;
40 int d, i;
41
42 if (y.val <= 0) {
43 z.val = 0;
44 } else {
45 // y * ln(x)
46 t = (x - 1) / (x + 1);
47 t2 = t * t;
48 d = 1;
49 lnx = 0;
50 do {
51 lnx0 = lnx;
52 lnx += t / d;
53 t *= t2;
54 d += 2;
55 } while (::abs(lnx.val - lnx0.val) > 2);
56 lnx.val <<= 1;
57 t = y * lnx;
58 // exp(y * ln(x))
59 t2 = t;
60 d = 1;
61 i = 1;
62 z = 1;
63 do {
64 z0 = z;
65 z += t2 / d;
66 t2 *= t;
67 ++i;
68 d *= i;
69 } while (::abs(z.val - z0.val) > 2 && d < (1 << fixptShift));
70 }
71 return z;
72 }
73
mul(int x,int y)74 int FixedPoint::mul(int x, int y) {
75 #if 1 //~tmp
76 return ((FixPtInt64)x * y) >> fixptShift;
77 #else
78 int ah0, ah, bh, al, bl;
79 ah0 = x & fixptMaskH;
80 ah = x >> fixptShift;
81 al = x - ah0;
82 bh = y >> fixptShift;
83 bl = y - (bh << fixptShift);
84 return ah0 * bh + ah * bl + al * bh + ((al * bl) >> fixptShift);
85 #endif
86 }
87
div(int x,int y)88 int FixedPoint::div(int x, int y) {
89 #if 1 //~tmp
90 return ((FixPtInt64)x << fixptShift) / y;
91 #else
92 #endif
93 }
94
95 #endif // USE_FIXEDPOINT
96