载入中...
搜索中...
未找到
inv.c 文件参考
inv.c 的引用(Include)关系图:函数 | |
| mp_limb_t | lmmp_inv_1_ (mp_limb_t x) |
| 1阶逆元计算 (inv1) | |
| mp_limb_t | lmmp_inv_2_1_ (mp_limb_t xh, mp_limb_t xl) |
| 2-1阶逆元计算 (inv21) | |
函数说明
◆ lmmp_inv_1_()
1阶逆元计算 (inv1)
- 参数
-
x 输入的64位无符号整数，最高位为1(MSB(x)=1)
- 返回
- 计算结果：(B^2-1)/x - B
- 警告
- MSB(x)=1, 即x>=2^63
107 {
108 mp_limb_t r, m;
109 {
110 mp_limb_t p, ql;
111 unsigned ul, uh, qh;
112
113 ul = x & LLIMB_MASK;
114 uh = x >> (LIMB_BITS / 2);
115 qh = (x ^ LIMB_MAX) / uh;
116
118 p = (mp_limb_t)qh * ul;
119 if (r < p) {
120 qh--;
121 r += x;
122 if (r >= x)
123 if (r < p) {
124 qh--;
125 r += x;
126 }
127 }
128 r -= p;
129 p = (r >> (LIMB_BITS / 2)) * qh + r;
130 ql = (p >> (LIMB_BITS / 2)) + 1;
133 ql--;
134 r += x;
135 }
137 if (r >= x) {
138 m++;
139 r -= x;
140 }
141 }
142 return m;
143}
引用了 LIMB_BITS, LIMB_MAX , 以及 LLIMB_MASK.
被这些函数引用 lmmp_div_1_(), lmmp_div_1_s_(), lmmp_inv_basecase_(), lmmp_mod_1_() , 以及 lmmp_mulmod_ulong_().
这是这个函数的调用关系图:◆ lmmp_inv_2_1_()
2-1阶逆元计算 (inv21)
- 参数
-
xh 输入数的高64位部分 xl 输入数的低64位部分
- 返回
- 计算结果：(B^3-1)/(xh*B+xl) - B
- 警告
- MSB(xh)=1, 即xh>=2^63
10 {
11 mp_limb_t r, m;
12 {
13 mp_limb_t p, ql;
14 unsigned ul, uh, qh;
15
16 /* For notation, let b denote the half-limb base, so that B = b^2.
17 Split u1 = b uh + ul. */
18 ul = xh & LLIMB_MASK;
19 uh = xh >> (LIMB_BITS / 2);
20
21 /* Approximation of the high half of quotient. Differs from the 2/1
22 inverse of the half limb uh, since we have already subtracted
23 u0. */
24 qh = (xh ^ LIMB_MAX) / uh;
25
26 /* Adjust to get a half-limb 3/2 inverse, i.e., we want
27
28 qh' = floor( (b^3 - 1) / u) - b = floor ((b^3 - b u - 1) / u
29 = floor( (b (~u) + b-1) / u),
30
31 and the remainder
32
33 r = b (~u) + b-1 - qh (b uh + ul)
34 = b (~u - qh uh) + b-1 - qh ul
35
36 Subtraction of qh ul may underflow, which implies adjustments.
37 But by normalization, 2 u >= B > qh ul, so we need to adjust by
38 at most 2.
39 */
40
42
43 p = (mp_limb_t)qh * ul;
44 /* Adjustment steps taken from udiv_qrnnd_c */
45 if (r < p) {
46 qh--;
47 r += xh;
48 if (r >= xh) /* i.e. we didn't get carry when adding to r */
49 if (r < p) {
50 qh--;
51 r += xh;
52 }
53 }
54 r -= p;
55
56 /* Low half of the quotient is
57
58 ql = floor ( (b r + b-1) / u1).
59
60 This is a 3/2 division (on half-limbs), for which qh is a
61 suitable inverse. */
62
63 p = (r >> (LIMB_BITS / 2)) * qh + r;
64 /* Unlike full-limb 3/2, we can add 1 without overflow. For this to
65 work, it is essential that ql is a full mp_limb_t. */
66 ql = (p >> (LIMB_BITS / 2)) + 1;
67
68 /* By the 3/2 trick, we don't need the high half limb. */
70
72 ql--;
73 r += xh;
74 }
76 if (r >= xh) {
77 m++;
78 r -= xh;
79 }
80 }
81
82 /* Now m is the 2/1 inverse of u1. If u0 > 0, adjust it to become a
83 3/2 inverse. */
84 if (xl > 0) {
85 mp_limb_t th, tl;
86 r = ~r;
87 r += xl;
88 if (r < xl) {
89 m--;
90 if (r >= xh) {
91 m--;
92 r -= xh;
93 }
94 r -= xh;
95 }
96 _umul64to128_(xl, m, &tl, &th);
97 r += th;
98 if (r < th) {
99 m--;
100 m -= ((r > xh) | ((r == xh) & (tl > xl)));
101 }
102 }
103
104 return m;
105}
static void _umul64to128_(uint64_t a, uint64_t b, uint64_t *low, uint64_t *high)
Definition longlong.h:31
引用了 _umul64to128_(), LIMB_BITS, LIMB_MAX , 以及 LLIMB_MASK.
被这些函数引用 lmmp_bninv_appr_newton_(), lmmp_div_(), lmmp_div_2_(), lmmp_div_2_s_(), lmmp_div_s_(), lmmp_inv_basecase_() , 以及 lmmp_mod_2_().
函数调用图: 这是这个函数的调用关系图: