github.com/ericlagergren/ctb@v0.0.0-20220810041818-96749d9c394d/lll/f64.go (about)

     1  package lll
     2  
     3  import (
     4  	"math"
     5  
     6  	"gonum.org/v1/gonum/floats"
     7  )
     8  
     9  // Reduction64 computes the Lenstra–Lenstra–Lovász
    10  // lattice basis reduction algorithm.
    11  //
    12  // B is a lattice basis
    13  //    b0, b1, ... bn in Z^m
    14  // delta must be in (1/4, 1), typically 3/4.
    15  func Reduction64(delta float64, B [][]float64) [][]float64 {
    16  	if delta < 1/4 || delta > 1 {
    17  		panic("delta out of range")
    18  	}
    19  	Bstar := gramSchmidt64(nil, B)
    20  	mu := func(i, j int) float64 {
    21  		return projCoff64(Bstar[j], B[i])
    22  	}
    23  	n := len(B)
    24  	k := 1
    25  	for k < n {
    26  		for j := k - 1; j >= 0; j-- {
    27  			mukj := mu(k, j)
    28  			if math.Abs(mukj) > 0.5 {
    29  				bj := scale64(nil, B[j], math.Round(mukj))
    30  				B[k] = sub64(B[k], B[k], bj)
    31  				Bstar = gramSchmidt64(Bstar, B)
    32  			}
    33  		}
    34  		if sdot64(Bstar[k]) >= (delta-math.Pow(mu(k, k-1), 2))*sdot64(Bstar[k-1]) {
    35  			k++
    36  		} else {
    37  			B[k], B[k-1] = B[k-1], B[k]
    38  			Bstar = gramSchmidt64(Bstar, B)
    39  			k--
    40  			if k < 1 {
    41  				k = 1
    42  			}
    43  		}
    44  	}
    45  	return B
    46  }
    47  
    48  func gramSchmidt64(u, v [][]float64) [][]float64 {
    49  	u = u[:0]
    50  	for _, vi := range v {
    51  		ui := vi
    52  		for _, uj := range u {
    53  			// ui -= uj*vi
    54  			uj = proj64(nil, uj, vi)
    55  			ui = sub64(nil, ui, uj)
    56  		}
    57  		if len(ui) > 0 {
    58  			u = append(u, ui)
    59  		}
    60  	}
    61  	return u
    62  }
    63  
    64  func scale64(z, x []float64, c float64) []float64 {
    65  	z = zmake64(z, len(x))
    66  	return floats.ScaleTo(z, c, x)
    67  }
    68  
    69  func mul64(z, x, y []float64) []float64 {
    70  	z = zmake64(z, len(x))
    71  	return floats.MulTo(z, x, y)
    72  }
    73  
    74  func sub64(z, x, y []float64) []float64 {
    75  	z = zmake64(z, len(x))
    76  	return floats.SubTo(z, x, y)
    77  }
    78  
    79  func proj64(z, x, y []float64) []float64 {
    80  	z = zmake64(z, len(x))
    81  	return scale64(z, x, projCoff64(x, y))
    82  }
    83  
    84  func projCoff64(x, y []float64) float64 {
    85  	return dot64(x, y) / sdot64(x)
    86  }
    87  
    88  func dot64(x, y []float64) float64 {
    89  	return floats.Dot(x, y)
    90  }
    91  
    92  func sdot64(x []float64) float64 {
    93  	return dot64(x, x)
    94  }
    95  
    96  func zmake64(z []float64, n int) []float64 {
    97  	if n <= cap(z) {
    98  		return z[:n]
    99  	}
   100  	return make([]float64, n)
   101  }