在有限的 16 字节字符串上将 IEEE 754-1985 双倍写入 ASCII

Writing IEEE 754-1985 double as ASCII on a limited 16 bytes string

这是我 的后续。但为了清楚起见,我会重复一遍:

根据 DICOM 标准,可以使用 Decimal 字符串的值表示来存储一种浮点数。见 :

Decimal String: A string of characters representing either a fixed point number or a floating point number. A fixed point number shall contain only the characters 0-9 with an optional leading "+" or "-" and an optional "." to mark the decimal point. A floating point number shall be conveyed as defined in ANSI X3.9, with an "E" or "e" to indicate the start of the exponent. Decimal Strings may be padded with leading or trailing spaces. Embedded spaces are not allowed.

"0"-"9", "+", "-", "E", "e", "." and the SPACE character of Default Character Repertoire. 16 bytes maximum

标准说文本表示是定点与浮点。该标准仅涉及值在 DICOM 数据集本身中的表示方式。因此,不需要将定点文本表示加载到定点变量中。

现在很明显 DICOM 标准暗示 citely 推荐 double (IEEE 754-1985) 来表示类型 Decimal StringValue Representation(最大16 位有效数字)。我的问题是如何使用标准 C I/O 库将此二进制表示从内存转换回 ASCII 到此有限大小的字符串?

来自互联网上的随机来源,这很重要,但通常 accepted solution 是:

printf("%1.16e\n", d); // Round-trippable double, always with an exponent

printf("%.17g\n", d); // Round-trippable double, shortest possible

当然,这两个表达式在我的例子中都是无效的,因为它们产生的输出比我限制的最大 16 字节 长得多。那么,当将任意双精度值写入有限的 16 字节字符串时,最小化 precision 中的损失的解决方案是什么?


编辑:如果这个不清楚,请我按照标准来。我不能使用 hex/uuencode 编码。

编辑2:我是运行使用travis-ci进行比较见:here

到目前为止,建议的代码是:

  1. Serge Ballesta
  2. chux
  3. Mark Dickinson
  4. chux

我在这里看到的结果是:

所以 compute4.c 导致最好的 precision (0.001287056579548422 < 4.050031792674619),但整体执行时间增加三倍 (x3)(仅在调试模式下使用 time 测试命令)。

C 库格式化程序没有满足您要求的直接格式。简单来说,如果你能接受标准%g格式的字符浪费e20写成e+020:浪费2个字符) , 你可以:

  • 生成 %.17g 格式的输出
  • 如果大于 16 个字符,计算导致 16 的精度
  • 生成该格式的输出。

代码可能如下所示:

void encode(double f, char *buf) {
    char line[40];
    char format[8];
    int prec;
    int l;

    l = sprintf(line, "%.17g", f);
    if (l > 16) {
        prec = 33 - strlen(line);
        l = sprintf(line, "%.*g", prec, f);
        while(l > 16) {
            /* putc('.', stdout);*/
            prec -=1;
            l = sprintf(line, "%.*g", prec, f);
        }
    }
    strcpy(buf, line);
}

如果你真的想做到最优(意思是写 e30 而不是 e+030),你可以尝试使用 %1.16e 格式并 post 处理输出。基本原理(正数):

  • %1.16e 格式允许您将尾数和指数(以 10 为底)分开
  • 如果指数在size-2(含)和size(不含)之间:只需将尾数正确舍入到int部分并显示即可
  • 如果指数介于 0 和 size-2 之间(均包含在内):显示圆角尾数并正确放置圆点
  • 如果指数介于 -1 和 -3 之间(均包含在内):以点开头,添加最终的 0 并填充圆角尾数
  • 否则使用指数部分最小的 e 格式并用四舍五入的尾数填充

极端案例:

  • 对于负数,放一个起始-,加上对数和size-1的显示
  • 舍入:如果第一个被拒绝的数字是 >=5,则增加前面的数字并迭代它是否是 9。处理 9.9999999999... 作为一个特殊情况四舍五入到 10

可能的代码:

void clean(char *mant) {
    char *ix = mant + strlen(mant) - 1;
    while(('0' == *ix) && (ix > mant)) {
        *ix-- = '[=11=]';
    }
    if ('.' == *ix) {
        *ix = '[=11=]';
    }
}

int add1(char *buf, int n) {
    if (n < 0) return 1;
    if (buf[n] == '9') {
        buf[n] = '0';
        return add1(buf, n-1);
    }
    else {
        buf[n] += 1;
    }
    return 0;
}

int doround(char *buf, unsigned int n) {
    char c;
    if (n >= strlen(buf)) return 0;
    c = buf[n];
    buf[n] = 0;
    if ((c >= '5') && (c <= '9')) return add1(buf, n-1);
    return 0;
}

int roundat(char *buf, unsigned int i, int iexp) {
    if (doround(buf, i) != 0) {
        iexp += 1;
        switch(iexp) {
            case -2:
                strcpy(buf, ".01");
                break;
            case -1:
                strcpy(buf, ".1");
                break;
            case 0:
                strcpy(buf, "1.");
                break;
            case 1:
                strcpy(buf, "10");
                break;
            case 2:
                strcpy(buf, "100");
                break;
            default:
                sprintf(buf, "1e%d", iexp);
        }
        return 1;
    }
    return 0;
}

