// Copyright ©2016 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

//+build !noasm,!appengine

#include "textflag.h"

#define MOVDDUP_XPTR__X3    LONG $0x1E120FF2 // MOVDDUP (SI), X3
#define MOVDDUP_XPTR_INCX__X5    LONG $0x120F42F2; WORD $0x062C // MOVDDUP (SI)(R8*1), X5
#define MOVDDUP_XPTR_INCX_2__X7    LONG $0x120F42F2; WORD $0x463C // MOVDDUP (SI)(R8*2), X7
#define MOVDDUP_XPTR_INCx3X__X9    LONG $0x120F46F2; WORD $0x0E0C // MOVDDUP (SI)(R9*1), X9

#define MOVDDUP_8_XPTR__X2    LONG $0x56120FF2; BYTE $0x08 // MOVDDUP 8(SI), X2
#define MOVDDUP_8_XPTR_INCX__X4    LONG $0x120F42F2; WORD $0x0664; BYTE $0x08 // MOVDDUP 8(SI)(R8*1), X4
#define MOVDDUP_8_XPTR_INCX_2__X6    LONG $0x120F42F2; WORD $0x4674; BYTE $0x08 // MOVDDUP 8(SI)(R8*2), X6
#define MOVDDUP_8_XPTR_INCx3X__X8    LONG $0x120F46F2; WORD $0x0E44; BYTE $0x08 // MOVDDUP 8(SI)(R9*1), X8

#define ADDSUBPD_X2_X3    LONG $0xDAD00F66 // ADDSUBPD X2, X3
#define ADDSUBPD_X4_X5    LONG $0xECD00F66 // ADDSUBPD X4, X5
#define ADDSUBPD_X6_X7    LONG $0xFED00F66 // ADDSUBPD X6, X7
#define ADDSUBPD_X8_X9    LONG $0xD00F4566; BYTE $0xC8 // ADDSUBPD X8, X9

#define X_PTR SI
#define Y_PTR DI
#define LEN CX
#define TAIL BX
#define SUM X0
#define P_SUM X1
#define INC_X R8
#define INCx3_X R9
#define INC_Y R10
#define INCx3_Y R11
#define NEG1 X15
#define P_NEG1 X14

// func DotcInc(x, y []complex128, n, incX, incY, ix, iy uintptr) (sum complex128)
TEXT ·DotcInc(SB), NOSPLIT, $0
	MOVQ   x_base+0(FP), X_PTR       // X_PTR = &x
	MOVQ   y_base+24(FP), Y_PTR      // Y_PTR = &y
	MOVQ   n+48(FP), LEN             // LEN = n
	PXOR   SUM, SUM                  // SUM = 0
	CMPQ   LEN, $0                   // if LEN == 0 { return }
	JE     dot_end
	PXOR   P_SUM, P_SUM              // P_SUM = 0
	MOVQ   ix+72(FP), INC_X          // INC_X = ix * sizeof(complex128)
	SHLQ   $4, INC_X
	MOVQ   iy+80(FP), INC_Y          // INC_Y = iy * sizeof(complex128)
	SHLQ   $4, INC_Y
	LEAQ   (X_PTR)(INC_X*1), X_PTR   // X_PTR = &(X_PTR[ix])
	LEAQ   (Y_PTR)(INC_Y*1), Y_PTR   // Y_PTR = &(Y_PTR[iy])
	MOVQ   incX+56(FP), INC_X        // INC_X = incX
	SHLQ   $4, INC_X                 // INC_X *=  sizeof(complex128)
	MOVQ   incY+64(FP), INC_Y        // INC_Y = incY
	SHLQ   $4, INC_Y                 // INC_Y *=  sizeof(complex128)
	MOVSD  $(-1.0), NEG1
	SHUFPD $0, NEG1, NEG1            // { -1, -1 }
	MOVQ   LEN, TAIL
	ANDQ   $3, TAIL                  // TAIL = n % 4
	SHRQ   $2, LEN                   // LEN = floor( n / 4 )
	JZ     dot_tail                  // if n <= 4 { goto dot_tail }
	MOVAPS NEG1, P_NEG1              // Copy NEG1 to P_NEG1 for pipelining
	LEAQ   (INC_X)(INC_X*2), INCx3_X // INCx3_X = 3 * incX * sizeof(complex128)
	LEAQ   (INC_Y)(INC_Y*2), INCx3_Y // INCx3_Y = 3 * incY * sizeof(complex128)

