CFRAC and SQUFOF methods to completely factorise any number into its prime roots. [SLOW AFTER 30+ digits]

cfrac.go

package main

import (
	"fmt"
	"log"
	"math/big"
	"os"
	"time"
)

var mone = big.NewInt(int64(-1))
var zero = big.NewInt(int64(0))
var one = big.NewInt(int64(1))
var two = big.NewInt(int64(2))
var three = big.NewInt(int64(3))
var four = big.NewInt(int64(4))

type kv struct {
	Key   string
	Value int
}

type Ring struct {
	Items [2]*big.Int
}

// Shift items left and add replace last index with n
func (r *Ring) Push(n *big.Int) {
	r.Items[0] = r.Items[1]
	r.Items[1] = n
}

// Return items in ring individually (n-1, n)
func (r Ring) Get() (*big.Int, *big.Int) {
	return r.Items[0], r.Items[1]
}

func main() {
	if len(os.Args[1:]) < 1 {
		log.Fatal("Not enough arguments")
	}
	if len(os.Args[1:]) > 2 {
		log.Fatal("Too many arguments")
	}

	N := new(big.Int)
	k := big.NewInt(int64(1))
	N.SetString(os.Args[1], 10)
	if len(os.Args[1:]) == 2 {
		k.SetString(os.Args[2], 10)
	}

	start := time.Now()
	r := CFRAC(N, k)
	//r := SQUFOF(N)
	elapsed := time.Since(start)
	r = Quicksort(r)
	fmt.Printf("%v = %v (%v) (%d)\n", N, r, elapsed, len(N.String()))
}

/*
Morrison-Brillhart method which is an Extension of
the continuos fraction prime factorisation
Expands sqrt(kN) to find congruences A^2/B^2
Using a series of equations determine A-Q pairs
Using these A-Q pairs find a suitable match such that
GCD(A[n-1]-sqrt(Q[n]), N) gives a non trivial factor of N
*/
func CFRAC(N, k *big.Int) []*big.Int {
	if prime := N.ProbablyPrime(10); prime {
		return []*big.Int{N}
	}

	// Initialise factors array
	factors := []*big.Int{}

	// Take out factor of 2 on even numbers > 2
	Num := new(big.Int)
	Num = Num.Set(N)
	even := new(big.Int)
	even = even.Mod(Num, two)
	if N.Cmp(two) > 0 && even.Cmp(zero) == 0 {
		factors = append(factors, two)
		Num = Num.Div(Num, two)
	}

	kN := new(big.Int)
	kN = kN.Mul(k, Num)
	tsqkN := new(big.Int)
	tsqkN = tsqkN.Sqrt(kN).Mul(tsqkN, two)

	// Store expansion data Ring structs
	var A, Q, P, r, q Ring
	A.Push(zero)
	A.Push(one)
	Q.Push(kN)
	Q.Push(one)
	P.Push(zero)
	g := new(big.Int)
	g = g.Sqrt(kN)
	r.Push(g)
	q.Push(g)

	// Initialise arrays
	D, B := new(big.Int), new(big.Int)
	FD, FB := []*big.Int{}, []*big.Int{}

	// Period of sqrt(kN) is usually sqrt(kN)
	for i := big.NewInt(int64(0)); i.Cmp(g) <= 0; i.Add(i, one) {
		// g+P[n]=q[n]Q[n]+r[n]
		_, Pn := P.Get()
		gPn := new(big.Int)
		gPn = gPn.Add(g, Pn)

		Qm, Qn := Q.Get()

		// Catch Q error and exit
		if i.Cmp(zero) > 0 && (Qn.Cmp(zero) <= 0 || Qn.Cmp(tsqkN) >= 0) {
			break
		}

		qp, rp := new(big.Int), new(big.Int)
		qp, rp = qp.QuoRem(gPn, Qn, rp)
		q.Push(qp)
		r.Push(rp)

		// A[n]=q[n]A[n-1]+A[n-2]
		Am, An := A.Get()
		qA := new(big.Int)
		qA = qA.Mul(qp, An).Add(qA, Am).Mod(qA, N)
		A.Push(qA)

		// g+P[n+1]=2g-r[n]
		tg := big.NewInt(int64(2))
		rm, rn := r.Get()
		tg = tg.Mul(tg, g).Sub(tg, rn).Sub(tg, g)
		P.Push(tg)

		// Catch P error and exit
		if i.Cmp(zero) > 0 && (tg.Cmp(zero) < 0 || tg.Cmp(g) >= 0) {
			break
		}

		// Q[n+1]=Q[n-1]+q[n](r[n]-r[n-1])
		Qp := new(big.Int)
		Qp = Qp.Sub(rn, rm).Mul(Qp, qp).Add(Qp, Qm)
		Q.Push(Qp)

		// A^2_[n-1] (mod N) = (-1)^nQ[n] (mod N)
		sign := new(big.Int)
		sign = sign.Exp(mone, i, nil).Mul(sign, Qn)
		smod := new(big.Int)
		smod = smod.Mod(sign, N)
		ASq := new(big.Int)
		ASq = ASq.Mul(An, An).Mod(ASq, N)

		if (i.Cmp(zero) > 0 && smod.Cmp(ASq) != 0) || sign.Cmp(one) <= 0 {
			continue
		}

		// If Q is a perfect square check if A-Q pair leads
		// to a factor of N
		square, QSq := PerfectSquare(sign)
		if !square {
			continue
		}

		AmSqQ := new(big.Int)
		AmSqQ = AmSqQ.Sub(An, QSq)

		// If D=GCD(A-sqrt(Q), N) and
		// 1 < D < N then D is a non trivial factor
		D = D.GCD(nil, nil, AmSqQ, Num)
		if D.Cmp(one) > 0 && D.Cmp(Num) < 0 {
			B = B.Div(Num, D)
			break
		}
	}

	if B.Cmp(zero) != 0 && D.Cmp(zero) != 0 {
		// Recursion step to find all factors
		FD = CFRAC(D, k)
		FB = CFRAC(B, k)
	} else {
		// If no factors found use Shanks method to factorise N
		FD = SQUFOF(N)
	}

	factors = append(factors, FD...)
	factors = append(factors, FB...)

	return factors
}

