1 |
/* |
2 |
* fpu/mathlib.cpp - Floating-point math support library |
3 |
* |
4 |
* Basilisk II (C) 1997-2008 Christian Bauer |
5 |
* |
6 |
* MC68881/68040 fpu emulation |
7 |
* |
8 |
* Original UAE FPU, copyright 1996 Herman ten Brugge |
9 |
* Rewrite for x86, copyright 1999-2001 Lauri Pesonen |
10 |
* New framework, copyright 2000-2001 Gwenole Beauchesne |
11 |
* Adapted for JIT compilation (c) Bernd Meyer, 2000-2001 |
12 |
* |
13 |
* This program is free software; you can redistribute it and/or modify |
14 |
* it under the terms of the GNU General Public License as published by |
15 |
* the Free Software Foundation; either version 2 of the License, or |
16 |
* (at your option) any later version. |
17 |
* |
18 |
* This program is distributed in the hope that it will be useful, |
19 |
* but WITHOUT ANY WARRANTY; without even the implied warranty of |
20 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
21 |
* GNU General Public License for more details. |
22 |
* |
23 |
* You should have received a copy of the GNU General Public License |
24 |
* along with this program; if not, write to the Free Software |
25 |
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
26 |
*/ |
27 |
|
28 |
/* NOTE: this file shall be included only from fpu/fpu_*.cpp */ |
29 |
#undef PRIVATE |
30 |
#define PRIVATE static |
31 |
|
32 |
#undef PUBLIC |
33 |
#define PUBLIC /**/ |
34 |
|
35 |
#undef FFPU |
36 |
#define FFPU /**/ |
37 |
|
38 |
#undef FPU |
39 |
#define FPU fpu. |
40 |
|
41 |
#if defined(FPU_IEEE) && defined(USE_X87_ASSEMBLY) |
42 |
|
43 |
PRIVATE fpu_extended fp_do_pow(fpu_extended x, fpu_extended y) |
44 |
{ |
45 |
fpu_extended value, exponent; |
46 |
uae_s64 p = (uae_s64)y; |
47 |
|
48 |
if (x == 0.0) { |
49 |
if (y > 0.0) |
50 |
return (y == (double) p && (p & 1) != 0 ? x : 0.0); |
51 |
else if (y < 0.0) |
52 |
return (y == (double) p && (-p & 1) != 0 ? 1.0 / x : 1.0 / fp_fabs (x)); |
53 |
} |
54 |
|
55 |
if (y == (double) p) { |
56 |
fpu_extended r = 1.0; |
57 |
if (p == 0) |
58 |
return 1.0; |
59 |
if (p < 0) { |
60 |
p = -p; |
61 |
x = 1.0 / x; |
62 |
} |
63 |
while (1) { |
64 |
if (p & 1) |
65 |
r *= x; |
66 |
p >>= 1; |
67 |
if (p == 0) |
68 |
return r; |
69 |
x *= x; |
70 |
} |
71 |
} |
72 |
|
73 |
__asm__ __volatile__("fyl2x" : "=t" (value) : "0" (x), "u" (1.0) : "st(1)"); |
74 |
__asm__ __volatile__("fmul %%st(1) # y * log2(x)\n\t" |
75 |
"fst %%st(1)\n\t" |
76 |
"frndint # int(y * log2(x))\n\t" |
77 |
"fxch\n\t" |
78 |
"fsub %%st(1) # fract(y * log2(x))\n\t" |
79 |
"f2xm1 # 2^(fract(y * log2(x))) - 1\n\t" |
80 |
: "=t" (value), "=u" (exponent) : "0" (y), "1" (value)); |
81 |
value += 1.0; |
82 |
__asm__ __volatile__("fscale" : "=t" (value) : "0" (value), "u" (exponent)); |
83 |
return value; |
84 |
} |
85 |
|
86 |
PRIVATE fpu_extended fp_do_log1p(fpu_extended x) |
87 |
{ |
88 |
// TODO: handle NaN and +inf/-inf |
89 |
fpu_extended value; |
90 |
// The fyl2xp1 can only be used for values in |
91 |
// -1 + sqrt(2) / 2 <= x <= 1 - sqrt(2) / 2 |
92 |
// 0.29 is a safe value. |
93 |
if (fp_fabs(x) <= 0.29) |
94 |
__asm__ __volatile__("fldln2; fxch; fyl2xp1" : "=t" (value) : "0" (x)); |
95 |
else |
96 |
__asm__ __volatile__("fldln2; fxch; fyl2x" : "=t" (value) : "0" (x + 1.0)); |
97 |
return value; |
98 |
} |
99 |
|
100 |
#endif |