天元逆算

#include <gmp.h>

#include <gmp-impl.h>

static const uint8_t INV_TABLE_8[128] = {

1, 171, 205, 183, 57, 163, 197, 239, 241, 27, 61, 167, 41, 19, 53, 223, 225, 139, 173, 151, 25, 131, 165, 207, 209, 251, 29, 135, 9, 243, 21, 191, 193, 107, 141, 119, 249, 99, 133, 175, 177, 219, 253, 103, 233, 211, 245, 159, 161, 75, 109, 87, 217, 67, 101, 143, 145, 187, 221, 71, 201, 179, 213, 127, 129, 43, 77, 55, 185, 35, 69, 111, 113, 155, 189, 39, 169, 147, 181, 95, 97, 11, 45, 23, 153, 3, 37, 79, 81, 123, 157, 7, 137, 115, 149, 63, 65, 235, 13, 247, 121, 227, 5, 47, 49, 91, 125, 231, 105, 83, 117, 31, 33, 203, 237, 215, 89, 195, 229, 15, 17, 59, 93, 199, 73, 51, 85, 255};

typedef union {

mp_limb_t u64; // 64位整数

uint32_t u32[2]; // 两个32位整数

uint16_t u16[4]; // 四个16位整数

uint8_t u8[8]; // 八个8位整数

} LimbUnion;

mp_limb_t inv_n(LimbUnion n){ // a ≡ (-n)⁻¹(mod 2⁶⁴)

LimbUnion a, t

a.u8[0] = INV_TABLE_8[(256 - n.u8[0]) >> 1]; // 查表获得8位初值， // 牛顿迭代精准变形：a <= a + a(an + 1)

t.u32[0]=   a.u8[0] * n.u16[0]; a.u8[1] =(uint8_t) (a.u8[0] *(t.u8[1]  + 1)); // 倍增高半截；致敬牛顿，超越牛顿迭代！

t.u64 =   a.u16[0] * n.u32[0]; a.u16[1]=(uint16_t)(a.u16[0]*(t.u16[1] + 1)); // an的低半截是全1，右移都省了！

t.u64 = (mp_limb_t)(a.u32[0] * n.u64)  ; a.u32[1]=(uint32_t)(a.u32[0]*(t.u32[1] + 1)); // a的低端稳定，所以直接算高半截！

return a.u64

}