Some NTL Routines I have written
I've written a few NTL routines to
investigate Number Theory, and will list some of them
here for reference. Hopefully you will find them helpful,
if not for your own use, then for understanding how they work.
IsSquare()
Often time you will need to test if a number
is square. The obvious approach is to take the integer square root,
square it, and test if you get back the same number. However, this
can really hurt performance, especially if you are doing this thousands
if not millions of times. Fortunately there is an easy trick that can
help you identify most non-squares - quadratic residues. If you
test various squares modulo 3, you will notice that the result is
never 2. This is because 2 is not a 'quadratic residue' modulo 3.
With further tests you will notice that 2 and 3 are not quadratic
residues modulo 5. Using sets of non-residues, we can quickly
eliminate numbers that aren't square. Here is the routine I
typically use in my NTL programs...
// Returns 1 if x is a perfect square. Otherwise 0.
int IsSquare(ZZ &x) {
int z = x % 4849845; // 4849845 = 3*5*7*11*13*17*19
// Do some quick tests on x to see if we
// can quickly eliminate it as a square
// using quadratic-residues.
if(z%3 == 2) return 0;
int t;
t = z%5; if((t==2)||(t==3)) return 0;
t = z%7; if((t==3)||(t==5)||(t==6)) return 0;
t = z%11; if((t==2)||(t==6)||(t==7)||(t==8)||(t==10)) return 0;
t = z%13; if((t==2)||(t==5)||(t==6)||(t==7)||(t==8)||(t==11)) return 0;
t = z%17; if((t==3)||(t==5)||(t==6)||(t==7)||(t==10)||(t==11)||(t==12)||(t==14)) return 0;
t = z%19; if((t==2)||(t==3)||(t==8)||(t==10)||(t==12)||(t==13)||(t==14)||(t==15)||(t==18)) return 0;
// If we get here, we'll have to just try taking
// square-root & compare its square...
static ZZ r;
SqrRoot(r, abs(x));
if(r*r == abs(x)) return 1;
return 0;
}
PellSolve()
Solves the Brouncker/Pell Equation x^2 - D*y^2 = +1 given D.
// Solves x^2-D*y^2 = +1, given D.
// Uses traditional continued-fraction approach.
void PellSolve(const ZZ &D, ZZ &x, ZZ &y) {
// Return trivial solution for D <= 1...
if(D <= 1) {
x = 1;
y = 0;
return;
}
static ZZ pn[5];
static ZZ qn[5];
static ZZ an[5];
static ZZ Pn[5];
static ZZ Qn[5];
ZZ a1;
ZZ t;
int i;
// Initialize continued fraction sequence...
Pn[0] = 0;
Qn[0] = 1;
SqrRoot(a1, D);
an[0] = a1;
Pn[1] = a1;
Qn[1] = D - a1*a1;
// On perfect-squares, return trivial solution.
if(Qn[1] == 0) {
x = 1;
y = 0;
return;
}
qn[0] = 1;
pn[0] = a1;
an[1] = (a1 + Pn[1])/Qn[1];
qn[1] = an[1];
pn[1] = a1*an[1]+1;
for(i=2;;i++) {
// Traverse chain...
Pn[i%4] = an[(i-1)%4] * Qn[(i-1)%4] - Pn[(i-1)%4];
Qn[i%4] = (D - Pn[i%4]*Pn[i%4])/Qn[(i-1)%4];
an[i%4] = (a1 + Pn[i%4])/Qn[i%4];
pn[i%4] = an[i%4]*pn[(i-1)%4] + pn[(i-2)%4];
qn[i%4] = an[i%4]*qn[(i-1)%4] + qn[(i-2)%4];
// Found repetition?
if(an[i%4] == a1*2) {
// Obtain solution...
x = pn[(i-1)%4];
y = qn[(i-1)%4];
t = x*x - D*y*y;
// If this is a -1 solution, transform into a +1...
if(t == -1) {
t = x;
x = x*x + D*y*y;
y = 2*t*y;
}
// exit...
break;
}
}
}
Tries to find a factor of a number using the Pollard p-1 method.
Before calling you should verify the number isn't prime using ProbPrime(), and check for
small divisors (e.g. < 100,000) using trial-division.
ZZ pollard(ZZ &n) {
PrimeSeq seq;
ZZ p, t;
ZZ b;
b = 7;
// Init with lots of small primes...
for(int i=0; i<20; i++) {
p = seq.next();
for(int j=0; j<4; j++) {
PowerMod(b, b, p, n);
}
}
// Continue with one of every other prime...
for(i=0;;i++) {
p = seq.next();
PowerMod(b, b, p, n);
if(i%50 == 0) {
GCD(t, b-1, n);
if( (t!=1) && (t!=n) ) {
return t;
}
i = 0;
// Throw in some extra products for good measure and luck...
PowerMod(b, b, b, n);
}
}
}
Tries to find a factor of a number using the rho method.
Before calling you should verify the number isn't prime using ProbPrime(), and check for
small divisors (e.g. < 100,000) using trial-division.
ZZ rho(ZZ &n) {
ZZ x=5, y=7, Q=1, t;
int a = 3;
for(int i=0;;i++) {
x = (x*x-a)%n;
y = (y*y-a)%n;
y = (y*y-a)%n;
Q = (Q*(y-x))%n;
if(i%100 == 0) {
GCD(t, Q, n);
if( (t!=1) && (t!=n) ) return t;
}
}
}
KthRoot()
Computes the Kth integer root using simple binary search.
inline ZZ KthRoot(ZZ &n, int k) {
static ZZ high, low, mid, tmp;
high = n;
low = 1;
mid = (high + low)/2;
while(high > low+1) {
power(tmp, mid, k);
if(tmp > n) high = mid;
else if(tmp < n) low = mid;
else return mid;
mid = (high + low)/2;
}
power(tmp, high, k);
if(tmp < n) return high;
return low;
}
IsCube()
Works just like IsSquare(), but for cubes. Falls back
on the KthRoot() routine if the number isn't eliminated by cubic non-residues.
// Returns 1 if x is a perfect cube. Otherwise 0.
int IsCube(ZZ &x) {
int z = x % 53599; // 53599 = 7*13*19*31
// Do some quick tests on x to see if we
// can quickly eliminate it as a cube
// using cubic-residues.
int t;
t = z%7; if( (t>=2)&&(t<=5) ) return 0;
t = z%13; if( (t<=11)&&(t>=2)&&(t!=5)&&(t!=8) ) return 0;
t = z%19; if( (t>=2)&&(t<=6) ) return 0;
if( (t==9) || (t==10) ) return 0;
if( (t>=13) && (t<=17) ) return 0;
t = z%31; if(! ((t==0)||(t==1)||(t==2)||(t==4)||(t==8)||(t==15)||(t==16)||(t==23)||(t==27)||(t==29)||(t==30)||(t==31)) ) return 0;
// If we get here, we'll have to just try taking
// cube root & compare its cube...
static ZZ r;
r = KthRoot(abs(x), 3);
if(r*r*r == abs(x)) return 1;
return 0;
}
LinearCRT()
A lot of times in number theory you can find two different
linear forms for a number. For instance, you might determine in one way that
a number you are searching for is of the form 3x+1, and in another way determine
it must also be of the form 5y+2. Based on these two restrictions you can conclude
the number must be of the form 15z+7. This routine allows you to make this kind
of derivation, taking two such linear forms and returning one combined linear form.
// Take two linear forms ax+b & cy+d and create one
// form that assumes both values mz+r. We set a=m, b=r.
// This routine is advantageous to typical CRT in that
// we allow gcd(a,c)>1.
// Returns 0 if contradiction or if cy+d doesn't help
// restrict ax+b. Otherwise returns 1.
int LinearCRT(ZZ &a, ZZ &b, ZZ c, ZZ d) {
// Swap greater modulus in case one doesn't help
// we'll stick with the most restrictive
int swapped = 0;
if(c > a) {
ZZ t;
t = a;
a = c;
c = t;
t = b;
b = d;
d = t;
swapped = 1;
}
// Reduce residues...
if(a>0) b %= a;
if(c>0) d %= c;
// Throw out contradictions and useless forms...
if( (a<=1) || (c<=1) ) return swapped;
ZZ g;
GCD(g, a, c);
while(g > 1) {
if(b%g != d%g) return 0;
c /= g;
d %= c;
if(c<=1) return swapped;
GCD(g, a, c);
}
// m is easy to calculate now
ZZ m;
m = a*c;
// Find r such that r = b (mod a) and r = d (mod c)
ZZ y, r;
y = 1;
CRT(r, y, b, a);
CRT(r, y, d, c);
r %= m;
// Set results and return
a = m;
b = r;
return 1;
}
SquareSum()
Finds representation of given prime as a sum of squares.
Returns 0 if there exists no such representation. It's a good algorithm, able to handle primes with a thousand digits or more!
int SquareSum(ZZ p, ZZ &A, ZZ &B) {
ZZ u, v, t, r, x;
// Trivial case...
if(p==2) {
A = 1;
B = 1;
return 1;
}
// All p>2 must equal 1 mod 4
if(p%4 != 1) return 0;
// Otherwise we can find solution to x^2 = -1 (mod p)
// Btw: Explicit such solution is +/- [(p-1)/2]!
SqrRootMod(x, p-1, p);
// Here we do Euclidean algorithm until quotient is less than sqrt(p)
SqrRoot(r, p);
u = p;
v = x;
do {
t = u % v;
u = v;
v = t;
if(t < r) {
A = t;
SqrRoot(B, p - t*t);
if(A*A + B*B == p) {
A = u;
B = v;
return 1;
}
else return 0;
}
} while (v != 0);
return 0;
}
ChooseFirst(), ChooseNext()
During investigations, you may need to iteratively
obtain a subset of values from a larger group. This can be confusing
to implement, so I have some handy routines ready. Here they are, along
with an example routine demonstrating how they can be used.
// Variables used...
int g_chooseN = 0;
ZZ *g_chooseArr = 0;
ZZ *g_chooseTgt = 0;
int g_chooseK = 0;
int *g_loopCounter = 0;
// Choose k values from arr[n] and put them in tgt[k]
int ChooseFirst(ZZ *arr, int n, ZZ *tgt, int k) {
g_chooseArr = arr;
g_chooseTgt = tgt;
if(g_loopCounter) {
free(g_loopCounter);
g_loopCounter = 0;
}
g_loopCounter = (int *)malloc(sizeof(int) * k);
if(!g_loopCounter) return 0;
for(int i=0; i< k; i++) {
g_loopCounter[i] = i;
}
g_chooseK = k;
g_chooseN = n;
// set indices
for(i=0; i< k; i++) {
g_chooseTgt[i] = g_chooseArr[ g_loopCounter[i] ];
}
return 1;
}
// Used internally (do not call), returns TRUE if successful, otherwise ZERO.
int IncrementLoop(int pos) {
if(pos < 0) return 0;
if((g_loopCounter[pos] < g_chooseN-1) && (g_loopCounter[pos] + g_chooseK-pos < g_chooseN)) {
g_loopCounter[pos]++;
// Now re-base subsequent loop counters...
for(int i=1; pos+i< g_chooseK; i++) {
if(g_loopCounter[pos+i] < g_loopCounter[pos]+i) {
g_loopCounter[pos+i] = g_loopCounter[pos]+i;
// Went to far?
if(g_loopCounter[pos+i] >= g_chooseN) return 0;
// todo: Above test could be done without iteration
}
}
}
else {
g_loopCounter[pos] = pos;
if(!IncrementLoop(pos-1)) return 0;
}
return 1;
}
// Return TRUE if new value stored, otherwise ZERO.
int ChooseNext(void) {
if(!IncrementLoop(g_chooseK-1)) return 0;
// set indices
for(int i=0; i< g_chooseK; i++) {
g_chooseTgt[i] = g_chooseArr[ g_loopCounter[i] ];
}
return 1;
}
// Example of how it is used...
void ChooseTest(void) {
ZZ val[8];
ZZ tgt[8];
for(int i=0; i<8; i++) val[i] = i+1;
int k = 4;
ChooseFirst(val, 8, tgt, k);
do {
cout << "{" << tgt[0] << flush;
outs << "{" << tgt[0] << flush;
for(int i=1; i< k; i++) {
cout << "," << tgt[i] << flush;
outs << "," << tgt[i] << flush;
}
cout << "}\n" << flush;
outs << "}\n" << flush;
} while(ChooseNext());
}
Joe K. Crump |
Number Theory Home