var multipliers []*big.Int = []*big.Int{
	big.NewInt(int64(1)),
	big.NewInt(int64(3)),
	big.NewInt(int64(5)),
	big.NewInt(int64(7)),
	big.NewInt(int64(11)),
	big.NewInt(int64(15)),
	big.NewInt(int64(21)),
	big.NewInt(int64(33)),
	big.NewInt(int64(35)),
	big.NewInt(int64(55)),
	big.NewInt(int64(77)),
	big.NewInt(int64(105)),
	big.NewInt(int64(165)),
	big.NewInt(int64(231)),
	big.NewInt(int64(385)),
	big.NewInt(int64(1155)),
}

/*
Square Forms Factorisation technique (Shanks)
*/
func SQUFOF(N *big.Int) []*big.Int {
	if prime := N.ProbablyPrime(10); prime {
		return []*big.Int{N}
	}

	s := new(big.Int)
	s = s.Sqrt(N)

	square, sq := PerfectSquare(N)

	// Perfect square
	if square {
		o := new(big.Int)
		o = sq
		for square {
			square, o = PerfectSquare(sq)
			if square {
				sq = o
			}
		}
		factors := SQUFOF(sq)
		factors = append(factors, factors...)
		return factors
	}

	for i := 0; i < len(multipliers); i++ {
		k := multipliers[i]

		kN := new(big.Int)
		kN = kN.Mul(k, N)
		D := new(big.Int)
		D = D.Sqrt(kN)

		// Initialise rings for P Q and b values
		var P, Q, b Ring

		// P[0] = floor(sqrt(kN))
		Po := new(big.Int)
		Po = D
		P.Push(D)
		PoPo := new(big.Int)
		PoPo = PoPo.Mul(Po, Po)

		// Q[0] = 1, Q[1] = kN-P[0]^2
		Q.Push(one)
		Nms := new(big.Int)
		Nms = Nms.Sub(kN, PoPo)
		Q.Push(Nms)

		// Set up limits
		B := new(big.Int)
		B = B.Mul(D, two).Sqrt(B).Mul(B, two)
		n := new(big.Int)

		for n = big.NewInt(int64(1)); n.Cmp(B) <= 0; n.Add(n, one) {
			_, Pn := P.Get()
			Qm, Qn := Q.Get()

			// b[i] = floor((P[0]+P[i-1])Q[i])
			bi := new(big.Int)
			bi = bi.Add(Po, Pn).Div(bi, Qn)
			b.Push(bi)

			// P[i] = b[i]Q[i]-P[i-1]
			Pi := new(big.Int)
			Pi = Pi.Mul(bi, Qn).Sub(Pi, Pn)
			P.Push(Pi)

			// Q[i+1]=Q[i-1]+b[i](P[i-1]-P[i])
			bPP := new(big.Int)
			bPP = bPP.Sub(Pn, Pi).Mul(bPP, bi)
			Qi := new(big.Int)
			Qi = Qi.Add(Qm, bPP)
			Q.Push(Qi)

			//fmt.Printf("i: %v, P: %v, Q: %v, b: %v\n", n, Pi, Qn, bi)
			// When 2|(i+1) and Q[i+1] is a perfect square exit
			square, _ := PerfectSquare(Qi)
			even := new(big.Int)
			even = even.Add(n, one).Mod(even, two)
			if even.Cmp(zero) == 0 && square {
				break
			}
		}

		_, Pn := P.Get()
		_, Qn := Q.Get()
		square, Qs := PerfectSquare(Qn)
		if n.Cmp(B) >= 0 && !square {
			continue
		}

		// b[0] = floor((P[0]-P[i])/sqrt(Q[i+1]))
		bo := new(big.Int)
		bo = bo.Sub(Po, Pn).Div(bo, Qs)
		b.Push(bo)

		// P[0]=b[0]sqrt(Q[i+1])+P[i]
		Pp := new(big.Int)
		Pp = Pp.Mul(bo, Qs).Add(Pp, Pn)
		P.Push(Pp)

		// Q[0] = sqrt(Q[i+1]), Q[1] = (kN-P[0]^2)/Q[0]
		Pp2 := new(big.Int)
		Pp2 = Pp2.Mul(Pp, Pp)
		NmP := new(big.Int)
		NmP = NmP.Sub(kN, Pp2).Div(NmP, Qs)
		Q.Push(Qs)
		Q.Push(NmP)

		for {
			_, Pn := P.Get()
			Qm, Qn := Q.Get()

			// b[i] = floor((P[0]+P[i-1])Q[i])
			bi := new(big.Int)
			bi = bi.Add(Po, Pn).Div(bi, Qn)
			b.Push(bi)

			// P[i] = b[i]Q[i]-P[i-1]
			Pi := new(big.Int)
			Pi = Pi.Mul(bi, Qn).Sub(Pi, Pn)
			P.Push(Pi)

			// Q[i+1]=Q[i-1]+b[i](P[i-1]-P[i])
			bPP := new(big.Int)
			bPP = bPP.Sub(Pn, Pi).Mul(bPP, bi)
			Qi := new(big.Int)
			Qi = Qi.Add(Qm, bPP)
			Q.Push(Qi)

			//fmt.Printf("P: %v, Q: %v, b: %v\n", Pi, Qn, bi)
			// When P[i] = P[i-1] exit
			if Pn.Cmp(Pi) == 0 {
				break
			}
		}

		// f = GCD(N, P[i]) if f != 1 or f != N then f is a non
		// trivial factor of N
		_, Pn = P.Get()
		fac := new(big.Int)
		fac = fac.GCD(nil, nil, N, Pn)
		if fac.Cmp(one) == 0 || fac.Cmp(N) == 0 {
			continue
		}
		fac2 := new(big.Int)
		fac2 = fac2.Div(N, fac)

		FA := SQUFOF(fac)
		FB := SQUFOF(fac2)

		factors := []*big.Int{}
		factors = append(factors, FA...)
		factors = append(factors, FB...)

		return factors
	}

	return []*big.Int{N}
}

