C#程序 計算圓周率

關於圓周率，可參閱“維基百科-圓周率”。我前天在博客園也發表了一篇隨筆：“計算圓周率的C程序” 。這個C#版的計算圓周率程序就是在C程序的基礎上改寫的。C#版的程序必須使用C#2.0編譯，算法和C程序是一樣的，都是利用圓周率的反正切展式的泰勒級數來計算，但C#程序充分使用面象對象的編程方法，並且程序中有適當的注釋，比C程序容易理解多了。C#程序從配置文件中讀取計算所用的公式，允許自己增加計算公式。

using System;
using System.IO;
using System.Text.RegularExpressions;
using System.Collections.Generic;
using System.Diagnostics;

namespace Skyiv.Util
{
// 用反正切展式計算圓周率
// 例如: pi= + 16 * arctan(1/5) - 4 * arctan(1/239) [Machin]
sealed class ThePi
{
// 表示 (positive ? + : -) coefficient * arctan(1/denominator)
public sealed class Term
{
bool positive; // 系數是否正數
int coefficient; // 系數
int denominator; // 分母

public bool Positive { get { return positive; } }
public int Coefficient { get { return coefficient; } }
public int Denominator { get { return denominator; } }

public Term(bool positive, int coefficient, int denominator)
{
this.positive = positive;
this.coefficient = coefficient;
this.denominator = denominator;
}

public override string ToString()
{
return string.Format("{0} {1} * arctan(1/{2}) ", positive ? "+" : "-", coefficient, denominator);
}
}

const int overDigits = 3; // 額外計算的位數，以消除誤差
List<Term> list = new List<Term>(); // 反正切展式
string name; // 反正切展式的發現者

public ThePi(string name)
{
this.name = name;
}

public void Add(Term term)
{
list.Add(term);
}

// 返回反正切展式，例如: pi= + 16 * arctan(1/5) - 4 * arctan(1/239) [Machin]
public override string ToString()
{
string s = "pi= ";
foreach (Term term in list) s += term.ToString();
return s + "[" + name + "]";
}

// 將圓周率精確到小數點後 digits 位的值以及所用時間輸出到 tw
public void Out(TextWriter tw, int digits)
{
const int digitsPerGroup = 5;
const int groupsPerLine = 13;
const int digitsPerLine = groupsPerLine * digitsPerGroup;
tw.WriteLine(this);
TimeSpan elapsed;
int [] piValue = Compute(digits, out elapsed);
tw.Write("pi= {0}.", piValue[piValue.Length - 1]);
int position = 6;
for (int i = piValue.Length - 2; i >= overDigits; i--, position++)
{
tw.Write(piValue[i]);
if (position % digitsPerLine == 0) tw.WriteLine();
else if (position % digitsPerGroup == 0) tw.Write(" ");
}
if ((position - 1) % digitsPerLine != 0) tw.WriteLine();
tw.WriteLine("[{0}] DIGITS:{1:N0} ELAPSED:{2}", name, digits, elapsed);
}

// 計算圓周率到小數點後 digits 位, 並統計所用時間 elapsed
int [] Compute(int digits, out TimeSpan elapsed)
{
Stopwatch stopWatch = new Stopwatch();
stopWatch.Start();
int [] piValue = Compute(digits + overDigits);
Format(piValue);
stopWatch.Stop();
elapsed = stopWatch.Elapsed;
return piValue;
}

// 計算圓周率到小數點後 digits 位, 結果未格式化
int [] Compute(int digits)
{
int [] pi = new int [digits + 1]; // 圓周率的值, 反序存放，如: 951413
int [] tmp = new int [digits + 1]; // 中間計算結果，也是反序存放
foreach (Term term in list)
{
// arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + x^9/9 - x^11/11 + ..
int validDigits = digits;
int divisor = term.Denominator;
bool positive = term.Positive;
Array.Clear(tmp, 0, tmp.Length);
tmp[digits] = term.Coefficient;
Divide(true, positive, true, ref validDigits, pi, tmp, divisor);
positive = !positive;
divisor *= divisor;
for (int step = 3; validDigits > 0; positive = !positive, step += 2)
{
Divide(false, true, true, ref validDigits, null, tmp, divisor);
Divide(true, positive, false, ref validDigits, pi, tmp, step);
}
}
return pi;
}

// 計算 sum += 或 -= (dividend /= 或 / divisor)
void Divide(
bool updateSum, // 是否更新sum
bool positive, // 系數是否正數
bool updateDividend, // 是否更新被除數
ref int digits, // 被除數的有效位數
int [] sum, // 和數
int [] dividend, // 被除數
int divisor // 除數
)
{
for (int remainder = 0, i = digits; i >= 0; i--)
{
int quotient = 10 * remainder + dividend[i];
remainder = quotient % divisor;
quotient /= divisor;
if (updateDividend) dividend[i] = quotient;