codevs3119 (大整数开根 牛顿迭代法)

http://codevs.cn/problem/3119/

Description

给出一个正整数n，求n开根号后的整数部分的值。n的位数不超过1000位。

Input Description

读入一个不超过1000位的正整数n。

Output Description

Sample Input

17

Sample Output

4

Data Size & Hint

n为不超过1000位的正整数 满足 1<=n <101000

思路

对于大整数我们开根号时，若采用double 精度有限 且效率不是很快

首先我们发现其具有单调性，即可以二分

二分法$O(n^3)$

其中二分次数为$O(3.32*n)$ 每一次二分为$O(2n^2+2n)$

对于本题n1000时间勉强

有没有更快逼近这样点的方法呢？

关于牛顿迭代法

https://en.wikipedia.org/wiki/Newton%27s_method

https://www.zhihu.com/question/20690553

牛顿迭代法 $O(n^2logn)$

其中迭代次数为$log(n)$ 每一次迭代为$O(2n^2+2n)$

代码示例

//高精度(压四位)
#include<bits/stdc++.h>
using namespace std;

const int ten[4]={1,10,100,1000};//辅助压位 压四位 这样乘法不会超int
const int maxl=1500; //最大位数

struct BigNumber{
    int d[maxl];//每4位压位 d[0]为压位后的位数 d[1]为原先的最低四位
    BigNumber (string s){
        int len=s.size();
        d[0]=(len-1)/4+1;
        int i,j,k;
        for(i=1;i<maxl;++i) d[i]=0;
        for(i=len-1;i>=0;--i){
            j=(len-i-1)/4+1;
            k=(len-i-1)%4;
            d[j]+=ten[k]*(s[i]-'0');
        }
        while(d[0]>1 && d[d[0]]==0) --d[0];//去除前导零
    }
    BigNumber(){
        *this=BigNumber(string("0"));
    }
    string toString(){
        string s("");
        int i,j,temp;
        //先处理高1位
        for(i=3;i>=1;--i) if(d[d[0]]>=ten[i]) break;
        temp=d[d[0]];
        for(j=i;j>=0;--j){
            s=s+(char)(temp/ten[j]+'0');
            temp%=ten[j];
        }
        //再处理剩下的位
        for(i=d[0]-1;i>0;--i){
            temp=d[i];
            for(j=3;j>=0;--j){
                s=s+(char)(temp/ten[j]+'0');
                temp%=ten[j];
            }
        }
        return s;
    }
}zero("0"),d,mid1[15];//d为作除法时的余数（可看做是高精度模高精度） mid1辅助作除法


bool operator <(const BigNumber &a,const BigNumber &b){
    if(a.d[0]!=b.d[0]) return a.d[0]<b.d[0];
    int i;
    for(i=a.d[0];i>0;--i) if(a.d[i]!=b.d[i]) return a.d[i] < b.d[i];
    return false;
}

//一定没有前导零
BigNumber operator +(const BigNumber &a,const BigNumber &b){
    BigNumber c;
    c.d[0]=max(a.d[0],b.d[0]);
    int i,x=0;
    for(i=1;i<=c.d[0];++i){
        x=a.d[i]+b.d[i]+x;
        c.d[i]=x%10000;
        x/=10000;
    }
    while(x!=0)
    {
        c.d[++c.d[0]]=x%10000;
        x/=10000;
    }
    return c;
}

BigNumber operator - (const BigNumber &a,const BigNumber &b){
    BigNumber c;
    c.d[0]=a.d[0];
    int i,x=0;
    for(i=1;i<=c.d[0];++i){
        x=10000+a.d[i]-b.d[i]+x;
        c.d[i]=x%10000;
        x=x/10000-1;
    }
    while((c.d[0]>1)&&(c.d[c.d[0]]==0)) --c.d[0];
    return c;
}

BigNumber operator *(const BigNumber &a,const int &k){//乘以较小数
    BigNumber c;
    c.d[0]=a.d[0];
    int i,x=0;
    for(i=1;i<=a.d[0];++i){
        x=a.d[i]*k+x;
        c.d[i]=x%10000;
        x/=10000;
    }
    while(x>0)
    {
        c.d[++c.d[0]]=x%10000;
        x/=10000;
    }
    while((c.d[0]>1) && (c.d[c.d[0]]==0) ) --c.d[0];
    return c;
}


