emulated_arm64.go

     1  // Copyright 2026 The Go Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  //go:build goexperiment.simd && arm64
     6  
     7  package bridge
     8  
     9  import "simd/archsimd"
    10  
    11  func (xx Uint64x2nclm) CarrylessMultiplyEven(yy Uint64x2nclm) Uint64x2nclm {
    12  	x := archsimd.Uint64x2(xx)
    13  	y := archsimd.Uint64x2(yy)
    14  
    15  	return Uint64x2nclm(carrylessMultiply_(x, y))
    16  }
    17  
    18  func (xx Uint64x2nclm) CarrylessMultiplyOdd(yy Uint64x2nclm) Uint64x2nclm {
    19  	x := archsimd.Uint64x2(xx)
    20  	y := archsimd.Uint64x2(yy)
    21  
    22  	x = x.SetElem(0, x.GetElem(1))
    23  	y = y.SetElem(0, y.GetElem(1))
    24  	return Uint64x2nclm(carrylessMultiply_(x, y))
    25  }
    26  
    27  func new64x2(lo, hi uint64) archsimd.Uint64x2 {
    28  	return archsimd.Uint64x2{}.SetElem(0, lo).SetElem(1, hi)
    29  }
    30  
    31  // These masks all have 4 zeroes between 1s.
    32  var m0_ = new64x2(0x1084210842108421, 0x2108421084210842)
    33  var m1_ = new64x2(0x2108421084210842, 0x4210842108421084)
    34  var m2_ = new64x2(0x4210842108421084, 0x8421084210842108)
    35  var m3_ = new64x2(0x8421084210842108, 0x0842108421084210)
    36  var m4_ = new64x2(0x0842108421084210, 0x1084210842108421)
    37  
    38  // Selects the middle 64 bits of a 128-bit simd value
    39  var middle = new64x2(0xffffffff00000000, 0x00000000ffffffff)
    40  
    41  // mwl_ is a 64x64 into 128 multiply that is missing
    42  // some carries that we don't need for CLMUL emulation.
    43  // The high 64 bits of each input are ignored.
    44  // Also just for fun, accumulate sums with Xor.
    45  func mwl(x, y archsimd.Uint64x2) archsimd.Uint64x2 {
    46  	// reshape input into Uint32x4
    47  	// input is  {a b _ _}.mwl_{c d _ _}
    48  	// need the sum of
    49  	// ac0_ac1
    50  	//   0 ad0_ad1
    51  	//   0 bc0_bc1
    52  	//   0   0 bd0_bd1
    53  	// This "sum" is where the carries (not propagated
    54  	// across lanes) are lost.
    55  	ab__ := x.ReshapeToUint32s()
    56  	cd__ := y.ReshapeToUint32s()
    57  	ac0_ac1_bd0_bd1 := ab__.MulWidenLo(cd__)
    58  
    59  	dc__ := y.RotateAllLeft(32).ReshapeToUint32s()
    60  	ad0_ad1_bc0_bc1 := ab__.MulWidenLo(dc__)
    61  	//
    62  	// have        ad0, ad1, bc0, bc1
    63  	// want        0, ad0+bc0, ad1+bc1, 0
    64  	// to add to    ac0_ac1_bd0_bd1
    65  	//
    66  	// swap 64-bit halves of ad0_ad1_bc0_bc1
    67  	// to get   bc0_bc1_ad0_ad1
    68  	bc0_bc1_ad0_ad1 := archsimd.Uint64x2{}.SetElem(0, ad0_ad1_bc0_bc1.GetElem(1)).SetElem(1, ad0_ad1_bc0_bc1.GetElem(0))
    69  
    70  	// added to ad0_ad1_bc0_bc1 yields
    71  	//   bc0+ad0, bc1+ad1, bc0+ad0, bc1+ad1
    72  	// rotate 32 (within the two 64-bit elements) yields
    73  	//   bc1+ad1, bc0+ad0, bc1+ad1, bc0+ad0
    74  	// and then intersect with mask:
    75  	//   0      , bc0+ad0, bc1+ad1, 0
    76  	//
    77  	// use xor to make it a worse multiply
    78  	zzz_adPbc0_adPbc1_zzz := bc0_bc1_ad0_ad1.Xor(ad0_ad1_bc0_bc1).RotateAllLeft(32).And(middle)
    79  	return ac0_ac1_bd0_bd1.Xor(zzz_adPbc0_adPbc1_zzz)
    80  }
    81  
    82  // carrylessMultiply is constant time carrless multiply implemented with an
    83  // absurd number of multiplication given that the emulation platforms only have
    84  // 32x32 into 64, it might make sense to rework this into that primitive, but,
    85  // for now this works and is easily tested in scalar Go.
    86  func carrylessMultiply_(x, y archsimd.Uint64x2) archsimd.Uint64x2 {
    87  
    88  	// This by masking the two inputs into 5 thinned inputs, with
    89  	// 4 zeroes separating any 2 set bits.  Multiply will potentially
    90  	// set more bits with addition of overlapping terms, however this
    91  	// technique allows as many as 31 additions (filling all 4 separation
    92  	// positions with 1) without perturbing the bits we care about.  Since
    93  	// there's at most 13 set bits in a thinned input, 31 is not a problem.
    94  	// If there were only 3 set bits, there are 16 1s per thinned input and
    95  	// only 15 additions can be tolerated -- so that's not possible.
    96  
    97  	// This is also discussed at
    98  	// https://timtaubert.de/blog/2017/06/verified-binary-multiplication-for-ghash/
    99  
   100  	x0 := x.And(m0_)
   101  	x1 := x.And(m1_)
   102  	x2 := x.And(m2_)
   103  	x3 := x.And(m3_)
   104  	x4 := x.And(m4_)
   105  
   106  	y0 := y.And(m0_)
   107  	y1 := y.And(m1_)
   108  	y2 := y.And(m2_)
   109  	y3 := y.And(m3_)
   110  	y4 := y.And(m4_)
   111  
   112  	var z archsimd.Uint64x2
   113  	// for a given line, combining (xI).mwl_(yJ) terms, I+J == K mod 5; mask index = K
   114  	z = (mwl(x0, y0)).Xor(mwl(x1, y4)).Xor(mwl(x4, y1)).Xor(mwl(x2, y3)).Xor(mwl(x3, y2)).And(m0_)
   115  	z = (mwl(x3, y3)).Xor(mwl(x2, y4)).Xor(mwl(x4, y2)).Xor(mwl(x0, y1)).Xor(mwl(x1, y0)).And(m1_).Or(z)
   116  	z = (mwl(x1, y1)).Xor(mwl(x3, y4)).Xor(mwl(x4, y3)).Xor(mwl(x0, y2)).Xor(mwl(x2, y0)).And(m2_).Or(z)
   117  	z = (mwl(x4, y4)).Xor(mwl(x0, y3)).Xor(mwl(x3, y0)).Xor(mwl(x1, y2)).Xor(mwl(x2, y1)).And(m3_).Or(z)
   118  	z = (mwl(x2, y2)).Xor(mwl(x0, y4)).Xor(mwl(x4, y0)).Xor(mwl(x1, y3)).Xor(mwl(x3, y1)).And(m4_).Or(z)
   119  
   120  	return z
   121  }
   122
View as plain text