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); 27 extern struct fp_ext *fp_fmul(struct fp_ext *dest, const struct fp_ext *src); 28 29 struct fp_ext * 30 fp_fsqrt(struct fp_ext *dest, struct fp_ext *src) 31 { 32 struct fp_ext tmp, src2; 33 int i, exp; 34 35 dprint(PINSTR, "fsqrt\n"); 36 37 fp_monadic_check(dest, src); 38 39 if (IS_ZERO(dest)) 40 return dest; 41 42 if (dest->sign) { 43 fp_set_nan(dest); 44 return dest; 45 } 46 if (IS_INF(dest)) 47 return dest; 48 49 /* 50 * sqrt(m) * 2^(p) , if e = 2*p 51 * sqrt(m*2^e) = 52 * sqrt(2*m) * 2^(p) , if e = 2*p + 1 53 * 54 * So we use the last bit of the exponent to decide wether to 55 * use the m or 2*m. 56 * 57 * Since only the fractional part of the mantissa is stored and 58 * the integer part is assumed to be one, we place a 1 or 2 into 59 * the fixed point representation. 60 */ 61 exp = dest->exp; 62 dest->exp = 0x3FFF; 63 if (!(exp & 1)) /* lowest bit of exponent is set */ 64 dest->exp++; 65 fp_copy_ext(&src2, dest); 66 67 /* 68 * The taylor row around a for sqrt(x) is: 69 * sqrt(x) = sqrt(a) + 1/(2*sqrt(a))*(x-a) + R 70 * With a=1 this gives: 71 * sqrt(x) = 1 + 1/2*(x-1) 72 * = 1/2*(1+x) 73 */ 74 fp_fadd(dest, &fp_one); 75 dest->exp--; /* * 1/2 */ 76 77 /* 78 * We now apply the newton rule to the function 79 * f(x) := x^2 - r 80 * which has a null point on x = sqrt(r). 81 * 82 * It gives: 83 * x' := x - f(x)/f'(x) 84 * = x - (x^2 -r)/(2*x) 85 * = x - (x - r/x)/2 86 * = (2*x - x + r/x)/2 87 * = (x + r/x)/2 88 */ 89 for (i = 0; i < 9; i++) { 90 fp_copy_ext(&tmp, &src2); 91 92 fp_fdiv(&tmp, dest); 93 fp_fadd(dest, &tmp); 94 dest->exp--; 95 } 96 97 dest->exp += (exp - 0x3FFF) / 2; 98 99 return dest; 100 } 101 102 struct fp_ext * 103 fp_fetoxm1(struct fp_ext *dest, struct fp_ext *src) 104 { 105 uprint("fetoxm1\n"); 106 107 fp_monadic_check(dest, src); 108 109 if (IS_ZERO(dest)) 110 return dest; 111 112 return dest; 113 } 114 115 struct fp_ext * 116 fp_fetox(struct fp_ext *dest, struct fp_ext *src) 117 { 118 uprint("fetox\n"); 119 120 fp_monadic_check(dest, src); 121 122 return dest; 123 } 124 125 struct fp_ext * 126 fp_ftwotox(struct fp_ext *dest, struct fp_ext *src) 127 { 128 uprint("ftwotox\n"); 129 130 fp_monadic_check(dest, src); 131 132 return dest; 133 } 134 135 struct fp_ext * 136 fp_ftentox(struct fp_ext *dest, struct fp_ext *src) 137 { 138 uprint("ftentox\n"); 139 140 fp_monadic_check(dest, src); 141 142 return dest; 143 } 144 145 struct fp_ext * 146 fp_flogn(struct fp_ext *dest, struct fp_ext *src) 147 { 148 uprint("flogn\n"); 149 150 fp_monadic_check(dest, src); 151 152 return dest; 153 } 154 155 struct fp_ext * 156 fp_flognp1(struct fp_ext *dest, struct fp_ext *src) 157 { 158 uprint("flognp1\n"); 159 160 fp_monadic_check(dest, src); 161 162 return dest; 163 } 164 165 struct fp_ext * 166 fp_flog10(struct fp_ext *dest, struct fp_ext *src) 167 { 168 uprint("flog10\n"); 169 170 fp_monadic_check(dest, src); 171 172 return dest; 173 } 174 175 struct fp_ext * 176 fp_flog2(struct fp_ext *dest, struct fp_ext *src) 177 { 178 uprint("flog2\n"); 179 180 fp_monadic_check(dest, src); 181 182 return dest; 183 } 184 185 struct fp_ext * 186 fp_fgetexp(struct fp_ext *dest, struct fp_ext *src) 187 { 188 dprint(PINSTR, "fgetexp\n"); 189 190 fp_monadic_check(dest, src); 191 192 if (IS_INF(dest)) { 193 fp_set_nan(dest); 194 return dest; 195 } 196 if (IS_ZERO(dest)) 197 return dest; 198 199 fp_conv_long2ext(dest, (int)dest->exp - 0x3FFF); 200 201 fp_normalize_ext(dest); 202 203 return dest; 204 } 205 206 struct fp_ext * 207 fp_fgetman(struct fp_ext *dest, struct fp_ext *src) 208 { 209 dprint(PINSTR, "fgetman\n"); 210 211 fp_monadic_check(dest, src); 212 213 if (IS_ZERO(dest)) 214 return dest; 215 216 if (IS_INF(dest)) 217 return dest; 218 219 dest->exp = 0x3FFF; 220 221 return dest; 222 } 223 224