代码验证斯特林公式的准确性

fliter

发布于 2024-02-26 15:45:02

660

发布于 2024-02-26 15:45:02

关于斯特林公式[1]

斯特林公式（Stirling's approximation或Stirling's formula）是一个用于近似计算阶乘（n!）的公式。当要为某些极大的n求阶乘时，直接计算n!的计算量会随着n的增大而快速增长，导致计算变得不实际，尤其是在计算机程序中。斯特林公式提供了一种有效的方式来近似这种大数的阶乘，能够将求解阶乘的复杂度降低到对数级。

公式如下：

[

n! \approx \sqrt{2 \pi n} \left( \frac{n}{e} \right)^n

]

其中，(e)是自然对数的底数（约等于2.71828），(

\pi

)是圆周率。

斯特林是一位苏格兰数学家,生活在大概清朝康雍乾时期. 不是现在在英超踢快乐足球的那位

但实际上这个公式最早是棣莫弗发现的,就是提出概率论中棣莫弗-拉普拉斯中心极限定理的那位,法国裔英国籍的数学家

更多可参考:

大数定律与中心极限定理[2]

揭示函数渐近展开之奥秘——棣莫弗-拉普拉斯定理的应用与实践[3]

作用

斯特林公式在数学、物理学、工程学和计算机科学等领域有广泛的应用，主要用于：

近似计算大数的阶乘。
在概率论、统计学中估算组合数和排列数。
分析算法的复杂度，特别是那些涉及到阶乘计算的算法。

使用Go代码验证斯特林公式的准确性

如下编写一个简单的Go程序来计算斯特林公式的近似值，并与实际的阶乘值进行比较，以此来验证斯特林公式的准确性

package main

import (
 "fmt"
 "math"
)

// 计算阶乘 n!
func factorial(n int64) int64 {
 if n == 0 {
  return 1
 }
 return n * factorial(n-1)
}

// 斯特林公式近似计算 n!
func stirlingApproximation(n float64) float64 {
 return math.Sqrt(2*math.Pi*n) * math.Pow(n/math.E， n)
}

func main() {
 var a int64 = 5 // 示例：计算5!
 exact := factorial(a)
 approx := stirlingApproximation(float64(a))

 fmt.Printf("Exact factorial of %d is: %d\n"， a， exact)
 fmt.Printf("Stirling's approximation of %d is: %f\n"， a， approx)
 fmt.Printf("Difference: %f\n"， math.Abs(float64(exact)-approx))

 fmt.Printf("误差率: %f\n"， math.Abs(float64(exact)-approx)/float64(exact))

 fmt.Println()

 var b int64 = 10 // 示例：计算10!
 exact = factorial(b)
 approx = stirlingApproximation(float64(b))

 fmt.Printf("Exact factorial of %d is: %d\n"， b， exact)
 fmt.Printf("Stirling's approximation of %d is: %f\n"， b， approx)
 fmt.Printf("Difference: %f\n"， math.Abs(float64(exact)-approx))

 fmt.Printf("误差率: %f\n"， math.Abs(float64(exact)-approx)/float64(exact))

 fmt.Println()

 var c int64 = 20 // 示例：计算20!
 exact = factorial(c)
 approx = stirlingApproximation(float64(c))

 fmt.Printf("Exact factorial of %d is: %d\n"， c， exact)
 fmt.Printf("Stirling's approximation of %d is: %f\n"， c， approx)
 fmt.Printf("Difference: %f\n"， math.Abs(float64(exact)-approx))

 fmt.Printf("误差率: %f\n"， math.Abs(float64(exact)-approx)/float64(exact))

}

在线运行[4]

输出:

Exact factorial of 5 is: 120
Stirling's approximation of 5 is: 118.019168
Difference: 1.980832
误差率: 0.016507

Exact factorial of 10 is: 3628800
Stirling's approximation of 10 is: 3598695.618741
Difference: 30104.381259
误差率: 0.008296

Exact factorial of 20 is: 2432902008176640000
Stirling's approximation of 20 is: 2422786846761136640.000000
Difference: 10115161415503360.000000
误差率: 0.004158

上面程序中，factorial函数递归地计算了一个数的阶乘，而stirlingApproximation函数则根据斯特林公式计算了阶乘的近似值。通过比较两者的结果，可以看到斯特林公式给出的近似值与实际阶乘值之间的差异。

看起来，n越大，斯特林公式计算的结果，和实际n的阶乘值之间的误差会越小。在实际应用中，通常只用斯特林公式来近似计算大数的阶乘。

如果对于非常大的n值，直接计算阶乘可能会导致整数溢出。

2的64次方: 9223372036854775807 20的阶乘: 2432902008176640000

如果n=21，就会超出int64的范围了，需要使用bigint

如下:

package main

import (
 "fmt"
 "math"
 "math/big"
)

func factorial(n int) *big.Int {
 bigN := big.NewInt(int64(n))
 if n == 0 {
  return big.NewInt(1)
 }
 result := big.NewInt(1)
 for i := big.NewInt(2); i.Cmp(bigN) <= 0; i.Add(i， big.NewInt(1)) {
  result.Mul(result， i)
 }
 return result
}

func stirlingApproximation(n float64) float64 {
 return math.Sqrt(2*math.Pi*n) * math.Pow(n/math.E， n)
}

