编译时 std::ratio 的 std::ratio 功率？

Posted 2023-02-22

【中文标题】编译时 std::ratio 的 std::ratio 功率？

【英文标题】：std::ratio power of a std::ratio at compile-time?

【发布时间】：2013-11-06 21:30:28

【问题描述】：

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

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

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

三个问题：

你认为在不爆炸编译的情况下这是可能的吗 时间？ 使用什么算法？ 如何实现这个操作？

【问题讨论】：

已经有一个ratio_multiply 可以做到这一点。如果会发生整数溢出，我不知道您所说的“失败”是什么意思。你希望它如何失败？

ratio_multiply 是两个比率的乘积。在这里，我寻找两个比率的幂，例如：(num1/den1)^(num2/den2)。

计算机通常难以处理无理数。请举一些例子 - R2 是像1/2 这样的简单分数还是更复杂，像1/1337？

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

@MichaelSimbirsky：在模板元编程中，素数分解实际上并不难；这只是(X % prime == 0) 的部分专业化。如果为真，您有一个因素并继续X/prime；如果没有，请继续使用X % nextPrime。素数列表可以硬编码，这不太可能改变;)

【参考方案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

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

这个算法给出了满足 s^N 的最大整数 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;
;

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

【参考方案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, ... 和a%2 != 0 -> 不可能。 (1/9)^(1/2) -> a=0, b=-2, b%2 = 0, 可能, 结果 = 3^(-2/2)。