POJ1845

题意:求A^B的约数和mod 9901.(0 <= A,B <= 50000000)

分析:记f(n)为n的约数和.
求f(a^b) mod c.

f(n)=∏(pi^(qi+1)-1)/(pi-1)
pi为质因子,qi为质因子个数.

(pi^(qi+1)-1)/(pi-1)=1+pi+pi^2+......+pi^qi

转化为等比数列的和.
可以用二分.例如:

1+p+p^2+p^3+p^4=p^2+(1+p^3)(1+p)
1+p+p^2+p^3+p^4+p^5=(1+p^3)(1+p+p^2)

递归进行.

code:

const mo=9901;
var   p,q:array[0..20] of longint;
      ans,tmp,num,a,b,i,j,n:longint;


      function pow(x,y,z:longint):longint;
      var   g,h:int64;
      begin
            g:=1;
            h:=x;
            while y>0 do
            begin
                  if y and 1=1 then g:=g*h mod z;
                  h:=h*h mod z;
                  y:=y>>1;
            end;
            exit(g);
      end;

      function work(a,b:longint):int64;
      begin
            if b=0 then exit(1);
            if b and 1=1 then
              exit((pow(a,b>>1+1,mo)+1)*work(a,b>>1) mod mo)
            else
              exit((pow(a,b>>1,mo)+(pow(a,b>>1+1,mo)+1)*work(a,b>>1-1)) mod mo);
      end;

begin
      readln(a,b);
      num:=a;
      for i:=2 to trunc(sqrt(a)) do
      if num mod i=0 then
      begin
            inc(n);
            p[n]:=i;
            while num mod i=0 do
            begin
                  inc(q[n]);
                  num:=num div i;
            end;
      end;
      if num>1 then
      begin
            inc(n);
            p[n]:=num;
            q[n]:=1;
      end;

      for i:=1 to n do q[i]:=q[i]*b;
      ans:=1;
      for i:=1 to n do
      begin
            tmp:=work(p[i],q[i]);
            ans:=ans*tmp mod mo;
      end;
      writeln(ans);
end.