新来在书上看到壹段运用多项式逼近计算平方根的代码,于今邑没拥有搞皓白干者是怎么铰算出产阿谁公式的。但在尝试处理效实的经过中,学到了不微少东方正西,于是便拥有了此雕刻篇心得,写出产到来和父亲家共享。就中拥有错漏的中,还请父亲家多多赐教养。

  确实,正如此多人所说的这么,当今拥有拥有FPU,拥有3DNow,拥有SIMD,讨论绵软件算法如同不符时宜。关于sqrt的话题实则早在2003年便已在 GameDev.net上违反掉落了普遍的讨论(却见我真实什分火星了,天然不扫摒除还拥有其他尚在冥王星的人,嘿嘿)。而尝诈寻求该话题则完整顿是出产于己己己的志趣和猎零数心(换句子话说坚硬是蒙昧)。

  我条是个beginner,因此此雕刻种父亲是父亲匪的效实我也说不清楚(在GameDev.net上也拥有很多相像的争议)。但无论何以,Carmack在DOOM3中还是运用了绵软件算法,而多知道壹点数学知对3D编程到来说也条要更加处没拥有变质处。3D图形编程实则坚硬是数学,数学,还是数学。

  =========================================================

  在3D图形编程中,日日要寻求平方根或平方根的倒腾数,比如:寻求向量的长度或将向量归壹募化。C数学函数库中的sqrt具拥有雄心的稀度,但关于3D游玩程式到来说快度太缓。我们期望却以在保障趾够的稀度的同时,进壹步提迅快度。

  Carmack在QUAKE3中运用了下面的算法,它第壹次在帮群场合出产即兴的时分,信直震住了所拥局部人。耳闻该算法实则并不是Carmack发皓的,它真正的干者是Nvidia的Gary Tarolli(不经证皓)。

  -----------------------------------

  //

  // 计算参数x的平方根的倒腾数

  //

  float InvSqrt (float x)

  {

  float xhalf=0.5f*x;

  int i=*(int*)&x;

  i=0x5f3759df – (i >> 1); // 计算第壹个相近根

  x=*(float*)&i;

  x=x*(1.5f – xhalf*x*x); // 牛顿迭代法

  return x;

  }

  ----------------------------------

  该算法的淡色实则坚硬是牛顿迭代法(Newton-Raphson Method,信称NR),而NR的基础则是泰勒级数(Taylor Series)。NR是壹种寻求方程的相近根的方法。比值先要估计壹个与方程的根比较接近的数值,然后根据公式铰算下壹个更其相近的数值,时时重骈直到却以得到满意的稀度。其公式如次:

  -----------------------------------

  函数:y=f(x)

  其壹阶带数为:y’=f'(x)

  则方程:f(x)=0 的第n+1个相近根为

  x[n+1]=x[n] – f(x[n]) / f'(x[n])

  -----------------------------------

  ?NR最关键的中在于估计第壹个相近根。假设该相近根与真根趾够接近的话,这么条需寻求微少半几次迭代,就却以违反掉落满意的松。

  当今回度过火到来看看何以使用牛顿法到来处理我们的效实。寻求平方根的倒腾数,还愿坚硬是寻求方程1/(x^2)-a=0的松。将该方程按牛顿迭代法的公式展开为:

  x[n+1]=1/2*x[n]*(3-a*x[n]*x[n])

  ?将1/2放到括号外面面,就违反掉落了下面阿谁函数的倒腾数第二行。

  接着,我们要想法估计第壹个相近根。此雕刻亦下面的函数最神物零数的中。它经度过某种方法算出产了壹个与真根什分接近的相近根,故此它条需寻求运用壹次迭代经过就得到了较满意的松。它是怎么做到的呢?所拥局部微妙就在于此雕刻壹行:

  i=0x5f3759df – (i >> 1); // 计算第壹个相近根

  超级莫皓其妙的语句子,不是吗?但细心想壹下的话,还是却以了松的。我们知道,IEEE规范下,float典型的数据在32位体系上是此雕刻么体即兴的(父亲体到来说坚硬是此雕刻么,但节微了很多底细,拥有志趣却以GOOGLE):

  -------------------------------

  bits:31 30 … 0

  31:标记位

  30-23:共8位,管指数(E)

  22-0:共23位,管条数(M)

  -------------------------------

  因此,32位的浮点数用什进制次数体即兴坚硬是:M*2^E。开根然后倒腾数坚硬是:M^(-1/2)*2^(-E/2)。当今就什鲜皓晰了。语句子i> >1其工干坚硬是将指数摒除以2,完成2^(E/2)的片断。而前面用壹个日数减去它,目的坚硬是违反掉落M^(1/2)同时反转所拥有指数的标记。

  到于阿谁0x5f3759df,呃,我不得不说,确实是壹个超级的Magic Number。

  阿谁Magic Number是却以铰带出产到来的,但我并不规划在此雕刻边讨论,鉴于真实太万端琐了。骈杂到来说,其规律如次:鉴于IEEE的浮点数中,条数M节微了最前面的1,因此还愿的条数是1+M。假设你在父亲学上数学课没拥有拥有打瞌睡的话,这么当你看到(1+M)^(-1/2)此雕刻么的方法时,应当会立雕刻联想的到它的泰勒级数展开,而该展开式的第壹项坚硬是日数。下面给出产骈杂的铰带经过:

  -------------------------------

  关于次数R>0,假定其在IEEE的浮点体即兴中,

  指数为E,条数为M,则:

  R^(-1/2)

  =(1+M)^(-1/2) * 2^(-E/2)

  将(1+M)^(-1/2)按泰勒级数展开,取第壹项,得:

  原式

  =(1-M/2) * 2^(-E/2)

  =2^(-E/2) – (M/2) * 2^(-E/2)

  假设不考虑指数的标记的话,

  (M/2)*2^(E/2)正是(R>>1),

  而在IEEE体即兴中,指数的标记条需骈杂地加以上壹个偏移即却,

  而式儿子的前半片断方好是个日数,因此原式却以转募化为:

  原式=C – (M/2)*2^(E/2)=C – (R>>1),就中C为日数

  因此条需寻求松方程:

  R^(-1/2)

  =(1+M)^(-1/2) * 2^(-E/2)

  =C – (R>>1)

  寻求出产令到对立误差最小的C值就却以了

  -------------------------------

  下面的铰带经过条是我团弄体的了松,并不违反掉落证皓。而Chris Lomont则在他的论文中详细讨论了最末阿谁方程的松法,并尝试在还愿的机具上寻摸最佳的日数C。拥有志趣的对象却以在文末了找到他的论文的链接。

  因此,所谓的Magic Number,并不是从N元宇宙的某个星系鉴于时空诬蔑而掉落到地球上的,而是几佰年前就拥局部数学即兴实。条需熟识NR和泰勒级数,你我异样拥有才干干出产相像的优募化。

  在GameDev.net上拥有人做度过测试,该函数的对立误差条约为0.177585%,快度比C规范库的sqrt提高皓度过20%。假设添加以壹次迭代经过,对立误差却以投降低到e-004 的级数,但快度也会投降到和sqrt差不多。耳闻在DOOM3中,Carmack经度过查找表进壹步优募化了该算法,稀度接近完备,同时快度也比原版提高了壹截(正竭力弄源码,谁拥有发我壹份)。

  犯得着剩意的是,在Chris Lomont的演算中,即兴实上最优秀的日数(稀度最高)是0x5f37642f,同时在还愿测试中,假设条运用壹次迭代的话,其效实亦最好的。但零数异的是,经度过两次NR后,在该日数下松的稀度将投降低得什分剧凶(天知道是怎么回事!)。经度过还愿的测试,Chris Lomont认为,最优秀的日数是0x5f375a86。假设换成64位的double版本的话,算法还是壹样的,而最优日数则为0x5fe6ec85e7de30da(又壹个令人冒汗的Magic Number – -b)。

  此雕刻个算法依顶赖于浮点数的外面部体即兴和字节以次,因此是不具移栽性的。假设放到Mac上跑就会挂掉落。假设想具拥有却移栽性,还是乖乖用sqrt好了。但算法思惟是畅通用的。父亲家却以尝试铰算壹下相应的平方根算法。

  下面给出产Carmack在QUAKE3中运用的平方根算法。Carmack曾经将QUAKE3的所拥有源代码捐给开源了,因此父亲家却以担心运用,不用担心会受到律师信。

  ---------------------------------

  //

  // Carmack在QUAKE3中运用的计算平方根的函数

  //

  float CarmSqrt(float x){

  union{

  int intPart;

  float floatPart;

  } convertor;

  union{

  int intPart;

  float floatPart;

  } convertor2;

  convertor.floatPart=x;

  convertor2.floatPart=x;

  convertor.intPart=0x1FBCF800 + (convertor.intPart >> 1);

  convertor2.intPart=0x5f3759df – (convertor2.intPart >> 1);

  return 0.5f*(convertor.floatPart + (x * convertor2.floatPart));

  }

  ---------------------------------------