--- /dev/null
+// 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 MOVSHDUP_X3_X2 LONG $0xD3160FF3 // MOVSHDUP X3, X2
+#define MOVSHDUP_X5_X4 LONG $0xE5160FF3 // MOVSHDUP X5, X4
+#define MOVSHDUP_X7_X6 LONG $0xF7160FF3 // MOVSHDUP X7, X6
+#define MOVSHDUP_X9_X8 LONG $0x160F45F3; BYTE $0xC1 // MOVSHDUP X9, X8
+
+#define MOVSLDUP_X3_X3 LONG $0xDB120FF3 // MOVSLDUP X3, X3
+#define MOVSLDUP_X5_X5 LONG $0xED120FF3 // MOVSLDUP X5, X5
+#define MOVSLDUP_X7_X7 LONG $0xFF120FF3 // MOVSLDUP X7, X7
+#define MOVSLDUP_X9_X9 LONG $0x120F45F3; BYTE $0xC9 // MOVSLDUP X9, X9
+
+#define ADDSUBPS_X2_X3 LONG $0xDAD00FF2 // ADDSUBPS X2, X3
+#define ADDSUBPS_X4_X5 LONG $0xECD00FF2 // ADDSUBPS X4, X5
+#define ADDSUBPS_X6_X7 LONG $0xFED00FF2 // ADDSUBPS X6, X7
+#define ADDSUBPS_X8_X9 LONG $0xD00F45F2; BYTE $0xC8 // ADDSUBPS 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
+
+// func DotuInc(x, y []complex64, n, incX, incY, ix, iy uintptr) (sum complex64)
+TEXT ·DotuInc(SB), NOSPLIT, $0
+ MOVQ x_base+0(FP), X_PTR // X_PTR = &x
+ MOVQ y_base+24(FP), Y_PTR // Y_PTR = &y
+ PXOR SUM, SUM // SUM = 0
+ PXOR P_SUM, P_SUM // P_SUM = 0
+ MOVQ n+48(FP), LEN // LEN = n
+ CMPQ LEN, $0 // if LEN == 0 { return }
+ JE dotu_end
+ MOVQ ix+72(FP), INC_X
+ MOVQ iy+80(FP), INC_Y
+ LEAQ (X_PTR)(INC_X*8), X_PTR // X_PTR = &(X_PTR[ix])
+ LEAQ (Y_PTR)(INC_Y*8), Y_PTR // Y_PTR = &(Y_PTR[iy])
+ MOVQ incX+56(FP), INC_X // INC_X = incX * sizeof(complex64)
+ SHLQ $3, INC_X
+ MOVQ incY+64(FP), INC_Y // INC_Y = incY * sizeof(complex64)
+ SHLQ $3, INC_Y
+
+ MOVQ LEN, TAIL
+ ANDQ $3, TAIL // TAIL = LEN % 4
+ SHRQ $2, LEN // LEN = floor( LEN / 4 )
+ JZ dotu_tail // if TAIL == 0 { goto dotu_tail }
+
+ LEAQ (INC_X)(INC_X*2), INCx3_X // INCx3_X = INC_X * 3
+ LEAQ (INC_Y)(INC_Y*2), INCx3_Y // INCx3_Y = INC_Y * 3
+
+dotu_loop: // do {
+ MOVSD (X_PTR), X3 // X_i = { imag(x[i]), real(x[i]) }
+ MOVSD (X_PTR)(INC_X*1), X5
+ MOVSD (X_PTR)(INC_X*2), X7
+ MOVSD (X_PTR)(INCx3_X*1), X9
+
+ // X_(i-1) = { imag(x[i]), imag(x[i]) }
+ MOVSHDUP_X3_X2
+ MOVSHDUP_X5_X4
+ MOVSHDUP_X7_X6
+ MOVSHDUP_X9_X8
+
+ // X_i = { real(x[i]), real(x[i]) }
+ MOVSLDUP_X3_X3
+ MOVSLDUP_X5_X5
+ MOVSLDUP_X7_X7
+ MOVSLDUP_X9_X9
+
+ // X_j = { imag(y[i]), real(y[i]) }
+ MOVSD (Y_PTR), X10
+ MOVSD (Y_PTR)(INC_Y*1), X11
+ MOVSD (Y_PTR)(INC_Y*2), X12
+ MOVSD (Y_PTR)(INCx3_Y*1), X13
+
+ // X_i = { imag(y[i]) * real(x[i]), real(y[i]) * real(x[i]) }
+ MULPS X10, X3
+ MULPS X11, X5
+ MULPS X12, X7
+ MULPS X13, X9
+
+ // X_j = { real(y[i]), imag(y[i]) }
+ SHUFPS $0xB1, X10, X10
+ SHUFPS $0xB1, X11, X11
+ SHUFPS $0xB1, X12, X12
+ SHUFPS $0xB1, X13, X13
+
+ // X_(i-1) = { real(y[i]) * imag(x[i]), imag(y[i]) * imag(x[i]) }
+ MULPS X10, X2
+ MULPS X11, X4
+ MULPS X12, X6
+ MULPS X13, X8
+
+ // X_i = {
+ // imag(result[i]): imag(y[i]) * real(x[i]) + real(y[i]) * imag(x[i]),
+ // real(result[i]): real(y[i]) * real(x[i]) - imag(y[i]) * imag(x[i]) }
+ ADDSUBPS_X2_X3
+ ADDSUBPS_X4_X5
+ ADDSUBPS_X6_X7
+ ADDSUBPS_X8_X9
+
+ // SUM += X_i
+ ADDPS X3, SUM
+ ADDPS X5, P_SUM
+ ADDPS X7, SUM
+ ADDPS X9, P_SUM
+
+ LEAQ (X_PTR)(INC_X*4), X_PTR // X_PTR = &(X_PTR[INC_X*4])
+ LEAQ (Y_PTR)(INC_Y*4), Y_PTR // Y_PTR = &(Y_PTR[INC_Y*4])
+
+ DECQ LEN
+ JNZ dotu_loop // } while --LEN > 0
+
+ ADDPS P_SUM, SUM // SUM = { P_SUM + SUM }
+ CMPQ TAIL, $0 // if TAIL == 0 { return }
+ JE dotu_end
+
+dotu_tail: // do {
+ MOVSD (X_PTR), X3 // X_i = { imag(x[i]), real(x[i]) }
+ MOVSHDUP_X3_X2 // X_(i-1) = { imag(x[i]), imag(x[i]) }
+ MOVSLDUP_X3_X3 // X_i = { real(x[i]), real(x[i]) }
+ MOVUPS (Y_PTR), X10 // X_j = { imag(y[i]), real(y[i]) }
+ MULPS X10, X3 // X_i = { imag(y[i]) * real(x[i]), real(y[i]) * real(x[i]) }
+ SHUFPS $0x1, X10, X10 // X_j = { real(y[i]), imag(y[i]) }
+ MULPS X10, X2 // X_(i-1) = { real(y[i]) * imag(x[i]), imag(y[i]) * imag(x[i]) }
+
+ // X_i = {
+ // imag(result[i]): imag(y[i])*real(x[i]) + real(y[i])*imag(x[i]),
+ // real(result[i]): real(y[i])*real(x[i]) - imag(y[i])*imag(x[i]) }
+ ADDSUBPS_X2_X3
+ ADDPS X3, SUM // SUM += X_i
+ ADDQ INC_X, X_PTR // X_PTR += INC_X
+ ADDQ INC_Y, Y_PTR // Y_PTR += INC_Y
+ DECQ TAIL
+ JNZ dotu_tail // } while --TAIL > 0
+
+dotu_end:
+ MOVSD SUM, sum+88(FP) // return SUM
+ RET