lll.go raw

   1  package gnarlring
   2  
   3  import "math"
   4  
   5  // lllReduce reduces a lattice basis in-place using the LLL algorithm
   6  // with parameter delta (typically 0.99). The input basis is a slice
   7  // of d columns, each a vector of length d.
   8  func lllReduce(basis [][]int64, delta float64) {
   9  	d := len(basis)
  10  	if d == 0 {
  11  		return
  12  	}
  13  
  14  	mu := make([][]float64, d)
  15  	for i := range mu {
  16  		mu[i] = make([]float64, d)
  17  	}
  18  	bStar := make([][]float64, d)
  19  	bStarNorm := make([]float64, d)
  20  
  21  	computeGS(basis, mu, bStar, bStarNorm)
  22  
  23  	k := 1
  24  	for k < d {
  25  		for j := k - 1; j >= 0; j-- {
  26  			r := int(math.Round(mu[k][j]))
  27  			if r != 0 {
  28  				addScaled(basis[k], basis[j], -r)
  29  				updateMuRow(mu, basis, bStar, bStarNorm, k, j)
  30  				updateBStar(bStar, mu, basis, bStarNorm, k)
  31  			}
  32  		}
  33  
  34  		if bStarNorm[k] >= (delta-mu[k][k-1]*mu[k][k-1])*bStarNorm[k-1] {
  35  			k++
  36  		} else {
  37  			basis[k], basis[k-1] = basis[k-1], basis[k]
  38  			for i := k - 1; i < d; i++ {
  39  				updateBStar(bStar, mu, basis, bStarNorm, i)
  40  				updateMuCol(mu, basis, bStar, bStarNorm, i, d)
  41  			}
  42  			if k > 1 {
  43  				k--
  44  			}
  45  		}
  46  	}
  47  }
  48  
  49  func computeGS(basis [][]int64, mu [][]float64, bStar [][]float64, bStarNorm []float64) {
  50  	d := len(basis)
  51  	for i := 0; i < d; i++ {
  52  		bStar[i] = make([]float64, d)
  53  		for j := 0; j < d; j++ {
  54  			bStar[i][j] = float64(basis[i][j])
  55  		}
  56  	}
  57  	for i := 0; i < d; i++ {
  58  		for j := 0; j < i; j++ {
  59  			mu[i][j] = dotFloat(bStar[i], bStar[j]) / bStarNorm[j]
  60  			for k := 0; k < d; k++ {
  61  				bStar[i][k] -= mu[i][j] * bStar[j][k]
  62  			}
  63  		}
  64  		bStarNorm[i] = dotFloat(bStar[i], bStar[i])
  65  	}
  66  }
  67  
  68  func updateBStar(bStar [][]float64, mu [][]float64, basis [][]int64, bStarNorm []float64, i int) {
  69  	d := len(basis)
  70  	for k := 0; k < d; k++ {
  71  		bStar[i][k] = float64(basis[i][k])
  72  	}
  73  	for j := 0; j < i; j++ {
  74  		mu[i][j] = dotFloat(bStar[i], bStar[j]) / bStarNorm[j]
  75  		for k := 0; k < d; k++ {
  76  			bStar[i][k] -= mu[i][j] * bStar[j][k]
  77  		}
  78  	}
  79  	bStarNorm[i] = dotFloat(bStar[i], bStar[i])
  80  }
  81  
  82  func updateMuCol(mu [][]float64, basis [][]int64, bStar [][]float64, bStarNorm []float64, i, d int) {
  83  	for t := i + 1; t < d; t++ {
  84  		if bStarNorm[i] > 1e-20 {
  85  			mu[t][i] = int64DotFloat(basis[t], bStar[i]) / bStarNorm[i]
  86  		} else {
  87  			mu[t][i] = 0
  88  		}
  89  	}
  90  }
  91  
  92  func updateMuRow(mu [][]float64, basis [][]int64, bStar [][]float64, bStarNorm []float64, k, j int) {
  93  	for jj := 0; jj <= j; jj++ {
  94  		if bStarNorm[jj] > 1e-20 {
  95  			mu[k][jj] = int64DotFloat(basis[k], bStar[jj]) / bStarNorm[jj]
  96  		} else {
  97  			mu[k][jj] = 0
  98  		}
  99  	}
 100  }
 101  
 102  func dotFloat(a, b []float64) float64 {
 103  	var sum float64
 104  	for i := range a {
 105  		sum += a[i] * b[i]
 106  	}
 107  	return sum
 108  }
 109  
 110  func int64DotFloat(a []int64, b []float64) float64 {
 111  	var sum float64
 112  	for i := range a {
 113  		sum += float64(a[i]) * b[i]
 114  	}
 115  	return sum
 116  }
 117  
 118  func addScaled(dst, src []int64, r int) {
 119  	for i := range dst {
 120  		dst[i] += int64(r) * src[i]
 121  	}
 122  }
 123