原題意是, 1+2+...+n >= 10^20,求最小之 n。

 

這問題像是學校作業跑出來的,由於和提問的人之前互動沒很良好,故也沒打算給解答,留參作為紀錄。

 

下面錯誤的想法先提出來

 

錯誤一

int i=0, sum=0;

while(sum < 100000000000000000000) ++i, sum+=i;

printf("i = %d\n", i);

即使用 unsigned long long 都不會過,好的 compiler 在編譯期時就會提出錯誤或警告。

 

錯誤二

double i=0.0, sum=0.0;

while(sum < 1.0E20) ++i, sum+=i;

printf("i = %lf\n", i);

用  double 有效位數只有 52+1 位,到後面 +1 會完全加不上去,卡死。

 

其實不少人有想到,n*(n+1) / 2 >= 10^20 ,二邊取 log

log (n*(n+1) / 2 ) >= 20,

log(n) + log(n+1) >= 20 - log2,

最後問題還是要去解 log(n),實際上可能更麻煩,且不準確。

 

解這題其實不難,但無論如何都會用到大數演算法,若真是用 1+2+...+n > 10 ** 20,其實也才 6 分多鐘 (好吧,我承認老師或助教不會接受這六分多鐘..),另有用 n=1 時,代入 n*(n+1)/2 公式,n 從 1 開始算,一直算到 > 10 ** 20 為止,其實沒什麼兩樣。

 

數值分析。

 

數值分析裡面有個章節叫「非線性方程式求解」,裡面有個方法算是效率最差,但普遍性而言還是非常好、可廣為接受的,叫二分法 (Bisector),即使用二分法暴力破解,最差的情況也只是 log2(up - low) 而已。所以 10**20 ,最差迭代次數也才 20 * log2(10) = 67 次而已。至於效率比較好的牛頓法在這題收斂速度會更快,只是我懶得再用大數實作 non-linear solution。

 

以二分法來做時,只有做幾個部份

 

(1) 設 low = 0 , up= 10**20

(2) 算 n = mid = 0.5* (low + up),算 sum = n * (n+1) / 2,也就是 sum = 1 + 2 +... + n

(3) sum > 10 ** 20 時, up = mid ;反之則 low = mid

(4) 重覆 2, 3 步驟,進行  log2(10**20) 次,即 20 * log2(10) = 67 次。

 

以下程式碼很多都沒優化過,隨便寫寫,但這題很容易跑出 0 secs。

 

/*
 * 
 * filename : BigNum_BiSectorRoot.c
 * author : edison.shih/edisonx
 * compiler : Visual C++ 2008
 * date : 2011.11.26
 *
 * A.L.L. R.I.G.H.T.S. R.E.S.E.R.V.E.
 * 
*/ #include <stdio.h> #include <time.h> #include <math.h> #include <string.h> // ----------------------------------------- #define MAX 12 // 6*2 = 12 #define CAP 10000 #define LEN 4 typedef int Big[MAX]; // ----------------------------------------- // rst = a*b void big_mul(Big rst, Big a, Big b) { int i=0, j, carry, sum; memset(rst, 0, sizeof(int)*MAX); for(i=0; i<MAX; ++i) for(j=0; i+j<MAX; ++j) //if(i+j < MAX) rst[i+j]+=a[i]*b[j]; for(carry=i=0; i<MAX; ++i){ sum = carry + rst[i]; carry = sum / CAP; rst[i] = sum % CAP; } } // ----------------------------------------- // rst = a/b void big_div(Big rst, Big a, int b) { int i=0; int x, remainder=0; for(i=MAX-1; i>=0; --i){ x = remainder * CAP + a[i]; rst[i] = x / b; remainder = a[i] % b; } } // ----------------------------------------- // rst = a+b void big_add(Big rst, Big a, Big b) { int i, carry=0, sum=0; for(i=0; i<MAX; ++i){ sum = a[i] + b[i] + carry; rst[i] = sum % CAP; carry = sum / CAP; } } // ----------------------------------------- // a = a+1 void big_inc(Big a) { int i=0; ++a[i]; while(a[i]==CAP && i<MAX) a[i]=0, ++a[++i]; } // ----------------------------------------- // rst = a+1 void big_add1(Big rst, Big a) // rst = a+1 { memcpy(rst, a, sizeof(int)*MAX); big_inc(rst); } // ----------------------------------------- // big number compare int big_compare(Big a, Big b) { int i=MAX-1; while(a[i]==b[i] && i>=0) --i; if(i==0) return 0; // a==b else if(a[i]>b[i]) return 1; // a>b else return -1; // a<b } // ----------------------------------------- // show the big number void big_show(Big a) { int i=MAX-1; while(a[i]==0 && i>0) --i; printf("%d", a[i]); --i; while(i>=0){ printf("%0*d", LEN, a[i]); --i; } } // ----------------------------------------- // rst = x * (x+1) / 2 void fitness_func(Big rst, Big x) { Big x2, tmp, tmp2; big_add1(x2, x); big_mul(tmp, x, x2); big_div(rst, tmp, 2); } // ----------------------------------------- // main function int main() { int i, j, ret, max_iterator; Big up={0}, low={0}, mid={0}; Big fitness={0}, tmp, ans={0}; Big target={0}; // target value clock_t t = clock(); up[5]=target[5]=1; // set target value and up value max_iterator = (int)(1.0 + log(1.0E20)/log(2.0)); // max itera printf("\n > max_iterator = %d", max_iterator); for(i=0; i<max_iterator; ++i) {// log2(10^20) = 66.xxx big_add(tmp, up, low); // tmp1 = up + low big_div(mid, tmp, 2); // mid = (up+low) / 2 = n fitness_func(fitness, mid); // cal fitness ret = big_compare(fitness, target); if(ret==0) break; else if(ret>0) memcpy(up, mid, sizeof(int)*MAX); else memcpy(low, mid, sizeof(int)*MAX); } // check the last answer if(ret==0) memcpy(ans, mid,sizeof(int)*MAX); // right. else if(ret==1) memcpy(ans, low,sizeof(int)*MAX); // too big else memcpy(ans, up,sizeof(int)*MAX);// too small // calculate fitness again. fitness_func(fitness, ans); // output the result printf("\n > n = "), big_show(ans); printf("\n > 1+2+...+n = "), big_show(fitness); printf(" > "), big_show(target); printf("\n > elaspe : %ld\n", clock()-t); return 0; }

 

總之我認為真正瓶頸不在於大數,而在於求根的演算法。

arrow
arrow
    全站熱搜

    edisonx 發表在 痞客邦 留言(0) 人氣()