From 1a57267a17c2fc17fb6e104846fabc3e363c326c Mon Sep 17 00:00:00 2001 From: Emile Date: Fri, 16 Aug 2024 19:50:26 +0200 Subject: initial commit --- vendor/modernc.org/mathutil/poly.go | 248 ++++++++++++++++++++++++++++++++++++ 1 file changed, 248 insertions(+) create mode 100644 vendor/modernc.org/mathutil/poly.go (limited to 'vendor/modernc.org/mathutil/poly.go') diff --git a/vendor/modernc.org/mathutil/poly.go b/vendor/modernc.org/mathutil/poly.go new file mode 100644 index 0000000..c15e696 --- /dev/null +++ b/vendor/modernc.org/mathutil/poly.go @@ -0,0 +1,248 @@ +// Copyright (c) 2016 The mathutil Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package mathutil // import "modernc.org/mathutil" + +import ( + "fmt" + "math/big" +) + +func abs(n int) uint64 { + if n >= 0 { + return uint64(n) + } + + return uint64(-n) +} + +// QuadPolyDiscriminant returns the discriminant of a quadratic polynomial in +// one variable of the form a*x^2+b*x+c with integer coefficients a, b, c, or +// an error on overflow. +// +// ds is the square of the discriminant. If |ds| is a square number, d is set +// to sqrt(|ds|), otherwise d is < 0. +func QuadPolyDiscriminant(a, b, c int) (ds, d int, _ error) { + if 2*BitLenUint64(abs(b)) > IntBits-1 || + 2+BitLenUint64(abs(a))+BitLenUint64(abs(c)) > IntBits-1 { + return 0, 0, fmt.Errorf("overflow") + } + + ds = b*b - 4*a*c + s := ds + if s < 0 { + s = -s + } + d64 := SqrtUint64(uint64(s)) + if d64*d64 != uint64(s) { + return ds, -1, nil + } + + return ds, int(d64), nil +} + +// PolyFactor describes an irreducible factor of a polynomial in one variable +// with integer coefficients P, Q of the form P*x+Q. +type PolyFactor struct { + P, Q int +} + +// QuadPolyFactors returns the content and the irreducible factors of the +// primitive part of a quadratic polynomial in one variable with integer +// coefficients a, b, c of the form a*x^2+b*x+c in integers, or an error on +// overflow. +// +// If the factorization in integers does not exists, the return value is (0, +// nil, nil). +// +// See also: +// https://en.wikipedia.org/wiki/Factorization_of_polynomials#Primitive_part.E2.80.93content_factorization +func QuadPolyFactors(a, b, c int) (content int, primitivePart []PolyFactor, _ error) { + content = int(GCDUint64(abs(a), GCDUint64(abs(b), abs(c)))) + switch { + case content == 0: + content = 1 + case content > 0: + if a < 0 || a == 0 && b < 0 { + content = -content + } + } + a /= content + b /= content + c /= content + if a == 0 { + if b == 0 { + return content, []PolyFactor{{0, c}}, nil + } + + if b < 0 && c < 0 { + b = -b + c = -c + } + if b < 0 { + b = -b + c = -c + } + return content, []PolyFactor{{b, c}}, nil + } + + ds, d, err := QuadPolyDiscriminant(a, b, c) + if err != nil { + return 0, nil, err + } + + if ds < 0 || d < 0 { + return 0, nil, nil + } + + x1num := -b + d + x1denom := 2 * a + gcd := int(GCDUint64(abs(x1num), abs(x1denom))) + x1num /= gcd + x1denom /= gcd + + x2num := -b - d + x2denom := 2 * a + gcd = int(GCDUint64(abs(x2num), abs(x2denom))) + x2num /= gcd + x2denom /= gcd + + return content, []PolyFactor{{x1denom, -x1num}, {x2denom, -x2num}}, nil +} + +// QuadPolyDiscriminantBig returns the discriminant of a quadratic polynomial +// in one variable of the form a*x^2+b*x+c with integer coefficients a, b, c. +// +// ds is the square of the discriminant. If |ds| is a square number, d is set +// to sqrt(|ds|), otherwise d is nil. +func QuadPolyDiscriminantBig(a, b, c *big.Int) (ds, d *big.Int) { + ds = big.NewInt(0).Set(b) + ds.Mul(ds, b) + x := big.NewInt(4) + x.Mul(x, a) + x.Mul(x, c) + ds.Sub(ds, x) + + s := big.NewInt(0).Set(ds) + if s.Sign() < 0 { + s.Neg(s) + } + + if s.Bit(1) != 0 { // s is not a square number + return ds, nil + } + + d = SqrtBig(s) + x.Set(d) + x.Mul(x, x) + if x.Cmp(s) != 0 { // s is not a square number + d = nil + } + return ds, d +} + +// PolyFactorBig describes an irreducible factor of a polynomial in one +// variable with integer coefficients P, Q of the form P*x+Q. +type PolyFactorBig struct { + P, Q *big.Int +} + +// QuadPolyFactorsBig returns the content and the irreducible factors of the +// primitive part of a quadratic polynomial in one variable with integer +// coefficients a, b, c of the form a*x^2+b*x+c in integers. +// +// If the factorization in integers does not exists, the return value is (nil, +// nil). +// +// See also: +// https://en.wikipedia.org/wiki/Factorization_of_polynomials#Primitive_part.E2.80.93content_factorization +func QuadPolyFactorsBig(a, b, c *big.Int) (content *big.Int, primitivePart []PolyFactorBig) { + content = bigGCD(bigAbs(a), bigGCD(bigAbs(b), bigAbs(c))) + switch { + case content.Sign() == 0: + content.SetInt64(1) + case content.Sign() > 0: + if a.Sign() < 0 || a.Sign() == 0 && b.Sign() < 0 { + content = bigNeg(content) + } + } + a = bigDiv(a, content) + b = bigDiv(b, content) + c = bigDiv(c, content) + + if a.Sign() == 0 { + if b.Sign() == 0 { + return content, []PolyFactorBig{{big.NewInt(0), c}} + } + + if b.Sign() < 0 && c.Sign() < 0 { + b = bigNeg(b) + c = bigNeg(c) + } + if b.Sign() < 0 { + b = bigNeg(b) + c = bigNeg(c) + } + return content, []PolyFactorBig{{b, c}} + } + + ds, d := QuadPolyDiscriminantBig(a, b, c) + if ds.Sign() < 0 || d == nil { + return nil, nil + } + + x1num := bigAdd(bigNeg(b), d) + x1denom := bigMul(_2, a) + gcd := bigGCD(bigAbs(x1num), bigAbs(x1denom)) + x1num = bigDiv(x1num, gcd) + x1denom = bigDiv(x1denom, gcd) + + x2num := bigAdd(bigNeg(b), bigNeg(d)) + x2denom := bigMul(_2, a) + gcd = bigGCD(bigAbs(x2num), bigAbs(x2denom)) + x2num = bigDiv(x2num, gcd) + x2denom = bigDiv(x2denom, gcd) + + return content, []PolyFactorBig{{x1denom, bigNeg(x1num)}, {x2denom, bigNeg(x2num)}} +} + +func bigAbs(n *big.Int) *big.Int { + n = big.NewInt(0).Set(n) + if n.Sign() >= 0 { + return n + } + + return n.Neg(n) +} + +func bigDiv(a, b *big.Int) *big.Int { + a = big.NewInt(0).Set(a) + return a.Div(a, b) +} + +func bigGCD(a, b *big.Int) *big.Int { + a = big.NewInt(0).Set(a) + b = big.NewInt(0).Set(b) + for b.Sign() != 0 { + c := big.NewInt(0) + c.Mod(a, b) + a, b = b, c + } + return a +} + +func bigNeg(n *big.Int) *big.Int { + n = big.NewInt(0).Set(n) + return n.Neg(n) +} + +func bigMul(a, b *big.Int) *big.Int { + r := big.NewInt(0).Set(a) + return r.Mul(r, b) +} + +func bigAdd(a, b *big.Int) *big.Int { + r := big.NewInt(0).Set(a) + return r.Add(r, b) +} -- cgit 1.4.1