数学库函数的加速实现

吃点汉堡
4,458 阅读19分钟

基于C语言数学库函数的加速实现

查表法的基本原理

查表法（Lookup Table, LUT）是一种通过预先计算和存储函数值，在需要时直接查找表中的值以实现快速计算的方法。这种方法广泛应用于计算密集型任务中，如数学函数计算、图像处理、信号处理等。查表法的核心思想是用空间换取时间，通过预先计算和存储结果，减少实时计算的负担。

原理详述

  1. 预计算和存储：

    • 将一个复杂函数的值在某个范围内进行预先计算。
    • 将这些计算结果存储在一个数组（查找表）中，数组的索引对应函数的输入值。

  2. 查找和插值：

    • 当需要计算函数值时，通过查找表直接获取预先计算好的值。
    • 如果输入值不在查找表的索引上，可以使用插值方法（如线性插值）在邻近的表值之间进行近似计算。

数学函数查表法

以计算三角函数 sin⁡(x)\sin(x) 为例，使用查表法的过程如下：

  1. 构建查找表：

    • 在 [0,2π][0, 2\pi] 范围内，选择等间隔的点，如 0,2πN,4πN,…,2π0, \frac{2\pi}{N}, \frac{4\pi}{N}, \dots, 2\pi，其中 NN 是分割的区间数。
    • 计算这些点上的 sin⁡\sin 值，并存储在数组中，如 sinTable[i]=sin⁡(2πiN)\text{sinTable}[i] = \sin\left(\frac{2\pi i}{N}\right)。

  2. 查找和插值：

    • 对于任意输入 xx，首先将 xx 归一化到 [0,2π][0, 2\pi] 范围内（利用周期性），然后找到最接近的两个表值 xix_i 和 xi+1x_{i+1}。

    • 使用这些表值进行插值计算 sin⁡(x)\sin(x) 的近似值，如线性插值：

      sin⁡(x)≈sin⁡(xi)+sin⁡(xi+1)−sin⁡(xi)xi+1−xi(x−xi)\sin(x) \approx \sin(x_i) + \frac{\sin(x_{i+1}) - \sin(x_i)}{x_{i+1} - x_i} (x - x_i)

优缺点

  • 优点：

    • 快速：通过查找和简单的插值计算，极大地减少了复杂函数的计算时间。
    • 简单：实现简单，只需存储和查找，不需要复杂的数学运算。

  • 缺点：

    • 存储空间：需要预先存储大量的计算结果，消耗较多的内存。
    • 精度受限：受表格分辨率限制，插值方法可能带来误差。

实例：三角函数查表法实现

以 cosf 函数为例，使用查表法计算余弦值：

float cosf(float x)
{
	  int K = (int)(x * TABLE_SIZEDivTwoPi);
	  x = x * TABLE_SIZEDivTwoPi - (float)K;
	  x *= TwoPiDivTABLE_SIZE;
	  K &= TABLE_MASK;
	  float SinK = sincosTable[K];
	  float CosK = sincosTable[K + 32];
	  float x2 = x * x, x3 = x2 * x, x4 = x2 * x2, x5 = x4 * x;
	  float cos = CosK - x * SinK + x2 * Coef0 * CosK + x3 * Coef1_pos * SinK + x4 * Coef2 * CosK + x5 * Coef3_neg * SinK;
	  return cos;
}

解释：

  1. 归一化和索引计算：

    • 将输入 xx 转换为对应的查表索引 KK，并将 xx 归一化到一个小范围内。

  2. 查找表值：

    • 根据索引 KK 从查找表中获取对应的 sin⁡\sin 和 cos⁡\cos 值。

  3. 多项式展开：

    • 使用多项式展开方法，通过查找的表值和预先定义的多项式系数进行计算，得到近似的余弦值。

通过查表法，可以大大提高函数计算的效率，特别是在需要频繁计算相同函数的大量输入时，查表法的优势更加明显。

多项式展开与级数展开和快速平方根倒数算法的原理

1. 多项式展开