BigNumber operator *(const BigNumber &a,const BigNumber &b){
    BigNumber c;
    c.d[0]=a.d[0]+b.d[0];
    int i,j,x;
    for(i=1;i<=a.d[0];++i){
        x=0;
        for(j=1;j<=b.d[0];++j){
            x=a.d[i]*b.d[j]+x+c.d[i+j-1];
            c.d[i+j-1]=x%10000;
            x/=10000;
        }
        c.d[i+b.d[0]]=x;//保留进位 可能产生前导零
    }
    while((c.d[0]>1)&&(c.d[c.d[0]]==0)) --c.d[0];
    return c;
}

//辅助除法
bool smaller(const BigNumber & a,const BigNumber &b ,int delta){
    if(a.d[0]+delta!=b.d[0]) return a.d[0]+delta<b.d[0];
    //a为试减的数 b为余数
    //如果b的位数和a+delta是相等的
    int i;
    //从最高位开始
    for(i=a.d[0];i>0;--i) if(a.d[i]!=b.d[i+delta]) return a.d[i]<b.d[i+delta];//i+delta为b高位对应a的数
    return true;//都相等返回true
}

//辅助除法
void Minus(BigNumber &a,const BigNumber &b,int delta){
    //引用a 即余数   b为要减去的数 即除数的倍数
    int i,x=0;
    //a最高的a.d[0]-delta位减去b
    for(i=1;i<=a.d[0]-delta;++i){
        x=10000+a.d[i+delta]-b.d[i]+x;
        a.d[i+delta]=x%10000;
        x=x/10000-1;
    }
    while((a.d[0]>1) && (a.d[a.d[0]]==0)) --a.d[0];
}

BigNumber operator / (const BigNumber &a,const BigNumber &b){
    BigNumber c;
    d=a;//余数
    int i,j,temp;
    mid1[0]=b;
    for(i=1;i<=13;++i){
        mid1[i]=mid1[i-1]*2;
    }
    for(i=a.d[0]-b.d[0];i>=0;--i){
        temp=8192;//2^13   倍增思想
        for(j=13;j>=0;--j){ //mid1[j] 为 除数b*2^j
            if(smaller(mid1[j],d,i)){
                Minus(d,mid1[j],i);
                c.d[i+1]+=temp;//i+1是商的位权 因为i是delta
            }
            temp/=2;
        }
    }
    c.d[0]=max(1,a.d[0]-b.d[0]+1);
    while((c.d[0]>1) && (c.d[c.d[0]]==0)) --c.d[0];
    return c;
}

BigNumber tp;
string s;
BigNumber f(BigNumber x)
{
    return x*x-tp;
}

BigNumber deri(BigNumber x)
{
    return x*2;
}

BigNumber newton_sqrt(BigNumber a)
{
    int num=a.d[0]*4;
    string tmp="1";
    for(int i=0;i<num/2+1;++i){
        tmp+="0";
    }//开始的a大约是原位数的一半
    a=BigNumber(tmp);
    num=log(num)+20;
    while(num--)
    {
        a=a-f(a)/deri(a);
        //cout<<num<<"******"<<a.toString()<<endl;
    }
    while(tp<a*a){
        a=a-BigNumber("1");
    }
    return a;
}

BigNumber bina_sqrt(BigNumber a)
{
    BigNumber l("0");
    BigNumber r(s);
    BigNumber cp("1");
    while(cp<r-l)
    {
        BigNumber mid=(l+r)/BigNumber("2");
        //cout<<"******"<<mid.toString()<<endl;
        if(mid*mid<tp) l=mid;
        else r=mid;
    }
    while(tp<r*r){
        r=r-BigNumber("1");
    }
    return r;
}

int main()
{
    freopen("in.txt","r",stdin);
    ios::sync_with_stdio(false);
    cin>>s;
    BigNumber a(s);
    tp=a;
    BigNumber ans=newton_sqrt(a);
    cout<<ans.toString()<<endl;
    cout<<clock()<<endl;
    cout<<endl<<"**********************************"<<endl;

    BigNumber bns=bina_sqrt(a);
    cout<<bns.toString()<<endl;
    cout<<clock()<<endl;

    cout<<endl<<"**********************************"<<endl;

    char temp[3000];
    strcpy(temp,s.c_str());
    double d;
    d=atof(temp);
    cout<<floor(sqrt(d))<<endl;
    cout<<clock()<<endl;
    return 0;
}