# curve448.py

```import core
from core import *

##################
# Isogeny switch #
##################

# Traditionally, Ed448 doesn't use X448's birationally equivalent curve,
# but an isogeny thereof.  This introduces a slight complication when
# trying to lean on existing Ed448 code.
#
# Here I decided to describe both alternatives:
# - If isogeny is True, we use the standard Ed448 code.
# - If isogeny is False, we use the birationally equivalent Edwards curve.
#
# Both methods generate the exact same results.
isogeny = True

####################
# field parameters #
####################
GF.p = 2**448 - 2**224 - 1

def is_negative(self):
"""True iff self is odd"""
return (self.val % GF.p) % 2 == 1

GF.is_negative = is_negative

#########################
# Square root functions #
#########################
def sqrt(n):
"""Principal square root of n

If n is not a square, the behaviour is undefined.

To know how it works, refer to https//elligator.org/sqrt
"""
square = n**((GF.p+1) // 4)
square = cmove(square, -square, square.is_negative())
return square

def inv_sqrt(x):
"""Inverse square root of x

Returns (0         , True ) if x is zero.
Returns (sqrt( 1/x), True ) if x is non-zero square.
Returns (sqrt(-1/x), False) if x is not a square.

The return value is *not* guaranteed to be non-negative.
"""
isr       = x**((GF.p - 3) // 4)
legendre  = x * isr**2         # x**((p-1)/2) == -1, 0, or 1
is_square = legendre != GF(-1) # use constant time comparison
return isr, is_square

core.sqrt     = sqrt
core.inv_sqrt = inv_sqrt

################################
# Birational map (and isogeny) #
################################
def mt_to_edwards(u):
y = (u + GF(1)) / (u - GF(1))
x = sqrt((y**2 - GF(1)) / (Ed.d * y**2 - GF(1)))
return (x, y, GF(1))

def isogeny_to_ed(point):
x, y, z = point
x2 = x**2
y2 = y**2
du = z**2*GF(2) - x2 - y2
v  = y2 - x2 # dv, actually. We use v to save space
d  = du * v
u  = v * x * y * GF(2)*sqrt(Ed.d)
v  = (y2 + x2) * du
return (u, v, d)

def edwards_to_mt(point):
x, y, z = point  # in projective coordinates
return (y + z) / (y - z)

####################
# Curve parameters #
####################

# Montgomery constants (We already assume B = 1)
Mt.A = GF(156326)

# Edwards constants
Ed.a = GF(1) # 1 -> not twisted
if isogeny: Ed.d = GF(-39081)
else      : Ed.d = GF(39082) / GF(39081)

# curve order and cofactor
Mt.cofactor = 4

# Standard base point, that generates the prime order sub-group
Mt.base = GF(5)
if isogeny:
# From RFC 8032
Ed.base = (GF(224580040295924300187604334099896036246789641632564134246125461686950415467406032909029192869357953282578032075146446173674602635247710),
GF(298819210078481492676017930443930673437544040154080242095928241372331506189835876003536878655418784733982303233503462500531545062832660),
GF(1))
else:
# Base point of the (non-standard) birational curve
Ed.base = mt_to_edwards(Mt.base)

# Low order point (of order 4), used to add the cofactor component
# There are 2 such points: (1, 0) and (-1, 0)
# We chose (1, 0) somewhat arbitrarily
Ed.lop = (GF(1), GF(0), GF(1))

# "Dirty" Base point, that generates the whole curve.
# mt_base_c = mt_base + (lop * co_clear)
co_clear  = Mt.order % Mt.cofactor # 3
lop_c     = Ed.scalarmult(Ed.lop, co_clear)
ed_base   = isogeny_to_ed(Ed.base) if isogeny else Ed.base
Mt.base_c = edwards_to_mt(Ed.base_c)

"""Adding a low order point, fast

# Equivalent to the following, except constant time
x, y, z = point
if i == 0: return  x,  y, z
if i == 1: return  y, -x, z
if i == 2: return -x, -y, z
if i == 3: return -y,  x, z
"""
x, y, z = point
l    = (i//1) % 2 == 1
h    = (i//2) % 2 == 1
x, y = cswap(x, y, l)
x    = cmove(x, -x, h)
y    = cmove(y, -y, l != h) # use XOR instead of !=
return x, y, z

# Replacing the generic Ed.co_scalarmult()
# with a custom method that uses the fast add_lop()
# This saves a full point addition.
def co_scalarmult(scalar, c):
main_point  = Ed.scalarmult(Ed.base, clamp(scalar))
if isogeny:
main_point = isogeny_to_ed(main_point)
low_order_p = Ed.scalarmult(Ed.lop, c)
montgomery1 = edwards_to_mt(Ed.add(main_point, low_order_p)) # slow
montgomery2 = edwards_to_mt(add_lop(main_point, c))          # fast
if montgomery1 != montgomery2:
raise ValueError('Incoherent low order point selection')
return montgomery2

Ed.co_scalarmult = co_scalarmult

########################
# Elligator parameters #
########################
core.Z       = GF(-1)
core.ufactor = -core.Z            # ufactor ==  1
core.vfactor = sqrt(core.ufactor) # vfactor == -1
```