The Largest Semiprime Palindrome

# Problem:
#
# nnnnn * mmmmm = P
# nnnnn, mmmmm - five-digit primes
# P - the maximum palindrome
# (P - the largest semiprime palindrome)
# P - ?


# If the palindrome contains 10 digits, then it is divided by 11. That won't pass!
typealias TNumber_9   UInt32 # nine-digit unsigned numbers
typealias TNumber_5   UInt32 # five-digit ...

typealias TIndex      TNumber_5
typealias TIteration  TNumber_9

# signed version of numbers (for type casting)
typealias TINumber_5  Int32
typealias TINumber_9  Int32

typealias TIIteration Int32

type TProduct
    value::TNumber_9 # nnnnn * mmmmm
    index::TIndex
end

type TBitSet
    # BitArrays are space-efficient "packed" boolean arrays, which store one bit per boolean value.
    # One bit - one number (optimization)
    data::BitArray
    k::TIndex
end

const MAX_5  = 99989     # the largest five-digit prime number, but 99997 is good too
const MIN_5  = 10007     # the smallest ...,                    but 10001 is good too
const MAX_9  = 999999997 # 999,999,999 and 999,999,998 are not primes

SQRT_I = 2171            # primes[2171] == 31607 ~ sqrt(MAX_9) and isPrime(31607)==true
                         # (use eratosthenesSieve())


function eratosthenesSieve() ::Array{TNumber_5}
    # In cryptography we use sieve of Atkin, but this algorithm is the best for big numbers
    # The idea:
    # 1. sift composit numbers on 1..sqrt(MAX_5)
    # 2. sift numbers on MIN_5..MAX_5 that are divisible by primes on 1..sqrt(MAX_5)

    bitSet ::TBitSet = TBitSet(trues(MAX_5), 0) # bitSet.k == 0
    _end   ::TIndex  = MAX_5
    start  ::TIndex  = MIN_5

    # sift even numbers to the square root only (optimization)
    primaryEnd ::TIndex = TIndex(floor(sqrt(_end)))
    for i::TIndex= 4 : 2 : primaryEnd
        bitSet.data[i]= false
    end

    # sift even five-digit only (optimization)
    for i::TIndex= div(start, 2) * 2 : 2 : _end
        bitSet.data[i]=false
    end

    # sift odd five-digit
    for i::TIndex= 3 : 2 : primaryEnd # to the primaryEnd, not to the _end!
        if bitSet.data[i] # isPrime?
            # optimization: preliminary small sieve (from sqr(prime))
            for j::TIndex= i*i : i : primaryEnd
                bitSet.data[j]= false
            end
            # basic large sieve (from sqr(prime))
            for m::TIndex= div(start, i) * i : i : _end
                bitSet.data[m]= false
            end
        end
    end

    bitSet.data[1]= false # not a prime!
    bitSet.data[2]= true  # prime!

    # ñalculate the number of all five-digit primes
    # bitSet.k == 0
    oddStart ::TIndex = start + (mod(start, 2) == 0)
    for i= oddStart : 2 : _end
        if bitSet.data[i]      # isPrime?
            bitSet.k += 1      # the number of all five-digit primes
        end
    end

    # form primes into Array (for quick access)
    primes ::Array{TNumber_5} = zeros(Int, bitSet.k)
    m= 1
    for i= oddStart : 2 : _end
        if bitSet.data[i]      # isPrime?
            primes[m]= i
            m += 1
        end
    end

    return primes

end #eratosthenesSieve


function isPalindrome(number::TNumber_9) ::Bool

    tempNumber ::TNumber_9 = number;
    reversedNumber ::TNumber_9 = 0;
    rest::TNumber_9 = 0;

    # accumulate the reverse number
    while tempNumber != 0
        rest= mod(tempNumber, 10)
        tempNumber = div(tempNumber, 10)
        reversedNumber = reversedNumber*10 + rest
    end

    return (number == reversedNumber) # really a palindrome?

end #isPalindrome


function getProduct(primes::Array{TNumber_5}, diagram::Array{TIndex}, i::TIndex) ::TNumber_9
    return (primes[i] * primes[diagram[i]]) # nnnnn * mmmmm = PPPPPPPPP
