PKU1945 Power Hungry Cows

【题目大意】

开始时有两个变量（A=x^a,B=x^b），每次可以做以下爱几种运算：

A*B，A*A，B*B，A/B，B/A，A/A，B/B

并把结果放到某个变量中，问最少多少步可以得到目标X^P

【题目分析】

我们可以只记录当前两个变量的指数，这样各种操作就变成了加和减，可以直接进行BFS，但是显然会超时，于是我们先从状态入手，规定a<=b，这样就把状态数减少了一半，而操作中的除法也就只剩下b/a这一种情况了，再继续考虑优化。

对于各种操作，其实A*A这个变化时与一些操作重复的，所以可以去掉这个变换，这样就减少了搜索树结点的度，使能达到的节点数更少了。

状态树上已经不能再优化了，继续考虑较少扩展量，有如下两个剪枝

①如果p mod gcd(A,B)<>0 那么状态<A,B>一定必能出解

②如果A>=100且B>=50000那么该状态无解

以上两个剪枝笔者都不会证明。。。。

以上两个剪枝加入后程序性能有所改善，但是仍然不能在时间限制内出解。我们只能考虑搜索扩展规则上的优化，使用A*算法。

 

介绍一下A*算法，给每一个状态一个估价值，表示该状态达到目标的期望水平，函数值越大，则该节点越有可能达到解，因此，我们在进行A*算法时需要维护两个表，表一代替了朴素BFS中的队列，能够在较短时间内取出并删除函数值最大的结点，还能在较短时间内插入新节点，而第二个表则是用来判重的，能够支持插入节点和查询节点是否已经存在，对于数据结构了解比较清楚的同学可能已经知道他们应该用什么实现了，表一用一个二叉堆来实现，在logn时间内完成插入操作或取出最小值并维护形态，表二则用一个平衡树或hash表来实现，个人推荐平衡树，比较稳定，不会像hash那样出问题。这样，A*的数据结构就选好了。

接下来确定的是估价函数，最快情况下变量按2的指数级别增长，所以定义H(state)为变量B不断自乘2最少多少次能大于等于P，定义G(state)为该节点的深度，定义估价函数F(state)=G(state)*2+H(state)

这里给G(state)乘2是为了保证估价函数的相容性，即保证先扩展出来的结点的步数一定是最优的，这样才能得到最优解。同时，利用刚刚的朴素BFS程序，我们还可以得到数据范围内的最大步数大概在20左右，所以如果某结点函数F(state)>20*2，那么我们完全可以舍弃这个节点。

如果还是不能出解，我们就要考虑牺牲解的正确性了，我的做法是在数据较强时估价函数G不乘2，这样时间大幅下降。

这样加入所有的优化后，程序效率明显上升，能够在较短时间内出解。而且结点总数在20万以下，可以说已经优化的很好了。

 

{$inline on}
program usaco_power(input,output);
const
    maxpower= 200000;
type
    nodes=record
        x,y:longint;
    end;
    size_balance_tree=object
        root,tot:longint;
        s,left,right:array[0..200000] of longint;
        key:array[0..300000] of nodes;
        procedure clean;
        procedure left_rotate(var t:longint);
        procedure right_rotate(var t:longint);
        procedure maintain(var t:longint;flag:boolean);
        procedure insert(var t:longint;xx:nodes);
        function find(t:longint;xx:nodes):boolean;
    end;
    node=record
        x,y,w,step:longint;
    end;
    Theap=object
        tot:longint;
        list:array[0..100000] of node;
        procedure clean;
        procedure swap(var aa,bb:longint);
        procedure up(now:longint);
        procedure down(now:longint);
        procedure insert(now:node);
        function delete():node;
    end;
procedure size_balance_tree.clean(); 
begin
    root:=0;
    tot:=0;
    fillchar(left,sizeof(left),0);
    fillchar(right,sizeof(right),0);
    fillchar(s,sizeof(s),0);
    fillchar(key,sizeof(key),0);
end;
procedure size_balance_tree.left_rotate(var t:longint); 
var
    k:longint;
begin
    k:=right[t];
    right[t]:=left[k];
    left[k]:=t;
    s[k]:=s[t];
    s[t]:=s[left[t]]+s[right[t]]+1;
    t:=k;
end;
procedure size_balance_tree.right_rotate(var t:longint); 
var
    k:longint;
begin
    k:=left[t];
    left[t]:=right[k];
    right[k]:=t;
    s[k]:=s[t];
    s[t]:=s[left[t]]+s[right[t]]+1;
    t:=k;
end;
procedure size_balance_tree.maintain(var t:longint;flag:boolean);
begin
    if not flag then
    begin
        if s[left[left[t]]]>s[right[t]] then
            size_balance_tree.right_rotate(t)
        else
            if s[right[left[t]]]>s[right[t]] then
            begin
                size_balance_tree.left_rotate(left[t]);
                size_balance_tree.right_rotate(t);
            end
            else
                exit;
    end
    else
    begin
        if s[right[right[t]]]>s[left[t]] then
            size_balance_tree.left_rotate(t)
        else
            if s[left[right[t]]]>s[left[t]] then
            begin
                size_balance_tree.right_rotate(right[t]);
                size_balance_tree.left_rotate(t);
            end
            else
                exit;
    end;
    size_balance_tree.maintain(left[t],false);
    size_balance_tree.maintain(right[t],true);
    size_balance_tree.maintain(t,false);
    size_balance_tree.maintain(t,true);
end;
procedure size_balance_tree.insert(var t:longint;xx:nodes); 
begin
    if t=0 then
    begin
        inc(tot);
        t:=tot;
        left[t]:=0;
        right[t]:=0;
        s[t]:=1;
        key[t]:=xx;
    end
    else
    begin
        inc(s[t]);
        if (xx.x*maxpower+xx.y)<(key[t].x*maxpower+key[t].y) then
            size_balance_tree.insert(left[t],xx)
        else
            size_balance_tree.insert(right[t],xx);
        size_balance_tree.maintain(t,xx.x*maxpower+xx.y>=key[t].x*maxpower+key[t].y);
    end;
end;
function size_balance_tree.find(t:longint;xx:nodes):boolean; 
begin
    if t=0 then
        exit(false);
    if (key[t].x=xx.x)and(key[t].y=xx.y) then
        exit(true);
    if (xx.x*maxpower+xx.y)<(key[t].x*maxpower+key[t].y) then
        exit(size_balance_tree.find(left[t],xx))
    else
        exit(size_balance_tree.find(right[t],xx));
end;


procedure Theap.clean(); 
begin
    fillchar(list,sizeof(list),0);
    tot:=0;
end;
procedure Theap.swap(var aa,bb:longint); 
var
    tt:longint;
begin
    tt:=aa;
    aa:=bb;
    bb:=tt;
end;
procedure Theap.up(now:longint); 
begin
    while (now>1)and(list[now].w<list[now>>1].w) do
    begin
        Theap.swap(list[now].x,list[now>>1].x);
        Theap.swap(list[now].y,list[now>>1].y);
        Theap.swap(list[now].w,list[now>>1].w);
        Theap.swap(list[now].step,list[now>>1].step);
        now:=now>>1;
    end;
end;
procedure Theap.down(now:longint); 
var
    tmp:longint;
begin
    while (now*2<=tot) do
    begin
        tmp:=now*2;
        if (tmp+1<=tot)and(list[tmp+1].w<list[tmp].w) then
            inc(tmp);
        if list[tmp].w<list[now].w then
        begin
            Theap.swap(list[now].x,list[tmp].x);
            Theap.swap(list[now].y,list[tmp].y);
            Theap.swap(list[now].w,list[tmp].w);
            Theap.swap(list[now].step,list[tmp].step);
            now:=tmp;
        end
        else
            break;
    end;
end;
procedure Theap.insert(now:node); 
begin
    inc(tot);
    list[tot]:=now;
    Theap.up(tot);
end;
function Theap.delete():node; 
begin
    delete:=list[1];
    Theap.swap(list[1].x,list[tot].x);
    Theap.swap(list[1].y,list[tot].y);
    Theap.swap(list[1].w,list[tot].w);
    Theap.swap(list[1].step,list[tot].step);    
    dec(tot);
    Theap.down(1);
