I have seen this piece of code years ago somewhere in the internet, but can't remember exactly where. So, I will be glad if anyone can comment/mail me the original source. However, for spoj, I need to write a much much better version of this in C++ :( Exam sucks life...
import java.math.BigInteger;
import java.security.SecureRandom;
import java.io.*;
import java.util.*;
public class PollardRho {
private final static BigInteger ZERO = new BigInteger("0");
private final static BigInteger ONE = new BigInteger("1");
private final static BigInteger TWO = new BigInteger("2");
private final static SecureRandom random = new SecureRandom();
static Vector<BigInteger> v = new Vector<BigInteger>();
public static BigInteger rho(BigInteger N) {
BigInteger divisor;
BigInteger c = new BigInteger(N.bitLength(), random);
BigInteger x = new BigInteger(N.bitLength(), random);
BigInteger xx = x;
if (N.mod(TWO).compareTo(ZERO) == 0) return TWO;
do {
x = x.multiply(x).mod(N).add(c).mod(N);
xx = xx.multiply(xx).mod(N).add(c).mod(N);
xx = xx.multiply(xx).mod(N).add(c).mod(N);
divisor = x.subtract(xx).gcd(N);
} while((divisor.compareTo(ONE)) == 0);
return divisor;
}
public static void factor(BigInteger N) {
if (N.compareTo(ONE) == 0) return;
if (N.isProbablePrime(20)) {
v.add(N);
return;
}
BigInteger divisor = rho(N);
factor(divisor);
factor(N.divide(divisor));
}
public static void main(String[] args) throws Exception {
String string = "";
InputStreamReader input = new InputStreamReader(System.in);
BufferedReader reader = new BufferedReader(input);
while(null != (string = reader.readLine())) {
BigInteger N = new BigInteger(string);
v.clear();
factor(N);
Collections.sort(v);
for(int i = 0; i < v.size(); i++) System.out.println(v.get(i));
System.out.println();
}
}
}
Thursday, January 13, 2011
Pollard's Rho in Java
Thursday, February 11, 2010
Primitive Root
Primitive root modulo n
This time, lets talk about some sort of modular arithmetic. Suppose, you are asked a peculiar question like the following:
"For a prime number p, check whether g, g2, g3, ... , gp-1 are all distinct (mod p)." [ No doubt, here, p is sufficiently large to prevent you from following naive approach ]
In other words, you are asked:
"Given any prime number p, and another positive number g, you are to tell whether g is a primitive root modulo p."
Now the question is, what is a Primitive Root?
A positive integer g is said to be a primitive root modulo n, if for every integer h such that h and n are co-prime, i.e. GCD(h,n) = 1, there is an integer k such that gk ≡ h (mod n). This also means, g is a generator of the multiplicative group of integers modulo n.
According to the above problem, here, n is actually any prime number p, which makes the condition GCD(h,p) = 1 very obvious.
Example:
The number 3 is a primitive root in (mod 7) because
31=3 ≡ 3 (mod 7)
32=9 ≡ 2 (mod 7)
33=27 ≡ 6 (mod 7)
34=81 ≡ 4 (mod 7)
35=243 ≡ 5 (mod 7)
36=729 ≡ 1 (mod 7)
The remainders above which are 3, 2, 6, 4, 5, and 1 are simply a permutation of 1, 2, ..., (p-1). In this case, p-1=6 because the prime number is 7.
Now, how to solve the problem easily?
As far as I know, there is no simple formula known yet to find primitive roots, however, there is a faster way to test whether a number is a primitive root rather than simply trying out all candidates. This technique can be derived from the definition of primitive roots. From its alternate definition, we know, if the multiplicative order of a number g modulo n is equal to φ(n), then g is a primitive root. In fact, the converse is also true, g is a primitive root modulo n if the multiplicative order of g is φ(n). From this axiom, we get a way to test whether, for any given g and n, g is a primitive root modulo n, as follows:
- First we compute φ(n), [ In case n is a prime p, φ(n) = p-1 ]
- Factorize φ(n), let it has k different prime factors p1, p2, ... , pk, [ Even if m is as large as 231, k will be around 32 ]
- Compute gφ(n)/pi mod n; { for i = 1, 2, 3, ... ,k }. Using a successive squaring algorithm, this can be done very easily and fast. If all the k results are different than 1, then g is a primitive root modulo n.
hu la la... , really easy. From programming approach, φ(n) can be calculated easily by standard division algorithm, dividing by the primes less than √n. [ obviously we can use Sieve of Eratosthenes to get a prime table up to √max, and if n is a prime, φ(n) = n-1 ]. Then factorize φ(n) with the same standard division algorithm and get k different primes pi, { for i = 1, 2, 3, ... , k }. Then the rest can be done in linear time to k where each k modular exponentiation is of logarithmic complexity.
Have fun ...
Thursday, December 10, 2009
Drawing Regular n-gon
How to draw a regular pentagon with just ruler and compass?
How about a regular hexagon?
Primes of the form 22n + 1 are known as Fermat primes.
A regular n-gon is constructible using straightedge (ruler) and compass if and only if n = 2i · m where m is a product of any number of distinct Fermat Prime and i is any natural number, including zero.
Only five Fermat primes are known: 3, 5, 17, 257, and 65,537.
Related study: Compass and straightedge constructions.
Monday, July 13, 2009
Extended Euclidean Algorithm
This method computes expressions of the form ri = axi + byi for the remainder in each step i of the Euclidean algorithm. Each modulus can be written in terms of the previous two remainders and their whole quotient as follows:
By substitution, this gives:
The first two values are the initial arguments to the algorithm:
- r1 = a = a(1) + b(0)
- r2 = b = a(0) + b(1)
The expression for the last non-zero remainder gives the desired results since this method computes every remainder in terms of a and b, as desired.
Example: Compute the GCD of 120 and 23. Or, more formally, compute: x, y, g for 120x + 23y = g; where x, y are two integers and g is the gcd of 120 and 23...
The computation proceeds as follows:- Initial values:
Step 1: Reminder = 120;
Combine terms: 120 = 120 x 1 + 23 x 0
Step 2: Reminder = 23;
Combine terms: 23 = 120 x 0 + 23 x 1
Iterative steps:
Step 3: Quotient = 120 / 23 = 5; Reminder = 120 % 23 = 5;
5 = 120 - 23 x 5
=> 5 = (120 x 1 + 23 x 0) - (120 x 0 + 23 x 1) x 5 ;[from Step 1 and 2]
=> 5 = 120 x 1 + 23 x -5
Step 4: Quotient = 23 / 5 = 4; Reminder = 23 % 5 = 3;
3 = 23 - 5 x 4
=> 3 = (120 x 0 + 23 x 1) - (120 x 1 + 23 x -5) x 4 ;[from Step 2 and 3]
=> 3 = 120 x -4 + 23 x 21
Step 5: Quotient = 5 / 3 = 1; Reminder = 5 % 3 = 2;
2 = 5 - 3 x 1
=> 2 = (120 x 1 + 23 x -5) - (120 x -4 + 23 x 21) x 1 ;[from Step 3 and 4]
=> 2 = 120 x 5 + 23 x -26
Step 6: Quotient = 3 / 2 = 1; Reminder = 3 % 2 = 1;
1 = 3 - 2 x 1
=> 1 = (120 x -4 + 23 x 21) - (120 x 5 + 23 x -26) x 1 ;[from Step 4 and 5]
=> 1 = 120 x -9 + 23 x 47
End of Algorithm.
The last line reads 1 = −9×120 + 47×23, which is the required solution: x = −9 and y = 47, and obviously g = gcd(120,23) = 1
This also means that −9 is the multiplicative inverse of 120 modulo 23, and that 47 is the multiplicative inverse of 23 modulo 120.
- −9 × 120 ≡ 1 mod 23 and also 47 × 23 ≡ 1 mod 120.
By routine algebra of expanding and grouping like terms (refer to the previous example), the following algorithm for iterative method is obtained:
- Apply Euclidean algorithm, and let qn(n starts from 1) be a finite list of quotients in the division.
- Initialize x0, x1 as 1, 0, and y0, y1 as 0,1 respectively.
- Then for each i so long as qi is defined,
- Compute xi+1= xi-1- qixi
- Compute yi+1= yi-1- qiyi
- Repeat the above after incrementing i by 1.
- The answers are the second-to-last of xn and yn.
Here is a sample program written in C++ which implements the Extended Euclidean Algorithm:
#include <stdio.h>
/*
Takes a, b as input.
Returns gcd(a, b).
Updates x, y via pointer reference.
*/
int Extended_Euclid(int A, int B, int *X, int *Y)
{
int x, y, u, v, m, n, a, b, q, r;
/* B = A(0) + B(1) */
x = 0; y = 1;
/* A = A(1) + B(0) */
u = 1; v = 0;
for (a = A, b = B; 0 != a; b = a, a = r, x = u, y = v, u = m, v = n)
{
/* b = aq + r and 0 <= r < a */
q = b / a;
r = b % a;
/* r = Ax + By - aq = Ax + By - (Au + Bv)q = A(x - uq) + B(y - vq) */
m = x - (u * q);
n = y - (v * q);
}
/* Ax + By = gcd(A, B) */
*X = x; *Y = y;
return b;
}
int main()
{
int a, b, x, y, g;
scanf("%d %d", &a, &b);
g = Extended_Euclid(a, b, &x, &y);
printf("X = %d; Y = %d; G = %d\n", x, y, g);
return 0;
}
Python implementation:
def Extended_Euclid(a, b):
x, last_x = 0, 1
y, last_y = 1, 0
while b:
quotient = a / b
a, b = b, a % b
x, last_x = last_x - quotient*x, x
y, last_y = last_y - quotient*y, y
return (last_x, last_y, a)
For complete reading, click here. It is also useful to have a look at the Chinese Remainder Theorem.
Hope this will help...
Sunday, July 12, 2009
Largest impossible score
If two possible points in a game is p and q, what is the largest impossible score of the game which can't be achieved? Or it can't be determined? In other words, for integers p, q ≥ 0, what is the maximum number m which can't be written as a linear combination of p, q, i.e. px + qy, where, x, y ≥ 0.
Solution:
The score can only be determined when p and q are co-primes, i.e. g = GCD(a,b) = 1. Otherwise, we can only make the multiples of g and the remaining numbers can't be made which means the maximum number can't be determined, or infinite! This also covers the case if exactly one of p or q is 0.
One other required condition is p,q > 1. If anyone of p or q is 1, there is no such score which can't be achieved, i.e. answer will be 0. If p and q both are 0, we can't even make any score...
Now, we can arrange the numbers from 0 to p*q-1 in a rectangular array as follows:
0 1 ....... q-1
q q+1 ....... 2q-1
... ... ....... ...
... ... ....... ...
(p-1)q (p-1)(q-1) ....... pq-1
For example, if p = 3 and q = 7, we have the following table:
(0) 1 2 (3) 4 5 (6)
7 8 (9) 10 11 (12) 13
14 (15) 16 17 (18) 19 20
We have circled the multiples of p in the above table and claim these observations:
- There will be q circled numbers. As we have exactly p numbers in each of q columns, so each column will have exactly one circled number, i.e. multiple of p, which implies, there will be a total of q circled number in this table.
- No two circled numbers are in same same column and there will be only one circled number in any column. For the same reason stated above, exactly one circled number in each column.
- Any circled number is possible. This is trivial, because, the circled numbers are multiple of one of the two given numbers, i.e. in the form p*x + q*0 or p*0 + q*x.
- Any number below the circled one in the same column is possible. We have seen any circled number is possible. A circled number is the first multiple of p in the respective column, i.e. n*p, where n ≥1, so the following numbers in the column will be: n*p + q, n*p + 2*q. n*p + 3*q, ... ... ... which are always possible.
- Any number above the circled one in the same column is impossible. This is also explained in the previous step. As any circled number in a column is the first one which is a multiple of one element, the number above this one will be neither the multiple of any of the two numbers. So it is impossible.
- Any number ≥ p*q is possible. As we complete p rows of the table with q columns, we already have one circled number in each column, so, any number > p*q-1 (i.e. the last element of the table) will be surely following any of the q circled numbers, so any number starting from the p+1th row, starting from p*q will be possible.
Now considering all these observations, we can see that, largest impossible number is the number just above the last circled number. The last circled number is p*q-p, i.e. 21-3 = 18 and number above it is (p*q-p)-q, i.e. 18-7 = 11. So in general, largest impossible number is: p*q - p - q.
Conclusion:
If you are given p, q :: two positive numbers, the largest number, m, which can't be expressed as a linear combination of p, q i.e. in the form px + qy, where x and y are two non-negative numbers, is p * q - p - q, when p, q > 1 and GCD(p, q) is 1 i.e. p, q are co-primes.
Wednesday, July 8, 2009
Factorial of negative number
Is it possible to find factorial of negative numbers?
By definition,
N! = 1 x 2 x .... x N
⇒ N! = 1 x 2 x .... x N x (N+1) / (N+1)
⇒ N! = (N+1)! / (N+!)
Now, let's try some values:
3! = 4! / 4 = 6
2! = 3! / 3 = 2
1! = 2! / 2 = 1
0! = 1! / 1 = 1
So, that's how 0! = 1. Let's not stop at 0!, if we continue this fashion, we get:
(-1)! = 0! / 0 = Undetermined
Many people still keeps arguing about whether x / 0 = +∞ or Undefined, it really doesn't matter though, because it is well established that anything divided by zero is undefined.
However, in this post, we are going to talk about a problem from UVa Online Judge, UVa-10323 (Factorial! You Must be Kidding!!!). This problem actually requires you to calculate even the factorials of negative numbers, which does not really exist in real world.
In this context, we are forced to assume that x / 0 = +∞. There was a mistake on my previous post, this section is updated as per the comment from a fellow reader আহনাফ তাহমিদ:
(-1)! = 0! / 0 = +∞
(-2)! = (-1)! / -1 = -∞
(-3)! = (-2)! / -2 = +∞
(-4)! = (-3)! / -3 = -∞
(-5)! = (-4)! / -4 = +∞
and so on..
Hence, we can see that, if N is odd, then (-N)! = +∞, otherwise (-N)! = -∞. Note that, this specific problem statement calls +∞ as Overflow and -∞ as Underflow.
Hope it helps!
Monday, July 6, 2009
Negative Base Number System
The common names for negative-base positional numeral systems are formed by prefixing nega- to the name of the corresponding positive-base system; for example, negadecimal (base -10) corresponds to decimal (base 10), negaternary (base -3) to ternary (base 3), and negabinary (base -2) to binary (base 2).
Denoting the base as − r, every integer a can be written uniquely as:
where each digit dk is an integer from 0 to r − 1 and the leading digit dn is > 0 (unless n=0). The base − r expansion of a is then given by the string dndn-1.....d1d0.
Calculation:
The base − r expansion of a number can be found by repeated division by − r, recording the non-negative remainders of 0,1,2,.........,r-1; and concatenating those remainders, starting with the last. Note that if a / b = c, remainder d, then bc + d = a. For example, 146 in negaternary:
146 / -3 = -48; reminder = 2
-48 / -3 = 16; reminder = 0
16 / -3 = -5; reminder = 1
-5 / -3 = 2; reminder = 1
2 / -3 = 0; reminder = 2
Therefore, the negaternary expansion of 146 is 21102.
Note that in most programming languages, the result (in integer arithmetic) of dividing a negative number by a negative number is rounded towards 0, usually leaving a negative remainder; to get the correct result in such case, computer implementations of the above algorithm should add 1 and r to the quotient and remainder respectively.The numbers -15 to 15 with their expansions in a number of positive and corresponding negative bases are:
| Decimal | Negadecimal | Binary | Negabinary | Ternary | Negaternary |
|---|---|---|---|---|---|
| -15 | 25 | -1111 | 110001 | -120 | 1220 |
| -14 | 26 | -1110 | 110110 | -112 | 1221 |
| -13 | 27 | -1101 | 110111 | -111 | 1222 |
| -12 | 28 | -1100 | 110100 | -110 | 1210 |
| -11 | 29 | -1011 | 110101 | -102 | 1211 |
| -10 | 10 | -1010 | 1010 | -101 | 1212 |
| -9 | 11 | -1001 | 1011 | -100 | 1200 |
| -8 | 12 | -1000 | 1000 | -22 | 1201 |
| -7 | 13 | -111 | 1001 | -21 | 1202 |
| -6 | 14 | -110 | 1110 | -20 | 20 |
| -5 | 15 | -101 | 1111 | -12 | 21 |
| -4 | 16 | -100 | 1100 | -11 | 22 |
| -3 | 17 | -11 | 1101 | -10 | 10 |
| -2 | 18 | -10 | 10 | -2 | 11 |
| -1 | 19 | -1 | 11 | -1 | 12 |
| 0 | 0 | 0 | 0 | 0 | 0 |
| 1 | 1 | 1 | 1 | 1 | 1 |
| 2 | 2 | 10 | 110 | 2 | 2 |
| 3 | 3 | 11 | 111 | 10 | 120 |
| 4 | 4 | 100 | 100 | 11 | 121 |
| 5 | 5 | 101 | 101 | 12 | 122 |
| 6 | 6 | 110 | 11010 | 20 | 110 |
| 7 | 7 | 111 | 11011 | 21 | 111 |
| 8 | 8 | 1000 | 11000 | 22 | 112 |
| 9 | 9 | 1001 | 11001 | 100 | 100 |
| 10 | 190 | 1010 | 11110 | 101 | 101 |
| 11 | 191 | 1011 | 11111 | 102 | 102 |
| 12 | 192 | 1100 | 11100 | 110 | 220 |
| 13 | 193 | 1101 | 11101 | 111 | 221 |
| 14 | 194 | 1110 | 10010 | 112 | 222 |
| 15 | 195 | 1111 | 10011 | 120 | 210 |
Note that the base − r expansions of negative integers have an even number of digits, while the base − r expansions of the non-negative integers have an odd number of digits.
Program for calculating NegabinaryIn python:
def negabinary(i):
digits = []
while i != 0:
i, remainder = divmod(i, -2)
if remainder < 0:
i, remainder = i + 1, remainder + 2
digits.insert(0, str(remainder))
return ''.join(digits)
In C++:
string negabinary(int i)
{
int reminder;
string digits;
while(i != 0)
{
reminder = i % -2;
i /= -2;
if(reminder < 0)
{
i++;
reminder += 2;
}
digits.push_back(reminder+'0');
}
reverse(digits.begin(),digits.end());
return digits;
}
The same procedure it to be followed for calculating other negative number system.
Source: Wikipedia
0.99999999................... = 1
Let, x = 0.99999999999999....................... => 10x = 9.99999999999999....................... => 10x = 9 + 0.99999999999999................... => 10x = 9 + x => 9x = 9 => x = 1 So, 0.99999999999999............................ = 1Is anything going wrong here?