dot_loop: // do {
	MOVDDUP_XPTR__X3        // X_(i+1) = { real(x[i], real(x[i]) }
	MOVDDUP_XPTR_INCX__X5
	MOVDDUP_XPTR_INCX_2__X7
	MOVDDUP_XPTR_INCx3X__X9

	MOVDDUP_8_XPTR__X2        // X_i = { imag(x[i]), imag(x[i]) }
	MOVDDUP_8_XPTR_INCX__X4
	MOVDDUP_8_XPTR_INCX_2__X6
	MOVDDUP_8_XPTR_INCx3X__X8

	// X_i = { -imag(x[i]), -imag(x[i]) }
	MULPD NEG1, X2
	MULPD P_NEG1, X4
	MULPD NEG1, X6
	MULPD P_NEG1, X8

	// X_j = { imag(y[i]), real(y[i]) }
	MOVUPS (Y_PTR), X10
	MOVUPS (Y_PTR)(INC_Y*1), X11
	MOVUPS (Y_PTR)(INC_Y*2), X12
	MOVUPS (Y_PTR)(INCx3_Y*1), X13

	// X_(i+1) = { imag(a) * real(x[i]), real(a) * real(x[i])  }
	MULPD X10, X3
	MULPD X11, X5
	MULPD X12, X7
	MULPD X13, X9

	// X_j     = { real(y[i]), imag(y[i]) }
	SHUFPD $0x1, X10, X10
	SHUFPD $0x1, X11, X11
	SHUFPD $0x1, X12, X12
	SHUFPD $0x1, X13, X13

	// X_i     = { real(a) * imag(x[i]), imag(a) * imag(x[i])  }
	MULPD X10, X2
	MULPD X11, X4
	MULPD X12, X6
	MULPD X13, X8

	// X_(i+1) = {
	//	imag(result[i]):  imag(a)*real(x[i]) + real(a)*imag(x[i]),
	//	real(result[i]):  real(a)*real(x[i]) - imag(a)*imag(x[i])
	//  }
	ADDSUBPD_X2_X3
	ADDSUBPD_X4_X5
	ADDSUBPD_X6_X7
	ADDSUBPD_X8_X9

	// psum += result[i]
	ADDPD X3, SUM
	ADDPD X5, P_SUM
	ADDPD X7, SUM
	ADDPD X9, P_SUM

	LEAQ (X_PTR)(INC_X*4), X_PTR // X_PTR = &(X_PTR[incX*4])
	LEAQ (Y_PTR)(INC_Y*4), Y_PTR // Y_PTR = &(Y_PTR[incY*4])

	DECQ  LEN
	JNZ   dot_loop   // } while --LEN > 0
	ADDPD P_SUM, SUM // sum += psum
	CMPQ  TAIL, $0   // if TAIL == 0 { return }
	JE    dot_end

dot_tail: // do {
	MOVDDUP_XPTR__X3      // X_(i+1) = { real(x[i], real(x[i]) }
	MOVDDUP_8_XPTR__X2    // X_i = { imag(x[i]), imag(x[i]) }
	MULPD  NEG1, X2       // X_i     = { -imag(x[i])          , -imag(x[i])           }
	MOVUPS (Y_PTR), X10   // X_j     = {  imag(y[i])          ,  real(y[i])           }
	MULPD  X10, X3        // X_(i+1) = {  imag(a) * real(x[i]),  real(a) * real(x[i]) }
	SHUFPD $0x1, X10, X10 // X_j     = {  real(y[i])          ,  imag(y[i])           }
	MULPD  X10, X2        // X_i     = {  real(a) * imag(x[i]),  imag(a) * imag(x[i]) }

	// X_(i+1) = {
	//	imag(result[i]):  imag(a)*real(x[i]) + real(a)*imag(x[i]),
	//	real(result[i]):  real(a)*real(x[i]) - imag(a)*imag(x[i])
	//  }
	ADDSUBPD_X2_X3
	ADDPD X3, SUM      // sum += result[i]
	ADDQ  INC_X, X_PTR // X_PTR += incX
	ADDQ  INC_Y, Y_PTR // Y_PTR += incY
	DECQ  TAIL
	JNZ   dot_tail     // }  while --TAIL > 0

dot_end:
	MOVUPS SUM, sum+88(FP)
	RET