Pythagorean Factorization Plagiarized

"""
Quick ~1000 line proof to show Pythagorean Factorization is just masked Fermat's
Factorization Algorithm.

In the absence of reference implementation, here is the best-effort replication
of the Pythagorean Factorization, based on the author Robert Edward Grant's
illustration (https://archive.is/nN5Ai):
"""
import math

def pythagorean_factorization(semiprime: int):
    B = math.sqrt(semiprime)
    r = math.ceil(B) - B
    x = 0
    while True:
        # Solve leg A with the formula in Grant's image
        A = math.sqrt((x + r) * (2 * B + (x + r)))

        # Note: In Grant's illustration, the 'ready-to-use'
        # formula for solving A is missing the square root
        # used above, but it is present in the bottom right
        # of the image, where A=24 is obtained with sqrt(576).

        C = round(math.sqrt(A ** 2 + B ** 2))

        PRIME1 = round(C - A)
        PRIME2 = round(C + A)

        if PRIME1 * PRIME2 == semiprime:
            return PRIME1, PRIME2

        x += 1


"""
Below is a 2013 implementation of Fermat's Factorization Algorithm from the
1700th century. The implementation precedes Grant's supposed algorithm by seven
years. It can be found at https://stackoverflow.com/a/20465181
"""
def isqrt(n):
    x = n
    y = (x + n // x) // 2
    while y < x:
        x = y
        y = (x + n // x) // 2
    return x

def fermat(n, verbose = True):
    a = isqrt(n)  # int(ceil(n**0.5))
    b2 = a * a - n
    b = isqrt(n)  # int(b2**0.5)
    count = 0
    while b * b != b2:
        if verbose:
            print('Trying: a=%s b2=%s b=%s' % (a, b2, b))
        a = a + 1
        b2 = a * a - n
        b = isqrt(b2)  # int(b2**0.5)
        count += 1
    p = a + b
    q = a - b
    return p, q


"""
To show Grant's Pythagorean Factorization is the exact same algorithm as
Fermat's factorization algorithm, we can compare the number of steps and see
that they indeed take an equal number for each input, but we can also do this
mathematically and just move things around and show they're equivalent.

We're going to do that in four parts.
    1. Clean unnecessary stuff from Fermat implementation's code
    2. Switch to using the same variable names
    3. Make the Fermat's algorithm more readable
    4. Show that the algorithms are the same, and that Grant is just obfuscating
       his plagiarism
"""


# Part I - Cleaning the code


"""
Let's start by

* removing the `isqrt` function from above, as it's just an implementation of
  math.isqrt, which we can import from the math library,

* removing the code comments (trailing stuff that starts with '#'),

* removing the
    if verbose:
        print('Trying: a=%s b2=%s b=%s' % (a, b2, b))
  part as unnecessary,

* removing the lines
      count = 0
  and
      count +=1
  for now, as we're not analyzing the algorithm's efficiency, just yet.

After every change description like this, the code below will reflect the
changes.
"""
from math import isqrt

def version1(n):
    a = isqrt(n)
    b2 = a * a - n
    b = isqrt(n)

    while b * b != b2:
        a = a + 1
        b2 = a * a - n
        b = isqrt(b2)

    p = a + b
    q = a - b
    return p, q


"""
The line
    from math import isqrt
doesn't need to be done again and again so we'll leave it here at the top of
this document.
"""


# Part II - Switch to using the same variable names


"""
Next, let's swap some names of variables so that they match those of the
Pythagorean Factorization.

Since finding the prime factors `p` and `q` have identical final steps in both
algorithms, we can identify some variable names to change.

      Fermat  |     Grant
    ---------------------------
    p = a+b   |  PRIME2 = C + A
    q = a-b   |  PRIME1 = C - A

So, let's change:
    p → PRIME2
    q → PRIME1
    a → C
    b → A

The parameter `n` is the semiprime that Grant incorrectly calls the
"public key".

An RSA public key consists of two parts
    1) the modulus (n), that is the semiprime (product of two large primes
       numbers), and
    2) a public exponent, usually 65537.

But since breaking the public key is about finding the factors of the modulus,
we'll just ignore this as another instance of "Grant has never read a book on
the topic".

In Grant's illustration, the side B of the right triangle is the square root of
the modulus `n`. So,
    B² = n
and thus,
    B = √(n)

To make the variables in the algorithm easier to recognize, we'll change them
to lowercase.

We'll also agree on a notation that capitalized letters represent squares of
the legs and the hypotenuse, i.e.
    A = a²
    B = b²
    C = c²


    ^                                           ^
    | \                                         | \
    |  \                                        |  \
    |   \                                       |   \
    |    \                                      |    \
    |     \                                     |     \
    |      \                                    |      \
    |       \  C                                |       \  c
 B  |        \                   🡆          b  |        \
    |         \                                 |         \
    |          \                                |          \
    |           \                               |           \
    |            \                              |            \
    |             \                             |             \
    |              \                            |              \
    |_______________\                           |_______________\
            A                                           a

    Grant's triangle                            New variable names


Pythagorean theorem respectively

      A² + B² = C²                🡆               a² + b² = c²

                                            (or)   A  + B  = C    (as upper case
                                                                   represents
                                                                   square value)

With the new syntax, now
    b² = n
that is,
    B  = n

We thus know `n` represents a square of the longer leg in Grant's illustration,
and we can replace `n` with `B` in the algorithm.

That together with demoting the original `A` and `C` to `a` and `c`, gives us
"""
def version2(B):
    c = isqrt(B)
    b2 = c * c - B
    a = isqrt(B)

    while a * a != b2:
        c = c + 1
        b2 = c * c - B
        a = isqrt(b2)

    PRIME2 = c + a
    PRIME1 = c - a
    return PRIME2, PRIME1


"""
Next, let's fix the wonky variable name `b2`.

As per Pythagoras,
    a² + b² = c²
we can subtract b² from both sides and get
    a² = c² - b²

The middle statement inside the while loop says
    b2 = c * c - B
or
    b2 =   c²  - b²
we thus now know that
    b2 = a² = A,

and we can replace `b2` with `A`.

Let's also swap the order of PRIME2 and PRIME1 as it's annoying PRIME2 is listed
first.
"""
def version3(B):
    c = isqrt(B)
    A = c * c - B
    a = isqrt(B)

    while a * a != A:
        c = c + 1
        A = c * c - B
        a = isqrt(A)

    PRIME1 = c - a
    PRIME2 = c + a
    return PRIME1, PRIME2


"""
All variable names are now replaced.
"""


# PART III - Make the Fermat's algorithm more readable

"""
Next, let's swap the places of the two first lines in the function
    A = c * c - B
    a = isqrt(B)
to better show `a` and `c` are set to the same variable.
"""
def version4(B):
    c = isqrt(B)
    a = isqrt(B)
    A = c * c - B

    while a * a != A:
        c = c + 1
        A = c * c - B
        a = isqrt(A)

    PRIME1 = c - a
    PRIME2 = c + a
    return PRIME1, PRIME2


"""
Next, let's ensure we can use capitalized letters without adding lines. This is
done by setting a new variable on the same line. We can separate a new statement
from the previous one in Python with the `;` symbol.

This allows us to define `C` immediately when we know `c`, by doing
         ; C = c*c
on the same line, we define `c`, e.g.,
    c = isqrt(B)                        ; C = c*c

We'll use that below, and we can then use C in place of c*c. We'll just have to
reassign `C` every time `c` updates.
"""
def version5(B):
    c = isqrt(B)                        ; C = c * c
    a = isqrt(B)
    A = C - B

    while a * a != A:
        c = c + 1                       ; C = c * c
        A = C - B
        a = isqrt(A)

    PRIME1 = c - a
    PRIME2 = c + a
    return PRIME1, PRIME2


"""
Next, we can combine the lines
    A = C - B
    a = isqrt(A)
to
    a = isqrt(C - B)

This allows us to see the connection to the Pythagorean theorem, where
    a = √(C - B)
or
    a² = c² - b²      (for reference, the Pythagorean theorem is a² + b² = c²)
"""
def version6(B):
    c = isqrt(B)                        ; C = c * c
    a = isqrt(B)
    A = C - B

    while a * a != A:
        c = c + 1                       ; C = c * c
        a = isqrt(C - B)

    PRIME1 = c - a
    PRIME2 = c + a
    return PRIME1, PRIME2


"""
Finally, we can replace the while loop condition's `A` with `C-B` to help us
understand the solution detection condition. This also allows us to remove the
third line of the function
    A = C - B
as `A` is no longer needed.
"""
def version7(B):
    c = isqrt(B)                        ; C = c * c
    a = isqrt(B)

    while a * a != C - B:
        c = c + 1                       ; C = c * c
        a = isqrt(C - B)

    PRIME1 = c - a
    PRIME2 = c + a
    return PRIME1, PRIME2


"""
At this point, we have made Fermat's algorithm as readable as possible.

The algorithm performs the following steps:
    1. Assume leg `a` and hypotenuse `c` are both the integer square root of B,
       the modulus.
           An integer square root will first calculate a normal square root to
       get an accurate decimal approximation, and it then rounds it down to the
       next integer. To compare with normal square root:
           sqrt(5) = 2.2360679...
           isqrt(5) = 2
    2. Raise the integer square root `a` to the second power and see if it
       completes the Pythagorean theorem
           a*a = C  - B
           a²  = c² - b²

    3. Enter a `while` -loop with `x != y` condition. It means the loop
       continues as long as `x` and `y` are not equal. Inside the loop
        a) Increase the hypotenuse `c` by 1 to get a new hypotenuse candidate.
        b) Uses isqrt to get new leg `a` candidate from known leg `b`
              a = isqrt(c² - b²)
        c) Jump back to check if the new `a`, when squared, completes the
           Pythagorean theorem for the triangle. The check
              a² == c² - b² ?
           will remain False, i.e. still not equal, for quite some time,
           so it jumps back to step a) to increase `c` by one. Trying all values
           (of `c`) is called brute forcing, and it's the worst efficiency an
           attack can have. It's essentially the "try all keys" method.

        Once the a*a will equal c² - b², the loop exits. At this point, it knows
          * `c`, the hypotenuse, is the middle point between the two primes.
          * `a`, the bottom side leg, is the distance of primes from the center.

    4. Calculate the prime factors with
           PRIME2 = c + a
           PRIME1 = c - a

    6. Return the prime factors to the user.
"""


# Part IV - Plagiarism


"""

Fermat's algorithm describes the factoring method beautifully. But being so
simple, plagiarizing is trivial to detect. Its nature of brute forcing the
hypotenuse one integer at a time immediately shows its time complexity is 𝒪(√n).

But we need to plagiarize the algorithm to attract fame and suckers. So, we
can turn to our buddy, Nassim Haramein for advice. Nassim has a long history
of obfuscating existing math by replacing variables and moving things around:
    https://youtu.be/_W2WBeqGNM0?si=v_QiR5um05k1xaYp&t=1884

So let's obfuscate the algorithm by brute-forcing the leg `a` instead of brute-
forcing the hypotenuse `c`. Since leg `b` is already known, we can't do much
else anyway.

To do that, we need to remove our guess for `a`, i.e. the second line
    a = isqrt(B)

(At this stage the algorithm won't work, so we'll denote these intermediate,
temporarily broken factoring implementations with a "_broken"-suffix.)

We'll also comment out the broken versions with '#' chars so they won't cause
issues if this file is opened in a code editor.
"""
# def version8_broken(B):
#     c = isqrt(B)                      ; C = c*c
#
#     while a * a != C - B:
#         c = c + 1                     ; C = c*c
#         a = isqrt(C - B)
#
#     PRIME1 = c - a
#     PRIME2 = c + a
#     return PRIME1, PRIME2


"""
The first trouble is, we need `a` immediately in the while loop condition. So
let's move the brute-force-is-done check
    a * a != C - B
inside the loop. Since our while loop continued as long as there was an
inequality, an if-statement checking for the condition inside the while loop
must check for equality, to know when to `break` from the loop.

The while loop's condition is changed to `True` so that the while loop's
condition check itself never terminates the loop.
"""
# def version9_broken(B):
#     c = isqrt(B)                      ; C = c*c
#
#     while True:
#         c = c + 1                     ; C = c*c
#         a = isqrt(C - B)
#
#         if a * a == C - B:
#             break
#
#     PRIME2 = c + a
#     PRIME1 = c - a
#     return PRIME2, PRIME1

"""
Since we're no longer using the first `C` in first ;-separated `C = c*c`, we can
delete that.
"""
# def version10_broken(B):
#     c = isqrt(B)
#
#     while True:
#         c = c + 1                     ; C = c*c
#         a = isqrt(C - B)
#
#         if a * a == C - B:
#             break
#
#     PRIME1 = c - a
#     PRIME2 = c + a
#     return PRIME1, PRIME2


"""
Next, to simplify, we can move the prime-forming to below the if-condition, and
return factors from the loop immediately after obtaining them.
"""
# def version11_broken(B):
#     c = isqrt(B)
#
#     while True:
#         c = c + 1                     ; C = c*c
#         a = isqrt(C - B)
#
#         if a * a == C - B:
#             PRIME1 = c - a
#             PRIME2 = c + a
#             return PRIME1, PRIME2


"""
We've now solved the problem of `a` missing from the initial loop break
condition check. But we have another problem.

This algorithm immediately increases `c` by 1. Thus, it will jump over some
solutions if prime factors are very close to one another.

But thanks to the Pythagorean theorem, we can find `c` from `b` with `a`, as
    c = sqrt(a² + b²)

That means we can delete the entire line
    c = c + 1                     ; C = c*c

And replace it with
    c = sqrt(a² + B).
Raising to exponents in Python is done with '**' operator, so we do
    c = sqrt(a**2 + B)

"""
# from math import sqrt
#
# def version12_broken(B):
#     c = isqrt(B)
#
#     while True:
#         a = isqrt(C - B)
#         c = sqrt(a**2 + B)
#
#         if a * a == C - B:
#             PRIME1 = c - a
#             PRIME2 = c + a
#             return PRIME1, PRIME2


"""
Since we can't get `a` candidate with
    a = isqrt(C - B)
any more, we'll delete the line.
"""
# def version13_broken(B):
#     c = isqrt(B)
#
#     while True:
#         c = sqrt(a ** 2 + B)
#
#         if a * a == C - B:
#             PRIME1 = c - a
#             PRIME2 = c + a
#             return PRIME1, PRIME2


"""
And now the top row
    c = isqrt(B)
has become useless, so we remove it too.
"""
# def version14_broken(B):
#
#     while True:
#         c = sqrt(a ** 2 + B)
#
#         if a * a == C - B:
#             PRIME1 = c - a
#             PRIME2 = c + a
#             return PRIME1, PRIME2


"""
Next, let's figure out how to solve `a` without `c`.

The algorithm's exit condition is already
    A  = C  - B
    a² = c² - b²

Now, as per difference of squares ( https://en.wikipedia.org/wiki/Difference_of_two_squares ):
    a²-b² = (a+b)(a-b)
Thus, we can substitute
    c² - b²
on the right with
    (b+c)(c-b)
and get
    a² = (b+c)(c-b)
and since multiplication is commutative, we can change the order to
    a² = (c-b)(b+c)

Let's then add some spaces
    a² = (c-b)(b  +  c   )

And then use a bit of mathematical
intuition, and add `b-b` (which is 0) to the right side.
                                             ↓
    a² = (c-b)(b  + b-b + c   )

Let's then rearrange the terms to get another c-b out:
    a² = (c-b)(b+b       -b +c )
    a² = (c-b)(2b        + c-b )
                           ---

Now, let's replace the c-b with two new variables, Blonk and Blorp:
    https://youtu.be/_W2WBeqGNM0?si=v_QiR5um05k1xaYp&t=1884

Let's define
    c-b = x+r

Now
    a² =    (c-b)(2b + (c-b))
becomes
    a² =    (x+r)(2b + (x+r))
    a  = √( (x+r)(2b + (x+r)) )

And we're done. We have solved for `a` without `c`.

We can add this as
    a = sqrt( (x+r) * (2*b + (x+r)) )
"""
# def version15_broken(B):
#
#     while True:
#         a = sqrt( (x+r) * (2*b + (x+r)) )
#         c = sqrt(a ** 2 + B)
#
#         if a * a == C - B:
#             PRIME1 = c - a
#             PRIME2 = c + a
#             return PRIME1, PRIME2


"""
Now we notice that we are missing three variables, `x`, `r`, and `b`. Of these,
`b` is the easy one, it's just √B. So we'll start with that, and add
    b = sqrt(B)
to the beginning.
"""
# def version16_broken(B):
#     b = sqrt(B)
#
#     while True:
#         a = sqrt( (x+r) * (2*b + (x+r)) )
#         c = sqrt(a ** 2 + B)
#
#         if a * a == C - B:
#             PRIME1 = c - a
#             PRIME2 = c + a
#             return PRIME1, PRIME2


"""
Next, let's change the way we check for completeness. We'll first get the prime
candidates without the check, and then see if the two prime factors, when
rounded as integers produce the semiprime B.
"""
# def version17_broken(B):
#     b = sqrt(B)
#
#     while True:
#         a = sqrt((x + r) * (2 * b + (x + r)))
#         c = sqrt(a ** 2 + B)
#
#         PRIME1 = round(c - a)
#         PRIME2 = round(c + a)
#         if PRIME1 * PRIME2 == B:
#             return PRIME1, PRIME2


"""
Next, let's solve r. Since the
    a² + b² = c²
only works with integers when the solution is correct, `r` will simply serve as
a value that rounds the irrational value b = √B, to the next integer.

For this, we need what's called a ceiling function, which rounds to the next
integer. The ceiling function is denoted with ⌈ and ⌉, and e.g. ⌈5.1⌉ = 6.
In Python, the equivalent is math.ceil(5.1).

So
    √B + r = ⌈√B⌉
is thus
         r = ⌈√B⌉ - √B

Or in Python,
    from math import ceil
#  (r =     ⌈   √ B ⌉ -    √ B)
    r = ceil(sqrt(B)) - sqrt(b)

And since we already calculated
    b = sqrt(B)
on the first line of the function, we can get r with
    r = ceil(b) - b

Which is what Grant did in his illustration.
"""
# def version18_broken(B):
#     b = sqrt(B)
#     r = ceil(b) - b
#
#     while True:
#         a = sqrt((x + r) * (2 * b + (x + r)))
#         c = sqrt(a ** 2 + B)
#
#         PRIME1 = round(c - a)
#         PRIME2 = round(c + a)
#         if PRIME1 * PRIME2 == B:
#             return PRIME1, PRIME2


"""
Now, because we're incompetent or malicious (pick your poison), we will remove
the square root from our formula for solving `a`:
    a = sqrt((x + r) * (2 * b + (x + r)))
becomes
    a =      (x + r) * (2 * b + (x + r))

This is the state in which Grant was peddling the algorithm to his viewers.
"""
# from math import sqrt, ceil
#
# def version19_broken(B):
#     b = sqrt(B)
#     r = ceil(b) - b
#
#     while True:
#         a = (x + r) * (2 * b + (x + r))
#         c = sqrt(a ** 2 + B)
#
#         PRIME1 = round(c - a)
#         PRIME2 = round(c + a)
#         if PRIME1 * PRIME2 == B:
#             return PRIME1, PRIME2


"""
Intentional or not, removing the square root means the algorithm isn't working.

Also, he didn't explain where x comes from, he just baked its value (2) into the
image:
    https://archive.is/nN5Ai
The reason x is so small is he intentionally used an invalid RSA key where the
prime factors are practically next to one another.

The reader who manages to follow the steps in Grant's illustration with no
beginning or end will fail to notice that the description "x is a whole number"
is doing a LOT of lifting, that is, trillions and trillions of super-computing
years' worth of 'lifting'. Finding the right `x` is the key to factoring numbers
with this algorithm. It's the number of brute force steps. And unlike the
illustration implies, it's never, ever 2. For average-size RSA keys (3072 bits),
it's verified by key generation algorithms to be at least

    2^((3072/2)-100)) = \
        19014012429663796324419459331224055645026991688903
        74749247369562864802585451004209832696515053742316
        84659417885497720016514658819276907584086181741652
        52126912455130430925740762843286365008491375209869
        58696226722405264137600414778995896732994290595012
        11878942986520479340011031387572153845034397193548
        21112051579461663487198093927294533414503946670529
        27912914855801753338166857982316138938427044449697
        165870468846196415120486732660736

To put that number into perspective, Frontier, the world's most powerful
super-computer could run the x only up until roughly
    265000000000000000000000000000000000
before our Sun runs out of fuel, in 5 billion years.

---
Before we close out, all that's left for us to do now is put back the square
root to solve for `a`, and add the `x` brute force to the algorithm, and count
the steps.

The Fermat implementation also has a variable `count` that counts the number of
iterations of the while loop.

In Grant's variation, `x` will also conveniently count the number of steps for
us. Since we're checking for exit-condition inside the loop, we must make sure
to increment `x` before the if-statement to get correct results, and since we're
starting from 0, we also return x + 1 to get the steps.

This allows us to verify the performance is identical to that of
the original.
"""
from math import sqrt, ceil

def grant(B):
    b = sqrt(B)
    r = ceil(b) - b
    x = 0
    while True:
        a = sqrt((x + r) * (2 * b + (x + r)))
        c = sqrt(a ** 2 + B)
        x = x + 1
        PRIME1 = round(c - a)
        PRIME2 = round(c + a)
        if PRIME1 * PRIME2 == B:
            return PRIME1, PRIME2, x + 1

def fermat2(n):
    a = isqrt(n)
    b2 = a * a - n
    b = isqrt(n)
    count = 0
    while b * b != b2:
        a = a + 1
        b2 = a * a - n
        b = isqrt(b2)  # int(b2**0.5)
        count += 1
    p = a + b
    q = a - b
    return p, q, count + 1


if __name__ == '__main__':
    import random

    # The first 1000 primes to help us generate moduli below. Keep scrolling :)
    primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
              61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131,
              137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197,
              199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271,
              277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353,
              359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433,
              439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509,
              521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
              607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677,
              683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769,
              773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859,
              863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953,
              967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033,
              1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103,
              1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193,
              1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279,
              1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361,
              1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447,
              1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
              1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597,
              1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667,
              1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753,
              1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861,
              1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933,
              1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017,
              2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099,
              2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203,
              2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281,
              2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
              2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437,
              2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543,
              2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647,
              2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707,
              2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789,
              2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861,
              2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963,
              2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061,
              3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169,
              3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257,
              3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343,
              3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449,
              3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529,
              3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607,
              3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691,
              3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779,
              3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, 3877,
              3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947,
              3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051,
              4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139,
              4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, 4241,
              4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337,
              4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 4441,
              4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519,
              4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637,
              4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721,
              4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801,
              4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931,
              4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999,
              5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087,
              5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, 5189,
              5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, 5281, 5297,
              5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, 5407,
              5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471, 5477, 5479,
              5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569,
              5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659,
              5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749,
              5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849,
              5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927,
              5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053,
              6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, 6143,
              6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247,
              6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, 6311, 6317, 6323,
              6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379, 6389, 6397,
              6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521, 6529, 6547,
              6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637,
              6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709, 6719,
              6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827,
              6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911,
              6917, 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997,
              7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, 7109,
              7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213,
              7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, 7307, 7309, 7321,
              7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457,
              7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517, 7523, 7529, 7537,
              7541, 7547, 7549, 7559, 7561, 7573, 7577, 7583, 7589, 7591, 7603,
              7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691, 7699,
              7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817,
              7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907,
              7919]

    rng = random.SystemRandom()
    N   = rng.choice(primes) * rng.choice(primes)

    p1, q1, count1 = fermat2(N)
    p2, q2, count2 = grant(N)

    print(f"Fermat's algo:                   p={p1}, q={q1}, steps={count1}")
    print(f"Fermat's algo (Grant's variant): p={p2}, q={q2}, steps={count2}")

    # No matter how many times this code is executed, the
    # application will always yield the same number of steps.


"""
The plagiarism
--------------

The only difference between these two algorithms was that Grant went out of his
way to solve for leg `a`, instead of just increasing `c`, the hypotenuse, by one
at a time.

Fermat's method is much more readable and understandable. Grant never explained
how he got his weird formula to solve `a`:
    a = √( (x+r) (2b + (x+r)) )

The most likely answer is that he derived it the way we showed above, or one of
his grifter buddies, our guess being Nassim Haramein or Talal Ghannam, helped
him obfuscate Fermat's algorithm.

This is, again, exactly what Nassim Haramein has been doing, and considering
Grant and Haramein are acquainted from another new-age grift called Resonance
Science Foundation, it wouldn't be surprising if this was intentional,
collaborated obfuscated plagiarism.

Further evidence of this behavior is, that Grant also earlier, in 2019, hosted a
press conference where he used a camouflaged version of CADO-NFS, an open-source
factoring algorithm. This couldn't be even considered forking, because they did
not state they were using the General Number Field Sieve algorithm that CADO-NFS
uses, they claimed the software used their algorithm, called
    "CS BiPrime Algorithm (ver. 1.1.)".

Again, Grant never described where he got the formula to derive a with. If one
would argue he's so competent it just comes to him, then how was he also so
incompetent he forgot to include the outermost square root when solving `a`,
and why did he claim this algorithm runs in constant time when it famously runs
in 𝒪(√n) time? And why wasn't the algorithm ever released as code?

The answer should be obvious by now.
"""