"what the fuck ?"

[ed9er]

/*
** 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;
}


牛顿法求倒数平方，我自己看不懂，谁能讲解清楚？

[Smile_Tiger]

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

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

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

[Smile_Tiger]

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

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

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

[lczddd]

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

[DaNiao]

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

既然你问了我就花点时间说一下
#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的地址里面,结果是迥然不同的
他实际上是模拟了浮点运算

[DaNiao]

累死我了

[elaiyo]

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

[ed9er]

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

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

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

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

[genesis]

可将常数换为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

[Kevin_qing]

:)
快速的近似float运算.

[ed9er]

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

[genesis]

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

[Kevin_qing]

/21
szn真像是苍蝇

[ed9er]

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

[genesis]

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内的函数绝大多数是用级数的。


这里只是抛砖引玉。

[ed9er]

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

我看你的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的时候不要太精确哦~~

继续跟在砖头后面拣玉

[Kevin_qing]

精彩～
收藏～～～^_^

[genesis]

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

太忙了，说不了几句。
应该说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)，所以肯定有偏差，修正的方法要查书了，我想是统计方面的书，我没时间，
有兴趣自己解决吧。
知道误差的原因，可以找的误差最大的浮点数类型。

[Eros]

关注

[ed9er]

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

[Smile_Tiger]

to DaNiao(鸿雁):

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

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

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