include("twihes.jl")

############################
#
# Parameters
#
############################

#Prime and field
const la = big(2)
const lb = big(3)
const ea = 305
const eb = 192
const f = big(1)
const p = la^ea*lb^eb*f-1
const F = Fp2(p)

#cuberoot of One
const w = F(big(1319470205536631335981810349644143385091780115593611158375203778232819918005279075081612765167581266740524628378608982325666905211062864268703698577903039608673459647224627806555078655), big(1286451788956606916776387017696978305185994757974417726547130267369253169617279161336419642919745234101683127375727357361442589772170347650505223652853823514074141695255872287918030044))
if !isone(w^3)
    error("w is not a cuberoot...")
end
#Starting curve: a = 1 for speed, solve for d to get j(E) = some supersingular j-invariant. See sage-script
#Starting j-invariant = 831837561004103883868671088207491838344992437494159978595307997269883652199194625608995144292694837047770068423678474079025058740865040338182857297729622238883560881131572628427170113 + 771813078065723508849643803754084945896394345666593331261515932318350081581366794480345255301600363646900179413481795555101011900773965712484190126442309614812527433851551995117691457*i
#See sage script for construction of d0 from j-Invariant and a0
const a0 = one(F)
const d0 = F(big(526466053434329606348792315213226961488878919273384911349255872484662011698700897679785221256842528945834974679658632918959678155926738589102758106233735036087191018029528283224394826), big(322622151547826919747264788491387568702182872571850474819101171989083545313493484238990598522402415188382526845235802789803276920938914974974558706304117646595374335297821576916453987))


#Generates torsion basises from random points
function genTorsionBasis(alicePX, alicePY, aliceQX, aliceQY, bobPX, bobPY, bobQX, bobQY)
    Z = one(F)
    aliceP = EC(alicePX, alicePY, Z)
    aliceQ = EC(aliceQX, aliceQY, Z)

    if !isValid(aliceP, a0, d0)
        error("Alice' P is not a valid point!")
    end

    if !isValid(aliceQ, a0, d0)
        error("Alice' Q is not a valid point!")
    end

    Pa = normalized(mul(f*lb^eb, aliceP, a0))
    Qa = normalized(mul(f*lb^eb, aliceQ, a0))

    Patest = doublePow(Pa, a0, ea-1)
    Qatest = doublePow(Qa, a0, ea-1)

    if !isIdentity(double(Patest, a0))
        error("Pa : Something does not look right with the order of the curve...")
    end

    if !isIdentity(double(Qatest, a0))
        error("Qa : Something does not look right with the order of the curve...")
    end

    if isIdentity(Patest)
        error("Alice' point Pa does not have the right order!")
    end
    if isIdentity(Qatest)
        error("Alice' point Qa does not have the right order!")
    end

    if Patest == Qatest
        error("Alice' points are not independent!")
    end

    bobP = EC(bobPX, bobPY, Z)
    bobQ = EC(bobQX, bobQY, Z)

    if !isValid(bobP, a0, d0)
        error("Bob's P is not a valid point!")
    end

    if !isValid(bobQ, a0, d0)
        error("Bob's Q is not a valid point!")
    end

    Pb = normalized(mul(f*la^ea, bobP, a0))
    Qb = normalized(mul(f*la^ea, bobQ, a0))

    Pbtest = triplePow(Pb, a0, d0, eb-1)
    Qbtest = triplePow(Qb, a0, d0, eb-1)

    if !isIdentity(triple(Pbtest, a0, d0))
        error("Pb : Something does not look right with the order of the curve...")
    end

    if !isIdentity(triple(Qbtest, a0, d0))
        error("Qb : Something does not look right with the order of the curve...")
    end

    if isIdentity(Pbtest)
        error("Bob's point Pb does not have the right order!")
    end

    if isIdentity(Qbtest)
        error("Bob's point Qb does not have the right order!")
    end

    if Pbtest == Qbtest
        error("Bob's points are not independent!")
    end

    return Pa, Qa, Pb, Qb
end

#See sage script for generating random points
PX1 = F(big(2492006491871983676452451839020564198147973496177388329872836720660635968389647718878527574597574695774796088332985623054551169799526739171989607694478941255607446159858803583994964802), big(1106461312022518869103531097009806634520011809887491682622823938979378314609312783672420779078369845610390072881748838641128935933666342275393037280723333398214385754330629136835945615))
PY1 = F(big(568934704890564342552989237837373499705027664549344609870746088351623756422430685800281891502113550175496052961717514133893321580563982611360102558084198226260195186543398872742166956), big(1719319066282017378159238737975467015356375587815788392734006098067712564812437626919347516402785703024563499386845184938500959018676538501237827797206911141692035339948528241047463062))
QX1 = F(big(2517628532562506074245481407647129152141059415595858026166161117305445674436190840142297568565926716458956505364881489417252719704505645622987495305020004927497493134683221021490239959), big(1433263917761027493460060659243243110082141390402415821474572364169645458452572424878937372751197553782045370168747546443025397983411476342832168979291996203472762061817060320239732514))
QY1 = F(big(1268486271669556135524139358247425138875339494127874128975087540722306664933940157281323943300180893898834305285190175720304656654248228384119897648002761039198456950277535984384435960), big(870717938688314703034698794917729951973369728070119897658763707633485814633592046219916287841159157678092192037611181964671011367024700114197655484075793009679462073178107401411271146))
PX2 = F(big(1087251779026402284780690514599816326569275126413001776550994011768620358684555955678133912577235126758668170004393171497761389593621955822052207714093613781329468475697631867387474483), big(2469211885106793446942611087302463642283587716018469215650568305955990628425421892841320028890624561652810865126897501340849958642894467779823256283556264901094677686850284342042454424))
PY2 = F(big(2042684845558020702068009872943604896265659382991530782022341450605024204133331153581602581183717652980627685160557204787602548705930777904538999384048875962096503211454749544788971900), big(1952407962161458071632111111545938309344853088503297399264610051714768899866468987194505212765542390314727724800985868656215669828095236952855844761115022653991396268869521633770184697))
QX2 = F(big(1948557883143239148216542337109344092480664013194035763418796831802277689220591717171651875237624876590110990487339142795957213635061021311223506194896335710514651491187248696544839890), big(1561987540323249359134988383755115761443665949729569112934702895896188229271478800280961840021334352088226396987001996010282799910199376720872821614759744240212908203945883767928868921))
QY2 = F(big(2434803430031357671431182080606412625586548768754456399219302240585675185680813127143693009941726516214200502993859529022071607523524471311775260681161873786254554974167028608047308730), big(1246986093659258949315414845812705415673341880493879348339825574550001003997591027583897079561051808540633480867564839740453013245250501290325921658949100205817371785062094851400235056))
const Pa, Qa, Pb, Qb = genTorsionBasis(PX1, PY1, QX1, QY1, PX2, PY2, QX2, QY2)

# Compute the optimal strategy of size n, given weights p and q
function computeStrat(n, p, q)
    S = Dict{Int64, Array{Int64, 1}}()
    C = Dict{Int64, Int64}()
    S[1] = []
    C[1] = 0
    for i in 2:n
        b = argmin([C[i-b] + C[b] + b*p + (i-b)*q for b in 1:i-1])
        S[i] = append!(append!([b], S[i-b]), S[b])
        C[i] = C[i-b] + C[b] + b*p + (i-b)*q
    end
    return S[n]
end

#The weights are estimated using @benchmark
const strat_A = Vector(computeStrat(ea, 23358, 32132))

const strat_B = Vector(computeStrat(eb, 45620, 34278))
