7

从数学、算法和元编程递归的角度来看，我有一个具有挑战性的问题。考虑以下声明：

template<class R1, class R2>
using ratio_power = /* to be defined */;

基于std::ratio类似的操作示例std::ratio_add。给定，二std::ratio R1和R2这个操作应该计算R1^R2当且仅当R1^R2是一个有理数。如果它是不合理的，那么实现应该会失败，比如当一个人试图将两个非常大的比率相乘时，编译器会说存在整数溢出。

三个问题：

  1. 您认为在不增加编译时间的情况下这是可能的吗？
  2. 使用什么算法？
  3. 如何实现这个操作？
4

2 回答 2

14

此计算需要两个构建块：

  • 编译时整数的 n 次方
  • 编译时整数的第 n 个根

注意：我使用 int 作为分子和分母的类型以节省一些打字，我希望重点能得到理解。我从一个工作实现中提取了以下代码，但我不能保证我不会在某处打错字；）

第一个相当简单：您使用 x^(2n) = x^n * x^n 或 x^(2n+1) = x^n * x^n * x 这样，您实例化的模板最少，例如x^39 计算如下： x39 = x19 * x19 * x x19 = x9 * x9 * x x9 = x4 * x4 * x x4 = x2 * x2 x2 = x1 * x1 x1 = x0 * x x0 = 1

template <int Base, int Exponent>
struct static_pow
{
  static const int temp = static_pow<Base, Exponent / 2>::value;
  static const int value = temp * temp * (Exponent % 2 == 1 ? Base : 1);
};

template <int Base>
struct static_pow<Base, 0>
{
  static const int value = 1;
};

第二个有点棘手，它使用包围算法：给定 x 和 N，我们想找到一个数字 r 使得 r^N = x

  • 将包含解的区间 [low, high] 设置为 [1, 1 + x / N]
  • 计算中点平均值 = (low + high) / 2
  • 确定，如果 mean^N >= x
    • 如果是，将间隔设置为 [low, mean]
    • 如果不是，则将间隔设置为 [mean+1, high]
  • 如果区间只包含一个数字，则计算完成
  • 否则，再次迭代

该算法给出了满足 s^N <= x 的最大整数 s

所以检查是否 s^N == x。如果是，则 x 的 N 次根是整数，否则不是。

现在让我们把它写成编译时程序：

基本界面：

template <int x, int N>
struct static_root : static_root_helper<x, N, 1, 1 + x / N> { };

帮手：

template <int x, int N, int low, int high>
struct static_root_helper
{
  static const int mean = (low + high) / 2;
  static const bool is_left = calculate_left<mean, N, x>::value;
  static const int value = static_root_helper<x, N, (is_left ? low : mean + 1), (is_left ? mean, high)>::value;
};

区间仅包含一个条目的递归端点：

template <int x, int N, int mid>
struct static_root_helper<x, N, mid, mid>
{
  static const int value = mid;
};

检测乘法溢出的助手（我认为，您可以将 boost-header 交换为 c++11 constexpr-numeric_limits）。如果乘法 a * b 溢出，则返回 true。

#include "boost/integer_traits.hpp"

template <typename T, T a, T b>
struct mul_overflow
{
  static_assert(std::is_integral<T>::value, "T must be integral");
  static const bool value = (a > boost::integer_traits<T>::const_max / b);
};

现在我们需要实现 calculate_left 来计算 x^N 的解是在均值的左侧还是在均值的右侧。我们希望能够计算任意根，因此像 static_pow > x 这样的简单实现会很快溢出并给出错误的结果。因此我们使用以下方案：我们要计算 x^N > B

  • 设置 A = x 和 i = 1
  • 如果 A >= B 我们已经完成 -> A^N 肯定会大于 B
  • A * x 会溢出吗？
    • 如果是 -> A^N 肯定会大于 B
    • 如果不是 -> A *= x 和 i += 1
  • 如果 i == N，我们就完成了，我们可以和 B 做一个简单的比较

现在让我们把它写成元程序

template <int A, int N, int B>
struct calculate_left : calculate_left_helper<A, 1, A, N, B, (A >= B)> { };

