In[]:=
(*=====================================================================*)(*FASTZETATRANSFORMOVERGF(2)ANDSETSEQUENCEGENERATOR*)(*=====================================================================*)ClearAll[FastZetaTransform,BitwiseORConvolution,ComputeSequence];(*1.FastZetaTransform(FZT)Modulo2*)(*Performsanin-placetransformationinO(NlogN)time*)FastZetaTransform[arr_List]:=Module[{a=arr,n=Length[arr],step=1,i,j},While[step<n,i=1;While[i<=n,For[j=0,j<step,j++,(*a[i+j+step]=a[i+j+step]\[OPlus]a[i+j]*)a[[i+j+step]]=BitXor[a[[i+j+step]],a[[i+j]]]];i+=2*step;];step*=2;];a];(*2.BitwiseORConvolutionusingFZT*)(*Convolutiontheorem:InverseFZT(FZT(A)·FZT(B))*)BitwiseORConvolution[AList_,BList_]:=FastZetaTransform[BitAnd[FastZetaTransform[AList],FastZetaTransform[BList]]];(*3.CoreSequenceGenerator*)(*ComputesthebooleanarrayforS_muptoamaximumintegervalueR*)ComputeSequence[m_Integer,R_Integer]:=Module[{k,NSize,A,B,C,nextArr,i},(*Padarraystothenextpowerof2forFZTcompatibility*)k=Ceiling[Log2[R+1]];NSize=2^k;(*BaseCases:S_1={0},S_2={1}.NoteMathematicais1-indexed.*)A=ConstantArray[0,NSize];A[[1]]=1;B=ConstantArray[0,NSize];B[[2]]=1;If[m==1,Return[A]];If[m==2,Return[B]];(*Evaluaterecurrenceuptostepm*)For[i=3,i<=m,i++,(*ComputeA\[OTimes]B*)C=BitwiseORConvolution[A,B];(*ComputeA\[OPlus]B\[OPlus]C*)nextArr=BitXor[A,BitXor[B,C]];(*ApplyInc()byright-shiftingthebooleanarray*)nextArr=Prepend[Drop[nextArr,-1],0];(*Updatehistory*)A=B;B=nextArr;];B];(*=====================================================================*)(*HELPERFUNCTIONSFORUSAGE*)(*=====================================================================*)(*PredicateFunction:ReturnsTrueifr∈S_m*)IsInSetQ[m_Integer,r_Integer]:=Module[{arr},If[r<0,Return[False]];arr=ComputeSequence[m,r];arr[[r+1]]==1];(*SetExtractionFunction:ReturnstheactuallistofelementsinS_muptoMaxR*)GetSetElements[m_Integer,MaxR_Integer]:=Flatten[Position[ComputeSequence[m,MaxR],1]]-1;(*ExampleExecution*)Print["S_8 up to 20: ",GetSetElements[8,20]];Print["Is 16 in S_8? ",IsInSetQ[8,16]];Print["Is 13 in S_8? ",IsInSetQ[8,13]];Print["Is 500 in S_1000? ",IsInSetQ[1000,500]];Print["Is 501 in S_1000? ",IsInSetQ[1000,501]];
S_8 up to 20: {4,6,8,9,12,14,16}
Is 16 in S_8? True
Is 13 in S_8? False
Is 500 in S_1000? True
Is 501 in S_1000? False