编译时如何计算std::ratio的幂？

Question

我有一个具有数学、算法和元编程递归观点的具有挑战性的问题。考虑以下声明：

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. 如何实现此操作？

1. 执行这个操作而不增加编译时间是不可能的。 2. 可以参考二进制分解幂运算的算法。 3. 可以使用模板元编程技术（template metaprogramming）来实现此操作，具体而言，可以使用递归模板函数计算幂，同时使用编译时if（constexpr if）避免不必要的计算和类型转换。

- Vincent

---

已经有一个 ratio_multiply 可以实现这个功能。但是，如果整数溢出会发生什么“失败”我不知道你的意思。你想如何让它失败？ - David

ratio_multiply 是两个比率的乘积。在这里，我寻求两个比率的指数运算，例如：(num1/den1)^(num2/den2)。 - Vincent

计算机通常难以处理无理数。请举一些例子 - R2是否总是一个简单的分数，如 1/2，还是可以更加复杂，比如 1/1337？ - anatolyg

R1和R2可以是任何有效的std::ratio（例如1/1337）。如果结果是无理数（例如2^(1/2)），则元函数在实例化期间应崩溃。 - Vincent

众所周知，元编程是图灵机。这意味着它是可能的，并且应该沿着实现运行时解决方案的方式进行实现。由于涉及到质因数分解，不要期望有简单的解决方案。 - Michael Simbirsky

@MichaelSimbirsky：在模板元编程中，质因数分解实际上并不难；它只是对(X % prime == 0)进行部分特化。如果为真，则您有一个因子，并继续使用X/prime；如果不是，则继续使用X % nextPrime。质数列表可以硬编码，这很不可能改变 ;) - MSalters

Answer 1

你需要两个构建块来进行此计算：

• 在编译时计算整数的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。

1. 设置包含解的区间[low, high]为[1, 1 + x / N] 2. 计算中点 mean = (low + high) / 2 3. 确定 mean^N 是否大于等于 x - 如果是，则将区间设置为 [low, mean] - 如果不是，则将区间设置为 [mean+1, high] 4. 如果区间只包含一个数字，则计算完成 5. 否则，重复上述步骤

该算法可以得出最大的整数s，满足 s^N <= x

因此，检查 s^N 是否等于 x。如果是，则x的N次方根是整数，否则不是。

现在让我们将其编写为编译时程序：

基本接口：

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

helper:

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;
};

检测乘法溢出的辅助程序（您可以使用c++11 constexpr-numeric_limits替换boost-header，我认为）。如果乘法 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;
};

所以，我希望至少基本思想能够传达出来。

- MadScientist

Answer 2

是的，这是可能的。

让我们定义 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, ... and a%2 != 0 -> 不可能。 (1/9)^(1/2) -> a=0, b=-2, b%2 = 0, 可能, 结果为3^(-2/2)。

- MSalters