【经验总结】不能实施并行处理的情况	浅谈并行编程中的任务分解模式

"what the fuck ?"

楼主ed9er（始祖鸟）2000-12-28 20:23:00 在专题开发/技术/项目 / 数据结构与算法提问

/*
** float q_rsqrt( float number )
*/
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;

x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // what the fuck?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

return y;
}

牛顿法求倒数平方，我自己看不懂，谁能讲解清楚？问题点数：100、回复次数：34Top

1 楼ed9er（始祖鸟）回复于 2000-12-28 20:25:00 得分 0

敲错了，是倒数开方，开平方Top

2 楼Wingsun（孙春阳）回复于 2000-12-28 21:09:00 得分 0

i = 0x5f3759df - ( i >> 1 ); // what the fuck?
这句话不就等于
i=0x5f3759df - ( i/2 );
0x5f3759df 肯定是一个常数啊。
Top

3 楼ed9er（始祖鸟）回复于 2000-12-29 10:06:00 得分 0

楼上的，怎么不顺便给我讲解一下long是什么意思？float是什么意思？

不会要我先把牛顿法讲解一遍吧：
x(k+1)=x(k)-f(x(k))/f'(x(k))
求初始点x0的过程看不懂，求得的x0已经很靠近解点，两
次迭代已经很精确了，要求不高时甚至不需要迭代。
Top

4 楼The_east_key（东方一键）回复于 2000-12-29 14:23:00 得分 0

建议您访问www.etechbase.net/tech,里面有很多资料，也许可以解决您的问题。
访问http://168.168.18.11:81/etechbase/advsearch.php将您的问题输入查询内容框，选择不同的精确程度，即可以找到你所需要的答案。效果还是可以的。Top

5 楼ed9er（始祖鸟）回复于 2001-01-16 15:28:00 得分 0

翻出来
Top

6 楼Kevin_qing（）回复于 2001-02-26 13:20:00 得分 0

你明白这句吗
i = * ( long * ) &y;

Top

7 楼ed9er（始祖鸟）回复于 2001-02-26 14:05:00 得分 0

这句应该是为下面关键的一句做准备的，下面一句完了以后马上又转回来了
这个跟float的机内表示关系很大Top

8 楼elaiyo（小牛头）回复于 2001-02-26 16:51:00 得分 0

这里只是我的理解，提出来探讨一下。：）
i = * ( long * ) &y; /*强行将一个数从float型转换为长整型。*/
i = 0x5f3759df - ( i >> 1 ); /* 关键部分*/
y = y * ( threehalfs - ( x2 * y * y ) ); /* 重新从long型转换为长整型*/
我觉得这个算法是非常奇怪的，真的要弄懂的话，必须要想清楚从float型转化为长整型发生了什么事。这涉及到一个数字在计算机中的表达。还有相关特有的算法。
ed9er,很佩服你。可否做个朋友？？？
我的email : elaiyo@etang.com
Top

9 楼elaiyo（小牛头）回复于 2001-02-26 17:11:00 得分 0

另外，ed9er兄，对下面程序的printf语句，我还是有点不明白
#include <stdio.h>
main(_){char*x="=#6#694.21..#62.*.*2-)1#62.*..*.=#=1..29..#";
while(_=*x/4)_-=8,printf("\n%*s"+!!_,_+_,"_/_/_/"+*x++%4*2);}
问题如下：
（1）%*s 是否真的还代表一个字符串输出？？
（2）其中有两个逗号，应是给出了两个参数，为什么只有一个%s？？？
谢谢
Top

10 楼holyfire（谁最衰啊你最衰，谁最帅啊我最帅）回复于 2001-02-26 17:36:00 得分 0

i = 0x5f3759df - ( i >> 1 ); =>> i = 0x5f3759df - i /2;
=>> i = 0x5f3759df *2 / 3 ;
Top

11 楼holyfire（谁最衰啊你最衰，谁最帅啊我最帅）回复于 2001-02-26 17:41:00 得分 0

有点怪。看来它是利用了机器的误差。我看不出什么道道来。Top

12 楼holyfire（谁最衰啊你最衰，谁最帅啊我最帅）回复于 2001-02-26 17:43:00 得分 0

我弄错了i = 0x5f3759df *2 / 3 ;
不成立。
不好意思。Top