func main() {
 var x = 50
 exact := factorial(x)
 approx := stirlingApproximation(float64(x))

 fmt.Println("Exact factorial of 50 is: "， exact.String())
 fmt.Printf("Stirling's approximation of 50 is: %f\n"， approx)

 floatExact， _ := exact.Float64()
 fmt.Printf("Difference: %f\n"， approx-floatExact)
 fmt.Printf("误差率: %f\n"， math.Abs(floatExact-approx)/floatExact)
 fmt.Println()
 fmt.Println()

 var y = 100
 exact = factorial(y)
 approx = stirlingApproximation(float64(y))

 fmt.Println("Exact factorial of 100 is: "， exact.String())
 fmt.Printf("Stirling's approximation of 100 is: %f\n"， approx)

 floatExact， _ = exact.Float64()
 fmt.Printf("Difference: %f\n"， approx-floatExact)
 fmt.Printf("误差率: %f\n"， math.Abs(floatExact-approx)/floatExact)
 fmt.Println()
 fmt.Println()

 var z = 150
 exact = factorial(z)
 approx = stirlingApproximation(float64(z))

 fmt.Println("Exact factorial of 150 is: "， exact.String())
 fmt.Printf("Stirling's approximation of 150 is: %f\n"， approx)

 floatExact， _ = exact.Float64()
 fmt.Printf("Difference: %f\n"， approx-floatExact)
 fmt.Printf("误差率: %f\n"， math.Abs(floatExact-approx)/floatExact)

}

输出:

Exact factorial of 50 is:  30414093201713378043612608166064768844377641568960512000000000000
Stirling's approximation of 50 is: 30363445939381662539822092816663010554176972669373577771666636800.000000
Difference: -50647262331712294376073382985838127571388053403738488862408704.000000
误差率: 0.001665


Exact factorial of 100 is:  93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
Stirling's approximation of 100 is: 93248476252694086082741480126603978924814282139093327106422842133464418571195672826024908864958629150354343629008809013568909876278442701697442603499744395264.000000
Difference: -77739191250067288427407758913009195131184458378162476302178214323666797065713990911211707211836434286969363646771402787323207418418652286983705207165157376.000000
误差率: 0.000833


Exact factorial of 150 is:  57133839564458545904789328652610540031895535786011264182548375833179829124845398393126574488675311145377107878746854204162666250198684504466355949195922066574942592095735778929325357290444962472405416790722118445437122269675520000000000000000000000000000000000000
Stirling's approximation of 150 is: 57102107404794002531486784484402246707941441176612129409152808169561298568174402178138464969562539359822454467015048872053054189489998393451341243701758138113713430727594194645443385915763391969083296804782229203527249066910085527728440667512127808337057891221504.000000
Difference: -31732159664543345709716527525897500549797926077936066874981976179631513538039613253514105887862044104270936341040554986498420072620437122053086846109932370973210446316192698228766975256370780291479993145852823806527483089444575765382819081148258251866386726912.000000
误差率: 0.000555

二者性能差距

先测算下同样计算100! 共计1000次的总耗时:

package main

import (
 "testing"
 "time"
)

func TestFactorial(t *testing.T) {

 n := 100

 t.Run("Factorial"， func(t *testing.T) {
  start := time.Now()
  for i := 0; i < 1000; i++ {
   factorial(n)
  }
  dur := time.Since(start)
  t.Logf("Factorial duration: %v"， dur)
 })

 t.Run("Stirling"， func(t *testing.T) {
  start := time.Now()
  for i := 0; i < 1000; i++ {
   stirlingApproximation(float64(n))
  }
  dur := time.Since(start)
  t.Logf("Stirling duration: %v"， dur)
 })

}

执行后输出:

=== RUN   TestFactorial
=== RUN   TestFactorial/Factorial
    bench_test.go:18: Factorial duration: 2.608916ms
=== RUN   TestFactorial/Stirling
    bench_test.go:27: Stirling duration: 23.958µs
--- PASS: TestFactorial (0.00s)
    --- PASS: TestFactorial/Factorial (0.00s)
    --- PASS: TestFactorial/Stirling (0.00s)
PASS

再进行benchmark:

package main

import (
 "testing"
)

// 基准测试factorial函数
func BenchmarkFactorial(b *testing.B) {
 for i := 0; i < b.N; i++ {
  factorial(100) // 使用固定的n值进行测试，这里n假设为100
 }
}

// 基准测试stirlingApproximation函数
func BenchmarkStirlingApproximation(b *testing.B) {
 for i := 0; i < b.N; i++ {
  stirlingApproximation(100.0) // 同样使用固定的n值进行测试
 }
}

执行 go test -bench=".*" -benchmem

BenchmarkFactorial-8                      590892              2029 ns/op             240 B/op          6 allocs/op
BenchmarkStirlingApproximation-8        68640559                17.48 ns/op            0 B/op          0 allocs/op

单次执行相差将近100倍，而且n越大，二者的差距越大

「阶乘秘方」斯特林公式：计算大数阶乘的神奇逼近算法[5]

参考资料

[1]

斯特林公式: https://baike.baidu.com/item/%E6%96%AF%E7%89%B9%E6%9E%97%E5%85%AC%E5%BC%8F

[2]

大数定律与中心极限定理: https://zhuanlan.zhihu.com/p/94140620

[3]

揭示函数渐近展开之奥秘——棣莫弗-拉普拉斯定理的应用与实践: https://baijiahao.baidu.com/s?id=1774269727173589266

[4]

在线运行: https://go.dev/play/p/X_FUmDaN84o

[5]