filifa
/
mathtools
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

move to lib

This commit is contained in:
filifa 2025-09-29 21:52:52 -04:00
parent d2275cd69f
commit b55fe61820
16 changed files with 452 additions and 296 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
cmd/discreteLog.go
View File

@ -17,11 +17,11 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
package cmd
import (
"errors"
"fmt"
"math/big"
"github.com/spf13/cobra"
"scm.dairydemon.net/filifa/mathtools/internal/lib"
)
var discreteLogModulus string
@ -29,62 +29,6 @@ var discreteLogBase string
var discreteLogElement string
var discreteLogOrder string
// whyyyy doesn't math/big have a ceil functionnnn
func ceilSqrt(x *big.Int) *big.Int {
z := new(big.Int).Sqrt(x)
s := new(big.Int).Exp(z, big.NewInt(2), nil)
if s.Cmp(x) != 0 {
z.Add(z, big.NewInt(1))
}
return z
}
// TODO: this can be extended to work with n, b not coprime
// https://cp-algorithms.com/algebra/discrete-log.html
func babyStepGiantStep(n, b, x, order *big.Int) (*big.Int, error) {
z := new(big.Int).GCD(nil, nil, b, n)
if z.Cmp(big.NewInt(1)) != 0 {
return nil, fmt.Errorf("base %v and modulus %v are not coprime", b, n)
}
var m *big.Int
if order == nil {
// m = ceil(sqrt(n - 1))
z := big.NewInt(1)
z.Sub(n, z)
m = ceilSqrt(z)
} else {
m = ceilSqrt(order)
}
table := make(map[string]*big.Int)
for j := big.NewInt(1); j.Cmp(m) <= 0; j.Add(j, big.NewInt(1)) {
a := new(big.Int).Exp(b, j, n)
table[a.String()] = new(big.Int).Set(j)
}
// p = b^-m modulo n
p := new(big.Int).Neg(m)
p.Exp(b, p, n)
gamma := new(big.Int).Set(x)
for i := big.NewInt(0); i.Cmp(m) == -1; i.Add(i, big.NewInt(1)) {
j, ok := table[gamma.String()]
if ok {
i.Mul(i, m)
i.Add(i, j)
return i, nil
}
gamma.Mul(gamma, p)
gamma.Mod(gamma, n)
}
return nil, errors.New("no solution")
}
func discreteLog(cmd *cobra.Command, args []string) {
m, ok := new(big.Int).SetString(discreteLogModulus, 10)
if !ok {
@ -109,7 +53,7 @@ func discreteLog(cmd *cobra.Command, args []string) {
}
}
k, err := babyStepGiantStep(m, b, x, order)
k, err := lib.BabyStepGiantStep(m, b, x, order)
if err != nil {
cobra.CheckErr(err)
}

22
cmd/divisor.go
View File

@ -21,27 +21,9 @@ import (
"math/big"
"github.com/spf13/cobra"
"scm.dairydemon.net/filifa/mathtools/internal/lib"
)
// TODO: expand to work for different sigmas
func divisorSummatory(n *big.Int) *big.Int {
// employing Dirichlet's hyperbola method
sqrt := new(big.Int).Sqrt(n)
total := big.NewInt(0)
for x := big.NewInt(1); x.Cmp(sqrt) <= 0; x.Add(x, big.NewInt(1)) {
z := new(big.Int).Div(n, x)
total.Add(total, z)
}
total.Mul(total, big.NewInt(2))
sqrt.Exp(sqrt, big.NewInt(2), nil)
total.Sub(total, sqrt)
return total
}
func divisorSum(cmd *cobra.Command, args []string) {
for _, arg := range args {
n, ok := new(big.Int).SetString(arg, 10)
@ -50,7 +32,7 @@ func divisorSum(cmd *cobra.Command, args []string) {
continue
}
d := divisorSummatory(n)
d := lib.DivisorSummatory(n)
fmt.Println(d)
}
}

61
cmd/divisors.go
View File

@ -20,71 +20,14 @@ import (
"fmt"
"github.com/spf13/cobra"
"scm.dairydemon.net/filifa/mathtools/internal/lib"
)
var divisorsN uint
var divisorsE uint
func pow(base uint, exp uint) uint {
result := uint(1)
for exp > 0 {
if exp%2 == 1 {
result *= base
}
exp >>= 1
base *= base
}
return result
}
func updateMultiples(sieve []uint, x uint, p uint, n uint) {
for q := p; ; q *= p {
// sigma_x(a*b) = sigma_x(a) * sigma_x(b) if gcd(a,b) = 1
for i := 2 * q; i < n; i += q {
if i%(p*q) != 0 {
sieve[i] *= sieve[q]
}
}
if p*q >= n {
break
}
// sigma_x(p^k) = p^(kx) + sigma_x(p^(k-1))
sieve[p*q] = pow(p*q, x) + sieve[q]
}
}
func divisorsSieve(n uint, x uint) chan uint {
sieve := make([]uint, n)
sieve[0] = 0
for i := uint(1); i < n; i++ {
sieve[i] = 1
}
ch := make(chan uint)
go func() {
for i := uint(0); i < n; i++ {
if i == 0 || i == 1 || sieve[i] != 1 {
ch <- sieve[i]
continue
}
sieve[i] = pow(i, x) + 1
updateMultiples(sieve, x, i, n)
ch <- sieve[i]
}
close(ch)
}()
return ch
}
func divisors(cmd *cobra.Command, args []string) {
ch := divisorsSieve(divisorsN, divisorsE)
ch := lib.DivisorsSieve(divisorsN, divisorsE)
for i := 0; ; i++ {
v, ok := <-ch
if !ok {

34
cmd/mobius.go
View File

@ -20,43 +20,13 @@ import (
"fmt"
"github.com/spf13/cobra"
"scm.dairydemon.net/filifa/mathtools/internal/lib"
)
var mobiusN uint
func mobiusSieve(n uint) chan int {
sieve := make([]int, n)
for i := 0; i < int(n); i++ {
sieve[i] = i
}
ch := make(chan int)
go func() {
for i := 0; i < int(n); i++ {
if i == 0 || i == 1 || sieve[i] != i {
ch <- sieve[i]
continue
}
ch <- -1
for j := 2 * i; j < int(n); j += i {
if j%(i*i) == 0 {
sieve[j] = 0
}
sieve[j] /= -i
}
}
close(ch)
}()
return ch
}
func mobius(cmd *cobra.Command, args []string) {
ch := mobiusSieve(mobiusN)
ch := lib.MobiusSieve(mobiusN)
for i := 0; ; i++ {
v, ok := <-ch
if !ok {

35
cmd/partitions.go
View File

@ -18,46 +18,15 @@ package cmd
import (
"fmt"
"math/big"
"strconv"
"github.com/spf13/cobra"
"scm.dairydemon.net/filifa/mathtools/internal/lib"
)
var partitionsN string
var partitionsK string
// given a slice where vals[i] = p_{k-1}(i) for a given k, nextPartition updates the slice so vals[i] = p_k(i) using the property that p_k(n) = p_k(n-k) + p_{k-1}(n)
func nextPartition(k int, vals []*big.Int) {
for i := 1; i < len(vals); i++ {
if i-k >= 0 {
vals[i].Add(vals[i], vals[i-k])
}
}
}
func p(n, k int) *big.Int {
if n < 0 {
return big.NewInt(0)
}
if k > n {
k = n
}
vals := make([]*big.Int, n+1)
for i := range vals {
vals[i] = big.NewInt(0)
}
vals[0] = big.NewInt(1)
for i := 1; i <= k; i++ {
nextPartition(i, vals)
}
return vals[n]
}
func partitions(cmd *cobra.Command, args []string) {
n, err := strconv.Atoi(partitionsN)
if err != nil {
@ -83,7 +52,7 @@ func partitions(cmd *cobra.Command, args []string) {
cobra.CheckErr("k must be nonnegative")
}
z := p(n, k)
z := lib.Partitions(n, k)
fmt.Println(z)
}

47
cmd/primeOmega.go
View File

@ -20,57 +20,14 @@ import (
"fmt"
"github.com/spf13/cobra"
"scm.dairydemon.net/filifa/mathtools/internal/lib"
)
var primeOmegaN uint
var primeOmegaMul bool
func primeOmegaUpdateMultiples(sieve []uint, p uint, n uint, multiplicity bool) {
for q := p; ; q *= p {
// omega(a*b) = omega(a) + omega(b) if gcd(a,b) = 1
for i := 2 * q; i < n; i += q {
if i%(p*q) != 0 {
sieve[i] += sieve[q]
}
}
if p*q >= n {
break
}
if multiplicity {
sieve[p*q] = 1 + sieve[q]
}
}
}
func primeOmegaSieve(n uint) chan uint {
sieve := make([]uint, n)
for i := uint(0); i < n; i++ {
sieve[i] = 0
}
ch := make(chan uint)
go func() {
for i := uint(0); i < n; i++ {
if i == 0 || i == 1 || sieve[i] != 0 {
ch <- sieve[i]
continue
}
sieve[i] = 1
primeOmegaUpdateMultiples(sieve, i, n, primeOmegaMul)
ch <- sieve[i]
}
close(ch)
}()
return ch
}
func primeOmega(cmd *cobra.Command, args []string) {
ch := primeOmegaSieve(primeOmegaN)
ch := lib.PrimeOmegaSieve(primeOmegaN, primeOmegaMul)
for i := 0; ; i++ {
v, ok := <-ch
if !ok {

27
cmd/shoelace.go
View File

@ -24,15 +24,11 @@ import (
"strings"
"github.com/spf13/cobra"
"scm.dairydemon.net/filifa/mathtools/internal/lib"
)
var shoelaceFile string
type point struct {
x float64
y float64
}
func splitLines(file *os.File) []string {
slice := make([]string, 0)
scanner := bufio.NewScanner(file)
@ -43,8 +39,8 @@ func splitLines(file *os.File) []string {
return slice
}
func readFromFile(filepath string) ([]point, error) {
points := make([]point, 0)
func readFromFile(filepath string) ([]lib.Point, error) {
points := make([]lib.Point, 0)
file, err := os.Open(filepath)
if err != nil {
return points, err
@ -65,25 +61,14 @@ func readFromFile(filepath string) ([]point, error) {
return points, err
}
points = append(points, point{x, y})
points = append(points, lib.Point{x, y})
}
return points, nil
}
func area(points []point) float64 {
total := float64(0)
n := len(points)
for i, p := range points {
q := points[(i+1)%n]
total += (p.y + q.y) * (p.x - q.x)
}
return total / 2
}
func shoelace(cmd *cobra.Command, args []string) {
var coordinates []point
var coordinates []lib.Point
var err error
if shoelaceFile != "" {
coordinates, err = readFromFile(shoelaceFile)
@ -94,7 +79,7 @@ func shoelace(cmd *cobra.Command, args []string) {
cobra.CheckErr("filename required")
}
fmt.Println(area(coordinates))
fmt.Println(lib.Area(coordinates))
}
// shoelaceCmd represents the shoelace command

30
cmd/totient.go
View File

@ -20,39 +20,13 @@ import (
"fmt"
"github.com/spf13/cobra"
"scm.dairydemon.net/filifa/mathtools/internal/lib"
)
var totientN uint
func totientSieve(n uint) chan uint {
totients := make([]uint, n)
totients[0] = 0
totients[1] = 1
for i := uint(2); i < n; i++ {
totients[i] = i - 1
}
ch := make(chan uint)
go func() {
for i := uint(0); i < n; i++ {
ch <- totients[i]
if i == 0 || i == 1 || totients[i] != i-1 {
continue
}
for j := uint(2 * i); j < n; j += i {
totients[j] -= totients[j] / i
}
}
close(ch)
}()
return ch
}
func totient(cmd *cobra.Command, args []string) {
for v := range totientSieve(totientN) {
for v := range lib.TotientSieve(totientN) {
if v == 0 {
continue
}

79
internal/lib/discreteLog.go Normal file
View File

@ -0,0 +1,79 @@
/*
Copyright © 2025 filifa
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package lib
import (
"errors"
"fmt"
"math/big"
)
// whyyyy doesn't math/big have a ceil functionnnn
func ceilSqrt(x *big.Int) *big.Int {
z := new(big.Int).Sqrt(x)
s := new(big.Int).Exp(z, big.NewInt(2), nil)
if s.Cmp(x) != 0 {
z.Add(z, big.NewInt(1))
}
return z
}
// TODO: this can be extended to work with n, b not coprime
// https://cp-algorithms.com/algebra/discrete-log.html
func BabyStepGiantStep(n, b, x, order *big.Int) (*big.Int, error) {
z := new(big.Int).GCD(nil, nil, b, n)
if z.Cmp(big.NewInt(1)) != 0 {
return nil, fmt.Errorf("base %v and modulus %v are not coprime", b, n)
}
var m *big.Int
if order == nil {
// m = ceil(sqrt(n - 1))
z := big.NewInt(1)
z.Sub(n, z)
m = ceilSqrt(z)
} else {
m = ceilSqrt(order)
}
table := make(map[string]*big.Int)
for j := big.NewInt(1); j.Cmp(m) <= 0; j.Add(j, big.NewInt(1)) {
a := new(big.Int).Exp(b, j, n)
table[a.String()] = new(big.Int).Set(j)
}
// p = b^-m modulo n
p := new(big.Int).Neg(m)
p.Exp(b, p, n)
gamma := new(big.Int).Set(x)
for i := big.NewInt(0); i.Cmp(m) == -1; i.Add(i, big.NewInt(1)) {
j, ok := table[gamma.String()]
if ok {
i.Mul(i, m)
i.Add(i, j)
return i, nil
}
gamma.Mul(gamma, p)
gamma.Mod(gamma, n)
}
return nil, errors.New("no solution")
}

40
internal/lib/divisor.go Normal file
View File

@ -0,0 +1,40 @@
/*
Copyright © 2025 filifa
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package lib
import (
"math/big"
)
// TODO: expand to work for different sigmas
func DivisorSummatory(n *big.Int) *big.Int {
// employing Dirichlet's hyperbola method
sqrt := new(big.Int).Sqrt(n)
total := big.NewInt(0)
for x := big.NewInt(1); x.Cmp(sqrt) <= 0; x.Add(x, big.NewInt(1)) {
z := new(big.Int).Div(n, x)
total.Add(total, z)
}
total.Mul(total, big.NewInt(2))
sqrt.Exp(sqrt, big.NewInt(2), nil)
total.Sub(total, sqrt)
return total
}

75
internal/lib/divisors.go Normal file
View File

@ -0,0 +1,75 @@
/*
Copyright © 2025 filifa
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package lib
func pow(base uint, exp uint) uint {
result := uint(1)
for exp > 0 {
if exp%2 == 1 {
result *= base
}
exp >>= 1
base *= base
}
return result
}
func updateMultiples(sieve []uint, x uint, p uint, n uint) {
for q := p; ; q *= p {
// sigma_x(a*b) = sigma_x(a) * sigma_x(b) if gcd(a,b) = 1
for i := 2 * q; i < n; i += q {
if i%(p*q) != 0 {
sieve[i] *= sieve[q]
}
}
if p*q >= n {
break
}
// sigma_x(p^k) = p^(kx) + sigma_x(p^(k-1))
sieve[p*q] = pow(p*q, x) + sieve[q]
}
}
func DivisorsSieve(n uint, x uint) chan uint {
sieve := make([]uint, n)
sieve[0] = 0
for i := uint(1); i < n; i++ {
sieve[i] = 1
}
ch := make(chan uint)
go func() {
for i := uint(0); i < n; i++ {
if i == 0 || i == 1 || sieve[i] != 1 {
ch <- sieve[i]
continue
}
sieve[i] = pow(i, x) + 1
updateMultiples(sieve, x, i, n)
ch <- sieve[i]
}
close(ch)
}()
return ch
}

48
internal/lib/mobius.go Normal file
View File

@ -0,0 +1,48 @@
/*
Copyright © 2025 filifa
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package lib
func MobiusSieve(n uint) chan int {
sieve := make([]int, n)
for i := 0; i < int(n); i++ {
sieve[i] = i
}
ch := make(chan int)
go func() {
for i := 0; i < int(n); i++ {
if i == 0 || i == 1 || sieve[i] != i {
ch <- sieve[i]
continue
}
ch <- -1
for j := 2 * i; j < int(n); j += i {
if j%(i*i) == 0 {
sieve[j] = 0
}
sieve[j] /= -i
}
}
close(ch)
}()
return ch
}

52
internal/lib/partitions.go Normal file
View File

@ -0,0 +1,52 @@
/*
Copyright © 2025 filifa
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package lib
import (
"math/big"
)
// given a slice where vals[i] = p_{k-1}(i) for a given k, nextPartition updates the slice so vals[i] = p_k(i) using the property that p_k(n) = p_k(n-k) + p_{k-1}(n)
func nextPartition(k int, vals []*big.Int) {
for i := 1; i < len(vals); i++ {
if i-k >= 0 {
vals[i].Add(vals[i], vals[i-k])
}
}
}
func Partitions(n, k int) *big.Int {
if n < 0 {
return big.NewInt(0)
}
if k > n {
k = n
}
vals := make([]*big.Int, n+1)
for i := range vals {
vals[i] = big.NewInt(0)
}
vals[0] = big.NewInt(1)
for i := 1; i <= k; i++ {
nextPartition(i, vals)
}
return vals[n]
}

61
internal/lib/primeOmega.go Normal file
View File

@ -0,0 +1,61 @@
/*
Copyright © 2025 filifa
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package lib
func primeOmegaUpdateMultiples(sieve []uint, p uint, n uint, multiplicity bool) {
for q := p; ; q *= p {
// omega(a*b) = omega(a) + omega(b) if gcd(a,b) = 1
for i := 2 * q; i < n; i += q {
if i%(p*q) != 0 {
sieve[i] += sieve[q]
}
}
if p*q >= n {
break
}
if multiplicity {
sieve[p*q] = 1 + sieve[q]
}
}
}
func PrimeOmegaSieve(n uint, multiplicity bool) chan uint {
sieve := make([]uint, n)
for i := uint(0); i < n; i++ {
sieve[i] = 0
}
ch := make(chan uint)
go func() {
for i := uint(0); i < n; i++ {
if i == 0 || i == 1 || sieve[i] != 0 {
ch <- sieve[i]
continue
}
sieve[i] = 1
primeOmegaUpdateMultiples(sieve, i, n, multiplicity)
ch <- sieve[i]
}
close(ch)
}()
return ch
}

33
internal/lib/shoelace.go Normal file
View File

@ -0,0 +1,33 @@
/*
Copyright © 2025 filifa
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package lib
type Point struct {
X float64
Y float64
}
func Area(points []Point) float64 {
total := float64(0)
n := len(points)
for i, p := range points {
q := points[(i+1)%n]
total += (p.Y + q.Y) * (p.X - q.X)
}
return total / 2
}

44
internal/lib/totient.go Normal file
View File

@ -0,0 +1,44 @@
/*
Copyright © 2025 filifa
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package lib
func TotientSieve(n uint) chan uint {
totients := make([]uint, n)
totients[0] = 0
totients[1] = 1
for i := uint(2); i < n; i++ {
totients[i] = i - 1
}
ch := make(chan uint)
go func() {
for i := uint(0); i < n; i++ {
ch <- totients[i]
if i == 0 || i == 1 || totients[i] != i-1 {
continue
}
for j := uint(2 * i); j < n; j += i {
totients[j] -= totients[j] / i
}
}
close(ch)
}()
return ch
}