📄 lib.txt
字号:
02341
02342 /* add lower pair of 32 bits */
02343 carry = ((unsigned long) 0xFFFFFFFF - e1->l_32 < e2->l_32);
02344 e1->l_32 += e2->l_32;
02345 if ((carry) && (++e1->h_32 == 0))
02346 return(1); /* had a 64 bit overflow */
02347 else
02348 return(overflow); /* return status from higher add */
02349 }
02351 /* The following tables can be computed with the following bc(1)
02352 program:
02353
02354 obase=16
02355 scale=0
02356 define t(x){
02357 auto a, b, c
02358 a=2;b=1;c=2^32;n=1
02359 while(a<x) {
02360 b=a;n+=n;a*=a
02361 }
02362 n/=2
02363 a=b
02364 while(b<x) {
02365 a=b;b*=c;n+=32
02366 }
02367 n-=32
02368 b=a
02369 while(a<x) {
02370 b=a;a+=a;n+=1
02371 }
02372 n-=1
02373 x*=16^16
02374 b=x%a
02375 x/=a
02376 if(a<=(2*b)) x+=1
02377 obase=10
02378 n
02379 obase=16
02380 return(x)
02381 }
02382 for (i=1;i<28;i++) {
02383 t(10^i)
02384 }
02385 0
02386 for (i=1;i<20;i++) {
02387 t(10^(28*i))
02388 }
02389 0
02390 define r(x){
02391 auto a, b, c
02392 a=2;b=1;c=2^32;n=1
02393 while(a<x) {
02394 b=a;n+=n;a*=a
02395 }
02396 n/=2
02397 a=b
02398 while(b<x) {
02399 a=b;b*=c;n+=32
02400 }
02401 n-=32
02402 b=a
02403 while(a<x) {
02404 b=a;a+=a;n+=1
02405 }
02406 a=b
02407 a*=16^16
02408 b=a%x
02409 a/=x
02410 if(x<=(2*b)) a+=1
02411 obase=10
02412 -n
02413 obase=16
02414 return(a)
02415 }
02416 for (i=1;i<28;i++) {
02417 r(10^i)
02418 }
02419 0
02420 for (i=1;i<20;i++) {
02421 r(10^(28*i))
02422 }
02423 0
02424
02425 */
02426 static struct EXTEND ten_powers[] = { /* representation of 10 ** i */
02427 { 0, 0, 0x80000000, 0 },
02428 { 0, 3, 0xA0000000, 0 },
02429 { 0, 6, 0xC8000000, 0 },
02430 { 0, 9, 0xFA000000, 0 },
02431 { 0, 13, 0x9C400000, 0 },
02432 { 0, 16, 0xC3500000, 0 },
02433 { 0, 19, 0xF4240000, 0 },
02434 { 0, 23, 0x98968000, 0 },
02435 { 0, 26, 0xBEBC2000, 0 },
02436 { 0, 29, 0xEE6B2800, 0 },
02437 { 0, 33, 0x9502F900, 0 },
02438 { 0, 36, 0xBA43B740, 0 },
02439 { 0, 39, 0xE8D4A510, 0 },
02440 { 0, 43, 0x9184E72A, 0 },
02441 { 0, 46, 0xB5E620F4, 0x80000000 },
02442 { 0, 49, 0xE35FA931, 0xA0000000 },
02443 { 0, 53, 0x8E1BC9BF, 0x04000000 },
02444 { 0, 56, 0xB1A2BC2E, 0xC5000000 },
02445 { 0, 59, 0xDE0B6B3A, 0x76400000 },
02446 { 0, 63, 0x8AC72304, 0x89E80000 },
02447 { 0, 66, 0xAD78EBC5, 0xAC620000 },
02448 { 0, 69, 0xD8D726B7, 0x177A8000 },
02449 { 0, 73, 0x87867832, 0x6EAC9000 },
02450 { 0, 76, 0xA968163F, 0x0A57B400 },
02451 { 0, 79, 0xD3C21BCE, 0xCCEDA100 },
02452 { 0, 83, 0x84595161, 0x401484A0 },
02453 { 0, 86, 0xA56FA5B9, 0x9019A5C8 },
02454 { 0, 89, 0xCECB8F27, 0xF4200F3A }
02455 };
02456 static struct EXTEND big_ten_powers[] = { /* representation of 10 ** (28*i) */
02457 { 0, 0, 0x80000000, 0 },
02458 { 0, 93, 0x813F3978, 0xF8940984 },
02459 { 0, 186, 0x82818F12, 0x81ED44A0 },
02460 { 0, 279, 0x83C7088E, 0x1AAB65DB },
02461 { 0, 372, 0x850FADC0, 0x9923329E },
02462 { 0, 465, 0x865B8692, 0x5B9BC5C2 },
02463 { 0, 558, 0x87AA9AFF, 0x79042287 },
02464 { 0, 651, 0x88FCF317, 0xF22241E2 },
02465 { 0, 744, 0x8A5296FF, 0xE33CC930 },
02466 { 0, 837, 0x8BAB8EEF, 0xB6409C1A },
02467 { 0, 930, 0x8D07E334, 0x55637EB3 },
02468 { 0, 1023, 0x8E679C2F, 0x5E44FF8F },
02469 { 0, 1116, 0x8FCAC257, 0x558EE4E6 },
02470 { 0, 1209, 0x91315E37, 0xDB165AA9 },
02471 { 0, 1302, 0x929B7871, 0xDE7F22B9 },
02472 { 0, 1395, 0x940919BB, 0xD4620B6D },
02473 { 0, 1488, 0x957A4AE1, 0xEBF7F3D4 },
02474 { 0, 1581, 0x96EF14C6, 0x454AA840 },
02475 { 0, 1674, 0x98678061, 0x27ECE4F5 },
02476 { 0, 1767, 0x99E396C1, 0x3A3ACFF2 }
02477 };
02478
02479 static struct EXTEND r_ten_powers[] = { /* representation of 10 ** -i */
02480 { 0, 0, 0x80000000, 0 },
02481 { 0, -4, 0xCCCCCCCC, 0xCCCCCCCD },
02482 { 0, -7, 0xA3D70A3D, 0x70A3D70A },
02483 { 0, -10, 0x83126E97, 0x8D4FDF3B },
02484 { 0, -14, 0xD1B71758, 0xE219652C },
02485 { 0, -17, 0xA7C5AC47, 0x1B478423 },
02486 { 0, -20, 0x8637BD05, 0xAF6C69B6 },
02487 { 0, -24, 0xD6BF94D5, 0xE57A42BC },
02488 { 0, -27, 0xABCC7711, 0x8461CEFD },
02489 { 0, -30, 0x89705F41, 0x36B4A597 },
02490 { 0, -34, 0xDBE6FECE, 0xBDEDD5BF },
02491 { 0, -37, 0xAFEBFF0B, 0xCB24AAFF },
02492 { 0, -40, 0x8CBCCC09, 0x6F5088CC },
02493 { 0, -44, 0xE12E1342, 0x4BB40E13 },
02494 { 0, -47, 0xB424DC35, 0x095CD80F },
02495 { 0, -50, 0x901D7CF7, 0x3AB0ACD9 },
02496 { 0, -54, 0xE69594BE, 0xC44DE15B },
02497 { 0, -57, 0xB877AA32, 0x36A4B449 },
02498 { 0, -60, 0x9392EE8E, 0x921D5D07 },
02499 { 0, -64, 0xEC1E4A7D, 0xB69561A5 },
02500 { 0, -67, 0xBCE50864, 0x92111AEB },
02501 { 0, -70, 0x971DA050, 0x74DA7BEF },
02502 { 0, -74, 0xF1C90080, 0xBAF72CB1 },
02503 { 0, -77, 0xC16D9A00, 0x95928A27 },
02504 { 0, -80, 0x9ABE14CD, 0x44753B53 },
02505 { 0, -84, 0xF79687AE, 0xD3EEC551 },
02506 { 0, -87, 0xC6120625, 0x76589DDB },
02507 { 0, -90, 0x9E74D1B7, 0x91E07E48 }
02508 };
02509
02510 static struct EXTEND r_big_ten_powers[] = { /* representation of 10 ** -(28*i) */
02511 { 0, 0, 0x80000000, 0 },
02512 { 0, -94, 0xFD87B5F2, 0x8300CA0E },
02513 { 0, -187, 0xFB158592, 0xBE068D2F },
02514 { 0, -280, 0xF8A95FCF, 0x88747D94 },
02515 { 0, -373, 0xF64335BC, 0xF065D37D },
02516 { 0, -466, 0xF3E2F893, 0xDEC3F126 },
02517 { 0, -559, 0xF18899B1, 0xBC3F8CA2 },
02518 { 0, -652, 0xEF340A98, 0x172AACE5 },
02519 { 0, -745, 0xECE53CEC, 0x4A314EBE },
02520 { 0, -838, 0xEA9C2277, 0x23EE8BCB },
02521 { 0, -931, 0xE858AD24, 0x8F5C22CA },
02522 { 0, -1024, 0xE61ACF03, 0x3D1A45DF },
02523 { 0, -1117, 0xE3E27A44, 0x4D8D98B8 },
02524 { 0, -1210, 0xE1AFA13A, 0xFBD14D6E },
02525 { 0, -1303, 0xDF82365C, 0x497B5454 },
02526 { 0, -1396, 0xDD5A2C3E, 0xAB3097CC },
02527 { 0, -1489, 0xDB377599, 0xB6074245 },
02528 { 0, -1582, 0xD91A0545, 0xCDB51186 },
02529 { 0, -1675, 0xD701CE3B, 0xD387BF48 },
02530 { 0, -1768, 0xD4EEC394, 0xD6258BF8 }
02531 };
02532
02533 #define TP (int)(sizeof(ten_powers)/sizeof(ten_powers[0]))
02534 #define BTP (int)(sizeof(big_ten_powers)/sizeof(big_ten_powers[0]))
02535 #define MAX_EXP (TP * BTP - 1)
02536
02537 static
02538 add_exponent(struct EXTEND *e, int exp)
02539 {
02540 int neg = exp < 0;
02541 int divsz, modsz;
02542 struct EXTEND x;
02543
02544 if (neg) exp = -exp;
02545 divsz = exp / TP;
02546 modsz = exp % TP;
02547 if (neg) {
02548 mul_ext(e, &r_ten_powers[modsz], &x);
02549 mul_ext(&x, &r_big_ten_powers[divsz], e);
02550 }
02551 else {
02552 mul_ext(e, &ten_powers[modsz], &x);
02553 mul_ext(&x, &big_ten_powers[divsz], e);
02554 }
02555 }
02557 _str_ext_cvt(const char *s, char **ss, struct EXTEND *e)
02558 {
02559 /* Like strtod, but for extended precision */
02560 register int c;
02561 int dotseen = 0;
02562 int digitseen = 0;
02563 int exp = 0;
02564
02565 if (ss) *ss = (char *)s;
02566 while (isspace(*s)) s++;
02567
02568 e->sign = 0;
02569 e->exp = 0;
02570 e->m1 = e->m2 = 0;
02571
02572 c = *s;
02573 switch(c) {
02574 case '-':
02575 e->sign = 1;
02576 case '+':
02577 s++;
02578 }
02579 while (c = *s++, isdigit(c) || (c == '.' && ! dotseen++)) {
02580 if (c == '.') continue;
02581 digitseen = 1;
02582 if (e->m1 <= (unsigned long)(0xFFFFFFFF)/10) {
02583 struct mantissa a1;
02584
02585 a1 = e->mantissa;
02586 b64_sft(&(e->mantissa), -3);
02587 b64_sft(&a1, -1);
02588 b64_add(&(e->mantissa), &a1);
02589 a1.h_32 = 0;
02590 a1.l_32 = c - '0';
02591 b64_add(&(e->mantissa), &a1);
02592 }
02593 else exp++;
02594 if (dotseen) exp--;
02595 }
02596 if (! digitseen) return;
02597
02598 if (ss) *ss = (char *)s - 1;
02599
02600 if (c == 'E' || c == 'e') {
02601 int exp1 = 0;
02602 int sign = 1;
02603 int exp_overflow = 0;
02604
02605 switch(*s) {
02606 case '-':
02607 sign = -1;
02608 case '+':
02609 s++;
02610 }
02611 if (c = *s, isdigit(c)) {
02612 do {
02613 int tmp;
02614
02615 exp1 = 10 * exp1 + (c - '0');
02616 if ((tmp = sign * exp1 + exp) > MAX_EXP ||
02617 tmp < -MAX_EXP) {
02618 exp_overflow = 1;
02619 }
02620 } while (c = *++s, isdigit(c));
02621 if (ss) *ss = (char *)s;
02622 }
02623 exp += sign * exp1;
02624 if (exp_overflow) {
02625 exp = sign * MAX_EXP;
02626 if (e->m1 != 0 || e->m2 != 0) errno = ERANGE;
02627 }
02628 }
02629 if (e->m1 == 0 && e->m2 == 0) return;
02630 e->exp = 63;
02631 while (! (e->m1 & 0x80000000)) {
02632 b64_sft(&(e->mantissa),-1);
02633 e->exp--;
02634 }
02635 add_exponent(e, exp);
02636 }
02638 #include <math.h>
02639
02640 static
02641 ten_mult(struct EXTEND *e)
02642 {
02643 struct EXTEND e1 = *e;
02644
02645 e1.exp++;
02646 e->exp += 3;
02647 add_ext(e, &e1, e);
02648 }
02650 #define NDIGITS 128
02651 #define NSIGNIFICANT 19
02652
02653 char *
02654 _ext_str_cvt(struct EXTEND *e, int ndigit, int *decpt, int *sign, int ecvtflag)
02655 {
02656 /* Like cvt(), but for extended precision */
02657
02658 static char buf[NDIGITS+1];
02659 struct EXTEND m;
02660 register char *p = buf;
02661 register char *pe;
02662 int findex = 0;
02663
02664 if (ndigit < 0) ndigit = 0;
02665 if (ndigit > NDIGITS) ndigit = NDIGITS;
02666 pe = &buf[ndigit];
02667 buf[0] = '\0';
02668
02669 *sign = 0;
02670 if (e->sign) {
02671 *sign = 1;
02672 e->sign = 0;
02673 }
02674
02675 *decpt = 0;
02676 if (e->m1 != 0) {
02677 register struct EXTEND *pp = &big_ten_powers[1];
02678
02679 while(cmp_ext(e,pp) >= 0) {
02680 pp++;
02681 findex = pp - big_ten_powers;
02682 if (findex >= BTP) break;
02683 }
02684 pp--;
02685 findex = pp - big_ten_powers;
02686 mul_ext(e,&r_big_ten_powers[findex],e);
02687 *decpt += findex * TP;
02688 pp = &ten_powers[1];
02689 while(pp < &ten_powers[TP] && cmp_ext(e, pp) >= 0) pp++;
02690 pp--;
02691 findex = pp - ten_powers;
02692 *decpt += findex;
02693
02694 if (cmp_ext(e, &ten_powers[0]) < 0) {
02695 pp = &r_big_ten_powers[1];
02696 while(cmp_ext(e,pp) < 0) pp++;
02697 pp--;
02698 findex = pp - r_big_ten_powers;
02699 mul_ext(e, &big_ten_powers[findex], e);
02700 *decpt -= findex * TP;
02701 /* here, value >= 10 ** -28 */
02702 ten_mult(e);
02703 (*decpt)--;
02704 pp = &r_ten_powers[0];
02705 while(cmp_ext(e, pp) < 0) pp++;
02706 findex = pp - r_ten_powers;
02707 mul_ext(e, &ten_powers[findex], e);
02708 *decpt -= findex;
02709 findex = 0;
02710 }
02711 (*decpt)++; /* because now value in [1.0, 10.0) */
02712 }
02713 if (! ecvtflag) {
02714 /* for fcvt() we need ndigit digits behind the dot */
02715 pe += *decpt;
02716 if (pe > &buf[NDIGITS]) pe = &buf[NDIGITS];
02717 }
02718 m.exp = -62;
02719 m.sign = 0;
02720 m.m1 = 0xA0000000;
02721 m.m2 = 0;
02722 while (p <= pe) {
02723 struct EXTEND oneminm;
02724
02725 if (p - pe > NSIGNIFICANT) {
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -