'''
Copyright (c) 2021 Serguei Zolotov

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as
published by the Free Software Foundation, either version 3 of the
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

'''
import numpy
import sympy
from fractions import Fraction 

UPPER_K = 64 #64

# generate list of primes
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n
def primes(n):
    # Returns a array of primes, 2 <= p < n
    m = numpy.ones(n // 2, dtype = numpy.bool)
    for i in range(3, int(n**0.5) + 1, 2):
        if m[i//2]:
            m[i*i//2::i] = False
    return [2] + list(2*numpy.nonzero(m)[0][1::]+1)

# when UPPER_K is large - it takes lots of system memory
# check your system configuration before changing UPPER_K
_known_primes_long = primes(2 ** (UPPER_K // 2) + 1)

# check if number is prime - direct method
def is_prime_s(n):
    
    if n & 1 == 0:
        return False

    d = 3
    while d * d <= n:
        if n % d == 0:
            return False
        d = d + 2

    return True

# check if number is prime - Miller-Rabin
def is_prime(n):

# https://www.ams.org/journals/mcom/1993-61-204/S0025-5718-1993-1192971-8/S0025-5718-1993-1192971-8.pdf
# https://oeis.org/A014233
    arr = []
    if n < 2047:
        arr = [2]
    elif n < 1373653: 
        arr = [2,3]
    elif n < 9080191:
        arr = [31,73]
    elif n < 25326001:
        arr = [2,3,5]
    elif n < 3215031751: 
        arr = [2,3,5,7]
    elif n < 4759123141:
        arr = [2,7,61]
    elif n < 1122004669633:
        arr = [2,13,23,1662803]
    elif n < 2152302898747: 
        arr = [2,3,5,7,11]
    elif n < 3474749660383:
        arr = [2,3,5,7,11,13]
    elif n < 341550071728321:
        arr = [2,3,5,7,11,13,17]
    elif n < 3825123056546413051:
        arr = [2,3,5,7,11,13,17,19,23]
    # elif n < 18446744073709551616:
    #     arr = [2,3,5,7,11,13,17,19,23,29,31,37]
    elif n < 318665857834031151167461:
        arr = [2,3,5,7,11,13,17,19,23,29,31,37]
    elif n < 3317044064679887385961981:
        arr = [2,3,5,7,11,13,17,19,23,29,31,37,41]
    else:
        return is_prime_s(n)

    # Miller-Rabin primality test
    d = n - 1
    s = 0
    while d % 2 == 0:
        d //= 2
        s += 1

    for a in arr:

        x = pow(a, d, n)
        if x == 1 or x == n - 1: 
            continue

        for _ in range(s - 1):
            x = pow(x, 2, n)
            if x == 1: 
                return False
            if x == n - 1: 
                break
        else: 
            return False

    return True

# find largest prime below N from list - binary search
def getMaxListPrime(N, nolist, PP):

    # edge case - last or above all
    if N > PP[-1]:
        # we do not know how much bigger it is
        return -1

    # edge case - first or below all
    if N < 2:
        # low case shall be handled already
        return -1

    left, right = 0, len(PP) - 1
    mid = 0

    while left < right:
        mid = (left + right) // 2  # find the mid
        if N < PP[mid]:
            right = mid
        elif N > PP[mid]:
            left = mid + 1
        else:
            break

    while (PP[mid] in nolist or PP[mid] > N) and mid > 0:
        mid = mid-1

    return PP[mid]

# find largest prime withon boundaries
def largestPrimeBetween(U, L, nolist):

    if U < 2:
        return -1
    if U == 2:
        # check if 2 was already used
        if U not in nolist:
            return 2
        else:
            return -1

    # odd prime low boundary     
    if U%2 == 0:
        U = U-1

    # get the prime candidate
    p = getMaxListPrime(U, nolist, _known_primes_long)
    if p in nolist:
        return -1

    if p > 1:
        # found it
        return p

    # list is too small
    while U >= L:
        # iterate down
        if is_prime(U) and U not in nolist:
            return U
        if U%3 == 2:
            U = U-4
        else:
            U = U-2

    return -1

# get largest number from within boundares with the given number of dividers
def getCandidate_s(U, L, K):

    for a in range(U, L, -1):
        if sympy.divisor_count(a) == K:
            return a
    return L

# compute candidate based on exponent array
def getCandidate(U, L, vexp, nolist):

    if len(vexp) == 1:
        # single prime
        e = vexp[0]
        rexp = Fraction(1, e)
        lprime = largestPrimeBetween(int(U**rexp), int(L**rexp), nolist)
        if lprime > 0:
            L = max(L, (lprime**e))

    else:
        prime1 = int(U**(1.0/(sum(vexp))))
        for p in _known_primes_long:
            if p > prime1:
                break            
            if p in nolist:
                continue       

            # recursive
            localnolist = nolist[:] + [p]
            for i, e in enumerate(vexp):
                pp = p**e
                # update exponents
                nvexp = vexp[:]
                nvexp.pop(i)
                m = getCandidate(U//pp, L//pp, nvexp, localnolist)
                if m > 0:
                    L = max(L, pp * m)

    return L

# return exponents up to 2 and flag if more exponents are possible
def computekd(K):
    kd = []
    kd.append([K - 1,])
    triples = False

    # add pairs
    for c in range(2, int(K ** 0.5) + 1):
        if K % c == 0:
            M = K // c
            kd.append([c - 1, M - 1])
            # check triplets
            if not triples:
                for b in range(c, int(M ** 0.5) + 1):
                    if M % b == 0:
                        triples = True
                        break

    return kd, triples

# compute largest number for all exponent arrays
def get_num(L, U, K):

    # get exponents
    exponents, domore = computekd(K)

    W = L
    # start searching
    for c in range(len(exponents)):
        vexp = exponents[c]
        W = max(W, getCandidate(U, W, vexp, []))

    if domore:
        W = max(W, getCandidate_s(U, W, K))
        # do not need to do the rest

    return W

# compute b-sequence
def compute(M):

    f8 = open("b341339.txt",'w')

    A = [[0 for i in range(M + 1)] for j in range(M + 1)]
    n = 2

    # antidiagonals
    for i in range(1, M + 1):
        for j in range(1, i + 1):
            x = i - j + 1
            y = j

            L = 0
            if x > 0:
                L = A[y][x - 1]            
            
            val = get_num(L, 2**x, y + 1)
            A[y][x] = int(val)

            print ("[%d:%d] = {%d}"%(x, y, val), end=" ", flush=" ")
            print ("%d %d"%(n, val), file=f8)
            f8.flush()
            n += 1

    f8.close()      

if __name__ == '__main__':
    # ...
    compute(UPPER_K)