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 whether 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 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 220