牛顿法分形生成中θ函数计算加速

Question

我一直在尝试生成牛顿迭代法分形的雅可比theta函数 - 我用mpmath尝试了很长时间，所以我尝试用C编写它。

用于生成以下图像的源代码在这里：http://owen.maresh.info/allegra.c，可以使用gcc allegra.c -o allegra -lm进行编译，然后应该调用./allegra > jacobi.pnm
_{(来源：maresh.info)}
所以： * 有没有办法加速评估 - 这需要半个小时的墙时间才能产生这个图像？（我想能够快速生成不同名称的这些图像，以便我可以制作电影） * 我知道我在θ函数定义中犯了一个错误，但我很难找到间断的原因。

供参考，此图像是通过对ϑ₃(z,0.001-0.3019*i)进行标准牛顿迭代得到的。

- graveolensa

---

你介意我把你的代码传递给苹果公司作为一个实例应用程序，演示使用他们新的clang C编译器编译会比使用gcc编译速度变慢的情况吗？ - Brian Swift

没问题。我稍后会尝试使用icc。 - graveolensa

"icc" 是一个很好的想法。我很好奇它相比于 "gcc" 会有多少改进。 - Brian Swift

icc » /usr/bin/time ./allegra2 >> chack.pnm 164.21用户 1.35系统 2:46.03经过的时间 99%CPU (0avgtext+0avgdata 14864maxresident)k 0输入+4920输出 (0major+985minor)页面错误 0交换 owen@orrery ~/math/thenewt » md5sum chack.pnm 03b3ed17194d6e77970310bce55c967d chack.pnm gcc » /usr/bin/time ./allegra >> geese.pnm 1257.28用户 23.44系统 21:32.45经过的时间 99%CPU (0avgtext+0avgdata 15376maxresident)k 32输入+4920输出 (0major+1010minor)页面错误 0交换 owen@orrery ~/math/thenewt » md5sum geese.pnm 03b3ed17194d6e77970310bce55c967d geese.pnm - graveolensa

哇，icc 生成的代码比 gcc 运行速度快了7倍以上，这真是令人印象深刻。（另外感谢提到 icc，我之前不知道英特尔在 Linux 上免费提供非商业用途的开发工具）。 - Brian Swift

附言：我已经完成了渲染，可以在YouTube上查看：http://www.youtube.com/watch?v=f0ZGfCmPjWA。下一步是尝试使用英特尔数学核心库（假设它比math.h更快）。 - graveolensa

Answer 1

首先尝试使用-O3和/或-fast启用编译器优化。在我的系统上进行快速测试，显示了三倍的性能提升。

另外，在尝试改进性能的代码更改时，将主循环更改为for(a=0;a<10 /* 512*/ ;a++)可能有助于缩短运行时间。

还要注意：GCC支持复数，请参阅手册中的complex、cpow和cexp，以及包含文件/usr/include/complex.h。

我对应用程序进行了剖析，并发现它大部分时间都花在了powc()上。不幸的是，当我将powc()更改为使用数学库中的cpow()时，它的运行速度比您的实现要慢。

如果您所在的系统具有多个核心，则可以通过OpenMP并行化外部主循环，从而很容易地降低墙钟时间。但是，当您为动画生成图像帧时，最有效的方法可能是使每个帧都由单独的进程生成（我喜欢使用xargs -P # -n 1进行这种粗粒度并行化）。

- Brian Swift

---

顺便说一下，我有一种直觉，认为没有人真正使用c99标准中的复杂数据类型做任何事情，只是因为当时有人认为这是个好主意才将其添加到语言中。（另外，我正在推动专用函数ASIC的开发，因为这些函数（Jacobi）theta在软件中被重新实现了无数次，而且所有的软件实现都非常慢。） - graveolensa

我认为在C99之前，gcc已经有了对复数的支持。然而，如果大多数重度用户仍在使用Fortran，我也不会感到惊讶。 - Brian Swift

Answer 2

当你在IRC上提到这个时，我情绪奇怪，花了一段时间进行优化。现在在我的Mac上至少快了4倍，不包括编译器优化标志，在其他一些平台上更快。

当涉及到高等数学时，我很无知，但我确实知道一些关于优化的事情。我认为这个计算与原始计算相同，除了用系统cexp()替换你的expc()实现，它产生相同的输出。你可以决定它是否仍然对你足够稳定。

正如Brian Swift所指出的，powc()是昂贵的，因为log()和pow()函数。

那些取得巨大成功的事情：

- pjtheta()和pjtheta3()中的计算可以合并 - 该计算可以成为newt()中的内部循环，并且其中的一些内容可以移出内部或两个循环 - cpow()可能对Brian（和我）来说速度较慢，但cexp()肯定比你的代码快，至少在我的机器上是这样。尝试两种方式 - 编译器标志中的-ffast-math会消除对不良行为数字的标准支持，并大大加快速度

另一个重大的胜利是将cexp()和cpow()中的算法转换为单精度，但这产生了略微不同的结果，您可能会或可能不会关心。

您可能无法再认出该程序，但它位于：

https://github.com/cgull/allegra.git

我注意到了一些事情，将其进一步优化了25%-33%（哇，这是一个收敛的迭代函数！）

我相信，比我更懂高等数学的人可以在其中找到另外2-4倍的性能。

- cgull