公布源码和优化报告了

[huanyun]

兑现我的承诺

在T2050笔记本上测试结果双精度0.671秒 单精度:0.312秒

希望能看到其他高手的代码 学习学习! :)

下载地址如下

http://www.pudn.com/downloads68/sourcecode/windows/other/detail244074.html

[huanyun]

我是用2005编译的 我把程序放上去 麻烦你找个core 2跑一下
我的程序还没有在core2上跑过 :)

http://download.csdn.net/source/160802

[flyingdog]

我的E6600 和 xeno 5130都是借来的。
现在再要跑，有困难啊。

[flyingdog]

我的pow部分，是靠icc优化的，没有用汇编优化。用vc2005编译吃亏了。

[huanyun]

昏 你能借到xeno? 我从来没有用过那些处理器 寒

[huanyun]

我用2005 + icc编译的

[flyingdog]

其实xeno也没怎么样的。E6600 和 xeno 5130基于同样的内核。E6600 2.4G xeno 5130 2.0G
结果还是E6600 快。
那个xeno就前端总线快了一点，对于这个程序没什么用处。


我借到的还是双xeno 5130的，还特地屏蔽掉了一个才做程序的。

其实关键是酷锐的核。否则p4 xeno多的是。

[huanyun]

估计是你使用了两个r数组
我的t2050 cache 不够大

[flyingdog]

t2050？很神奇的型号。
intel网站上查不到。
网上查也没看明白是什么版本的。
t2050的L1 cache是多大的？

[huanyun]

笔记本啊 dell 640m 就是了

[huanyun]

昨天查了精度的原因

__inline double InvSqrt_Nr(double x)
{
double xhalf=0.5f*x;
__int64 i= *(__int64*) &x;
i = 0x5fe6ec85e7de30da - (i>>1);
x = *(double*)&i;
x = x*(1.5f-xhalf*x*x);
x = x*(1.5f-xhalf*x*x);
x = x*(1.5f-xhalf*x*x);
return x;
}
少了一次迭代 可改为
__inline double InvSqrt_Nr(double x)
{
double xhalf=0.5f*x;
__int64 i= *(__int64*) &x;
i = 0x5fe6ec85e7de30da - (i>>1);
x = *(double*)&i;
x = x*(1.5f-xhalf*x*x);
x = x*(1.5f-xhalf*x*x);
x = x*(1.5f-xhalf*x*x);
x = x*(1.5f-xhalf*x*x);
return x;
}

由于该函数使用次数极少 仅在XGetComputePotPoint中使用
__inline double XGetComputePotPoint(int i, int j)
{
int ii = (i >> 1) << 1;
int jj = (j >> 1) << 1;
int mm = i & 1;
int nn = j & 1;
return InvSqrt_Nr((r[ii][mm] - r[jj][nn]) * (r[ii][mm] - r[jj][nn]) +
(r[ii][mm+2] - r[jj][nn+2]) * (r[ii][mm+2] - r[jj][nn+2]) +
(r[ii][mm+4] - r[jj][nn+4]) * (r[ii][mm+4] - r[jj][nn+4]));
}
所以并不影响速度

[huanyun]

如flyingdog(flyingdog) 所说的 CHRIS LOMONT的值需要四次迭代才可以达到7位精度的要求

[flyingdog]

你试试用这个值能否三次迭代达到精度要求。
0x5fe6d250b0000000

你的程序我没试过。我当初在测精度时发现，如果全部都用那种方法估计初值，那使用0x5fe6d250b0000000会好一点。

[huanyun]

其实
return InvSqrt_Nr((r[ii][mm] - r[jj][nn]) * (r[ii][mm] - r[jj][nn]) +
(r[ii][mm+2] - r[jj][nn+2]) * (r[ii][mm+2] - r[jj][nn+2]) +
(r[ii][mm+4] - r[jj][nn+4]) * (r[ii][mm+4] - r[jj][nn+4]));
可以该为
return invsqrt((r[ii][mm] - r[jj][nn]) * (r[ii][mm] - r[jj][nn]) +
(r[ii][mm+2] - r[jj][nn+2]) * (r[ii][mm+2] - r[jj][nn+2]) +
(r[ii][mm+4] - r[jj][nn+4]) * (r[ii][mm+4] - r[jj][nn+4]));

