最少乘法次数【zjoi2005DAY2】

admin

由 开始，通过最少的乘法次数得出
提交文件名：XN.PAS
输入输出格式：
输入文件名：XN.IN。该文件共含10行，每行1个测试数据n（1≤n≤2000）
输出文件名：OUTPUT.TXT。该文件共含10组结果，第i组结果为第i个测试数据对应的乘法方案和最少乘法次数
一组测试数据的输入输出样例：
输入：
XN.IN
  ……
191
……
输出：
OUTPUT.TXT
  ……
X1
X2=X1*X1
X3=X2*X2
X4=X3*X3
X5=X4*X4
X6=X5*X5
X7=X4*X6
X8=X2*X7
X9=X6*X8
X10=X9*X9
X11=X8*X10
X12=X1*X11
Times of multiplies=11
……

admin

因为两式相乘等于方幂相加，所以本题可以等效的表示为：构造一个数列 ，满足

要求 ，并且使 最小。
我们选择回溯法作为本程序的主体结构：当搜索到第 层时， 的取值范围在 到 之间，为了尽快接近目标 ，应当从 开始从大到小为 取值，当然，选取的数都不能大于 。当搜索到 出现时，就得到了一个解，与当前保存的解比较取优。最终搜索结束之后，即得到最终的最优解。
如果按上述结构直接进行搜索，效率显然很低。因此，我们可以考虑使用最优性剪枝进行优化：
优化之一：当然首先要建立估价函数。由于使数列中的最大数加倍无疑是最快的增长方式，所以一种最简单的估价函数为（设数列中当前的最大者是 ，即当前搜索深度为 ）：

然而，这个估价函数的准确性太差，特别是当 大于 时， 只能等于1，根本不能发挥剪枝的作用。因此，无论是深度优先的回溯法还是宽度优先的A*搜索方法，只要还使用这个估价函数，其优化效果都比较有限。
下面，我们再讨论一些进一步的优化手段——
优化之二：着眼于估价函数的“生命”——准确性。我们可以利用加法在奇偶性上的特点的推广，使估价函数在某些情况下的准确性得到一定的提高（具体改进请参见附录）。
优化之三：我们新建立的这个估价函数虽然准确性稍有提高，但时间复杂度也相应提高。为了提高剪枝的效率，我们可以使用一种“逐步细化”的技巧，即先使用一次原先的估价函数（快而不准）进行“过滤”，再使用新的估价函数（稍准但较慢）减少“漏网之鱼”，以求准确性和运行效率的平衡。
优化之四：我们可以在搜索之前将 分解质因数，对每个质因数先使用上述搜索方法求解（如某个质因数仍太大，还可以将其减1，再分解），即相当于将构造数列的过程，划分成若干阶段处理。这样得到的结果虽然不一定是全局的最优解，却可以作为初始设定的优度下界，使后面的全局搜索少走很多弯路。事实上，该估计方法相当准确，在很多情况下都能直接得到最优解，使程序效率得到极大提高。
优化之五：当数列中的当前最大数 超过 后，原估价函数就难以发挥作用了。但是，此后的最优方案，实际上就等价于从 到 这 个数中，选出尽可能少的数（可重复），使它们的和等于 。这个问题已经可以使用动态规划方法来解决。这样得到的“估价函数”不但完全准确，甚至直接就可以代替后面的搜索过程。这里也体现出了搜索和动态规划的优势互补。
测试情况
N值        乘法次数        运行时间
                程序1        程序2        程序3        程序4        程序5        程序6        程序7
127        10         0.22         0.22         0.05         0.06         0.06        <0.01         0.16
191        11         3.52         3.90         0.83         0.55         0.55         0.27         1.75
382        11        13.02        13.29         1.05         0.94         0.66         0.33         2.69
511        12        Very long        Very long         0.82         5.66         0.33         0.27         1.81
635        13        46.96        40.37        20.65        19.61        10.06         8.13        Overflow
719        13        45.59        39.11        45.64        13.45        39.11         6.64        27.35
1002        13        11.81        10.65         2.74         3.46         1.10         0.99         6.26
1023        13        Very long        Very long         2.66        Very long         1.04         0.88         6.15
1357        13         9.51         6.98         5.93         5.44         3.40         2.74        13.35
1894        14        44.87        37.24        13.89        15.49         6.37         4.61        23.68

