This is a follow up of the prior blog post on primes that are congruent to 1 mod 4 being uniquely the sum of squares. In other words, such primes can be written as x2 + y2 with x and y positive integer. There is a natural follow up question, are there any more ways to combine squares of integers to produce primes. And yes, there is one more easily computable example. If a prime is congruent to 1 mod 3, then it can be written uniquely as the three times the square of one integer plus the square of another. In other words, such a prime can be written as 3x2 + y2 with x and y positive integers. And like for primes congruent to 1 mod 4, there is an efficient way to compute the break down of the prime into the squares. For more on this general topic see Fermat Sum of Two Squares.
For the primes congruent to 1 mod 4, we used the Gaussian integers and applied the Euclidean algorithm. For the primes congruent to 1 mod 3, we use the Cyclotomic 3rd root of unity integers and still use the Euclidean algorithm. Note, the Gaussian integers are actually the Cyclotomic 4th root of unity integers. See Cyclotomic Fields Wiki for a definition of Cyclotomic fields and the definition of a ring of integers inside the field.
The key fact is that Cyclotomic 3 and Cyclotomic 4 fields are quadratic extensions of the rational numbers. And the general theory of Quadratic Forms says that finding solutions mod a prime will, for enough primes and powers of primes, will give you a solution in the integers. For our two simple cases, finding a solution mod p is sufficient for a mod p case to be lifted to integer solutions.
The code for the case of primes congruent to 1 mod 3 is similar enough to the code for primes congruent to 1 mod 4, that instead of providing code just to do the 1 mod 3 case, the code does both cases simultaneously. In particular if a prime is congruent to 1 mod 12, you will see a break down into a sum of squares and the breakdown into three times a square plus a square.
As in the prior blog post, you can execute the code at Groovy Web Console. In the results, the mod 4 case is prefixed with M4 and the mod 3 case is prefixed with M3. We write out the break down in ascii math. So 3*2^2 + 5^2 should be read as three times the square of the number 2 plus the square of the number 5 giving a result of 37. The code, given a value for start, will look for primes greater or equal to start that are congruent to one mod 4 or one mod 3. When it finds the primes, it will then break them down into squares. If you click on the Result tab at the bottom of the web console, it will show you the primes and break down into squares with the number of primes being determined by the value of the count variable.
// Used to quickly show an integer is not divisible by any prime 13 or less.
def gridN = 4 * 3 * 5 * 7 * 11 * 13
boolean[] mod4Array = new boolean[gridN]
boolean[] mod3Array = new boolean[gridN]
boolean[] fullArray = new boolean[gridN]
def testN = [2, 3, 5, 7, 11, 13]
for (int i = 5; i < gridN; i++) {
boolean foundIt = false
for (def n : testN) {
if ((i % n) == 0) {
foundIt = true
break
}
}
if (!foundIt) {
mod3Array[i] = (i % 3) == 1
mod4Array[i] = (i % 4) == 1
fullArray[i] = true
}
}
// Find all primes up to *endRange*. We test primes greater than 13 directly against
// any potential prime *p*.
int endRange = 50_000
List<Integer> additionalTestN = []
// Compute primes up to *endRange*.
boolean[] primeArray = new boolean[endRange]
for (int i = 17; i < endRange; i++) {
if (!primeArray[i] && fullArray[i]) {
additionalTestN.add(i)
int lastMult = endRange.intdiv(i)
for (int j = 2; j < lastMult; j++) {
primeArray[i * j] = true
}
}
}
// Example start numbers.
BigInteger strongPseudoPrime = 1373653 // If you wish to use this, set eliminateSmallFactors to false
BigInteger largerPseudoPrime = 1195_068768_795265_792518_361315_725116_351898_245581 // Pseudo to bases up to 36.
BigInteger smallPrimes = 17
BigInteger superLargePrimes = 1E200 // Primes with 200 digits.
boolean printStartingExamples = false
def startingExamples = [strongPseudoPrime, largerPseudoPrime, smallPrimes, superLargePrimes]
if (printStartingExamples) {
println(startingExamples)
}
// Change start and count as desired.
BigInteger start = smallPrimes
int count = 10 // Number of primes to find starting at *start*.
// Whether to eliminate candidates who have prime factors less than *endRange*. Set to false when hunting for pseudo primes.
boolean eliminateSmallFactors = true
int i = 0
BigInteger cur = start
BigInteger end = start + count * 2000
boolean doMod4 = true
boolean doMod3 = true
boolean printThoseWithFactors = false
List<S> results = []
while (cur < end && i < count) {
boolean canDoMod3 = mod3Array[cur % gridN] && doMod3
boolean canDoMod4 = mod4Array[cur % gridN] && doMod4
// Make sure the number is not divisible by prime 13 or less.
if (canDoMod3 || canDoMod4) {
// Eliminate potential primes that have divisors less than *endRange*.
Integer foundFactor = null
if (eliminateSmallFactors) {
for (def n : additionalTestN) {
if ((cur % n) == 0) {
if (cur == n) {
break
}
foundFactor = n
}
}
}
def p = cur
if (foundFactor == null) {
List<Boolean> l = []
if (canDoMod3) {
l.add(false)
}
if (canDoMod4) {
l.add(true)
}
// Any factor of *p* is now guaranteed to be larger than *endRange* if *eliminateSmallFactors* was true
println("--[count=${i + 1}] Potential prime ${f(p)}")
S result = new S()
l.each {isGaussian ->
String prefix = isGaussian ? "Mod4 " : "Mod3 "
CS cs = squares(p, isGaussian)
if (cs) {
def gcd = cs.gcd
println("${prefix}Result ${gcd}[squares-norm=${f(gcd.norm())}]")
if (isGaussian) result.mod4 = cs else result.mod3 = cs
} else {
// For very large *p*, *p* is likely to still not be a prime, even though not divisible
// by any prime less than *endRange*.
println("${prefix}Cannot find generator for ${f(p)}")
}
}
if (result.mod3 || result.mod4) {
results.add(result)
i++
}
} else {
if (printThoseWithFactors) {
println("--Potential prime ${f(p)} has factor ${foundFactor}")
}
}
}
cur++
}
class CS {
C gcd
Integer notPrimeGen
String toValString() {
String prefix = gcd.isGaussian ? "M4=" : "M3="
String notPrimeGenStr = notPrimeGen ? "[notPrimeGen=${notPrimeGen}]" : ""
return prefix + notPrimeGenStr + gcd.sqVariant()
}
}
class S {
CS mod4
CS mod3
String toString() {
if (!mod4 && !mod3) {
return ""
}
List<CS> l = []
if (mod4) {
l.add(mod4)
}
if (mod3) {
l.add(mod3)
}
def norms = l.collect {it.gcd.norm()}
def n = norms[0]
if (norms.find {it != n} != null) {
throw new Exception("Mismatched norms ${norms}")
}
"${F.f(n)} -> ${l.collect {it.toValString()}.join(" ")}"
}
}
// Find number field integers whose norm is congruent to 1 mod p. We do this looking for potential generators for the cyclic
// group multiplicative group of integers mod p. We also, use this to potentially disprove that `p` is a prime.
// Note, if isGaussian is false, then doing 3rd root of unity ring integers.
CS squares(BigInteger p, boolean isGaussian) {
int testInt = isGaussian ? 4 : 6
def r = (p - 1).intdiv(testInt)
if (r * testInt != p - 1) {
println("Potential prime not congruent to 1 mod ${testInt}")
return null
}
Integer notPrimeGen = null
List<Integer> generators = []
List<BigInteger> powers = []
// We find all potential generators from a short list. If none of them disprove that *p* is a prime, then *p* is extremely
// likely to be a prime.
def toTry = [2, 3, 5, 7, 11, 13, 17, 37, 325, 9375, 28178, 450775]
for (def gen : toTry) {
if (gen >= p) {
break
}
BigInteger pow = powN(gen, r, p)
def mult = (pow * pow) % p
if (!isGaussian) {
mult = (mult * pow) % p
}
BigInteger tPow = mult
if (tPow == p - 1 && pow != p - 1) {
generators.add(gen)
powers.add(pow)
} else if (tPow != 1 && pow != p - 1) {
notPrimeGen = gen
break
}
}
String noPrimeMsg = notPrimeGen != null ? ", not prime because of generator ${notPrimeGen}" : ""
println("Generators = " + generators + noPrimeMsg + ".")
CS retval = null
if (powers.size() > 0) {
retval = new CS()
retval.gcd = gcd(new C(p, 0, isGaussian), new C(powers[0], 1, isGaussian))
retval.notPrimeGen = notPrimeGen
}
return retval
}
// The x to the power r mod n, but do so in C*log(r) steps.
static BigInteger powN(Number xx, Number r, Number n) {
BigInteger x = xx as BigInteger
if (x < 0) {
BigInteger t = -x % n
x = t == 0 ? 0 : n - t
}
// Unwrap the more typical recursion algorithm. When the numbers get too large, the recursion exceeds the max
// Java stack.
List<Boolean> computeChain = []
def curR = r
while (curR > 2) {
long b = curR % 2
Number s = b == 1 ? curR - 1 : curR
computeChain.add(b == 1)
curR = s.intdiv(2)
if (computeChain.size() > 2000) {
throw new Exception("Too many descent steps on taking exponent mod ${f(n)}")
}
}
BigInteger startSeed
if (curR <= 0) {
startSeed = 1
} else if (curR == 1) {
startSeed = x % n
} else {
startSeed = (x * x) % n
}
computeChain = computeChain.reverse()
def curVal = startSeed
for (def b : computeChain) {
curVal = (curVal * curVal) % n
if (b) {
curVal = (curVal * x) % n
}
}
return curVal
}
// Is either a gaussian or cyclotomic-3 integer (the integers in the numbers containing a 3rd root of unity)
class C {
boolean isGaussian // If false, then 3rd root of unity ring of integers. If true, then integers inside complex numbers.
BigInteger x // Real part
BigInteger y // Coefficient for root of unity (imaginary part of complex number if doing Gaussian integers).
C(Number x, Number y, boolean isGaussian) {
this.isGaussian = isGaussian
this.x = x as BigInteger
this.y = y as BigInteger
}
// Actually the square of the norm.
BigInteger norm() {
def n = x * x + y * y
if (!isGaussian) {
n = n - x * y
}
return n
}
C conj() {
def nx = x
def ny = -y
if (!isGaussian) {
nx = nx - y
}
new C(nx, ny, isGaussian)
}
C mult(C o) {
def nx = x * o.x - y * o.y
def ny = x * o.y + y * o.x
if (!isGaussian) {
ny = ny - y * o.y
}
return new C(nx, ny, isGaussian)
}
// Divide by integer and round to nearest integer for each coordinate
C round(BigInteger divisor) {
def xd = round(x, divisor)
def yd = round(y, divisor)
return new C(xd, yd, isGaussian)
}
static BigInteger round(BigInteger n, BigInteger divisor) {
BigInteger halfDivisor = divisor.intdiv(2)
boolean isNeg = n < 0
def nn = isNeg ? -n : n
def nnd = nn.intdiv(divisor)
def remainder = nn % divisor
if (remainder > halfDivisor) {
// Round up.
nnd += 1
}
return isNeg ? -nnd : nnd
}
// Apply quotient remainder division in number field integers
C remainder(C d) {
def p = mult(d.conj())
def intPart = p.round(d.norm())
def m = d.mult(intPart)
return new C(x - m.x, y - m.y, isGaussian)
}
C bestVersion() {
if (isGaussian) {
return new C(x.abs(), y.abs(), true)
}
// We currently have a version that satisfied x*x - x*y + y*y = p (where p is the prime).
// There is a simple algorithm for producing a solution to 3*x*x + y*y for our current
// x and y, but it does not always produce integers. In order to produce integers
// x and y need to be both odd.
// In what follows, we transform x and y to a different solution of x*x - x*y + y*y = p, but
// where x and y are both odd, x is positive and the absolute value of x is bigger than
// the absolute value of y. There is only one such solution. From that we
// can produce our solution to 3*x*x + y*y = p.
def x1 = x
def y1 = y
if (y1 % 2 == 0) {
x1 = y
y1 = x
}
if (x1 % 2 == 0) {
x1 = x1 - y1
y1 = -y1
}
if (y1.abs() > x1.abs()) {
def t = x1
x1 = y1
y1 = t
}
if (x1 < 0) {
x1 = -x1
y1 = -y1
}
return new C(x1, y1, false)
}
boolean isZero() {
return x == 0 && y == 0
}
String toString() {
return "(" + F.f(x) + "," + F.f(y) + ")"
}
// Produce the display version of the sum of squares combination.
String sqVariant() {
def best = bestVersion()
String n1
String n2
if (isGaussian) {
n1 = F.f(best.x)
n2 = F.f(best.y)
} else {
// Convert to solution to 3*x*x + y*y = p.
n1 = F.f((best.x - best.y).intdiv(2))
n2 = F.f((best.x + best.y).intdiv(2))
}
return toSqs(n1, n2)
}
String toSqs(String n1, String n2) {
String cf = isGaussian ? "" : "3*"
return cf + n1 + "^2" + "+" + n2 + "^2"
}
}
// Do the Euclidean algorithm for determining Greatest Common Divisor in Cyclotomic 3 and 4 root ring integers.
def gcd(C x, C y) {
def curX = x
def curY = y
int maxLoop = 10_000
int count = 0
// Unwrap the more typical recursion algorithm. When the numbers get to large, the recursion exceeds the max
// Java stack.
while (true) {
def rem = curX.remainder(curY)
if (rem.isZero()) {
return curY.bestVersion()
}
println("x = " + curX + "; y = " + curY + "; remainder=" + rem)
if (rem.norm() >= curY.norm()) {
throw new Exception("Remainder failed for ${curX} and ${curY}, and remainder=${rem}")
}
if (count++ > maxLoop) {
throw new Exception("Too many iterations of algorithm for GCD")
}
curX = curY
curY = rem
}
}
// Class used so other classes can access function *f*.
class F {
// Format big integer with underscores every six characters (instead of commas every three characters).
// Generally those who present very large numbers use a separator every six characters instead of every three.
static String f(Number inN) {
boolean isNeg = inN < 0
def n = isNeg ? -inN : inN
String s = n.toString()
def l = s.length()
def beginOffset = l % 6
if (beginOffset == 0) {
beginOffset = 6
}
StringBuilder out = new StringBuilder(s.length() * 2)
if (isNeg) {
out.append("-")
}
out.append(s, 0, beginOffset)
int offset = beginOffset
while (offset < l) {
out.append("_")
out.append(s, offset, offset + 6)
offset += 6
}
return out.toString()
}
}
static String f(Number inN) {
return F.f(inN)
}
results.collect {r -> r.toString()}.join("\n")