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