CMPQ    BX, $0x2000000
	JAE	loop_preheader_avx2_huge

loop_avx2:
	VMOVDQU	Y0, 0(DI)
	VMOVDQU	Y0, 32(DI)
	VMOVDQU	Y0, 64(DI)
	VMOVDQU	Y0, 96(DI)
	SUBQ	$128, BX
	ADDQ	$128, DI
	CMPQ	BX, $128
	JAE	loop_avx2
	VMOVDQU  Y0, -32(DI)(BX*1)
	VMOVDQU  Y0, -64(DI)(BX*1)
	VMOVDQU  Y0, -96(DI)(BX*1)
	VMOVDQU  Y0, -128(DI)(BX*1)
	VZEROUPPER
	RET

loop_preheader_erms:
#ifndef hasAVX2

	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

	VGF2P8AFFINEINVQB $27, Y7, Y12, K3, Y0             // 62f39d2bcfc71b
	VGF2P8AFFINEINVQB $27, (BX), Y12, K3, Y0           // 62f39d2bcf031b
	VGF2P8AFFINEINVQB $27, -17(BP)(SI*1), Y12, K3, Y0  // 62f39d2bcf8435efffffff1b
	VGF2P8AFFINEINVQB $27, Y5, Y3, K3, Y0              // 62f3e52bcfc51b
	VGF2P8AFFINEINVQB $27, Y28, Y3, K3, Y0             // 6293e52bcfc41b
	VGF2P8AFFINEINVQB $27, Y7, Y3, K3, Y0              // 62f3e52bcfc71b

		y6 = (y8 - w3mw5*y6) >> 3
		y7 = (y8 - w3pw5*y7) >> 3

		// Stage 2.
		y8 = y0 + y1
		y0 -= y1
		y1 = w6*(y3+y2) + 4
		y2 = (y1 - w2pw6*y2) >> 3
		y3 = (y1 + w2mw6*y3) >> 3
		y1 = y4 + y6
		y4 -= y6
		y6 = y5 + y7
		y5 -= y7

		// Stage 3.
		y7 = y8 + y3
		y8 -= y3
		y3 = y0 + y2
		y0 -= y2
		y2 = (r2*(y4+y5) + 128) >> 8
		y4 = (r2*(y4-y5) + 128) >> 8

		// Stage 4.

	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
}

		println("BUG: got", f, "want -0.0")
	}

	test32(-0.0)
	test32(x0)
	test32(x1)
	test32(x2)
	test32(x3)

	if f := -y0; f != 0 || !math.Signbit(f) {
		println("BUG: got", f, "want -0.0")
	}

	test64(-0.0)
	test64(y0)
	test64(y1)
	test64(y2)
	test64(y3)

		return 1 + z*(-0.25+(r/s)) // |x| < 1.00
	}
	u := 0.5 * x
	return (1+u)*(1-u) + z*(r/s) // 1.0 < |x| < 2.0
}

// Y0 returns the order-zero Bessel function of the second kind.
//
// Special cases are:
//
//	Y0(+Inf) = 0
//	Y0(0) = -Inf
//	Y0(x < 0) = NaN
//	Y0(NaN) = NaN
func Y0(x float64) float64 {
	const (
		TwoM27 = 1.0 / (1 << 27)             // 2**-27 0x3e40000000000000