用兩個二分查詢實現的海倫信封程式

lt發表於2017-01-01

原問題採用sqrt函式判斷完全平方數。
參考帖子採用字典查詢。
在結構中新增一個成員，用來儲存子序列長度。第1個二分查詢第一個成員，第2個查詢第2個成員。

#include<cstdio>
#include<cmath>
#include<cstdlib>
#define MAXP 10000000
int gcd(int a,int b)
{
    int a1,b1,r1;
    if(a==1||b==1) return 1;
    a1=a;
    b1=b;
    r1=a1%b1;
    while(r1)
    {
        a1=b1;
        b1=r1;
        r1=a1%b1;
    }
    return b1;
}

typedef struct
{
    int A=0;
    int B=0;
    int dis=0;
} ST ;
static ST st[MAXP*2];
int cnt=0;

//sort
int init()
{
    for (int a = 2; ; a++)
    {
        int done = 1;
        for (int x, y, b = 1; b < a; b++)
        {
            if ((x = a*b << 1) + (y = a*a - b*b) >= MAXP/2) break;
            if ((1 & (a + b)) == 0 || gcd(a, b) != 1) continue;
            for (int u = x, v = y; u + v < MAXP/2; u += x, v += y)
                //{ Add(dict, u, v); Add(dict, v, u); }
            {
                st[cnt].A=u;
                st[cnt].B=v;
                cnt++;
                //swap
                st[cnt].A=v;
                st[cnt].B=u;
                cnt++;
            }
            done = 0;
        }
        if (done==1) break;
    }
    printf("cnt=%d\n",cnt);
    return 0;
}

int cmp2 ( const void *a, const void *b )
{

    ST *ia = (ST *)a; //asc order
    ST *ib = (ST *)b; //asc order
    long long x=(ia->A - ib->A);
    if (x>0)
        return 1;
    else if(x<0)
        return -1;
    else //(x==0)
        return (ia->B > ib->B)?1:-1;

}
int sort_ggs(int cnt)
{
    qsort(st,cnt,sizeof(st[0]),cmp2);
    return 0;
}
int bsech(int a)
{

    int lo=0,hi=cnt-1;
    while(lo<=hi )
    {
        int mid= (lo+hi)/2;
        if(st[mid].A<a)
        {
            lo=mid+1;
        }
        else if (st[mid].A>a)
        {
            hi=mid -1;
        }
        else
        {
            return mid;
        }
    }
    return -1;


}
int bsech2(int low,int high,int b)
{

    size_t lo=low,hi=high;
    while(lo<=hi)
    {
        size_t mid= (lo+hi)/2;
        if(st[mid].B<b)
        {
            lo=mid+1;
        }
        else if (st[mid].B>b)
        {
            hi=mid -1;
        }
        else
        {
            return mid;
        }
    }
    return -1;
}
int main()
{
    long long a,b,c,d,e,b_2,p;
    int c1=0;
    long long h=0;
    long long sum=0;
    long long a0,b0,a1,b1;
    init();
    sort_ggs(cnt);
    int dis=1;
    for(int i=0; i<cnt-1; i++)
    {
        if(st[i].A==st[i+1].A)
        {
            st[i].dis=dis;
            dis++;
        }
        else
        {
            st[i].dis=dis;
            st[i-dis+1].dis=-dis;
            dis=1;
        }
    }
    if(st[cnt-1].A!=st[cnt-2].A)st[cnt-1].dis=-1;

    int rep=1;
    #pragma omp parallel for reduction(+:sum,c1, a0,b0,a1,b1,a,b,b_2,h,p)
    for(int i=0; i<cnt; i++)
    {
        a0=st[i].A;
        b0=st[i].B;
        if(b0%2==1)
            continue;

        b_2=b0/2;

        int po=bsech(b_2);

        if (po==-1)
            continue;

        //for(; po>0 && st[po].A==b_2; po--); //bsech may return mid,get first match
        int startpos=0;
        int len=0;
        if(st[po].dis<0)
            startpos=po;
        else
            startpos=po-st[po].dis+1;

        len=-st[startpos].dis;

        if(len<=1)continue;
        //if(po==cnt-2 || st[po+1].A!=st[po+2].A) continue;
        for(int j=startpos; j<startpos+len/*j<cnt && st[j].A==b_2&& (st[j].tg>st[i].tg*2.0)*/; j++)
        {
            b1=st[j].A; //swap a b
            a1=st[j].B; //h

            b=b0;//lcm(b0,b1);

            if(b_2*2 +(a0+a1)*2 >=MAXP) //b+ 2a+2h
                break; //a1 asc order
            //continue;

            h=a1; // 29.309
            //h=a0; //33.441
            a=a0; // 29.309

            if(h<a)
                //for(int j2=j+1; j2<cnt && st[j2].A==b_2; j2++) //search a+h
                //for(int j2=j+1; j2<startpos+len; j2++) //search a+h
                    //if(st[j2].B==a+h)//found
                         if(bsech2(j+1,startpos+len-1,a+h)!=-1)
                        // if(fabs(round(sqrt(a*a+4*b_2*b_2))-sqrt(a*a+4*b_2*b_2)) <1e-7)
                    {
                        if((p=(2*a+2*b_2+2*round(sqrt(h*h+b_2*b_2))))<=MAXP)
                        {
                            //int n=MAXP/p;
                            sum+=p;//p*n*(1+n)/2;
                            //printf("%lld %lld %lld %lld %lld\n",a,b,b_2,h,p);
                            c1++;//=n;
                        }
                    }
        }

    }
    printf("%lld %d\n",sum,c1);
    return 0;
}

D:\>g++ p583tg9.cpp -O2

D:\>timer a
Timer 9.01 : Igor Pavlov : Public domain : 2009-05-31
cnt=17205400
1174137929000 233887

Kernel Time  =     0.405 =    1%
User Time    =    34.257 =   97%
Process Time =    34.663 =   98%
Global Time  =    35.235 =  100%

D:\>g++ p583tg9.cpp -O2 -march=native -fopenmp -lm -o ap.exe

D:\>timer ap
Timer 9.01 : Igor Pavlov : Public domain : 2009-05-31
cnt=17205400
1174137929000 233887

Kernel Time  =     0.296 =    1%
User Time    =    45.209 =  206%
Process Time =    45.505 =  207%
Global Time  =    21.888 =  100%

相關文章