apfloat源码分析经验交流

[liangbch]

总的来说，apfloat 写的非常好，其大数乘法使用快速数论变换实现，速度很快。但FNT的实现也很复杂，不容易弄懂。特开此贴，方便大家交流 apfloat源码分析的心得，便于共同提高。

[shines77]

看得还不是很懂，但FNT流程已经很清晰了，apfloat的代码还是写得很清爽的，先up一下

CRT我已经弄懂，但FNT具体的转换细节还有些不太明白的地方

[wenminghu]

看得头大,楼上的两位可以留下msn么?

[gxqcn]

别急于看代码，应先看文档，把脉络先领会透。
（我就是这样过来的，写出了比apfloat效率更高的代码）

[yaos]

一般是类似2重循环的一个结构
内部实现a + b
(a + b)c^k的运算在数组上

[liangbch]

最核心最复杂的是6步FNT变换。

[liangbch]

按照 IEEE754 标准，double 型数为 8byte(64bits)，其中尾数位 52+1 (因为最高bit 一定为1，故实际中部存储，可节省1bit)，但是在FPU内部，浮点处理器是80位，故可设为64位精度的表示。


Sign: 1 Bit; Exponent: 11 Bits; Mantissa: 52 Bits

Value = (-1)^S X 2^(E-1023) X (1.M)

Range: +/- 2^(-1022) to (2^1023)x(2 – (2^-52))
+/- 2.23 X (10^-308) to 1.8 x (10^308)

Example:
0 = 0 00000000000 00000……0
-1.5 = 1 01111111111 10000..….0
3.75 = 0 10000000000 11100……0

[shines77]

谢谢liangbch(宝宝)，刚才我也想到了

因为http://www.pediy.com/tutorial/chap2/Chap2-4.htm
有讲mem80，另外好像在MSDN还是哪里也提到过80bit，由于对FPU不是很了解没有注意

刚才我测试了一下，SSE的float的复数乘法运算比c写的快了30%，不知道double会快多少，

对于FNT的6步FNT变换，也是我现在写FFT的参照，不过我最近可能需要放慢一下脚步，争取早日走上FNT的道路

[yaos]

_control87 (RC_CHOP | PC_64, MCW_RC | MCW_PC);似乎是x87控制字吧 ：）

[yaos]

SSE的float的复数乘法运算比c写的快了30%，不知道double会快多少

：），似乎不止吧，最少50% ：）

[wenminghu]

不过我更喜欢先看实现.

[shines77]

_control87 (RC_CHOP | PC_64, MCW_RC | MCW_PC);

这个的意思是不是把FPU当作64bit大整数做乘除法，因为我实在无法理解64位精度的浮点怎么表示，_PC_24 可以认为是float，_PC_53 应该是double，如果是这样，那是不是用SSE2会更快呢，可惜我的CPU不支持SSE2，所以我写的SSE2的复数乘法(用于FFT,double精度)也测试不了

[shines77]

http://www.77studio.net/files/cmplx_test.rar

我做了一下16字节数据对齐，结果还是快了30%左右，里面有double版的SSE2代码，顺便帮我测一下，我的CPU不支持SSE2，跟C的比一下

[yaos]

C的：

BASE = 10000

Whether display result on screen? [0=NO, 1=YES]: 0

BigInt Length = [1-131072, 0 = exit] ? 131072

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

cmplx_multiply_double [C version]:

cmplx_multiply length: 262144

cmplx_multiply used time: 0.00818647 s.

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Total time: 1.09034483 s.

FFTSquareWorstError: 0.000000000

SSE2的：
BASE = 10000

Whether display result on screen? [0=NO, 1=YES]: 0

BigInt Length = [1-131072, 0 = exit] ? 131072

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

cmplx_multiply_double [SSE version]:
cmplx_multiply length: 262144

cmplx_multiply used time: 0.01037766 s.

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Total time: 1.19107649 s.

FFTSquareWorstError: 0.000000000

[yaos]

嘿嘿，你这是怎么优化的？ ：）

[yaos]

C:
BASE = 10000

Whether display result on screen? [0=NO, 1=YES]: 0

BigInt Length = [1-131072, 0 = exit] ? 131072

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

cmplx_multiply_float [C version]:

cmplx_multiply length: 524288

cmplx_multiply used time: 0.39171647 s.

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Total time: 10.05411307 s.

FFTSquareWorstError: 0.000000000

SSE:
BASE = 10000

Whether display result on screen? [0=NO, 1=YES]: 0

BigInt Length = [1-131072, 0 = exit] ? 131072

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

cmplx_multiply_float [SSE version]:
cmplx_multiply length: 524288

cmplx_multiply used time: 0.28453780 s.

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Total time: 9.58576214 s.

