數學成分濃厚,有關程式計算整數平方根問題(所得位數要很長) |
尚未結案
|
ENIX007
高階會員 發表:28 回覆:274 積分:185 註冊:2003-11-27 發送簡訊給我 |
各位好,遇到一個難題...
比如說要得到根號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 ; i程式迷人之處,在於邏輯思考,然而卻也是惱人之處~~
------
程式迷人之處,在於邏輯思考,然而卻也是惱人之處~~ |
richtop
資深會員 發表:122 回覆:646 積分:468 註冊:2003-06-10 發送簡訊給我 |
ENIX007 您好: 依稀記得這是我那個年代國中時學過用筆算求平方根的方法,今日再次看到其重現江湖(應該一直都存在著才對)後不禁有種時不我予的感慨!好了言歸正傳!
錯誤的結果來自於整數的精確度太短!
ENIX007 您的方法似乎是正確的(老實說,是我看不出來哪裡有錯誤)!唯一的問題應該出現在整數的有效位數上。所以只要設法找到能算長整數加減乘的程式碼,取代程式中work,Num,temp的計算,相信您就能得到您需要的結果,理論上應該要幾位就能有幾位!
希望您早日破案!並請公佈真相,昭告天下!
|
fusung
中階會員 發表:26 回覆:169 積分:99 註冊:2003-11-25 發送簡訊給我 |
哈囉,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! |
fusung
中階會員 發表:26 回覆:169 積分:99 註冊:2003-11-25 發送簡訊給我 |
哈囉,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! |
richtop
資深會員 發表:122 回覆:646 積分:468 註冊:2003-06-10 發送簡訊給我 |
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 ; iRichTop 敬上 =====***** 把數學當工具,可以解決問題;將數學變能力,能夠發現並解決問題! =====##### |
fusung
中階會員 發表:26 回覆:169 積分:99 註冊:2003-11-25 發送簡訊給我 |
各位好: 在使用過大數運算函式庫的時候,因為自己第一次使用,發現一些很容易掉入的陷阱,特地整理如下: 【心得 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! |
richtop
資深會員 發表:122 回覆:646 積分:468 註冊:2003-06-10 發送簡訊給我 |
ENIX007、fusung 您們好: 其實我很早之前就對計算圓周率PI的值感興趣!不過這同樣需要長整數的加減乘除法運算。
因為看了這個問題,想一起將先前的願望實現,所以寫了類別LongInt來做計算。
多虧ENIX007的程式啟發,目前看來似乎能解決求平方根的問題。不過我是一邊改一邊測試(除法尚未實作),顯然還存在很多問題才對!
先提供大家參考,也請把發現的問題告訴我以便修改!
程式連結如右:http://delphi.ktop.com.tw/loadfile.php?TOPICID=25690427&CC=574553 我曾經用Maple計算相同的結果與10000!,彈指間即完成!其速度之快令人嘆為觀止!
因此想知道pwipwi大大的算法,會不會比我的快一點,以便直接拿來用!所以想拜託fusung:
由於您所提供的程式碼似乎沒有將pwipwi大大的程式碼加入?我嘗試下載pwipwi大大的程式碼,但不能連結!
可否請您在下次版本更新時一併附上?我想看看這個版本會不會快一點! RichTop 敬上 =====*****
把數學當工具,可以解決問題;將數學變能力,能夠發現並解決問題!
=====#####
|
fusung
中階會員 發表:26 回覆:169 積分:99 註冊:2003-11-25 發送簡訊給我 |
各位好: 我又來了,特地附上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! |
richtop
資深會員 發表:122 回覆:646 積分:468 註冊:2003-06-10 發送簡訊給我 |
|
ENIX007
高階會員 發表:28 回覆:274 積分:185 註冊:2003-11-27 發送簡訊給我 |
真是抱歉,最近忙著測試,就忘了來看文章了< >
不過看到兩位大大的回應,真是令小弟感動< >
後來看到這篇文章
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
------
程式迷人之處,在於邏輯思考,然而卻也是惱人之處~~ |
bugmans
高階會員 發表:95 回覆:322 積分:188 註冊:2003-04-12 發送簡訊給我 |
參考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 iAnd the final result sqrt(A) = A x . nThat 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 |
fusung
中階會員 發表:26 回覆:169 積分:99 註冊:2003-11-25 發送簡訊給我 |
謝謝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! |
本站聲明 |
1. 本論壇為無營利行為之開放平台,所有文章都是由網友自行張貼,如牽涉到法律糾紛一切與本站無關。 2. 假如網友發表之內容涉及侵權,而損及您的利益,請立即通知版主刪除。 3. 請勿批評中華民國元首及政府或批評各政黨,是藍是綠本站無權干涉,但這裡不是政治性論壇! |