GMPを使ってDBL_MAXの文字列をdoubleに変換するルーチン

/*
 
  dbl_max_gmp.c
 
  GMPを用いて、DBL_MAXの値(10進数文字列)をdoubleに変換する。
 
  コンパイル:
  % cl /Ic:\include dbl_max_gmp.c gmp.lib /link /LIBPATH:c:\lib
 
  実行例:
  % ./dbl_max_gmp
     mpf = 0.17976931348623158e309
  double =  1.79769313486231570e+308
 
 */
 
#include <stdio.h>
#include <float.h>
#include <string.h>
#include "gmp.h"
 
#define LOG_TEN_TWO  3.32192809488736234789
#define bprec(n) (int)(((n+10)*LOG_TEN_TWO)+2)
 
void atompf(char s[], mpf_t mp);
 
main()
{
    mpf_t f;
    long prec, dprec;
 
    dprec = 1000L;  /* decimal precision */
    prec = bprec(dprec);    /* binary precision (plus alpha) */
    mpf_set_default_prec(prec);
 
    mpf_init(f);
    atompf("1.7976931348623158e+308", f);   /* DBL_MAX in Visual C++ */
 
    printf("   mpf = ");
    mpf_out_str(stdout, 10, dprec + 10, f);
    putchar('\n');
    printf("double = %25.17e\n", mpf_get_d(f));
 
    return 0;
}
 
/* atompf:  convert a string into multi-precision floating-point */
void atompf(char s[], mpf_t f)
{
    int i, c, exp;
    int frc_sign, frc_exp, exp_sign, exp_exp;
    mpz_t mp_frc, mp_var, mp_exp;
    mpf_t mp_fp1, mp_fp2, mp_fp3;
 
    /* init */
    mpz_init(mp_frc);
    mpz_init(mp_var);
    mpz_init(mp_exp);
    mpf_init(mp_fp1);
    mpf_init(mp_fp2);
    mpf_init(mp_fp3);
 
    /* skip leading blanks */
    for (i=0; isspace(s[i]); i++)
        ;
 
    /* sign */
    frc_sign = 1;
    if ((c = s[i]) == '-')
        frc_sign = -1;
    if (c == '-' || c == '+')
        i++;
 
    /* integral part */
    mpz_set_ui(mp_frc, 0);
    for (; isdigit(s[i]); i++) {
        mpz_mul_ui(mp_var, mp_frc, 10);
        mpz_add_ui(mp_frc, mp_var, s[i] - '0');
    }
 
    /* fractional part */
    if (s[i] == '.') {
        frc_exp = 0;
        for (i++; isdigit(s[i]); i++) {
            mpz_mul_ui(mp_var, mp_frc, 10);
            mpz_add_ui(mp_frc, mp_var, s[i] - '0');
            frc_exp++;
        }
    }
 
    /* exponential part */
    if ((c = toupper(s[i])) == 'E' || c == 'D')
        i++;
    exp_sign = 1;
    if ((c = s[i]) == '-')
        exp_sign = -1;
    if (c == '+' || c == '-')
        i++;
    exp_exp = 0;
    for (; isdigit(s[i]); i++) {
        exp_exp = 10 * exp_exp + (s[i] - '0');
    }
 
    /* add exponent and store its power in mp_exp */
    exp = exp_sign * exp_exp - frc_exp;
    if (exp > 0)
        exp_sign = 1;
    else {
        exp_sign = -1;
        exp = -exp;
    }
    mpz_set_ui(mp_exp, 1);
    for (; exp > 0; exp--) {
        mpz_mul_ui(mp_var, mp_exp, 10);
        mpz_swap(mp_exp, mp_var);
    }
 
    /* build a mpf_t value */
    mpf_set_z(mp_fp1, mp_frc);
    mpf_set_z(mp_fp2, mp_exp);
    if (exp_sign > 0) {
        mpf_mul(mp_fp3, mp_fp1, mp_fp2);
    } else {
        mpf_div(mp_fp3, mp_fp1, mp_fp2);
    }
 
    /* copy */
    mpf_set(f, mp_fp3);
 
    /* clear */
    mpz_clear(mp_frc);
    mpz_clear(mp_exp);
    mpz_clear(mp_var);
    mpf_clear(mp_fp1);
    mpf_clear(mp_fp2);
    mpf_clear(mp_fp3);
}