#! Given two GF(p^5) elements on stack, this routine computes modular
#! addition over extension field GF(p^5) s.t. p = 2^64 - 2^32 + 1
#!
#! Expected stack state :
#!
#! [a0, a1, a2, a3, a4, b0, b1, b2, b3, b4, ...]
#!
#! After application of routine stack :
#!
#! [c0, c1, c2, c3, c4, ...] s.t. c = a + b
#!
#! See section 3.2 of https://eprint.iacr.org/2022/274.pdf
#!
#! For reference implementation in high level language, see
#! https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L607-L616
export.add
    repeat.5
        movup.5
        add
        movdn.4
    end
end

#! Given two GF(p^5) elements on stack, this routine subtracts second
#! element from first one, over extension field GF(p^5) s.t. p = 2^64 - 2^32 + 1
#!
#! Expected stack state :
#!
#! [a0, a1, a2, a3, a4, b0, b1, b2, b3, b4, ...]
#!
#! After application of routine stack :
#!
#! [c0, c1, c2, c3, c4, ...] s.t. c = a - b
#!
#! See section 3.2 of https://eprint.iacr.org/2022/274.pdf
#!
#! For reference implementation in high level language, see
#! https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L629-L638
export.sub
    repeat.5
        movup.5
        sub
        movdn.4
    end
end

#! Given two GF(p^5) elements on stack, this routine computes modular
#! multiplication ( including reduction by irreducible polynomial )
#! over extension field GF(p^5) s.t. p = 2^64 - 2^32 + 1
#!
#! Expected stack state :
#!
#! [a0, a1, a2, a3, a4, b0, b1, b2, b3, b4, ...]
#!
#! After application of routine stack :
#!
#! [c0, c1, c2, c3, c4, ...] s.t. c = a * b
#!
#! See section 3.2 of https://eprint.iacr.org/2022/274.pdf
#!
#! For reference implementation in high level language, see
#! https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L676-L689
export.mul
    # compute {c0, c1, c2, c3, c4} - five coefficients of resulting
    # degree-4 polynomial

    # compute c4
    dup.9
    dup.1
    mul

    dup.9
    dup.3
    mul

    add

    dup.8
    dup.4
    mul

    add

    dup.7
    dup.5
    mul

    add

    dup.6
    dup.6
    mul

    add

    # compute c3
    dup.9
    dup.2
    mul

    dup.9
    dup.4
    mul

    add

    dup.8
    dup.5
    mul

    add

    dup.7
    dup.6
    mul

    add

    dup.11
    dup.7
    mul
    mul.3

    add

    # compute c2
    dup.9
    dup.3
    mul

    dup.9
    dup.5
    mul

    add

    dup.8
    dup.6
    mul

    add

    dup.12
    dup.7
    mul
    mul.3

    add

    dup.11
    dup.8
    mul
    mul.3

    add

    # compute c1
    dup.9
    dup.4
    mul

    dup.9
    dup.6
    mul

    add

    dup.13
    dup.7
    mul
    mul.3

    add

    dup.12
    dup.8
    mul
    mul.3

    add

    dup.11
    dup.9
    mul
    mul.3

    add

    # compute c0
    movup.9
    movup.5
    mul

    movup.12
    movup.6
    mul
    mul.3

    add

    movup.10
    movup.6
    mul
    mul.3

    add

    movup.8
    movup.6
    mul
    mul.3

    add

    movup.6
    movup.6
    mul
    mul.3

    add
end

#! Given one GF(p^5) element on stack, this routine computes modular
#! squaring ( including reduction by irreducible polynomial )
#! over extension field GF(p^5) s.t. p = 2^64 - 2^32 + 1
#!
#! This routine has same effect as calling mul(a, a) | a ∈ GF(p^5)
#!
#! Expected stack state :
#!
#! [a0, a1, a2, a3, a4, ...]
#!
#! After application of routine stack :
#!
#! [b0, b1, b2, b3, b4, ...] s.t. b = a * a
#!
#! See section 3.2 of https://eprint.iacr.org/2022/274.pdf
#!
#! For reference implementation in high level language, see
#! https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L709-L715
export.square
    # compute {b0, b1, b2, b3, b4} - five coefficients of resulting
    # degree-4 polynomial

    # compute b4
    dup.2
    dup.3
    mul

    dup.5
    dup.2
    mul
    mul.2

    add

    dup.4
    dup.3
    mul
    mul.2

    add

    # compute b3
    dup.4
    dup.2
    mul
    mul.2

    dup.4
    dup.4
    mul
    mul.2

    add

    dup.6
    dup.7
    mul
    mul.3

    add

    # compute b2
    dup.3
    dup.4
    mul

    dup.5
    dup.4
    mul
    mul.2

    add

    dup.7
    dup.7
    mul
    mul.6

    add

    # compute b1
    dup.4
    dup.4
    mul
    mul.2

    dup.7
    dup.8
    mul
    mul.3

    add

    dup.8
    dup.7
    mul
    mul.6

    add

    # compute b0
    dup.4
    movup.5
    mul

    movup.8
    movup.6
    mul
    mul.6

    add

    movup.6
    movup.6
    mul
    mul.6

    add
end

#! Given an element a ∈ GF(p^5), this routine applies Frobenius operator
#! once, raising the element to the power of p | p = 2^64 - 2^32 + 1.
#!
#! Expected stack state :
#!
#! [a0, a1, a2, a3, a4, ...]
#!
#! Final stack state :
#!
#! [b0, b1, b2, b3, b4, ...]
#!
#! See https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L723-L737
#! for reference implementation in high-level language.
proc.frobenius_once
    movup.4
    mul.1373043270956696022

    movup.4
    mul.211587555138949697

    movup.4
    mul.15820824984080659046

    movup.4
    mul.1041288259238279555

    movup.4
end

#! Given an element a ∈ GF(p^5), this routine applies Frobenius operator
#! twice, raising the element to the power of p^2 | p = 2^64 - 2^32 + 1.
#!
#! Expected stack state :
#!
#! [a0, a1, a2, a3, a4, ...]
#!
#! Final stack state :
#!
#! [b0, b1, b2, b3, b4, ...]
#!
#! See https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L739-L749
#! for reference implementation in high-level language.
proc.frobenius_twice
    movup.4
    mul.211587555138949697

    movup.4
    mul.1041288259238279555

    movup.4
    mul.1373043270956696022

    movup.4
    mul.15820824984080659046

    movup.4
end

#! Given one GF(p^5) element on stack, this routine computes multiplicative
#! inverse over extension field GF(p^5) s.t. p = 2^64 - 2^32 + 1
#!
#! Expected stack state :
#!
#! [a0, a1, a2, a3, a4, ...]
#!
#! After application of routine stack :
#!
#! [b0, b1, b2, b3, b4, ...] s.t. b = 1 / a
#!
#! See section 3.2 of https://eprint.iacr.org/2022/274.pdf
#!
#! For reference implementation in high level language, see
#! https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L751-L775
#!
#! Note, this routine will not panic even when operand `a` is zero.
export.inv
    repeat.5
        dup.4
    end

    exec.frobenius_once # = t0

    repeat.5
        dup.4
    end

    exec.frobenius_once # = t0.frobenius_once()
    exec.mul            # = t1

    repeat.5
        dup.4
    end

    exec.frobenius_twice # = t1.frobenius_twice()
    exec.mul             # = t2

    movup.5
    dup.1
    mul

    movup.6
    dup.6
    mul
    mul.3

    add

    movup.6
    dup.5
    mul
    mul.3

    add

    movup.6
    dup.4
    mul
    mul.3

    add

    movup.6
    dup.3
    mul
    mul.3

    add                    # = t3

    dup
    push.0
    eq
    add
    inv                    # = t4

    movup.5
    dup.1
    mul

    movup.5
    dup.2
    mul

    movup.5
    dup.3
    mul

    movup.5
    dup.4
    mul

    movup.5
    movup.5
    mul
end

#! Given two GF(p^5) elements ( say a, b ) on stack, this routine computes
#! modular division over extension field GF(p^5) s.t. p = 2^64 - 2^32 + 1
#!
#! Expected stack state :
#!
#! [a0, a1, a2, a3, a4, b0, b1, b2, b3, b4, ...]
#!
#! After application of routine stack :
#!
#! [c0, c1, c2, c3, c4, ...] s.t. c = a / b
#!
#! See section 3.2 of https://eprint.iacr.org/2022/274.pdf
#!
#! For reference implementation in high level language, see
#! https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L777-L781
export.div
    repeat.5
        movup.9
    end

    exec.inv
    exec.mul
end

#! Given an element v ∈ Z_q | q = 2^64 - 2^32 + 1, and n on stack, this routine
#! raises it to the power 2^n, by means of n successive squarings
#!
#! Expected stack stack
#!
#! [v, n, ...] | n >= 0
#!
#! After finishing execution stack
#!
#! [v', ...] s.t. v' = v ^ (2^n)
#!
#! See https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L461-L469
#! for reference implementation in higher level language
proc.base_msquare
    swap
    dup
    neq.0

    while.true
        sub.1
        swap
        dup
        mul

        swap
        dup
        neq.0
    end

    drop
end

#! Given an element v ∈ Z_q | q = 2^64 - 2^32 + 1, this routine attempts to compute
#! square root of v, if that number is a square.
#!
#! Expected stack state :
#!
#! [v, ...]
#!
#! After finishing execution stack looks like :
#!
#! [v', flg, ...]
#!
#! If flg = 1, it denotes v' is square root of v i.e. v' * v' = v ( mod q )
#! If flg = 0, then v' = 0, denoting v doesn't have a square root
#!
#! See https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L349-L446
#! for reference implementation in higher level language.
proc.base_sqrt
    dup # = x

    push.31
    swap
    exec.base_msquare # = u

    dup
    dup
    mul # = u^2

    movup.2
    dup
    eq.0
    add

    div # = v

    # j = 1
    # i = 32 - j = 31
    dup
    push.30
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.4614640910117430873
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.1753635133440165772
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 2
    # i = 32 - j = 30
    dup
    push.29
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.9123114210336311365
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.4614640910117430873
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 3
    # i = 32 - j = 29
    dup
    push.28
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.16116352524544190054
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.9123114210336311365
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 4
    # i = 32 - j = 28
    dup
    push.27
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.6414415596519834757
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.16116352524544190054
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 5
    # i = 32 - j = 27
    dup
    push.26
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.1213594585890690845
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.6414415596519834757
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 6
    # i = 32 - j = 26
    dup
    push.25
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.17096174751763063430
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.1213594585890690845
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 7
    # i = 32 - j = 25
    dup
    push.24
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.5456943929260765144
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.17096174751763063430
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 8
    # i = 32 - j = 24
    dup
    push.23
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.9713644485405565297
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.5456943929260765144
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 9
    # i = 32 - j = 23
    dup
    push.22
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.16905767614792059275
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.9713644485405565297
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 10
    # i = 32 - j = 22
    dup
    push.21
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.5416168637041100469
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.16905767614792059275
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 11
    # i = 32 - j = 21
    dup
    push.20
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.17654865857378133588
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.5416168637041100469
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 12
    # i = 32 - j = 20
    dup
    push.19
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.3511170319078647661
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.17654865857378133588
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 13
    # i = 32 - j = 19
    dup
    push.18
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.18146160046829613826
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.3511170319078647661
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 14
    # i = 32 - j = 18
    dup
    push.17
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.9306717745644682924
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.18146160046829613826
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 15
    # i = 32 - j = 17
    dup
    push.16
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.12380578893860276750
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.9306717745644682924
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 16
    # i = 32 - j = 16
    dup
    push.15
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.6115771955107415310
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.12380578893860276750
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 17
    # i = 32 - j = 15
    dup
    push.14
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.17776499369601055404
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.6115771955107415310
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 18
    # i = 32 - j = 14
    dup
    push.13
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.16207902636198568418
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.17776499369601055404
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 19
    # i = 32 - j = 13
    dup
    push.12
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.1532612707718625687
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.16207902636198568418
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 20
    # i = 32 - j = 12
    dup
    push.11
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.17492915097719143606
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.1532612707718625687
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 21
    # i = 32 - j = 11
    dup
    push.10
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.455906449640507599
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.17492915097719143606
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 22
    # i = 32 - j = 10
    dup
    push.9
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.11353340290879379826
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.455906449640507599
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 23
    # i = 32 - j = 9
    dup
    push.8
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.1803076106186727246
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.11353340290879379826
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 24
    # i = 32 - j = 8
    dup
    push.7
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.13797081185216407910
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.1803076106186727246
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 25
    # i = 32 - j = 7
    dup
    push.6
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.17870292113338400769
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.13797081185216407910
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 26
    # i = 32 - j = 6
    dup
    push.5
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.549755813888
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.17870292113338400769
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 27
    # i = 32 - j = 5
    dup
    push.4
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.70368744161280
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.549755813888
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 28
    # i = 32 - j = 4
    dup
    push.3
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.17293822564807737345
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.70368744161280
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 29
    # i = 32 - j = 3
    dup
    push.2
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.18446744069397807105
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.17293822564807737345
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 30
    # i = 32 - j = 2
    dup
    push.1
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.281474976710656
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.18446744069397807105
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap

    # j = 31
    # i = 32 - j = 1
    dup
    push.0
    swap
    exec.base_msquare # = w

    eq.18446744069414584320 # = cc

    dup.1
    mul.18446744069414584320
    movup.2
    swap
    dup.2
    cdrop # = v'

    dup.2
    mul.281474976710656
    movup.3
    swap
    movup.3
    cdrop # = u'

    swap # On stack [v, u, ...]

    dup
    eq.0
    swap
    eq.1
    or # = cc

    swap
    dup.1
    mul # On stack [u * cc, cc, ...]
end