多项式展开（Polynomial Expansion）是一种用多项式近似函数的方法，通过有限项的多项式来逼近复杂函数，从而提高计算效率。

泰勒展开：

  • 泰勒展开是多项式展开的一种常见形式，它将函数在某一点附近展开成多项式形式：

    f(x)=f(a)+f′(a)(x−a)+f′′(a)2!(x−a)2+f′′′(a)3!(x−a)3+⋯f(x) = f(a) + f'(a)(x-a) + \frac{f''(a)}{2!}(x-a)^2 + \frac{f'''(a)}{3!}(x-a)^3 + \cdots

  • 通过截断高阶项，可以得到函数的近似值。

举例：余弦函数的泰勒展开：

cos⁡(x)≈1−x22!+x44!−x66!\cos(x) \approx 1 - \frac{x^2}{2!} + \frac{x^4}{4!} - \frac{x^6}{6!}

在计算机实现中，可以通过预定义系数和幂次进行多项式计算，从而快速得到近似值。

float cosf(float x) {
    float x2 = x * x;
    float x4 = x2 * x2;
    return 1 - (x2 / 2) + (x4 / 24);
}

通过有限项的多项式展开，可以在不牺牲太多精度的情况下，显著提高计算速度。

2. 级数展开

级数展开（Series Expansion）是另一种逼近函数的方法，它通过无穷级数的和来表示函数。在实际计算中，通常截断到有限项以达到足够的精度。

麦克劳林级数：

  • 麦克劳林级数是泰勒级数在 a=0a = 0 点的展开：

    f(x)=f(0)+f′(0)x+f′′(0)2!x2+f′′′(0)3!x3+⋯f(x) = f(0) + f'(0)x + \frac{f''(0)}{2!}x^2 + \frac{f'''(0)}{3!}x^3 + \cdots

举例：双曲正弦函数的级数展开：

sinh⁡(x)=x+x33!+x55!+⋯\sinh(x) = x + \frac{x^3}{3!} + \frac{x^5}{5!} + \cdots

通过计算有限项，可以近似计算双曲正弦值。

实现：

float sinhf(float x) {
    float sum = x;
    float term = x;
    for (int n = 1; term > 0.000001; n++) {
        term *= (x * x) / ((2 * n) * (2 * n + 1));
        sum += term;
    }
    return sum;
}

通过逐项累加，使用有限项的级数展开可以在保证精度的同时提高计算效率。

3. 快速平方根倒数算法

快速平方根倒数算法（Fast Inverse Square Root, Q_rsqrt）是一种高效计算浮点数平方根倒数的方法。最著名的实现是用于计算机图形学中的 Q_rsqrt 算法。

原理：

  • 利用浮点数的位级表示，通过魔数和位操作快速得到平方根倒数的近似值，然后通过牛顿迭代法进行修正。

步骤：

  1. 通过位操作取得浮点数的初始近似值：

    i=magic−(i>>1)i = \text{magic} - (i >> 1)

    其中，magic 是一个特定的常数，对于 32 位浮点数，该常数通常为 0x5f3759df 或 0x5f375a86。

  2. 将位操作后的值转换回浮点数，并进行一到两次牛顿迭代，以提高精度：

    y=y×(1.5−0.5×x×y2)y = y \times (1.5 - 0.5 \times x \times y^2)

实现：

float isqrtf(float x) {
    float xhalf = 0.5f * x;
    int i = *(int*)&x;
    i = 0x5f375a86 - (i >> 1);
    x = *(float*)&i;
    x = x * (1.5f - xhalf * x * x); // 第一次迭代
    x = x * (1.5f - xhalf * x * x); // 第二次迭代
    return x;
}

float sqrtf(float x) {
    return x * isqrtf(x);
}

通过这个算法，可以快速计算平方根的倒数，然后再通过一次乘法得到平方根。

原理总结

