# hw3.txt, 2014-02-02, Jeff Ames

# Problem 2
with(linalg):

# for d from 10 to 100 by 10 do
#     Digits := d:
#     ed := evalf(det(Hil(50))):
#     de := det(evalf(Hil(50))):
#     printf("%3d digits: %16e, %16e\n", d, ed, de):
# od:

# produces this output:
#
#      10 digits:   1.392616e-1466,   -5.096100e-455
#      20 digits:   1.392616e-1466,    3.967187e-784
#      30 digits:   1.392616e-1466,   5.671928e-1038
#      40 digits:   1.392616e-1466,   3.673654e-1225
#      50 digits:   1.392616e-1466,  -8.393330e-1353
#      60 digits:   1.392616e-1466,   2.708020e-1430
#      70 digits:   1.392616e-1466,  -2.304665e-1463
#      80 digits:   1.392616e-1466,   1.392615e-1466
#      90 digits:   1.392616e-1466,   1.392616e-1466
#     100 digits:   1.392616e-1466,   1.392616e-1466
#
# The values agree somewhere between Digits = 70 and Digits = 80.

# Investigating this range more closely:

# for d from 70 to 80 do
#     Digits := d:
#     ed := evalf(det(Hil(50))):
#     de := det(evalf(Hil(50))):
#     printf("%2d digits: %16e, %16e\n", d, ed, de):
# od:

# gives these results:
#
#    70 digits:   1.392616e-1466,  -2.304665e-1463
#    71 digits:   1.392616e-1466,  -8.337626e-1465
#    72 digits:   1.392616e-1466,   9.103217e-1466
#    73 digits:   1.392616e-1466,  -1.709360e-1467
#    74 digits:   1.392616e-1466,   1.468210e-1466
#    75 digits:   1.392616e-1466,   1.379502e-1466
#    76 digits:   1.392616e-1466,   1.390410e-1466
#    77 digits:   1.392616e-1466,   1.392531e-1466
#    78 digits:   1.392616e-1466,   1.392612e-1466
#    79 digits:   1.392616e-1466,   1.392616e-1466
#    80 digits:   1.392616e-1466,   1.392615e-1466
#
# showing agreement at around Digits = 74.

# Generalizing this a bit, albeit inefficiently:

# evalDetDiff := proc(n, perc)
# local d, de, ed:
#     for d from 10 to 500 do
#         Digits := d:
#         ed := evalf(det(Hil(n))):
#         de := det(evalf(Hil(n))):
#         if abs((ed - de) / de) < perc then
#             return d:
#         fi:
#     od:
#     return -1:
# end:

# for n from 10 to 100 by 10 do
#     printf("%d: %d\n", n, evalDetDiff(n, 0.10)):
# od:

# produces
#
#    10: 14
#    20: 27
#    30: 44
#    40: 59
#    50: 74
#    60: 90
#    70: 104
#    80: 119
#    90: 134
#    100: 149
#
# which appears to be approximately (3/2 * n).

# Problem 3
inv1 := proc(p)
local i, j, n, ninv:
    ninv := 0:
    n := nops(p):
    for i from 1 to n do
        for j from i+1 to n do
            if p[i] > p[j] then
                ninv := ninv + 1:
            fi:
        od:
    od:
    ninv:
end:

gf := proc(n, q)
local sum, p:
   sum := 0:
   for p in permute(n) do
       sum := sum + q^inv1(p)
   od:
   sum:
end:

for n from 1 to 8 do
    gf(n, q);
od;

# produces
#
# 1
# 1+q
# q^3+2*q^2+2*q+1
# q^6+3*q^5+5*q^4+6*q^3+5*q^2+3*q+1
# q^10+4*q^9+9*q^8+15*q^7+20*q^6+22*q^5+20*q^4+15*q^3+9*q^2+4*q+1
# q^15+5*q^14+14*q^13+29*q^12+49*q^11+71*q^10+90*q^9+101*q^8+101*q^7+90*q^6+71*q^5+49*q^4+29*q^3+14*q^2+5*q+1
# q^21+6*q^20+20*q^19+49*q^18+98*q^17+169*q^16+259*q^15+359*q^14+455*q^13+531*q^12+573*q^11+573*q^10+531*q^9+455*q^8+359*q^7+259*q^6+169*q^5+98*q^4+49*q^3+20*q^2+6*q+1
# q^28+7*q^27+27*q^26+76*q^25+174*q^24+343*q^23+602*q^22+961*q^21+1415*q^20+1940*q^19+2493*q^18+3017*q^17+3450*q^16+3736*q^15+3836*q^14+3736*q^13+3450*q^12+3017*q^11+2493*q^10+1940*q^9+1415*q^8+961*q^7+602*q^6+343*q^5+174*q^4+76*q^3+27*q^2+7*q+1
#
# cofficients only:
#
# 1
# 1, 1
# 1, 2, 2, 1
# 1, 3, 5, 6, 5, 3, 1
# 1, 4, 9, 15, 20, 22, 20, 15, 9, 4, 1
# 1, 5, 14, 29, 49, 71, 90, 101, 101, 90, 71, 49, 29, 14, 5, 1
# 1, 6, 20, 49, 98, 169, 259, 359, 455, 531, 573, 573, 531, 455, 359, 259, 169, 98, 49, 20, 6, 1
# 1, 7, 27, 76, 174, 343, 602, 961, 1415, 1940, 2493, 3017, 3450, 3736, 3836, 3736, 3450, 3017, 2493, 1940, 1415, 961, 602, 343, 174, 76, 27, 7, 1
#
# N(r,c) = N(r-1, c) + N(r, c-1) almost works here, but fails for some central elements (r for row, c for column)