目的

円周率（π）と指数（e）を計算する関数，ユーザが定義するクラスの print method，二項演算子の関数定義の応用など。

使用法

calc.pi(length)
calc.e(length)

引数

length   計算したい小数点以下の桁数

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/multibyte.R", encoding="euc-jp")

# 多倍長計算を行う
"%add%" <- function(ans, b)                          # 足し算の演算子　ans %add% b を行い結果を返す
{                                                       # ans, b は "multibyte" クラスの多倍長整数
        ans <- ans+b                                 # 各桁の足し算を行う
        for (i in length(ans):1) {                      # 各桁について下の桁から，
                if (ans[i] >= 10000000000) {            # 繰り上がり処理を行う
                        ans[i] <- ans[i]-10000000000
                        ans[i-1] <- ans[i-1]+1
                }
        }
        return(ans)                                     # 結果を返す
}
#
"%sub%" <- function(ans, b)                          # 引き算の演算子　ans %sub% b を行い結果を返す
{                                                       # ans, b は "multibyte" クラスの多倍長整数
        ans <- ans-b                                 # 各桁の引き算を行う
        for (i in length(ans):1) {                      # 各桁について下の桁から，
                if (ans[i] < 0) {                    # 繰り下がり処理を行う
                        ans[i] <- ans[i]+10000000000
                        ans[i-1] <- ans[i-1]-1
                }
        }
        return(ans)                                     # 結果を返す
}
#
"%div%" <- function(ans, n)                          # 割り算の演算子　ans %div% n を行い結果を返す
{                                                       # 注：n は "multibyte" クラスではなく普通の整数値
        r <- 0                                               # 剰余
        for (i in 1:length(ans)) {                      # 各桁について上の桁から，
                x <- ans[i]+r*10000000000            # より上の位での剰余を考慮した，被除算数
                ans[i] <- x%/%n                              # 割り算を行い結果を格納
                r <- x-n*ans[i]                              # 今回の剰余
        }
        return(ans)                                     # 結果を返す
}
#
calc.pi <- function(len)                             # πの計算例　小数点以下 len 桁まで求める
{
        len <- len %/% 10+3                          # "multibyte" クラスの多倍長整数の必要桁数
        pi <- a <- b <- numeric(len)                   # 多倍長変数の領域確保　pi <- a <- b <- 0
        a[1] <- 16*5                                 # a[1] が最上位桁（a は小数点以下の数値を格納）
        n <- 1                                               # 除数
        repeat {                                        # 十分な精度を持つまで繰り返し計算
                a <- a %div% 25                              # a <- a / 25
                b <- a %div% n                               # b <- a / n
                pi <- pi %add% b                     # pi <- pi + b
                n <- n+2                             # n <- n + 2
                a <- a %div% 25                              # a <- a / 25
                b <- a %div% n                               # b <- a / n
                pi <- pi %sub% b                     # pi <- pi - b
                n <- n+2                             # n <- n + 2
                if (sum(b) == 0) break                  # pi の値が変化しなくなったらループ終了
        }
        a <- b <- numeric(len)                            # 多倍長変数の領域確保　a <- b <- 0
        a[1] <- 4*239                                        # a <- 4*239
        n <- 1                                               # 除数
        repeat {                                        # 十分な精度を持つまで繰り返し計算
                a <- a %div% 57121                   # a <- a / 57121
                b <- a %div% n                               # b <- a / n
                pi <- pi %sub% b                     # pi <- pi - b
                n <- n+2                             # n <- +2
                a <- a %div% 57121                   # a <- a / 57121
                b <- a %div% n                               # b <- a / n
                pi <- pi %add% b                     # pi <- pi + b
                n <- n+2                             # n <- n + 2
                if (sum(b) == 0) break                  # pi の値が変化しなくなったらループ終了
        }
        class(pi) <- "multibyte"                     # "multibyte" クラスにする（プリント・メソッドを使うため）
        return(pi)                                      # 結果を返す
}
#
calc.e <- function(len)                                      # e の計算例　小数点以下 len 桁まで求める
{
        len <- len %/% 10+3                          # "multibyte" クラスの多倍長整数の必要桁数
        e <- t <- numeric(len)                            # 多倍長変数の領域確保　e <- t <- 0
        e[1] <- t[1] <- 1                         # e <- t <- 1
        i <- 0                                               # 除数
        repeat {                                        # 十分な精度を持つまで繰り返し計算
                i <- i+1                             # i = 1, 2, ...
                t <- t %div% i                               # t = 1 / i!
                if (sum(t) == 0) break                  # t が 0 になるまで
                e <- e %add% t                               # e <- e + t
        }
        class(e) <- "multibyte"                              # "multibyte" クラスにする（プリント・メソッドを使うため）
        return(e)                                       # 結果を返す
}
#
print.multibyte <- function(ans)                     # プリント・メソッド
{
        if (ans[1] == 3) {                              # πの計算結果
                cat("π = 3.\n")
        }
        else {                                          # e の計算結果
                cat("e = 2.\n")
        }
        for (i in 2:(length(ans)-2)) {                  # 小数点以下の答えを出力する（ちょっと冗長だが）
                out <- paste("0000000000", as.character(ans[i]), sep="")
                len <- nchar(out)
                cat(sprintf(" %10s", substring(out, len-9, len)))
                if ((i-1) %% 5 == 0) cat("\n")          # 1行に50桁ずつ出力する
        }
        cat("\n")
}


