squareRoot 1 에서 99 의 사이의 평방근의 값을 다자리수로 산출한다.
int squareRoot(int *data, int keta, int n);
data (출력) 평방근의 값이 10진 4자리수씩 들어간 정수의 배열 keta (입력) 소수점 이하의 자리수의 지정 n (입력) 계산하려고 하는 평방근의 2승의 값
정상적으로 계산할 수 있었을 때에, 계산 결과가 되는 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;
}
pi+1 = pi2 - 2
qi+1 = pi qi
(을)를 계산해, 마지막으로, pi+1 / qi+1 로 평방근을 요구한다.
덧붙여 초기치 p0, q0 은 펠 방정식 p02 - M q02 = 4 를 채운다 정수해이다. 다만, M 는 계산하려고 하는 평방근의 2승의 값.
수학적 귀납법으로 증명할 수 있지만, 위의 점화식에서 얻을 수 있던 pi+1, qi+1 에 의해,
pi+1 / qi+1 = (M + 4 / (pi qi)2)1 / 2
된다. 따라서 평방근의 계산에 대해, 본알고리즘의 방법에서는 뉴턴법등 다른 방법보다 효율이 좋고, 몇차례의 계산으로 필요로 하는 정밀도의 것을 얻을 수 있다.