數學/手開平方/sgu111 Very simple problem

題意

　　給出一個數n，求sqrt(n) (1≤n≤10¹⁰⁰⁰)

分析

　　題意很簡單，就是開一個數的平方

　　在網上看了看一些方法，一下摘自“風中落葉”hi.baidu.com/xiamengy

1.舉例

上式意為65536的開平方為256。手開方過程類似於除法計算。為了方便表述，以下仍稱類似位置的數為“被除數”、“除數”、“商”。

以65536為例，其具體計算過程如下：

Step1：將被開方數(為了形象，表述成“被除數”，此例中即為65536)從個位往高位每兩位一斷寫成6,55,35的形式，為了方便表述，以下每一個“,”稱為一步

Step2：從高位開始計算開方。例如第一步為6，由於2²=4<6<9=3²，因此只能商2(這就是和除法不同的地方，“除數”和“商”的計算位必須相同)。於是將2寫在根號上方，計算開方余項。即高位余項加一步低位，此例中，即為高位余項2和低位一步55，余項即為255。

Step3：將Step2得到的第一步開方得數2乘以20(原理在后面證明)作為第二步除數的高位。即本步除數是4x(四十幾)。按照要求，本步的商必須是x。因為45×5=225<255<46×6=276，所以本步商5。

Step4：按照類似方法，繼續計算以后的各步。其中，每一步的除數高位都是20×已求出的部分商。例如第三步的除數高位就是25×20=500，所以第三步除數為50x。本例中，506×6=3036恰好能整除，所以256就是最終計算結果。

 

2.字母表示和手開方公式的證明：

既然要證明，必須先把公式一般化。簡言之，用字母而不是特殊值來表示計算過程和結果。

任意正整數均可表示成

則正整數M開方計算得到的就是A。根據手開方公式的思路，應該寫成：

失一般性，對A進行推廣。前面A表示正整數，現在A可以表示任意實數。因為計算開平方問題上，對於數值，正負是無所謂的。

因此不妨假設A為任意正實數。即可記

(即用科學計數法表示，例如134.87可以表示為

1.3487×10²=(1+3×0.1+4×0.01+8×0.001+7×0.0001)×10²)

如此，每一步的開方余項都用該步的“除數”和“商”表示出來，因此，手開方公式是精確的。

再進一步推廣，對於更一般的情況，即使開方結果是無理數，或循環小數的，只需令n→∞即可。由於以上證明對任意n均成立，可以推得對於n→∞也相應成立。

 

上文的摘抄介紹了手開平方的方法以及證明，下面來說一下sgu111具體實現

這個我們主要其實是要解決“除數”

根據上文的證明可知，對於開平方，我們需要找到某個數(1~9)，使得上一步的“商”×20+這個數，能盡可能地接近“開方余項”

設“這個數”為y，上一步的商為x

 

如上圖的45=2*20+5

當(2*20+y)*y盡可能逼近255時，這個y可取。y=5時，(2*20+5)*5=225<255；y=6時，(2*20+6)*6=276>255，不符合，所以y取5

類似地，上圖中的506也可以由(25*20+y)*y得到，只不過此時當y=6時，能正好取到3036

從高位求到低位的開方算法其實是取值范圍的划定。

如果當前求當前位，設前面開好的一個大整數x，新數就是x*10+y，設當前的被開方數是a，必有：
a <= (x*10+y)² = x²*100+20xy+y² = x²+y*(20x+y)

根據這個算法的思想，我們要在上述不等式中使y求到最大，如果滿足條件的情況下，y不是最大的，不論后面開出來的數有多大，結果都是有在x的數位前有很大的數字無法開盡，更談不上逐步求解，縮小誤差了。

具體實現：
首先我們給數字分位，要考慮奇數位的特殊情況。
一般的高精度a[1]是最低位，我這里反着存。
b[i-1]:=b[i-1]*2 這句話其實是b=20x

b[i]=j是b=20x+y

calc函數用來判斷此時的y*(20x+y)能否取到最大值

Accepted Code

由於雲昊哥最近旅游去了,所以...沒人給我改C++的程序了....先用pascal寫吧~

{
        PROBLEM:sgu111
        AUTHER:Rinyo
        MEMO:手開平方
}
Program sgu111;
Const
  Infile = 'sgu111.in';
  Outfile = 'sgu111.out';
Var
  a,b,c:Array[-1..1030]Of Longint;
  len,i,j:Longint;
  ch:Char;
Function calc(x,y:Longint):Boolean;
Var
  t,i:Longint;
Begin
  t:=0;
  For i:=x DOwnto 0 DO Begin
    c[i]:=b[i]*y+t;
    t:=c[i] Div 10;
    c[i]:=c[i] Mod 10;
  End;
  c[-1]:=t;
  For i:=-1 To x Do
    If c[i]>a[i+x] Then Begin
      calc:=false;
      Exit;
    End Else If c[i]<a[i+x] Then Begin
      calc:=true;
      Exit;
    End;
  calc:=true;
ENd;

Begin
  Assign(input,infile);Reset(input);
  Assign(output,outfile);Rewrite(output);
  Fillchar(a,sizeof(a),0);
  Fillchar(b,sizeof(b),0);
  len:=0;
  While not eoln Do Begin
    Inc(len);Read(ch);a[len]:=ord(ch)-ord('0');
  End;
  If odd(len) Then Begin
    For i:=len+1 Downto 1 Do a[i]:=a[i-1];
    Inc(len);
  End;
  For i:=1 To len Div 2 Do Begin
    b[i-1]:=b[i-1]*2;
    If b[i-1]>9 Then Begin b[i-1]:=b[i-1]-10;b[i-2]:=b[i-2]+1; End;
    For j:=9 Downto 0 Do Begin
      b[i]:=j;
      If calc(i,j) Then Begin Write(j);break; End;
    End;
    For j:=i*2 Downto i-1 Do Begin
      Dec(a[j],c[j-i]);
      If a[j]<0 Then Begin Dec(a[j-1]);Inc(a[j],10); End;
    End;
  End;
  WriteLn;
  Close(input);Close(output);
End.