使用例

# πを1000桁計算

 > calc.pi(1000)
 π = 3.
  1415926535 8979323846 2643383279 5028841971 6939937510
  5820974944 5923078164 0628620899 8628034825 3421170679
  8214808651 3282306647 0938446095 5058223172 5359408128
  4811174502 8410270193 8521105559 6446229489 5493038196
  4428810975 6659334461 2847564823 3786783165 2712019091
  4564856692 3460348610 4543266482 1339360726 0249141273
  7245870066 0631558817 4881520920 9628292540 9171536436
  7892590360 0113305305 4882046652 1384146951 9415116094
  3305727036 5759591953 0921861173 8193261179 3105118548
  0744623799 6274956735 1885752724 8912279381 8301194912
  9833673362 4406566430 8602139494 6395224737 1907021798
  6094370277 0539217176 2931767523 8467481846 7669405132
  0005681271 4526356082 7785771342 7577896091 7363717872
  1468440901 2249534301 4654958537 1050792279 6892589235
  4201995611 2129021960 8640344181 5981362977 4771309960
  5187072113 4999999837 2978049951 0597317328 1609631859
  5024459455 3469083026 4252230825 3344685035 2619311881
  7101000313 7838752886 5875332083 8142061717 7669147303
  5982534904 2875546873 1159562863 8823537875 9375195778
  1857780532 1712268066 1300192787 6611195909 2164201989

# e を1000桁計算
> calc.e(1000)
 e = 2.
  7182818284 5904523536 0287471352 6624977572 4709369995
  9574966967 6277240766 3035354759 4571382178 5251664274
  2746639193 2003059921 8174135966 2904357290 0334295260
  5956307381 3232862794 3490763233 8298807531 9525101901
  1573834187 9307021540 8914993488 4167509244 7614606680
  8226480016 8477411853 7423454424 3710753907 7744992069
  5517027618 3860626133 1384583000 7520449338 2656029760
  6737113200 7093287091 2744374704 7230696977 2093101416
  9283681902 5515108657 4637721112 5238978442 5056953696
  7707854499 6996794686 4454905987 9316368892 3009879312
  7736178215 4249992295 7635148220 8269895193 6680331825
  2886939849 6465105820 9392398294 8879332036 2509443117
  3012381970 6841614039 7019837679 3206832823 7646480429
  5311802328 7825098194 5581530175 6717361332 0698112509
  9618188159 3041690351 5988885193 4580727386 6738589422
  8792284998 9208680582 5749279610 4841984443 6346324496
  8487560233 6248270419 7862320900 2160990235 3043699418
  4914631409 3431738143 6405462531 5209618369 0888707016
  7683964243 7814059271 4563549061 3031072085 1038375051
  0115747704 1718986106 8739696552 1267154688 9570350354