提高计算效率和精度的原因

  • 预计算和查表：通过预先计算和存储复杂函数值，查表法可以在运行时快速查找和插值，减少了实时计算的复杂度。
  • 多项式和级数展开：这些方法通过有限项的多项式或级数近似函数，避免了复杂的函数计算，大大提高了计算速度。在保证足够精度的情况下，截断高阶项可以显著减少计算量。
  • 快速平方根倒数算法：利用浮点数的位级表示，通过简单的位操作和少量迭代，实现了平方根倒数的高效计算。这个算法利用了计算机硬件的高效整数运算，极大地提高了计算速度。

通过以上方法，可以在计算复杂数学函数时，显著提高效率和精度，满足高性能计算的需求。

具体实现

正弦函数 (sinf)

float sinf(float x)
{
	  int K = (int)(x * TABLE_SIZEDivTwoPi);
	  x = x * TABLE_SIZEDivTwoPi - (float)K;
	  x *= TwoPiDivTABLE_SIZE;
	  K &= TABLE_MASK;
	  K += 32;
	  float CosK = sincosTable[K];
	  float SinK = sincosTable[K - 32];
	  float x2 = x * x, x3 = x2 * x, x4 = x2 * x2, x5 = x4 * x;
	  float sin = SinK + x * CosK + x2 * Coef0 * SinK + x3 * Coef1 * CosK + x4 * Coef2 * SinK + x5 * Coef3 * CosK;
	  return sin;
}

数学原理：

  • sinf 函数计算给定值的正弦值，即 sin⁡(x)\sin(x)。

  • 算法步骤：

    1. 通过查表法获取近似的正弦和余弦值。
    2. 使用多项式展开计算正弦值。

实现细节：

  • 使用查表法获取近似的正弦和余弦值。
  • 通过多项式展开进行精确计算，使用预定义的系数进行计算。

双曲正弦函数 (sinhf)

float sinhf(float x)
{
	float fVal = x;
	float Xn = x, temp = 1.0, den = 6.0;
	float count = 3.0;
	float sum = x;
	while(fabsf(temp) > TINY_VALUE)
	{
		Xn = Xn * fVal * fVal;
		temp  = Xn / den;
		sum += temp;
		count = count + 2.0f;
		den  = den * (count - 1.0f) * count;
	}
	return sum;
}

数学原理：

  • sinhf 函数计算给定值的双曲正弦值，即 sinh⁡(x)\sinh(x)。
  • 使用级数展开计算双曲正弦值： sinh⁡(x)=∑n=0∞x2n+1(2n+1)!\sinh(x) = \sum_{n=0}^{\infty} \frac{x^{2n+1}}{(2n+1)!}

实现细节：

  • 使用循环和累加计算级数展开的各项，直到达到精度要求（TINY_VALUE）。
  • 每次循环计算当前项的分子和分母，并将当前项加到和中。

余弦函数 (cosf)

float cosf(float x)
{
	  int K = (int)(x * TABLE_SIZEDivTwoPi);
	  x = x * TABLE_SIZEDivTwoPi - (float)K;
	  x *= TwoPiDivTABLE_SIZE;
	  K &= TABLE_MASK;
	  float SinK = sincosTable[K];
	  float CosK = sincosTable[K + 32];
	  float x2 = x * x, x3 = x2 * x, x4 = x2 * x2, x5 = x4 * x;
	  float cos = CosK - x * SinK + x2 * Coef0 * CosK + x3 * Coef1_pos * SinK + x4 * Coef2 * CosK + x5 * Coef3_neg * SinK;
	  return cos;
}

数学原理：

  • cosf 函数计算给定值的余弦值，即 cos⁡(x)\cos(x)。
  • 通过查表法获取近似值，并结合多项式展开进行精确计算。

实现细节：

  • 使用查表法获取近似的正弦和余弦值。
  • 通过多项式展开进行精确计算，使用预定义的系数进行计算。

双曲余弦函数 (coshf)

float coshf(float x)
{
	float fVal = x;
	float Xn = 1.0, temp = 1.0, den = 2.0;
	float count = 2.0;
	float sum = 1.0;
	while(fabsf(temp) > TINY_VALUE)
	{
		Xn = Xn * fVal * fVal;
		temp  = Xn / den;
		sum += temp;
		count = count + 2.0f;
		den  = den * (count - 1.0f) * count;
	}
	return sum;
}

