各位好,遇到一個難題... 比如說要得到根號2小數點100位的值(不是第100位數喔) 其結果越精確越好,精確條件如下 假設求取根號2小數後4位,取其平方後小於該整數之最大值 SOL1:1.4142 < 2 SOL2:1.4143 > 2 所以後4位是取1.4142為最接近值    小弟所知的求取法則有2...    1.夾擠法:可以很精確的求值,原理也簡單,問題是100位數應該是無法以目前的 型別去做運算...    2.開方法:小弟依稀記得有這麼一個方法,後來詢問數學系朋友後得知算法, 過程有點繁複,不過都能以整數直接運算,因此小弟以型別INT進行運算,最終值 再用字串接起來,雖然無法達到無限位數,不過根號2是可以算到小數後224位, 很興奮的使用小算盤驗證結果,小算盤可以顯示到小數後31位,結果在第9位數 就開始出現誤差了,而且位數越多誤差越大...    列出根號2小數以下10位結果 我的程式(開方法):1.4142135601 小算盤:1.4142135623 很明顯的小算盤是非常精確的,因為只要再加0.0000000001平方結果正好大於2!    小算盤只列出31位,或許也是受到型別限制,因此它能使用夾擠法來求值(猜測) 那麼如何做到求取小數後100位呢?    請有興趣的大大來討論討論吧< > 以下是小弟程式,也就是開方法,如果程式寫錯也請不吝指正< >

   
  int Times = Edit3->Text.ToIntDef(0);  //求取的位數
  Memo1->Clear();
  int Input = Edit1->Text.ToIntDef(0);  //欲開方的整數
  int Num = sqrt(Input);        //開方的整數部分
  int work = Input - pow(Num,2);
  int temp=0;
  AnsiString Output;  //最終輸出結果
  Output = Num;
  if (work != 0)
  {
    Output = Output   ".";
    Num  = Num;
    for (int i=0 ; iLines->Add(Output);

程式迷人之處,在於邏輯思考,然而卻也是惱人之處~~

------
程式迷人之處,在於邏輯思考,然而卻也是惱人之處~~

ENIX007 您好：    依稀記得這是我那個年代國中時學過用筆算求平方根的方法，今日再次看到其重現江湖(應該一直都存在著才對)後不禁有種時不我予的感慨！好了言歸正傳！ 錯誤的結果來自於整數的精確度太短！ ENIX007 您的方法似乎是正確的(老實說，是我看不出來哪裡有錯誤)！唯一的問題應該出現在整數的有效位數上。所以只要設法找到能算長整數加減乘的程式碼，取代程式中work,Num,temp的計算，相信您就能得到您需要的結果，理論上應該要幾位就能有幾位！ 希望您早日破案！並請公佈真相，昭告天下！

哈囉,ENIX007:    看完你的解法,我直覺聯想到之前pwipwi先進所發表的一篇文章    大數運算＋分數運算＋多項式運算＋四則運算http://delphi.ktop.com.tw/topic.php?topic_id=49456 原理大致上跟你的一樣,不過有遇到一些困難,我先把我做好的結果先post上來, 有興趣,有時間的人再繼續做下去囉!! 【檔案下載,可以執行但還沒完全解決問題】< href="http://delphi.ktop.com.tw/loadfile.php?TOPICID=25675403&CC=574217">http://delphi.ktop.com.tw/loadfile.php?TOPICID=25675403&CC=574217 目前我只針對根號2,求得testM = 14142135623730,(ps我整數跟小數還沒拆開), testM^2 = 199999999999973116139112900 < 200000000000000000000000000 (嘿嘿,還沒爆掉) ( >, 真可惜,真希望有人可以繼續走下去< >。 有空我會再回來挑戰...< > <> <> >

------


The first step toward proving things for yourself is to understand how others have done it before!

哈囉,ENIX007    終於知道我的方法為什麼會有位數的限制了 針對根號 href="http://delphi.ktop.com.tw/loadfile.php?TOPICID=25677281&CC=574259">http://delphi.ktop.com.tw/loadfile.php?TOPICID=25677281&CC=574259 有問題再討論吧,如有錯誤歡迎賜教。

 請輸入一整數,no = 2
請輸入小數位數長度(例如length=12),length = 17
input = 2
temp = 1
temp3 = 14
sqrt(整數,小數點位數) ≡ sqrt(2, 1) = 1.4
===============================    temp3 = 141
sqrt(整數,小數點位數) ≡ sqrt(2, 2) = 1.41
===============================    temp3 = 1414
sqrt(整數,小數點位數) ≡ sqrt(2, 3) = 1.414
===============================    temp3 = 14142
sqrt(整數,小數點位數) ≡ sqrt(2, 4) = 1.4142
===============================    temp3 = 141421
sqrt(整數,小數點位數) ≡ sqrt(2, 5) = 1.41421
===============================    temp3 = 1414213
sqrt(整數,小數點位數) ≡ sqrt(2, 6) = 1.414213
===============================    temp3 = 14142135
sqrt(整數,小數點位數) ≡ sqrt(2, 7) = 1.4142135
===============================    temp3 = 141421356
sqrt(整數,小數點位數) ≡ sqrt(2, 8) = 1.41421356
===============================    temp3 = 1414213562
sqrt(整數,小數點位數) ≡ sqrt(2, 9) = 1.414213562
===============================    temp3 = 14142135623
sqrt(整數,小數點位數) ≡ sqrt(2, 10) = 1.4142135623
===============================    temp3 = 141421356237
sqrt(整數,小數點位數) ≡ sqrt(2, 11) = 1.41421356237
===============================    temp3 = 1414213562373
sqrt(整數,小數點位數) ≡ sqrt(2, 12) = 1.414213562373
===============================    temp3 = 14142135623730
sqrt(整數,小數點位數) ≡ sqrt(2, 13) = 1.4142135623730
===============================    temp3 = 141421356237308
===極限===
x = 19999999999999574355611086864
sqrt(整數,小數點位數) ≡ sqrt(2, 14) = 1.41421356237309
===============================    temp3 = 1414213562373088
===極限===
x = 1999999999999980062978106655744
sqrt(整數,小數點位數) ≡ sqrt(2, 15) = 1.414213562373099
===============================    temp3 = 14142135623730888
===極限===
x = 199999999999998232571980645268544
sqrt(整數,小數點位數) ≡ sqrt(2, 16) = 1.4142135623730999
===============================    請按任意鍵繼續 . . .

The first step toward proving things for yourself is to understand how others have done it before!

------


The first step toward proving things for yourself is to understand how others have done it before!

ENIX007 您好： 我跑出來的結果： sqrt(2)=1.414213562373095048801688724209698078569671875376948073176679
7379907324784621070388503875343276415727
Maple: 1.4142135623730950488016887242096980785696718753769480731766797379907324784621
07038850387534327641573
您的程式有部份要修改，如紅色部分所標示。我寫了一個類別：LongInt，以便完成上述計算。算出的結果我已與Maple比較，
沒有錯誤！您參考一下！

  int Times = Edit3->Text.ToIntDef(0);  //求取的位數
  Memo1->Clear();
  int Input = Edit1->Text.ToIntDef(0);  //欲開方的整數
  int Num = sqrt(Input);        //開方的整數部分
  int work = Input - pow(Num,2);
  int temp=0;
  AnsiString Output;  //最終輸出結果
  Output = Num;
  if (work != 0)
  {
    Output = Output   ".";
    Num  = Num;
    for (int i=0 ; iNum)  // Num 1
      {
        Output = Output   "0";
        //work *= 100; //多乘一次 
      }
      else
      {
        for (int j=1 ; j<=10 ; j  )
        {
          temp = (Num j)*j;
          if (work < temp)
          {
            temp = (Num (j-1))*(j-1);
            Num  = ((j-1)*2);
            work -= temp;
            Output = Output   IntToStr(j-1);
            break;
          }
        }
      }
    }
  }
  Memo1->Lines->Add(Output);

RichTop 敬上 =====***** 把數學當工具，可以解決問題；將數學變能力，能夠發現並解決問題！ =====#####

各位好: 在使用過大數運算函式庫的時候,因為自己第一次使用,發現一些很容易掉入的陷阱,特地整理如下: 【心得 1】兩個大數比較大小

arbitrary x;
arbitrary y;
if ( x>y ) // 有時候會錯誤

比較保險的作法是用一個中間變數z

arbitrary z;
z = x;
z -= y;

然後再利用(z>0)是否成立去做判斷,真是傑克,太神奇了!! 【心得 2】某大數取平方

    cout << "(1) 詭異的地方(陷阱)==========================================" << endl;
    arbitrary test(" 14142135623731");
    cout << "test = " << test << endl;
    test *= test;
    cout << "(錯誤值)test*=test等於  " << test << endl;
    arbitrary test1(" 14142135623731");
    arbitrary test2(" 14142135623731");
    cout << "test1 = " << test1 << endl;
    cout << "test2 = " << test2 << endl;
    test1 *= test2;
    cout << "(正確值)test1*=test2等於  " << test1 << endl << endl;

經過一番debug,我把結果先貼上來。 sqrt(整數,小數點位數) ≡ sqrt(2, 120) = 1.41421356237309504880168872420969807856967187537694807317667973799073247846210
7038850387534327641572
735013846230912297024 richtop的答案 1.41421356237309504880168872420969807856967187537694807317667973799073247846210
7038850387534327641572
7 真慶幸,經過比較之後: 我的計算結果小數點後100位的數值都跟richtop前輩算的結果相同 驗證部分: x為求得平方根的平方(放大至整數),noScaled為對應x的上界 noScaled = 200000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000
00000000000000000000000000000000000000000 x = 199999999999999999999999999999999999999999999999999999999999999999999999999999
9999999999999999999999
999999999999999999997384168613688492465578826760641151533279013384126427746181
1948353969026636684840
78999911765307750557929677524899999256576 目前測試過到1000位小數點都可以喔 如果大數運算函式庫沒有限制,加上演算法是正確的,應該就可以達到richtop所說的 "理論上應該要幾位就能有幾位！"
【檔案下載-改版2】http://delphi.ktop.com.tw/loadfile.php?TOPICID=25685732&CC=574448
The first step toward proving things for yourself is to understand how others have done it before!

------


The first step toward proving things for yourself is to understand how others have done it before!

ENIX007、fusung 您們好：    其實我很早之前就對計算圓周率PI的值感興趣！不過這同樣需要長整數的加減乘除法運算。 因為看了這個問題，想一起將先前的願望實現，所以寫了類別LongInt來做計算。 多虧ENIX007的程式啟發，目前看來似乎能解決求平方根的問題。不過我是一邊改一邊測試(除法尚未實作)，顯然還存在很多問題才對！ 先提供大家參考，也請把發現的問題告訴我以便修改！ 程式連結如右：http://delphi.ktop.com.tw/loadfile.php?TOPICID=25690427&CC=574553 我曾經用Maple計算相同的結果與10000!，彈指間即完成！其速度之快令人嘆為觀止！ 因此想知道pwipwi大大的算法，會不會比我的快一點，以便直接拿來用！所以想拜託fusung： 由於您所提供的程式碼似乎沒有將pwipwi大大的程式碼加入？我嘗試下載pwipwi大大的程式碼，但不能連結！ 可否請您在下次版本更新時一併附上？我想看看這個版本會不會快一點！ RichTop 敬上 =====***** 把數學當工具，可以解決問題；將數學變能力，能夠發現並解決問題！ =====#####

各位好:    我又來了,特地附上pwipwi的函式庫下載網址。    記得當初是從pwipwi範例中的關鍵字arbi.h從GOOGLE查到下面這個網址。 【C template library。可任意結合大數運算、分數運算、矩陣運算】http://rt.openfoundry.org/Foundry/Project/Download/?Queue=169 蠻期待richtop的測試結果... pwipwi v.s. richtop <> <> >

------


The first step toward proving things for yourself is to understand how others have done it before!

ENIX007、fusung 您們好：    謝謝fusung告知pwipwi大大的程式碼下載點！ 怎麼變成fusung與我在討論？幸好ENIX007還沒回來！ 我試過了之前程式使用 > 看來我可以直接用 > 最後再次謝謝 >

真是抱歉,最近忙著測試,就忘了來看文章了< > 不過看到兩位大大的回應,真是令小弟感動< > 後來看到這篇文章 http://delphi.ktop.com.tw/topic.php?topic_id=23866 於是利用無限位數加法與乘法,配合最單純的夾擠法完成 利用這樣的方法,根號2的100位數需要20秒左右,也測試過999位數需要 3小時多...因此還有很大的成長空間 目前仍忙於測試工作,無暇細看兩位大大的程式,待我有空時再來研究< > 給分方面,實在很難決斷,只好採用亂數選取(花了 > 先謝謝兩位大大的精采討論 < > 以下是程式部分 < class="code"> int OPAdd(AnsiString &s,int id,int value) { //遞回加法運算 AnsiString ss,sv; int Result = 0; if (id <= 0) { sv = IntToStr(value); s = sv s; Result = sv.Length(); } else { // ss = IntToStr(ord(s[id])-48 value); ss = IntToStr(StrToInt(s[id]) value); if (ss.Length() > 1) { s[id] = ss[2]; // Result = OPAdd(s,id-1,ord(ss[1])-48); Result = OPAdd(s,id-1,StrToInt(ss[1])); } else s[id] = ss[1]; } return Result; } //--------------------------------------------------------------------------- AnsiString InfinitAdd(AnsiString s1,AnsiString s2) { //整數版無窮位數加法運算 int i,n1,n2; AnsiString Result; n1 = s1.Length(); n2 = s2.Length(); if (n2 > n1) { Result = InfinitAdd(s2,s1); } else { Result = s1; for (int i=1 ; i<=n2 ; i ) n1 = n1 OPAdd(Result,n1-n2 i,StrToInt(s2[i])); } return Result; } //--------------------------------------------------------------------------- AnsiString InfinitAddPlus(AnsiString s1,AnsiString s2) { //小數點版無窮位數加法運算 int pos1 = s1.Pos("."); int pos2 = s2.Pos("."); int pos; //輸出的小數點位置 AnsiString Result; AnsiString s11=NULL,s12=NULL,s21=NULL,s22=NULL; int n11=0,n12=0,n21=0,n22=0; if (pos1 != 0) { s11 = s1.SubString(1,pos1-1); //小數點前 n11 = s11.Length(); n12 = s1.Length() - s11.Length() - 1; //小數點後 s12 = s1.SubString(pos1 1,n12); } else { s11 = s1; n12 = 1; s1 = s1 "."; } if (pos2 != 0) { s21 = s2.SubString(1,pos2-1); //小數點前 n21 = s21.Length(); n22 = s2.Length() - s21.Length() - 1; //小數點後 s22 = s2.SubString(pos2 1,n22); } else { s21 = s2; n22 = 1; s2 = s2 "."; } //補零 if (n12 < n22) { for (int i=0 ; i n1) Result = InfinitMul(s2,s1); else { n = n1; Result = AnsiString::StringOfChar('0',n1); for ( i=1 ; i<=n2 ; i ) { for ( j=n1 ; j>=1 ; j--) { n = n OPAdd(Result,n-n1-n2 i j,(StrToInt(s2[i]) * StrToInt(s1[j]))); } } } return Result; } //--------------------------------------------------------------------------- AnsiString InfinitMulPlus(AnsiString s1,AnsiString s2) { //小數點版無窮位數乘法運算 int pos1 = s1.Pos("."); int pos2 = s2.Pos("."); int pos; //輸出的小數點位置 AnsiString Result; AnsiString s11=NULL,s12=NULL,s21=NULL,s22=NULL; int n12=0,n22=0; if (pos1 != 0) { s11 = s1.SubString(1,pos1-1); //小數點前 n12 = s1.Length() - s11.Length() - 1; //小數點後 s12 = s1.SubString(pos1 1,n12); } else { s11 = s1; n12 = 1; s1 = s1 "."; } if (pos2 != 0) { s21 = s2.SubString(1,pos2-1); //小數點前 n22 = s2.Length() - s21.Length() - 1; //小數點後 s22 = s2.SubString(pos2 1,n22); } else { s21 = s2; n22 = 1; s2 = s2 "."; } //計算前將小數點拿掉,整數部分為零也不列入計算 if (s11 == "0") s1 = s12; else s1 = s11 s12; if (s21 == "0") s2 = s22; else s2 = s21 s22; Result = InfinitMul(s1,s2); //輸出值得小數點位置 pos = Result.Length()-n12-n22 1; Result.Insert(".",pos); //無整數部分...補零 if (pos <= 1) Result.Insert("0",1); //將小數點後的零去除 int Length = Result.Length(); for (int i=0 ; iText.ToIntDef(0); //求取的位數 Memo1->Clear(); int Input = Edit1->Text.ToIntDef(0); //欲開方的整數 int Num = sqrt(Input); //開方的整數部分 int work = Input - pow(Num,2); AnsiString Output; //最終輸出結果 Output = Num; if (work != 0) { int Pos=0; //小數點位置 AnsiString temp,result; //使用十分逼近法求取平方根 Output = Output "."; for (int i=1 ; i<=Times ; i ) { for (int j=1 ; j<=9 ; j ) { Application->ProcessMessages(); //使系統不鎖住 temp = Output IntToStr(j); result = InfinitMulPlus(temp,temp); Pos = result.Pos("."); result = result.SubString(1,Pos-1); if (result.ToInt() >= Input) //判定小數點前數字大小 { Output = Output IntToStr(j-1); break; } else { if (j == 9) //如果取到9還無法比輸入值大,取9 Output = Output IntToStr(j); } } } } Memo1->Lines->Add(Output); DWORD lasttime = timeGetTime(); Caption = "完成,共花費" IntToStr((lasttime-pretime)/1000) "秒"; } 其中無限位數加乘是取自ccchen大大,小弟翻成了C版本 另外多加了小數點版本,因應程式需要 採用的是十分逼近法,如果改用二分逼近法應該會更快... 不過小弟最近實在沒時間呵呵 再次感謝兩位大大 程式迷人之處,在於邏輯思考,然而卻也是惱人之處~~

------
程式迷人之處,在於邏輯思考,然而卻也是惱人之處~~

參考http://www.swox.com/gmp/ 下載GMP 4.1.4 http://ftp.sunet.se/pub/gnu/gmp/gmp-4.1.4.tar.gz 解開後在gmp-4.1.4/doc/projects.htm說明文件有提到利用迭代的方法求平方根 Faster sqrt     The current code uses divisions, which are reasonably fast, but it'd be possible to use only multiplications by computing 1/sqrt(A) using this formula: 

                                    2
                   x   = x  (3 - A x )/2.
                    i 1   i         i

And the final result

                     sqrt(A) = A x .
                                  n

That final multiply might be the full size of the input (though it might only need the high half of that), so there may or may not be any speedup overall. 另外gmp的說明文件http://www.swox.com/gmp/gmp-man-4.1.4.pdf中第16.5.1 Square Root提到 Square roots are taken using the "Karatsuba Square Root" algorithm by Paul Zimmermann(see Appendix B[References],page 113). This is expressed in a divide and conquer form,... 我另外在http://www.koders.com/java/fidD30CD3D79A65E1D12D57BFC8ABC988E56A861969.aspx?s=karatsubaSquare 找到java的程式碼

    public static BigInteger karatsubaSquare(BigInteger x) {
        BigInteger product = recursiveKaratsubaSquare(x.abs());
        return product;
    }        private static BigInteger recursiveKaratsubaSquare(BigInteger x) {
        assert x.signum() >= 0;            if (x.bitLength() < karatshubaThreshold) {
            return x.pow(2); // faster than x.multiply(x)
        } else {
            BigInteger[] halves = new BigInteger[2];                int byteSize = x.bitLength() >> 4;
            int n = byteSize << 3;                split(x, byteSize, halves);
            BigInteger a = halves[0], b = halves[1];                BigInteger squareA = recursiveKaratsubaSquare(a);
            BigInteger squareB = recursiveKaratsubaSquare(b);                BigInteger aPlusB = a.add(b);
            BigInteger squareAPlusB = recursiveKaratsubaSquare(aPlusB);
            aPlusB = null;                BigInteger product = squareA.add(
                squareAPlusB.subtract(squareA).subtract(squareB).shiftLeft(n)).add(
                    squareB.shiftLeft(n << 1));                return product;
        }

相信這兩種方法都比二分法或是直式開平方法來得快 發表人 - bugmans 於 2005/12/05 21:50:23

謝謝bugmans的補充資料。    A proof of GMP square root using the Coq assistantftp://ftp.inria.fr/INRIA/publication/publi-pdf/RR/RR-4475.pdf Abstract : We present a formal proof (at the implementation level) of an efficient algorithm proposed in to compute square roots of arbitrarily large integers. This program, which is part of the GNU Multiple Precision Arithmetic Library (GMP), is completely proven within the system. Proofs are developed using the Correctness tool to deal with imperative features of the program. The formalization is rather large (more than 13000 lines) and requires some advanced techniques for proof management and reuse. 有興趣可以看一下這個方法的證明 發表人 -

------


The first step toward proving things for yourself is to understand how others have done it before!