just a test

ibm.py

from fractions import gcd

M = 19*41*47*61
totient = 18*40*46*60

def p(n):
  return (n*n*n*n + n*n + n + 1) % M

# compute inverse of  a  mod  M
#  or None if no inverse exists

def inv(a):
	if gcd(a,M) != 1: return None
	return pow(a,totient-1,M)

#  given (a,b,c), test if  p(n) == q(n) for all n
#  where  q(n) = (n^2 + an + b)(n^2 - an + c)

def testpoly(a,b,c):
	q = lambda n: ((n*n + a*n + b) * (n*n - a*n + c)) % M
	return all([p(n) == q(n) for n in range(M)])

#  search for solution to the following system of (mod M) equations
#  bc = 1
#  a(c-b) = 1
#  b+c-aa = 1

onehalf = inv(2)

for a in range(1,M):
	if inv(a) is None: continue
	c = ((inv(a) + a*a + 1) * onehalf) % M
	b = (1 + a*a - c) % M
	if (b*c % M) == 1:
		print a,b,c,
		print "*" if testpoly(a,b,c) else " "

print "~fini~"