package bignum import ( "fmt" "math" "math/big" "github.com/ALTree/bigfloat" ) const pi = "3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989" const log2 = "0.693147180559945309417232121458176568075500134360255254120680009493393621969694715605863326996418687542001481020570685733685520235758130557032670751635075961930727570828371435190307038623891673471123350115364497955239120475172681574932065155524734139525882950453007095326366642654104239157814952043740430385500801944170641671518644712839968171784546957026271631064546150257207402481637773389638550695260668341137273873722928956493547025762652098859693201965058554764703306793654432547632744951250406069438147104689946506220167720424524529612687946546193165174681392672504103802546259656869144192871608293803172714367782654877566485085674077648451464439940461422603193096735402574446070308096085047486638523138181676751438667476647890881437141985494231519973548803751658612753529166100071053558249879414729509293113897155998205654392871700072180857610252368892132449713893203784393530887748259701715591070882368362758984258918535302436342143670611892367891923723146723217205340164925687274778234453534764811494186423867767744060695626573796008670762571991847340226514628379048830620330611446300737194890027436439650025809365194430411911506080948793067865158870900605203468429736193841289652556539686022194122924207574321757489097706753" // Pi returns Pi with prec bits of precision. func Pi(prec uint) *big.Float { pi, _ := new(big.Float).SetPrec(prec).SetString(pi) return pi } func Log2(prec uint) *big.Float { log2, _ := new(big.Float).SetPrec(prec).SetString(log2) return log2 } // NewFloat creates a new big.Float element with "prec" bits of precision. // Valide types for x are: int, int64, uint, uint64, float64, *big.Int or *big.Float. func NewFloat(x interface{}, prec uint) (y *big.Float) { y = new(big.Float) y.SetPrec(prec) // decimal precision if x == nil { return } switch x := x.(type) { case int: y.SetInt64(int64(x)) case int64: y.SetInt64(x) case uint: y.SetUint64(uint64(x)) case uint64: y.SetUint64(x) case float64: y.SetFloat64(x) case *big.Int: y.SetInt(x) case *big.Float: y.Set(x) default: panic(fmt.Errorf("invalid x.(type): valide types are int, int64, uint, uint64, float64, *big.Int or *big.Float but is %T", x)) } return } // Round returns round(x). func Round(x *big.Float) (r *big.Float) { r = new(big.Float).Set(x) if r.Cmp(new(big.Float)) >= 0 { r.Add(r, new(big.Float).SetFloat64(0.5)) } else { r.Sub(r, new(big.Float).SetFloat64(0.5)) } tmp := new(big.Int) r.Int(tmp) r.SetInt(tmp) return } // Cos is an iterative arbitrary precision computation of Cos(x) // Iterative process with an error of ~10^{−0.60206*k} = (1/4)^k after k iterations. // ref : Johansson, B. Tomas, An elementary algorithm to evaluate trigonometric functions to high precision, 2018 func Cos(x *big.Float) (cosx *big.Float) { tmp := new(big.Float) t := NewFloat(0.5, x.Prec()) half := new(big.Float).Copy(t) for i := uint(1); i < (x.Prec()>>1)-1; i++ { t.Mul(t, half) } s := new(big.Float).Mul(x, t) s.Mul(s, x) s.Mul(s, t) four := NewFloat(4.0, x.Prec()) for i := uint(1); i < x.Prec()>>1; i++ { // (1/4)^k = (1/2)^(2*k) tmp.Sub(four, s) s.Mul(s, tmp) } cosx = new(big.Float).Quo(s, NewFloat(2.0, x.Prec())) cosx.Sub(NewFloat(1.0, x.Prec()), cosx) return } func Sin(x *big.Float) (sinx *big.Float) { halfPi := Pi(x.Prec()) halfPi.Quo(halfPi, new(big.Float).SetInt64(2)) return Cos(new(big.Float).Sub(x, halfPi)) } // Log return ln(x) with 2^precisions bits. func Log(x *big.Float) (ln *big.Float) { return bigfloat.Log(x) } // Exp returns exp(x) with 2^precisions bits. func Exp(x *big.Float) (exp *big.Float) { return bigfloat.Exp(x) } // Pow returns x^y func Pow(x, y *big.Float) (pow *big.Float) { return bigfloat.Pow(x, y) } // SinH returns hyperbolic sin(x) with 2^precisions bits. func SinH(x *big.Float) (sinh *big.Float) { sinh = new(big.Float).Set(x) sinh.Add(sinh, sinh) sinh.Neg(sinh) sinh = Exp(sinh) sinh.Neg(sinh) sinh.Add(sinh, NewFloat(1, x.Prec())) tmp := new(big.Float).Set(x) tmp.Neg(tmp) tmp = Exp(tmp) tmp.Add(tmp, tmp) sinh.Quo(sinh, tmp) return } // TanH returns hyperbolic tan(x) with 2^precisions bits. func TanH(x *big.Float) (tanh *big.Float) { tanh = new(big.Float).Set(x) tanh.Add(tanh, tanh) tanh = Exp(tanh) tmp := new(big.Float).Set(tanh) tmp.Add(tmp, NewFloat(1, x.Prec())) tanh.Sub(tanh, NewFloat(1, x.Prec())) tanh.Quo(tanh, tmp) return } // ArithmeticGeometricMean returns the arithmetic–geometric mean of x and y with 2^precisions bits. func ArithmeticGeometricMean(x, y *big.Float) *big.Float { precision := x.Prec() a := new(big.Float).Set(x) g := new(big.Float).Set(y) tmp := new(big.Float) half := NewFloat(0.5, x.Prec()) for i := 0; i < int(math.Log2(float64(precision))); i++ { tmp.Mul(a, g) a.Add(a, g) a.Mul(a, half) g.Sqrt(tmp) } return a } func Sign(x *big.Float) (y *big.Float) { return NewFloat(float64(x.Cmp(NewFloat(0.0, x.Prec()))), x.Prec()) }