squareRoot  1 から 99 の間の平方根の値を多桁数で算出する。

int squareRoot(int *data, int keta, int n);

data  （出力）平方根の値が10進4桁ずつ入った整数の配列
keta  （入力）小数点以下の桁数の指定
n     （入力）計算しようとする平方根の２乗の値

正常に計算できたときに、計算結果となる data の長さ。エラーのとき 0。

int squareRoot(int *data, int keta, int n)
{
    int  *pbuf, *qbuf, *tbuf;
    int  *pp, *qq, *tt;
    int  pketa, qketa;
    int  i;
    int  ret;
    static long p0[101] = {    /* p0*p0 = 4+n*q0*q0 */
        0,-1, 6726, 2702,   -2, 5778, 9602, 4048, 6726,   -3, 1442,
        7940, 2702, 1298,  898, 3842,   -4, 4354, 1154,  340, 5778,
        2525,  394, 2302, 9602,   -5,  102, 2702, 4048,  727,  482,
        3040, 6726, 2114, 4898, 1692,   -6,  146, 5474, 2498, 1442,
        4098,  674, 6964, 7940, 2207,48670L,9214, 2702,   -7,  198,
        9998, 1298, 2599,  970,  178,   898,  302,39206L,1060,3842,
        1523,  126, 4048,   -8,  258,   130,97684L,4354, 623,  502,
        6960, 1154,4562498L,7398,2702,  340, 6239,  106, 160, 5778,
          -9,  326,  164,  110, 6887, 20810, 3134,  394,1000002L,1442,
        3148, 2302,  839,4286590L,6082,9602,125619266L,198,7940,-10 };
    static long q0[101] = {
        0,-1, 4756, 1560,   -2, 2584, 3920, 1530, 2378,   -3,  456,
        2394,  780,  360,  240,  992,   -4, 1056,  272,   78, 1292,
         551,   84,  480, 1960,   -5,   20,  520,  765,  135,   88,
         546, 1189,  368,  840,  286,   -6,   24,  888,  400,  228,
         640,  104, 1062, 1197,  329, 7176, 1344,  390,   -7,   28,
        1400,  180,  357,  132,   24,  120,   40, 5148,  138,  496,
         195,   16,  510,   -8,   32,   16,11934,  528,   75,   60,
         826,  136,534000L,860,  312,   39,  711,   12,   18,  646,
          -9,   36,   18,   12,  747, 2244,  336,   42,106000L,152,
         330,  240,   87,442128L,624,  980,12754704L,20, 798,  -10 };

    if (n < 1 || n > 100 || keta < 1) return 0;
    
    if (p0[n] < 0) {
        *data = -p0[n];
        return 1;
    }

    if (p0[n]*p0[n] != 4 + q0[n]*q0[n]*n)
        return 0;				/* Oh, my God ! */
    ret = 0;
    pbuf = qbuf = tbuf = NULL;
    pbuf = (int *)calloc(2*keta+3, sizeof(int));
    qbuf = (int *)calloc(keta+3, sizeof(int));
    tbuf = (int *)calloc(keta+3, sizeof(int));

    if (pbuf != NULL && qbuf != NULL && tbuf != NULL) {

        if (p0[n] > 10000)
             *pbuf = p0[n]%10000, *(pbuf+1) = p0[n]/10000, pketa = 2;
        else *pbuf = p0[n], pketa = 1;
        if (q0[n] > 10000)
             *qbuf = q0[n]%10000, *(qbuf+1) = q0[n]/10000, qketa = 2;
        else *qbuf = q0[n], qketa = 1;

        while (qketa*4 < keta) {
            /* q = p * q */
            qketa = mulf(tbuf, pbuf, pketa, qbuf, qketa);
            for (tt = tbuf, qq = qbuf, i = qketa; i > 0; i--) *qq++ = *tt++;

            /* p = p * p - 2 */
            pketa = mulf(tbuf, pbuf, pketa, pbuf, pketa);
            for (tt = tbuf, pp = pbuf, i = pketa; i > 0; i--) *pp++ = *tt++;
            *pbuf -= 2;
        }
        ret = divf(data, pbuf, pketa, qbuf, qketa, tbuf);
    }
    if (tbuf != NULL) free(tbuf);
    if (pbuf != NULL) free(qbuf);
    if (qbuf != NULL) free(pbuf);
    return ret;
}

