最少乘法次数【zjoi2005DAY2】 - 华师一附中OI组

{$A+,B-,D-,E+,F-,G+,I-,L-,N+,O-,P-,Q-,R-,S-,T-,V+,X+,Y-}
{$M 65520,0,655360}
program LeastMultiply; {说明: 1.为了测试计时的方便,本程序从命令行读入n值。 2.程序结束后，在output.txt中给出执行乘法的方式，并给出总的乘法次数。
3.在搜索过程中，由于与动态规划结合，所以在没有搜索到底层的时候，就可以得到最优解的数列长度(但此时没有得到完整的最优幂次数列)，所以在搜索结束后，程序中调用formkeep过程在keep中生成完整的构造数列。
4.由于程序中的搜索过程生成的是最优幂次数列，而没有直接给出乘法的进行方式，所以在输出结果的过程output中，对其进行了转换。
5.为了尽可能的提高程序的时间效率，程序中有几处细节上的优化，请参见程序内的注释。
6.程序中，逻辑上相对独立的程序段用空行分开。
}
const
max=20;{数列最大长度}
maxr=2000;{动态规划计数数组的最大长度(输入的n不能超过maxr的2倍)}
mp=100;{预处理估价时,可以直接搜索处理的数的范围上限}
power2:array[0..12]of integer= (1,2,4,8,16,32,64,128,256,512,1024,2048,4096); {2的方幂}
type
atype=array[0..max]of integer;{用于记录构造的幂次数列的数组类型，0号元素记录数列长度}
var
n:integer;{读入的目标数字}
time:array[0..maxr]of integer;
{动态规划计数数组，time[i]表示在当前构造数列的基础上，组成数i至少需要的加数个数}
range:integer;{使用动态规划处理的范围：range=[n/2]}
a:atype;{搜索中记录数列}
kp:atype;{预处理估界时，记录结果数列的临时数组}
keep:atype;{记录最优结果数列的数组}
best:integer;{当前最优解}
f:text;{输出文件}
procedure init;{初始化}
var
temp:integer;{临时变量}
begin
val(paramstr(1),n,temp);{从命令行读入n}
keep[0]:=1;
keep[1]:=1;
best:=maxint;{最优数列长度的初值}
assign(f,'output.txt');{连接输出文件}
end;
procedure search(n:integer;var best:integer;var keep:atype);{搜索主过程}
{搜索之前，给出的搜索目标为n；
在best中存放搜索前已经给出的优度下界；
在keep中存放初始优度下界对应的最优幂次数列。
搜索结束之后，在best中给出的是构造的最优幂次数列的长度,即最少乘法次数加1；
在keep中给出所构造的最优幂次数列。
}
var
kn:integer;{n所含的2的方幂的次数}
i:integer;{循环变量}
function getk(num:integer):integer;{求num所含的2的方幂次数,即论文中所设的k(num)}
var
i:integer;{循环变量}
begin
i:=0;
while not odd(num) do begin
num:=num shr 1;
inc(i);
end;
getk:=i-1;
end;
procedure find(step:integer);{递归搜索过程}
var
i:integer;{循环变量}
k:integer;{本层搜索的循环范围上限}
function ok(num:integer):boolean;{判断数num能否在当前被构造出来}
{为了提高程序的效率，这里利用了动态规划的结果}
var
i,j:integer;{循环变量}
begin
if num<=range then begin{待判断数num在[n/2]以内}
ok:=(time[num]=2);{直接利用最少需要的加数是否为2来判断}
exit;
end;
for i:=step-1 downto 1 do begin
if a[i]+a[i]<num then break;
if time[num-a[i]]=1 then begin
{time[t]=1表明数t在已有数列中出现过，这样可以避免使用循环判断}
ok:=true;
exit;
end;
end;
ok:=false;
end;
procedure evltime;{动态规划子过程}
var
i,j:integer;{循环变量}
p:integer;{本次动态规划递推的上限}
begin
p:=k;
if p>range then p:=range;
for i:=a[step-1]+1 to p do begin
time[i]:=time[i-a[step-1]]+1;{目标函数赋初值}
for j:=step-2 downto 2 do{决策枚举}
if time[i-a[j]]+1<time[i] then time[i]:=time[i-a[j]]+1;
end;
end;
function h(num:integer):byte;{最初的简单估价函数h}
var
i:integer;{循环计数变量}
begin
i:=0;
while num<n do begin
num:=num shl 1;
inc(i);
end;
h:=i;
end;
function cut1:boolean;{最初的剪枝判断，直接调用估价函数h}
begin
if h(a[step])+step>=best then cut1:=true else cut1:=false;
end;
function cut2:boolean;{使用改进的估价函数的剪枝判断}
var
pt:integer;{k(a[step])的值，即a[step]中所含的2的幂的次数}
rest:integer;{新构造的数列中，满足k(rest)=k(n)的数}
i:integer;{循环计数变量}
begin
if h(a[step]+a[step-1])+step+1<best then begin
{如果新构造数列的第一步选择a[step]+a[step-1]，而不是2*a[step]，就必然导致剪枝，
这是新的估价函数起作用的必要条件。
}
cut2:=false;
exit;
end;
pt:=getk(a[step]);{求k(a[step])}
if pt>kn then rest:=a[step-1]
else
if pt=kn then rest:=a[step]
else rest:=maxint;
{给rest赋初值，以防新构造的数列中的所有数的k值都小于k(n)}
i:=0;
repeat
a[step+i+1]:=a[step+i] shl 1;
inc(i);
if pt+i=kn then rest:=a[step+i];
until a[step+i]>n;
dec(i);
{以上程序段为构造数列的过程}
if (n-a[step+i]>rest)and(step+i+2>=best) then cut2:=true else cut2:=false;{剪枝判断}
end;
begin{Find}
if a[step-1]+a[step-1]>=n then begin
{数列中的当前最大数已经超过[n/2]，则直接引用动态规划的结果}
if time[n-a[step-1]]+step-1<best then begin{判断出解}
best:=time[n-a[step-1]]+step-1;
keep:=a;
end;
exit;
end;
k:=a[step-1]+a[step-1];{计算a[step]的可选范围的上限}
evltime;{调用动态规划子过程}
inc(a[0]);
for i:=k downto a[step-1]+1 do
if ok(i) then begin
if i<=range then time[i]:=1;
a[step]:=i;
if cut1 then break;{调用剪枝判断1}
if cut2 then continue;{调用剪枝判断2}
find(step+1);{递归调用下一层搜索}
end;
dec(a[0]);
end;
procedure formkeep;{生成最优数列keep}
{由于在搜索时，如果a[step]>[n/2]则直接引用动态规划的结果，所以最优结果数列的最后一
部分实际上并未得到，所以需要在本过程中将其补充完整。
这里还需要使用递归和回溯（实际上就是一个小规模的搜索），不过，由于动态规划给出的结
果表示的是从已有数列中，选出最少的数求和而得到n，所以对这些数是可以定序的。
}
var
found:boolean;{找到最优数列的标志}
procedure check(from,left:integer);{回溯法生成最优数列的最后一部分}
{from表示当前层内循环的上界(用于定序)，left表示剩余的需要通过求和而得到的数}
var
i:integer;{循环变量}
begin
if keep[0]=best then begin
if left=0 then found:=true;{找到最优数列}
exit;
end;
inc(keep[0]);
for i:=from downto 1 do{循环枚举数列中的数}
if keep[i]<=left then begin
keep[keep[0]]:=keep[keep[0]-1]+keep[i];
check(i,left-keep[i]);{调用下一层搜索，但所使用的数不能再超过keep[i](定序)}
if found then exit;
end;
dec(keep[0]);
end;
begin{FromKeep}
found:=false;
check(keep[0],n-keep[keep[0]]);{调用搜索过程}
end;
begin{Search}
kn:=getk(n);{计算k(n)}
time[0]:=0;
time[1]:=1;
a[0]:=1;
a[1]:=1;
range:=n div 2;
{以上程序段为搜索前的初始化}
find(2);{调用搜索过程}
formkeep;{搜索结束后，在keep中生成完整的构造数列}
end;
function guess(n:integer):integer;{递归方式实现的预处理估界}
{说明：
1.子程序guess运行结束后，返回的值即为对n估价的结果；同时，在keep数组中得到对应的数列。
2.为了能够使估价尽可能准确，guess中同时考虑了两种分解策略，即使用了回溯结构。
3.由于每次递归调用guess，其最终结果都必须记入keep数组，所以每次guess运行都只操作keep数组的一部分，同时还要借助于kp,kpt等临时数组。
}
var
num:integer;{将n中的所有2的因子都提出后剩下的数}
pfact:integer;{表示num的素因子的变量}
best:integer;{向下调用时记录返回的估价值}
b2:integer;{记录num中的2的因子的数目}
s:integer;{对n的估价步数}
i:integer;{循环变量}
sq:integer;{num的平方根的整数部分，分解素因子时作为上界}
s1,k2,k:integer;{临时变量}
kpt:atype;{临时记录最优结果数列}
begin
num:=n;
b2:=0;
while not odd(num) do begin
num:=num div 2;
inc(b2);
inc(keep[0]);
keep[keep[0]]:=power2[b2];
end;
s:=b2+1;
k2:=keep[0];
{以上程序段的作用是将n中所含的2的因子全部提出，剩下的部分即num}
if num<=mp then begin{如果num足够小(小于等于mp)，则直接调用搜索过程处理}
best:=maxint;
search(num,best,kp);{对n调用搜索过程，得到的数列存放在kp中}
inc(s,best-1);
for i:=2 to best do keep[i+keep[0]-1]:=kp[i];
inc(keep[0],best-1);
end
else begin{使用估价方法处理}
k:=keep[0];
best:=guess(num-1);{递归调用guess}
keep[keep[0]+1]:=keep[keep[0]]+1;
inc(keep[0]);
inc(s,best);
{以上程序段为估价的第一种策略：将num减1后，对num-1进行估价}
sq:=trunc(sqrt(num));
pfact:=3;
while (pfact<=sq)and(num mod pfact>0) do inc(pfact,2);
if pfact<=sq then begin
kpt:=keep;
keep[0]:=k2;
s1:=b2+1;
{以上程序段为使用第二种策略前的初始化}
k:=keep[0];
best:=guess(pfact);
inc(s1,best-1);
{以上程序段中递归调用guess，对质因数pfact进行估价}
k:=keep[0];
best:=guess(num div pfact);
for i:=k+1 to keep[0] do keep[i]:=pfact*keep[i];
inc(s1,best-1);
{以上程序段对[num/pfact]进行估价}
if s1<s then s:=s1 else keep:=kpt;{比较两种策略结果的优劣}
end;
{以上程序段是估价的第二种策略：将num分解出一个较小的质因数后再处理。
估价的结果存放在s1中，使用kpt临时存放第一种策略所构造的数列。
}
end;
for i:=k2+1 to keep[0] do keep[i]:=keep[i]*power2[b2];
guess:=s;{返回估价结果}
end;
procedure output;{输出结果}
var
i,j,k:integer;{循环变量}
begin
rewrite(f);
writeln(f,'X1');
for i:=2 to best do begin
for j:=1 to i-1 do begin
for k:=j to i-1 do
if keep[j]+keep[k]=keep[i] then break;
if keep[j]+keep[k]=keep[i] then break;
end;
writeln(f,'X',i,'=X',j,'*X',k);
end;
writeln(f,'Times of multiplies=',best-1);
close(f);
end;
begin{主程序}
init;{读入数据，初始化}
best:=guess(n);{调用预处理估界}
search(n,best,keep);{调用搜索主过程}
output;{输出结果}
end.

复制代码