#Nathan Fox
#Homework 25
#I give permission for this work to be posted online

#Read procedures from class/last homework
read(`C25.txt`):
read(`C6.txt`):

#Help procedure
Help:=proc():
    print(` Vg(m, p, MaxLoss) , Vxg(m, p, x, MaxLoss) `):
    print(` GuessPol1(L, n, d) , GuessPol(L, n) `):
    print(` BdP(m0, p) , Cpoly(m, p) , ProveBdP(M0, p) `):
end:

##PROBLEM 1##

#Conjectures:
#Expectation: V(n, n), which is asymptotically c*sqrt(n), where c is approximately 0.521
#Variance: asymptotically c*n, where c is approximately 0.145
#The standardized moments seems to decrease, and the odd ones at least converge to zero.  But, it's not asymptoticlaly normal.

##PROBLEM 2##

#Vg(m, p, MaxLoss): computes the expected gain where the player
#refuses to ever have m-p more than MaxLoss (in other words, he or
#she do not want to take a chance, however small, of losing more
#than MaxLoss dollars).
Vg:=proc(m, p, MaxLoss)
option remember:
    if m = 0 then
        return p:
    elif p = 0 then
        return 0:
    elif m - p > MaxLoss then
        return 0:
    else
        return max(0, m/(m+p)*(-1+Vg(m-1, p, MaxLoss))+p/(m+p)*(1+Vg(m, p-1, MaxLoss))):
    fi:
end:

#The sequences evalf([seq(Vg(i,i,MaxLoss),i=1..100)]), for MaxLoss from 1 to 5? are
#[.5000000000, .6666666667, .8500000000, 1., 1.119047619, 1.214285714, 1.291666667, 1.355555556, 1.409090909, 1.454545455, 1.493589744, 1.527472527, 1.557142857, 1.583333333, 1.606617647, 1.627450980, 1.646198830, 1.663157895, 1.678571429, 1.692640693, 1.705533597, 1.717391304, 1.728333333, 1.738461538, 1.747863248, 1.756613757, 1.764778325, 1.772413793, 1.779569892, 1.786290323, 1.792613636, 1.798573975, 1.804201681, 1.809523810, 1.814564565, 1.819345661, 1.823886640, 1.828205128, 1.832317073, 1.836236934, 1.839977852, 1.843551797, 1.846969697, 1.850241546, 1.853376503, 1.856382979, 1.859268707, 1.862040816, 1.864705882, 1.867269985, 1.869738752, 1.872117400, 1.874410774, 1.876623377, 1.878759398, 1.880822747, 1.882817066, 1.884745763, 1.886612022, 1.888418826, 1.890168971, 1.891865079, 1.893509615, 1.895104895, 1.896653098, 1.898156277, 1.899616368, 1.901035197, 1.902414487, 1.903755869, 1.905060883, 1.906330989, 1.907567568, 1.908771930, 1.909945318, 1.911088911, 1.912203830, 1.913291139, 1.914351852, 1.915386932, 1.916397296, 1.917383821, 1.918347339, 1.919288646, 1.920208500, 1.921107628, 1.921986721, 1.922846442, 1.923687424, 1.924510272, 1.925315568, 1.926103866, 1.926875700, 1.927631579, 1.928371993, 1.929097412, 1.929808287, 1.930505051, 1.931188119, 1.931857892]
#[.5000000000, .6666666667, .8500000000, 1., 1.119047619, 1.229437229, 1.335372960, 1.434110334, 1.524701769, 1.607227911, 1.682220257, 1.750369431, 1.812388035, 1.868949986, 1.920666995, 1.968082757, 2.011675258, 2.051862478, 2.089009160, 2.123433558, 2.155413673, 2.185192797, 2.212984351, 2.238976036, 2.263333399, 2.286202856, 2.307714279, 2.327983185, 2.347112615, 2.365194727, 2.382312166, 2.398539232, 2.413942892, 2.428583646, 2.442516276, 2.455790501, 2.468451535, 2.480540586, 2.492095281, 2.503150044, 2.513736429, 2.523883406, 2.533617622, 2.542963624, 2.551944063, 2.560579872, 2.568890424, 2.576893673, 2.584606284, 2.592043742, 2.599220457, 2.606149854, 2.612844456, 2.619315957, 2.625575292, 2.631632692, 2.637497745, 2.643179444, 2.648686231, 2.654026038, 2.659206327, 2.664234124, 2.669116048, 2.673858342, 2.678466898, 2.682947283, 2.687304759, 2.691544305, 2.695670636, 2.699688217, 2.703601285, 2.707413857, 2.711129750, 2.714752588, 2.718285818, 2.721732719, 2.725096411, 2.728379868, 2.731585921, 2.734717272, 2.737776499, 2.740766061, 2.743688308, 2.746545485, 2.749339740, 2.752073125, 2.754747606, 2.757365064, 2.759927301, 2.762436046, 2.764892953, 2.767299612, 2.769657547, 2.771968221, 2.774233039, 2.776453352, 2.778630458, 2.780765605, 2.782859992, 2.784914776]
#[.5000000000, .6666666667, .8500000000, 1., 1.119047619, 1.229437229, 1.335372960, 1.434110334, 1.524701769, 1.607227911, 1.687397226, 1.765411833, 1.840809473, 1.913207824, 1.982391061, 2.048281891, 2.110900672, 2.170331558, 2.226697835, 2.280144834, 2.330828399, 2.378907283, 2.424538231, 2.467872919, 2.509056134, 2.548224797, 2.585507538, 2.621024642, 2.654888231, 2.687202580, 2.718064529, 2.747563930, 2.775784113, 2.802802347, 2.828690287, 2.853514407, 2.877336396, 2.900213539, 2.922199064, 2.943342465, 2.963689801, 2.983283969, 3.002164951, 3.020370049, 3.037934090, 3.054889622, 3.071267087, 3.087094985, 3.102400017, 3.117207224, 3.131540106, 3.145420740, 3.158869878, 3.171907048, 3.184550634, 3.196817965, 3.208725379, 3.220288300, 3.231521292, 3.242438124, 3.253051818, 3.263374700, 3.273418443, 3.283194114, 3.292712205, 3.301982677, 3.311014988, 3.319818124, 3.328400630, 3.336770634, 3.344935876, 3.352903724, 3.360681203, 3.368275008, 3.375691529, 3.382936863, 3.390016834, 3.396937006, 3.403702698, 3.410318996, 3.416790769, 3.423122675, 3.429319178, 3.435384553, 3.441322899, 3.447138146, 3.452834066, 3.458414277, 3.463882254, 3.469241336, 3.474494730, 3.479645519, 3.484696668, 3.489651030, 3.494511350, 3.499280270, 3.503960336, 3.508553999, 3.513063622, 3.517491482]
#[.5000000000, .6666666667, .8500000000, 1., 1.119047619, 1.229437229, 1.335372960, 1.434110334, 1.524701769, 1.607227911, 1.687397226, 1.765411833, 1.840809473, 1.913207824, 1.982391061, 2.048281891, 2.112213746, 2.174760756, 2.235998731, 2.295858627, 2.354246421, 2.411081973, 2.466309395, 2.519897229, 2.571835246, 2.622130460, 2.670803326, 2.717884456, 2.763411894, 2.807428937, 2.849982397, 2.891121250, 2.930895589, 2.969355819, 3.006552052, 3.042533652, 3.077348905, 3.111044779, 3.143666758, 3.175258741, 3.205862971, 3.235520014, 3.264268750, 3.292146395, 3.319188528, 3.345429133, 3.370900650, 3.395634028, 3.419658781, 3.443003048, 3.465693654, 3.487756166, 3.509214955, 3.530093252, 3.550413200, 3.570195915, 3.589461531, 3.608229252, 3.626517400, 3.644343458, 3.661724116, 3.678675309, 3.695212257, 3.711349500, 3.727100938, 3.742479855, 3.757498960, 3.772170410, 3.786505837, 3.800516382, 3.814212711, 3.827605044, 3.840703175, 3.853516493, 3.866054005, 3.878324348, 3.890335815, 3.902096365, 3.913613644, 3.924894993, 3.935947472, 3.946777865, 3.957392697, 3.967798243, 3.978000544, 3.988005413, 3.997818448, 4.007445040, 4.016890386, 4.026159491, 4.035257182, 4.044188115, 4.052956779, 4.061567508, 4.070024483, 4.078331741, 4.086493183, 4.094512575, 4.102393556, 4.110139645]
#[.5000000000, .6666666667, .8500000000, 1., 1.119047619, 1.229437229, 1.335372960, 1.434110334, 1.524701769, 1.607227911, 1.687397226, 1.765411833, 1.840809473, 1.913207824, 1.982391061, 2.048281891, 2.112213746, 2.174760756, 2.235998731, 2.295858627, 2.354246421, 2.411081973, 2.466309395, 2.520043731, 2.572700079, 2.624475549, 2.675437302, 2.725591031, 2.774916777, 2.823386557, 2.870972791, 2.917652079, 2.963406640, 3.008224575, 3.052099558, 3.095030272, 3.137019764, 3.178074787, 3.218205195, 3.257423387, 3.295743822, 3.333182592, 3.369757062, 3.405485554, 3.440387088, 3.474481160, 3.507787555, 3.540326190, 3.572116994, 3.603179795, 3.633534235, 3.663199706, 3.692195287, 3.720539704, 3.748251294, 3.775347980, 3.801847251, 3.827766150, 3.853121263, 3.877928716, 3.902204178, 3.925962853, 3.949219494, 3.971988404, 3.994283441, 4.016118033, 4.037505180, 4.058457469, 4.078987083, 4.099105812, 4.118825062, 4.138155870, 4.157108912, 4.175694519, 4.193922679, 4.211803058, 4.229345007, 4.246557569, 4.263449495, 4.280029252, 4.296305031, 4.312284760, 4.327976111, 4.343386509, 4.358523142, 4.373392971, 4.388002731, 4.402358950, 4.416467946, 4.430335841, 4.443968565, 4.457371867, 4.470551313, 4.483512304, 4.496260071, 4.508799688, 4.521136075, 4.533274005, 4.545218106, 4.556972870]
#These all start out with initial segments of evalf([seq(V(i,i),i=1..100)]), but then they start being less than the latter
#This is as expected

##PROBLEM 3##

#Vxg(m, p, x, MaxLoss): finds the probability generating funcion
#with the above restiction
Vxg:=proc(m, p, x, MaxLoss)
option remember:
    if m = 0 then
        return x^p:
    elif Vg(m, p, MaxLoss) = 0 then
        return 1:
    else
        return expand(m/(m+p)*(Vxg(m-1, p, x, MaxLoss)/x)+p/(m+p)*(Vxg(m, p-1, x, MaxLoss)*x)):
    fi:
end:

##PROBLEM 4##

#GuessPol1(L, n, d): inputs a list L, a symbol n, and a positive
#integer d (<=nops(L)-3), and outputs a polynomial of degree d in
#n, let's call it POL, such that L[i]=POL(i)
#(for i=1,2, ..., nops(L)). If it fails, it returns FAIL.
GuessPol1:=proc(L, n, d) local POL, a, eq, var, i:
    if d > nops(L)-3 then
        return FAIL:
    fi:
    POL:=add(a[i]*n^i, i=0..d):
    var:={seq(a[i], i=0..d)}:
    eq:={seq(subs(n=i, POL)=L[i], i=1..nops(L))}:
    var:=solve(eq, var):
    if var = NULL then
        return FAIL:
    fi:
    return subs(var, POL):
end:

##PROBLEM 5##

#GuessPol(L, n): inputs a list of numbers L, and outputs a
#polynomial of degree <=nops(L)-3 that interpolates it,
#if it exits, or returns FAIL.
GuessPol:=proc(L, n) local P, d:
    for d from 0 to nops(L) - 3 do
        P:=GuessPol1(L, n, d):
        if P <> FAIL then
            return P:
        fi:
    od:
    return FAIL:
end:

##PROBLEM 6##

#NOTE: I switched the arguments for my own sanity

#BdP(m0, p): inputs a SYMBOL p, a positive integer m0, and outputs
#the smallest integer, p0, such that Bd(m0, p0) is non-zero, and
#conjectures an explicit polynomial expressions for Bd(m0, p)
#(valid for p>=p0).
BdP:=proc(m0, p) local p0, i, P, L, last:
    for p0 from 1 while Bd(m0, p0) = 0 do od:
    last:=2*p0+1:
    L:=[seq(Bd(m0, i), i=p0..last)]:
    while true do
        P:=GuessPol(L, p):
        if P <> FAIL then
            return m0, sort(simplify(subs(p=p-p0+1, P))):
        fi:
        last:=2*last+1:
        L:=[seq(Bd(m0, i), i=p0..last)]:
    od:
end:

#BdP(m0, p) for m0 from 0 to 20:
#0: p
#1: p^2
#2: 1/2*p^3+1/2*p^2-p
#3: 1/6*p^4+1/2*p^3-2/3*p^2-2*p+2
#4: 1/24*p^5+1/4*p^4-1/24*p^3-9/4*p^2-p+6
#5: 1/120*p^6+1/12*p^5+1/8*p^4-13/12*p^3-47/15*p^2+3*p+12
#6: 1/720*p^7+1/48*p^6+11/144*p^5-13/48*p^4-101/45*p^3-9/4*p^2+32/3*p+34
#7: 1/5040*p^8+1/240*p^7+19/720*p^6-1/48*p^5-319/360*p^4-179/60*p^3+571/420*p^2+73/2*p+82
#8: 1/40320*p^9+1/1440*p^8+19/2880*p^7+1/90*p^6-1297/5760*p^5-2261/1440*p^4-27029/10080*p^3+1867/120*p^2+979/10*p+133
#9: 1/362880*p^10+1/10080*p^9+79/60480*p^8+1/180*p^7-127/3456*p^6-751/1440*p^5-189961/90720*p^4+1039/630*p^3+66943/1260*p^2+2728/15*p+208
#10: 1/3628800*p^11+1/80640*p^10+13/60480*p^9+61/40320*p^8-3209/1209600*p^7-1417/11520*p^6-62525/72576*p^5-31943/20160*p^4+798583/50400*p^3+96353/840*p^2+32027/105*p+306
#11: 1/39916800*p^12+1/725760*p^11+11/362880*p^10+37/120960*p^9+23/44800*p^8-1033/48384*p^7-8891/36288*p^6-19015/18144*p^5+1704049/907200*p^4+3033/70*p^3+11662697/55440*p^2+199391/420*p+4082
#12: 1/479001600*p^13+1/7257600*p^12+163/43545600*p^11+73/1451520*p^10+1157/4838400*p^9-6379/2419200*p^8-2290511/43545600*p^7-538037/1451520*p^6-6485161/10886400*p^5+18499759/1814400*p^4+56015471/623700*p^3+1762373/5040*p^2+1370671/315*p-2564
#13: 1/6227020800*p^14+1/79833600*p^13+197/479001600*p^12+17/2419200*p^11+53/967680*p^10-143/806400*p^9-55037/6220800*p^8-97121/1036800*p^7-4400449/10886400*p^6+343459/302400*p^5+64745231/2494800*p^4+269384497/1663200*p^3+1067098453/450450*p^2-70037/210*p+4947
#14: 1/87178291200*p^15+1/958003200*p^14+13/319334400*p^13+821/958003200*p^12+8969/958003200*p^11+379/29030400*p^10-715151/609638400*p^9-1606907/87091200*p^8-26263/201600*p^7-2060489/10886400*p^6+150662417/29937600*p^5+116144783/2217600*p^4+132596685211/151351200*p^3+292122679/277200*p^2+23890079/4620*p+3271
#15: 1/1307674368000*p^16+1/12454041600*p^15+137/37362124800*p^14+89/958003200*p^13+18967/14370048000*p^12+307/45619200*p^11-215833/1828915200*p^10-256463/87091200*p^9-19639561/653184000*p^8-72517/544320*p^7+365435809/718502400*p^6+69924743/5702400*p^5+6688607831513/27243216000*p^4+104203216213/129729600*p^3+17981100701/5405400*p^2+11933737/1980*p+4845
#16: 1/20922789888000*p^17+1/174356582400*p^16+317/1046139494400*p^15+113/12454041600*p^14+239693/1494484992000*p^13+163/119750400*p^12-596549/80472268800*p^11-236017/609638400*p^10-803720849/146313216000*p^9-6891929/174182400*p^8-307994521/5748019200*p^7+241720883/119750400*p^6+23859303237923/435891456000*p^5+47460682337/145297152*p^4+2884742608217/1816214400*p^3+35112583021/7207200*p^2+3027562793/360360*p+6700
#17: 1/355687428096000*p^18+1/2615348736000*p^17+11/475517952000*p^16+211/261534873600*p^15+180203/10461394944000*p^14+18979/93405312000*p^13+6751/44706816000*p^12-841117/20118067200*p^11-121301057/146313216000*p^10-153992939/18289152000*p^9-7809355847/201180672000*p^8+262792477/1437004800*p^7+1171245612647/118879488000*p^6+1264486165367/13621608000*p^5+293074114013/504504000*p^4+1145346626771/454053600*p^3+1506103645643/214414200*p^2+2107836457/180180*p+3072991
#18: 1/6402373705728000*p^19+1/41845579776000*p^18+103/62768369664000*p^17+53/804722688000*p^16+9487/5706215424000*p^15+526531/20922789888000*p^14+1548373/11769069312000*p^13-960023/268240896000*p^12-1019379761/9656672256000*p^11-419599489/292626432000*p^10-51517395767/4828336128000*p^9-5720054173/402361344000*p^8+413704031449/294226732800*p^7+52842580206329/2615348736000*p^6+19794544808329/118879488000*p^5+2005829228693/2095632000*p^4+585124633287629/154378224000*p^3+75517916401/7567560*p^2+27747615437/9009*p-14948582
#19: 1/121645100408832000*p^20+1/711374856192000*p^19+29/266765571072000*p^18+1/201180672000*p^17+9161/62768369664000*p^16+56699/20922789888000*p^15+2311259/94152554496000*p^14-731957/3487131648000*p^13-80261/7023034368*p^12-658181899/3218890752000*p^11-5080539361/2414168064000*p^10-8218731403/804722688000*p^9+13702528728101/94152554496000*p^8+9144236592973/2615348736000*p^7+50107690063019/1307674368000*p^6+4754861633917/16765056000*p^5+5529033779033227/3705077376000*p^4+94898477076223/17153136000*p^3+13740753250887701/8888443200*p^2-1610443978817/120120*p+30799051
#20: 1/2432902008176640000*p^21+1/12804747411456000*p^20+173/25609494822912000*p^19+1489/4268249137152000*p^18+83663/7113748561920000*p^17+32663/125536739328000*p^16+2473781/753220435968000*p^15-403829/376610217984000*p^14-7874887753/7532204359680000*p^13-43869943/1755758592000*p^12-95905879/286123622400*p^11-50202587263/19313344512000*p^10+1505580598261/235381386240000*p^9+5679995605307/11769069312000*p^8+679494534920413/94152554496000*p^7+118648732749737/1743565824000*p^6+50911381192082627/111152321280000*p^5+8311774500415951/3705077376000*p^4+1521022368345217841/2933186256000*p^3-2773741327715497/467812800*p^2+149190447259853/6126120*p-26851891

##PROBLEM 7##

#Need binomial coefficients as polynomials
Cpoly:=proc(m, p) local i:
    #return binomial(m+p, m):
    return expand(product(i, i=1..m+p)/(product(i, i=1..m)*product(i, i=1..p))):
end:

#Again, I switched the arguments intentionally

#ProveBdP(M0, p): first uses the above BdP to GUESS polynomial
#equations for Bd(p, m) and m0<=M0, and then goes on to prove
#these guesses rigorously (by induction), by using the recurrence
#(10) in William M. Boyce's article "On a simple optimal
#stopping problem."
#
#Returns the guesses, the polynomials from the recurrence, and
#either true or false (hopefully true) if the two match
ProveBdP:=proc(M0, p) local m0, PL, PL2, nxt:#p0L:
    #p0L:=[]:
    PL:=[p]:
    PL2:=[p]:
    for m0 from 1 to M0 do
        nxt:=BdP(m0, p):
        #p0L:=[op(m0L), nxt[1]]:
        PL:=[op(PL), nxt[2]]:
        PL2:=[op(PL2), simplify(Cpoly(m0, p-1)-Cpoly(m0-1, p)+PL2[-1]+subs(p=p-1, PL[-1]))]:
    od:
    return PL, PL2, evalb({op(PL2-PL)}={0}):
end:

#It successfully proved up to 20