/*
a254077.go calculates sequence oeis.org/A254077

Flag arguments likely to be useful include:

    -num=<integer>   [20]
	Attempt to generate this many terms.
	The default of 20 is useful for simple testing.

    -prog=<integer>   [0]
	If positive, show progress on standard err.
	Useful for comparing performance of different implementations.
	Progress will be reported when the number of terms is a multiple.
	For example, if -prog=100000 is specified,
	progress is reported at every multiple of 100000

	A typical progress report might look like

	    a[200000]=204569: 2.849688009s 1m35.221598503s

	    showing
		the sequence value for the 200000th term
		the elapsed time to compute the last 100000 terms
		the total elapsed time computing terms
		[Earlier versions of this program total elapsed time
		 since the job began, rather than time computing terms.
		 This conflated initialization with time computing terms.]

    -headroom=<integer>   [5]
        A sieve, Cbits, is used to record which composite integers
	have not yet appeared in the sequence.  There is no known
	way of predicting how large the sieve must be to assure
	that all composite candidates can be interrogated.
	For now we use the number of terms requested
	(see the -num flag) times this overhead factor.
	It is adequate for 2.5 billion terms.
	Overestimating it may exhaust main memory.
	Underestimating it will cause the program to terminate
	with an array bounds error.

    -seento=<integer>   [2003]
        For primes less than or equal to this, keep track of the
	least unseen multiple of that prime.  Important for reducing
	repeated skipping over smaller multiples.

    -gcd1max=<integer>   [unlimited]
	Insist that gcd(a[n],a[n-1]) be no greater than this.

    -time=<bool>   [true]
	Report major times and some event counts on standard err.
	For example, for the default of computing 20 elements,
	you might see

	    Making sieves for 51103 elements: 604.755us
	    Finding 20 elements: 75.568us
	    Total run time: 755.822us
	    42 gcd() calls, 17 wrong order, 103 loops, 0 back2back

	The event counts on the last line are probably of limited interest
	unless you are modifying the program, but do no harm.

     -study=<bool>   [false]
 	Report considerable detail about the generation of each term.
	This will slow things down a little, and expand the output a lot.
	Some typical output on the way to 10 million terms for an
	earlier version of this program included

        18 20:2^2*5:4>1:1:3^3:11:20:4:1:1:0
	19 33:3*11:3>1:1:2^2*5:11:22:4:2:8:6
	20 24:2^3*3:4>3:3:11p:11:22:4:1:2:1
	21 11:11P:11>1:1:2^3*3:11:22:5:1:0:0
	22 26:2*13:2>1:1:11p:13:22:5:1:3:2
	23 22:2*11:11>2:2:13p:0:22:5:0:1:0
	    ... following split after field 4 to avoid very long lines
	9999571 10175857:1213*8389:1213>1:1:5*13^2*12041:
	    5087323:10171703:354156:1:1079:3076
	9999572 10174655:5*71*28661:5>1:1:1213*8389:
	    5087323:10171703:354156:1:51:2902
	9999573 10177070:2*5*839*1213:1213>5:5:71*28661:
	    5087323:10171703:354156:2:2178:3190

        10 colon-delimited fields follow n and a[n] on each line.

	    For brevity, we will refer to
		gcd2 = gcd(a[n], a[n-2])
		gcd1 = gcd(a[n], a[n-1])

	The fields are:

	1) The factors of a[n].  They are followed by a "P" iff a[n] is prime.
	2) gcd2 > gcd1
	3) The factors of gcd1.
	4) The factors of a[n-1]/gcd1.  They are followed by a "p" iff prime.
	5) The least prime not yet seen
	6) The least composite not yet seen
	7) The total number of primes seen so far
	8) The number of primes considered as candidates for a[n]
	9) The number of composites considered as candidates for a[n]
	10) The number of composite candidates skipped.
	    A candidate may be skipped because it is, in fact, a prime,
	    or because it is a composite that has already appeared.

	Some observations about the terms above.

	For n=21, fields 9 and 10 are both 0.  A prime candidate
	less than the least composite candidate proved satisfactory.

	For n=23, fields 5 and 8 are both 0.  Because a[n-2] is prime,
	no other prime can have gcd2 other than 1, so no primes are considered.

	For n=9999573, fields 9 and 10 are quite large.  That is because
	a[n-2] is the product of two large primes.  Composite candidates
	in the version of the program that produced the output were only
	searched sequentially.  We encountered many (3190) candidates
	that were either prime or already seen, and many (2178) that
	were unseen composites only the last of which was satisfactory.
	This observation led to changes in the way composite candidates
	were considered.
*/

package main

import (
	"flag"
	"fmt"
	// "github.com/davecheney/profile"
	"math"
	"os"
	"time"
)

var Huge int64 = (1 << 63) - 1 // Big as a 64-bit int can hold

var num = flag.Int64("num", 20, "Number of elements to report")
var growth = flag.Int64("grow", 1000, "Growth when no candidates remain")
var showtimes = flag.Bool("time", true, "Show processing times on standard err")
var progress = flag.Int64("prog", 0, "Show progress at multiples of this")
var gcd1max = flag.Int64("gcd1max", Huge, "Maximum value of gcd1")
var headroom = flag.Int64("headroom", 5, "Control size of composites sieve")
var study = flag.Bool("study", false, "Show considerable detail about terms")
var fseen = flag.Int64("seento", 2003, "Reduce skipping for primes through this")

/*
Program history (in case someone plans to modify it):

Originally, we kept just a single "candidate pool" in an Unused structure.
It was observed that the pool tended to fill up with unseen primes.
When the proper candidate was a composite (which is most of the time),
we had to step through all the primes before getting to the first
composite candidate.

The program was changed to maintain separate pools of prime and
composite candidates, but still in Unused structures.  This is a
very satisfactory mechanism for maintaining prime candidates,
since it has been observed (but not yet proved) that the primes
are added to the sequence in order of their appearance as integers.
The candidates are removed in order until the structure is empty,
then replenished.

This is *not* a satisfactory mechanism for maintaining composite
candidates!  When a[n-2] is prime, only composite multiples of that prime
can appear as a[n].  As larger and larger primes appear in the sequence,
stepping through the composite candidates sequentially, the only method
available with the Unused structure, is inordinately expensive.
For prime p on the order of a million or more, we visit almost that
many composite candidates before encountering a multiple of p.
The composite pool grows for every *growth (1000) new candidates added,
probably forcing the pool to be reallocated and copied.
With 2p usually (but not yet provably) already seen, enabling
the appearance of p two steps later in the squence, we don't get
a serious composite until we reach 3p.  So there was a lot of
"flailing" whenever we encounter a prime, something that was visible
watching the sequence grow, "stalling" whenever a prime was the
penultimate term.

So the program was modified to maintain composite candidates as a
"sieve", initially identical to the sieve used to identify the primes.
Although that allowed us to directly investigate whether any given
integer had already been seen, sequential access to find the next
unseen integer was rather slow, as discussed in the comment above
about the study flag.  So we "flipped the emphasis".  Instead of
walking through unseen candidates, looking for those that were not
coprime with a[n-2], we instead used the factors of a[n-2] to
to generate, in order, integers that could not be coprime with a[n-2],
and used the sieve to determine if they were unseen.
*/

// Keep track of unused integers
type Unused struct {
	next   int64            // never been seen at all
	seen   []int64          // seen but not used
	last   int64            // where we left off in seen
	zthru  int64            // how many leading 0 elements in seen
	grow   int64            // how many new terms to add when filling
	filter func(int64) bool // filter for new terms if non-nil
}

func fillup(u *Unused) {
	if u.filter == nil {
		for len(u.seen) < cap(u.seen) {
			u.seen = append(u.seen, u.next)
			u.next++
		}
		return
	}
	nxt := u.next
	for ; ; nxt++ {
		if u.filter(nxt) {
			u.seen = append(u.seen, nxt)
			if len(u.seen) == cap(u.seen) {
				break
			}
		}
	}
	u.next = nxt + 1
}

func NewU(startAt, growBy int64, filter func(int64) bool) *Unused {
	u := new(Unused)
	u.next = startAt
	u.grow = growBy
	u.filter = filter
	u.seen = make([]int64, 0, u.grow)
	u.last = -1
	u.zthru = 0
	fillup(u)
	return u
}

func squeeze(from, to *[]int64) {
	for _, i := range *from {
		if i != 0 {
			*to = append(*to, i)
		}
	}
}

func grow(u *Unused) {
	var size int64 = int64(cap(u.seen)) + *growth
	bigger := make([]int64, 0, size)
	squeeze(&u.seen, &bigger)
	u.zthru = 0
	u.seen = bigger
	u.last = int64(len(u.seen)) - 1
	fillup(u)
}

