aboutsummaryrefslogtreecommitdiff
path: root/src/base/fmt/float.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/base/fmt/float.c')
-rw-r--r--src/base/fmt/float.c507
1 files changed, 3 insertions, 504 deletions
diff --git a/src/base/fmt/float.c b/src/base/fmt/float.c
index ba2637e..c99630c 100644
--- a/src/base/fmt/float.c
+++ b/src/base/fmt/float.c
@@ -157,507 +157,6 @@ sub1(char *a, int n)
rt·exit(1);
}
-// -----------------------------------------------------------------------
-// strtod
-
-#define Nbits 28
-#define Nmant 53
-#define Prec ((Nmant+Nbits+1)/Nbits)
-
-#define Sigbit (1<<(Prec*Nbits-Nmant)) /* first significant bit of Prec-th word */
-#define Ndig 1500
-#define One (ulong)(1<<Nbits)
-#define Half (ulong)(One>>1)
-#define Maxe 310
-
-#define Fsign (1<<0) /* found - */
-#define Fesign (1<<1) /* found e- */
-#define Fdpoint (1<<2) /* found . */
-
-#define S0 0 /* _ _S0 +S1 #S2 .S3 */
-#define S1 1 /* _+ #S2 .S3 */
-#define S2 2 /* _+# #S2 .S4 eS5 */
-#define S3 3 /* _+. #S4 */
-#define S4 4 /* _+#.# #S4 eS5 */
-#define S5 5 /* _+#.#e +S6 #S7 */
-#define S6 6 /* _+#.#e+ #S7 */
-#define S7 7 /* _+#.#e+# #S7 */
-
-typedef struct Tab Tab;
-struct Tab
-{
- int bp;
- int siz;
- char *cmp;
-};
-
-static ulong
-umuldiv(ulong a, ulong b, ulong c)
-{
- double d;
-
- d = ((double)a * (double)b) / (double)c;
- if(d >= 4294967295.)
- d = 4294967295.;
- return (ulong)d;
-}
-
-static void
-frnorm(ulong *f)
-{
- int i, c;
-
- c = 0;
- for(i=Prec-1; i>0; i--) {
- f[i] += c;
- c = f[i] >> Nbits;
- f[i] &= One-1;
- }
- f[0] += c;
-}
-
-static int
-fpcmp(char *a, ulong* f)
-{
- ulong tf[Prec];
- int i, d, c;
-
- for(i=0; i<Prec; i++)
- tf[i] = f[i];
-
- for(;;) {
- /* tf *= 10 */
- for(i=0; i<Prec; i++)
- tf[i] = tf[i]*10;
- frnorm(tf);
- d = (tf[0] >> Nbits) + '0';
- tf[0] &= One-1;
-
- /* compare next digit */
- c = *a;
- if(c == 0) {
- if('0' < d)
- return -1;
- if(tf[0] != 0)
- goto cont;
- for(i=1; i<Prec; i++)
- if(tf[i] != 0)
- goto cont;
- return 0;
- }
- if(c > d)
- return +1;
- if(c < d)
- return -1;
- a++;
- cont:;
-}
-}
-
-static void
-divby(char *a, int *na, int b)
-{
- int n, c;
- char *p;
-
- p = a;
- n = 0;
- while(n>>b == 0){
- c = *a++;
- if(c == 0) {
- while(n) {
- c = n*10;
- if(c>>b)
- break;
- n = c;
- }
- goto xx;
- }
- n = n*10 + c-'0';
- (*na)--;
- }
- for(;;){
- c = n>>b;
- n -= c<<b;
- *p++ = c + '0';
- c = *a++;
- if(c == 0)
- break;
- n = n*10 + c-'0';
- }
- (*na)++;
- xx:
- while(n){
- n = n*10;
- c = n>>b;
- n -= c<<b;
- *p++ = c + '0';
- (*na)++;
- }
- *p = 0;
-}
-
-static Tab tab1[] =
-{
- 1, 0, "",
- 3, 1, "7",
- 6, 2, "63",
- 9, 3, "511",
- 13, 4, "8191",
- 16, 5, "65535",
- 19, 6, "524287",
- 23, 7, "8388607",
- 26, 8, "67108863",
- 27, 9, "134217727",
-};
-
-static void
-divascii(char *a, int *na, int *dp, int *bp)
-{
- int b, d;
- Tab *t;
-
- d = *dp;
- if(d >= (int)(arrlen(tab1)))
- d = (int)(arrlen(tab1))-1;
- t = tab1 + d;
- b = t->bp;
- if(mem·compare(a, t->siz, t->cmp) > 0)
- d--;
- *dp -= d;
- *bp += b;
- divby(a, na, b);
-}
-
-static void
-mulby(char *a, char *p, char *q, int b)
-{
- int n, c;
-
- n = 0;
- *p = 0;
- for(;;) {
- q--;
- if(q < a)
- break;
- c = *q - '0';
- c = (c<<b) + n;
- n = c/10;
- c -= n*10;
- p--;
- *p = c + '0';
- }
- while(n) {
- c = n;
- n = c/10;
- c -= n*10;
- p--;
- *p = c + '0';
- }
-}
-
-static Tab tab2[] =
-{
- 1, 1, "", /* dp = 0-0 */
- 3, 3, "125",
- 6, 5, "15625",
- 9, 7, "1953125",
- 13, 10, "1220703125",
- 16, 12, "152587890625",
- 19, 14, "19073486328125",
- 23, 17, "11920928955078125",
- 26, 19, "1490116119384765625",
- 27, 19, "7450580596923828125", /* dp 8-9 */
-};
-
-static void
-mulascii(char *a, int *na, int *dp, int *bp)
-{
- char *p;
- int d, b;
- Tab *t;
-
- d = -*dp;
- if(d >= (int)(arrlen(tab2)))
- d = (int)(arrlen(tab2))-1;
- t = tab2 + d;
- b = t->bp;
- if(mem·compare(a, t->siz, t->cmp) < 0)
- d--;
- p = a + *na;
- *bp -= b;
- *dp += d;
- *na += d;
- mulby(a, p+d, p, b);
-}
-
-static int
-cmp(char *a, char *b)
-{
- int c1, c2;
-
- while((c1 = *b++) != '\0') {
- c2 = *a++;
- if(utf8·isupper(c2))
- c2 = utf8·tolower(c2);
- if(c1 != c2)
- return 1;
- }
- return 0;
-}
-
-double
-fmtstrtod(char *as, char **aas)
-{
- int na, ex, dp, bp, c, i, flag, state;
- ulong low[Prec], hig[Prec], mid[Prec];
- double d;
- char *s, a[Ndig];
-
- flag = 0; /* Fsign, Fesign, Fdpoint */
- na = 0; /* number of digits of a[] */
- dp = 0; /* na of decimal point */
- ex = 0; /* exonent */
-
- state = S0;
- for(s=as;;s++){
- c = *s;
- if('0' <= c && c <= '9'){
- switch(state){
- case S0: case S1: case S2:
- state = S2;
- break;
- case S3: case S4:
- state = S4;
- break;
- case S5: case S6: case S7:
- state = S7;
- ex = ex*10 + (c-'0');
- continue;
- }
-
- if(na == 0 && c == '0'){
- dp--;
- continue;
- }
- if(na < Ndig-50)
- a[na++] = c;
- continue;
- }
- switch(c){
- case '\t': case '\n': case '\v': case '\f': case '\r': case ' ':
- if(state == S0)
- continue;
- break;
- case '-':
- if(state == S0)
- flag |= Fsign;
- else
- flag |= Fesign;
- case '+':
- if(state == S0)
- state = S1;
- else
- if(state == S5)
- state = S6;
- else
- break; /* syntax */
- continue;
- case '.':
- flag |= Fdpoint;
- dp = na;
- if(state == S0 || state == S1){
- state = S3;
- continue;
- }
- if(state == S2){
- state = S4;
- continue;
- }
- break;
- case 'e': case 'E':
- if(state == S2 || state == S4){
- state = S5;
- continue;
- }
- break;
- }
- break;
- }
-
- /* clean up return char-pointer */
- switch(state) {
- case S0:
- if(cmp(s, "nan") == 0){
- if(aas != nil)
- *aas = s+3;
- goto retnan;
- }
- case S1:
- if(cmp(s, "infinity") == 0){
- if(aas != nil)
- *aas = s+8;
- goto retinf;
- }
- if(cmp(s, "inf") == 0){
- if(aas != nil)
- *aas = s+3;
- goto retinf;
- }
- case S3:
- if(aas != nil)
- *aas = as;
- goto ret0; /* no digits found */
- case S6:
- s--; /* back over +- */
- case S5:
- s--; /* back over e */
- break;
- }
- if(aas != nil)
- *aas = s;
-
- if(flag & Fdpoint)
- while(na > 0 && a[na-1] == '0')
- na--;
- if(na == 0)
- goto ret0; /* zero */
- a[na] = 0;
- if(!(flag & Fdpoint))
- dp = na;
- if(flag & Fesign)
- ex = -ex;
- dp += ex;
- if(dp < -Maxe){
- /* errno = ERANGE; */
- goto ret0; /* underflow by exp */
- } else
- if(dp > +Maxe)
- goto retinf; /* overflow by exp */
-
- /*
- * normalize the decimal ascii number
- * to range .[5-9][0-9]* e0
- */
- bp = 0; /* binary exponent */
- while(dp > 0)
- divascii(a, &na, &dp, &bp);
- while(dp < 0 || a[0] < '5')
- mulascii(a, &na, &dp, &bp);
-
- /* close approx by naive conversion */
- mid[0] = 0;
- mid[1] = 1;
- for(i=0; (c=a[i]) != '\0'; i++) {
- mid[0] = mid[0]*10 + (c-'0');
- mid[1] = mid[1]*10;
- if(i >= 8)
- break;
- }
- low[0] = umuldiv(mid[0], One, mid[1]);
- hig[0] = umuldiv(mid[0]+1, One, mid[1]);
- for(i=1; i<Prec; i++) {
- low[i] = 0;
- hig[i] = One-1;
- }
-
- /* binary search for closest mantissa */
- for(;;) {
- /* mid = (hig + low) / 2 */
- c = 0;
- for(i=0; i<Prec; i++) {
- mid[i] = hig[i] + low[i];
- if(c)
- mid[i] += One;
- c = mid[i] & 1;
- mid[i] >>= 1;
- }
- frnorm(mid);
-
- /* compare */
- c = fpcmp(a, mid);
- if(c > 0) {
- c = 1;
- for(i=0; i<Prec; i++)
- if(low[i] != mid[i]) {
- c = 0;
- low[i] = mid[i];
- }
- if(c)
- break; /* between mid and hig */
- continue;
- }
- if(c < 0) {
- for(i=0; i<Prec; i++)
- hig[i] = mid[i];
- continue;
- }
-
- /* only hard part is if even/odd roundings wants to go up */
- c = mid[Prec-1] & (Sigbit-1);
- if(c == Sigbit/2 && (mid[Prec-1]&Sigbit) == 0)
- mid[Prec-1] -= c;
- break; /* exactly mid */
- }
-
- /* normal rounding applies */
- c = mid[Prec-1] & (Sigbit-1);
- mid[Prec-1] -= c;
- if(c >= Sigbit/2) {
- mid[Prec-1] += Sigbit;
- frnorm(mid);
- }
- goto out;
-
-ret0:
- return 0;
-
-retnan:
- return NaN();
-
-retinf:
- /* Unix strtod requires these. Plan 9 would return Inf(0) or Inf(-1). */
- /* errno = ERANGE; */
- if(flag & Fsign)
- return -HUGE_VAL;
- return HUGE_VAL;
-
-out:
- d = 0;
- for(i=0; i<Prec; i++)
- d = d*One + mid[i];
- if(flag & Fsign)
- d = -d;
- d = ldexp(d, bp - Prec*Nbits);
- if(d == 0) /* underflow */
- /* errno = ERANGE; */
-
- return d;
-}
-
-#undef Nbits
-#undef Nmant
-#undef Prec
-
-#undef Sigbit
-#undef Ndig
-#undef One
-#undef Half
-#undef Maxe
-
-#undef Fsign
-#undef Fesign
-#undef Fdpoint
-
-#undef S0
-#undef S1
-#undef S2
-#undef S3
-#undef S4
-#undef S5
-#undef S6
-#undef S7
-
static void
fmtexp(char *p, int e, int ucase)
{
@@ -712,7 +211,7 @@ dtoa(double f, char *s, int *exp, int *neg, int *len)
return;
}
- frexp(f, &e2);
+ math·frexp(f, &e2);
e = (int)(e2 * .301029995664);
g = f * fpow10(-e);
while(g < 1) {
@@ -736,7 +235,7 @@ dtoa(double f, char *s, int *exp, int *neg, int *len)
fmtexp(s+NSIGNIF, e, 0);
for(i=0; i<10; i++) {
- g=fmtstrtod(s, nil);
+ g=str·asfloat(s, nil);
if(f > g) {
if(add1(s, NSIGNIF)){
/* gained a digit */
@@ -763,7 +262,7 @@ dtoa(double f, char *s, int *exp, int *neg, int *len)
c = s[i];
if(c != '0'){
s[i] = '0';
- g=fmtstrtod(s, nil);
+ g=str·asfloat(s, nil);
if(g != f){
s[i] = c;
break;
generated by cgit v1.2.1 (git 2.18.0) at 2024-07-03 13:52:40 +0000