数学原理：

  • coshf 函数计算给定值的双曲余弦值，即 cosh⁡(x)\cosh(x)。
  • 使用级数展开方法计算双曲余弦值，即 cosh⁡(x)=∑n=0∞x2n(2n)!\cosh(x) = \sum_{n=0}^{\infty} \frac{x^{2n}}{(2n)!}。

实现细节：

  • 使用循环和累加计算级数展开的各项，直到达到精度要求。
  • 通过条件判断确保计算精度。

反余弦函数 (acosf)

float acosf(float x)
{
	float input = x;
	if (x < 0.0f)
	{
		input = -x;
	}
	float m = input * input;
	m = sqrtf(1 - m);
	float acos_v = m / input;
	float a, b;
	a = atanf(acos_v) / (2.0f * pi);
	b = atanf(1.0f / acos_v) / (2.0f * pi);
	if (acos_v >= 1.0f)
		a = 0.125f;
	else if (acos_v <= -1.0f)
		a = -0.125f;
	if (1.0f / acos_v >= 1.0f)
		b = 0.125f;
	else if (1.0f / acos_v <= -1.0f)
		b = -0.125f;
	a = 2.0f * pi * a;
	b = 2.0f * pi * b;
	b = pi / 2.0f - b;
	if (acos_v > 1.0f)
		a = b;
	if (x < 0.0f)
		a = pi - b;
	return a;
}

数学原理：

  • acosf 函数计算给定值的反余弦，即 cos⁡−1(x)\cos^{-1}(x)，返回值的范围是 [0, π]。
  • 首先将输入值取绝对值，如果输入值小于0，则将其变为正值。
  • 计算 m=1−x2m = \sqrt{1 - x^2}，即反余弦的对称轴部分。
  • 使用反正切函数 atanf 来计算角度，并进行适当的缩放和调整。
  • 通过条件判断调整结果，以确保结果在正确的范围内。

实现细节：

  • sqrtf 函数用于计算平方根，atanf 函数用于计算反正切。
  • 2.0f * pi 和 pi / 2.0f 用于将弧度转换为角度。
  • 最后，通过条件判断调整角度结果，以确保结果的准确性和范围。

反正弦函数 (asinf)

float asinf(float X)
{
	float input = X;
	if (X < 0.0f) {
		input = -X;
	}
	float x = input * input;
	x = sqrtf(1.0f - x);
	float acos_v = input / x;
	float R0H, R3H;
	R0H = atanf(acos_v) / (2.0f * pi);
	R3H = atanf(1.0f / acos_v) / (2.0f * pi);
	if (acos_v >= 1.0f)
		R0H = 0.125f;
	else if (acos_v <= -1.0f)
	    R0H = -0.125f;
	if (1 / acos_v >= 1.0f)
	    R3H = 0.125f;
	else if (1.0f / acos_v <= -1.0f)
	    R3H = -0.125f;
	R0H = 2.0f * pi * R0H;
	R3H = 2.0f * pi * R3H;
	R3H = pi / 2.0f - R3H;
	if (acos_v > 1.0f)
	    R0H = R3H;
	if (X < 0.0f)
	    R0H = pi - R3H;
	return R0H;
}

数学原理：

  • asinf 函数计算给定值的反正弦，即 sin⁡−1(x)\sin^{-1}(x)，返回值的范围是 [-π/2, π/2]。
  • 类似于 acosf，首先取绝对值，然后计算 1−x2\sqrt{1 - x^2}。
  • 通过计算反余弦的对称轴部分，再使用反正切函数来计算角度。
  • 根据输入值的符号和范围进行调整，以确保结果在正确的范围内。

实现细节：

  • 与 acosf 类似，使用 sqrtf 和 atanf 计算，并通过常量调整结果。
  • 最后，通过条件判断调整角度结果，以确保准确性和范围。

反正切函数 (atanf)