admin

1. {$A+,B-,D-,E+,F-,G+,I-,L-,N+,O-,P-,Q-,R-,S-,T-,V+,X+,Y-}
   
2. {$M 65520,0,655360}
   
3. program LeastMultiply;  {说明:      1.为了测试计时的方便,本程序从命令行读入n值。      2.程序结束后，在output.txt中给出执行乘法的方式，并给出总的乘法次数。
   
4.      3.在搜索过程中，由于与动态规划结合，所以在没有搜索到底层的时候，就可以得到最优解的数列长度(但此时没有得到完整的最优幂次数列)，所以在搜索结束后，程序中调用formkeep过程在keep中生成完整的构造数列。
   
5.      4.由于程序中的搜索过程生成的是最优幂次数列，而没有直接给出乘法的进行方式，所以在输出结果的过程output中，对其进行了转换。
   
6. 5.为了尽可能的提高程序的时间效率，程序中有几处细节上的优化，请参见程序内的注释。
   
7. 6.程序中，逻辑上相对独立的程序段用空行分开。
   
8. }
   
9. 
   
10. const
    
11.   max=20;{数列最大长度}
    
12.   maxr=2000;{动态规划计数数组的最大长度(输入的n不能超过maxr的2倍)}
    
13.   mp=100;{预处理估价时,可以直接搜索处理的数的范围上限}
    
14.   power2:array[0..12]of integer= (1,2,4,8,16,32,64,128,256,512,1024,2048,4096); {2的方幂}
    
15. 
    
16. type
    
17.   atype=array[0..max]of integer;{用于记录构造的幂次数列的数组类型，0号元素记录数列长度}
    
18. 
    
19. var
    
20.   n:integer;{读入的目标数字}
    
21.   time:array[0..maxr]of integer;
    
22.   {动态规划计数数组，time[i]表示在当前构造数列的基础上，组成数i至少需要的加数个数}
    
23.   range:integer;{使用动态规划处理的范围：range=[n/2]}
    
24.   a:atype;{搜索中记录数列}
    
25.   kp:atype;{预处理估界时，记录结果数列的临时数组}
    
26.   keep:atype;{记录最优结果数列的数组}
    
27.   best:integer;{当前最优解}
    
28.   f:text;{输出文件}
    
29. 
    
30. procedure init;{初始化}
    
31.   var
    
32.     temp:integer;{临时变量}
    
33.   begin
    
34.     val(paramstr(1),n,temp);{从命令行读入n}
    
35.     keep[0]:=1;
    
36.     keep[1]:=1;
    
37.     best:=maxint;{最优数列长度的初值}
    
38.     assign(f,'output.txt');{连接输出文件}
    
39.   end;
    
40. 
    
41. procedure search(n:integer;var best:integer;var keep:atype);{搜索主过程}
    
42. {搜索之前，给出的搜索目标为n；
    
43.            在best中存放搜索前已经给出的优度下界；
    
44.            在keep中存放初始优度下界对应的最优幂次数列。
    
45. 搜索结束之后，在best中给出的是构造的最优幂次数列的长度,即最少乘法次数加1；
    
46.                在keep中给出所构造的最优幂次数列。
    
47. }
    
48.   var
    
49.     kn:integer;{n所含的2的方幂的次数}
    
50.     i:integer;{循环变量}
    
51. 
    
52.   function getk(num:integer):integer;{求num所含的2的方幂次数,即论文中所设的k(num)}
    
53.     var
    
54.       i:integer;{循环变量}
    
55.     begin
    
56.       i:=0;
    
57.       while not odd(num) do begin
    
58.         num:=num shr 1;
    
59.         inc(i);
    
60.       end;
    
61.       getk:=i-1;
    
62.     end;
    
63. 
    
64.   procedure find(step:integer);{递归搜索过程}
    
