VPGATHERDQ Y0, 664(X4*8), Y6            // c4e2fd9034e598020000
	VPGATHERDQ Y0, 8(X14*1), Y6             // c4a2fd90343508000000
	VPGATHERDQ Y0, -8(X14*1), Y6            // c4a2fd903435f8ffffff
	VPGATHERDQ Y0, 0(X14*1), Y6             // c4a2fd90343500000000
	VPGATHERDQ Y0, 664(X14*1), Y6           // c4a2fd90343598020000
	VPGATHERDQ Y0, 8(X14*8), Y6             // c4a2fd9034f508000000

	LDP	0*16(a_ptr), (x0, x1)
	LDP	1*16(a_ptr), (x2, x3)
	LDP	0*16(b_ptr), (y0, y1)
	LDP	1*16(b_ptr), (y2, y3)

	// y[0] * x
	MUL	y0, x0, acc0
	UMULH	y0, x0, acc1

	MUL	y0, x1, t0
	ADDS	t0, acc1
	UMULH	y0, x1, acc2

	MUL	y0, x2, t0
	ADCS	t0, acc2
	UMULH	y0, x2, acc3

	MUL	y0, x3, t0
	ADCS	t0, acc3
	UMULH	y0, x3, acc4
	ADC	$0, acc4
	// First reduction step
	MUL	acc0, hlp1, hlp0

	LXVD2X (R19)(P1ptr), Y0
	LXVD2X (R20)(P1ptr), Y1
	XXPERMDI Y0, Y0, $2, Y0
	XXPERMDI Y1, Y1, $2, Y1

	CALL p256MulInternal<>(SB)

	// Leave T0, T1 as is.
	XXPERMDI T0, T0, $2, TT0
	XXPERMDI T1, T1, $2, TT1
	STXVD2X TT0, (R19)(P3ptr)
	STXVD2X TT1, (R20)(P3ptr)

	// X-  ; Y=X ; MUL; T-    // Y3 = Y3²
	VOR  X0, X0, Y0
	VOR  X1, X1, Y1
	CALL p256MulInternal<>(SB)

	ADDQ R15, R11
	ADDQ DI, R11

	VPSRLQ $8, Y1, Y2
	VPSLLQ $(64-8), Y1, Y1
	VPOR   Y2, Y1, Y1

	VPXOR Y8, Y3, Y3
	VPXOR Y1, Y3, Y1

	VPADDQ Y1, Y0, Y0

	VPERM2F128 $0x0, Y0, Y0, Y4

	VPAND MASK_YMM_LO<>(SB), Y0, Y0

	VPERM2F128 $0x11, Y7, Y7, Y2
	VPSRLQ     $6, Y2, Y8

	MOVQ  R11, DI
	RORXQ $41, R8, R13
	RORXQ $18, R8, R14
	ADDQ  1*8+frame_YFER(SP), R10
	ORQ   BX, DI

// On supported platforms this is an intrinsic lowered by the compiler.
func Mul64(x, y uint64) (hi, lo uint64) {
	const mask32 = 1<<32 - 1
	x0 := x & mask32
	x1 := x >> 32
	y0 := y & mask32
	y1 := y >> 32
	w0 := x0 * y0
	t := x1*y0 + w0>>32
	w1 := t & mask32
	w2 := t >> 32
	w1 += x0 * y1
	hi = x1*y1 + w2 + w1>>32
	lo = x * y
	return
}

#define X0    V0
#define X1    V1
#define Y0    V2
#define Y1    V3

TEXT p256SqrInternal<>(SB), NOFRAME|NOSPLIT, $0
	VLR X0, Y0
	VLR X1, Y1
	BR  p256MulInternal<>(SB)

#undef X0
#undef X1
#undef Y0
#undef Y1

#define p256SubInternal(T1, T0, X1, X0, Y1, Y0) \
	VZERO   ZER                \
	VSCBIQ  Y0, X0, CAR1       \
	VSQ     Y0, X0, T0         \
	VSBCBIQ X1, Y1, CAR1, SEL1 \

	RORXL   $13, a, T1;                  \ // T1 = a >> 13					// S0B
	XORL    y1, y0;                      \ // y0 = (e>>25) ^ (e>>11)		// S1
	XORL    g, y2;                       \ // y2 = f^g						// CH
	;                                    \
	RORXL   $6, e, y1;                   \ // y1 = (e >> 6)				// S1
	XORL    y1, y0;                      \ // y0 = (e>>25) ^ (e>>11) ^ (e>>6)	// S1

	x0, x1, dx := r.Min.X, r.Max.X, 1
	y0, y1, dy := r.Min.Y, r.Max.Y, 1
	if r.Overlaps(r.Add(sp.Sub(r.Min))) {
		if sp.Y < r.Min.Y || sp.Y == r.Min.Y && sp.X < r.Min.X {
			x0, x1, dx = x1-1, x0-1, -1
			y0, y1, dy = y1-1, y0-1, -1
		}
	}

	sy := sp.Y + y0 - r.Min.Y
	my := mp.Y + y0 - r.Min.Y
	sx0 := sp.X + x0 - r.Min.X
	mx0 := mp.X + x0 - r.Min.X

	}
	x, y = curve.Add(xG, yG, x0, y0)
	if x.Cmp(xG) != 0 || y.Cmp(yG) != 0 {
		t.Errorf("x+∞ != x")
	}
	x, y = curve.Add(x0, y0, xG, yG)
	if x.Cmp(xG) != 0 || y.Cmp(yG) != 0 {
		t.Errorf("∞+x != x")
	}

	if curve.IsOnCurve(x0, y0) {
		t.Errorf("IsOnCurve(∞) == true")
	}

	if xx, yy := Unmarshal(curve, Marshal(curve, x0, y0)); xx != nil || yy != nil {

	}

	// compute yd (or the negative value if underflow occurs)
	yd := z[2*n+n2 : 3*n]
	if subVV(yd, y0, y1) != 0 { // y0-y1
		s = -s
		subVV(yd, y1, y0) // y1-y0
	}

	// p = (x1-x0)*(y0-y1) == x1*y0 - x1*y1 - x0*y0 + x0*y1 for s > 0
	// p = (x0-x1)*(y0-y1) == x0*y0 - x0*y1 - x1*y0 + x1*y1 for s < 0
	p := z[n*3:]
	karatsuba(p, xd, yd)

	// save original z2:z0