zInvSq := new(p256Element)

	for j := 0; j < 32; j++ {
		p.Add(p, base)

		// Convert p to affine coordinates.
		p256Inverse(zInv, &p.z)
		p256Sqr(zInvSq, zInv, 1)
		p256Mul(zInv, zInv, zInvSq)

		p256Mul(&p.x, &p.x, zInvSq)
		p256Mul(&p.y, &p.y, zInv)
		p.z = p256One

		if p256Equal(&table[j].x, &p.x) != 1 || p256Equal(&table[j].y, &p.y) != 1 {
			t.Fatalf("incorrect table entry at index %d", j)
		}
	}

	var t1 = new(p256Element)

	p256Sqr(z, in, 1)
	p256Mul(z, in, z)
	p256Sqr(z, z, 1)
	p256Mul(z, in, z)
	p256Sqr(t0, z, 3)
	p256Mul(t0, z, t0)
	p256Sqr(t1, t0, 6)
	p256Mul(t0, t0, t1)
	p256Sqr(t0, t0, 3)
	p256Mul(z, z, t0)
	p256Sqr(t0, z, 1)
	p256Mul(t0, in, t0)
	p256Sqr(t1, t0, 16)
	p256Mul(t0, t0, t1)
	p256Sqr(t0, t0, 15)
	p256Mul(z, z, t0)
	p256Sqr(t0, t0, 17)
	p256Mul(t0, in, t0)

DATA p256mul<>+0x38(SB)/8, $0x1c1d1e1f0c0d0e0f // SEL d0  0 d1 d0
DATA p256mul<>+0x40(SB)/8, $0x040506071c1d1e1f // SEL  0 d1 d0 d1
DATA p256mul<>+0x48(SB)/8, $0x0c0d0e0f1c1d1e1f // SEL  0 d1 d0 d1
DATA p256mul<>+0x50(SB)/8, $0x0405060704050607 // SEL  0  0 d1 d0
DATA p256mul<>+0x58(SB)/8, $0x1c1d1e1f0c0d0e0f // SEL  0  0 d1 d0
DATA p256mul<>+0x60(SB)/8, $0x0c0d0e0f1c1d1e1f // SEL d1 d0 d1 d0

DATA p256mul<>+0x00(SB)/8, $0xffffffff00000001 // P256
DATA p256mul<>+0x08(SB)/8, $0x0000000000000000 // P256
DATA p256mul<>+0x10(SB)/8, $0x00000000ffffffff // P256
DATA p256mul<>+0x18(SB)/8, $0xffffffffffffffff // P256
DATA p256mul<>+0x20(SB)/8, $0x1c1d1e1f00000000 // SEL d0  0  0 d0
DATA p256mul<>+0x28(SB)/8, $0x000000001c1d1e1f // SEL d0  0  0 d0
DATA p256mul<>+0x30(SB)/8, $0x0001020304050607 // SEL d0  0 d1 d0

	return e
}

// Sub sets e = t1 - t2, and returns e.
func (e *P256Element) Sub(t1, t2 *P256Element) *P256Element {
	p256Sub(&e.x, &t1.x, &t2.x)
	return e
}

// Mul sets e = t1 * t2, and returns e.
func (e *P256Element) Mul(t1, t2 *P256Element) *P256Element {
	p256Mul(&e.x, &t1.x, &t2.x)
	return e
}

// Square sets e = t * t, and returns e.
func (e *P256Element) Square(t *P256Element) *P256Element {

	x2 := ((x1 & arg3) | ((^x1) & arg2))
	*out1 = x2
}

// p256Mul multiplies two field elements in the Montgomery domain.
//
// Preconditions:
//
//	0 ≤ eval arg1 < m
//	0 ≤ eval arg2 < m
//
// Postconditions:
//
//	eval (from_montgomery out1) mod m = (eval (from_montgomery arg1) * eval (from_montgomery arg2)) mod m
//	0 ≤ eval out1 < m

	MOVD	y1, x1
	MOVD	y2, x2
	MOVD	y3, x3
	CBNZ	b_ptr, sqrLoop

	STP	(y0, y1), 0*16(res_ptr)
	STP	(y2, y3), 1*16(res_ptr)
	RET
/* ---------------------------------------*/
// func p256Mul(res, in1, in2 *p256Element)
TEXT ·p256Mul(SB),NOSPLIT,$0
	MOVD	res+0(FP), res_ptr
	MOVD	in1+8(FP), a_ptr
	MOVD	in2+16(FP), b_ptr

	MOVD	p256const0<>(SB), const0
	MOVD	p256const1<>(SB), const1

	LDP	0*16(a_ptr), (x0, x1)

	MOVQ acc2, (8*2)(res_ptr)
	MOVQ acc3, (8*3)(res_ptr)
	MOVQ res_ptr, x_ptr
	DECQ BX
	JNE  sqrLoop

	RET
/* ---------------------------------------*/
// func p256Mul(res, in1, in2 *p256Element)
TEXT ·p256Mul(SB),NOSPLIT,$0
	MOVQ res+0(FP), res_ptr
	MOVQ in1+8(FP), x_ptr
	MOVQ in2+16(FP), y_ptr
	// x * y[0]
	MOVQ (8*0)(y_ptr), t0

	MOVQ (8*0)(x_ptr), AX
	MULQ t0
	MOVQ AX, acc0

func (e *P224Element) Sub(t1, t2 *P224Element) *P224Element {
	p224Sub(&e.x, &t1.x, &t2.x)
	return e
}

// Mul sets e = t1 * t2, and returns e.
func (e *P224Element) Mul(t1, t2 *P224Element) *P224Element {
	p224Mul(&e.x, &t1.x, &t2.x)
	return e
}

// Square sets e = t * t, and returns e.
func (e *P224Element) Square(t *P224Element) *P224Element {
	p224Square(&e.x, &t.x)
	return e
}

func (e *P521Element) Sub(t1, t2 *P521Element) *P521Element {
	p521Sub(&e.x, &t1.x, &t2.x)
	return e
}

// Mul sets e = t1 * t2, and returns e.
func (e *P521Element) Mul(t1, t2 *P521Element) *P521Element {
	p521Mul(&e.x, &t1.x, &t2.x)
	return e
}

// Square sets e = t * t, and returns e.
func (e *P521Element) Square(t *P521Element) *P521Element {
	p521Square(&e.x, &t.x)
	return e
}