FFTSquareWorstError: 0.000000000

[shines77]

SSE2的写得有好大问题，我改一下，待会另外开贴讨论，顺便把代码共享一下

[shines77]

帖子地址：
http://community.csdn.net/Expert/topic/3206/3206642.xml?temp=.5171167

代码以及测试exe下载：
http://www.77studio.net/files/cmplx_src.rar （那个帖子主题的下载地址是错位的，一时疏忽）

[liangbch]

最朴素的fnt正变换和逆变换代码。
void fnt (modint data[], modint pr, int isign, size_t nn, int s)
{
size_t m, i, j, istep, mmax;
modint w, z, a, b;

if (nn < 2) return;

if (isign > 0)
w = pow (pr, modint::modulus - 1 - (modint::modulus - 1) / nn);
else
w = pow (pr, (modint::modulus - 1) / nn);

mmax = nn >> 1;
while (mmax)
{
istep = mmax << 1;

// Optimize first step when z = 1

for (i = 0; i < nn; i += istep)
{
j = i + mmax;
a = data[i];
b = data[j];
data[i] = a + b;
data[j] = a - b;
}
z = w;

for (m = 1; m < mmax; m++)
{
for (i = m; i < nn; i += istep)
{
j = i + mmax;
a = data[i];
b = data[j];
data[i] = a + b;
data[j] = z * (a - b);
}
z *= w;
}
w *= w;
mmax >>= 1;
}

if (s) scramble (data, nn);
}


void ifnt (modint data[], modint pr, int isign, size_t nn, int s)
{
size_t m, i, j, istep, mmax;
modint w, wt, wr, wtemp;

if (nn < 2) return;

if (isign > 0)
w = pow (pr, modint::modulus - 1 - (modint::modulus - 1) / nn);
else
w = pow (pr, (modint::modulus - 1) / nn);

if (s) scramble (data, nn);

mmax = 1;
while (nn > mmax)
{
istep = mmax << 1;
wr = wt = pow (w, nn / istep);

// Optimize first step when wr = 1

for (i = 0; i < nn; i += istep)
{
j = i + mmax;
wtemp = data[j];
data[j] = data[i] - wtemp;
data[i] += wtemp;
}

for (m = 1; m < mmax; m++)
{
for (i = m; i < nn; i += istep)
{
j = i + mmax;
wtemp = wr * data[j];
data[j] = data[i] - wtemp;
data[i] += wtemp;
}
wr *= wt;
}
mmax = istep;
}
}


下面是使用最普通的 fnt 计算 大数相乘的代码，代码中调用的函数请参照apfloat。

void fnt_DivN (modint data[], size_t nn)
{
for (int i=0;i<nn;i++)
data[i] /= nn;
}

void convolution1 (const DWORD *buf1, const DWORD *buf2, DWORD *prod,size_t nn, int *i)
{
modint *rst1;
modint *rst2;
modint *rst3;
modint *tmp;

size_t n2 = (nn & -nn);
assert(nn == n2);
//----------------------------------------------
tmp = new modint[nn];
//----------------------------------------------
setmodulus (moduli[2]);

rst1 = new modint[nn];
memcpy(rst1,buf1,sizeof(modint)*nn);
fnt (rst1, primitiveroots[2], 1, nn);

memcpy(tmp,buf2,sizeof(modint)*nn);
fnt (tmp, primitiveroots[2], 1, nn);

multiplyinplace (rst1, tmp,nn);
fnt (rst1, primitiveroots[2], -1, nn);
fnt_DivN (rst1, nn);
//====================================
setmodulus (moduli[1]);

rst2 = new modint[nn];
memcpy(rst2,buf1,sizeof(modint)*nn);
fnt (rst2, primitiveroots[1], 1, nn);

memcpy(tmp,buf2,sizeof(modint)*nn);
fnt (tmp, primitiveroots[1], 1, nn);

multiplyinplace (rst2, tmp,nn);
fnt (rst2, primitiveroots[1], -1, nn);
fnt_DivN (rst2, nn);
//====================================

setmodulus (moduli[0]);

rst3 = new modint[nn];
memcpy(rst3,buf1,sizeof(modint)*nn);
fnt (rst3, primitiveroots[0], 1, nn);

memcpy(tmp,buf2,sizeof(modint)*nn);
fnt (tmp, primitiveroots[0], 1, nn);

multiplyinplace (rst3, tmp,nn);
fnt (rst3, primitiveroots[0], -1, nn);
fnt_DivN (rst3, nn);
//-------------------------------------------------
*i = carrycrt (rst3, rst2, rst1, nn);
memcpy(prod,rst3,sizeof(modint)*nn);

delete[] rst1;
delete[] rst2;
delete[] rst3;
delete[] tmp;
}