import time
import numpy
from gmpy2 import is_square

def proofandcounting():
    if is_square(term2 - y4):
        global counter
        counter += 1
        global stringterm
        stringterm += (str(term) + "^2 = " + str(int((term2 - y4)**0.5)).rjust(11) +
                       "^2 + " + str(int(y4**0.25)).rjust(9) +
                       "^4    solution " + str(counter) + "\n")

timestart = time.time()
#█████████████ Mainparameter ██████████████████████████████
                                                        #██
begin  = 1                                              #██
howfar = 1000000   # 1000000 should be about 3 Minutes. #██
solutions_to_screen = [3, 4, 5, 6, 7, 8, 9, 10, 11, 12] #██
                                                        #██
#██████████████████████████████████████████████████████████
yyyy = (numpy.arange(0, int(pow(howfar * howfar, 0.25)) + 1, dtype=object))**4

for term in range (begin, howfar + 1): 
    term2 = term*term
    counter = 0
    stringterm = ""
    if is_square(term): ylimit = int(term ** 0.5)
    else: ylimit = int(term ** 0.5)+1
        
    if term2 % 24 == 1:                                     #  Some thoughts about accelerating the code: 
        for y4 in yyyy[1:ylimit]:      proofandcounting()   #  y^4 mod 24       (z^2 can only be [0,1,4,9,12,16] mod 24)
                                                            #     |
    elif term2 % 24 in [4, 16]:                             #  16 | 16 (17) (20)   1    4  ( 8)
        for y4 in yyyy[2:ylimit][::2]: proofandcounting()   #   9 |  9 (10) (13) (18) (21)   1  
                                                            #   1 |  1 ( 2) ( 5) (10) (13) (17)                 
    elif term2 % 24 == 9:                                   #   0 |  0   1    4    9   12   16                 
        for y4 in yyyy[3:ylimit][::3]: proofandcounting()   #  ---+------------------------------> x^2 mod 24  
                                                            #     |  0   1    4    9   12   16   
    else:   # ( = if term2 % 24 in [0, 12]:)
        for y4 in yyyy[6:ylimit][::6]: proofandcounting()   #  Because of x^2 + y^4 is an addition, we can add the residues mod 24 too.
                                                            #  Thus we have only to search in the non brackets ranges.
    if counter > 0:                                         #  (primes^2 allways have 1 mod 24, for all primes >= 5)
        f = open(str(counter) + "_Solutions.txt", "a") 
        if counter == 1: f.write(stringterm)
        if counter >  1: f.write(stringterm + "\n")
        f.close() 
        if counter in solutions_to_screen: print(stringterm)

print("Time = %3.3f " % (time.time() - timestart))