end


function initDiagram(primes::Array{TNumber_5}) ::Array{TIndex}
    # nnnnn * mmmmm = PPPPPPPPP
    # form the structure (diagram) for searching for factors (nnnnn * mmmmm; MIN_5 <= nnnnn, mmmmm <= MAX_5)

    # nnnnn >= mmmmm
    # mmmmm <= sqrt(MAX_9), because: mmmmm > sqrt(MAX_9) => mmmmm > nnnnn
    # (optimization)
    diagram ::Array{TIndex} = zeros(Int, SQRT_I)
    j::TIndex= length(primes)
    for i::TIndex= 1:SQRT_I
        _max ::TNumber_5 = div(MAX_9, primes[i]) # mmmmm <= MAX_9 / nnnnn
        while j >= i
            if primes[j] <= _max
                diagram[i]= j
                break
            end
        j -= 1
        end
    end

    return diagram

end #initDiagram


function initSorted(primes::Array{TNumber_5}, diagram::Array{TIndex}) ::Array{TProduct}
    # transform diagram into the sorted structure (sorted::Array{TProduct})
    # (optimization)

    # initialization
    sorted ::Array{TProduct} = Array(TProduct, SQRT_I)
    for i::TIndex= 1:SQRT_I      # the system function fill() have problems with references:
                                 # https://github.com/JuliaLang/julia/issues/17796
        sorted[i]= TProduct(0,0)
    end
    for i::TIndex= 1:SQRT_I
        sorted[i].value= getProduct(primes, diagram, i)
        sorted[i].index= i
    end

    sort!(sorted; alg= QuickSort, lt= (x, y) -> x.value > y.value)

    return sorted

end #initSorted


function _sort!(sorted::Array{TProduct}, border::TIndex)
    # Insertion sort (one iteration) for Array.
    # Array is not good for optimization, but...
    # ...Julia does not support full-value pointers, full-value recursive data types and full-value references.
    # ...IntSet (sorted set; implemented as a bit string) have the problem with big numbers:
    #    https://docs.julialang.org/en/stable/stdlib/collections/?highlight=intset#Base.IntSet

    j::TIndex= 2
    key::TProduct= sorted[1]
    while (j <= border) && (key.value < sorted[j].value)
        sorted[j-1]= sorted[j] # shift all of sorted[j], where sorted[j].value < key.value (for all j)
        j += 1
    end
    sorted[j-1]= key

    return

end # _sort


function findPalindrome(N::TIteration) ::Dict{String, TINumber_9}

    # cascade inicialization
    primes  ::Array{TNumber_5} = eratosthenesSieve(                  ) # primes on MIN_5...MAX_5
    diagram ::Array{TIndex}    = initDiagram(        primes          ) # the structure for searching for factors
    sorted  ::Array{TProduct}  = initSorted(         primes, diagram ) # sorted diagram

    # main ñycle
    border ::TIndex = SQRT_I
    i::TIteration= 1
    while ((N == 0) && (border >= 1)) || (i <= N) # N iterations or to the end (N == 0)

        if isPalindrome(sorted[1].value)
            return (Dict([ # found!
                ("product"   , TINumber_9( sorted[1].value)),
                ("factor-1"  , TINumber_5( primes[sorted[1].index])),
                ("factor-2"  , TINumber_5( primes[diagram[sorted[1].index]])),
                ("iterations", TIIteration(i)) # numbers of iterations
            ]))
        end

        # delete the pair (nnnnn, mmmmm) from sorted (and diagram)
        # nnnnn * mmmmm == sorted[1].value == maximum at this iteration
        diagram[sorted[1].index] -= 1
        sorted[1].value= getProduct(primes, diagram, sorted[1].index)

        _sort!(sorted, border) # sort in place!

        # correct border (optimization)
        if primes[border] > primes[diagram[border]]
            border -= 1
            # println("border: ", border)
        end

        i += 1 # next iteration

    end # main ñycle

    return Dict()

end #findPalindrome


findPalindrome()= findPalindrome(TIteration(0)) # N == 0


@time result= findPalindrome()