Primes One Mod Four

// See *start* and *count* later in code.

// Used to quickly show an integer is not divisible by any prime 13 or less.
def gridN = 4 * 3 * 5 * 7 * 11 * 13
boolean[] array = 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) {
        array[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 = superLargePrimes
int count = 2 // 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 printThoseWithFactors = false

List<S> results = []

while (cur < end && i < count) {
   // Make sure the number is not divisible by prime 13 or less.
   if (array[cur % gridN]) {
       // 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) {
           // 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 = squares(p)
           if (result) {
                def gcd = result.gcd
                println("Result ${result.gcd}[sum-squares=${f(gcd.norm())}] for potential prime ${f(p)}")
                results.add(result)
                i++
           } 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("Cannot find generator for ${f(p)}")
           }
       } else {
           if (printThoseWithFactors) {
                println("--Potential prime ${f(p)} has factor ${foundFactor}")
           }
       }
   }
    cur++
}
class S {
    C gcd
    Integer notPrimeGen
    String toString() {
        String notPrimeGenStr = notPrimeGen ? "[notPrimeGen=${notPrimeGen}]" : ""
       "${F.f(gcd.norm())}${notPrimeGenStr} -> " + gcd.toString()
   }
}

// Find Gaussian 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.
S squares(BigInteger p) {
    def r = (p - 1).intdiv(4)
   if (r * 4 != p - 1) {
        println("Potential prime not congruent to 1 mod 4")
       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)
        BigInteger tPow = (pow * pow) % p
       if (tPow == p - 1) {
           // Note, gen is not necessarily a generator. It is just a non-square mod p, which is sufficient.
           // However, it is likely to be a generator, hence our abuse of language.
           generators.add(gen)
            powers.add(pow)
       } else if (tPow != 1) {
            notPrimeGen = gen
           break
       }
   }
    String noPrimeMsg = notPrimeGen != null ? ", not prime because of generator ${notPrimeGen}" : ""
    println("Generators = " + generators + noPrimeMsg + ".")
    S retval = null
   if (powers.size() > 0) {
        retval = new S()
        retval.gcd = gcd(new C(p, 0), new C(powers[0], 1))
        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) {
    def 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 = []
    BigInteger curR = r as BigInteger
   while (curR > 2) {
       long b = (curR % 2) as Long

        BigInteger 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)}")
       }
   }
    BigDecimal 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
}

// Gaussian integers (the integers in the complex numbers)
class C {
    BigInteger x // Real part
   BigInteger y // Imaginary part
   C(Number x, Number y) {
       this.x = x as BigInteger
       this.y = y as BigInteger
   }

   // Actually the square of the norm.
   BigInteger norm() { return x * x + y * y}

    C conj() {
       return new C(x, -y)
   }

    C mult(C o) {
       return new C(x * o.x - y * o.y, x * o.y + y * o.x)
   }

   // 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)
   }

   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 Gaussian 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)
   }

    C abs() {
       return new C(x.abs(), y.abs())
   }

   boolean isZero() {
       return x == 0 && y == 0
   }

    String toString() {
       return "(" + F.f(x) + "," + F.f(y) + ")"
   }
}

// Do the Euclidean algorithm for determining Greatest Common Divisor between two Gaussian integers.
def gcd(C x, C y) {
    def curX = x
    def curY = y
   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.abs()
       }
        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++ > 10_000) {
           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")