@@ -51,7 +51,7 @@ public static long[][] reduce(long[] a, long[] m) {
5151
5252 // Split each modulus into prime-power factors
5353 for (int i = 0 ; i < a .length ; i ++) {
54- List <Long > factors = PrimeFactorization . primeFactorization (m [i ]);
54+ List <Long > factors = primeFactorization (m [i ]);
5555 Collections .sort (factors );
5656
5757 int j = 0 ;
@@ -132,4 +132,107 @@ static long[] egcd(long a, long b) {
132132 ret [2 ] = tmp ;
133133 return ret ;
134134 }
135+
136+ private static List <Long > primeFactorization (long n ) {
137+ if (n <= 0 )
138+ throw new IllegalArgumentException ();
139+ List <Long > factors = new ArrayList <>();
140+ factor (n , factors );
141+ return factors ;
142+ }
143+
144+ private static void factor (long n , List <Long > factors ) {
145+ if (n == 1 )
146+ return ;
147+ if (isPrime (n )) {
148+ factors .add (n );
149+ return ;
150+ }
151+ long d = pollardRho (n );
152+ factor (d , factors );
153+ factor (n / d , factors );
154+ }
155+
156+ private static long pollardRho (long n ) {
157+ if (n % 2 == 0 )
158+ return 2 ;
159+ long x = 2 + (long ) (999999 * Math .random ());
160+ long c = 2 + (long ) (999999 * Math .random ());
161+ long y = x ;
162+ long d = 1 ;
163+ while (d == 1 ) {
164+ x = mulMod (x , x , n ) + c ;
165+ if (x >= n )
166+ x -= n ;
167+ y = mulMod (y , y , n ) + c ;
168+ if (y >= n )
169+ y -= n ;
170+ y = mulMod (y , y , n ) + c ;
171+ if (y >= n )
172+ y -= n ;
173+ d = gcd (Math .abs (x - y ), n );
174+ if (d == n )
175+ break ;
176+ }
177+ return d ;
178+ }
179+
180+ /**
181+ * Deterministic Miller-Rabin primality test, correct for all long values. Uses 12 witnesses
182+ * guaranteed correct for n < 3.317 × 10^24.
183+ */
184+ private static boolean isPrime (long n ) {
185+ if (n < 2 )
186+ return false ;
187+ if (n < 4 )
188+ return true ;
189+ if (n % 2 == 0 || n % 3 == 0 )
190+ return false ;
191+
192+ long d = n - 1 ;
193+ int r = Long .numberOfTrailingZeros (d );
194+ d >>= r ;
195+
196+ for (long a : new long [] {2 , 3 , 5 , 7 , 11 , 13 , 17 , 19 , 23 , 29 , 31 , 37 }) {
197+ if (a >= n )
198+ continue ;
199+ long x = powMod (a , d , n );
200+ if (x == 1 || x == n - 1 )
201+ continue ;
202+ boolean composite = true ;
203+ for (int i = 0 ; i < r - 1 ; i ++) {
204+ x = mulMod (x , x , n );
205+ if (x == n - 1 ) {
206+ composite = false ;
207+ break ;
208+ }
209+ }
210+ if (composite )
211+ return false ;
212+ }
213+ return true ;
214+ }
215+
216+ private static long powMod (long base , long exp , long mod ) {
217+ long result = 1 ;
218+ base %= mod ;
219+ while (exp > 0 ) {
220+ if ((exp & 1 ) == 1 )
221+ result = mulMod (result , base , mod );
222+ exp >>= 1 ;
223+ base = mulMod (base , base , mod );
224+ }
225+ return result ;
226+ }
227+
228+ private static long mulMod (long a , long b , long mod ) {
229+ return java .math .BigInteger .valueOf (a )
230+ .multiply (java .math .BigInteger .valueOf (b ))
231+ .mod (java .math .BigInteger .valueOf (mod ))
232+ .longValue ();
233+ }
234+
235+ private static long gcd (long a , long b ) {
236+ return b == 0 ? a : gcd (b , a % b );
237+ }
135238}
0 commit comments