13 楼ed9er（始祖鸟）回复于 2001-02-26 17:45:00 得分 0

感兴趣越来越多了~~关键问题是0x5f3759df是如何来的，其他都是小事，IEEE的float规范：

Single Precision
The IEEE single precision floating point standard representation requires a 32 bit word, which may be represented as numbered from 0 to 31, left to right. The first bit is the sign bit, S, the next eight bits are the exponent bits, 'E', and the final 23 bits are the fraction 'F':

S EEEEEEEE FFFFFFFFFFFFFFFFFFFFFFF
0 1 8 9 31

The value V represented by the word may be determined as follows:

If E=255 and F is nonzero, then V=NaN ("Not a number")
If E=255 and F is zero and S is 1, then V=-Infinity
If E=255 and F is zero and S is 0, then V=Infinity
If 0<E<255 then V=(-1)**S * 2 ** (E-127) * (1.F) where "1.F" is intended to represent the binary number created by prefixing F with an implicit leading 1 and a binary point.
If E=0 and F is nonzero, then V=(-1)**S * 2 ** (-126) * (0.F) These are "unnormalized" values.
If E=0 and F is zero and S is 1, then V=-0
If E=0 and F is zero and S is 0, then V=0
In particular,

0 00000000 00000000000000000000000 = 0
1 00000000 00000000000000000000000 = -0

0 11111111 00000000000000000000000 = Infinity
1 11111111 00000000000000000000000 = -Infinity

0 11111111 00000100000000000000000 = NaN
1 11111111 00100010001001010101010 = NaN

0 10000000 00000000000000000000000 = +1 * 2**(128-127) * 1.0 = 2
0 10000001 10100000000000000000000 = +1 * 2**(129-127) * 1.101 = 6.5
1 10000001 10100000000000000000000 = -1 * 2**(129-127) * 1.101 = -6.5

0 00000001 00000000000000000000000 = +1 * 2**(1-127) * 1.0 = 2**(-126)
0 00000000 10000000000000000000000 = +1 * 2**(-126) * 0.1 = 2**(-127)
0 00000000 00000000000000000000001 = +1 * 2**(-126) *
0.00000000000000000000001 =
2**(-149) (Smallest positive value)

// *s要两个参数，第一个是整个串打出来的长度，后一个串的内容Top

14 楼Smile_Tiger（笑面虎）回复于 2001-02-27 16:29:00 得分 0

请教各位牛顿迭代法的具体描述Top

15 楼Smile_Tiger（笑面虎）回复于 2001-02-27 16:41:00 得分 0

另外 elaiyo(牛气冲天)关于 i = * ( long * ) &y;的说法是错误的

i = * ( long * ) &y; 的意思是将float类型的y变量地址以long类型解释，然后赋值给i变量
此处注意到 long类型和float类型都是四字节的，

如果是"强行将一个数从float型转换为长整型",那么应该是这样的：
i = (long)y; 或者 i = float(y);(c++特有的表示法)
Top

16 楼Smile_Tiger（笑面虎）回复于 2001-02-27 16:42:00 得分 0

17 楼elaiyo（小牛头）回复于 2001-02-28 10:21:00 得分 0

楼上的，你有没有弄清楚我的意思？？
我的意思是如果i用float在计算机中表示为XXXX，那么i = * ( long * ) &y转换后i在计算集中还表示为XXXX，只不过要把它当作long型的数来理解。Top

18 楼lczddd（李找乐）回复于 2001-02-28 14:13:00 得分 0

唉！看来是要出人命了，我逃！Top

19 楼DaNiao（鸿雁）回复于 2001-02-28 19:39:00 得分 0

牛兄:
我不知道你对这些搞笑代码为什么这么有兴趣?
这种东西纯粹是一些人为了炫耀自己的一点点歪才搞出来吓人的
以前还有些吃饱饭没事干的人搞过一个模糊代码比赛,比的是谁的代码更难懂
对这些东西你看了一乐就完了,何必这么深究呢?

既然你问了我就花点时间说一下
#include <stdio.h>
main(_)//_是一个int变量,或者叫形参,其值由OS给出是命令行参数的个数(不同的OS稍有不同)
//在PDP-11的c里int型的声明可以省去int三个字
//在c里_被看作字母,与英文字母等效
{
//一个char型指针
//"=#6#694.21..#62.*.*2-)1#62.*..*.=#=1..29..#"是一个字符串常量
char*x="=#6#694.21..#62.*.*2-)1#62.*..*.=#=1..29..#";

while(_=*x/4)//把*x/4的值付给_(前面提到的那个变量)等于0时跳出
//*x是x所指字符串的第一个字

//循环体,逗号表达式,
//当_为0时!!_为0,否则!!_为1
//当_为非0时"\n%*s"+!!_相当于"%*s" (注意\n只是0X0A一个字符与DOS的回车不同)
//当_为0时还是"\n%*s"
//这里的%*s与scanf里的不同,在这里*表示后面紧接着的一个参数为宽度
//比如printf("%*d,3,x)就相当于printf("%3d,x)
//x的值被引用后再加1
//*x++%4*2就是x所指的字符的ASCII值对4取mod再乘2
//这决定了_/的个数
_-=8,printf("\n%*s"+!!_,_+_,"_/_/_/"+*x++%4*2);
}
他在一开始给出的那个字符串是刻意凑出来的

i = * ( long * ) &y;
i = 0x5f3759df - ( i >> 1 );
y = * ( float * ) &i;

i=(long)y是不能被编译器接受的,所以他用了上面的方法
i=y;所得的结果是不同的
i=y;是把y的值给i,也就是说如果y=3则i=3,y=0.1则i=0
这中间有复杂的转化(不仅仅是舍入的问题)
而i = * ( long * ) &y;则是把内存中y的内容(在内存里float的格式始祖鸟已经给出了)
直接搬到i的地址里面,结果是迥然不同的
他实际上是模拟了浮点运算
Top

20 楼DaNiao（鸿雁）回复于 2001-02-28 19:41:00 得分 0

累死我了Top

21 楼ed9er（始祖鸟）回复于 2001-02-28 23:03:00 得分 0

我想把问题再明了一点，少些于事无益的回应

//what the fuck那行其实已经进行了开方并求倒数的运算，可能是用了某些级数的公式或是什么其他的数学算式，也结合了float规范的定义，但是误差比较大，然后下面进行牛顿法迭代，关键的问题是这个融汇在一行里面的算法究竟是什么？是从什么数学式子里推导出来的？又为什么能表示成这么个形式？

把出处也说明一下，这是在id公布的源代码里面的，应该不是Carmark自己写的，可能是早期的hacker们搞出来的东西，我不知道对于一个3D引擎来说，这样的代码是不是搞笑，是不是某种炫耀，是不是吃饱了没事干

希望能有老大能帮我解决这个问题，我真的很想知道究竟是怎么回事Top

22 楼genesis（）回复于 2001-03-01 13:33:00 得分 0

可将常数换为0x5f3504f3,对比结果如下：
第一行为新常量的结果
第二行为原常量的结果
原常量对2^n的结果，较精确！其他稍差。
我只知道新常量的算法，不知是否有人感兴趣，感兴趣的请举手…… ：）
4:
0.498640
0.499154
9:
0.333238
0.332953
16:
0.249320
0.249577
25:
0.199872
0.199690
36:
0.166619
0.166477
49:
0.142852
0.142762
64:
0.124660
0.124788
81:
0.111110
0.111086
8.4752:
0.343468
0.343272
Top

23 楼Kevin_qing（）回复于 2001-03-01 13:37:00 得分 0

:)
快速的近似float运算.Top

24 楼ed9er（始祖鸟）回复于 2001-03-01 14:00:00 得分 0

知道你就说嘛，你不说我怎么能知道呢……Top

25 楼genesis（）回复于 2001-03-01 14:14:00 得分 0

…… ……
太麻烦了，自己一时写不完，明天贴！
不过，我想知道源程序从哪来的？Top

26 楼Kevin_qing（）回复于 2001-03-01 14:23:00 得分 0

/21
szn真像是苍蝇Top

27 楼ed9er（始祖鸟）回复于 2001-03-01 14:32:00 得分 0

quake源代码，id software
给个思路先~~~Top

28 楼genesis（）回复于 2001-03-01 17:50:00 得分 0