我担心编译过不了 才改用InvSqrt_Nr

[huanyun]

-----------------0x5fe6d250b0000000-----------------
0: Potential: 686445.9565352
10: Potential: 280129.9852357
20: Potential: 206437.9449030
30: Potential: 170828.5177899
40: Potential: 147831.4599590
50: Potential: 133650.1392165
60: Potential: 122234.9888389-----?
70: Potential: 114473.8684963
80: Potential: 106560.9887754
90: Potential: 99199.0808710
100: Potential: 93318.9772015
110: Potential: 89056.6733336
120: Potential: 85229.7706144
130: Potential: 82200.7831734
140: Potential: 79693.6008187
150: Potential: 77797.8899493
160: Potential: 75235.8488249
170: Potential: 73208.5403272
180: Potential: 71404.9412059
190: Potential: 69980.3005699
200: Potential: 68477.5907029
Seconds = 0.687000000

使用0x5fe6d250b0000000 还是一样


-----------------Begin-----------------
0: Potential: 686445.9565352
10: Potential: 280129.9852357
20: Potential: 206437.9449030
30: Potential: 170828.5177899
40: Potential: 147831.4599590
50: Potential: 133650.1392165
60: Potential: 122234.9888390
70: Potential: 114473.8684963
80: Potential: 106560.9887754
90: Potential: 99199.0808710
100: Potential: 93318.9772015
110: Potential: 89056.6733336
120: Potential: 85229.7706144
130: Potential: 82200.7831734
140: Potential: 79693.6008187
150: Potential: 77797.8899493
160: Potential: 75235.8488249
170: Potential: 73208.5403272
180: Potential: 71404.9412059
190: Potential: 69980.3005699
200: Potential: 68477.5907029
Seconds = 0.687000000

这个完全正确 是四次迭代的结果

[flyingdog]

埃，看来不同情况下，幻数的选择还是不一样。
你有兴趣的话，可以多试试其他值，肯定能找到合适的幻数的。
我本来想穷举，找到最好的值的。但是后来发现rsqrtps更好用。就放弃穷举了。

[huanyun]

我的算法99%的计算是用__mm_rsqrt_ps 用单精度牛顿下山法提高精度 然后再用双精度牛顿下山法提高精度 所以也不用在意魔法数的选择

可以看到四次迭代的时间并没有增加

另外 我在我的机器上编译测试flyingdog 和housisong 的程序
都使用我设置的编译选项 估计是编译选项我设置的不好

flyingdog
200: Potential: 68477.5907029
Seconds = 0.796000000

housisong
200: Potential: 68477.5907029
Seconds = 0.734000000.

我的
200: Potential: 68477.5907029
Seconds = 0.687000000


我估计是我使用了 SetPriorityClass(hProcess, REALTIME_PRIORITY_CLASS);
提高了成绩

[huanyun]

都加入SetPriorityClass(hProcess, REALTIME_PRIORITY_CLASS) 后

flyingdog
200: Potential: 68477.5907029
Seconds = 0.781000000

housisong
200: Potential: 68477.5907029
Seconds = 0.703000000

[housisong]

我也用穷举的方法想找到一个更好的魔法数，也没有成功；
（但也没有排除，如果针对这个竞赛问题的值域来找到一个更好的魔法数；
从我的测试中针对我的一组特定数据，还是能找到比0x5fe6ec85e7de30da更好的数的)
）

[huanyun]

使用魔法数需要至少3次迭代
而__mm_rsqrt_ps只需要两次 一次float 一次double

我是为了处理边角数据才使用自己写得InvSqrt_Nr
其实intel库里的invsqrt 就可以了

[阿浩No_1]

收藏