diff --git a/cmd/convergents.go b/cmd/convergents.go index 73ed939..e5b9bda 100644 --- a/cmd/convergents.go +++ b/cmd/convergents.go @@ -21,67 +21,9 @@ import ( "math/big" "github.com/spf13/cobra" + "scm.dairydemon.net/filifa/mathtools/internal/lib" ) -func cycle(seq []*big.Int) <-chan *big.Int { - ch := make(chan *big.Int) - n := len(seq) - go func() { - for i := 0; true; i = (i + 1) % n { - ch <- seq[i] - } - }() - - return ch -} - -func gaussianBrackets(ch <-chan *big.Int) <-chan *big.Int { - out := make(chan *big.Int) - - xprev := big.NewInt(0) - x := big.NewInt(1) - - go func() { - tmp := new(big.Int) - for a := range ch { - out <- x - - tmp.Mul(a, x) - tmp.Add(tmp, xprev) - - xprev.Set(x) - x = new(big.Int).Set(tmp) - } - }() - - return out -} - -func cFracConvergents(a0 *big.Int, denoms []*big.Int) <-chan *big.Rat { - hc := cycle(denoms) - _ = <-hc - hch := gaussianBrackets(hc) - - kc := cycle(denoms) - kch := gaussianBrackets(kc) - _ = <-kch - - a := new(big.Rat).SetInt(a0) - - out := make(chan *big.Rat) - - go func() { - for { - h, k := <-hch, <-kch - r := new(big.Rat).SetFrac(h, k) - r.Add(r, a) - out <- r - } - }() - - return out -} - func convergents(cmd *cobra.Command, args []string) { a0, ok := new(big.Int).SetString(args[0], 10) if !ok { @@ -97,7 +39,7 @@ func convergents(cmd *cobra.Command, args []string) { } } - for r := range cFracConvergents(a0, denoms) { + for r := range lib.CFracConvergents(a0, denoms) { fmt.Println(r) } } diff --git a/cmd/pell.go b/cmd/pell.go index 6251da6..9124071 100644 --- a/cmd/pell.go +++ b/cmd/pell.go @@ -47,8 +47,9 @@ func pell(cmd *cobra.Command, args []string) { period = 2 * r } - ch := cFracConvergents(a0, repetend) + ch := lib.CFracConvergents(a0, repetend) + // TODO: consider breaking after finding the fundamental solution and using the recurrence relation to find the following solutions for i := 1; true; i = (i + 1) % period { r := <-ch if i == period-1 { diff --git a/internal/lib/continuedFrac.go b/internal/lib/continuedFrac.go index 66d2eb6..1f96e8d 100644 --- a/internal/lib/continuedFrac.go +++ b/internal/lib/continuedFrac.go @@ -36,3 +36,62 @@ func SqrtRepetend(x *big.Int) ([]*big.Int, error) { return repetend, nil } + +func cycle(seq []*big.Int) <-chan *big.Int { + ch := make(chan *big.Int) + n := len(seq) + go func() { + for i := 0; true; i = (i + 1) % n { + ch <- seq[i] + } + }() + + return ch +} + +func GaussianBrackets(ch <-chan *big.Int) <-chan *big.Int { + out := make(chan *big.Int) + + xprev := big.NewInt(0) + x := big.NewInt(1) + + go func() { + tmp := new(big.Int) + for a := range ch { + out <- x + + tmp.Mul(a, x) + tmp.Add(tmp, xprev) + + xprev.Set(x) + x = new(big.Int).Set(tmp) + } + }() + + return out +} + +func CFracConvergents(a0 *big.Int, denoms []*big.Int) <-chan *big.Rat { + hc := cycle(denoms) + _ = <-hc + hch := GaussianBrackets(hc) + + kc := cycle(denoms) + kch := GaussianBrackets(kc) + _ = <-kch + + a := new(big.Rat).SetInt(a0) + + out := make(chan *big.Rat) + + go func() { + for { + h, k := <-hch, <-kch + r := new(big.Rat).SetFrac(h, k) + r.Add(r, a) + out <- r + } + }() + + return out +}