y = y * ( threehalfs - ( x2 * y * y ) )
为泰勒级数的表达式：
y=x**(-1/2)
y=f(x0)+f'(x0)*(x-x0)
=(x0)**(-1/2)-(1/2)*(x0)**(-3/2)*(x-x0)
=(x0)**(-1/2)-(1/2)*(x0)**(-3/2)*x + (1/2)*x0**(-1/2)
=1.5*(x0)**(-1/2)-(1/2)*x*(x0)**(-3/2)
若y0就非常接近(x0)**(-1/2),公式可迭代为
y=1.5*y - x/2 * y * y * y
即一开始的公式。

问题就剩下y0的求解：
If 0<E<255 then V=(-1)**S * 2 ** (E-127) * (1.F) where "1.F" is intended to represent the binary number created by prefixing F with an implicit leading 1 and a binary point.
If E=0 and F is nonzero, then V=(-1)**S * 2 ** (-126) * (0.F)=(-1)**S * 2 ** (E-127) * (G), G=(0.F)左移一位。
所以，只讨论0<=E<255,S=0时,x= 2 ** (E-127) * (G)
已知，y=x**(-1/2)=2**(63.5-E/2)*(G**(-1/2));
G**(-1/2)接近1,可以认为1.F中的F=0x0000;
x>>1=2**(E/2-127)*(G'),同理G'也可认为近似1；
常数C: 2 ** (E'-127) * (G"),
C-(x>>1)=2 ** (E'-E/2-127) * (Z),同理Z也可认为近似1;
所以，E'-E/2-127=63.5-E/2
E'=190.5,但E'为整数，E'=190,G"=2**(1/2);
所以可认为C: 2**(190-127)*2**(1/2),
C=0x5f3504f3。

不知，各位看官以为如何？
math.h内的函数绝大多数是用级数的。

这里只是抛砖引玉。Top

29 楼ed9er（始祖鸟）回复于 2001-03-01 19:23:00 得分 0

啪啦啦~~啪啦啦~~精彩啊

我看你的G"是估计的，不是求出来的，是吧？
你估的是1.414214
它的是1.432430
然后我看懂了你这个帖子后猜了个1.5，做出常数：5f400000
运行结果：
4 : 0.500000
9 : 0.330202
16 : 0.250000
25 : 0.198307
36 : 0.165101
49 : 0.141446
64 : 0.125000
81 : 0.110682
100 : 0.099153
啊……吐血，2^n的时候不要太精确哦~~

继续跟在砖头后面拣玉
Top

30 楼Kevin_qing（）回复于 2001-03-01 19:57:00 得分 0

精彩～
收藏～～～^_^Top

31 楼Smile_Tiger（笑面虎）回复于 2001-03-02 12:05:00 得分 0

to DaNiao(鸿雁):

你的 "i=(long)y是不能被编译器接受的,所以他用了上面的方法"的说法是错误的

首先 i=(long)y 是绝对能够编译通过的，这是ANSI C语言标准，意思是"将浮点数y的值转换成整型赋值给long型数i",它和 i = y 是等价的

其次，当然是你这句话的因果关系不成立啦
Top

32 楼genesis（）回复于 2001-03-02 13:44:00 得分 100

先兑现部分分数吧…… ：）

太忙了，说不了几句。
应该说G"不是估计值。
对Z的描述有点误差。
常数C: 2 ** (E'-127) * (G")，G"=1.F"
x>>1: 2 **(E/2 -127) * (G'),G'=1.F' or 0.F'
C - (x>>1) = 2 ** (E'-E/2-127) * (Z), Z=1.(F"-F')
因为E'-E/2-127为整数，对比x**(-1/2)的表达式
x**(-1/2)=2**(63.5-E/2)*(G**(-1/2))
因为G=1.F or 0.F,所以G**(-1/2)可认为为1。
x**(-1/2)=2**(63-E/2)*(2**0.5)
E'-E/2-127=63-E/2,
Z=2**0.5=1.(F"-F')
因为F'近似F>>1,如果忽略，则1.F"=2**0.5=G"
由于忽略了F',近似了G**(-1/2)，所以肯定有偏差，修正的方法要查书了，我想是统计方面的书，我没时间，
有兴趣自己解决吧。
知道误差的原因，可以找的误差最大的浮点数类型。
Top

33 楼Eros（）回复于 2001-03-02 16:06:00 得分 0

关注Top

34 楼ed9er（始祖鸟）回复于 2001-03-02 20:50:00 得分 0

差不多了，方法和思路肯定是对的了，再搞下去我也受不了了Top