yukicoder No.313 π

問題
No.313 π - yukicoder
20万桁ある円周率のうち，一箇所だけ間違いがある．
どの数字が間違っているか指摘する問題．

解法1
（解法1では解けません）
円周率を20万桁まで計算して比較するという方法を取ろうとしてしまった．
しかし5秒以内に計算するとなったら1万桁が限界だったので，諦めて解説を見たら出てくる数を数えるだけだったのでショックを受けている．
頭が固かったです．
折角なので1万桁を計算する方法を紹介します．
Machin の公式

を使います．arctanは

と表されるので，次数の低いものから足しあわせていけば円周率に近づけることができる．
しかしこれでは収束が遅いので
高野喜久雄の公式

を使ったが20万桁は達成できなかった．

解法2
ソースコードに0~9の出現する回数を埋め込んでおきそれと比較することで間違いを検出する．

ソースコード1

#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;

#define RREP(i,s,e) for (i = s; i >= e; i--)
#define rrep(i,n) RREP(i,(int)(n)-1,0)
#define REP(i,s,e) for (i = s; i <= e; i++)
#define rep(i,n) REP(i,0,(int)(n)-1)
#define INF 100000000
#define MAX (ll)1e8
#define SIZE 1254

typedef long long ll;
typedef vector<ll> vl;

void div(vl& a, int d, vl& b) {
    int i;
    ll down = 0;
    rep (i,min(a.size(),b.size())) {
        ll tmp = (down * MAX + a[i]) / d;
        down = (down * MAX + a[i]) % d;
        b[i] = tmp;
    }
}

void add(vl& a, vl& b, vl& c) {
    int i, up = 0;
    rrep (i,min({a.size(),b.size(),c.size()})) {
        c[i] = a[i] + b[i] + up;
        up = 0;
        while (c[i] < 0 && a[0] + b[0] > 0 || c[i] <= -MAX) {
            up -= 1;
            c[i] += MAX;
        }
        while (c[i] > 0 && a[0] + b[0] < 0 || c[i] >= MAX) {
            up += 1;
            c[i] -= MAX;
        }
    }
}

void sub(vl& a, vl& b, vl& c) {
    int i, up = 0;
    rrep (i,min({a.size(),b.size(),c.size()})) {
        c[i] = a[i] - b[i] + up;
        up = 0;
        while (c[i] < 0 && a[0] + b[0] > 0 || c[i] <= -MAX) {
            up -= 1;
            c[i] += MAX;
        }
        while (c[i] > 0 && a[0] + b[0] < 0 || c[i] >= MAX) {
            up += 1;
            c[i] -= MAX;
        }
    }
}

void print(vl a) {
    int i;
    cout << a[0] << ".";
    REP (i,1,a.size()-4)
        printf("%08lld",a[i]);
    printf("\n");
}

int main() {
    int i, j;
    vl pi(SIZE), a(SIZE), b(SIZE), c(SIZE), d(SIZE), e(SIZE);
    bool endd = false;
    a[0] = 48*49; b[0] = 128*57; c[0] = 20*239; d[0] = 48*110443;
    //a[0] = 16*5; b[0] = 4*239;
    for (i = 0; ; i++) {
        int j;
        int n = 2 * i + 1;
        /*
        div(a,25,a);
        div(b,239*239,b);
        sub(a,b,c);
        div(c,n,c);
        if (i % 2 == 0)
            add(pi,c,pi);
        else
            sub(pi,c,pi);
        bool end = true;
        rep (j,c.size()) end &= ~c[j];
        // */
        // /*
        div(a,49*49,a);
        div(b,57*57,b);
        div(c,239*239,c);
        div(d,110443,d);
        div(d,110443,d);
        add(a,b,e);
        sub(e,c,e);
        add(e,d,e);
        div(e,n,e);
        if (i % 2 == 0)
            add(pi,e,pi);
        else
            sub(pi,e,pi);
        bool end = true;
        rep (j,e.size()) end &= ~e[j];
        // */
        if (end)
            break;
    }
    print(pi);
    return 0;
}

ソースコード2

#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;

#define RREP(i,s,e) for (i = s; i >= e; i--)
#define rrep(i,n) RREP(i,(int)(n)-1,0)
#define REP(i,s,e) for (i = s; i <= e; i++)
#define rep(i,n) REP(i,0,(int)(n)-1)
#define INF 100000000

typedef long long ll;

int main() {
    int i, fault, ok;
    int currect[10] = {20104,20063,19892,20011,19874,20199,19898,20163,19956,19841};
    int cnt[10] {};
    string s;
    cin >> s;
    for (auto c : s) {
        if (c != '.')
            cnt[c-'0']++;
    }
    rep (i,10) {
        if (cnt[i] == currect[i] - 1)
            ok = i;
        else if (cnt[i] == currect[i] + 1)
            fault = i;
    }
    cout << fault << " " << ok << endl;
    return 0;
}