new; cls; /* Calculates a^n%m where powermod(a,n,m) uses mmod(a,n,m) for more digits. */
print powermod(1234,56789012,3456789);
print powermod2(1234,56789012,3456789);


proc powermod(a,n,m);
   local v;
   v=1;
   do while n>=1;
      if mod(n,2)==1;
         v=mod(a*v,m);
      endif;
      a=mod(a^2,m);
      n=floor(n/2);
   endo;
   retp(v);
endp;

proc powermod2(a,n,m);
   local v;
   a=mod(a,m);
   if a<0;
      a=a+m;
   endif;
   v=1;
   do while n>=1;
      if mod(n,2)==1;
         v=mmod(a,v,m);
      endif;
      a=mmod(a,a,m);
      n=floor(n/2);
   endo;
   retp(v);
endp;

proc mmod(a,n,m);
   local mul;
   a=mod(a,m);
   if a<0;
      a=a+m;
   endif;
   mul=0;
   do while n>=1;
      if mod(n,2)==1;
         mul=mod(a+mul,m);
      endif;
      a=a*mod(2,m);
      n=floor(n/2);
   endo;
   retp(mul);
endp;

proc mod(x,m);
   retp( (x.>=0).*(x%m)+(x.<0).*((x+m.*ceil(abs(x)/m))%m) );
endp;