fast_muls.c
void mul_va128128_sse(mzd_t *res, mzd_t *V, mzd_t *A, int clear)
{
//	unsigned arows = A->nrows;
	unsigned acols = A->ncols;
	unsigned vcols = V->ncols;
	unsigned vrows = V->nrows;
	unsigned vcolhi = vcols - 64;
	unsigned acolhi = acols - 64;
//	assert(vcols <= 128);
//	assert(arows == vcols);

	for (int j = 0; j < vrows; j++)
	{
//		__m128i v  = _mm_set_epi64x(mzd_read_bits(V, j, 64, vcolhi), mzd_read_bits(V, j, 0, 64));
		__m128i v  = _mm_loadu_si128((__m128i*)V->rows[j]);

		__m128i m2, m3, rr2, rr3;
		__m128i acc = _mm_setzero_si128();
		__m128i m1  = _mm_set1_epi64x(1);

		for (int i = 0; i < 64; i++)
		{
//			rr2  = _mm_set_epi64x(mzd_read_bits(A, i, 64, acolhi), mzd_read_bits(A, i, 0, 64));
			rr2 = _mm_loadu_si128((__m128i*)A->rows[i]); // FIXME: always correct?
			m2  = _mm_and_si128(v, m1);
			m2  = _mm_cmpeq_epi64(m1, m2);
			if (i < vcolhi)
			{
//				rr3 = _mm_set_epi64x(mzd_read_bits(A, i+64, 64, acolhi), mzd_read_bits(A, i+64, 0, 64));
				rr3 = _mm_loadu_si128((__m128i*)A->rows[i+64]); // FIXME: always correct?
				m3  = _mm_shuffle_epi32(m2,0xFF);
				rr3 = _mm_and_si128(rr3, m3);
				acc = _mm_xor_si128(acc, rr3);
			}
			m2  = _mm_shuffle_epi32(m2,0x00);
			rr2 = _mm_and_si128(rr2, m2);
			acc = _mm_xor_si128(acc, rr2);
			m1  = _mm_slli_epi64(m1, 1);
		}

		if (clear)
		{
			mzd_and_bits(res, j, 0, 64, 0);
			mzd_and_bits(res, j, 64, acolhi, 0);
		}
		mzd_xor_bits(res, j, 0, 64, _mm_extract_epi64(acc, 0));
		mzd_xor_bits(res, j, 64, acolhi, _mm_extract_epi64(acc, 1));
	}

	return;
}
/*
 * All input allocated
 * if clear, reset res before adding the result
 * Computes VxA, where A is at most 128x128
 * Warning: also needs V to have > 64 columns (otherwise, a smaller implementation is selected)
 * Warning: also needs A to have > 64 columns (otherwise, a smaller implementation is selected)
 */
