with(gfun):

# The expected duration of the risk-adverse gambler drawing from an (m,
# p)-Shepp urn.
# (The risk-adverse gambler stops drawing as soon as he obtains a value of g,
# or if he cannot improve his value.)
E := proc(m, p)
    ExpD(m, p, p - m):
end:

# The expected duration of the risk-adverse gambler drawing from Shepp's urn
# given that his goal is g and his current value is 0.
# (The risk-adverse gambler stops drawing as soon as he obtains a value of g,
# or if he cannot improve his value.)
ExpD := proc(m, p, g)
    option remember:
    if g <= 0 then
        # Gamber is done.
        0
    elif m = 0 then
        # Gambler will draw until he reaches his goal.
        min(g, p)
    elif p = 0 then
        # Gambler cannot possibly improve his standing.
        0
    else
        # Gambler draws randomly and evaluates his position.
        m / (m + p) * ExpD(m - 1, p, g + 1) + p / (m + p) * ExpD(m, p - 1, g - 1) + 1:
    fi:
end:

# Guess the closed form for E(m, m + n) as a sequence in m, using N + 1
# datapoints.
# Try: guess_form_k(20, 1).
#      guess_form_k(20, 5).
guess_form_m := (N, k) -> rsolve(listtorec([seq(E(m, m + k), m=0..N)], a(m))[1], a(m)):

# Example outputs:

(*

                                               1/2
                                             Pi    GAMMA(m + 2)
                                        -1 + ------------------
                                               GAMMA(m + 3/2)

                                          1/2
                              2 (m + 2) Pi    GAMMA(m + 1)   2 (2 m + 3)
                              ---------------------------- - -----------
                                     GAMMA(m + 3/2)             m + 1

                                             1/2
                           3 (3 m + 7)   3 Pi    GAMMA(m + 1) (m + 3) (m + 2)
                         - ----------- + ------------------------------------
                              m + 1                 GAMMA(m + 5/2)

                          1/2   2                                  2
                      4 Pi    (m  + 7 m + 12) GAMMA(m + 1)   8 (2 m  + 11 m + 15)
                      ------------------------------------ - --------------------
                                 GAMMA(m + 5/2)                (m + 1) (m + 2)

                         2                    1/2                        2
                   5 (5 m  + 35 m + 62)   5 Pi    GAMMA(m + 1) (m + 3) (m  + 9 m + 20)
                 - -------------------- + --------------------------------------------
                     (m + 1) (m + 2)                     GAMMA(m + 7/2)

*)

# The term without the Gamma fucntions is O(1). I conjecture that the other
# term dominates.

# Grab the coefficient on sqrt(Pi) (the conjectured dominating term):

for k from 1 to 10 do coeff(guess_form_m(15, k), sqrt(Pi)); od;

(* Outputs:

                                             GAMMA(m + 2)
                                            --------------
                                            GAMMA(m + 3/2)

                                        2 GAMMA(m + 1) (m + 2)
                                        ----------------------
                                            GAMMA(m + 3/2)

                                    3 GAMMA(m + 1) (m + 3) (m + 2)
                                    ------------------------------
                                            GAMMA(m + 5/2)

                                                     2
                                    4 GAMMA(m + 1) (m  + 7 m + 12)
                                    ------------------------------
                                            GAMMA(m + 5/2)

                                                         2
                                5 GAMMA(m + 1) (m + 3) (m  + 9 m + 20)
                                --------------------------------------
                                            GAMMA(m + 7/2)

                                                3       2
                               6 GAMMA(m + 1) (m  + 15 m  + 74 m + 120)
                               ----------------------------------------
                                            GAMMA(m + 7/2)

                                                    3       2
                           7 GAMMA(m + 1) (m + 4) (m  + 18 m  + 107 m + 210)
                           -------------------------------------------------
                                            GAMMA(m + 9/2)

                                          4       3        2
                         8 GAMMA(m + 1) (m  + 26 m  + 251 m  + 1066 m + 1680)
                         ----------------------------------------------------
                                            GAMMA(m + 9/2)

                                              4       3        2
                     9 GAMMA(m + 1) (m + 5) (m  + 30 m  + 335 m  + 1650 m + 3024)
                     ------------------------------------------------------------
                                           GAMMA(m + 11/2)

                                     5       4        3         2
                   10 GAMMA(m + 1) (m  + 40 m  + 635 m  + 5000 m  + 19524 m + 30240)
                   -----------------------------------------------------------------
                                            GAMMA(m + 11/2)

*)

# For k >= 2, it looks like the "dominating" term is

#   sqrt(pi) * k * Gamma(m + 1) / Gamma(m + A(k) / 2) * (m^B(k) + O(m^{B(k) - 1})),

# where

#   A(k) = 2 * floor((k + 1) / 2) + 1
#   B(k) = floor((k + 1) / 2).

# From NIST, Gamma(m + 1) / Gamma(m + A(k) / 2) ~ m^(1 - A(k) / 2), thus the
# dominating term is

#       sqrt(pi) * k m^(1 + B(k) - A(k) / 2) (1 + O(1 / m))
#       sqrt(pi) * k m^{1 / 2} (1 + O(1 / m)).

# This is all *conjectural*, but seems right!