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 20static const struct fp_ext fp_one = 21{ 22 .exp = 0x3fff, 23}; 24 25extern struct fp_ext *fp_fadd(struct fp_ext *dest, const struct fp_ext *src); 26extern struct fp_ext *fp_fdiv(struct fp_ext *dest, const struct fp_ext *src); 27 28struct fp_ext * 29fp_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 101struct fp_ext * 102fp_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 114struct fp_ext * 115fp_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 124struct fp_ext * 125fp_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 134struct fp_ext * 135fp_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 144struct fp_ext * 145fp_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 154struct fp_ext * 155fp_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 164struct fp_ext * 165fp_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 174struct fp_ext * 175fp_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 184struct fp_ext * 185fp_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 205struct fp_ext * 206fp_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