<commit_before>package factorlib

import (
	"fmt"
	"math/big"
	"math/rand"
	"github.com/randall77/factorlib/matrix"
)

func init() {
	factorizers["qs"] = qs
}

func qs(n *big.Int, rnd *rand.Rand) []*big.Int {
	// Quadratic sieve:
	// Define f(x) = x^2 - n.  We wish to find x such that f(x) is smooth.
	// We will take x starting at ceil(sqrt(n)).

	// is f(x) divisible by p^k?
	// calculating mod p^k, f(x + p^k) = (x+p^k)^2 - n = x^2 - n = f(x)
	// so if f(x) is divisible by p^k, so is f(x + i * p^k) for all i.
	//
	// if x_1 and x_2 are the two solutions to f(x) == 0 mod p^k, then all x which
	// have f(x) divisible by p^k are x_1 + i * p^k and x_2 + i * p^k.
	//
	// f(x) == 0 mod p^k is x^2 == n mod p^k.

	// So the factor base contains small primes and powers of primes.
	// For each factor base element p^k, we store two offsets satisfying
	// f(start + off1) == 0 mod p^k and f(start + off2) == 0 mod p^k.
	// Then to sieve, we just add log_p to sieve[off1], sieve[off1+p^k], sieve[off1+2*p^k], ...
	// and the same for off2.  When we reach the end of the sieve, we update off1 and off2
	// for the next sieve.

	// the sieve returns x such that f(x) = x^2 - n = prod_i factors[i]
	// all factors <= B
	// factors are sorted in increasing order.

	// first, pick a factor base
	fb, a := makeFactorBase(n)
	if a != nil {
		return []*big.Int{a, new(big.Int).Div(n, a)}
	}

	maxp := fb[len(fb)-1]
	bigmaxp2 := new(big.Int).SetInt64(maxp*maxp)
	fmt.Printf("maxp=%d maxp2=%d len(fb)=%d\n", maxp, bigmaxp2, len(fb))

	// The x we start with
	start := sqrtCeil(n)

	si := makeSieveInfo(n, start, fb, rnd)

	sievelen := int(fb[len(fb)-1])
	if sievelen < 10000 {
		sievelen = 10000
	}

	sieve := make([]byte, sievelen)

	// temporary storage
	x := new(big.Int)
	y := new(big.Int)
	bigz := new(big.Int)
	r := new(big.Int)
	bigp := new(big.Int)
	t := new(big.Int)
	var factors []uint

	// matrix is used to do gaussian elimination on mod 2 exponents.
	m := matrix.NewBitMatrix(uint(len(fb)))

	// largeprimes records instances of the equation
	//   x^2 = prod(f) * p mod n
	// where f are small primes (indexes into factor base) and
	// p is a large prime.
	type largerecord struct {
		x *big.Int
		f []uint
	}
	largeprimes := map[int64]largerecord{}

	for {
		// clear sieve
		for i := 0; i < sievelen; i++ {
			sieve[i] = 0
		}
		
		// increment sieve entries for f(x) that are divisible
		// by each factor base prime.
		for i := range si {
			f := &si[i]
			pk := int(f.pk)
			lg_p := f.lg_p
			j := int(f.off1)
			for ; j < sievelen; j += pk {
				sieve[j] += lg_p
			}
			f.off1 = int32(j - sievelen) // for next time
			
			// same for off2
			j = int(f.off2)<commit_msg>Use renamed matrix.<commit_after>package factorlib

import (
	"fmt"
	"math/big"
	"math/rand"
	"github.com/randall77/factorlib/linear"
)

func init() {
	factorizers["qs"] = qs
}

func qs(n *big.Int, rnd *rand.Rand) []*big.Int {
	// Quadratic sieve:
	// Define f(x) = x^2 - n.  We wish to find x such that f(x) is smooth.
	// We will take x starting at ceil(sqrt(n)).

	// is f(x) divisible by p^k?
	// calculating mod p^k, f(x + p^k) = (x+p^k)^2 - n = x^2 - n = f(x)
	// so if f(x) is divisible by p^k, so is f(x + i * p^k) for all i.
	//
	// if x_1 and x_2 are the two solutions to f(x) == 0 mod p^k, then all x which
	// have f(x) divisible by p^k are x_1 + i * p^k and x_2 + i * p^k.
	//
	// f(x) == 0 mod p^k is x^2 == n mod p^k.

	// So the factor base contains small primes and powers of primes.
	// For each factor base element p^k, we store two offsets satisfying
	// f(start + off1) == 0 mod p^k and f(start + off2) == 0 mod p^k.
	// Then to sieve, we just add log_p to sieve[off1], sieve[off1+p^k], sieve[off1+2*p^k], ...
	// and the same for off2.  When we reach the end of the sieve, we update off1 and off2
	// for the next sieve.

	// the sieve returns x such that f(x) = x^2 - n = prod_i factors[i]
	// all factors <= B
	// factors are sorted in increasing order.

	// first, pick a factor base
	fb, a := makeFactorBase(n)
	if a != nil {
		return []*big.Int{a, new(big.Int).Div(n, a)}
	}

	maxp := fb[len(fb)-1]
	bigmaxp2 := new(big.Int).SetInt64(maxp*maxp)
	fmt.Printf("maxp=%d maxp2=%d len(fb)=%d\n", maxp, bigmaxp2, len(fb))

	// The x we start with
	start := sqrtCeil(n)

	si := makeSieveInfo(n, start, fb, rnd)

	sievelen := int(fb[len(fb)-1])
	if sievelen < 10000 {
		sievelen = 10000
	}

	sieve := make([]byte, sievelen)

	// temporary storage
	x := new(big.Int)
	y := new(big.Int)
	bigz := new(big.Int)
	r := new(big.Int)
	bigp := new(big.Int)
	t := new(big.Int)
	var factors []uint

	// matrix is used to do gaussian elimination on mod 2 exponents.
	m := linear.NewMatrix(uint(len(fb)))

	// largeprimes records instances of the equation
	//   x^2 = prod(f) * p mod n
	// where f are small primes (indexes into factor base) and
	// p is a large prime.
	type largerecord struct {
		x *big.Int
		f []uint
	}
	largeprimes := map[int64]largerecord{}

	for {
		// clear sieve
		for i := 0; i < sievelen; i++ {
			sieve[i] = 0
		}
		
		// increment sieve entries for f(x) that are divisible
		// by each factor base prime.
		for i := range si {
			f := &si[i]
			pk := int(f.pk)
			lg_p := f.lg_p
			j := int(f.off1)
			for ; j < sievelen; j += pk {
				sieve[j] += lg_p
			}
			f.off1 = int32(j - sievelen) // for next time
			
			// same for off2
			j = int(f.off2)
{'fraction_non_alphanumeric': 0.16309412861137, 'fraction_numerical': 0.019757688723205966, 'mean_word_length': 2.6703146374829, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3963494', 'n_tokens_mistral': 2392, 'n_tokens_neox': 2212, 'n_words': 723}