将任意精度有理数（OCaml，zarith）转换为近似浮点数-Java 学习之路

我正在使用Zarith库来执行任意精度的有理算术 . 假设我有一个 Q.t 类型的有理数 q ，它是两个大整数的比率（ Q 是Zarith的任意精度有理数模块） . 有时，为了便于阅读，我想将此数字打印为浮点数，有时我需要将此数字浮点数转换为以后的非任意精度计算 . 有没有办法将 q 转换为浮点数达到一定的精度？

我将 q 转换为浮点的方式现在无法保证，并且可以创建未定义的浮点数（ Z 是任意精度整数模块）：

let to_float q =
  let n, d = num q, den q in
  (* check if d is zero and raise an error if it is *)
  let nf, df = Z.to_float n, Z.to_float d in
  nf /. df

有没有更好的方法来处理这个问题，我可以获得最准确接近任何 q 的浮点数？

Edit

如果有人有兴趣，我很快就会在OCaml中写下Mark Dickinson的回答 . 它可能（绝对）可以改进和清理 . 如果我这样做，或者如果有人对改进有任何建议，我会编辑 . 但是现在这已经解决了我的问题！

let to_float q = 
  let n, d = num q, den q in
  let n_sign = Z.sign n in
  let d_sign = Z.sign d in (* always >= 0 *)
  if d_sign = 0 then raise Division_by_zero;
  let n = Z.abs n in
  if n_sign = 0 then 0. else
    let shift = (Z.numbits n) - (Z.numbits d) - 55 in
    let is_subnormal = shift < -1076 in
    let shift = if is_subnormal then -1076 else shift in
    let d = if shift >= 0 then Z.shift_left d shift else d in
    let n = if shift < 0 then Z.shift_left n (-shift)
      else n in
    let quotient, remainder = Z.div_rem n d in
    let quotient = if (Z.compare remainder (Z.zero)) = 0 && Z.is_even quotient then
        Z.add Z.one quotient else quotient in
    let quotient = if not is_subnormal then quotient else
        let round_select = Z.to_int @@ Z.rem quotient @@ Z.of_int 8 in
        Z.add quotient [|Z.zero;Z.minus_one;Z.of_int (-2);Z.one;Z.zero
                        ;Z.minus_one;Z.of_int 2;Z.one|].(round_select)
    in
    let unsigned_res = ldexp (Z.to_float quotient) shift in                                                                                                             
    if n_sign = 1 then unsigned_res else -.unsigned_res

我稍后会有'll look into writing an interface for GMP' s mpq_get_d功能，但我不完全确定如何做到这一点 . 我看到如何做到这一点的唯一方法是将 q : Q.t 转换为字符串并传递to:

int mpq_set_str (mpq_t rop, const char *str, int base)

有谁知道如何在OCaml中将 rop 传递给 mpq_get_d 或者有一个描述如何执行此操作的参考？我看了chapter 19 of RWO并没有看到这样的情况 .

2 回答

如果您有权访问

整数 log2 操作，和
能够将整数左移一个给定的位数

然后滚动自己正确的舍入转换相对容易 . 简而言之，该方法如下所示：

简化为案例 n > 0 ， d > 0 ;过滤掉明显的下溢/溢出
选择整数 shift ，以便 2^-shift*n/d 位于 2^54 和 2^56 之间 .
使用整数算术计算 x = 2^-shift*n/d ，使用round-to-odd舍入方法舍入为最接近的整数 .
将 x 转换为最接近的IEEE 754双精度值 dx ，采用通常的圆形连接到均匀舍入模式 .
返回 ldexp(dx, shift) .

我不熟悉OCaml，但下面的Python代码说明了积极输入的想法 . 我留给你对负输入进行明显的修改并将其除以零 . 您可能还希望尽早返回极端溢出和下溢的情况：通过查找下面的 shift 的额外大或小值，可以很容易地检测到这些情况 .

from math import ldexp

def to_float(numerator, denominator):
    """
    Convert numerator / denominator to float, correctly rounded.

    For simplicity, assume both inputs are positive.
    """
    # Shift satisfies 2**54 < (numerator / denominator) / 2**shift < 2**56
    shift = numerator.bit_length() - denominator.bit_length() - 55

    # Divide the fraction by 2**shift.
    if shift >= 0:
        denominator <<= shift
    else:
        numerator <<= -shift

    # Convert to the nearest integer, using round-to-odd.
    q, r = divmod(numerator, denominator)
    if r != 0 and q % 2 == 0:
        q += 1

    # Now convert to the nearest float and shift back.
    return ldexp(float(q), shift)

一些说明：

正整数 n 上的 bit_length 方法给出了表示 n 所需的位数，换句话说 1 + floor(log2(n)) .
divmod 是一个Python函数，它同时计算整数除法的商和余数 .
数量 q 适合（轻松）在64位整数
我们四舍五入：一次将移位的 numerator / denominator 转换为最接近的整数，再次将该整数四舍五入为浮点数 . 第一轮使用round-to-odd方法;这确保了第二轮（从int到float的转换中隐含）给出了相同的结果，就像我们将分数直接舍入到float一样 .
上述算法无法正确处理转换后的浮点值为低于正常的分数：在这种情况下， ldexp 操作可能会引入第三个舍入 . 可以小心处理这个问题 . 请参阅下面的一些代码 .

以上实际上是Python在将一个（大）整数除以另一个（大）整数以获得浮点结果时使用的算法的简化版本 . 你可以看到源here . long_true_divide 函数开头的注释概述了该方法 .

为了完整性，这里的变体也可以正确处理次正常结果 .

def to_float(numerator, denominator):
    """
    Convert numerator / denominator to float, correctly rounded.

    For simplicity, assume both inputs are positive.
    """
    # Choose shift so that 2**54 < numerator / denominator / 2**shift < 2**56
    shift = numerator.bit_length() - denominator.bit_length() - 55

    # The 'treat_as_subnormal' flag catches all cases of subnormal results,
    # along with some cases where the result is not subnormal but *is* still
    # smaller than 2**-1021. In all these cases, it's sufficient to find the
    # closest integer multiple of 2**-1074. We first round to the nearest
    # multiple of 2**-1076 using round-to-odd.
    treat_as_subnormal = shift < -1076
    if treat_as_subnormal:
        shift = -1076

    # Divide the fraction by 2**shift.
    if shift >= 0:
        denominator <<= shift
    else:
        numerator <<= -shift

    # Convert to the nearest integer, using round-to-odd.
    q, r = divmod(numerator, denominator)
    if r != 0 and q % 2 == 0:
        q += 1

    # Now convert to the nearest float and shift back.
    if treat_as_subnormal:
        # Round to the nearest multiple of 4, rounding ties to
        # the nearest multiple of 8. This avoids double rounding
        # from the ldexp call below.
        q += [0, -1, -2, 1, 0, -1, 2, 1][q%8]

    return ldexp(float(q), shift)

回复于 2024-04-30T00:01:46+08:00

2

这不是一个完整的答案，但环顾四周，我发现Zarith在内部使用GMP . 有一个名为 mpq_get_d 的GMP函数将理性转换为double . 如果它不能直接在Zarith中使用，那么应该（给定时间）为它添加一个接口 .

回复于 2024-04-30T00:01:46+08:00

将任意精度有理数（OCaml，zarith）转换为近似浮点数

2 回答

相关问题