void encode(double f, char *buf, int size) {
    char line[40];
    char *mant = line + 1;
    int iexp, lexp, i;
    char exp[6];

    if (f < 0) {
        f = -f;
        size -= 1;
        *buf++ = '-';
    }
    sprintf(line, "%1.16e", f);
    if (line[0] == '-') {
        f = -f;
    size -= 1;
    *buf++ = '-';
    sprintf(line, "%1.16e", f);
    }
    *mant = line[0];
    i = strcspn(mant, "eE");
    mant[i] = '[=11=]';
    iexp = strtol(mant + i + 1, NULL, 10);
    lexp = sprintf(exp, "e%d", iexp);
    if ((iexp >= size) || (iexp < -3)) {
        i = roundat(mant, size - 1 -lexp, iexp);
        if(i == 1) {
            strcpy(buf, mant);
            return;
        }
        buf[0] = mant[0];
        buf[1] = '.';
        strncpy(buf + i + 2, mant + 1, size - 2 - lexp);
        buf[size-lexp] = 0;
        clean(buf);
        strcat(buf, exp);
    }
    else if (iexp >= size - 2) {
        roundat(mant, iexp + 1, iexp);
        strcpy(buf, mant);
    }
    else if (iexp >= 0) {
        i = roundat(mant, size - 1, iexp);
        if (i == 1) {
            strcpy(buf, mant);
            return;
        }
        strncpy(buf, mant, iexp + 1);
        buf[iexp + 1] = '.';
        strncpy(buf + iexp + 2, mant + iexp + 1, size - iexp - 1);
        buf[size] = 0;
        clean(buf);
    }
    else {
        int j;
        i = roundat(mant, size + 1 + iexp, iexp);
        if (i == 1) {
            strcpy(buf, mant);
            return;
        }
        buf[0] = '.';
        for(j=0; j< -1 - iexp; j++) {
            buf[j+1] = '0';
        }
        if ((i == 1) && (iexp != -1)) {
            buf[-iexp] = '1';
            buf++;
        }
        strncpy(buf - iexp, mant, size + 1 + iexp);
        buf[size] = 0;
        clean(buf);
    }
}

对于有限浮点值,printf() 格式说明符 "%e" 非常匹配
"A floating point number shall be ... with an "E" 或 "e" 表示指数的开始

[−]d.ddd...ddde±dd

该符号以负数存在,可能 -0.0。指数至少为2位。

如果我们假设 DBL_MAX < 1e1000,(对于IEEE 754-1985 double是安全的),那么以下在所有情况下都有效:1个可选符号,1个前导数字, '.',8位,'e',符号,最多3位。

(注意:“最大16字节”似乎不是指C字符串空字符终止。如果需要调整1。)

// Room for 16 printable characters.
char buf[16+1];
int n = snprintf(buf, sizeof buf, "%.*e", 8, x);
assert(n >= 0 && n < sizeof buf);
puts(buf);

但这为可选符号和 2 到 3 个指数数字保留了空间。

诀窍是由于四舍五入,当一个数字使用 2 或使用 3 指数数字时的边界是模糊的。即使测试负数,-0.0 也是一个问题。

[编辑] 还需要对非常小的数字进行测试。

候选人:

// Room for 16 printable characters.
char buf[16+1];
assert(isfinite(x)); // for now, only address finite numbers

int precision = 8+1+1;
if (signbit(x)) precision--;  // Or simply `if (x <= 0.0) precision--;`
if (fabs(x) >= 9.99999999e99) precision--; // some refinement possible here.
else if (fabs(x) <= 1.0e-99) precision--;
int n = snprintf(buf, sizeof buf, "%.*e", precision, x);
assert(n >= 0 && n < sizeof buf);
puts(buf);

其他问题:

一些编译器至少打印 3 个指数数字。
IEEE 754-1985 double needed 的最大十进制有效数字位数因需求定义而异,但可能约为 15-17。 Printf width specifier to maintain precision of floating-point value

候选人2:一次测试输出太长

// Room for N printable characters.
#define N 16
char buf[N+1];
assert(isfinite(x)); // for now, only address finite numbers

int precision = N - 2 - 4;  // 1.xxxxxxxxxxe-dd
if (signbit(x)) precision--;
int n = snprintf(buf, sizeof buf, "%.*e", precision, x);
if (n >= sizeof buf) {
  n = snprintf(buf, sizeof buf, "%.*e", precision - (n - sizeof buf) - 1, x);
}
assert(n >= 0 && n < sizeof buf);
puts(buf);

无法以十进制打印,因为对于某些数字,需要 17 位尾数,这会在不打印指数的情况下用尽您的所有 space。更准确地说,打印十进制的双精度有时需要超过 16 个字符才能保证准确的往返。

相反,您应该使用十六进制打印底层二进制表示。假设不需要空终止符,这将正好使用 16 个字节。

如果你想使用少于 16 个字节打印结果,那么你基本上可以对它进行 uuencode。也就是说,使用超过 16 位的数字,以便您可以将更多的位压缩到每个数字中。如果您使用 64 个不同的字符(6 位),则可以用 11 个字符打印 64 位双精度数。可读性不是很好,但必须做出权衡。

我认为你最好的选择是使用 printf("%.17g\n", d);生成初始答案,然后 trim 它。 trim 的最简单方法是从尾数末尾删除数字,直到适合为止。这实际上工作得很好,但不会最小化错误,因为您正在截断而不是四舍五入到最近。

更好的解决方案是检查要删除的数字,将它们视为 0.0 到 1.0 之间的 n 位数字,因此“49”将是 0.49。如果它们的值小于 0.5,则删除它们。如果它们的值大于 0.50,则以十进制形式递增打印值。也就是说,将 1 添加到最后一位,并根据需要进行环绕和进位。创建的任何尾随零都应 trimmed.

这成为问题的唯一时间是进位一直传播到第一个数字并将其从 9 溢出到零。这可能是不可能的,但我不确定。在这种情况下 (+9.99999e17) 答案将是 +1e18,因此只要您对这种情况进行了测试就没问题。

因此,打印数字,将其拆分为 sign/mantissa 个字符串和一个指数整数,然后对它们进行字符串操作以获得结果。

这比最初想象的要棘手。

考虑到各种极端情况,最好尝试高精度,然后根据需要进行处理。

  1. 由于 '-'.

  2. ,任何负数都与正数打印相同,但精度降低 1
  3. '+' 符号不需要在字符串的开头,也不需要在 'e'.

  4. 之后
  5. '.' 不需要。

  6. 使用 sprintf() 以外的任何东西来完成给定的数学部分是危险的,所以 许多 极端情况。给定各种舍入模式,FLT_EVAL_METHOD 等,将繁重的编码留给完善的函数。

  7. 当一次尝试的长度超过 1 个字符时,可以保存迭代。例如。如果尝试使用精度 14,结果宽度为 20,则无需尝试精度 13 和 12,只需转到 11。

  8. 由于移除 '.' 而对指数进行缩放,必须在 sprintf() 之后完成,以便 1) 避免注入计算错误 2) 递减 double 低于其最小指数。

  9. -1.00000000049999e-200 一样,最大相对误差小于 2,000,000,000 分之一。平均相对误差约为 50,000,000,000 分之一。

  10. 14 位精度,最高,出现在像 12345678901234e1 这样的数字上,所以从 16-2 位开始。


static size_t shrink(char *fp_buffer) {
  int lead, expo;
  long long mant;
  int n0, n1;
  int n = sscanf(fp_buffer, "%d.%n%lld%ne%d", &lead, &n0, &mant, &n1, &expo);
  assert(n == 3);
  return sprintf(fp_buffer, "%d%0*llde%d", lead, n1 - n0, mant,
          expo - (n1 - n0));
}

int x16printf(char *dest, size_t width, double value) {
  if (!isfinite(value)) return 1;

  if (width < 5) return 2;
  if (signbit(value)) {
    value = -value;
    strcpy(dest++, "-");
    width--;
  }
  int precision = width - 2;
  while (precision > 0) {
    char buffer[width + 10];
    // %.*e prints 1 digit, '.' and then `precision - 1` digits
    snprintf(buffer, sizeof buffer, "%.*e", precision - 1, value);
    size_t n = shrink(buffer);
    if (n <= width) {
      strcpy(dest, buffer);
      return 0;
    }
    if (n > width + 1) precision -= n - width - 1;
    else precision--;
  }
  return 3;
}

测试代码

double rand_double(void) {
  union {
    double d;
    unsigned char uc[sizeof(double)];
  } u;
  do {
    for (size_t i = 0; i < sizeof(double); i++) {
      u.uc[i] = rand();
    }
  } while (!isfinite(u.d));
  return u.d;
}

void x16printf_test(double value) {
  printf("%-27.*e", 17, value);
  char buf[16+1];
  buf[0] = 0;
  int y = x16printf(buf, sizeof buf - 1, value);
  printf(" %d\n", y);
  printf("'%s'\n", buf);
}


int main(void) {
  for (int i = 0; i < 10; i++)
    x16printf_test(rand_double());
}

输出

-1.55736829786841915e+118   0
'-15573682979e108'
-3.06117209691283956e+125   0
'-30611720969e115'
8.05005611774356367e+175    0
'805005611774e164'
-1.06083057094522472e+132   0
'-10608305709e122'
3.39265065244054607e-209    0
'33926506524e-219'
-2.36818580315246204e-244   0
'-2368185803e-253'
7.91188576978592497e+301    0
'791188576979e290'
-1.40513111051994779e-53    0
'-14051311105e-63'
-1.37897140950449389e-14    0
'-13789714095e-24'
-2.15869805640288206e+125   0
'-21586980564e115'