All input allocated
if clear, reset res before adding the result
Computes VxA, where A is at most 128x128
Warning: also needs V to have > 64 columns (otherwise, a smaller implementation is selected)
Warning: also needs A to have > 64 columns (otherwise, a smaller implementation is selected)
All input allocated if clear reset res before adding the result Computes VxA where A is at most 128x128 Warning : also needs V to have > 64 columns ( otherwise a smaller implementation is selected ) Warning : also needs A to have > 64 columns ( otherwise a smaller implementation is selected )
void mul_va128128_sse(mzd_t *res, mzd_t *V, mzd_t *A, int clear)
{
	unsigned acols = A->ncols;
	unsigned vcols = V->ncols;
	unsigned vrows = V->nrows;
	unsigned vcolhi = vcols - 64;
	unsigned acolhi = acols - 64;
	for (int j = 0; j < vrows; j++)
	{
		__m128i v  = _mm_loadu_si128((__m128i*)V->rows[j]);
		__m128i m2, m3, rr2, rr3;
		__m128i acc = _mm_setzero_si128();
		__m128i m1  = _mm_set1_epi64x(1);
		for (int i = 0; i < 64; i++)
		{
			rr2 = _mm_loadu_si128((__m128i*)A->rows[i]); 
			m2  = _mm_and_si128(v, m1);
			m2  = _mm_cmpeq_epi64(m1, m2);
			if (i < vcolhi)
			{
				rr3 = _mm_loadu_si128((__m128i*)A->rows[i+64]); 
				m3  = _mm_shuffle_epi32(m2,0xFF);
				rr3 = _mm_and_si128(rr3, m3);
				acc = _mm_xor_si128(acc, rr3);
			}
			m2  = _mm_shuffle_epi32(m2,0x00);
			rr2 = _mm_and_si128(rr2, m2);
			acc = _mm_xor_si128(acc, rr2);
			m1  = _mm_slli_epi64(m1, 1);
		}
		if (clear)
		{
			mzd_and_bits(res, j, 0, 64, 0);
			mzd_and_bits(res, j, 64, acolhi, 0);
		}
		mzd_xor_bits(res, j, 0, 64, _mm_extract_epi64(acc, 0));
		mzd_xor_bits(res, j, 64, acolhi, _mm_extract_epi64(acc, 1));
	}
	return;
}
void mul_va128128_sse ( mzd_t * res , mzd_t * V , mzd_t * A , int clear ) { unsigned acols = A -> ncols ; unsigned vcols = V -> ncols ; unsigned vrows = V -> nrows ; unsigned vcolhi = vcols - 64 ; unsigned acolhi = acols - 64 ; for ( int j = 0 ; j < vrows ; j ++ ) { __m128i v = _mm_loadu_si128 ( ( __m128i * ) V -> rows [ j ] ) ; __m128i m2 , m3 , rr2 , rr3 ; __m128i acc = _mm_setzero_si128 ( ) ; __m128i m1 = _mm_set1_epi64x ( 1 ) ; for ( int i = 0 ; i < 64 ; i ++ ) { rr2 = _mm_loadu_si128 ( ( __m128i * ) A -> rows [ i ] ) ; m2 = _mm_and_si128 ( v , m1 ) ; m2 = _mm_cmpeq_epi64 ( m1 , m2 ) ; if ( i < vcolhi ) { rr3 = _mm_loadu_si128 ( ( __m128i * ) A -> rows [ i + 64 ] ) ; m3 = _mm_shuffle_epi32 ( m2 , 0xFF ) ; rr3 = _mm_and_si128 ( rr3 , m3 ) ; acc = _mm_xor_si128 ( acc , rr3 ) ; } m2 = _mm_shuffle_epi32 ( m2 , 0x00 ) ; rr2 = _mm_and_si128 ( rr2 , m2 ) ; acc = _mm_xor_si128 ( acc , rr2 ) ; m1 = _mm_slli_epi64 ( m1 , 1 ) ; } if ( clear ) { mzd_and_bits ( res , j , 0 , 64 , 0 ) ; mzd_and_bits ( res , j , 64 , acolhi , 0 ) ; } mzd_xor_bits ( res , j , 0 , 64 , _mm_extract_epi64 ( acc , 0 ) ) ; mzd_xor_bits ( res , j , 64 , acolhi , _mm_extract_epi64 ( acc , 1 ) ) ; } return ; }
All input allocated
 if clear, reset res before adding the result
 Computes VxA, where A is at most 128x128
 Warning: also needs V to have > 64 columns (otherwise, a smaller implementation is selected)
 Warning: also needs A to have > 64 columns (otherwise, a smaller implementation is selected)
All input allocated if clear reset res before adding the result Computes VxA where A is at most 128x128 Warning : also needs V to have > 64 columns ( otherwise a smaller implementation is selected ) Warning : also needs A to have > 64 columns ( otherwise a smaller implementation is selected )
//	unsigned arows = A->nrows; //	assert(vcols <= 128); //	assert(arows == vcols); //		__m128i v  = _mm_set_epi64x(mzd_read_bits(V, j, 64, vcolhi), mzd_read_bits(V, j, 0, 64)); //			rr2  = _mm_set_epi64x(mzd_read_bits(A, i, 64, acolhi), mzd_read_bits(A, i, 0, 64)); // FIXME: always correct? //				rr3 = _mm_set_epi64x(mzd_read_bits(A, i+64, 64, acolhi), mzd_read_bits(A, i+64, 0, 64)); // FIXME: always correct?