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 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 * 53 * So we use the last bit of the exponent to decide wether to 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 /* 67 * The taylor row around a for sqrt(x) is: 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 108 if (IS_ZERO(dest)) 109 return dest; 110 111 return dest; 112 } 113 114 struct fp_ext * 115 fp_fetox(struct fp_ext *dest, struct fp_ext *src) 116 { 117 uprint("fetox\n"); 118 119 fp_monadic_check(dest, src); 120 121 return dest; 122 } 123 124 struct fp_ext * 125 fp_ftwotox(struct fp_ext *dest, struct fp_ext *src) 126 { 127 uprint("ftwotox\n"); 128 129 fp_monadic_check(dest, src); 130 131 return dest; 132 } 133 134 struct fp_ext * 135 fp_ftentox(struct fp_ext *dest, struct fp_ext *src) 136 { 137 uprint("ftentox\n"); 138 139 fp_monadic_check(dest, src); 140 141 return dest; 142 } 143 144 struct fp_ext * 145 fp_flogn(struct fp_ext *dest, struct fp_ext *src) 146 { 147 uprint("flogn\n"); 148 149 fp_monadic_check(dest, src); 150 151 return dest; 152 } 153 154 struct fp_ext * 155 fp_flognp1(struct fp_ext *dest, struct fp_ext *src) 156 { 157 uprint("flognp1\n"); 158 159 fp_monadic_check(dest, src); 160 161 return dest; 162 } 163 164 struct fp_ext * 165 fp_flog10(struct fp_ext *dest, struct fp_ext *src) 166 { 167 uprint("flog10\n"); 168 169 fp_monadic_check(dest, src); 170 171 return dest; 172 } 173 174 struct fp_ext * 175 fp_flog2(struct fp_ext *dest, struct fp_ext *src) 176 { 177 uprint("flog2\n"); 178 179 fp_monadic_check(dest, src); 180 181 return dest; 182 } 183 184 struct fp_ext * 185 fp_fgetexp(struct fp_ext *dest, struct fp_ext *src) 186 { 187 dprint(PINSTR, "fgetexp\n"); 188 189 fp_monadic_check(dest, src); 190 191 if (IS_INF(dest)) { 192 fp_set_nan(dest); 193 return dest; 194 } 195 if (IS_ZERO(dest)) 196 return dest; 197 198 fp_conv_long2ext(dest, (int)dest->exp - 0x3FFF); 199 200 fp_normalize_ext(dest); 201 202 return dest; 203 } 204 205 struct fp_ext * 206 fp_fgetman(struct fp_ext *dest, struct fp_ext *src) 207 { 208 dprint(PINSTR, "fgetman\n"); 209 210 fp_monadic_check(dest, src); 211 212 if (IS_ZERO(dest)) 213 return dest; 214 215 if (IS_INF(dest)) 216 return dest; 217 218 dest->exp = 0x3FFF; 219 220 return dest; 221 } 222 223