Pari/GP 中的二进制拆分

Question

在不使用 FPU 的情况下尝试实现正弦函数时，我认识到所有输入都已经是有理数了，因此我决定尝试一种全有理数方法。可能更慢但是：为什么不呢？该级数是线性收敛的，因此有机会通过二进制拆分进行运行时间优化。甚至还有 highly detailed literature 如何做以及我用过的。

到目前为止，还不错。

我的数值算法原型设计工具是 Pari/GP，所以我将上述论文中的代码移植到 Pari/GP，正如您可能已经猜到的那样，我 post 在这里提问，没用。好吧，它确实有效，但无法将错误最小化。该论文还有其他几种不同功能的方法，但都显示出相同的行为。假设论文中有错字，我检查了作者在 CLN 中的实现。高度优化但基于论文中的代码，甚至是逐字逐句。

为了获得 MWE，我使用了他们的 exp(p/q) 配方（除了阶乘之外最简单的配方）并简化了 Pari/GP 代码。

exp_bin_split_rat_internal(n1, n2, x) = {
  \ R, L, r = [P, Q, B, T]
  \ a       = [p, q, b, a]
  local(diff, mn, L, R, r = vector(4));

  diff = n2 - n1;
  if(diff == 0,
     \ no actual error-handling here
     print("Error in bin_split_rat_internal: n2-n1 is zero.");
  );
  if( diff == 1,
    \ x = u/v
    if(n1 == 0,
      \ r.P = 1;
      r[1] = 1;
      \ r.Q = 1;
      r[2] = 1;
      \ r.B = b(0) = 1;
      r[3] = 1;
      \ r.T = a(0) * r.P = 1 * u;
      r[4] = 1 * r[1];
      return(r);
    , \ else
      \ r.P = u;
      r[1] = numerator(x);
      \ r.Q = n1 * v;
      r[2] = n1 * denominator(x);
      \ r.B = b(n) =  1;
      r[3] =  1;
      \ r.T = a(n) * r.P = 1 * u;
      r[4] = 1 * r[1];
      return(r);
    );
  );
  \ floor((n1 + n2)/2)
  nm = (n1 + n2);
  L = exp_bin_split_rat_internal(n1, nm, x);
  R = exp_bin_split_rat_internal(nm, n2, x);
  \            1  2  3  4
  \ R, L, r = [P, Q, B, T]
  \ r.P = L.P * R.P;
  r[1] = L[1] * R[1];
  \r.Q = L.Q * R.Q;
  r[2] = L[2] * R[2];
  \r.B = L.B * R.B;
  r[3] = L[3] * R[3];
  \r.T = R.B  * R.Q  * L.T   +  L.B  * L.P  * R.T;
  r[4] = (R[3] * R[2] * L[4]) + (L[3] * L[1] * R[4]);

  return(r);
}

exp_bin_split_rat(x, n) = {
  local(r, ret);
  r = exp_bin_split_rat_internal(0, n, x);
  \ r = [P, Q, B, T]
  \ S = T/(B*Q)
  ret = r[4]/(r[3] * r[2]);
  return(ret);
}

k = 1/1234;
tmp = exp_bin_split_rat(k, 10) * 1.0;print(tmp);tmp= exp(k);print(tmp);
tmp = exp_bin_split_rat(k, 100) * 1.0;print(tmp);tmp= exp(k);print(tmp);
tmp = exp_bin_split_rat(k, 1000) * 1.0;print(tmp);tmp= exp(k);print(tmp);
tmp = exp_bin_split_rat(k, 10000) * 1.0;print(tmp);tmp= exp(k);print(tmp);

（最后一个可能会有点长，想看的可以跳过运行）

如您所见，经过几十步后，误差不可能进一步减小。

所以这要么是我对算法的工作方式有误解，要么是 Pari/GP 的工作方式有误。它是哪一个，为什么？

Answer 1

你会因此讨厌 yourself/PARI（但喜欢堆栈溢出）。

在第一行你有：

local(diff, mn, L, R, r = vector(4));

而不是：

local(diff, nm, L, R, r = vector(4));

我建议使用 my 而不是 local。

my(diff, nm, L, R, r = vector(4));

Answer 2

请注意，从 2.13 版开始，它作为 Pari/GP 中的内置函数可用：请参阅 trans1.c:expQ()。计算小高度有理数的 exp(p/q) 时使用此函数。

通用驱动程序是私有的，但总有一天会导出，请参阅trans1.c：abpq_sum()。

Pari/GP 中的二进制拆分

Binary splitting in Pari/GP

rational-number

pari

exp

pari-gp