32位单精度浮点数Payne-Hanek范围缩减算法实现求助

我正在Java项目中开发一个32位单精度浮点数的小型数学函数库。在计算超大参数的正弦值时，无法成功实现Payne-Hanek范围缩减算法。

我研读了《ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit》等多篇文章，自认为理解双精度版本的算法，但多次尝试仍无法完成32位版本的适配。中等参数的正弦值已用Cody-Waite缩减法实现，但参数超过2^23（约109）后精度严重下降。

我的尝试与遇到的错误 #

我基于上述论文修改适配32位浮点数：

• 单精度浮点数尾数23位，指数范围-128至127，我计算得出需要存储至少181位的2/π值，最终用6个32位int存储了192位。
• 尝试拆分x和2/π为10位块进行精确乘法，再通过一系列步骤计算x%(π/2)的余数r。

但存在两处核心错误：

1. 不知如何计算2OverPi=0.6366...的位，应存储为整数形式还是浮点数的位？
2. 无法正确计算浮点数二进制点前的尾随零个数，我的countTrailingZeros()方法逻辑有误。

附上我的实现代码片段：

public static final float piOver2 = 1.57079637051f; //0x3fc90fdb;

/* 192 bits de Pi/2 for reduction : */
public static final int TwoOverPiBits[] = {
        0x00000000,
        0x28be60db,
        0x9391054a,
        0x7f09d5f4,
        0x7d4d3770,
        0x36d8a566,
        0x4f10e410
};

// return the exponant of a float in simple precision
public int getExp(float x) {
    int bits = Float.floatToIntBits(x);
    int exp = ((bits >> 23) & 0xFF) - 127;
    return exp;
}

// return the mantissa of a float in simple precision
public int getMant(float x) {
    int bits = Float.floatToIntBits(x);
    int mant = bits & 0x7fffff;
    return mant;
}

// return the sign of a float in simple precision
public int getSgn(float x) {
    int bits = Float.floatToIntBits(x);
    int sgn = (bits >> 31) & 1;
    return sgn;
}

// Method to extract 55 bits from TwoOverPiBits for the index k and stock them in 10 bits blocks
public int[] packBits(int k) {
    int numPackets = 6;
    int[] packets = new int[numPackets];

    int bitIndex = k; // Index of first bit
    for (int i = 0; i < numPackets; i++) {
        packets[i] = extract12Bits(bitIndex);
        bitIndex += 10;
    }

    return packets;
}

// Method to extract 10 bits from TwoOverPiBits at a given index
public int extract12Bits(int bitIndex) {
    int intIndex = bitIndex / 32; // Index de l'entier dans le tableau
    int bitOffset = bitIndex % 32; // Décalage du bit dans l'entier

    int result;
    if (bitOffset <= 22) {
        // All bits are in the same TwoOverPiBits block
        result = (deuxSurPiBits[intIndex] >> bitOffset) & 0x3FF; // Masqk to extract 10 bits
    } else {
        // 10 bits are in two cosecutive TwoOverPiBits blocks
        int bitsFirstPart = deuxSurPiBits[intIndex] >>> bitOffset;
        int bitsSecondPart = deuxSurPiBits[intIndex + 1] & ((1 << (10 - (32 - bitOffset))) - 1);
        result = (bitsSecondPart << (32 - bitOffset)) | bitsFirstPart;
    }

    return result;
}

//Methode to separate a float in 4 10 bits blocks
public int[] extract12BitPackets(float x) {

    int bits = Float.floatToIntBits(x);

    int[] packets = new int[4];

    // Extract the 4 blocks : 32 bits in 2-10-10-10
    // First block : 10 bits but only 2 that are informative : 1 sign bit and the first exposant bit
    // Second block : 7 exposant bits and 3 mantissa bits
    packets[1] = (bits >> 20) & 0x3FF;

    // third block : 10 mantissa bits
    packets[2] = (bits >> 10) & 0x3FF;

    // last block : 10 last mantissa bits
    packets[3] = (bits & 0x3FF);

    return packets;
}

//Method to count number of trailing zeros before the binary point
public int countTrailingZeros(int number) {
    if (number == 0) return 32;

    int count = 0;
    while ((number & 1) == 0) {
        count++;
        number >>= 1;
    }
    return count;
}


