/****************************************************************************** * * * COMPUTE MOBIUS FUNCTION (Deleglise & Rivat's algorithm) * * 08/31/15 (dkc) * * * ******************************************************************************/ #include <math.h> int newmob(unsigned int a, unsigned int b, int *out, unsigned int *table) { unsigned int i,count,p,ps,pmax,beg,rem,flag; int t; count=b-a; for (i=0; i<count; i++) out[i]=1; pmax=(unsigned int)sqrt((double)b); flag=0; for (i=0; i<17984; i++) { p=table[i]; if (p>pmax) { flag=1; break; } ps=p*p; rem=a-(a/ps)*ps; if (rem!=0) beg=ps-rem; else beg=0; while (beg<count) { out[beg]=0; beg=beg+ps; } rem=a-(a/p)*p; if (rem!=0) beg=p-rem; else beg=0; while (beg<count) { out[beg]=-out[beg]*p; beg=beg+p; } } if (flag==0) return(-1); for (i=0; i<count; i++) { if (out[i]==0) continue; t=out[i]; if (t<0) t=-t; if ((unsigned int)t<(i+a)) out[i]=-out[i]; } for (i=0; i<count; i++) { if (out[i]==0) continue; if (out[i]>0) out[i]=1; else out[i]=-1; } return(1); }