project-euler / Lua / 196.lua

-- First we're going to need a quick way to calculate the values
-- of a row without having to iterate up to the very high rows
-- we're asked to calculate.  Looking at the examples provided,
-- let X be a lower-triangular matrix that satisfies the given
-- conditions:
--
--      r       X[r,0]
--
--      1          1
--      2          2
--      3          4
--      4          7
--      5         11
--      6         16
--      7         22
--      8         29
--
-- By looking at these examples, we can come up with a recurrence
-- to describe the relationship between r and X[r, 0].  As it
-- turns out:
--
--      X[r, 0] = X[r-1, 0] + r - 1
--
-- Or written perhaps more familiarly:
--
--      T(n) = T(n - 1) + n - 1
--
-- Which is equivalent to:
--
--      T(n) = (n - 1) + (n - 2) + (n - 3) + ... + 1
--                    n
--           = n^2 -  Σ i + 1
--                   i=1
--                    n^2 + n
--           = n^2 - --------- + 1
--                       2
--
-- So now we have an easy equation to find out what number is
-- in the first column at any row.  And of course there are the
-- same number of columns in any given row as that row number.
--
-- You should never check the siblings of the current value;
-- they will always be even (of course) and therefore not prime.
--
-- Most cases of triplets can be drawn using the immediate
-- neighbors, using the grid below:
--
--      1  2  3
--      4  5  6
--      7  8  9
--
-- Where the number we're checking is in cell 5.  If the row
-- is even, only cells 1 3 and 8 must be checked.  If the row
-- is odd, only cells 2 7 and 9 must be checked.  This will
-- indicate the central cell is part of any triplet involving
-- both the upper and lower rows.
--
-- A triplet can also occur between two elements on the same
-- row, like in row 8:
--
--      x [x] x
--     [x] x [x]
--
-- The above method of checking immediate neighbors won't
-- find this.
--
-- There are two cases where forks can occur to create prime
-- triplets, which for example, occur in row 10.  These will
-- be entirely above or below the current row.  There are two
-- possibilities.
--
--         α               β
--
--     [x] x [x]       [x] x [x]
--     [x] x [x]        x [x] x
--      x [x] x         x [x] x         <- row n
--      x [x] x        [x] x [x]
--     [x] x [x]       [x] x [x]
--
-- Case α occurs when n is even, case β occurs when n is odd.
-- This observations will halve the number of calculations.
-- If one of these forks is found, we can immediately quit
-- working on that value and add it to the sum.

package.path = package.path .. ";" .. "/home/taylor/Programs/Libraries/Lua/?.lua"
package.cpath = package.cpath .. ";" .. "/home/taylor/Programs/Libraries/Lua/?.so"
require "list"; require "numerology"; require "prime"; require "bignum"