float atanf(float x)
{
	float input = x;
	float Numerator, Denominator;
	if (x < 0.0f)
		input = -x;
	if (1.0f >= input)
	{
		Numerator = input;
		Denominator = 1.0f;
	}
	else
	{
		Numerator = 1.0f;
		Denominator = input;
	}
	float Ratio = Numerator / Denominator;
	int index = (int)(Ratio * 64.0f);
	float A2 = atan2Table[3 * index + 2];
	float A1 = atan2Table[3 * index + 1];
	float A0 = atan2Table[3 * index + 0];
	float atan2 = A0 + Ratio * (A1 + A2 * Ratio);

	float angle = 0.0f;
	if (x >= 0.0f && 1.0f >= input)
		angle = atan2;
	else if (x >= 0.0f && 1.0f < input)
		angle = HalfPI - atan2;
	else if (x < 0.0f)
		angle = -atan2;
	return angle;
}

数学原理：

  • atanf 函数计算给定值的反正切，即 tan⁡−1(x)\tan^{-1}(x)，返回值的范围是 [-π/2, π/2]。
  • 通过取绝对值和比较，确定使用的分子和分母。
  • 计算比例 Ratio，然后通过查表法获取近似值，并进行插值计算。
  • 最后，根据输入值的符号和范围调整结果。

实现细节：

  • 使用 atan2Table 查表获取近似值，并通过插值计算得到结果。
  • 根据输入值的范围和符号调整结果，以确保准确性。

反正切函数（双参数）(atan2f)

float atan2f(float y, float x)
{
	float input1 = y, input2 = x;
	float Numerator, Denominator;
	if (y < 0.0f)
		input1 = -y;
	if (x < 0.0f)
		input2 = -x;
	if (input1 >= input2)
	{
		Numerator = input2;
		Denominator = input1;
	}
	else
	{
		Numerator = input1;
		Denominator = input2;
	}
	float Ratio = Numerator / Denominator;
	int index = (int)(Ratio * 64.0f);
	float A2 = atan2Table[3 * index + 2];
	float A1 = atan2Table[3 * index + 1];
	float A0 = atan2Table[3 * index + 0];
	float temp = A0 + Ratio * (A1 + A2 * Ratio);
	float angle = 0.0f;
	if (input2 >= input1)
		angle = temp;
	else if (input2 < input1)
		angle = HalfPI - temp;
	float flag = x * y;
	if(flag < 0.0f)
		angle = -angle;
	if(y >= 0.0f && x < 0.0f)
		angle = angle + pi;
	else if(y < 0.0f && x < 0.0f)
		angle = angle - pi;
	else if(y > 0.0f && x == 0.0f)
		angle = HalfPI;
	else if(y < 0.0f && x == 0.0f)
		angle = -HalfPI;
	if(x == 0.0f)
	{
		if(y == 0.0f)
			angle = 0.0f;
		else if(y == -0.0f)
			angle = -0.0f;
	}
	else if(x == -0.0f)
	{
		if(y == 0.0f)
			angle = pi;
		else if(y == -0.0f)
			angle = -pi;
	}
	return angle;
}

数学原理：

  • atan2f 函数计算两个浮点数的反正切，即 tan⁡−1(x/y)\tan^{-1}(x / y)，返回值的范围是 [-π, π]。
  • 处理特殊情况，如 y 为0时的情况。
  • 使用 atanf 计算反正切，并根据 y 的符号调整结果。

实现细节：

  • 通过条件判断处理特殊情况。
  • 使用 atanf 计算角度，并根据 y 的符号调整结果。

向上取整函数 (ceilf)

float ceilf(float x)
{
	int y = x;
	if(x > 0.0f && ((float)y - x) != 0)
		return ((float)y + 1.0f);
	return (float)y;
}

数学原理：

  • ceilf 函数返回大于或等于给定浮点数的最小整数值。
  • 通过取整后检查是否需要向上调整，以确保结果是向上取整。

实现细节：

  • 使用强制类型转换将浮点数转换为整数。
  • 通过条件判断确定是否需要向上调整结果。

指数函数 (expf)