/*
Check if N is a perfect square, and return true, sqrt(N)
else return false, 0
*/
func PerfectSquare(N *big.Int) (bool, *big.Int) {
	if N.Cmp(zero) < 0 {
		return false, zero
	}
	s := new(big.Int)
	s = s.Sqrt(N)
	ss := new(big.Int)
	ss = ss.Mul(s, s)

	if ss.Cmp(N) == 0 {
		return true, s
	}

	return false, zero
}

/*
Quick sort given slice by first finding a correctly placed
pivot point and then recursively sorting the left and right
subarrays until all subarrays are of length 1 or 0, and then
combining the sorted subarrays into a single slice
*/
func Quicksort(array []*big.Int) []*big.Int {
	if len(array) <= 1 {
		return array
	}
	p1 := 0
	p2 := len(array) - 1
	// Correctly place the pivot element
	for p1 != p2 {
		if (array[p1].Cmp(array[p2]) > 0 && p1 < p2) || (array[p1].Cmp(array[p2]) < 0 && p1 > p2) {
			array[p1], array[p2] = array[p2], array[p1]
			p1, p2 = p2, p1
		}
		if p1 < p2 {
			p1 += 1
		} else {
			p1 -= 1
		}
	}

	// Sort left and right subarrays
	left := Quicksort(array[:p1])
	right := Quicksort(array[p1+1:])

	// Combine into new sorted slice
	final := []*big.Int{}
	final = append(final, left...)
	final = append(final, array[p1])
	final = append(final, right...)

	return final
}