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