#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 7
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.

Help := proc(): print(`Gsomosk(Ini, k, n), GenHomPol(vars, d, a), GenNonHomPol(vars, d, a)`):

#############
# Problem 1 #
#############

Gsomosk := proc(Ini, k, n) option remember:
    if n <= 0 then:
        return FAIL:
    elif n <= k then:
        return Ini[n]:
    else:
        return add(Gsomosk(Ini, k, n-i)*Gsomosk(Ini, k, n-k+i), i=1..floor(k/2))/Gsomosk(Ini, k, n-k):
    fi:
end:

#############
# Problem 2 #
#############

# Running

#       seq(hastype(Gsomosk([1$i], i, 2*i+5), integer), i=4..100);

# yields a list starting with 4 "true"s, followed by a bunch of
# "false"s. Generating the Somos k-sequences for k from 4 through 7
# yields sequences that appear to be all integers, so the experimental
# evidence suggests that the Somos k-sequence consists of all integers
# iff 4 <= k <= 7.

#############
# Problem 3 #
#############

# I was not able to gather an adequate amount of data as maple
# repeatedly hung while attempting the calculations. However, I was able to
# discover that for 4 <= k <= 7, the sequence appears to consist
# entirely of Laurent polynomials, while for k = 8 and k = 9, some
# terms of the sequence are not Laurent polynomials. I would
# conjecture that the sequence eventually contains a rational function
# that is not a Laurent polynomial for k >= 8.

#############
# Problem 4 #
#############

GenHomPol := proc(vars,d,a, start := 1) local m, firstvar, newstart, i:
    m := nops(vars):
    if m <= 1 then:
        return a[start]*vars[1]^d:
    fi:
    firstvar := vars[1]:
    newstart := k -> add((i+m-2)!/i!/(m-2)!, i=0..k-1):
    return add(firstvar^(d-i)*GenHomPol(i, vars[2..m], a, start+newstart(i)), i=0..d):
end:

#############
# Problem 5 #
#############

GenNonHomPol := proc(vars,d, a) local m, start:
    m := nops(vars):
    start := k -> add(binomial(i+m-1, i), i=0..k-1)+1:
    return add(GenHomPol(i, vars, a, start(i)), i=0..d):
end: