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

#Read procedures from class
read(`C17.txt`):

Help:=proc():
    print(` SIrec(L,Ini,f,m,n) , dIBVP(L,Ini,Fini,f,m,N) , GR(p,N) , LGR(p,N) `):
end:

##PROBLEM 1##

#SIrec(L,Ini,f,m,n): same input as Srec but in addition an
#expression, f, in the discerte variable m, that solves
#Inhomogeneous linear recurrences 
#a(n)=L[1]*a(n-1)+L[2]*a(n-2)+ ... + L[-1]*a(n-nops(L))+f(n)
SIrec:=proc(L, Ini, f, m, n) local i:
option remember:
    if nops(L) <> nops(Ini) then
        print(Ini, L, `should be lists of the SAME length`):
        return FAIL:
    elif n < 0 then
        return 0:
    elif n<nops(Ini) then
        return Ini[n+1]:
    else
        return subs(m=n, f) +
            add(L[i]*SIrec(L,Ini,f,m,n-i), i=1..nops(L)):
    fi:
end:

##PROBLEM 2##

#dIBVP(L,Ini,Fini,f,m,N): the analogous thing to 
#dBVP(L,Ini,Fini,N):
dIBVP:=proc(L,Ini,Fini,f,m,N) local eq, var, a, i:
    if nops(L) <> nops(Ini) + nops(Fini) then
        return FAIL:
    fi:
    
    var:={seq(a[i], i=0..N)}:
    eq:={seq(a[i]=Ini[i+1], i=0..nops(Ini)-1), seq(a[N-i]=Fini[i+1], i=0..nops(Fini)-1), seq(a[i]-add(L[j]*a[i-j],j=1..nops(L))-subs(m=n,f), i=nops(L)..N)}:
    
    var:=solve(eq, var):
    if var = NULL then
        return FAIL:
    else
        return subs(var, [seq(a[i], i=0..N)]):
    fi:
end:

##PROBLEM 3##

#GR(p,N): inputs a real number p, between 0 and 1 (or symbolic),
#and a positive integer N, and outputting a list of length N+1,
#let's call it G, such G[i+1] is your chance of exiting a winner
#in a casino with the following rule:
#if you enter a casino where at each step you win a dollar with
#probability p and lose a dollar with probability 1-p, and you
#stop playing when you are either broke (0 dollars),
#or got N dollars
GR:=proc(p,N)
    return dBVP([1/p, (p-1)/p], [0], [1], N):
end:

##PROBLEM 4##

#LGR(p,N): inputs the same thing as above, and outputs the
#EXPECTED duration of the game until you end it
#(either winner or loser).
LGR:=proc(p,N) local m:
    return dIBVP([1/p, (p-1)/p], [0], [0], -1/p, m, N):
end:

##PROBLEM 5##

#Conjecture: Probability of winning with i dollars is i/N
#Conjecture: Expected duration with i dollars is i*(N-i)