@@ -54,7 +54,6 @@ public static uint ModOddInverse(uint[] m, uint[] x, uint[] z)
54
54
55
55
int bits = ( len32 << 5 ) - Integers . NumberOfLeadingZeros ( ( int ) m [ len32 - 1 ] ) ;
56
56
int len30 = ( bits + 29 ) / 30 ;
57
- int m0Inv30x4 = - ( int ) Inverse32 ( m [ 0 ] ) << 2 ;
58
57
59
58
int [ ] t = new int [ 4 ] ;
60
59
int [ ] D = new int [ len30 ] ;
@@ -69,32 +68,33 @@ public static uint ModOddInverse(uint[] m, uint[] x, uint[] z)
69
68
Array . Copy ( M , 0 , F , 0 , len30 ) ;
70
69
71
70
int eta = - 1 ;
71
+ int m0Inv32 = ( int ) Inverse32 ( ( uint ) M [ 0 ] ) ;
72
72
int maxDivsteps = GetMaximumDivsteps ( bits ) ;
73
73
74
74
for ( int divSteps = 0 ; divSteps < maxDivsteps ; divSteps += 30 )
75
75
{
76
76
eta = Divsteps30 ( eta , F [ 0 ] , G [ 0 ] , t ) ;
77
- UpdateDE30 ( len30 , D , E , t , m0Inv30x4 , M ) ;
77
+ UpdateDE30 ( len30 , D , E , t , m0Inv32 , M ) ;
78
78
UpdateFG30 ( len30 , F , G , t ) ;
79
79
}
80
80
81
81
int signF = F [ len30 - 1 ] >> 31 ;
82
- Debug . Assert ( - 1 == signF | 0 == signF ) ;
83
82
84
- CNegate30 ( len30 , signF , F ) ;
85
-
86
- int signD = CNegate30 ( len30 , signF , D ) ;
87
- Debug . Assert ( - 1 == signD | 0 == signD ) ;
88
-
89
- // TODO 'D' should already be in [P, -P), but absent a proof we support [-2P, 2P)
90
- signD = CSub30 ( len30 , ~ signD , D , M ) ;
91
- signD = CAdd30 ( len30 , signD , D , M ) ;
92
- signD = CAdd30 ( len30 , signD , D , M ) ;
93
- Debug . Assert ( 0 == signD ) ;
83
+ /*
84
+ * D is in the range (-2.M, M). First, conditionally add M if D is negative, to bring it
85
+ * into the range (-M, M). Then normalize by conditionally negating (according to signF)
86
+ * and/or then adding M, to bring it into the range [0, M).
87
+ */
88
+ int signD = D [ len30 - 1 ] >> 31 ;
89
+ signD = CAdd30 ( len30 , signD , D , M ) ;
90
+ CNormalize30 ( len30 , signF , D , M ) ;
94
91
95
92
Decode30 ( bits , D , 0 , z , 0 ) ;
96
93
Debug . Assert ( 0 != Nat . LessThan ( len32 , z , m ) ) ;
97
94
95
+ signF = CNegate30 ( len30 , signF , F ) ;
96
+ Debug . Assert ( 0 == signF ) ;
97
+
98
98
return ( uint ) ( EqualTo ( len30 , F , 1 ) & EqualToZero ( len30 , G ) ) ;
99
99
}
100
100
@@ -107,7 +107,6 @@ public static bool ModOddInverseVar(uint[] m, uint[] x, uint[] z)
107
107
108
108
int bits = ( len32 << 5 ) - Integers . NumberOfLeadingZeros ( ( int ) m [ len32 - 1 ] ) ;
109
109
int len30 = ( bits + 29 ) / 30 ;
110
- int m0Inv30x4 = - ( int ) Inverse32 ( m [ 0 ] ) << 2 ;
111
110
112
111
int [ ] t = new int [ 4 ] ;
113
112
int [ ] D = new int [ len30 ] ;
@@ -124,6 +123,7 @@ public static bool ModOddInverseVar(uint[] m, uint[] x, uint[] z)
124
123
int clzG = Integers . NumberOfLeadingZeros ( G [ len30 - 1 ] | 1 ) - ( len30 * 30 + 2 - bits ) ;
125
124
int eta = - 1 - clzG ;
126
125
int lenDE = len30 , lenFG = len30 ;
126
+ int m0Inv32 = ( int ) Inverse32 ( ( uint ) M [ 0 ] ) ;
127
127
int maxDivsteps = GetMaximumDivsteps ( bits ) ;
128
128
129
129
int divsteps = 0 ;
@@ -135,7 +135,7 @@ public static bool ModOddInverseVar(uint[] m, uint[] x, uint[] z)
135
135
divsteps += 30 ;
136
136
137
137
eta = Divsteps30Var ( eta , F [ 0 ] , G [ 0 ] , t ) ;
138
- UpdateDE30 ( lenDE , D , E , t , m0Inv30x4 , M ) ;
138
+ UpdateDE30 ( lenDE , D , E , t , m0Inv32 , M ) ;
139
139
UpdateFG30 ( lenFG , F , G , t ) ;
140
140
141
141
int fn = F [ lenFG - 1 ] ;
@@ -154,32 +154,30 @@ public static bool ModOddInverseVar(uint[] m, uint[] x, uint[] z)
154
154
}
155
155
156
156
int signF = F [ lenFG - 1 ] >> 31 ;
157
- Debug . Assert ( - 1 == signF | 0 == signF ) ;
158
-
159
- if ( 0 != signF )
160
- {
161
- Negate30 ( lenFG , F ) ;
162
- Negate30 ( lenDE , D ) ;
163
- }
164
-
165
- if ( ! IsOne ( lenFG , F ) )
166
- return false ;
167
157
158
+ /*
159
+ * D is in the range (-2.M, M). First, conditionally add M if D is negative, to bring it
160
+ * into the range (-M, M). Then normalize by conditionally negating (according to signF)
161
+ * and/or then adding M, to bring it into the range [0, M).
162
+ */
168
163
int signD = D [ lenDE - 1 ] >> 31 ;
169
- Debug . Assert ( - 1 == signD | 0 == signD ) ;
170
-
171
- // TODO 'D' should already be in [P, -P), but absent a proof we support [-2P, 2P)
172
164
if ( signD < 0 )
173
165
{
174
- signD = Add30 ( len30 , D , M ) ;
166
+ signD = Add30 ( lenDE , D , M ) ;
175
167
}
176
- else
168
+ if ( signF < 0 )
177
169
{
178
- signD = Sub30 ( len30 , D , M ) ;
170
+ signD = Negate30 ( lenDE , D ) ;
171
+ signF = Negate30 ( lenFG , F ) ;
179
172
}
173
+ Debug . Assert ( 0 == signF ) ;
174
+
175
+ if ( ! IsOne ( lenFG , F ) )
176
+ return false ;
177
+
180
178
if ( signD < 0 )
181
179
{
182
- signD = Add30 ( len30 , D , M ) ;
180
+ signD = Add30 ( lenDE , D , M ) ;
183
181
}
184
182
Debug . Assert ( 0 == signD ) ;
185
183
@@ -263,21 +261,22 @@ private static int CNegate30(int len30, int cond, int[] D)
263
261
return c ;
264
262
}
265
263
266
- private static int CSub30 ( int len30 , int cond , int [ ] D , int [ ] M )
264
+ private static void CNormalize30 ( int len30 , int condNegate , int [ ] D , int [ ] M )
267
265
{
268
266
Debug . Assert ( len30 > 0 ) ;
269
267
Debug . Assert ( D . Length >= len30 ) ;
270
268
Debug . Assert ( M . Length >= len30 ) ;
271
269
272
270
int c = 0 , last = len30 - 1 ;
271
+ int condAdd = ( D [ last ] >> 31 ) ^ condNegate ;
273
272
for ( int i = 0 ; i < last ; ++ i )
274
273
{
275
- c += D [ i ] - ( M [ i ] & cond ) ;
274
+ c += ( D [ i ] ^ condNegate ) - condNegate + ( M [ i ] & condAdd ) ;
276
275
D [ i ] = c & M30 ; c >>= 30 ;
277
276
}
278
- c += D [ last ] - ( M [ last ] & cond ) ;
279
- D [ last ] = c ; c >>= 30 ;
280
- return c ;
277
+ c += ( D [ last ] ^ condNegate ) - condNegate + ( M [ last ] & condAdd ) ;
278
+ D [ last ] = c ;
279
+ Debug . Assert ( c >> 30 == 0 ) ;
281
280
}
282
281
283
282
private static void Decode30 ( int bits , int [ ] x , int xOff , uint [ ] z , int zOff )
@@ -349,7 +348,7 @@ private static int Divsteps30Var(int eta, int f0, int g0, int[] t)
349
348
int f = f0 , g = g0 , m , w , x , y , z ;
350
349
int i = 30 , limit , zeros ;
351
350
352
- for ( ; ; )
351
+ for ( ; ; )
353
352
{
354
353
// Use a sentinel bit to count zeros only up to i.
355
354
zeros = Integers . NumberOfTrailingZeros ( g | ( - 1 << i ) ) ;
@@ -502,46 +501,46 @@ private static int Negate30(int len30, int[] D)
502
501
return c ;
503
502
}
504
503
505
- private static int Sub30 ( int len30 , int [ ] D , int [ ] M )
506
- {
507
- Debug . Assert ( len30 > 0 ) ;
508
- Debug . Assert ( D . Length >= len30 ) ;
509
- Debug . Assert ( M . Length >= len30 ) ;
510
-
511
- int c = 0 , last = len30 - 1 ;
512
- for ( int i = 0 ; i < last ; ++ i )
513
- {
514
- c += D [ i ] - M [ i ] ;
515
- D [ i ] = c & M30 ; c >>= 30 ;
516
- }
517
- c += D [ last ] - M [ last ] ;
518
- D [ last ] = c ; c >>= 30 ;
519
- return c ;
520
- }
521
-
522
- private static void UpdateDE30 ( int len30 , int [ ] D , int [ ] E , int [ ] t , int m0Inv30x4 , int [ ] M )
504
+ private static void UpdateDE30 ( int len30 , int [ ] D , int [ ] E , int [ ] t , int m0Inv32 , int [ ] M )
523
505
{
524
506
Debug . Assert ( len30 > 0 ) ;
525
507
Debug . Assert ( D . Length >= len30 ) ;
526
508
Debug . Assert ( E . Length >= len30 ) ;
527
509
Debug . Assert ( M . Length >= len30 ) ;
528
- Debug . Assert ( m0Inv30x4 * M [ 0 ] == - 1 << 2 ) ;
510
+ Debug . Assert ( m0Inv32 * M [ 0 ] == 1 ) ;
529
511
530
512
int u = t [ 0 ] , v = t [ 1 ] , q = t [ 2 ] , r = t [ 3 ] ;
531
- int di , ei , i , md , me ;
513
+ int di , ei , i , md , me , mi , sd , se ;
532
514
long cd , ce ;
533
515
516
+ /*
517
+ * We accept D (E) in the range (-2.M, M) and conceptually add the modulus to the input
518
+ * value if it is initially negative. Instead of adding it explicitly, we add u and/or v (q
519
+ * and/or r) to md (me).
520
+ */
521
+ sd = D [ len30 - 1 ] >> 31 ;
522
+ se = E [ len30 - 1 ] >> 31 ;
523
+
524
+ md = ( u & sd ) + ( v & se ) ;
525
+ me = ( q & sd ) + ( r & se ) ;
526
+
527
+ mi = M [ 0 ] ;
534
528
di = D [ 0 ] ;
535
529
ei = E [ 0 ] ;
536
530
537
531
cd = ( long ) u * di + ( long ) v * ei ;
538
532
ce = ( long ) q * di + ( long ) r * ei ;
539
533
540
- md = ( m0Inv30x4 * ( int ) cd ) >> 2 ;
541
- me = ( m0Inv30x4 * ( int ) ce ) >> 2 ;
534
+ /*
535
+ * Subtract from md/me an extra term in the range [0, 2^30) such that the low 30 bits of the
536
+ * intermediate D/E values will be 0, allowing clean division by 2^30. The final D/E are
537
+ * thus in the range (-2.M, M), consistent with the input constraint.
538
+ */
539
+ md -= ( m0Inv32 * ( int ) cd + md ) & M30 ;
540
+ me -= ( m0Inv32 * ( int ) ce + me ) & M30 ;
542
541
543
- cd += ( long ) M [ 0 ] * md ;
544
- ce += ( long ) M [ 0 ] * me ;
542
+ cd += ( long ) mi * md ;
543
+ ce += ( long ) mi * me ;
545
544
546
545
Debug . Assert ( ( ( int ) cd & M30 ) == 0 ) ;
547
546
Debug . Assert ( ( ( int ) ce & M30 ) == 0 ) ;
@@ -551,14 +550,12 @@ private static void UpdateDE30(int len30, int[] D, int[] E, int[] t, int m0Inv30
551
550
552
551
for ( i = 1 ; i < len30 ; ++ i )
553
552
{
553
+ mi = M [ i ] ;
554
554
di = D [ i ] ;
555
555
ei = E [ i ] ;
556
556
557
- cd += ( long ) u * di + ( long ) v * ei ;
558
- ce += ( long ) q * di + ( long ) r * ei ;
559
-
560
- cd += ( long ) M [ i ] * md ;
561
- ce += ( long ) M [ i ] * me ;
557
+ cd += ( long ) u * di + ( long ) v * ei + ( long ) mi * md ;
558
+ ce += ( long ) q * di + ( long ) r * ei + ( long ) mi * me ;
562
559
563
560
D [ i - 1 ] = ( int ) cd & M30 ; cd >>= 30 ;
564
561
E [ i - 1 ] = ( int ) ce & M30 ; ce >>= 30 ;
@@ -606,4 +603,4 @@ private static void UpdateFG30(int len30, int[] F, int[] G, int[] t)
606
603
G [ len30 - 1 ] = ( int ) cg ;
607
604
}
608
605
}
609
- }
606
+ }
0 commit comments