function s(n)
    local function first(n)
        return n^2 - ((n^2 + n) / 2) + 1
    end
    local function last(n)
        return first(n) + n - 1
    end
    -- given a value, a row, and a direction
    -- indicate what the value in that direction is
    -- checks boundaries automatically:
    --   * lower, lower-right, lower-lower-right always OK
    --   * upper-left, upper-upper-left, lower-left OK when x > first(n)
    --   * upper OK when x < last(n)
    --   * upper-right OK when x < last(n) - 1
    --   * upper-upper-right OK when x < last(n) - 2
    local function direction(x, n, dir)
        if dir == "upper-upper-left" then
            if x > first(n) then return x - 2 * n + 2 else return 1 end

        elseif dir == "upper-upper-right" then
            if x < last(n) - 2 then return x - 2 * n + 4 else return 1 end

        elseif dir == "upper-left" then
            if x > first(n) then return x - n else return 1 end

        elseif dir == "upper-right" then
            if x < last(n) - 1 then return x - n + 2 else return 1 end

        elseif dir == "upper" then
            if x < last(n) then return x - n + 1 else return 1 end

        elseif dir == "lower" then
            return x + n

        elseif dir == "lower-left" then
            if x > first(n) then return x + n - 1 else return 1 end

        elseif dir == "lower-right" then
            return x + n + 1

        elseif dir == "lower-lower-left" then
            if x > first(n) then return x + 2 * n else return 1 end

        elseif dir == "lower-lower-right" then
            return x + 2 * n + 2

        else
            io.stderr:write("Unknown direction: " .. dir)
            os.exit(1)
        end
    end
    local primes = {}
    local function isprime(x)
        return Prime.isprime("fermat", x, 10)
    end
    local function maybe_print(s)
        if false then
            print(s)
        end
    end
    local sum = Bignum.new(0)
    local goalstep = math.floor(n / 100)
    local goal = first(n) + goalstep
    local stepnum = 1
    for x = first(n), last(n) do
        if x > goal then
            if math.mod(stepnum, 10) == 0 then
                io.stderr:write(" " .. stepnum .. " ")
            else
                io.stderr:write(".")
            end
            goal = goal + goalstep
            stepnum = stepnum + 1
        end
        if isprime(x) then
            maybe_print("=== Checking prime " .. x)
            local total = 0
            local stopnow = false
            if math.mod(n, 2) == 0 then
                if isprime(direction(x, n, "upper-left")) then
                    if isprime(direction(x, n, "upper-upper-left")) then
                        maybe_print("(" .. n .. ", " .. x .. "): Triplet with upper-left and upper-upper-left.")
                        stopnow = true
                    elseif x > first(n) + 1 and isprime(x - 2) then
                        maybe_print("(" .. n .. ", " .. x .. "): Triplet with left sibling via upper route.")
                        stopnow = true
                    else
                        maybe_print("(" .. n .. ", " .. x .. "): Upper-left is prime.")
                        total = total + 1
                    end
                end
                if not stopnow and isprime(direction(x, n, "upper-right")) then
                    if isprime(direction(x, n, "upper-upper-right")) then
                        maybe_print("(" .. n .. ", " .. x .. "): Triplet with upper-right and upper-upper-right.")
                        stopnow = true
                    elseif x < last(n) - 1 and isprime(x + 2) then
                        maybe_print("(" .. n .. ", " .. x .. "): Triplet with right sibling via upper route.")
                        stopnow = true
                    else
                        maybe_print("(" .. n .. ", " .. x .. "): Upper-right is prime.")
                        total = total + 1
                    end
                end
                if not stopnow and isprime(direction(x, n, "lower")) then
                    if isprime(direction(x, n, "lower-lower-right")) or
                       isprime(direction(x, n, "lower-lower-left")) then
                        maybe_print("(" .. n .. ", " .. x .. "): Triplet with lower and some branch.")
                        stopnow = true
                    else
                        maybe_print("(" .. n .. ", " .. x .. "): Lower is prime.")
                        total = total + 1
                    end
                end
            else
                if isprime(direction(x, n, "upper")) then
                    if isprime(direction(x, n, "upper-upper-right")) or
                       isprime(direction(x, n, "upper-upper-left")) then
                        maybe_print("(" .. n .. ", " .. x .. "): Triplet with upper and some branch.")
                        stopnow = true
                    else
                        maybe_print("(" .. n .. ", " .. x .. "): Upper is prime.")
                        total = total + 1
                    end
                end
                if not stopnow and isprime(direction(x, n, "lower-left")) then
                    if isprime(direction(x, n, "lower-lower-left")) then
                        maybe_print("(" .. n .. ", " .. x .. "): Triplet with lower-left and lower-lower-left.")
                        stopnow = true
                    elseif x > first(n) + 1 and isprime(x - 2) then
                        maybe_print("(" .. n .. ", " .. x .. "): Triplet with left sibling via lower route.")
                        stopnow = true
                    else
                        maybe_print("(" .. n .. ", " .. x .. "): Lower-left is prime.")
                        total = total + 1
                    end
                end
                if not stopnow and isprime(direction(x, n, "lower-right")) then
                    if isprime(direction(x, n, "lower-lower-right")) then
                        maybe_print("(" .. n .. ", " .. x .. "): Triplet with lower-right and lower-lower-right.")
                        stopnow = true
                    elseif x < last(n) - 1 and isprime(x + 2) then
                        maybe_print("(" .. n .. ", " .. x .. "): Triplet with right sibling via lower route.")
                        stopnow = true
                    else
                        maybe_print("(" .. n .. ", " .. x .. "): Lower-right is prime.")
                        total = total + 1
                    end
                end
            end
            if stopnow == true or total > 1 then
                sum = sum + Bignum.new(x)
                if stopnow then maybe_print(">>> Sum is now " .. tostring(sum) .. " by stopnow.")
                else maybe_print(">>> Sum is now " .. tostring(sum) .. " by total = " .. total) end
            end
        end
    end
    return sum
end

print("\n" .. tostring(s(arg[1])))
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.