将任意精度有理数(OCaml、zarith)转换为近似浮点数
Converting an arbitrary-precision rational number (OCaml, zarith) to an approximate floating number
我正在使用 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
的浮点数?
编辑
如果有人感兴趣的话,我很快就用 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
我会考虑为 GMP 编写一个界面 mpq_get_d
function later, but I'm not entirely sure how to do that. The only way I see how to do that is to convert q : Q.t
to a string and pass that to:
int mpq_set_str (mpq_t rop, const char *str, int base)
有人知道如何在 OCaml 中将 rop
传递给 mpq_get_d
或者有描述如何执行此操作的参考吗?我翻了一下chapter 19 of RWO,没有看到这样的情况。
这不是一个完整的答案,但环顾四周,我发现 Zarith 在内部使用 GMP。有一个名为 mpq_get_d
的 GMP 函数可将有理数转换为双精度数。如果它不能直接在 Zarith 中使用,应该可以(给一些时间)为它添加一个接口。
如果您有权访问
- 整数
log2
运算,
- 将整数左移给定位数的能力
然后进行您自己的正确舍入转换相对容易。简而言之,该方法如下所示:
- 减少到
n > 0
、d > 0
的大小写;过滤掉明显的 underflow/overflow
- 选择一个整数
shift
,使 2^-shift*n/d
介于 2^54
和 2^56
之间。
- 使用整数算法计算
x = 2^-shift*n/d
,使用 round-to-odd 舍入方法舍入到最接近的整数。
- 将
x
转换为最接近的 IEEE 754 双精度值 dx
,使用通常的舍入到偶数舍入模式。
- Return
ldexp(dx, shift)
.
恐怕我对 OCaml 不是很流利,但下面的 Python 代码说明了正输入的想法。我留给你对负输入和除以零进行明显的修改。您可能还想为极端上溢和下溢的情况尽早设置 return:通过在下方查找 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
运算可能会引入 third 舍入。可以小心处理这个问题。请参阅下面的一些代码。
以上实际上是 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)
我正在使用 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
的浮点数?
编辑
如果有人感兴趣的话,我很快就用 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
我会考虑为 GMP 编写一个界面 mpq_get_d
function later, but I'm not entirely sure how to do that. The only way I see how to do that is to convert q : Q.t
to a string and pass that to:
int mpq_set_str (mpq_t rop, const char *str, int base)
有人知道如何在 OCaml 中将 rop
传递给 mpq_get_d
或者有描述如何执行此操作的参考吗?我翻了一下chapter 19 of RWO,没有看到这样的情况。
这不是一个完整的答案,但环顾四周,我发现 Zarith 在内部使用 GMP。有一个名为 mpq_get_d
的 GMP 函数可将有理数转换为双精度数。如果它不能直接在 Zarith 中使用,应该可以(给一些时间)为它添加一个接口。
如果您有权访问
- 整数
log2
运算, - 将整数左移给定位数的能力
然后进行您自己的正确舍入转换相对容易。简而言之,该方法如下所示:
- 减少到
n > 0
、d > 0
的大小写;过滤掉明显的 underflow/overflow - 选择一个整数
shift
,使2^-shift*n/d
介于2^54
和2^56
之间。 - 使用整数算法计算
x = 2^-shift*n/d
,使用 round-to-odd 舍入方法舍入到最接近的整数。 - 将
x
转换为最接近的 IEEE 754 双精度值dx
,使用通常的舍入到偶数舍入模式。 - Return
ldexp(dx, shift)
.
恐怕我对 OCaml 不是很流利,但下面的 Python 代码说明了正输入的想法。我留给你对负输入和除以零进行明显的修改。您可能还想为极端上溢和下溢的情况尽早设置 return:通过在下方查找 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
运算可能会引入 third 舍入。可以小心处理这个问题。请参阅下面的一些代码。
以上实际上是 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)