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

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

#Don't make a calc 3 mistake!
with(plots):

Help:=proc():
    print(` NumerBifOrbF(f, x, x0, K1, K2, d1, n, eps) `):
    print(` FirstBif(f, x) , PlotMS(reso, K, distfrom0:=5, faraway:=20) `):
end:

##PROBLEM 2##
 
#NumerBifOrbF(f,x,x0,K1,K2,d1,n,eps): inputs the same as
#OrbF(f,x,g,x0,K1,K2,d1), (note, f(x) has to be a function that vanishes at 0 and 1) as well as a positive integer n, and a small
#positive number eps, and outputs a list, let's call it B, of
#length n, such that B[i] is a small interval [ai,bi], with bi-ai < eps, and such that when g=ai the generalized logistic map has period 2^(i-1) and when g=bi it has period 2^i.
#Make sure K2 > 2^n
NumerBifOrbF:=proc(f,x,x0,K1,K2,d1,n,eps)
local A, A2, B, g, L, supernops, orbsize, first,
        FixFirstInterval, ok, j, midpt, hi, lo:
    supernops:=proc(L)
        #print(L):
        if L[-1] in {Float(infinity), Float(-infinity)} then
            return infinity:
        else
            return nops({op(L)}):
        fi:
    end:
    #Look for first interval in A1 containing at least one zero
    #Make one of those zeroes not a zero
    FixFirstInterval:=proc(A1)
    local bot, top, boti, topi, i, nex, orbsiz, L1,
            flog, glog:
        bot:=A1[1]:
        boti:=1:
        topi:=true:
        for i from 2 to nops(A1) do
            if i = boti+1 and A1[i] <> 0 then
                bot:=A1[i]:
                boti:=i:
            elif A1[i] <> 0 then
                top:=A1[i]:
                topi:=i:
                break:
            fi:
        od:
        if topi = true then
            #No zeroes!
            return A1:
        fi:
        #Bisect until replace one of the zeroes
        while true do
            nex:=(bot + top)/2:
            print(nex):
            L1:=OrbF(f,x,nex,x0,K1,K2,d1):
            orbsiz:=supernops(L1):
            if orbsiz = infinity then
                glog:=infinity:
                flog:=infinity:
            else
                glog:=log[2](orbsiz):
                flog:=floor(glog):
            fi:
            #print(glog, flog):
            if glog = flog and flog <= n and A1[flog+1] = 0 then
                return [op(1..flog, A1), nex,
                    op(flog+2..nops(A1), A1)]:
            elif flog <= n and A1[flog+1] <> 0 then
                bot:=nex:
            else
                top:=nex:
            fi:
            #print(bot, top):
        od:
    end:
    #A[i] has orbit size 2^(i-1) at the end
    #Except for A[n+2]
    A:=[0$(n+2)]:
    g:=1.1:
    first:=true:
    while true do
        L:=OrbF(f,x,g,x0,K1,K2,d1):
        orbsize:=supernops(L):
        ok:=evalb(type(log[2](orbsize), integer)):
        if ok and orbsize <= 2^n then
            A[log[2](orbsize) + 1]:=g:
            if A[1] <> 0 then
                first:=false:
            fi:
        elif orbsize > 2^n and not first then
            break:
        fi:
        if first then
            g:=g/2:
        else
            g:=g*2:
        fi:
    od:
    #This number is too big
    A[n+2]:=g:
    
    print(A):
    
    #There might still be zeroes in A, so get rid of them
    while true do
        A2:=FixFirstInterval(A):
        if A2 = A then
            break:
        fi:
        A:=A2:
        print(A):
    od:
    
    #Time to build B from A
    B:=[seq([A[j], A[j+1]], j=1..n)]:
    
    #Bisect the B intervals
    for j from 1 to nops(B) do
        hi:=B[j][2]:
        lo:=B[j][1]:
        while B[j][2] - B[j][1] >= eps do
            midpt:=(lo + hi) / 2:
            L:=OrbF(f,x,midpt,x0,K1,K2,d1):
            orbsize:=supernops(L):
            if orbsize > 2^(j-1) then
                hi:=midpt:
                if orbsize = 2^j then
                    B[j][2]:=hi:
                fi:
            else
                lo:=midpt:
                B[j][1]:=lo:
            fi:
        od:
    od:
    
    return B:
end:

#NumerBifOrbF(x^2*(1-x),x,.5,10000,1000,10,4,.02);
#runs forever, since x^2*(1-x) has a stable equilibrium at x=0
#until about 7.10581467145, at which point the process blows up
#to infinity

#NumerBifOrbF(x*(1-x)^2,x,.5,10000,1000,10,4,.02);
#returns [[2.595312500, 2.612500000], [4.795312500, 4.812500000], [5.160546875, 5.175585938], [5.264645386, 5.276982117]]

#NumerBifOrbF(x^2*(1-x)^2,x,.5,10000,1000,10,4,.02);
#runs forever, since x^2*(1-x)^2 has a stable equilibrium at x=0
#until about 19.31370850, at which point the process blows up
#to infinity

##PROBLEM 3##

#FirstBif(f,x): inputs an expression f in x (that vanishes at
#x=0 and x=1), a varibale x, and a positive integer g and finds
#the unique g1 such that for g < g1 the map x->g*f(x) has
#(eventual) period one, and when g > g1 it has eventual period 2.
FirstBif:=proc(f, x) local solvef, solvef2, g, i, j:
    solvef:=[op({solve(x=g*f, x)} minus {0})]:
    solvef2:=[op({solve(subs(x=g*f, g*f)=x, x)}
        minus {op(solvef), 0})]:
    print(solvef, solvef2):
    #return solve({solvef2=solvef}, g):
    return min(seq(seq(solve(solvef[i]=solvef2[j], g),
        j=1..nops(solvef2)), i=1..nops(solvef))):
end:

#FirstBif(x^2*(1-x),x); is 16/3
#FirstBif(x*(1-x)^2,x); is 4
#FirstBif(x^2*(1-x)^2,x); is infinity
#FirstBif(x^3*(1-x)^3,x); is 16807/432

##PROBLEM 4##

#PlotMS(reso,K): inputs a small resolution reso, and a positive
#integer K, and includes those c in the complex plane
#(approximated with a grid of resolution reso) such that iterating
#
#z->z^2+c
#
#(starting at z=0), K times does NOT go far away.
PlotMS:=proc(reso, K, distfrom0:=5, faraway:=5)
local x, y, z, c, i, M:
    M:={}:
    for x from -distfrom0 by reso to distfrom0 do
        for y from -distfrom0 by reso to distfrom0 do
            z:=0:
            c:=x+I*y:
            for i from 1 to K do
                z:=z^2+c:
            od:
            if evalf(abs(z)) < faraway then
                M:=M union {[x, y]}:
            fi:
        od:
    od:
    return pointplot(M):
end:

#See MS.jpg for plot