/* * Copyright (c) 2005-2014 Rich Felker, et al. * * Permission is hereby granted, free of charge, to any person obtaining * a copy of this software and associated documentation files (the * "Software"), to deal in the Software without restriction, including * without limitation the rights to use, copy, modify, merge, publish, * distribute, sublicense, and/or sell copies of the Software, and to * permit persons to whom the Software is furnished to do so, subject to * the following conditions: * * The above copyright notice and this permission notice shall be * included in all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
staticdouble decfloat(FFFILE *f, int c, int bits, int emin, int sign, int pok)
{
uint32_t x[KMAX]; staticconst uint32_t th[] = { LD_B1B_MAX }; int i, j, k, a, z; longlong lrp=0, dc=0; longlong e10=0; int lnz = 0; int gotdig = 0, gotrad = 0; int rp; int e2; int emax = -emin-bits+3; int denormal = 0; double y; double frac=0; double bias=0; staticconstint p10s[] = { 10, 100, 1000, 10000,
100000, 1000000, 10000000, 100000000 };
j=0;
k=0;
/* Don't let leading zeros consume buffer space */ for (; c=='0'; c = shgetc(f)) gotdig=1; if (c=='.') {
gotrad = 1; for (c = shgetc(f); c=='0'; c = shgetc(f)) gotdig=1, lrp--;
}
x[0] = 0; for (; c-'0'<10U || c=='.'; c = shgetc(f)) { if (c == '.') { if (gotrad) break;
gotrad = 1;
lrp = dc;
} elseif (k < KMAX-3) {
dc++; if (c!='0') lnz = dc; if (j) x[k] = x[k]*10 + c-'0'; else x[k] = c-'0'; if (++j==9) {
k++;
j=0;
}
gotdig=1;
} else {
dc++; if (c!='0') {
lnz = (KMAX-4)*9;
x[KMAX-4] |= 1;
}
}
} if (!gotrad) lrp=dc;
/* Align incomplete final B1B digit */ if (j) { for (; j<9; j++) x[k]*=10;
k++;
j=0;
}
a = 0;
z = k;
e2 = 0;
rp = lrp;
/* Optimize small to mid-size integers (even in exp. notation) */ if (lnz<9 && lnz<=rp && rp < 18) { int bitlim; if (rp == 9) return sign * (double)x[0]; if (rp < 9) return sign * (double)x[0] / p10s[8-rp];
bitlim = bits-3*(int)(rp-9); if (bitlim>30 || x[0]>>bitlim==0) return sign * (double)x[0] * p10s[rp-10];
}
/* Drop trailing zeros */ for (; !x[z-1]; z--);
/* Align radix point to B1B digit boundary */ if (rp % 9) { int rpm9 = rp>=0 ? rp%9 : rp%9+9; int p10 = p10s[8-rpm9];
uint32_t carry = 0; for (k=a; k!=z; k++) {
uint32_t tmp = x[k] % p10;
x[k] = x[k]/p10 + carry;
carry = 1000000000/p10 * tmp; if (k==a && !x[k]) {
a = (a+1 & MASK);
rp -= 9;
}
} if (carry) x[z++] = carry;
rp += 9-rpm9;
}
/* Upscale until desired number of bits are left of radix point */ while (rp < 9*LD_B1B_DIG || (rp == 9*LD_B1B_DIG && x[a]<th[0])) {
uint32_t carry = 0;
e2 -= 29; for (k=(z-1 & MASK); ; k=(k-1 & MASK)) {
uint64_t tmp = ((uint64_t)x[k] << 29) + carry; if (tmp > 1000000000) {
carry = tmp / 1000000000;
x[k] = tmp % 1000000000;
} else {
carry = 0;
x[k] = tmp;
} if (k==(z-1 & MASK) && k!=a && !x[k]) z = k; if (k==a) break;
} if (carry) {
rp += 9;
a = (a-1 & MASK); if (a == z) {
z = (z-1 & MASK);
x[z-1 & MASK] |= x[z];
}
x[a] = carry;
}
}
/* Downscale until exactly number of bits are left of radix point */ for (;;) {
uint32_t carry = 0; int sh = 1; for (i=0; i<LD_B1B_DIG; i++) {
k = (a+i & MASK); if (k == z || x[k] < th[i]) {
i=LD_B1B_DIG; break;
} if (x[a+i & MASK] > th[i]) break;
} if (i==LD_B1B_DIG && rp==9*LD_B1B_DIG) break; /* FIXME: find a way to compute optimal sh */ if (rp > 9+9*LD_B1B_DIG) sh = 9;
e2 += sh; for (k=a; k!=z; k=(k+1 & MASK)) {
uint32_t tmp = x[k] & (1<<sh)-1;
x[k] = (x[k]>>sh) + carry;
carry = (1000000000>>sh) * tmp; if (k==a && !x[k]) {
a = (a+1 & MASK);
i--;
rp -= 9;
}
} if (carry) { if ((z+1 & MASK) != a) {
x[z] = carry;
z = (z+1 & MASK);
} else x[z-1 & MASK] |= 1;
}
}
/* Assemble desired bits into floating point variable */ for (y=i=0; i<LD_B1B_DIG; i++) { if ((a+i & MASK)==z) x[(z=(z+1 & MASK))-1] = 0;
y = 1000000000.0L * y + x[a+i & MASK];
}
y *= sign;
/* Limit precision for denormal results */ if (bits > DBL_MANT_DIG+e2-emin) {
bits = DBL_MANT_DIG+e2-emin; if (bits<0) bits=0;
denormal = 1;
}
/* Calculate bias term to force rounding, move out lower bits */ if (bits < DBL_MANT_DIG) {
bias = copysign(scalbn(1, 2*DBL_MANT_DIG-bits-1), y);
frac = fmod(y, scalbn(1, DBL_MANT_DIG-bits));
y -= frac;
y += bias;
}
/* Process tail of decimal input so it can affect rounding */ if ((a+i & MASK) != z) {
uint32_t t = x[a+i & MASK]; if (t < 500000000 && (t || (a+i+1 & MASK) != z))
frac += 0.25*sign; elseif (t > 500000000)
frac += 0.75*sign; elseif (t == 500000000) { if ((a+i+1 & MASK) == z)
frac += 0.5*sign; else
frac += 0.75*sign;
} if (DBL_MANT_DIG-bits >= 2 && !fmod(frac, 1))
frac++;
}
y += frac;
y -= bias;
if ((e2+DBL_MANT_DIG & INT_MAX) > emax-5) { if (fabs(y) >= pow(2, DBL_MANT_DIG)) { if (denormal && bits==DBL_MANT_DIG+e2-emin)
denormal = 0;
y *= 0.5;
e2++;
} if (e2+DBL_MANT_DIG>emax || (denormal && frac))
errno = ERANGE;
}
return scalbn(y, e2);
}
staticdouble hexfloat(FFFILE *f, int bits, int emin, int sign, int pok)
{
uint32_t x = 0; double y = 0; double scale = 1; double bias = 0; int gottail = 0, gotrad = 0, gotdig = 0; longlong rp = 0; longlong dc = 0; longlong e2 = 0; int d; int c;
c = shgetc(f);
/* Skip leading zeros */ for (; c=='0'; c = shgetc(f))
gotdig = 1;
if (c=='.') {
gotrad = 1;
c = shgetc(f); /* Count zeros after the radix point before significand */ for (rp=0; c=='0'; c = shgetc(f), rp--) gotdig = 1;
}
for (; c-'0'<10U || (c|32)-'a'<6U || c=='.'; c = shgetc(f)) { if (c=='.') { if (gotrad) break;
rp = dc;
gotrad = 1;
} else {
gotdig = 1; if (c > '9') d = (c|32)+10-'a'; else d = c-'0'; if (dc<8) {
x = x*16 + d;
} elseif (dc < DBL_MANT_DIG/4+1) {
y += d*(scale/=16);
} elseif (d && !gottail) {
y += 0.5*scale;
gottail = 1;
}
dc++;
}
} if (!gotdig) {
shunget(f); if (pok) {
shunget(f); if (gotrad) shunget(f);
} else {
shlim(f, 0);
} return sign * 0.0;
} if (!gotrad) rp = dc; while (dc<8) x *= 16, dc++; if ((c|32)=='p') {
e2 = scanexp(f, pok); if (e2 == LLONG_MIN) { if (pok) {
shunget(f);
} else {
shlim(f, 0); return 0;
}
e2 = 0;
}
} else {
shunget(f);
}
e2 += 4*rp - 32;
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung und die Messung sind noch experimentell.