float expf(float x)
{
	float fVal = x;
	if(x < 0.0f)
		fVal = -fVal;
	int xm = (int)fVal;
	fVal = fVal - (float)xm;
	float en = 1.0;
	if(xm < 89)
		en = ExpTable[xm];
	else
	{
		int num = xm / 88;
		int man = xm - num * 88;
		for(int i = 0; i < num; i++)
			en *= ExpTable[88];
		en *= ExpTable[man];
	}
	float Xn = 1.0, temp = 1.0, den = 1.0;
	float count = 1.0;
	float sum = 1.0;
	while(temp > 0.000001f)
	{
		Xn = Xn * fVal;
		temp  = Xn / den;
		sum += temp;
		count = count + 1.0f;
		den *= count;
	}
	if(x < 0.0f)
		return 1.0f / (sum * en);
	return sum * en;
}

数学原理：

  • expf 函数计算给定值的指数值，即 exe^x。
  • 使用查表法获取整数部分的近似值，结合级数展开计算小数部分。

实现细节：

  • 使用 ExpTable 查表获取近似值，并结合循环和累加计算级数展开的各项。
  • 通过条件判断确保计算精度。

绝对值函数 (fabsf)

float fabsf(float x)
{
	return x < 0.0f ? -x : x;
}

数学原理：

  • fabsf 函数返回给定浮点数的绝对值。
  • 通过条件判断返回正值。

实现细节：

  • 使用简单的条件判断确定返回值。

向下取整函数 (floorf)

float floorf(float x)
{
	int y = x;
	if(x < 0.0f && ((float)y - x) != 0)
		return ((float)y - 1.0f);
	return (float)y;
}

数学原理：

  • floorf 函数返回小于或等于给定浮点数的最大整数值。
  • 通过取整后检查是否需要向下调整，以确保结果是向下取整。

实现细节：

  • 使用强制类型转换将浮点数转换为整数。
  • 通过条件判断确定是否需要向下调整结果。

浮点数取模函数 (fmodf)

float fmodf(float x, float y){
    float ind = x / y;
	float inter = (ind > 0.0f) ? floorf(ind) : ceilf(ind);
    float res = 0.0;
    res = x - inter * y;
    return res;
}

数学原理：

  • fmodf 函数返回浮点数除法的余数，即 xmod  yx \mod y。

  • 算法步骤：

    1. 计算 ind=xyind = \frac{x}{y}。
    2. 根据 indind 的符号选择向下取整或向上取整，得到整数部分 interinter。
    3. 计算 res=x−inter⋅yres = x - inter \cdot y。

实现细节：

  • 使用 floorf 或 ceilf 进行取整。
  • 减去整数部分的乘积，得到余数。

浮点数拆分函数 (frexpf)

float frexpf(float x, int *y)
{
    float fVal = x;
    if(x < 1.0f && x > -1.0f && x != 0)
        fVal =  1.0f / x;
    int intx = fVal;
    if(fVal < 0.0f)
        intx = -intx;
    while(intx > 0.0f)
    {
        *y = *y + 1;
        intx >>= 1;
    }
    float man;
    if(x < 1.0f && x > -1.0f && x != 0)
    {
        *y = -(*y) + 1;
        man = x * (1 << (-(*y)));
    }
    else
        man = x / (1 << *y);
    return man;
}

数学原理：

  • frexpf 函数将浮点数分解为尾数和指数，使得 x=man×2yx = \text{man} \times 2^y。

  • 算法步骤：

    1. 处理 xx 的绝对值在 (0, 1) 范围内的特殊情况。
    2. 使用位移操作计算指数部分 yy。
    3. 计算尾数部分 man\text{man}。

实现细节：

  • 使用位移操作进行计算。
  • 处理特殊情况，通过调整指数和尾数使得结果正确。

计算直角三角形斜边长 (hypotf)

float hypotf(float x,float y)
{
	float temp=x*x+y*y;
	float res=sqrtf(temp);
	return res;
}

数学原理：

  • hypotf 函数计算直角三角形的斜边长，即 x2+y2\sqrt{x^2 + y^2}。

  • 算法步骤：

    1. 计算 x2+y2x^2 + y^2。
    2. 使用平方根函数计算结果。