65.     var
    
66.       i:integer;{循环变量}
    
67.       k:integer;{本层搜索的循环范围上限}
    
68. 
    
69.     function ok(num:integer):boolean;{判断数num能否在当前被构造出来}
    
70.     {为了提高程序的效率，这里利用了动态规划的结果}
    
71.       var
    
72.         i,j:integer;{循环变量}
    
73.       begin
    
74.         if num<=range then begin{待判断数num在[n/2]以内}
    
75.           ok:=(time[num]=2);{直接利用最少需要的加数是否为2来判断}
    
76.           exit;
    
77.         end;
    
78.         for i:=step-1 downto 1 do begin
    
79.           if a[i]+a[i]<num then break;
    
80.           if time[num-a[i]]=1 then begin
    
81.           {time[t]=1表明数t在已有数列中出现过，这样可以避免使用循环判断}
    
82.             ok:=true;
    
83.             exit;
    
84.           end;
    
85.         end;
    
86.         ok:=false;
    
87.       end;
    
88. 
    
89.     procedure evltime;{动态规划子过程}
    
90.       var
    
91.         i,j:integer;{循环变量}
    
92.         p:integer;{本次动态规划递推的上限}
    
93.       begin
    
94.         p:=k;
    
95.         if p>range then p:=range;
    
96.         for i:=a[step-1]+1 to p do begin
    
97.           time[i]:=time[i-a[step-1]]+1;{目标函数赋初值}
    
98.           for j:=step-2 downto 2 do{决策枚举}
    
99.             if time[i-a[j]]+1<time[i] then time[i]:=time[i-a[j]]+1;
    
100.         end;
     
101.       end;
     
102. 
     
103.     function h(num:integer):byte;{最初的简单估价函数h}
     
104.       var
     
105.         i:integer;{循环计数变量}
     
106.       begin
     
107.         i:=0;
     
108.         while num<n do begin
     
109.           num:=num shl 1;
     
110.           inc(i);
     
111.         end;
     
112.         h:=i;
     
113.       end;
     
114. 
     
115.     function cut1:boolean;{最初的剪枝判断，直接调用估价函数h}
     
116.       begin
     
117.         if h(a[step])+step>=best then cut1:=true else cut1:=false;
     
118.       end;
     
119. 
     
120.     function cut2:boolean;{使用改进的估价函数的剪枝判断}
     
121.       var
     
122.         pt:integer;{k(a[step])的值，即a[step]中所含的2的幂的次数}
     
123.         rest:integer;{新构造的数列中，满足k(rest)=k(n)的数}
     
124.         i:integer;{循环计数变量}
     
125.       begin
     
126.         if h(a[step]+a[step-1])+step+1<best then begin
     
127.         {如果新构造数列的第一步选择a[step]+a[step-1]，而不是2*a[step]，就必然导致剪枝，
     
128.          这是新的估价函数起作用的必要条件。
     
129.         }
     
130.           cut2:=false;
     
131.           exit;
     
132.         end;
     
133. 
     
134.         pt:=getk(a[step]);{求k(a[step])}
     
135.         if pt>kn then rest:=a[step-1]
     
136.         else
     
137.           if pt=kn then rest:=a[step]
     
138.           else rest:=maxint;
     
139.              {给rest赋初值，以防新构造的数列中的所有数的k值都小于k(n)}
     
140. 
     
141.         i:=0;
     
142.         repeat
     
143.           a[step+i+1]:=a[step+i] shl 1;
     
144.           inc(i);
     
145.           if pt+i=kn then rest:=a[step+i];
     
146.         until a[step+i]>n;
     
147.         dec(i);
     
148.         {以上程序段为构造数列的过程}
     
149. 
     
150.         if (n-a[step+i]>rest)and(step+i+2>=best) then cut2:=true else cut2:=false;{剪枝判断}
     
151.       end;
     
152. 
     
153.     begin{Find}
     
154.       if a[step-1]+a[step-1]>=n then begin
     
155.       {数列中的当前最大数已经超过[n/2]，则直接引用动态规划的结果}
     