func reset(u *Unused) {
	u.last = u.zthru - 1
}

func saw(u *Unused) {
	u.seen[u.last] = 0
}

func next(u *Unused) int64 {
	for {
		possible := u.seen[u.last+1 : cap(u.seen)]
		var i int64 = 0
		for _, n := range possible {
			if n != 0 {
				if u.last+1 == u.zthru {
					u.zthru += i
				}
				u.last += i + 1
				return n
			}
			i++
		}
		grow(u)
	}
}

var Primes []int64 // enough primes to factor anything we might see
var Cbits []int64  // composite candidates packed into 64-bit words

func factor(n int64) [][2]int64 {
	factors := make([][2]int64, 0, 10)
	for _, p := range Primes {
		if p*p > n {
			break
		}
		if n%p == 0 {
			m := int64(1)
			n /= p
			for n%p == 0 {
				n /= p
				m++
			}
			factors = append(factors, [2]int64{p, m})
		}
	}
	if n != 1 {
		factors = append(factors, [2]int64{n, 1})
	}
	return factors
}

func isPrime(n int64) bool {
	if (n & 1) == 0 {
		return n == 2
	}
	for _, p := range Primes[1:] {
		if p*p > n {
			return true
		}
		if n%p == 0 {
			return false
		}
	}
	return true // should not happen
}

func isPrimeFact(factors [][2]int64) bool {
	return len(factors) == 1 && factors[0][1] == 1
}

func wordBit(index int64) (word, bit int64) {
	bit = 1 << uint(index&63)
	word = index >> 6
	return
}

func isCand(index int64) bool {
	word, bit := wordBit(index)
	return Cbits[word]&bit != 0
}

func makeCand(index int64) {
	word, bit := wordBit(index)
	Cbits[word] |= bit
}

func clearCand(index int64) {
	word, bit := wordBit(index)
	Cbits[word] &= ^bit
}

func MakeNonPrimeSieve(n int64) {
	// A bit b in Cbits should be true iff b is a composite,
	// and not already seen
	Cbits = make([]int64, 1+(n>>6))
	Primes = make([]int64, 0, int64(math.Log(float64(n))))
	n = 64 * int64(len(Cbits))
	var x int64 = 2

	for ; x < n; x++ {
		if !isCand(x) {
			Primes = append(Primes, x)
			for y := 2 * x; y < n; y += x {
				makeCand(y)
			}
		}
	}
}

func fS(factors [][2]int64) string {
	if len(factors) < 1 {
		return "1"
	}
	fs := ""
	sep := ""
	for _, pair := range factors {
		fs += sep
		sep = "*"
		fs += fmt.Sprintf("%d", pair[0])
		if pair[1] > 1 {
			fs += fmt.Sprintf("^%d", pair[1])
		}
	}
	return fs
}

var GcdCalls int64 = 0
var GcdDisorder int64 = 0
var GcdLoops int64 = 0

// Most of the variables that follow were originally local to main().
// Global variables are ugly, but, so too, is repeated code.
// Refactoring the repeated code into subroutines was easier
// if variables that were needed and potentially modified were global.
var n int64         // (0-based) index of next element in sequence
var Cbase int64 = 4 // First index in Cbits that hasn't been seen already
var Cnext int64     // Next index in Cbits that hasn't been seen or rejected
// Info about search for next element, not total
var Cskips int64        // Composite candidates skipped because seen
var Ccands int64        // Composite candidates considered
var Pcands int64        // Prime candidates considered
var back2, back1 int64  // a[n-2] and a[n-1]
var pback2, pback1 bool // are back2/back1 prime
var gcd2, gcd1 int64    // gcd(k, back2) and gcd(k, back1)
var u2 int64            // back2 / gcd1
var gcd1not1 bool       // was last gcd1 != 1,
var pu2 bool            // is u2 prime
var pthis bool          // is a[n] prime
var nxtP, nxtC int64    // next prime and composite candidates for a[n]
var back2factors, back1factors, u2factors [][2]int64
var CstepSize int64
var seenTo int64
var CSeenTo []int64

/*
It has been conjectured, but not (yet) proved, that
    if prime p has not appeared and a[n-2] is k*p
    then p does not divide a[n-1], so a[n] will be p.
This has been confirmed through 2,500,000,000 terms.
We count all instances where it fails to be true.
*/
var Back2backP = 0