实现细节：

  • 直接使用平方和的平方根计算斜边长。

浮点数乘以2的幂 (ldexpf)

float ldexpf(float x, int p){
    int fval = p;
    if(p < 0)
    	fval = -p;
    int exp = 1;
    exp <<= fval;
    float fexp = (float)exp;
    if(p < 0)
    	fexp = 1.0f / fexp;
    return x * fexp;
}

数学原理：

  • ldexpf 函数返回 x×2px \times 2^p 的值。

  • 算法步骤：

    1. 计算 2p2^p，并处理负指数情况。
    2. 返回 xx 乘以 2p2^p 的结果。

实现细节：

  • 使用位移操作计算 2p2^p。
  • 处理负指数，通过取倒数计算。

自然对数函数 (logf)

float logf(float a)
{
	float input = a;
	if (a < 0.0f)
		input = -a;
	int *ptr = (int *)&input;
	int exponent = *ptr >> 23;
	float f_exp = (float)exponent;
	int value = *ptr & LN_TABLE_MASK1;
	value |= LN_TABLE_MASK2;
	float *va_ptr = (float *)&value;
	float f_value = *va_ptr;
	float f_man = f_value - (float)((int)f_value);
	float ret1 = (f_exp - BIAS) * LNV2;
	int index = (int)(32.0f * f_man);
	float A2 = LnTable[3 * index + 2];
	float A1 = LnTable[3 * index + 1];
	float A0 = LnTable[3 * index + 0];
	float ret2 = A0 + f_man * (A1 + A2 * f_man);
	return ret1 + ret2;
}

数学原理：

  • logf 函数计算给定值的自然对数，即 ln⁡(a)\ln(a)。

  • 算法步骤：

    1. 将浮点数 aa 转换为整数位表示，提取指数部分和尾数部分。
    2. 通过查表法计算对数值。

实现细节：

  • 使用位操作提取浮点数的指数部分和尾数部分。
  • 通过查表法获取对数近似值，并进行插值计算。

常用对数函数 (log10f)

float log10f(float a)
{
	float input = a;
	if (a < 0.0f)
		input = -a;
	int *ptr = (int *)&input;
	int exponent = *ptr >> 23;
	float f_exp = (float)exponent;
	int value = *ptr & LN_TABLE_MASK1;
	value |= LN_TABLE_MASK2;
	float *va_ptr = (float *)&value;
	float f_value = *va_ptr;
	float f_man = f_value - (float)((int)f_value);
	float ret1 = (f_exp - BIAS) * LNV2;
	int index = (int)(32.0f * f_man);
	float A2 = LnTable[3 * index + 2];
	float A1 = LnTable[3 * index + 1];
	float A0 = LnTable[3 * index + 0];
	float ret2 = A0 + f_man * (A1 + A2 * f_man);
	return (ret1 + ret2)/LOG10;
}

数学原理：

  • log10f 函数计算给定值的常用对数，即 log⁡10(a)\log_{10}(a)。

  • 算法步骤：

    1. 将浮点数 aa 转换为整数位表示，提取指数部分和尾数部分。
    2. 通过查表法计算自然对数值，并转换为常用对数。

实现细节：

  • 与 logf 类似，通过位操作提取浮点数的指数部分和尾数部分。
  • 通过查表法获取自然对数近似值，最后除以常数 ln⁡(10)\ln(10) 转换为常用对数。

浮点数拆分函数 (modff)

float modff(float x, float *p){
    float i = (int)x;
    *p = i;
    float res = x - i;
    return res;
}

数学原理：

  • modff 函数将浮点数拆分为整数部分和小数部分，使得 x=int+fracx = \text{int} + \text{frac}。

  • 算法步骤：

    1. 计算整数部分 int\text{int}。
    2. 计算小数部分 frac=x−int\text{frac} = x - \text{int}。

实现细节：

  • 使用强制类型转换将浮点数转换为整数。
  • 减去整数部分，得到小数部分。

多项式计算函数 (polyf)