156.         if time[n-a[step-1]]+step-1<best then begin{判断出解}
     
157.           best:=time[n-a[step-1]]+step-1;
     
158.           keep:=a;
     
159.         end;
     
160.         exit;
     
161.       end;
     
162. 
     
163.       k:=a[step-1]+a[step-1];{计算a[step]的可选范围的上限}
     
164.       evltime;{调用动态规划子过程}
     
165. 
     
166.       inc(a[0]);
     
167.       for i:=k downto a[step-1]+1 do
     
168.         if ok(i) then begin
     
169.           if i<=range then time[i]:=1;
     
170.           a[step]:=i;
     
171.           if cut1 then break;{调用剪枝判断1}
     
172.           if cut2 then continue;{调用剪枝判断2}
     
173.           find(step+1);{递归调用下一层搜索}
     
174.         end;
     
175.       dec(a[0]);
     
176.     end;
     
177. 
     
178.   procedure formkeep;{生成最优数列keep}
     
179.   {由于在搜索时，如果a[step]>[n/2]则直接引用动态规划的结果，所以最优结果数列的最后一
     
180.    部分实际上并未得到，所以需要在本过程中将其补充完整。
     
181.    这里还需要使用递归和回溯（实际上就是一个小规模的搜索），不过，由于动态规划给出的结
     
182.    果表示的是从已有数列中，选出最少的数求和而得到n，所以对这些数是可以定序的。
     
183.   }
     
184.     var
     
185.       found:boolean;{找到最优数列的标志}
     
186. 
     
187.     procedure check(from,left:integer);{回溯法生成最优数列的最后一部分}
     
188. {from表示当前层内循环的上界(用于定序)，left表示剩余的需要通过求和而得到的数}
     
189.       var
     
190.         i:integer;{循环变量}
     
191.       begin
     
192.         if keep[0]=best then begin
     
193.           if left=0 then found:=true;{找到最优数列}
     
194.           exit;
     
195.         end;
     
196. 
     
197.         inc(keep[0]);
     
198.         for i:=from downto 1 do{循环枚举数列中的数}
     
199.           if keep[i]<=left then begin
     
200.             keep[keep[0]]:=keep[keep[0]-1]+keep[i];
     
201.             check(i,left-keep[i]);{调用下一层搜索，但所使用的数不能再超过keep[i](定序)}
     
202.             if found then exit;
     
203.           end;
     
204.         dec(keep[0]);
     
205.       end;
     
206. 
     
207.     begin{FromKeep}
     
208.       found:=false;
     
209.       check(keep[0],n-keep[keep[0]]);{调用搜索过程}
     
210.     end;
     
211. 
     
212.   begin{Search}
     
213.     kn:=getk(n);{计算k(n)}
     
214.     time[0]:=0;
     
215.     time[1]:=1;
     
216.     a[0]:=1;
     
217.     a[1]:=1;
     
218.     range:=n div 2;
     
219.     {以上程序段为搜索前的初始化}
     
220. 
     
221.     find(2);{调用搜索过程}
     
222.     formkeep;{搜索结束后，在keep中生成完整的构造数列}
     
223.   end;
     
224. 
     
225. function guess(n:integer):integer;{递归方式实现的预处理估界}
     
226. {说明：
     
227.      1.子程序guess运行结束后，返回的值即为对n估价的结果；同时，在keep数组中得到对应的数列。
     
228.      2.为了能够使估价尽可能准确，guess中同时考虑了两种分解策略，即使用了回溯结构。
     
229.      3.由于每次递归调用guess，其最终结果都必须记入keep数组，所以每次guess运行都只操作keep数组的一部分，同时还要借助于kp,kpt等临时数组。
     
230. }
     
231.   var
     
232.     num:integer;{将n中的所有2的因子都提出后剩下的数}
     
233.     pfact:integer;{表示num的素因子的变量}
     
234.     best:integer;{向下调用时记录返回的估价值}
     
235.     b2:integer;{记录num中的2的因子的数目}
     
236.     s:integer;{对n的估价步数}
     
