Changeset 1016
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/real.cpp
r1013 r1016 23 23 { 24 24 25 static int const BIGIT_BITS = 32; 26 25 27 real::real(float f) { *this = (double)f; } 26 28 real::real(int i) { *this = (double)i; } … … 47 49 } 48 50 49 m_mantissa[0] = u.x >> 36; 50 m_mantissa[1] = u.x >> 20; 51 m_mantissa[2] = u.x >> 4; 52 m_mantissa[3] = u.x << 12; 53 memset(m_mantissa + 4, 0, sizeof(m_mantissa) - 4 * sizeof(m_mantissa[0])); 51 m_mantissa[0] = u.x >> 20; 52 m_mantissa[1] = u.x << 12; 53 memset(m_mantissa + 2, 0, (BIGITS - 2) * sizeof(m_mantissa[0])); 54 54 } 55 55 … … 82 82 83 83 /* Store mantissa if necessary */ 84 u.x <<= 16;84 u.x <<= 32; 85 85 u.x |= m_mantissa[0]; 86 u.x <<= 16; 87 u.x |= m_mantissa[1]; 88 u.x <<= 16; 89 u.x |= m_mantissa[2]; 90 u.x <<= 4; 91 u.x |= m_mantissa[3] >> 12; 86 u.x <<= 20; 87 u.x |= m_mantissa[1] >> 12; 92 88 /* Rounding */ 93 u.x += (m_mantissa[ 3] >> 11) & 1;89 u.x += (m_mantissa[1] >> 11) & 1; 94 90 } 95 91 … … 132 128 int e2 = x.m_signexp - (1 << 30) + 1; 133 129 134 int bigoff = (e1 - e2) / (sizeof(uint16_t) * 8);135 int off = e1 - e2 - bigoff * (sizeof(uint16_t) * 8);130 int bigoff = (e1 - e2) / BIGIT_BITS; 131 int off = e1 - e2 - bigoff * BIGIT_BITS; 136 132 137 133 if (bigoff > BIGITS) … … 140 136 ret.m_signexp = m_signexp; 141 137 142 uint 32_t carry = 0;138 uint64_t carry = 0; 143 139 for (int i = BIGITS; i--; ) 144 140 { … … 147 143 carry += x.m_mantissa[i - bigoff] >> off; 148 144 149 if ( i - bigoff > 0)150 carry += (x.m_mantissa[i - bigoff - 1] << ( 16 - off)) & 0xffffu;145 if (off && i - bigoff > 0) 146 carry += (x.m_mantissa[i - bigoff - 1] << (BIGIT_BITS - off)) & 0xffffffffu; 151 147 else if (i - bigoff == 0) 152 carry += 0x0001u << (16- off);148 carry += (uint64_t)1 << (BIGIT_BITS - off); 153 149 154 150 ret.m_mantissa[i] = carry; 155 carry >>= 16;151 carry >>= BIGIT_BITS; 156 152 } 157 153 … … 162 158 for (int i = 0; i < BIGITS; i++) 163 159 { 164 uint 16_t tmp = ret.m_mantissa[i];165 ret.m_mantissa[i] = (carry << 15) | (tmp >> 1);166 carry = tmp & 0x0001u;160 uint32_t tmp = ret.m_mantissa[i]; 161 ret.m_mantissa[i] = (carry << (BIGIT_BITS - 1)) | (tmp >> 1); 162 carry = tmp & 1u; 167 163 } 168 164 ret.m_signexp++; … … 194 190 int e2 = x.m_signexp - (1 << 30) + 1; 195 191 196 int bigoff = (e1 - e2) / (sizeof(uint16_t) * 8);197 int off = e1 - e2 - bigoff * (sizeof(uint16_t) * 8);192 int bigoff = (e1 - e2) / BIGIT_BITS; 193 int off = e1 - e2 - bigoff * BIGIT_BITS; 198 194 199 195 if (bigoff > BIGITS) … … 202 198 ret.m_signexp = m_signexp; 203 199 204 int 32_t carry = 0;200 int64_t carry = 0; 205 201 for (int i = 0; i < bigoff; i++) 206 202 { 207 203 carry -= x.m_mantissa[BIGITS - i]; 208 carry = (carry & 0xffff0000u) | (carry >> 16); 209 } 210 carry -= x.m_mantissa[BIGITS - 1 - bigoff] & ((1 << off) - 1); 211 carry /= (1 << off); 204 /* Emulates a signed shift */ 205 carry >>= BIGIT_BITS; 206 carry |= carry << BIGIT_BITS; 207 } 208 carry -= x.m_mantissa[BIGITS - 1 - bigoff] & (((int64_t)1 << off) - 1); 209 carry /= (int64_t)1 << off; 212 210 213 211 for (int i = BIGITS; i--; ) … … 217 215 carry -= x.m_mantissa[i - bigoff] >> off; 218 216 219 if ( i - bigoff > 0)220 carry -= (x.m_mantissa[i - bigoff - 1] << ( 16 - off)) & 0xffffu;217 if (off && i - bigoff > 0) 218 carry -= (x.m_mantissa[i - bigoff - 1] << (BIGIT_BITS - off)) & 0xffffffffu; 221 219 else if (i - bigoff == 0) 222 carry -= 0x0001u << (16- off);220 carry -= (uint64_t)1 << (BIGIT_BITS - off); 223 221 224 222 ret.m_mantissa[i] = carry; 225 carry = (carry & 0xffff0000u) | (carry >> 16); 223 carry >>= BIGIT_BITS; 224 carry |= carry << BIGIT_BITS; 226 225 } 227 226 … … 238 237 if (!ret.m_mantissa[i]) 239 238 { 240 off += sizeof(uint16_t) * 8;239 off += BIGIT_BITS; 241 240 continue; 242 241 } 243 242 244 for (uint 16_t tmp = ret.m_mantissa[i]; tmp < 0x8000u; tmp <<= 1)243 for (uint32_t tmp = ret.m_mantissa[i]; tmp < 0x80000000u; tmp <<= 1) 245 244 off++; 246 245 break; 247 246 } 248 if (off == BIGITS * sizeof(uint16_t) * 8)247 if (off == BIGITS * BIGIT_BITS) 249 248 ret.m_signexp &= 0x80000000u; 250 249 else … … 253 252 ret.m_signexp -= off; 254 253 255 bigoff = off / (sizeof(uint16_t) * 8);256 off -= bigoff * sizeof(uint16_t) * 8;254 bigoff = off / BIGIT_BITS; 255 off -= bigoff * BIGIT_BITS; 257 256 258 257 for (int i = 0; i < BIGITS; i++) 259 258 { 260 uint 16_t tmp = 0;259 uint32_t tmp = 0; 261 260 if (i + bigoff < BIGITS) 262 261 tmp |= ret.m_mantissa[i + bigoff] << off; 263 if ( i + bigoff + 1 < BIGITS)264 tmp |= ret.m_mantissa[i + bigoff + 1] >> ( 16- off);262 if (off && i + bigoff + 1 < BIGITS) 263 tmp |= ret.m_mantissa[i + bigoff + 1] >> (BIGIT_BITS - off); 265 264 ret.m_mantissa[i] = tmp; 266 265 } … … 288 287 /* Accumulate low order product; no need to store it, we just 289 288 * want the carry value */ 290 uint64_t carry = 0 ;289 uint64_t carry = 0, hicarry = 0, prev; 291 290 for (int i = 0; i < BIGITS; i++) 292 291 { 293 292 for (int j = 0; j < i + 1; j++) 294 carry += (uint32_t)m_mantissa[BIGITS - 1 - j] 295 * (uint32_t)x.m_mantissa[BIGITS - 1 + j - i]; 296 carry >>= 16; 293 { 294 prev = carry; 295 carry += (uint64_t)m_mantissa[BIGITS - 1 - j] 296 * (uint64_t)x.m_mantissa[BIGITS - 1 + j - i]; 297 if (carry < prev) 298 hicarry++; 299 } 300 carry >>= BIGIT_BITS; 301 carry |= hicarry << BIGIT_BITS; 302 hicarry >>= BIGIT_BITS; 297 303 } 298 304 … … 300 306 { 301 307 for (int j = i + 1; j < BIGITS; j++) 302 carry += (uint32_t)m_mantissa[BIGITS - 1 - j] 303 * (uint32_t)x.m_mantissa[j - 1 - i]; 304 308 { 309 prev = carry; 310 carry += (uint64_t)m_mantissa[BIGITS - 1 - j] 311 * (uint64_t)x.m_mantissa[j - 1 - i]; 312 if (carry < prev) 313 hicarry++; 314 } 315 prev = carry; 305 316 carry += m_mantissa[BIGITS - 1 - i]; 306 317 carry += x.m_mantissa[BIGITS - 1 - i]; 307 ret.m_mantissa[BIGITS - 1 - i] = carry & 0xffffu; 308 carry >>= 16; 318 if (carry < prev) 319 hicarry++; 320 ret.m_mantissa[BIGITS - 1 - i] = carry & 0xffffffffu; 321 carry >>= BIGIT_BITS; 322 carry |= hicarry << BIGIT_BITS; 323 hicarry >>= BIGIT_BITS; 309 324 } 310 325 … … 315 330 for (int i = 0; i < BIGITS; i++) 316 331 { 317 uint 16_t tmp = ret.m_mantissa[i];318 ret.m_mantissa[i] = (carry << 15) | (tmp >> 1);319 carry = tmp & 0x0001u;332 uint32_t tmp = ret.m_mantissa[i]; 333 ret.m_mantissa[i] = (carry << (BIGIT_BITS - 1)) | (tmp >> 1); 334 carry = tmp & 1u; 320 335 } 321 336 e++; … … 473 488 /* Use the system's float inversion to approximate 1/x */ 474 489 union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f }; 475 v.x |= (uint32_t)x.m_mantissa[0] << 7; 476 v.x |= (uint32_t)x.m_mantissa[1] >> 9; 490 v.x |= x.m_mantissa[0] >> 9; 477 491 v.f = 1.0 / v.f; 478 492 479 493 real ret; 480 ret.m_mantissa[0] = (v.x >> 7) & 0xffffu; 481 ret.m_mantissa[1] = (v.x << 9) & 0xffffu; 494 ret.m_mantissa[0] = v.x << 9; 482 495 483 496 uint32_t sign = x.m_signexp & 0x80000000u; … … 490 503 /* FIXME: log2(BIGITS) steps of Newton-Raphson seems to be enough for 491 504 * convergence, but this hasn't been checked seriously. */ 492 for (int i = 1; i < real::BIGITS; i *= 2)505 for (int i = 2; i < real::BIGITS; i *= 2) 493 506 ret = ret * (real::R_2 - ret * x); 494 507 … … 518 531 union { float f; uint32_t x; } u = { 1.0f }, v = { 2.0f }; 519 532 v.x -= ((x.m_signexp & 1) << 23); 520 v.x |= (uint32_t)x.m_mantissa[0] << 7; 521 v.x |= (uint32_t)x.m_mantissa[1] >> 9; 533 v.x |= x.m_mantissa[0] >> 9; 522 534 v.f = 1.0 / sqrtf(v.f); 523 535 524 536 real ret; 525 ret.m_mantissa[0] = (v.x >> 7) & 0xffffu; 526 ret.m_mantissa[1] = (v.x << 9) & 0xffffu; 537 ret.m_mantissa[0] = v.x << 9; 527 538 528 539 uint32_t sign = x.m_signexp & 0x80000000u; … … 536 547 /* FIXME: log2(BIGITS) steps of Newton-Raphson seems to be enough for 537 548 * convergence, but this hasn't been checked seriously. */ 538 for (int i = 1; i < real::BIGITS; i *= 2)549 for (int i = 2; i < real::BIGITS; i *= 2) 539 550 { 540 551 ret = ret * (real::R_3 - ret * ret * x); … … 593 604 { 594 605 tmp.m_signexp = 0xffffffffu; 595 tmp.m_mantissa[0] = 0xffff u;606 tmp.m_mantissa[0] = 0xffffffffu; 596 607 return tmp; 597 608 } … … 658 669 if (exponent <= 0) 659 670 ret.m_mantissa[i] = 0; 660 else if (exponent < 8 * (int)sizeof(uint16_t))661 ret.m_mantissa[i] &= ~((1 << ( 16- exponent)) - 1);662 663 exponent -= 8 * sizeof(uint16_t);671 else if (exponent < BIGIT_BITS) 672 ret.m_mantissa[i] &= ~((1 << (BIGIT_BITS - exponent)) - 1); 673 674 exponent -= BIGIT_BITS; 664 675 } 665 676 -
trunk/src/real.h
r1003 r1016 109 109 * Newton-Raphson iterations) and real::print (the number of printed 110 110 * digits) */ 111 static int const BIGITS = 32;111 static int const BIGITS = 16; 112 112 113 113 uint32_t m_size; 114 114 uint32_t m_signexp; 115 uint 16_t m_mantissa[BIGITS];115 uint32_t m_mantissa[BIGITS]; 116 116 }; 117 117
Note: See TracChangeset
for help on using the changeset viewer.