float polyf(float x, int n, float c[])
{
	float sum=0;
	for(int i=0;i<=n;i++)
	{
		sum+=c[i]*int_powf(x,i);
	}
	return sum;
}

float int_powf(float x,int y)
{
	float sum=1;
	for(int i=0;i<y;i++)
	{
		sum*=x;
	}
	return sum;
}

数学原理：

  • polyf 函数计算多项式值，即 ∑i=0nc[i]⋅xi\sum_{i=0}^{n} c[i] \cdot x^i。

  • 算法步骤：

    1. 循环累加各项系数乘以相应的幂次。

实现细节：

  • 使用循环累加计算多项式的各项值。
  • int_powf 函数计算整数次幂，通过循环实现。

幂函数 (powf)

float powf(float x, float y)
{
	float ret, value = x;
	float r = 1.0;
	int p = (int)y;
	if (x == 0.0f && y > 0.0f)
		return 0.0;
	if (y == (float)p)
	{
	    if (p == 0)
	    	return 1.0f;
	    if (p < 0)
	    {
	    	p = -p;
	    	x = 1.0f / x;
	    }
	    while (1)
	    {
	    	if (p & 1)
	    		r *= x;
	    	p >>= 1;
	    	if (p == 0)
	    		return r;
	    	x *= x;
	    }
	}
	ret=expf(y*logf(value));
	return (float)ret;
}

数学原理：

  • powf 函数计算幂次，即 xyx^y。

  • 算法步骤：

    1. 处理整数次幂情况，使用二分法快速幂计算。
    2. 对于非整数次幂，使用对数和指数函数计算。

实现细节：

  • 使用二分法快速幂计算整数次幂。
  • 对于非整数次幂，使用 xy=ey⋅ln⁡(x)x^y = e^{y \cdot \ln(x)} 进行计算。

平方根函数 (sqrtf)

float isqrtf(float x)
{
	float xhalf = 0.5f * x;
	int i = *(int*)&x;
	i = 0x5f375a86 - (i>>1); //魔数，针对32位浮点float的魔数  https://zhuanlan.zhihu.com/p/445813662
	x = *(float*)&i;
	x = x*(1.5f-xhalf*x*x); //第一个迭代
	x = x*(1.5f-xhalf*x*x); //第二个迭代，可移除
	x = x*(1.5f-xhalf*x*x); //第三个迭代，可移除
	return x;
}

float sqrtf(float x)
{
	return x * isqrtf(x);
}

数学原理：

  • sqrtf 函数计算给定值的平方根，即 x\sqrt{x}。
  • 使用快速平方根倒数算法（Fast Inverse Square Root），基于牛顿迭代法优化计算。

实现细节：

  • isqrtf 函数使用魔数 0x5f375a86 和位操作计算平方根倒数的初始近似值，然后进行多次迭代以提高精度。
  • sqrtf 函数通过乘以原始值计算最终的平方根。

正切函数 (tanf)

float tanf(float x)
{
	return (sinf(x) / cosf(x));
}

数学原理：

  • tanf 函数计算给定值的正切值，即 tan⁡(x)\tan(x)。
  • 使用正弦和余弦函数计算

实现细节：

  • 直接调用 sinf 和 cosf 函数计算正弦和余弦值，然后计算其比值。

双曲正切函数 (tanhf)

float tanhf(float x)
{
    return (sinhf(x) / coshf(x));
}

数学原理：

  • tanhf 函数计算给定值的双曲正切值，即 tanh⁡(x)\tanh(x)。
  • 使用双曲正弦和双曲余弦函数计算： tanh⁡(x)=sinh⁡(x)cosh⁡(x)\tanh(x) = \frac{\sinh(x)}{\cosh(x)}

实现细节：

  • 直接调用 sinhf 和 coshf 函数计算双曲正弦和双曲余弦值，然后计算其比值。

总结

这些数学库函数实现了基本的数学运算和变换，使用了查表法、多项式展开、级数展开、快速平方根倒数算法等方法来提高计算的效率和精度。通过条件判断、位操作和循环迭代等手段，确保了计算结果的准确性和可靠性。

avatar
文章
阅读
粉丝