237.     i:integer;{循环变量}
     
238.     sq:integer;{num的平方根的整数部分，分解素因子时作为上界}
     
239.     s1,k2,k:integer;{临时变量}
     
240.     kpt:atype;{临时记录最优结果数列}
     
241. 
     
242.   begin
     
243.     num:=n;
     
244.     b2:=0;
     
245.     while not odd(num) do begin
     
246.       num:=num div 2;
     
247.       inc(b2);
     
248.       inc(keep[0]);
     
249.       keep[keep[0]]:=power2[b2];
     
250.     end;
     
251.     s:=b2+1;
     
252.     k2:=keep[0];
     
253.     {以上程序段的作用是将n中所含的2的因子全部提出，剩下的部分即num}
     
254. 
     
255.     if num<=mp then begin{如果num足够小(小于等于mp)，则直接调用搜索过程处理}
     
256.       best:=maxint;
     
257.       search(num,best,kp);{对n调用搜索过程，得到的数列存放在kp中}
     
258.       inc(s,best-1);
     
259.       for i:=2 to best do keep[i+keep[0]-1]:=kp[i];
     
260.       inc(keep[0],best-1);
     
261.     end
     
262.     else begin{使用估价方法处理}
     
263.       k:=keep[0];
     
264.       best:=guess(num-1);{递归调用guess}
     
265.       keep[keep[0]+1]:=keep[keep[0]]+1;
     
266.       inc(keep[0]);
     
267.       inc(s,best);
     
268.       {以上程序段为估价的第一种策略：将num减1后，对num-1进行估价}
     
269. 
     
270.       sq:=trunc(sqrt(num));
     
271.       pfact:=3;
     
272.       while (pfact<=sq)and(num mod pfact>0) do inc(pfact,2);
     
273.       if pfact<=sq then begin
     
274.         kpt:=keep;
     
275.         keep[0]:=k2;
     
276.         s1:=b2+1;
     
277.         {以上程序段为使用第二种策略前的初始化}
     
278. 
     
279.         k:=keep[0];
     
280.         best:=guess(pfact);
     
281.         inc(s1,best-1);
     
282.         {以上程序段中递归调用guess，对质因数pfact进行估价}
     
283. 
     
284.         k:=keep[0];
     
285.         best:=guess(num div pfact);
     
286.         for i:=k+1 to keep[0] do keep[i]:=pfact*keep[i];
     
287.         inc(s1,best-1);
     
288.         {以上程序段对[num/pfact]进行估价}
     
289. 
     
290.         if s1<s then s:=s1 else keep:=kpt;{比较两种策略结果的优劣}
     
291.       end;
     
292.       {以上程序段是估价的第二种策略：将num分解出一个较小的质因数后再处理。
     
293.        估价的结果存放在s1中，使用kpt临时存放第一种策略所构造的数列。
     
294.       }
     
295.     end;
     
296. 
     
297.     for i:=k2+1 to keep[0] do keep[i]:=keep[i]*power2[b2];
     
298.     guess:=s;{返回估价结果}
     
299.   end;
     
300. 
     
301. procedure output;{输出结果}
     
302.   var
     
303.     i,j,k:integer;{循环变量}
     
304.   begin
     
305.     rewrite(f);
     
306.     writeln(f,'X1');
     
307.     for i:=2 to best do begin
     
308.       for j:=1 to i-1 do begin
     
309.         for k:=j to i-1 do
     
310.           if keep[j]+keep[k]=keep[i] then break;
     
311.         if keep[j]+keep[k]=keep[i] then break;
     
312.       end;
     
313.       writeln(f,'X',i,'=X',j,'*X',k);
     
314.     end;
     
315.     writeln(f,'Times of multiplies=',best-1);
     
316.     close(f);
     
317.   end;
     
318. 
     
319. begin{主程序}
     
320.   init;{读入数据，初始化}
     
321.   best:=guess(n);{调用预处理估界}
     
322.   search(n,best,keep);{调用搜索主过程}
     
323.   output;{输出结果}
     
324. end.
     

复制代码