func gcd(a, b int64) int64 {
	GcdCalls++
	if a < b {
		GcdDisorder++
		a, b = b, a
	}
	for b != 0 {
		GcdLoops++
		a, b = b, a%b
	}
	return a
}

type advancer func() bool

var nextCfunc advancer

func fixedInit() {
	prime := u2factors[0][0]
	CstepSize = prime
	// The only possible gcds are prime**i, i >= 1
	// If gcd1 also has prime**k as a factor,
	// then we must *beat* k as a factor.
	baseline := nxtC
	if prime < seenTo && baseline < CSeenTo[prime] {
		baseline = CSeenTo[prime]
	}
	common := gcd1
	for common%prime == 0 {
		common /= prime
		CstepSize *= prime
	}
	nxtC = (baseline + CstepSize - 1) / CstepSize
	nxtC *= CstepSize
	for !isCand(nxtC) {
		nxtC += CstepSize
		Cskips++
	}
}

func nextCfixed() bool {
	// nxtC will always be a composite multiple of CstepSize.
	// On entry, nxtC failed to satisfy the gcd relationships.
	// On exit, nxtC will be the least unseen composite multiple.
	// On true exit, nxtC is less than nxtP.
	// On false exit, nxtC is greater than nxtP.
	nxtC += CstepSize
	for !isCand(nxtC) {
		nxtC += CstepSize
		Cskips++
	}
	return nxtC < nxtP
}

var FactorList [][2]int64

func fancyInit() {
	factors := u2factors
	nf := len(factors)
	FactorList = make([][2]int64, nf)
	var first int64
	for i := 0; i < nf; i++ {
		stepSize := factors[i][0] // prime factor
		FactorList[i][1] = stepSize
		if stepSize < int64(len(CSeenTo)) {
			first = CSeenTo[stepSize]
		} else {
			first = (nxtC + stepSize - 1) / stepSize
			first *= stepSize
		}
		for !isCand(first) {
			first += stepSize
			Cskips++
		}
		if stepSize < int64(len(CSeenTo)) {
			CSeenTo[stepSize] = first
		}
		FactorList[i][0] = first
	}
	nextCfancy() // set nxtC to the least multiple
	return
}

func nextCfancy() bool {
	least := Huge
	leastAt := -1
	for i, NxtStepSize := range FactorList {
		nxt, stepSize := NxtStepSize[0], NxtStepSize[1]
		if least > nxt {
			least = nxt
			leastAt = i
		} else if least == nxt {
			nxt += stepSize
			for !isCand(nxt) {
				nxt += stepSize
				Cskips++
			}
			FactorList[i][0] = nxt
		}
	}
	NxtStepSize := FactorList[leastAt]
	nxt, stepSize := NxtStepSize[0], NxtStepSize[1]
	nxt += stepSize
	for !isCand(nxt) {
		nxt += stepSize
		Cskips++
	}
	FactorList[leastAt][0] = nxt
	nxtC = least
	return nxtC < nxtP
}

func checkC() bool {
	// On entry, nxtC is the least composite that has not been seen.
	// On true exit, nxtC is a qualifying composite less than nxtP.
	// On false exit, nxtC is an unseen composite greater than nxtP.
	for cLow := true; cLow; cLow = nextCfunc() {
		Ccands++
		gcd2 = gcd(back2, nxtC)
		if gcd2 > 1 {
			gcd1 = gcd(back1, nxtC)
			if gcd2 > gcd1 && gcd1 <= *gcd1max {
				return true
			}
		}
	}
	return false
}

