C 语言中 CORDIC 对数的问题
ema*_*uts 3 c c99 numerical-methods
为了开始使用 CORDIC for ，我实现了此 PDF 第 6 页log10中派生的算法：
#include <stdio.h>\n#include <math.h>\n\n// https://www.mikrocontroller.net/attachment/31117/cordic1.pdf\nfloat logf_cordic (float x, int B)\n{\n float z = 0.0f;\n\n for (int i = 1; i <= B; i++)\n {\n if (x > 1.0f)\n {\n printf ("-");\n x = x - ldexpf (x, -i);\n z = z - log10f (1.0f - ldexpf (1.0f, -i));\n }\n else\n {\n printf ("+");\n x = x + ldexpf (x, -i);\n z = z - log10f (1.0f + ldexpf (1.0f, -i));\n }\n }\n return z;\n}\n这是论文中代码的直接 C99 实现，我在其中ldexpf(x,n)计算 x\xc2\xb72 n。该论文声称该方法收敛于0.4194 < x < 3.4627。该程序使用1 \xe2\x89\xa4 x \xe2\x89\xa4 2. 下面打印完整的 C99 代码：
+--+--+--+-+-+++--++ x = 1.00000: y = -0.000000, dy = -1.487272e-07\n-+++++++++++++++++++ x = 1.12500: y = +0.099773, dy = +4.862081e-02\n-+++++++++++++++++++ x = 1.25000: y = +0.099773, dy = +2.863325e-03\n-+++-+--+-+--+-++--+ x = 1.37500: y = +0.138302, dy = -4.023314e-07\n-++-+--++----+++--+- x = 1.50000: y = +0.176091, dy = -2.831221e-07\n-+-+++++-++-++-++++- x = 1.62500: y = +0.210854, dy = +2.831221e-07\n-+-+-+-+++++--+---++ x = 1.75000: y = +0.243038, dy = +2.235174e-07\n-+--++-+--+---+---+- x = 1.87500: y = +0.273001, dy = +0.000000e+00\n-+---+--+++--------- x = 2.00000: y = +0.301030, dy = -5.960464e-08\n所以它按预期工作，除了x = 1.125, 1.25误差很大并且在计算更多迭代时不会减少的地方。我现在盯着那个代码几个小时，但找不到我缺少的东西......
\n\n
完整的C99代码
\n#include <stdio.h>\n#include <math.h>\n\nfloat logf_cordic (float x, int B)\n{\n float z = 0.0f;\n\n for (int i = 1; i <= B; i++)\n {\n if (x > 1.0f)\n {\n printf ("-");\n x = x - ldexpf (x, -i);\n z = z - log10f (1.0f - ldexpf (1.0f, -i));\n }\n else\n {\n printf ("+");\n x = x + ldexpf (x, -i);\n z = z - log10f (1.0f + ldexpf (1.0f, -i));\n }\n }\n return z;\n}\n\nint main (int argc, char *argv[])\n{\n int ex = 3;\n int B = 20;\n\n if (argc >= 2)\n sscanf (argv[1], "%i", &ex);\n\n if (argc >= 3)\n sscanf (argv[2], "%i", &B);\n\n if (ex < 0) ex = 0;\n if (ex > 16) ex = 16;\n if (B > 100) B = 100;\n\n int n = 1 << ex;\n float dx = 1.0f / n;\n for (int i = 0; i <= n; ++i)\n {\n float x = 1.0f + i * dx;\n float y = logf_cordic (x, B);\n float dy = y - log10f (x);\n printf (" x = %.5f: y = %+f, dy = %+e\\n",\n (double) x, (double) y, (double) dy);\n }\n return 0;\n}\n作为参考，这里是论文中提出的算法：
\nlog10(x){\n z = 0;\n for ( i=1;i=<B;i++ ){\n if (x > 1)\n x = x - x*2^(-i);\n z = z - log10(1-2^(-i));\n else\n x = x + x*2^(-i);\n z = z - log10(1+2^(-i));\n }\nreturn(z)\n}\n这不是一个完整的答案。这个问题引起了我的注意，我决定潜水并享受一些乐趣。
\n我已经使用浮点数测试了Python中论文中给出的算法，并且还转换为定点，我的结果与你的结果相同。我还发现了一些其他有问题的点，例如x = 1.60000，，x = 2.00625和。x = 2.03750x = 2.06875
我将使用以下约定：
\nb_i+ = (1 + 2^-i)\nb_i- = (1 - 2^-i)\n因此该算法可以写成：
\nlog10(x) {\n z = 0;\n for (i = 1; i =< B; i++) {\n if (x > 1)\n x *= b_i-;\n z -= log10(b_i-);\n else\n x *= b_i+;\n z -= log10(b_i+);\n }\n return(z)\n}\n我不确定我们可以选择b_i形式的说法(1 \xc2\xb1 2^-i)本身是错误的，但我注意到选择每种形式的标准b_i可能是错误的。例如，当 时x == 1.125，当前序列是
x *= b_1- == 0.562500, z -= log10(b_1-) == 0.301030\nx *= b_2+ == 0.703125, z -= log10(b_2+) == 0.204120\nx *= b_3+ == 0.791016, z -= log10(b_3+) == 0.152967\nx *= b_4+ == 0.840454, z -= log10(b_4+) == 0.126639\nx *= b_5+ == 0.866718, z -= log10(b_5+) == 0.113275\nx *= b_6+ == 0.880261, z -= log10(b_6+) == 0.106541\nx *= b_7+ == 0.887138, z -= log10(b_7+) == 0.103161\nx *= b_8+ == 0.890603, z -= log10(b_8+) == 0.101468\n ...\nx永远不会达到1并且 的值z是错误的。b_1+但是，如果我们从(而不是)开始b_1-，则顺序为：
x *= b_1+ == 1.687500; z -= log10(b_1+) == -0.176091\nx *= b_2- == 1.265625; z -= log10(b_2-) == -0.051152\nx *= b_3- == 1.107421; z -= log10(b_3-) == 0.006839\nx *= b_4- == 1.038208; z -= log10(b_4-) == 0.034868\nx *= b_5- == 1.005764; z -= log10(b_5-) == 0.048656\nx *= b_6- == 0.990048; z -= log10(b_6+) == 0.055495\nx *= b_7+ == 0.997783; z -= log10(b_7+) == 0.052116\nx *= b_8+ == 1.001681; z -= log10(b_8-) == 0.050422\n ...\n收敛到正确的结果。
\n对于x == 2.00625，现在该算法使用
b_1- b_2- b_3+ b_4+ b_5+ b_6+ b_7+ b_8+, ...\n并x稳定在约x == 0.957. 但如果我们改成b_2+相反，顺序是
b_1- b_2+ b_3- b_4- b_5- b_6+ b_7- b_8-\n这使得x == 1.0002和z == 0.323317（再次正确）。
我的直觉[1]是我们总是可以选择一系列满足, [2]b_i形式的 s ，但使用(1 \xc2\xb1 2^-i)x * b_1 * ... * b_B == 1
b_i = 1 + 2^-i if x * b1 * b2 * ... * b_i-1 < 1\nb_i = 1 - 2^-i if x * b1 * b2 * ... * b_i-1 > 1\n不是合适的选择标准。这种情况(x > 1.0)应该被一种更聪明的情况所取代，但到目前为止我还没有想到。
\n
[1] 事实证明我的直觉是错误的（再次！）。请参阅OP的回答。
\n[2] 因此使用查找表和简单的移位加法运算来计算对数，这就是该算法的目的。
\n该算法是假的。问题是并不是所有的数字都能以这种形式呈现
\nx = \xce\xa0 (1 \xc2\xb1 2 \xe2\x88\x92k )
\n请参阅有关MSE 的讨论。\n甚至不存在区间 [a,b] 使得该区间中的所有 x 都可以以该形式表示。
\n\n
但确实存在的是形式的表示
\nx = \xce\xa0 (1 + k \xc2\xb72 \xe2\x88\x92k )
\n对于所有 1 \xe2\x89\xa4 x \xe2\x89\xa4 2.38 = x 0 ，在 { 0, 1 } 中有k，这足以修复算法：
\n要计算 log B (x)，首先对参数 x 进行标准化，使得\n1 / x 0 \xe2\x89\xa4 x \xe2\x89\xa4 1。
\n然后将以下算法应用于减少的 x：
\nlogB_cordic (x, N)\n z = 0\n xmax = 1 + 2^{\xe2\x88\x92N\xe2\x88\x921}\n FOR k = 1 ... N\n xk = x + x\xc2\xb72^{-k}\n IF xk \xe2\x89\xa4 xmax\n x = xk;\n z = z - logB (1 + 2^{\xe2\x88\x92k})\n return z\n实际上，N是固定的并且 NlogB值来自查找表。除此之外，该方法需要添加、比较以及通过可变偏移量进行移位。
比较IF xk \xe2\x89\xa4 xmax使得最终的绝对误差（几乎）在 0 附近对称。IF xk \xe2\x89\xa4 1绝对误差将是不对称的，并且高两倍。
结果
\n作为参考，下面是 C99 中的实现，它增加了固定算法实际工作的信心。B = 2，N = 10，该图显示了在 [0.5, 1] 上均匀分布的 16385 个 x 值。
（点击放大）
C99代码
\n#include <math.h>\n\ndouble log2_cordic (double x, int N)\n{\n double xmax = 1.0 + ldexp (1.0, -N - 1);\n double z = 0.0;\n\n for (int k = 1; k <= N; k++)\n {\n double x_bk = x + ldexp (x, -k);\n\n if (x_bk <= xmax)\n {\n x = x_bk;\n z = z - log2 (1.0 + ldexp (1.0, -k));\n }\n }\n return z;\n}\n\n
编辑：
\n刚刚发现我重新发明了轮子......根据维基百科，这是费曼在曼哈顿计划期间已经使用的方法。
\n