1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License, Version 1.0 only 6 * (the "License"). You may not use this file except in compliance 7 * with the License. 8 * 9 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 10 * or http://www.opensolaris.org/os/licensing. 11 * See the License for the specific language governing permissions 12 * and limitations under the License. 13 * 14 * When distributing Covered Code, include this CDDL HEADER in each 15 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 16 * If applicable, add the following below this CDDL HEADER, with the 17 * fields enclosed by brackets "[]" replaced with your own identifying 18 * information: Portions Copyright [yyyy] [name of copyright owner] 19 * 20 * CDDL HEADER END 21 */ 22 /* 23 * Copyright (c) 1994-1997, by Sun Microsystems, Inc. 24 * All rights reserved. 25 */ 26 27 /* 28 * This file contains __quad_mag_add and __quad_mag_sub, the core 29 * of the quad precision add and subtract operations. 30 */ 31 32 #include "quad.h" 33 34 /* 35 * __quad_mag_add(x, y, z, fsr) 36 * 37 * Sets *z = *x + *y, rounded according to the rounding mode in *fsr, 38 * and updates the current exceptions in *fsr. This routine assumes 39 * *x and *y are finite, with the same sign (i.e., an addition of 40 * magnitudes), |*x| >= |*y|, and *z already has its sign bit set. 41 */ 42 void 43 __quad_mag_add(const union longdouble *x, const union longdouble *y, 44 union longdouble *z, unsigned int *fsr) 45 { 46 unsigned int lx, ly, ex, ey, frac2, frac3, frac4; 47 unsigned int round, sticky, carry, rm; 48 int e, uflo; 49 50 /* get the leading significand words and exponents */ 51 ex = (x->l.msw & 0x7fffffff) >> 16; 52 lx = x->l.msw & 0xffff; 53 if (ex == 0) 54 ex = 1; 55 else 56 lx |= 0x10000; 57 58 ey = (y->l.msw & 0x7fffffff) >> 16; 59 ly = y->l.msw & 0xffff; 60 if (ey == 0) 61 ey = 1; 62 else 63 ly |= 0x10000; 64 65 /* prenormalize y */ 66 e = (int) ex - (int) ey; 67 round = sticky = 0; 68 if (e >= 114) { 69 frac2 = x->l.frac2; 70 frac3 = x->l.frac3; 71 frac4 = x->l.frac4; 72 sticky = ly | y->l.frac2 | y->l.frac3 | y->l.frac4; 73 } else { 74 frac2 = y->l.frac2; 75 frac3 = y->l.frac3; 76 frac4 = y->l.frac4; 77 if (e >= 96) { 78 sticky = frac4 | frac3 | (frac2 & 0x7fffffff); 79 round = frac2 & 0x80000000; 80 frac4 = ly; 81 frac3 = frac2 = ly = 0; 82 e -= 96; 83 } else if (e >= 64) { 84 sticky = frac4 | (frac3 & 0x7fffffff); 85 round = frac3 & 0x80000000; 86 frac4 = frac2; 87 frac3 = ly; 88 frac2 = ly = 0; 89 e -= 64; 90 } else if (e >= 32) { 91 sticky = frac4 & 0x7fffffff; 92 round = frac4 & 0x80000000; 93 frac4 = frac3; 94 frac3 = frac2; 95 frac2 = ly; 96 ly = 0; 97 e -= 32; 98 } 99 if (e) { 100 sticky |= round | (frac4 & ((1 << (e - 1)) - 1)); 101 round = frac4 & (1 << (e - 1)); 102 frac4 = (frac4 >> e) | (frac3 << (32 - e)); 103 frac3 = (frac3 >> e) | (frac2 << (32 - e)); 104 frac2 = (frac2 >> e) | (ly << (32 - e)); 105 ly >>= e; 106 } 107 108 /* add, propagating carries */ 109 frac4 += x->l.frac4; 110 carry = (frac4 < x->l.frac4); 111 frac3 += x->l.frac3; 112 if (carry) { 113 frac3++; 114 carry = (frac3 <= x->l.frac3); 115 } else { 116 carry = (frac3 < x->l.frac3); 117 } 118 frac2 += x->l.frac2; 119 if (carry) { 120 frac2++; 121 carry = (frac2 <= x->l.frac2); 122 } else { 123 carry = (frac2 < x->l.frac2); 124 } 125 lx += ly; 126 if (carry) 127 lx++; 128 129 /* postnormalize */ 130 if (lx >= 0x20000) { 131 sticky |= round; 132 round = frac4 & 1; 133 frac4 = (frac4 >> 1) | (frac3 << 31); 134 frac3 = (frac3 >> 1) | (frac2 << 31); 135 frac2 = (frac2 >> 1) | (lx << 31); 136 lx >>= 1; 137 ex++; 138 } 139 } 140 141 /* keep track of whether the result before rounding is tiny */ 142 uflo = (lx < 0x10000); 143 144 /* get the rounding mode, fudging directed rounding modes */ 145 /* as though the result were positive */ 146 rm = *fsr >> 30; 147 if (z->l.msw) 148 rm ^= (rm >> 1); 149 150 /* see if we need to round */ 151 if (round | sticky) { 152 *fsr |= FSR_NXC; 153 154 /* round up if necessary */ 155 if (rm == FSR_RP || (rm == FSR_RN && round && 156 (sticky || (frac4 & 1)))) { 157 if (++frac4 == 0) 158 if (++frac3 == 0) 159 if (++frac2 == 0) 160 if (++lx >= 0x20000) { 161 lx >>= 1; 162 ex++; 163 } 164 } 165 } 166 167 /* check for overflow */ 168 if (ex >= 0x7fff) { 169 /* store the default overflowed result */ 170 *fsr |= FSR_OFC | FSR_NXC; 171 if (rm == FSR_RN || rm == FSR_RP) { 172 z->l.msw |= 0x7fff0000; 173 z->l.frac2 = z->l.frac3 = z->l.frac4 = 0; 174 } else { 175 z->l.msw |= 0x7ffeffff; 176 z->l.frac2 = z->l.frac3 = z->l.frac4 = 0xffffffff; 177 } 178 } else { 179 /* store the result */ 180 if (lx >= 0x10000) 181 z->l.msw |= (ex << 16); 182 z->l.msw |= (lx & 0xffff); 183 z->l.frac2 = frac2; 184 z->l.frac3 = frac3; 185 z->l.frac4 = frac4; 186 187 /* if the pre-rounded result was tiny and underflow trapping */ 188 /* is enabled, simulate underflow */ 189 if (uflo && (*fsr & FSR_UFM)) 190 *fsr |= FSR_UFC; 191 } 192 } 193 194 /* 195 * __quad_mag_sub(x, y, z, fsr) 196 * 197 * Sets *z = *x - *y, rounded according to the rounding mode in *fsr, 198 * and updates the current exceptions in *fsr. This routine assumes 199 * *x and *y are finite, with opposite signs (i.e., a subtraction of 200 * magnitudes), |*x| >= |*y|, and *z already has its sign bit set. 201 */ 202 void 203 __quad_mag_sub(const union longdouble *x, const union longdouble *y, 204 union longdouble *z, unsigned int *fsr) 205 { 206 unsigned int lx, ly, ex, ey, frac2, frac3, frac4; 207 unsigned int guard, round, sticky, borrow, rm; 208 int e; 209 210 /* get the leading significand words and exponents */ 211 ex = (x->l.msw & 0x7fffffff) >> 16; 212 lx = x->l.msw & 0xffff; 213 if (ex == 0) 214 ex = 1; 215 else 216 lx |= 0x10000; 217 218 ey = (y->l.msw & 0x7fffffff) >> 16; 219 ly = y->l.msw & 0xffff; 220 if (ey == 0) 221 ey = 1; 222 else 223 ly |= 0x10000; 224 225 /* prenormalize y */ 226 e = (int) ex - (int) ey; 227 guard = round = sticky = 0; 228 if (e > 114) { 229 sticky = ly | y->l.frac2 | y->l.frac3 | y->l.frac4; 230 ly = frac2 = frac3 = frac4 = 0; 231 } else { 232 frac2 = y->l.frac2; 233 frac3 = y->l.frac3; 234 frac4 = y->l.frac4; 235 if (e >= 96) { 236 sticky = frac4 | frac3 | (frac2 & 0x3fffffff); 237 round = frac2 & 0x40000000; 238 guard = frac2 & 0x80000000; 239 frac4 = ly; 240 frac3 = frac2 = ly = 0; 241 e -= 96; 242 } else if (e >= 64) { 243 sticky = frac4 | (frac3 & 0x3fffffff); 244 round = frac3 & 0x40000000; 245 guard = frac3 & 0x80000000; 246 frac4 = frac2; 247 frac3 = ly; 248 frac2 = ly = 0; 249 e -= 64; 250 } else if (e >= 32) { 251 sticky = frac4 & 0x3fffffff; 252 round = frac4 & 0x40000000; 253 guard = frac4 & 0x80000000; 254 frac4 = frac3; 255 frac3 = frac2; 256 frac2 = ly; 257 ly = 0; 258 e -= 32; 259 } 260 if (e > 1) { 261 sticky |= guard | round | 262 (frac4 & ((1 << (e - 2)) - 1)); 263 round = frac4 & (1 << (e - 2)); 264 guard = frac4 & (1 << (e - 1)); 265 frac4 = (frac4 >> e) | (frac3 << (32 - e)); 266 frac3 = (frac3 >> e) | (frac2 << (32 - e)); 267 frac2 = (frac2 >> e) | (ly << (32 - e)); 268 ly >>= e; 269 } else if (e == 1) { 270 sticky |= round; 271 round = guard; 272 guard = frac4 & 1; 273 frac4 = (frac4 >> 1) | (frac3 << 31); 274 frac3 = (frac3 >> 1) | (frac2 << 31); 275 frac2 = (frac2 >> 1) | (ly << 31); 276 ly >>= 1; 277 } 278 } 279 280 /* complement guard, round, and sticky as need be */ 281 if (sticky) { 282 round = !round; 283 guard = !guard; 284 } else if (round) { 285 guard = !guard; 286 } 287 borrow = (guard | round | sticky); 288 289 /* subtract, propagating borrows */ 290 frac4 = x->l.frac4 - frac4; 291 if (borrow) { 292 frac4--; 293 borrow = (frac4 >= x->l.frac4); 294 } else { 295 borrow = (frac4 > x->l.frac4); 296 } 297 frac3 = x->l.frac3 - frac3; 298 if (borrow) { 299 frac3--; 300 borrow = (frac3 >= x->l.frac3); 301 } else { 302 borrow = (frac3 > x->l.frac3); 303 } 304 frac2 = x->l.frac2 - frac2; 305 if (borrow) { 306 frac2--; 307 borrow = (frac2 >= x->l.frac2); 308 } else { 309 borrow = (frac2 > x->l.frac2); 310 } 311 lx -= ly; 312 if (borrow) 313 lx--; 314 315 /* get the rounding mode */ 316 rm = *fsr >> 30; 317 318 /* handle zero result */ 319 if (!(lx | frac2 | frac3 | frac4 | guard)) { 320 z->l.msw = ((rm == FSR_RM)? 0x80000000 : 0); 321 z->l.frac2 = z->l.frac3 = z->l.frac4 = 0; 322 return; 323 } 324 325 /* postnormalize */ 326 if (lx < 0x10000) { 327 /* if cancellation occurred or the exponent is 1, */ 328 /* the result is exact */ 329 if (lx < 0x8000 || ex == 1) { 330 while ((lx | (frac2 & 0xfffe0000)) == 0 && ex > 32) { 331 lx = frac2; 332 frac2 = frac3; 333 frac3 = frac4; 334 frac4 = ((guard)? 0x80000000 : 0); 335 guard = 0; 336 ex -= 32; 337 } 338 while (lx < 0x10000 && ex > 1) { 339 lx = (lx << 1) | (frac2 >> 31); 340 frac2 = (frac2 << 1) | (frac3 >> 31); 341 frac3 = (frac3 << 1) | (frac4 >> 31); 342 frac4 <<= 1; 343 if (guard) { 344 frac4 |= 1; 345 guard = 0; 346 } 347 ex--; 348 } 349 if (lx >= 0x10000) 350 z->l.msw |= (ex << 16); 351 z->l.msw |= (lx & 0xffff); 352 z->l.frac2 = frac2; 353 z->l.frac3 = frac3; 354 z->l.frac4 = frac4; 355 356 /* if the result is tiny and underflow trapping is */ 357 /* enabled, simulate underflow */ 358 if (lx < 0x10000 && (*fsr & FSR_UFM)) 359 *fsr |= FSR_UFC; 360 return; 361 } 362 363 /* otherwise we only borrowed one place */ 364 lx = (lx << 1) | (frac2 >> 31); 365 frac2 = (frac2 << 1) | (frac3 >> 31); 366 frac3 = (frac3 << 1) | (frac4 >> 31); 367 frac4 <<= 1; 368 if (guard) 369 frac4 |= 1; 370 ex--; 371 } else { 372 sticky |= round; 373 round = guard; 374 } 375 376 /* fudge directed rounding modes as though the result were positive */ 377 if (z->l.msw) 378 rm ^= (rm >> 1); 379 380 /* see if we need to round */ 381 if (round | sticky) { 382 *fsr |= FSR_NXC; 383 384 /* round up if necessary */ 385 if (rm == FSR_RP || (rm == FSR_RN && round && 386 (sticky || (frac4 & 1)))) { 387 if (++frac4 == 0) 388 if (++frac3 == 0) 389 if (++frac2 == 0) 390 if (++lx >= 0x20000) { 391 lx >>= 1; 392 ex++; 393 } 394 } 395 } 396 397 /* store the result */ 398 z->l.msw |= (ex << 16) | (lx & 0xffff); 399 z->l.frac2 = frac2; 400 z->l.frac3 = frac3; 401 z->l.frac4 = frac4; 402 } 403