template <int x, int i, int A, int N, int B, bool short_circuit>
struct calulate_left_helper
{
  static const bool overflow = mul_overflow<int, x, A>::value;
  static const int next = calculate_next<x, A, overflow>::value;
  static const bool value = calculate_left_helper<next, i + 1, A, N, B, (overflow || next >= B)>::value;
};

i == N 的端点

template <int x, int A, int N, int B, bool short_circuit>
struct calculate_left_helper<x, N, A, N, B, short_circuit>
{
  static const bool value = (x >= B);
};

短路端点

template <int x, int i, int A, int N, int B>
struct calculate_down_helper<x, i, A, N, B, true>
{
  static const bool value = true;
};

template <int x, int A, int N, int B>
struct calculate_down_helper<x, N, A, N, B, true>
{
  static const bool value = true;
};

帮助器计算 x * A 的下一个值，考虑 x 溢出以消除编译器警告：

template <int a, int b, bool overflow>
struct calculate_next
{
  static const int value = a * b;
};

template <int a, int b>
struct calculate_next<a, b, true>
{
  static const int value = 0; // any value will do here, calculation will short-circuit anyway
};

所以，应该是这样的。我们需要一个额外的帮手

template <int x, int N>
struct has_integral_root
{
  static const int root = static_root<x, N>::value;
  static const bool value = (static_pow<root, N>::value == x);
};

现在我们可以实现 ratio_pow 如下：

template <typename, typename> struct ratio_pow;

template <int N1, int D1, int N2, int D2>
struct ratio_pow<std::ratio<N1, D1>, std::ratio<N2, D2>>
{
  // ensure that all roots are integral
  static_assert(has_integral_root<std::ratio<N1, D1>::num, std::ratio<N2, D2>::den>::value, "numerator has no integral root");
  static_assert(has_integral_root<std::ratio<N1, D1>::den, std::ratio<N2, D2>::den>::value, "denominator has no integral root");
  // calculate the "D2"-th root of (N1 / D1)
  static const int num1 = static_root<std::ratio<N1, D1>::num, std::ratio<N2, D2>::den>::value;
  static const int den1 = static_root<std::ratio<N1, D1>::den, std::ratio<N2, D2>::den>::value;
  // exchange num1 and den1 if the exponent is negative and set the exp to the absolute value of the exponent
  static const bool positive_exponent = std::ratio<N2, D2>::num >= 0;
  static const int num2 = positive_exponent ? num1 : den1;
  static const int den2 = positive_exponent ? den1 : num1;
  static const int exp = positive_exponent ? std::ratio<N2, D2>::num : - std::ratio<N2, D2>::num;
  //! calculate (num2 / den2) ^ "N2"
  typedef std::ratio<static_pow<num2, exp>::value, static_pow<den2, exp>::value> type;
};

所以，我希望至少能理解基本的想法。

于 2013-11-07T08:57:40.820 回答
5

是的，这是可能的。

让我们定义 R1 = P1/Q1、R2 = P2/Q2 和 R1^R2 = R3 = P3/Q3。进一步假设 P 和 Q 是互质数。

R1^R2 = R1^(P2/Q2) = R3
R1 ^ P2 = R3 ^ Q2.

R1^P2是已知的，并且有一个独特的因式分解到素数2^a * 3^b * 5^c * ...中注意它a, b, c可以是负数，因为 R1 是P1/Q1。现在第一个问题是是否所有a,b,c都是已知因子 Q2 的倍数。如果没有，那么你就直接失败了。如果是，那么R3 = 2^(a/Q2) * 3^(b/Q2) * 5^(c/Q2) ....

所有除法要么是精确的，要么结果不存在，因此我们可以在模板中使用纯整数数学。在模板中分解一个数字是相当简单的（部分专业化x%y==0）。

示例：2^(1/2) = R3 -> a=1, b=0, c=0, ... 和a%2 != 0-> 不可能。(1/9)^(1/2) -> a=0, b=-2, b%2 = 0, 可能, 结果 = 3^(-2/2)。

于 2013-11-06T23:09:43.993 回答