Pro_3.pas 3.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163
  1. {高精度计算问题,这里给出的是算法3-3}
  2. {$R-,S-,Q-}
  3. {$M 65520,0,655360}
  4. const
  5. Maxn =9999;{N的最大值}
  6. MaxK =99;{K的最大值}
  7. MaxValue =10000;
  8. Maxkk =(MaxK-1) div 4+1;
  9. type
  10. TNum =record{高精度类型}
  11. l :Integer;
  12. x :array[1..Maxkk+1] of Integer
  13. end;
  14. var
  15. Save :array[1..Maxn] of ^TNum;{保存计算过的信息}
  16. Split :array[1..Maxn] of Integer;{保存合数的分解}
  17. Result,Org :TNum;
  18. Time :LongInt;
  19. n,m,k,kk :Integer;
  20. Template :Longint;
  21. procedure InitSplit;{对合数进行分解}
  22. var
  23. i,j :Integer;
  24. begin
  25. FillChar(Split,Sizeof(Split),0);
  26. Split[1]:=1;
  27. for i:=2 to n do
  28. if Split[i]=0 then begin
  29. Split[i]:=1;
  30. j:=i+i;while j<=n do begin Split[j]:=i;Inc(j,i) end
  31. end
  32. end;
  33. procedure Mul(var x:TNum;y:TNum);{x:=x*y}
  34. var
  35. i,j :Integer;
  36. p,q,Tmp :LongInt;
  37. begin
  38. Org:=x;FillChar(x,Sizeof(x),0);
  39. for i:=1 to y.l do begin
  40. p:=y.x[i];
  41. if p<>0 then begin
  42. Tmp:=0;
  43. for j:=i to i+Org.l-1 do begin
  44. Inc(Tmp,p*Org.x[j-i+1]);
  45. q:=x.x[j];
  46. Inc(q,Tmp);
  47. Tmp:=q div MaxValue;
  48. Dec(q,Tmp*MaxValue);
  49. x.x[j]:=q
  50. end;
  51. Inc(x.x[j+1],Tmp)
  52. end;
  53. if Org.l+i=kk+1 then Dec(Org.l)
  54. end;
  55. x.l:=y.l+Org.l;if x.x[x.l]=0 then Dec(x.l)
  56. end;
  57. procedure Add(var x,y:TNum);{x:=x+y}
  58. var
  59. i :Integer;
  60. begin
  61. if x.l<y.l then x.l:=y.l;
  62. for i:=1 to y.l do begin
  63. Inc(x.x[i],y.x[i]);
  64. if x.x[i]>=MaxValue then begin
  65. Dec(x.x[i],MaxValue);
  66. Inc(x.x[i+1])
  67. end
  68. end;
  69. while (i<=x.l) and (x.x[i]>=MaxValue) do begin
  70. Dec(x.x[i],MaxValue);
  71. Inc(i);
  72. Inc(x.x[i])
  73. end;
  74. if (x.x[x.l+1]=1) and (x.l<kk) then Inc(x.l)
  75. end;
  76. procedure Power(var Res,x:TNum;m:Integer);{Res=x^m}
  77. begin
  78. if m=1 then Res:=x
  79. else begin
  80. Power(Res,x,m shr 1);
  81. Mul(Res,Res);
  82. if Odd(m) then Mul(Res,x)
  83. end
  84. end;
  85. function IntToStr(x:LongInt):string;
  86. var
  87. Res :string;
  88. begin
  89. Str(x,Res);
  90. if x<10 then IntToStr:='000'+Res else
  91. if x<100 then IntToStr:='00'+Res else
  92. if x<1000 then IntToStr:='0'+Res else IntToStr:=Res
  93. end;
  94. procedure Print;{输出}
  95. var
  96. i :Byte;
  97. Res :string;
  98. begin
  99. Res:='';
  100. for i:=Result.l downto 1 do Res:=Res+IntToStr(Result.x[i]);
  101. while Length(Res)>k do Delete(Res,1,1);
  102. while Res[1]='0' do Delete(Res,1,1);
  103. Writeln(Res)
  104. end;
  105. procedure Main;{计算主过程}
  106. var
  107. i :Integer;
  108. x,y :TNum;
  109. MustSave :array[1..Maxn] of Boolean;
  110. begin
  111. FillChar(MustSave,Sizeof(MustSave),False);
  112. for i:=1 to n do
  113. if Split[i]<>1 then begin
  114. MustSave[Split[i]]:=True;
  115. MustSave[i div Split[i]]:=True
  116. end;
  117. FillChar(Result,Sizeof(Result),0);
  118. Result.l:=1;Result.x[1]:=1;{1^m}
  119. for i:=2 to n do begin
  120. if Split[i]=1 then begin
  121. y.l:=1;y.x[1]:=i;
  122. Power(x,y,m)
  123. end else begin
  124. x:=Save[Split[i]]^;
  125. Mul(x,Save[i div Split[i]]^)
  126. end;
  127. if MustSave[i] then begin
  128. New(Save[i]);
  129. Save[i]^:=x
  130. end;
  131. Add(Result,x)
  132. end
  133. end;
  134. begin
  135. Write('Input n = ');Readln(n);
  136. Write('Input m = ');Readln(m);
  137. Write('Input k = ');Readln(k);
  138. kk:=(k-1) div 4+1;
  139. Time:=MemL[$40:$6C];
  140. InitSplit;
  141. Main;
  142. Time:=MemL[$40:$6C]-Time;
  143. Write('1^',m,'+2^',m,'+...+',n,'^',m,' = ...');
  144. Print;
  145. Writeln('Time Used = ',Time);
  146. end.