public float payneHanek(float x) {

    if (x<0){
        return -payneHanek(-x);
    }
    int q;
    float f, r;
    int k = getExp(x);
    int mant = getMant(x);
    int M = (k - 23) + countTrailingZeros(mant);

    int[] bits552Overpi = packBits(M - 6); //We begin at M-6 but we then erase the first five bits with the &0x3FF
    int 2Overpibits5 = bits552Overpi[5];
    int 2Overpibits4 = bits552Overpi[4];
    int 2Overpibits3 = bits552Overpi[3];
    int 2Overpibits2 = bits552Overpi[2];
    int 2Overpibits1 = bits552Overpi[1];
    int 2Overpibits0 = bits552Overpi[0] & 0x3FF;

    //Separate x in 10 bits block
    int[] bits32x = extract12BitPackets(x);
    int xbits3 = bits32x[3];
    int xbits2 = bits32x[2];
    int xbits1 = bits32x[1];
    int xbits0 = bits32x[0];

    //x*2OverPi with the blocks
    int y8 = xbits3 * 2Overpibits5;
    int y7 = xbits3 * 2Overpibits4 + xbits2 * 2Overpibits5;
    int y6 = xbits3 * 2Overpibits3 + xbits2 * 2Overpibits4 + xbits1 * 2Overpibits5;
    int y5 = xbits3 * 2Overpibits2 + xbits2 * 2Overpibits3 + xbits1 * 2Overpibits4 + xbits0 * 2Overpibits5;
    int y4 = xbits3 * 2Overpibits1 + xbits2 * 2Overpibits2 + xbits1 * 2Overpibits3 + xbits0 * 2Overpibits4;
    int y3 = xbits3 * 2Overpibits0 + xbits2 * 2Overpibits1 + xbits1 * 2Overpibits2 + xbits0 * 2Overpibits3;
    int y2 = xbits2 * 2Overpibits0 + xbits1 * 2Overpibits1 + xbits0 * 2Overpibits2;
    int y1 = xbits1 * 2Overpibits0 + xbits0 * 2Overpibits1;
    int y0 = xbits0 * 2Overpibits0;


    //Modulo of each blocks : those variables are exact
    y8 = y8 % 4;
    y7 = y7 % 4;
    y6 = y6 % 4;
    y5 = y5 % 4;
    y4 = y4 % 4;
    y3 = y3 % 4;
    y2 = y2 % 4;
    y1 = y1 % 4;
    y0 = y0 % 4;

    //final result :
    float y = y0 + (y1 + (y2 + (y3 + (y4 + (y5 + (y6 + (y7 + y8)))))));
    //error :
    int error = (((((((y1 - y) + y2) + y3) + y4) + y5) + y6) + y7) + y8;
    q = (int) Math.rint(y);
    f = (y - q) + error;
    r = f * piSur2;
    
    //Calculate sin(r)....

1. 2/π的位存储方式 #

Payne-Hanek算法需要的是2/π的二进制小数位序列，按32位一组打包为无符号整数存储：

• 2/π的二进制为0.1001001000011111101101010100010001000010110100011000...
• 存储时按从高位到低位的顺序，将前192位每32位转为一个十六进制整数，存入数组。比如第一个int存小数点后第1到32位，第二个存第33到64位，以此类推。
• 你当前的TwoOverPiBits数组值可能不准确，建议用高精度计算工具（如Python的decimal模块）生成准确的二进制位后再转换。

2. 二进制点前尾随零的计算修正 #

你的countTrailingZeros(mant)逻辑完全错误，因为mant只是浮点数的尾数部分，不是完整的有效数字。论文中M的正确计算方式是：

public int computeM(float x) {
    int exp = getExp(x);
    // 归一化尾数 m = 1 + mant/2^23
    float normalizedMant = 1.0f + getMant(x) / (1 << 23);
    // 计算 m*(2/π) 的值，用于找前导零
    float product = normalizedMant * (2.0f / (float)Math.PI);
    int leadingZeros = 0;
    // 统计二进制小数点后的前导零个数
    while (product < 0.5f) {
        product *= 2;
        leadingZeros++;
    }
    // log2(2/π) ≈ -0.65147，整数部分为-1
    return exp - 1 + leadingZeros;
}

其他关键修正点 #

• 变量名错误：extract12Bits中的deuxSurPiBits应改为TwoOverPiBits。
• x的拆分逻辑错误：不能直接拆分浮点数的原始32位，要将归一化后的x（(1+mant/2^23)*2^exp）拆分为带权重的10位块，而非整数块。
• 求和逻辑错误：取模后的y_i需要按2^-10, 2^-20等权重求和，不是直接相加整数。

替代算法参考 #

如果Payne-Hanek实现难度过高，可参考：

• fdlibm库的单精度正弦实现：其中包含成熟的超大参数范围缩减逻辑。
• 高阶多项式近似：用Remez算法生成拟合多项式，结合简单的范围缩减，精度略低于Payne-Hanek但实现更简单。

内容的提问来源于stack exchange，提问作者dananr