1 /* 1 /*
2 2
3 fp_log.c: floating-point math routines for t 3 fp_log.c: floating-point math routines for the Linux-m68k
4 floating point emulator. 4 floating point emulator.
5 5
6 Copyright (c) 1998-1999 David Huggins-Daines 6 Copyright (c) 1998-1999 David Huggins-Daines / Roman Zippel.
7 7
8 I hereby give permission, free of charge, to 8 I hereby give permission, free of charge, to copy, modify, and
9 redistribute this software, in source or bin 9 redistribute this software, in source or binary form, provided that
10 the above copyright notice and the following 10 the above copyright notice and the following disclaimer are included
11 in all such copies. 11 in all such copies.
12 12
13 THIS SOFTWARE IS PROVIDED "AS IS", WITH ABSO 13 THIS SOFTWARE IS PROVIDED "AS IS", WITH ABSOLUTELY NO WARRANTY, REAL
14 OR IMPLIED. 14 OR IMPLIED.
15 15
16 */ 16 */
17 17
18 #include "fp_arith.h" 18 #include "fp_arith.h"
19 #include "fp_emu.h" 19 #include "fp_emu.h"
20 #include "fp_log.h" 20 #include "fp_log.h"
21 21
22 static const struct fp_ext fp_one = { 22 static const struct fp_ext fp_one = {
23 .exp = 0x3fff, 23 .exp = 0x3fff,
24 }; 24 };
25 25
26 struct fp_ext *fp_fsqrt(struct fp_ext *dest, s 26 struct fp_ext *fp_fsqrt(struct fp_ext *dest, struct fp_ext *src)
27 { 27 {
28 struct fp_ext tmp, src2; 28 struct fp_ext tmp, src2;
29 int i, exp; 29 int i, exp;
30 30
31 dprint(PINSTR, "fsqrt\n"); 31 dprint(PINSTR, "fsqrt\n");
32 32
33 fp_monadic_check(dest, src); 33 fp_monadic_check(dest, src);
34 34
35 if (IS_ZERO(dest)) 35 if (IS_ZERO(dest))
36 return dest; 36 return dest;
37 37
38 if (dest->sign) { 38 if (dest->sign) {
39 fp_set_nan(dest); 39 fp_set_nan(dest);
40 return dest; 40 return dest;
41 } 41 }
42 if (IS_INF(dest)) 42 if (IS_INF(dest))
43 return dest; 43 return dest;
44 44
45 /* 45 /*
46 * sqrt(m) * 2^(p) 46 * sqrt(m) * 2^(p) , if e = 2*p
47 * sqrt(m*2^e) = 47 * sqrt(m*2^e) =
48 * sqrt(2*m) * 2^(p) 48 * sqrt(2*m) * 2^(p) , if e = 2*p + 1
49 * 49 *
50 * So we use the last bit of the expon 50 * So we use the last bit of the exponent to decide whether to
51 * use the m or 2*m. 51 * use the m or 2*m.
52 * 52 *
53 * Since only the fractional part of t 53 * Since only the fractional part of the mantissa is stored and
54 * the integer part is assumed to be o 54 * the integer part is assumed to be one, we place a 1 or 2 into
55 * the fixed point representation. 55 * the fixed point representation.
56 */ 56 */
57 exp = dest->exp; 57 exp = dest->exp;
58 dest->exp = 0x3FFF; 58 dest->exp = 0x3FFF;
59 if (!(exp & 1)) /* lowest bit 59 if (!(exp & 1)) /* lowest bit of exponent is set */
60 dest->exp++; 60 dest->exp++;
61 fp_copy_ext(&src2, dest); 61 fp_copy_ext(&src2, dest);
62 62
63 /* 63 /*
64 * The taylor row around a for sqrt(x) 64 * The taylor row around a for sqrt(x) is:
65 * sqrt(x) = sqrt(a) + 1/(2*sqrt( 65 * sqrt(x) = sqrt(a) + 1/(2*sqrt(a))*(x-a) + R
66 * With a=1 this gives: 66 * With a=1 this gives:
67 * sqrt(x) = 1 + 1/2*(x-1) 67 * sqrt(x) = 1 + 1/2*(x-1)
68 * = 1/2*(1+x) 68 * = 1/2*(1+x)
69 */ 69 */
70 /* It is safe to cast away the constne 70 /* It is safe to cast away the constness, as fp_one is normalized */
71 fp_fadd(dest, (struct fp_ext *)&fp_one 71 fp_fadd(dest, (struct fp_ext *)&fp_one);
72 dest->exp--; /* * 1/2 */ 72 dest->exp--; /* * 1/2 */
73 73
74 /* 74 /*
75 * We now apply the newton rule to the 75 * We now apply the newton rule to the function
76 * f(x) := x^2 - r 76 * f(x) := x^2 - r
77 * which has a null point on x = sqrt( 77 * which has a null point on x = sqrt(r).
78 * 78 *
79 * It gives: 79 * It gives:
80 * x' := x - f(x)/f'(x) 80 * x' := x - f(x)/f'(x)
81 * = x - (x^2 -r)/(2*x) 81 * = x - (x^2 -r)/(2*x)
82 * = x - (x - r/x)/2 82 * = x - (x - r/x)/2
83 * = (2*x - x + r/x)/2 83 * = (2*x - x + r/x)/2
84 * = (x + r/x)/2 84 * = (x + r/x)/2
85 */ 85 */
86 for (i = 0; i < 9; i++) { 86 for (i = 0; i < 9; i++) {
87 fp_copy_ext(&tmp, &src2); 87 fp_copy_ext(&tmp, &src2);
88 88
89 fp_fdiv(&tmp, dest); 89 fp_fdiv(&tmp, dest);
90 fp_fadd(dest, &tmp); 90 fp_fadd(dest, &tmp);
91 dest->exp--; 91 dest->exp--;
92 } 92 }
93 93
94 dest->exp += (exp - 0x3FFF) / 2; 94 dest->exp += (exp - 0x3FFF) / 2;
95 95
96 return dest; 96 return dest;
97 } 97 }
98 98
99 struct fp_ext *fp_fetoxm1(struct fp_ext *dest, 99 struct fp_ext *fp_fetoxm1(struct fp_ext *dest, struct fp_ext *src)
100 { 100 {
101 uprint("fetoxm1\n"); 101 uprint("fetoxm1\n");
102 102
103 fp_monadic_check(dest, src); 103 fp_monadic_check(dest, src);
104 104
105 return dest; 105 return dest;
106 } 106 }
107 107
108 struct fp_ext *fp_fetox(struct fp_ext *dest, s 108 struct fp_ext *fp_fetox(struct fp_ext *dest, struct fp_ext *src)
109 { 109 {
110 uprint("fetox\n"); 110 uprint("fetox\n");
111 111
112 fp_monadic_check(dest, src); 112 fp_monadic_check(dest, src);
113 113
114 return dest; 114 return dest;
115 } 115 }
116 116
117 struct fp_ext *fp_ftwotox(struct fp_ext *dest, 117 struct fp_ext *fp_ftwotox(struct fp_ext *dest, struct fp_ext *src)
118 { 118 {
119 uprint("ftwotox\n"); 119 uprint("ftwotox\n");
120 120
121 fp_monadic_check(dest, src); 121 fp_monadic_check(dest, src);
122 122
123 return dest; 123 return dest;
124 } 124 }
125 125
126 struct fp_ext *fp_ftentox(struct fp_ext *dest, 126 struct fp_ext *fp_ftentox(struct fp_ext *dest, struct fp_ext *src)
127 { 127 {
128 uprint("ftentox\n"); 128 uprint("ftentox\n");
129 129
130 fp_monadic_check(dest, src); 130 fp_monadic_check(dest, src);
131 131
132 return dest; 132 return dest;
133 } 133 }
134 134
135 struct fp_ext *fp_flogn(struct fp_ext *dest, s 135 struct fp_ext *fp_flogn(struct fp_ext *dest, struct fp_ext *src)
136 { 136 {
137 uprint("flogn\n"); 137 uprint("flogn\n");
138 138
139 fp_monadic_check(dest, src); 139 fp_monadic_check(dest, src);
140 140
141 return dest; 141 return dest;
142 } 142 }
143 143
144 struct fp_ext *fp_flognp1(struct fp_ext *dest, 144 struct fp_ext *fp_flognp1(struct fp_ext *dest, struct fp_ext *src)
145 { 145 {
146 uprint("flognp1\n"); 146 uprint("flognp1\n");
147 147
148 fp_monadic_check(dest, src); 148 fp_monadic_check(dest, src);
149 149
150 return dest; 150 return dest;
151 } 151 }
152 152
153 struct fp_ext *fp_flog10(struct fp_ext *dest, 153 struct fp_ext *fp_flog10(struct fp_ext *dest, struct fp_ext *src)
154 { 154 {
155 uprint("flog10\n"); 155 uprint("flog10\n");
156 156
157 fp_monadic_check(dest, src); 157 fp_monadic_check(dest, src);
158 158
159 return dest; 159 return dest;
160 } 160 }
161 161
162 struct fp_ext *fp_flog2(struct fp_ext *dest, s 162 struct fp_ext *fp_flog2(struct fp_ext *dest, struct fp_ext *src)
163 { 163 {
164 uprint("flog2\n"); 164 uprint("flog2\n");
165 165
166 fp_monadic_check(dest, src); 166 fp_monadic_check(dest, src);
167 167
168 return dest; 168 return dest;
169 } 169 }
170 170
171 struct fp_ext *fp_fgetexp(struct fp_ext *dest, 171 struct fp_ext *fp_fgetexp(struct fp_ext *dest, struct fp_ext *src)
172 { 172 {
173 dprint(PINSTR, "fgetexp\n"); 173 dprint(PINSTR, "fgetexp\n");
174 174
175 fp_monadic_check(dest, src); 175 fp_monadic_check(dest, src);
176 176
177 if (IS_INF(dest)) { 177 if (IS_INF(dest)) {
178 fp_set_nan(dest); 178 fp_set_nan(dest);
179 return dest; 179 return dest;
180 } 180 }
181 if (IS_ZERO(dest)) 181 if (IS_ZERO(dest))
182 return dest; 182 return dest;
183 183
184 fp_conv_long2ext(dest, (int)dest->exp 184 fp_conv_long2ext(dest, (int)dest->exp - 0x3FFF);
185 185
186 fp_normalize_ext(dest); 186 fp_normalize_ext(dest);
187 187
188 return dest; 188 return dest;
189 } 189 }
190 190
191 struct fp_ext *fp_fgetman(struct fp_ext *dest, 191 struct fp_ext *fp_fgetman(struct fp_ext *dest, struct fp_ext *src)
192 { 192 {
193 dprint(PINSTR, "fgetman\n"); 193 dprint(PINSTR, "fgetman\n");
194 194
195 fp_monadic_check(dest, src); 195 fp_monadic_check(dest, src);
196 196
197 if (IS_ZERO(dest)) 197 if (IS_ZERO(dest))
198 return dest; 198 return dest;
199 199
200 if (IS_INF(dest)) 200 if (IS_INF(dest))
201 return dest; 201 return dest;
202 202
203 dest->exp = 0x3FFF; 203 dest->exp = 0x3FFF;
204 204
205 return dest; 205 return dest;
206 } 206 }
207 207
208 208
Linux® is a registered trademark of Linus Torvalds in the United States and other countries.
TOMOYO® is a registered trademark of NTT DATA CORPORATION.