// C++ program from Daewon Kim, Aug 2026.
// Computes A048833(k) modulo a prime for 0 <= k <= N/2.
// Compile: g++ -O3 -march=native -mavx2 -fopenmp a048833.cpp.txt -x c++ -o a048833
// Usage:   ./a048833 N PRIME OUT [RSTART REND]
// The full exact table is reconstructed from runs modulo 17 primes by CRT.

#include <bits/stdc++.h>
#include <immintrin.h>
#include <omp.h>
using namespace std;
static uint64_t modpow(uint64_t a,uint64_t e,uint64_t m){uint64_t r=1;while(e){if(e&1)r=(uint64_t)((__uint128_t)r*a%m);a=(uint64_t)((__uint128_t)a*a%m);e>>=1;}return r;}
static inline uint32_t addm(uint32_t a,uint32_t b,uint32_t m){uint32_t x=a+b;return x>=m?x-m:x;}
static inline uint32_t subm(uint32_t a,uint32_t b,uint32_t m){return a>=b?a-b:a+m-b;}
int main(int argc,char**argv){
 if(argc<4){cerr<<"usage N PRIME OUT [RSTART REND]\n";return 2;}
 int N=atoi(argv[1]);uint32_t MOD=(uint32_t)strtoul(argv[2],nullptr,10);string outp=argv[3];
 int M=1,m=0;while((M<<1)<=N){M<<=1;m++;}if(N>=2*M)return 3;
 int reps=M/2;int rs=argc>=5?atoi(argv[4]):0;int re=argc>=6?atoi(argv[5]):reps;rs=max(rs,0);re=min(re,reps);
 int T=omp_get_max_threads();vector<vector<uint32_t>> partial(T,vector<uint32_t>(N/2+1));
 __m256i modv=_mm256_set1_epi32((int)MOD);
 double st=omp_get_wtime();
#pragma omp parallel
 {
  int tid=omp_get_thread_num();vector<uint32_t> f(N+8);
#pragma omp for schedule(static)
  for(int r=rs;r<re;r++){
   int u=2*r;fill(f.begin(),f.end(),0u);f[0]=1;
   for(int j=1;j<M;j++){
    bool neg=__builtin_parity((unsigned)(u&j));
    int n=j;
    if(j>=8){
      int end=N-7;
      if(!neg){
        for(;n<=end;n+=8){
          __m256i a=_mm256_loadu_si256((const __m256i*)&f[n]);
          __m256i b=_mm256_loadu_si256((const __m256i*)&f[n-j]);
          __m256i s=_mm256_add_epi32(a,b);
          __m256i sm=_mm256_sub_epi32(s,modv);
          __m256i v=_mm256_min_epu32(s,sm);
          _mm256_storeu_si256((__m256i*)&f[n],v);
        }
      }else{
        for(;n<=end;n+=8){
          __m256i a=_mm256_loadu_si256((const __m256i*)&f[n]);
          __m256i b=_mm256_loadu_si256((const __m256i*)&f[n-j]);
          __m256i d=_mm256_sub_epi32(a,b);
          __m256i dp=_mm256_add_epi32(d,modv);
          __m256i v=_mm256_min_epu32(d,dp);
          _mm256_storeu_si256((__m256i*)&f[n],v);
        }
      }
    }
    if(!neg){for(;n<=N;n++)f[n]=addm(f[n],f[n-j],MOD);}else{for(;n<=N;n++)f[n]=subm(f[n],f[n-j],MOD);}
   }
   auto &acc=partial[tid];for(int k=0;k<=N/2;k++)acc[k]=addm(acc[k],f[2*k],MOD);
  }
 }
 ofstream out(outp);if(!out)return 4;
 bool full=(rs==0&&re==reps);uint64_t inv=full?modpow(M/2,MOD-2,MOD):1;
 for(int k=0;k<=N/2;k++){uint64_t tot=0;for(int t=0;t<T;t++)tot+=partial[t][k];uint64_t v=tot%MOD;if(full)v=(uint64_t)((__uint128_t)v*inv%MOD);out<<k<<' '<<v<<'\n';}
 cerr<<"N="<<N<<" M="<<M<<" reps=["<<rs<<","<<re<<") elapsed="<<omp_get_wtime()-st<<"\n";
}
