函数

欧拉的一篇研究报告：关于整数因子和的一个非常奇特规律的发现

/ 12 条评论

    《数学与猜想》里引用了一段欧拉的这篇经典的研究报告，写的非常精彩。你可以从中看到一个数学家是如何进行发现、归纳、猜想和论证的。你可以看到两个完全不同的数学模型里出现了惊人的巧合，通过挖掘它们之间的内在联系，最终完成了伟大的统一。
    没扫描仪，拿相机拍的，效果非常不好，见谅了！
    另外，拜托大家不要盗链下面的图片。

  
  
  
  
  

旧闻一则：神秘的0x5f3759df 不可思议的Quake III源码

/ 35 条评论

    Quake III公开源码后，有人在game/code/q_math.c里发现了这样一段代码。它的作用是将一个数开平方并取倒，经测试这段代码比(float)(1.0/sqrt(x))快4倍：
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

  #ifndef Q3_VM
  #ifdef __linux__
    assert( !isnan(y) ); // bk010122 - FPE?
  #endif
  #endif
  return y;
}

    code/common/cm_trace.c中也出现了这样一段解释sqrt(x)的函数，与上面的代码唯一不同的就是这个函数返回的是number*y：
/*
================
SquareRootFloat
================
*/
float SquareRootFloat(float number) {
    long i;
    float x, y;
    const float f = 1.5F;

    x = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;
    i  = 0x5f3759df - ( i >> 1 );
    y  = * ( float * ) &i;
    y  = y * ( f - ( x * y * y ) );
    y  = y * ( f - ( x * y * y ) );
    return number * y;
}

    这样的代码速度肯定飞快，我就不用多说了；但算法的原理是什么呢？其实说穿了也不是很神，程序首先猜测了一个接近1/sqrt(number)的值，然后两次使用牛顿迭代法进行迭代。根号a的倒数实际上就是方程1/x^2 – a = 0的一个正实根，它的导数是-2/x^3。运用牛顿迭代公式x' = x – f(x)/f'(x)，式子化简为x' = x * (1.5 – 0.5a * x^2)。迭代几次后，x的值将趋于1/sqrt(a)。
    但这段代码真正牛B的是那个神秘的0x5f3759df，因为0x5f3759df – (i >> 1)出人意料地接近根号y的倒数。人们都不知道这个神秘的常数是怎么来的，只能把它当作神来膜拜。这个富有传奇色彩的常数到底咋回事，很少有人说得清楚。这篇论文比较科学地解释了这个常数。

牛顿迭代法快速寻找平方根

/ 37 条评论

    下面这种方法可以很有效地求出根号a的近似值：首先随便猜一个近似值x，然后不断令x等于x和a/x的平均数，迭代个六七次后x的值就已经相当精确了。
    例如，我想求根号2等于多少。假如我猜测的结果为4，虽然错的离谱，但你可以看到使用牛顿迭代法后这个值很快就趋近于根号2了：

(       4  + 2/   4     ) / 2 = 2.25
(    2.25  + 2/   2.25  ) / 2 = 1.56944..
( 1.56944..+ 2/1.56944..) / 2 = 1.42189..
( 1.42189..+ 2/1.42189..) / 2 = 1.41423..
….

      
    这种算法的原理很简单，我们仅仅是不断用(x,f(x))的切线来逼近方程x^2-a=0的根。根号a实际上就是x^2-a=0的一个正实根，这个函数的导数是2x。也就是说，函数上任一点(x,f(x))处的切线斜率是2x。那么，x-f(x)/(2x)就是一个比x更接近的近似值。代入f(x)=x^2-a得到x-(x^2-a)/(2x)，也就是(x+a/x)/2。
    同样的方法可以用在其它的近似值计算中。Quake III的源码中有一段非常牛B的开方取倒函数。