func main() {
	// defer profile.Start(profile.CPUProfile).Stop()

	flag.Parse()
	start := time.Now()
	progFirst := start
	progLast := start
	then := start

	var nxt, c1, p1 int64

	np := *num
	if np < *growth {
		np = *growth
	}
	np *= *headroom
	if int64(int(np)) != np {
		fmt.Fprintf(os.Stderr,
			"Warning! int != int64, results not reliable\n")
	}
	MakeNonPrimeSieve(np)
	seenTo = *fseen
	if seenTo < 19 {
		seenTo = 19 // Don't be silly
	}
	if seenTo > Primes[len(Primes)-1] {
		seenTo = Primes[len(Primes)-1]
	}
	seenTo++ // Allow direct indexing by prime number
	CSeenTo = make([]int64, seenTo)
	for _, p := range Primes {
		if p < seenTo {
			CSeenTo[p] = 2 * p // p is not composite
		} else {
			break
		}
	}
	if *showtimes {
		fmt.Fprintf(os.Stderr,
			"Making sieves for %d elements: %s\n",
			np, time.Since(then))
		then = time.Now()
	}
	progAt := *num + 1
	if *progress > 0 {
		progAt = *progress
		if *progress < 4 {
			progAt = (4 + *progress - 1) / *progress
			progAt *= *progress
		}
		progFirst = time.Now()
		progLast = progFirst
	}

	// report first 3 elements
	for n = 0; n < 3; n++ {
		fmt.Fprintf(os.Stdout, "%d %d\n", n+1, n+1)
	}
	primesSeen := 2 // 2 and 3
	pNums := NewU(n+1, *growth,
		func(n int64) bool {
			return isPrime(n)
		})

	Cbase = n + 1
	back2, pback2, back1, pback1 = 2, true, 3, true
	gcd1not1, gcd1, u2, pu2 = true, 1, back2, true
	back2factors, back1factors, u2factors =
		factor(back2), factor(back1), factor(u2)

	for n < *num {
		Pcands, Ccands, Cskips = 0, 0, 0
		for !isCand(Cbase) {
			Cbase++
			Cskips++
		}
		Cnext = Cbase
		nxtC = Cnext
		c1 = nxtC
		// Consider next prime candidate
		if pback2 {
			// All other primes will have gcd == 1
			nxtP = Huge
			p1 = 0
		} else {
			// It's possible the next term is a prime
			reset(pNums)
			nxtP = next(pNums)
			p1 = nxtP
		}
		var maxFact int64 = back2factors[len(back2factors)-1][0]
		if len(u2factors) == 1 {
			// Single factor
			fixedInit()
			nextCfunc = nextCfixed
		} else {
			// True composite
			fancyInit()
			nextCfunc = nextCfancy
		}
		for {
			for ; nxtP < nxtC; nxtP = next(pNums) {
				if nxtP > maxFact {
					nxtP = Huge
					break
				}
				Pcands++
				if (back2 % nxtP) != 0 {
					// gcd is 1, cannot qualify
					continue
				}
				gcd2 = nxtP
				if (back1 % nxtP) == 0 {
					// both have nxtP as a common factor
					// gcd2 and gcd1 are the same
					// so nxtP is rejected
					Back2backP++
					continue
				}
				gcd1 = 1
				nxt = nxtP
				saw(pNums)
				pthis = true
				goto report
			}
			for nxtC < nxtP {
				if checkC() {
					goto reportC
				}
			}
		}

	reportC: // Report a composite next term
		nxt = nxtC
		clearCand(nxt)
		pthis = false
	report:
		n++
		if n == progAt {
			progAt += *progress
			fmt.Fprintf(os.Stderr,
				"a[%d]=%d: %s %s\n",
				n, nxt,
				time.Since(progLast), time.Since(progFirst))
			progLast = time.Now()

		}
		gcd1not1 = (gcd1 != 1)
		if gcd1not1 {
			u2 = back1 / gcd1
			u2factors = factor(u2)
			pu2 = isPrimeFact(u2factors)
		} else {
			u2 = back1
			pu2 = pback1
			u2factors = back1factors
		}
		if pthis {
			primesSeen++
		}
		thisfactors := factor(nxt)
		if *study {
			Pthis, Pu2 := "", ""
			if pthis {
				Pthis = "P"
			}
			if pu2 {
				Pu2 = "p"
			}
			fmt.Fprintf(os.Stdout,
				"%d %d:%s%s:%d>%d:%s:%s%s:%d:%d:%d:%d:%d:%d\n",
				n, nxt,
				fS(thisfactors), Pthis, gcd2, gcd1,
				fS(factor(gcd1)), fS(u2factors),
				Pu2, p1, c1,
				primesSeen, Pcands, Ccands, Cskips)
		} else {
			fmt.Fprintf(os.Stdout, "%d %d\n", n, nxt)
		}
		back2, pback2, back1, pback1 =
			back1, pback1, nxt, pthis
		back2factors, back1factors = back1factors, thisfactors
	}
	if *showtimes {
		fmt.Fprintf(os.Stderr, "Finding %d elements: %s\n",
			*num, time.Since(progFirst))
		fmt.Fprintf(os.Stderr, "Total run time: %s\n",
			time.Since(start))
	}
	fmt.Fprintf(os.Stderr,
		"%d gcd() calls, %d wrong order, %d loops, %d back2back\n",
		GcdCalls, GcdDisorder, GcdLoops, Back2backP)
}