end;

var
    tree:size_balance_tree;
    heap:Theap;
    goal,answer:longint;
    answer_flag,check_goal:boolean;
function gcd(a,b:longint):longint; inline;
begin
    if b=0 then
        exit(a);
    exit(gcd(b,a mod b));
end;{ gcd }
procedure check();
var
    i,s:longint;
begin;
    s:=0;
    for i:=0 to 30 do
        if ((1<<i) and goal)>0 then
            inc(s);
    if s<=9 then
        check_goal:=true
    else
        check_goal:=false;
end;
function f(now:node):longint; inline;
var
    tmp:longint;
begin
    f:=0;
    tmp:=now.y;
    if tmp=0 then
        exit(maxlongint>>1);
    while tmp<goal do
    begin
        tmp:=tmp+tmp;
        inc(f);
    end;
    if (check_goal)and(goal>10000) then
        f:=f+now.step
    else
        f:=f+now.step<<1;
end;{ f }
procedure swap(var aa,bb:longint); inline;
var
    tt:longint;
begin
    tt:=aa;
    aa:=bb;
    bb:=tt;
end;{ swap }
procedure init; inline;
var
    now:node;
    new:nodes;
begin
    readln(goal);
    check();
    tree.clean();
    heap.clean();
    now.x:=0;
    now.y:=1;
    now.step:=0;
    now.w:=f(now);
    new.x:=0;
    new.y:=1;
    tree.insert(tree.root,new);
    heap.insert(now);
end;{ init }
function operand(now:node;flag:longint):node; inline;
begin
    case flag of
        1:begin
            operand.x:=now.x;
            operand.y:=now.x+now.y;
            operand.step:=now.step+1;
        end;
        2:begin
            operand.x:=now.x+now.y;
            operand.y:=now.y;
            operand.step:=now.step+1;
        end;
        3:begin
            operand.x:=now.y-now.x;
            operand.y:=now.y;
            operand.step:=now.step+1;
        end;
        4:begin
            operand.x:=now.x;
            operand.y:=now.y-now.x;
            operand.step:=now.step+1;
        end;
        {5:begin
            operand.x:=now.x+now.x;
            operand.y:=now.y;
            operand.step:=now.step+1;
        end;
        6:begin
            operand.x:=now.x;
            operand.y:=now.x+now.x;
            operand.step:=now.step+1;
        end;}
        5:begin
            operand.x:=now.y+now.y;
            operand.y:=now.y;
            operand.step:=now.step+1;
        end;
        6:begin
            operand.x:=now.x;
            operand.y:=now.y+now.y;
            operand.step:=now.step+1;
        end;
    end;
    if operand.y<operand.x then
        swap(operand.x,operand.y);
end;{ operand }
procedure expand(now:node); inline;
var
    newstate:node;
    pd:nodes;
    i:longint;
begin
    for i:=1 to 6 do
    begin
        newstate:=operand(now,i);
        if not ((newstate.x<=100)and(newstate.y<=50000)) then
            continue;
        if newstate.y=0 then
            continue;
        if (goal mod gcd(newstate.x,newstate.y))<>0 then
            continue;
        if (newstate.x=goal)or(newstate.y=goal) then
        begin
            answer_flag:=true;
            answer:=newstate.step;
            exit;
        end;
        newstate.w:=f(newstate);
        if newstate.w>40 then
            continue;
        pd.x:=newstate.x;
        pd.y:=newstate.y;
        if tree.find(tree.root,pd) then
            continue;
        tree.insert(tree.root,pd);
        heap.insert(newstate);
    end;
end;{ expand }
procedure A_star; inline;
var
    state:node;
begin
    answer_flag:=false;
    while heap.tot<>0 do
    begin
        state:=heap.delete();
        expand(state);
        if answer_flag then
            break;
    end;
end;{ A_star }
procedure print; inline;
begin
    writeln(answer);
end;{ print }
begin
    init;
    A_star;
    print;
end.