Almkvist-Giullera formula for pi: Difference between revisions
Add Kotlin implementation
m (Minor improvement to code.) |
(Add Kotlin implementation) |
||
| (5 intermediate revisions by 3 users not shown) | |||
Line 1,802:
π to 70 digits is 3.1415926535897932384626433832795028841971693993751058209749445923078164
Computer π is 3.1415926535897932384626433832795028841971693993751058209749445923078164
</pre>
=={{header|Kotlin}}==
{{trans|Java}}
<syntaxhighlight lang="Kotlin">
import java.math.BigDecimal
import java.math.BigInteger
import java.math.MathContext
import java.math.RoundingMode
object CodeKt{
@JvmStatic
fun main(args: Array<String>) {
println("n Integer part")
println("================================================")
for (n in 0..9) {
println(String.format("%d%47s", n, almkvistGiullera(n).toString()))
}
val decimalPlaces = 70
val mathContext = MathContext(decimalPlaces + 1, RoundingMode.HALF_EVEN)
val epsilon = BigDecimal.ONE.divide(BigDecimal.TEN.pow(decimalPlaces))
var previous = BigDecimal.ONE
var sum = BigDecimal.ZERO
var pi = BigDecimal.ZERO
var n = 0
while (pi.subtract(previous).abs().compareTo(epsilon) >= 0) {
val nextTerm = BigDecimal(almkvistGiullera(n)).divide(BigDecimal.TEN.pow(6 * n + 3), mathContext)
sum = sum.add(nextTerm)
previous = pi
n += 1
pi = BigDecimal.ONE.divide(sum, mathContext).sqrt(mathContext)
}
println("\npi to $decimalPlaces decimal places:")
println(pi)
}
private fun almkvistGiullera(aN: Int): BigInteger {
val term1 = factorial(6 * aN) * BigInteger.valueOf(32)
val term2 = BigInteger.valueOf(532L * aN * aN + 126 * aN + 9)
val term3 = factorial(aN).pow(6) * BigInteger.valueOf(3)
return term1 * term2 / term3
}
private fun factorial(aNumber: Int): BigInteger {
var result = BigInteger.ONE
for (i in 2..aNumber) {
result *= BigInteger.valueOf(i.toLong())
}
return result
}
private fun BigDecimal.sqrt(context: MathContext): BigDecimal {
var x = BigDecimal(Math.sqrt(this.toDouble()), context)
if (this == BigDecimal.ZERO) return BigDecimal.ZERO
val two = BigDecimal.valueOf(2)
while (true) {
val y = this.divide(x, context)
x = x.add(y).divide(two, context)
val nextY = this.divide(x, context)
if (y == nextY || y == nextY.add(BigDecimal.ONE.divide(BigDecimal.TEN.pow(context.precision), context))) {
break
}
}
return x
}
}
</syntaxhighlight>
{{out}}
<pre>
n Integer part
================================================
0 96
1 5122560
2 190722470400
3 7574824857600000
4 312546150372456000000
5 13207874703225491420651520
6 567273919793089083292259942400
7 24650600248172987140112763715584000
8 1080657854354639453670407474439566400000
9 47701779391594966287470570490839978880000000
pi to 70 decimal places:
3.1415926535897932384626433832795028841971693993751058209749445923078164
</pre>
Line 2,255 ⟶ 2,344:
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174503
────────────────────────────────────────────── ↑↑↑ the true pi to 160 fractional decimal digits (above) is ↑↑↑ ───────────────────────────────────────────────
</pre>
=={{header|RPL}}==
{{works with|HP|49}}
The starred formula can be implemented like this:
« → n
« 32 6 n * FACT *
532 n SQ * 126 n * + 9 + *
3 / n FACT 6 ^ /
» » '<span style="color:blue">ALMKV</span>' STO
or like that:
« → n '2^5*(6*n)!*(532*SQ(n)+126*n+9)/(3*n!^6)'
» '<span style="color:blue">ALMKV</span>' STO
which is more readable, although a little bit (0.6%) slower than the pure reverse Polish approach.
<code>LDIVN</code> is defined at [[Metallic ratios#RPL|Metallic ratios]]
« 0 → x t
« 0 1
'''WHILE''' DUP x ≤ '''REPEAT''' 4 * '''END'''
'''WHILE''' DUP 1 >
'''REPEAT''' 4 IQUOT
DUP2 + x SWAP - 't' STO
SWAP 2 IQUOT SWAP
'''IF''' t 0 ≥ '''THEN''' t 'x' STO SWAP OVER + SWAP '''END'''
'''END''' DROP
» » '<span style="color:blue">ISQRT</span>' STO
« -105 CF <span style="color:grey">@ set exact mode</span>
« n <span style="color:blue">ALMKV</span> » 'n' 0 9 1 SEQ
-1 → j
« 0 ""
'''DO''' SWAP 'j' INCR <span style="color:blue">ALMKV</span> 10 6 j * 3 + ^ / + EVAL
'''UNTIL''' DUP FXND <span style="color:blue">ISQRT</span> SWAP <span style="color:blue">ISQRT</span> 70 <span style="color:blue">LDIVN</span> ROT OVER ==
'''END'''
"π" →TAG NIP j "iterations" →TAG
» » '<span style="color:blue">TASK</span>' STO
{{out}}
<pre>
3: { 96 5122560 190722470400 7574824857600000 312546150372456000000 13207874703225491420651520 567273919793089083292259942400 24650600248172987140112763715584000 1080657854354639453670407474439566400000 47701779391594966287470570490839978880000000 }
2: π:"3.1415926535897932384626433832795028841971693993751058209749445923078164"
1: iterations:53
</pre>
=={{header|Scala}}==
{{trans|Java}}
<syntaxhighlight lang="Scala">
import java.math.{BigDecimal, BigInteger, MathContext, RoundingMode}
object AlmkvistGiulleraFormula extends App {
println("n Integer part")
println("================================================")
for (n <- 0 to 9) {
val term = almkvistGiullera(n).toString
println(f"$n%1d" + " " * (47 - term.length) + term)
}
val decimalPlaces = 70
val mathContext = new MathContext(decimalPlaces + 1, RoundingMode.HALF_EVEN)
val epsilon = BigDecimal.ONE.divide(BigDecimal.TEN.pow(decimalPlaces))
var previous = BigDecimal.ONE
var sum = BigDecimal.ZERO
var pi = BigDecimal.ZERO
var n = 0
while (pi.subtract(previous).abs.compareTo(epsilon) >= 0) {
val nextTerm = new BigDecimal(almkvistGiullera(n)).divide(BigDecimal.TEN.pow(6 * n + 3))
sum = sum.add(nextTerm)
previous = pi
n += 1
pi = BigDecimal.ONE.divide(sum, mathContext).sqrt(mathContext)
}
println("\npi to " + decimalPlaces + " decimal places:")
println(pi)
def almkvistGiullera(aN: Int): BigInteger = {
val term1 = factorial(6 * aN).multiply(BigInteger.valueOf(32))
val term2 = BigInteger.valueOf(532 * aN * aN + 126 * aN + 9)
val term3 = factorial(aN).pow(6).multiply(BigInteger.valueOf(3))
term1.multiply(term2).divide(term3)
}
def factorial(aNumber: Int): BigInteger = {
var result = BigInteger.ONE
for (i <- 2 to aNumber) {
result = result.multiply(BigInteger.valueOf(i))
}
result
}
}
</syntaxhighlight>
{{out}}
<pre>
n Integer part
================================================
0 96
1 5122560
2 190722470400
3 7574824857600000
4 312546150372456000000
5 13207874703225491420651520
6 567273919793089083292259942400
7 24650600248172987140112763715584000
8 1080657854354639453670407474439566400000
9 47701779391594966287470570490839978880000000
pi to 70 decimal places:
3.1415926535897932384626433832795028841971693993751058209749445923078164
</pre>
Line 2,367 ⟶ 2,565:
{{libheader|Wren-big}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="
import "./fmt" for Fmt
var factorial = Fn.new { |n|
| |||