用 ldexp 反转 frexp

Posted 2023-02-23

【中文标题】用 ldexp 反转 frexp

【英文标题】：Inverting frexp with ldexp

【发布时间】：2015-10-24 04:58:49

【问题描述】：

如果我理解documentationcorrectly，我们应该能够使用ldexp 来恢复一个浮点数，该浮点数分解为frexp 的有符号尾数和指数。我一直无法做到这一点。考虑以下代码：

#include <cmath>
#include <iostream>
#include <limits>

template <typename T>
void float_info() 
    std::cout << "max=" << std::numeric_limits<T>::max()          <<
        ", max_exp="    << std::numeric_limits<T>::max_exponent   <<
        ", max_10_exp=" << std::numeric_limits<T>::max_exponent10 <<
        ", min="        << std::numeric_limits<T>::min()          <<
        ", min_exp="    << std::numeric_limits<T>::min_exponent   <<
        ", min_10_exp=" << std::numeric_limits<T>::min_exponent10 <<
        ", dig="        << std::numeric_limits<T>::digits10       <<
        ", mant_dig="   << std::numeric_limits<T>::digits         <<
        ", epsilon="    << std::numeric_limits<T>::epsilon()      <<
        ", radix="      << std::numeric_limits<T>::radix          <<
        ", rounds="     << std::numeric_limits<T>::round_style    << std::endl;


template <typename T>
void compare(T a, T b) 
    std::cout << a << " " << b << " (" <<
        (a != b ? "un" : "") << "equal)" << std::endl;


template<typename T>
void test_ldexp() 
    float_info<T>();

    T x = 1 + std::numeric_limits<T>::epsilon();
    T y = ldexp(x, 0);
    int exponent;
    T mantissa = frexp(x, &exponent);
    T z = ldexp(mantissa, exponent);

    compare(x, y);
    compare(x, z);

    std::cout << std::endl;


int main() 
    std::cout.precision(25);
    test_ldexp<float>();
    test_ldexp<double>();
    test_ldexp<long double>();

在 Ubuntu 14.04.3 LTS 上使用 g++（版本 4.8.4）编译时，输出为：

max=3.402823466385288598117042e+38, max_exp=128, max_10_exp=38,
min=1.175494350822287507968737e-38, min_exp=-125, min_10_exp=-37, dig=6,
mant_dig=24, epsilon=1.1920928955078125e-07, radix=2, rounds=1
1.00000011920928955078125 1.00000011920928955078125 (equal)
1.00000011920928955078125 1.00000011920928955078125 (equal)

max=1.797693134862315708145274e+308, max_exp=1024, max_10_exp=308,
min=2.225073858507201383090233e-308, min_exp=-1021, min_10_exp=-307, dig=15,
mant_dig=53, epsilon=2.220446049250313080847263e-16, radix=2, rounds=1
1.000000000000000222044605 1.000000000000000222044605 (equal)
1.000000000000000222044605 1.000000000000000222044605 (equal)

max=1.189731495357231765021264e+4932, max_exp=16384, max_10_exp=4932,
min=3.362103143112093506262678e-4932, min_exp=-16381, min_10_exp=-4931, dig=18,
mant_dig=64, epsilon=1.084202172485504434007453e-19, radix=2, rounds=1
1.00000000000000000010842 1 (unequal)
1.00000000000000000010842 1 (unequal)

当使用long doubles 时，我们通过将x 分解为frexpr 似乎会丢失一些东西。如果我使用 python3（版本 3.4.3）运行以下脚本，我可以实现我期望的行为。

import math
import sys

def compare(a, b):
    print('a b (preequal)'.format(a=a, b=b,
        pre='un' if a != b else ''))

x = 1 + sys.float_info.epsilon
mantissa, exponent = math.frexp(x)

print(sys.float_info)
compare(x, math.ldexp(x, 0))
compare(x, math.ldexp(mantissa, exponent))

输出是：

sys.float_info(max=1.7976931348623157e+308, max_exp=1024, max_10_exp=308,
min=2.2250738585072014e-308, min_exp=-1021, min_10_exp=-307, dig=15,
mant_dig=53, epsilon=2.220446049250313e-16, radix=2, rounds=1)
1.0000000000000002 1.0000000000000002 (equal)
1.0000000000000002 1.0000000000000002 (equal)

请注意，这仅使用doubles。

我试图阅读cmath 头文件以了解frexpr 和ldexpr 是如何实现的，但我无法理解它。怎么回事？

【参考方案1】：

通过包含 typeinfo 标头并在 test_ldexp 函数中对 compare 的调用之前插入以下行来修改您的 C++ 程序。

std::cout << "types:" << std::endl;
std::cout << "  x         : " << typeid(x).name() << std::endl;
std::cout << "  mantissa  : " << typeid(frexp(x, &exponent)).name()
          << std::endl;
std::cout << "  ldexp(...): " << typeid(ldexp(x, 0)).name() << std::endl;
std::cout << "  ldexp(...): " << typeid(ldexp(mantissa, exponent)).name()
          << std::endl;

现在的输出是：

max=3.402823466385288598117042e+38, max_exp=128, max_10_exp=38,
min=1.175494350822287507968737e-38, min_exp=-125, min_10_exp=-37, dig=6,
mant_dig=24, epsilon=1.1920928955078125e-07, radix=2, rounds=1
types:
  x         : f
  mantissa  : d
  ldexp(...): d
  ldexp(...): d
1.00000011920928955078125 1.00000011920928955078125 (equal)
1.00000011920928955078125 1.00000011920928955078125 (equal)

max=1.797693134862315708145274e+308, max_exp=1024, max_10_exp=308,
min=2.225073858507201383090233e-308, min_exp=-1021, min_10_exp=-307, dig=15,
mant_dig=53, epsilon=2.220446049250313080847263e-16, radix=2, rounds=1
types:
  x         : d
  mantissa  : d
  ldexp(...): d
  ldexp(...): d
1.000000000000000222044605 1.000000000000000222044605 (equal)
1.000000000000000222044605 1.000000000000000222044605 (equal)

max=1.189731495357231765021264e+4932, max_exp=16384, max_10_exp=4932,
min=3.362103143112093506262678e-4932, min_exp=-16381, min_10_exp=-4931, dig=18,
mant_dig=64, epsilon=1.084202172485504434007453e-19, radix=2, rounds=1
types:
  x         : e
  mantissa  : d
  ldexp(...): d
  ldexp(...): d
1.00000000000000000010842 1 (unequal)
1.00000000000000000010842 1 (unequal)

frexpr 和 ldexpr 将返回 doubles 无论您输入什么类型！您似乎正在使用math.h 中定义的函数（参见here 和here），而不是cmath 中定义的函数。将您对frexpr 和ldexpr 的调用替换为对std::frexpr 和std::ldexpr 的调用，您的代码将按预期运行。