#! Given an element v ∈ Z_q | q = 2^64 - 2^32 + 1, this routine computes
#! legendre symbol, by raising that element to the power (p-1) / 2
#!
#! Expected stack state :
#!
#! [v, ...]
#!
#! After finishing execution stack looks like
#!
#! [v', ...] s.t. v' = legendre symbol of v
#!
#! See https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L448-L459
#! for reference implementation in higher level language.
proc.base_legendre
    repeat.31
        dup
        mul
    end

    dup

    repeat.32
        dup
        mul
    end

    swap
    dup
    eq.0
    add

    div
end

#! Given an element v ∈ GF(p^5), this routine computes its legendre symbol,
#! which is an element ∈ GF(p) | p = 2^64 - 2^32 + 1
#!
#! At beginning stack looks like
#!
#! [a0, a1, a2, a3, a4, ...]
#!
#! At end stack looks like
#!
#! [b, ...] s.t. b = legendre symbol of a
#!
#! See https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L857-L877
#! for reference implementation in higher level language.
export.legendre
    repeat.5
        dup.4
    end

    exec.frobenius_once

    repeat.5
        dup.4
    end

    exec.frobenius_once
    exec.mul

    repeat.5
        dup.4
    end

    exec.frobenius_twice
    exec.mul

    movup.5
    mul

    movup.5
    movup.5
    mul
    mul.3

    add

    movup.4
    movup.4
    mul
    mul.3

    add

    movup.3
    movup.3
    mul
    mul.3

    add

    movup.2
    movup.2
    mul
    mul.3

    add

    exec.base_legendre
end

#! Given an element v ∈ GF(p^5), this routine attempts to compute square root of v,
#! if that number is a square.
#!
#! At beginning stack looks like
#!
#! [a0, a1, a2, a3, a4, ...]
#!
#! At end stack looks like
#!
#! [b0, b1, b2, b3, b4, flg, ...]
#!
#! If flg = 1, it denotes v' = {b0, b1, b2, b3, b4} is square root of v i.e. v' * v' = v ( mod GF(p^5) )
#! If flg = 0, then v' = {0, 0, 0, 0, 0}, denoting v doesn't have a square root
#!
#! See https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L879-L910
#! for reference implementation in higher level language.
export.sqrt
    repeat.5
        dup.4
    end

    repeat.31
        repeat.5
            dup.4
        end

        exec.mul
    end # = v

    repeat.5
        dup.4
    end

    repeat.32
        repeat.5
            dup.4
        end

        exec.mul
    end

    exec.div

    repeat.5
        dup.9
    end

    exec.mul # = d

    repeat.5
        dup.4
    end

    exec.frobenius_twice
    exec.mul
    exec.frobenius_once # = e

    repeat.5
        dup.4
    end

    exec.square # = f

    movup.10
    mul

    swap
    movup.13
    mul
    mul.3
    add

    swap
    movup.11
    mul
    mul.3
    add

    swap
    movup.9
    mul
    mul.3
    add

    swap
    movup.7
    mul
    mul.3
    add # = g

    exec.base_sqrt # On stack [s, c, e0, e1, e2, e3, e4, ...]

    repeat.5
        movup.6
    end

    exec.inv # = e'

    repeat.5
        movup.4
        dup.5
        mul
    end

    movup.5
    drop # On stack [e0, e1, e2, e3, e4, c, ...]
end

#! Given two elements a, b ∈ GF(p^5), this routine produces single field element r,
#! denoting whether a == b.
#!
#! Expected stack state
#!
#! [a0, a1, a2, a3, a4, b0, b1, b2, b3, b4, ...]
#!
#! Final stack state
#!
#! [r, ...]
#!
#! If a == b { r = 1 } Else { r = 0 }
#!
#! See https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L797-L806
#! for reference implementation.
export.eq
    push.1

    swap
    movup.6
    eq
    and

    swap
    movup.5
    eq
    and

    swap
    movup.4
    eq
    and

    swap
    movup.3
    eq
    and

    swap
    movup.2
    eq
    and
end

#! Given two elements a, b ∈ GF(p^5), this routine produces single field element r,
#! denoting whether a != b.
#!
#! Expected stack state
#!
#! [a0, a1, a2, a3, a4, b0, b1, b2, b3, b4, ...]
#!
#! Final stack state
#!
#! [r, ...]
#!
#! If a != b { r = 1 } Else { r = 0 }
#!
#! See https://github.com/pornin/ecgfp5/blob/ce059c6/python/ecGFp5.py#L813-L822
#! for reference implementation.
export.neq
    push.0

    swap
    movup.6
    neq
    or

    swap
    movup.5
    neq
    or

    swap
    movup.4
    neq
    or

    swap
    movup.3
    neq
    or

    swap
    movup.2
    neq
    or
end