/* mm <-- aa * bb */
static int mulf(int *mm, int *aa, int la, int *bb, int lb)
{
    int  i, j;
    int  *a, *b, *m;
    int  lm, ca;
    long x;

    for (m = mm, i = la + lb; i > 0; i--) *m++ = 0;

    for (b = bb, j = 0; j < lb; j++, b++) {
        for (ca = 0, a = aa, i = 0; i < la; i++) {
            x = *a++;
            x = x * *b + *(mm + i + j) + ca;
            *(mm + i + j) = x % 10000;
            ca = x / 10000;
        }
        *(mm + i + j) = ca;
    }
    return la + lb - (ca == 0);
}         

/* qq <-- aa / bb */
static int divf(int *qq, int *aa, int la, int *bb, int lb, int *tt)
{
    int  i;
    int  *a, *b, *q, *t, *s;
    int  ca, k, Q, lt, lq, m;
    long x;

    for (a = aa + la + lb, s = aa + la, i = la; i > 0; i--) *--a = *--s;
    for (i = lb; i > 0; i--) *--a = 0;
    la += lb;

    if ((k = 4999 / *(bb + lb - 1) + 1) > 1) {
        for (ca = 0, a = aa, t = tt, i = 0; i < la; i++) {
            x = (long)k * *a++ + ca;
            *t++ = x % 10000, ca = x / 10000;
        }
        *t = ca;
        la = i + (ca > 0);
        for (t = tt, a = aa, i = la; i > 0; i--) *a++ = *t++;

        for (ca = 0, b = bb, t = tt, i = 0; i < lb; i++) {
            x = (long)k * *b++ + ca;
            *t++ = x % 10000, ca = x / 10000;
        }
        *t = ca;
        lb = i + (ca > 0);
        for (t = tt, b = bb, i = lb; i > 0; i--) *b++ = *t++;
    }

    lq = la - lb + 1;
    for (q = qq, i = lq; i > 0; i--) *q++ = 0;

    while ((m = la) >= lb) {
        if (*(aa + la - 1) >= *(bb + lb - 1)) {
            for (a = aa + la - 1, b = bb + lb - 1; b != bb; a--, b--) 
                if (*a != *b) break;
            if (b == bb) {
                *(qq + la - lb) = 1;
                la -= lb;
                continue;
            } else if (*a > *b) {
                for (s = b, ca = 0, a = aa + la - lb, b = bb; b <= s; a++) {
                    *a -= *b++ + ca;
                    ca = 0;
                    if (*a < 0) *a += 10000, ca = 1;
                }
                while (*--a == 0) ;
                *(qq + la - lb) = 1;
                la = a - aa + 1;
                continue;
            }
            Q = 9999;
        } else Q = (10000L * *(aa+la-1) + *(aa+la-2)) / *(bb+lb-1);
        if (la == lb) break;
        
        for ( ; ; ) {
            if (Q == 1) {
                *(bb + lb) = 0;
                for (ca=0, a=aa+(la-lb-1), b=bb, i=lb; i >= 0; i--, a++) {
                    *a -= *b++ + ca;
                    ca = 0;
                    if (*a < 0) *a += 10000, ca = 1;
                }
                while (*--a == 0) ;
                la = a - aa + 1;
                break;
            }
            for (ca = 0, t = tt, b = bb, i = 0; i < lb; i++) {
                x = (long)Q * *b++ + ca;
                *t++ = x % 10000, ca = x / 10000;
            }
            *t = ca;
            lt = i + 1;
            for (a = aa + la - 1, t = tt + lt - 1; t != tt; a--, t--)
                if (*a != *t) break;
            if (t == tt) {
                la -= lt;
                break;
            } else if (*a > *t) {
                for (s=t, ca=0, a = aa+(la-lt), t = tt; t <= s; a++) {
                    *a -= *t++ + ca;
                    ca = 0;
                    if (*a < 0) *a += 10000, ca = 1;
                }
                while (*--a == 0) ;
                la = a - aa + 1;
                break;
            } else Q--;
        }
        *(qq + m - lb - 1) = Q;
    }
    if (*(qq + lq - 1) == 0) lq--;

    return lq;
}

　　　　p_i+1 = p_i² - 2
　　　　q_i+1 = p_i q_i

を計算し、最後に、p_i+1 / q_i+1 で平方根を求める。

なお、初期値 p₀、q₀ は ペル方程式 p₀² - M q₀² = 4 を満たす 整数解である。ただし、M は計算しようとする平方根の２乗の値。

数学的帰納法で証明できるが、上の漸化式で得られた p_i+1、 q_i+1 により、

　　　　p_i+1 / q_i+1 = (M + 4 / (p_i q_i)²)^{1 / 2}

となる。 したがって平方根の計算において、本アルゴリズムの方法ではニュートン法等ほかの方法よりも 効率がよく、数回の計算で必要とする精度のものが得られる。