Re: [問題] 關於遞迴加快速度的迷思？
: 推 c2251393:公式解的話有個^n，要計算的話還是快不了吧
: → c2251393:(是說再深入下去感覺會離題XDD
不, 計算 n 次方的複雜度是 O(log N), 和 n 的 bit 數成正比,
在 N 大到某個程度以後, 會比前面討論的 O(N) 方法好太多了...
也就是說 x^n 次方的計算方法並非 for (r=1,c=0; c<n; c++) r*=x;
而是可以利用 x^(2n) = (x^n) * (x^n) 遞迴公式去計算...
令 m=n/2 先計算 x^m 再計算 x^n = x^m*x^m (n是偶數) 或 x^n=x^m*x^m*x (奇數)
(不用我說, 看完前文大家應該也知道這種簡單遞迴都是用單層迴圈實作的。)
公式法的問題出在高精確度浮點運算的速度和精確度, 不過這也不是不能改善...
* * *
那麼 Fibonacci 整數數列本身也是可以用類似的二分再二分法去計算的,
我一直在奇怪, 前面怎麼都沒有人提到過... 可能是目前用到的 n 太小的關係...
令函數 f(n) = (a, b, c, d) 傳回值是四個係數, 代表以 F[0]=x, F[1]=y 開頭的
Fibonacci 數列, 其最後兩項為 F[n-1]=ax+by, F[n]=cx+dy,
將兩條數列接起來, 則可推得 f(2n-1) = ((a^2+bc),(ab+bd),(ac+cd),(bc+d^2))
所以我們只需要把數列切成等長的兩半計算, 算完再用上面的公式接起來即可...
* * *
以下這段 code 計算 F[10^7] 不到 3 秒, 運算結果有 2,089,877 位數 (十進位)
使用前面討論的迴圈法要 10-20 分鐘, 你的 CPU 再快大概也不會小於 5 分鐘...
大數計算部份使用 GNU MP 的 mpz 整數運算,
MinGW 或任何版本的 Linux 應該都很容易抓到或者已內建 libgmp 這個套件...
那麼這個二分法能不能改寫成迴圈呢? 可以的, 我的寫法是 head recursive,
可以輕易改成單層迴圈... 有興趣的話各位不妨試試看...
// URL = http://codepad.org/Vn2lmxtz
// -------- >8 -------- >8 -------- >8 -------- >8 -------- >8 --------
#include <stdio.h>
#include <stdint.h>
#include <gmp.h>
#include <inttypes.h>
void f(uint64_t n, mpz_t a_1, mpz_t b_1, mpz_t a_0, mpz_t b_0)
{
uint64_t m;
mpz_t ma_1, mb_1, ma_0, mb_0;
if (n==2) {
mpz_set_ui(a_1, 0);
mpz_set_ui(b_1, 1);
mpz_set_ui(a_0, 1);
mpz_set_ui(b_0, 1);
} else if (n==3) {
mpz_set_ui(a_1, 1);
mpz_set_ui(b_1, 1);
mpz_set_ui(a_0, 1);
mpz_set_ui(b_0, 2);
} else {
mpz_init(ma_1);
mpz_init(mb_1);
mpz_init(ma_0);
mpz_init(mb_0);
m = (n+1)/2;
f(m, ma_1, mb_1, ma_0, mb_0);
// b_1 = (ma_1+mb_0)*mb_1
// a_0 = (ma_1+mb_0)*ma_0
// a_1 = mb_1*ma_0 + ma_1^2
// b_0 = mb_1*ma_0 + mb_0^2
mpz_add(b_1, ma_1, mb_0);
mpz_set(a_0, b_1);
mpz_mul(b_1, b_1, mb_1);
mpz_mul(a_0, a_0, ma_0);
mpz_mul(a_1, mb_1, ma_0);
mpz_set(b_0, a_1);
mpz_pow_ui(ma_1, ma_1, 2);
mpz_add(a_1, a_1, ma_1);
mpz_pow_ui(mb_0, mb_0, 2);
mpz_add(b_0, b_0, mb_0);
if (m*2-1 < n) {
mpz_swap(a_1, a_0);
mpz_swap(b_1, b_0);
mpz_add(a_0, a_1, a_0);
mpz_add(b_0, b_1, b_0);
}
mpz_add(ma_0, a_0, b_0);
mpz_clear(ma_1);
mpz_clear(mb_1);
mpz_clear(ma_0);
mpz_clear(mb_0);
}
}
int main(void)
{
uint64_t n;
mpz_t a_1, b_1, a_0, b_0;
double d;
signed long int exp;
n = 10000000;
mpz_init(a_1);
mpz_init(b_1);
mpz_init(a_0);
mpz_init(b_0);
f(n, a_1, b_1, a_0, b_0);
// Output number itself.
gmp_printf("F(%" PRIu64 ") = %Zd\n", n, b_0);
// Output in x*2^y format.
d = mpz_get_d_2exp(&exp, b_0);
printf("F(%" PRIu64 ") = %f * 2^%ld\n", n, d, exp);
mpz_clear(a_1);
mpz_clear(b_1);
mpz_clear(a_0);
mpz_clear(b_0);
return 0;
}
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 220.137.2.15
推
09/08 09:36, , 1F
09/08 09:36, 1F
討論串 (同標題文章)