1

在尝试在没有 FPU 的情况下实现正弦函数时，我意识到所有输入都已经是理性的，因此我决定尝试一种全理性的方法。可能更慢但是：为什么不呢？该系列是线性收敛的，因此有机会通过二进制拆分进行运行时优化。甚至还有关于如何做到这一点以及我使用的非常详细的文献。

到目前为止，一切都很好。

我的数值算法原型工具是 Pari/GP，所以我将上面提到的论文中的代码移植到了 Pari/GP 中，正如您可能已经从我在这里发布问题的事实中猜到的那样，它不起作用。好吧，它确实有效，但无法将错误最小化。该论文还有其他几种不同功能的配方，但都表现出相同的行为。假设论文中有错字，我检查了作者在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)\2;
  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 如何工作的误解。它是哪一个，为什么？

4

2 回答 2

2

你会为此讨厌自己/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));
于 2019-08-15T19:22:41.997 回答
1

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

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

于 2021-07-14T19:19:48.390 回答