Computing Mersenne 40

I ran across this program during one of my searches for obscure "bignum" code. It computes Mersenne 40 possibly using a Number Theoretic Transform.

I have no idea who wrote it and only a vague idea of how it works. Any information on either point would be appreciated.

The code was even more "obfuscated" when I found it. I haven't changed anything except to make it readable.

Depending on whether you're a purist or a pragmatist you will either be offended by its terseness or awed by its elegance.

Interestingly it compiles without warnings using Visual Studio under Windows or GCC under Linux. But it doesn't work properly under Windows. YMMV.

/* Computes Mersenne 40 2^20996011-1 */
int m=754974721,N,t[1<<24],a,*p,i,e=28046341,j,s,b,c,U;

f(d){

  for(s=1<<23;s;s/=2,d=d*1LL*d%m)
    if(s<N)
      for(p=t;p<t+N;p+=s)
        for(i=s,c=1;i;i--)
          b=*p+p[s],p[s]=(m+*p-p[s])*1LL*c%m,*p++=b%m,c=c*1LL*d%m;

  for(j=0;i<N-1;){
    for(s=N/2;!((j^=s)&s);s/=2);
    if(++i<j)a=t[i],t[i]=t[j],t[j]=a;
  }
           
}

main()
{
  *t=2;
  U=N=1;
  while(e/=2){
    N*=2;
    U=U*1LL*(m+1)/2%m;
    f(362);
    for(p=t;p<t+N;)
      *p++=*p*1LL**p%m*U%m;
    f(415027540);
    for(a=0,p=t;p<t+N;)
      a+=*p<<(e&1),*p++=a%10,a/=10;
   }
   while(!*--p);
   t[0]--;
   while(p>=t)
     printf("%d",*p--);
}