# 指定位數的圓周率

December 3, 2021

## 解法思路

PI = [16/5 - 16/(3*5³) + 16/(5*5⁵) - 16/(7*5⁷) + ……] - [4/239 - 4/(3*239³) + 4/(5*239⁵) - 4/(7*239⁷) + ……]

PI = [16/5 - 4/239] - [16/5³ - 4/239³]/3 + [16/5⁵ - 4/239⁵]/5 + ……

[16/5²ⁿ⁻¹ - 4/239²ⁿ⁻¹] / (2*n-1)

[16 /5²ⁿ⁻¹] / (2*n-1) = 10⁻ᶫ

n = L/(2log5) = L/1.39794

n = [L / 1.39794] + 1

[Wₙ₋₁/5²ⁿ - Vₙ₋₁/239²ⁿ] / (2*n-1)

``````div(w, 25, w);
divide(v, 57121, v); // 239 * 239 = 57121
sub(w, v, q);
div(q, 2 * k - 1, q);
``````

## 程式實作：

``````#include <stdio.h>
#include <stdlib.h>

#define L 1000
#define N L / 4 + 1

// L 為位數，N是array長度

// 只處理正數的大整數加、減、除
void subtract(int*, int*, int*);
void divide(int*, int, int*);

int main(void) {
int s[N] = {0};
int w[N] = {0};
int v[N] = {0};
int q[N] = {0};
int n = (int) (L / 1.39793 + 1);

w[0] = 16 * 5;
v[0] = 4 * 239;

int k;
for(k = 1; k <= n; k++) {
// 套用公式
divide(w, 25, w);
divide(v, 57121, v); // 239 * 239 = 57121
subtract(w, v, q);
divide(q, 2 * k - 1, q);

if(k % 2) {// 奇數項
} else {   // 偶數項
subtract(s, q, s);
}
}

printf("%d.", s[0]);
for(k = 1; k < N; k++) {
printf("%04d", s[k]);
}

return 0;
}

void add(int* a, int* b, int* c) {
int i, carry = 0;
for(i = N - 1; i >= 0; i--) {
c[i] = a[i] + b[i] + carry;
if(c[i] < 10000) {
carry = 0;
} else { // 進位
c[i] = c[i] - 10000;
carry = 1;
}
}
}

void subtract(int* a, int* b, int* c) {
int i, borrow = 0;
for(i = N - 1; i >= 0; i--) {
c[i] = a[i] - b[i] - borrow;
if(c[i] >= 0) {
borrow = 0;
} else { // 借位
c[i] = c[i] + 10000;
borrow = 1;
}
}
}

void divide(int* a, int b, int *c) {  // b 為除數
int i, tmp, remain = 0;
for(i = 0; i < N; i++) {
tmp = a[i] + remain;
c[i] = tmp / b;
remain = (tmp % b) * 10000;
}
}
``````
``````import java.math.BigInteger;

public class LongPI {
private final BigInteger PI;

public LongPI(int L) {
int n = (int) (L / 1.39793 + 1);

BigInteger b25 = BigInteger.valueOf(25);
BigInteger b57121 = BigInteger.valueOf(57121);
BigInteger scale = BigInteger.valueOf(10).pow(L);

BigInteger w = BigInteger.valueOf(16 * 5).multiply(scale);
BigInteger v = BigInteger.valueOf(4 * 239).multiply(scale);
BigInteger s = BigInteger.valueOf(0);
BigInteger q = null;
for(int k = 1; k <= n; k++) {
w = w.divide(b25);
v = v.divide(b57121);
q = w.subtract(v).divide(BigInteger.valueOf(2 * k - 1));
s = k % 2 == 1 ? s.add(q) : s.subtract(q);
}
PI = s;
}

public String toString() {
String str = PI.toString();
return String.format("%c.%s", str.charAt(0), str.substring(1));
}

public static void main(String[] args) {
System.out.println(new LongPI(1000));
}
}
``````
``````class LongPI:
def __init__(self, L):
n = int(L / 1.39793 + 1)
scale = 10 ** L

def pi(k, w, v):
if k == n + 1:
return 0
else:
wk = w // 25
vk = v // 57121
qk = (wk - vk) // (2 * k - 1)
return (qk if k % 2 == 1 else -qk) + pi(k + 1, wk, vk)

self.PI = pi(1, 16 * 5 * scale, 4 * 239 * scale)

def __str__(self):
s = str(self.PI)
return "%c.%s" % (s[0], s[1:])

print(LongPI(1000))
``````
``````import scala.math.BigInt

class LongPI private (pi: BigInt) {
override def toString = {
val str = pi.toString
"%c.%s".format(str(0), str.tail)
}
}

object LongPI {
def apply(l: Int) = {
val n = (l / 1.39793 + 1).toInt
val b25 = BigInt(25)
val b57121 = BigInt(57121)
val scale = BigInt(10).pow(l)

def pi(k: Int, w: BigInt, v: BigInt): BigInt = {
if(k == n + 1) BigInt(0)
else {
val wk = w / b25
val vk = v / b57121
val qk = (wk - vk) / BigInt(2 * k - 1)
(if(k % 2 == 1) qk else -qk) + pi(k + 1, wk, vk)
}
}

new LongPI(pi(1, BigInt(16 * 5) * scale, BigInt(4 * 239) * scale))
}
}

print(LongPI(1000))
``````
``````class LongPI
def initialize(l)
n = (l / 1.39793 + 1).to_i
scale = 10 ** l

p = ->(k, w, v) {
if k == n + 1
0
else
wk = w / 25
vk = v / 57121
qk = (wk - vk) / (2 * k - 1)
(if k % 2 == 1; qk else -qk end) + p.call(k + 1, wk, vk)
end
}

@pi = p.call(1, 16 * 5 * scale, 4 * 239 * scale)
end

def to_s
str = @pi.to_s
sprintf("%c.%s", str[0], str[1, str.size])
end
end

puts(LongPI.new(1000))
``````
``````// 用到大數運算中的 BigNumber

var LongPI = function() {
function apply(L) {
var n = parseInt(L / 1.39793 + 1);
var b25 = BigNumber('25');
var b57121 = BigNumber('57121');
var scale = BigNumber('1' + new Array(L + 1).join('0'));

var w = BigNumber('80').multiply(scale);
var v = BigNumber('956').multiply(scale);
var s = BigNumber('0');
for(var k = 1; k <= n; k++) {
w = w.divide(b25);
v = v.divide(b57121);
q = w.subtract(v).divide(BigNumber((2 * k - 1) + ''));
s = k % 2 === 1 ? s.add(q) : s.subtract(q);
}
return new LongPI(s);
}

function LongPI(PI) {
this.PI = PI;
}

LongPI.prototype.toString = function() {
var str = this.PI.toString();
return str[0] + '.' + str.substring(1);
};

return apply;
}();

print(LongPI(50));
``````
``````data LongPI = LongPI Integer

instance Show LongPI where
show (LongPI value) = (head s) : '.' : (tail s)
where s = show value

longPi l = LongPI \$ p 1 (16 * 5 * scale) (4 * 239 * scale)
where n = truncate (fromIntegral l / 1.39793 + 1)
scale = 10 ^ l
p k w v =
if k == n + 1 then 0
else
let wk = w `div` 25
vk = v `div` 57121
qk = (wk - vk) `div` (2 * k - 1)
in (if odd k then qk else -qk) + (p (k + 1) wk vk)

main = print \$ longPi 1000
``````