HDU5133(2014年广州现场赛G题)

2014年的两场比赛分别被两道G题卡成狗. 牡丹江的那场比完赛一周内就过了, 结果广州这场的G题拖到三年后的今天终于也过了. 可喜可贺...

题意: 给定一个九位以内的小数F (0<F<1), 求最小的P和Q满足 (P/Q)^2 在四舍五入之后等于F.

思路: 在14年的暑假集训中有遇到过类似的题目, 只不过那题没有开方. 这道题可以这么做. 首先, 四舍五入可以确定P/Q的范围, 比方说0.50=50/100, 那么P/Q应该在[sqrt(495/1000), sqrt(505/1000)) 之间. 我们称这个区间为我们的目标区间. 那么, 从[0/1, 1/1) 开始, 对于每个区间[ln/ld, rn/rd), 我们考虑 (ln+rn)/(ld+rd) 这个分数. 如果这个分数在目标区间里, 那么它就是我们的结果; 如果不在, 那么通过判断这个分数在目标区间的左侧还是右侧: 在左侧, 更新区间左端点, 也就是说现在我们的区间是[(ln+rn)/(ld+rd), rn/rd), 类似地在右侧, 更新区间右端点. 重复这个步骤一直到第一种情况, 即所得分数在目标区间中. 这个解法应该是来源于连分数.

好了, 有了思路之后就是程序实现了. 其实上面的思路我在比赛中大致想到了, 然而, 问题是写程序时在判断分数是否在目标区间那一步, 我直接用的是c++的sqrt函数, 这样做精度是不够的(唉又是精度). 而且, 对于类似0.999999999这样的样例, 如果按照上面的步骤去求解也实在是太慢了(循环10^9次, 每次都是加1...). 解决的方法是, 为了避免精度损失, 我们直接用整数运算做判断, 取平方后扔掉更号, 这里需要实现一个简单的高精度 (只需要乘法和比较大小). 而为了加快运算速度, 我们使用二分, 比方说, 如果(ln+rn)/(ld+rd)在目标区间的右侧并且ld和ln很小时, 那么我们就需要找最大的一个k, 使得(ln*k+rn)/(ld*k+rd)在目标区间的右侧, 以此减少上方循环的次数.

嗯其实很多细节我觉得要拿文字解释还是有些麻烦的. 所以直接上代码好了.

#include<bits/stdc++.h>

#define sqr(x) (x)*(x)

typedef long long ll;

const int maxn=55;
const ll mod=1000000007;

ll gcd(ll a,ll b)
{
    if(b==0) return a;
    return gcd(b,a%b);
}

ll m=1000000000;

struct bigint
{
    ll a[8];
    bigint(ll x=0){for(int i=0;i<8;i++) {a[i]=x%m;x/=m;}}
};

bigint operator *(bigint a,bigint b)
{
    bigint c;
    for(int i=0;i<8;i++) for(int j=0;i+j<8;j++) c.a[i+j]+=a.a[i]*b.a[j];
    for(int i=0;i<7;i++) c.a[i+1]+=c.a[i]/m,c.a[i]%=m;
    return c;
}

int cmp(ll a,ll b,ll c,ll d) // comapre c/d and sqrt(a/b)
{
    bigint lh=bigint(c)*bigint(c)*bigint(b);
    bigint rh=bigint(a)*bigint(d)*bigint(d);
    for(int i=7;i>=0;i--)
    {
        if(lh.a[i]!=rh.a[i])
        {
            return lh.a[i]<rh.a[i]?-1:1;
        }
    }
    return 0;
}

int n;
ll a,b,aa,bb;
char c[maxn];

ll ln,rn,ld,rd;

int main()
{
	int i,j,ca=1;
	while(scanf("%s",c)!=EOF)
	{
	    printf("Case #%d: ",ca++);
	    a=0,b=1;
        for(i=2;c[i];i++)
        {
            b*=10;a*=10;a+=c[i]-'0';
        }
        a*=10;b*=10;aa=a;bb=b;a-=5;aa+=5;
        ll d=gcd(a,b);a/=d;b/=d;
        d=gcd(aa,bb);aa/=d;bb/=d;
        ln=0;ld=rn=rd=1;
        for(;;)
        {
            ll c=ln+rn,d=ld+rd;
            //cout<<c<<","<<d<<endl;
            if(cmp(a,b,c,d)==-1)
            {
                if(rd<mod)
                {
                    ll lb=1,ub=mod*2;
                    while(ub-lb>1)
                    {
                        ll mid=(ub+lb)>>1;
                        if(cmp(a,b,ln+mid*rn,ld+mid*rd)>=0)
                        {
                            ub=mid;
                        }
                        else lb=mid;
                    }
                    ln=lb*rn+ln;ld=lb*rd+ld;
                    //printf("%I64d/%I64d\n",ub*rn+ln,ub*rd+ld);
                    //break;
                }
                else {ln=c;ld=d;}
            }
            else if(cmp(aa,bb,c,d)>=0)
            {
                if(ld<mod)
                {
                    ll lb=1,ub=mod*2;
                    while(ub-lb>1)
                    {
                        ll mid=(ub+lb)>>1;
                        if(cmp(aa,bb,ln*mid+rn,ld*mid+rd)<0)
                        {
                            ub=mid;
                        }
                        else lb=mid;
                    }
                    rn=lb*ln+rn;rd=lb*ld+rd;
                }
                else{rn=c;rd=d;}
            }
            else
            {
                printf("%I64d/%I64d\n",c,d);
                break;
            }
        }
	}

    return 0;
}