C语言求平方和倒数

文章目录

- 1. 代码实现
- - float类型数据
  - double类型数据
  - 使用 double 类型的调整
- 2. 魔数与位级别操作
- - 浮点数表示
  - 位级别魔数操作
- 3. 牛顿迭代
- 4. 复杂代码具体解释
- - 具体解释：
  - 目的：
  - 举例：
- 5.感谢

平方和倒数 广泛用于计算机图形学中，尤其是在处理大规模3D渲染时。这种方法的核心思想是利用了位级别的操作和牛顿迭代法，以非常高效的方式计算平方根的倒数。

平方根倒数（Reciprocal Square Root）计算的问题是，对于一个给定的数 n，我们希望快速找到 1/√n。在图形学和物理模拟中，这种计算非常常见，通常用于归一化向量。

传统的平方根计算通常使用浮点运算，较为耗时。为了提高效率，特别是在早期硬件性能有限的情况下，出现了一种基于“魔数”（magic number）的快速算法。

1. 代码实现

float类型数据

float Q_rsqrt(float n)
{long i;float x2, y;const float threehalfs = 1.5F;x2 = n * 0.5F;y = n;i = *(long *)&y;i = 0x5f3759df - (i >> 1);y = *(float *)&i;// First iterationy = y * (threehalfs - (x2 * y * y));// Second iteration for higher precisiony = y * (threehalfs - (x2 * y * y));// Third iteration for higher precisiony = y * (threehalfs - (x2 * y * y));// Fourth iteration for higher precisiony = y * (threehalfs - (x2 * y * y));return y;
}

double类型数据

double Q_rsqrt(double n)
{long long i;double x2, y;const double threehalfs = 1.5;x2 = n * 0.5;y = n;i = *(long long *)&y;i = 0x5fe6eb50c7b537a9 - (i >> 1); // 使用64位常量y = *(double *)&i;// First iterationy = y * (threehalfs - (x2 * y * y));// Second iteration for higher precisiony = y * (threehalfs - (x2 * y * y));// Third iteration for higher precisiony = y * (threehalfs - (x2 * y * y));// Fourth iteration for higher precisiony = y * (threehalfs - (x2 * y * y));return y;
}

使用 double 类型的调整

在处理 double 类型时，需要进行一些调整：

long 类型需要扩展为 long long 以适应 64 位浮点数的表示。
魔数 0x5f3759df 替换为适用于 64 位的 0x5fe6eb50c7b537a9。

2. 魔数与位级别操作

这个方法的核心是通过使用魔数对浮点数进行位级别的操作，得到平方根倒数的一个初始近似值。

浮点数表示

在现代计算机中，浮点数通常采用 IEEE 754 标准表示，包括符号位、指数部分和尾数部分。对于 32 位浮点数来说，结构如下：

1 位符号位
8 位指数部分
23 位尾数部分

位级别魔数操作

Q_rsqrt 方法的精妙之处在于，它将浮点数视作一个整数，通过操作位来得到平方根倒数的近似值。核心操作如下：

i = 0x5f3759df - (i >> 1);

在这个公式中，i 是浮点数 y 在内存中的位表示，0x5f3759df 是一个特定的常量，这个常量通过经验或实验计算得出，能够将位级别的操作映射到一个相对合理的平方根倒数的近似值。对于 double 类型，需要将这个常量替换为 64 位的版本，即 0x5fe6eb50c7b537a9。

这个操作的本质是近似计算了数值的对数，再通过魔数进行校正，使得结果在位操作的基础上接近实际的平方根倒数值。

3. 牛顿迭代

初始近似值得到后，使用牛顿迭代法进一步提高精度。牛顿迭代法是一种经典的数值方法，用于求解方程的根。对于平方根倒数问题，目标是找到 y，使得 y = 1/√n。

牛顿迭代的公式如下：

y = y * (threehalfs - (x2 * y * y));

其中：

x2 = n * 0.5，是数值的一半，用于简化计算。
threehalfs = 1.5，是一个常量，用于调整收敛速度和稳定性。

这个公式利用当前的 y 值，结合原数 n，逐步逼近实际的 1/√n。每次迭代后，y 的值会变得更接近真实的平方根倒数值。

为了提高精度，通常会进行多次迭代。在最常见的实现中，进行一次或两次迭代已经可以满足大部分应用场景中的精度要求。

4. 复杂代码具体解释

代码 i = *(long *)&y;和 y = *(float *)&i; 是一种类型惩断（type punning）的操作方式，用于在不改变底层数据的情况下，以另一种类型访问相同的内存位置。以i = *(long *)&y;为例：

具体解释：

y 是一个 float 类型的变量。
&y 获取的是变量 y 的内存地址，这个地址是指向 float 类型数据的指针。
(long *)&y 将这个 float 类型的指针强制转换为 long 类型的指针。注意，这个操作并不改变 y 的实际数据，只是告诉编译器，“虽然 &y 原本是 float* 类型，但我要以 long* 类型来看待它。”
*(long *)&y 是取出这个 long* 指针所指向的数据。这意味着它会将 y 所在内存位置的位模式当作一个 long 类型的整数来读取。