#include "WolframLibrary.h"
#include <stdlib.h>
#include <string.h>
#include <math.h>

/* fast classical-cipher search primitives for WolframInstitute/Puzzles.
   Index space is caller-defined (e.g. KRYPTOS-keyed): all letters are 0..25 ints. */

static double *QG = NULL;      /* 26^4 quadgram log-probs, flat */
static double  QFLOOR = -20.0;

/* xorshift128+ RNG (deterministic from seed) */
static unsigned long long RS0, RS1;
static void rseed(unsigned long long s){ RS0 = s?s:0x9E3779B97F4A7C15ULL; RS1 = RS0^0xD1B54A32D192ED03ULL; }
static inline unsigned long long rnext(void){
  unsigned long long x=RS0, y=RS1; RS0=y;
  x^=x<<23; x^=x>>17; x^=y^(y>>26); RS1=x; return x+y;
}
static inline int rr26(int n){ return (int)(rnext()%(unsigned)n); }
static inline double runif(void){ return (rnext()>>11)*(1.0/9007199254740992.0); }

DLLEXPORT int setQuadgrams(WolframLibraryData ld, mint argc, MArgument *args, MArgument res){
  MTensor t = MArgument_getMTensor(args[0]);
  mint n = ld->MTensor_getFlattenedLength(t);
  double *d = ld->MTensor_getRealData(t);
  if(QG) free(QG);
  QG = (double*)malloc(sizeof(double)*456976);
  if(!QG) return LIBRARY_FUNCTION_ERROR;
  for(mint i=0;i<456976 && i<n;i++) QG[i]=d[i];
  QFLOOR = MArgument_getReal(args[1]);
  MArgument_setInteger(res, n);
  return LIBRARY_NO_ERROR;
}

static inline double qscore(const int *p, int n){
  double s=0.0; int i;
  for(i=0;i<n-3;i++){
    int idx=((p[i]*26+p[i+1])*26+p[i+2])*26+p[i+3];
    s += QG[idx];
  }
  return s;
}

/* Vigenere/Quagmire hill-climb: plain[i] = (ct[i]-shift[i%p]) mod 26.
   args: ctIdx (Integer tensor), period, iters, restarts, seed, initShifts (Integer tensor, len p or 0).
   returns Real tensor {bestScore, shift0..shift_{p-1}} */
DLLEXPORT int hillClimbVig(WolframLibraryData ld, mint argc, MArgument *args, MArgument res){
  MTensor ctT = MArgument_getMTensor(args[0]);
  mint n = ld->MTensor_getFlattenedLength(ctT);
  mint *ctm = ld->MTensor_getIntegerData(ctT);
  int p = (int)MArgument_getInteger(args[1]);
  long iters = (long)MArgument_getInteger(args[2]);
  int restarts = (int)MArgument_getInteger(args[3]);
  rseed((unsigned long long)MArgument_getInteger(args[4]));
  MTensor initT = MArgument_getMTensor(args[5]);
  mint initLen = ld->MTensor_getFlattenedLength(initT);
  mint *initm = initLen>0 ? ld->MTensor_getIntegerData(initT) : NULL;

  int *ct=(int*)malloc(sizeof(int)*n);
  int *pl=(int*)malloc(sizeof(int)*n);
  int i; for(i=0;i<n;i++) ct[i]=(int)ctm[i];
  int *sh=(int*)malloc(sizeof(int)*p);
  int *bestsh=(int*)malloc(sizeof(int)*p);
  double best=-1e300;

  int r; for(r=0;r<restarts;r++){
    if(r==0 && initm){ for(i=0;i<p;i++) sh[i]=(int)initm[i]%26; }
    else { for(i=0;i<p;i++) sh[i]=rr26(26); }
    for(i=0;i<n;i++) pl[i]=(ct[i]-sh[i%p]+26)%26;
    double cur=qscore(pl,n);
    long it; for(it=0;it<iters;it++){
      int j=rr26(p), old=sh[j], nv=rr26(26);
      if(nv==old) continue;
      sh[j]=nv;
      /* recompute only affected positions? simpler: full recompute of column j positions */
      int k; for(k=j;k<n;k+=p) pl[k]=(ct[k]-nv+26)%26;
      double v=qscore(pl,n);
      double T=0.02+ (1.2)*(1.0-(double)it/iters); if(T<0.02)T=0.02;
      if(v>cur || runif() < exp((v-cur)/(T*n*0.02))){ cur=v; }
      else { sh[j]=old; for(k=j;k<n;k+=p) pl[k]=(ct[k]-old+26)%26; }
      if(cur>best){ best=cur; memcpy(bestsh,sh,sizeof(int)*p); }
    }
    /* greedy polish of this restart's best */
  }
  MTensor out; mint dims[1]; dims[0]=p+1;
  ld->MTensor_new(MType_Real,1,dims,&out);
  double *od=ld->MTensor_getRealData(out);
  od[0]=best; for(i=0;i<p;i++) od[i+1]=(double)bestsh[i];
  MArgument_setMTensor(res,out);
  free(ct);free(pl);free(sh);free(bestsh);
  return LIBRARY_NO_ERROR;
}

/* 5x5-square hill-climb for Playfair (mode 0) and Bifid (mode 1).  ct is English-letter
   indices 0..25 with J(9) already folded to I(8); the square is a permutation of the 25
   indices {0..25}\{9}.  Decrypt on the candidate square, score by quadgrams (QG, English-
   indexed).  args: ct, mode, period(bifid), iters, restarts, seed.
   returns {score, sq0..sq24} (English index at each of the 25 grid cells). */
static void squareDecrypt(const int *ct, int n, const int *sq, const int *pos, int mode, int period, int *pl){
  if(mode==0){ /* Playfair: digraphs */
    for(int i=0;i+1<n;i+=2){
      int pa=pos[ct[i]], pb=pos[ct[i+1]];
      int ra=pa/5,ca=pa%5, rb=pb/5,cb=pb%5, na,nb;
      if(ra==rb){ na=ra*5+(ca+4)%5; nb=rb*5+(cb+4)%5; }
      else if(ca==cb){ na=((ra+4)%5)*5+ca; nb=((rb+4)%5)*5+cb; }
      else { na=ra*5+cb; nb=rb*5+ca; }
      pl[i]=sq[na]; pl[i+1]=sq[nb];
    }
    if(n&1) pl[n-1]=ct[n-1];
  } else { /* Bifid: per-period block */
    for(int start=0;start<n;start+=period){
      int L=n-start<period?n-start:period;
      int rc[512]; /* coords in pair order: r0,c0,r1,c1,... */
      for(int k=0;k<L;k++){ int p=pos[ct[start+k]]; rc[2*k]=p/5; rc[2*k+1]=p%5; }
      /* first L of the flattened stream are plaintext rows, next L are cols */
      for(int k=0;k<L;k++){ int r=rc[k], c=rc[L+k]; pl[start+k]=sq[r*5+c]; }
    }
  }
}
DLLEXPORT int hillClimbSquare(WolframLibraryData ld, mint argc, MArgument *args, MArgument res){
  MTensor ctT=MArgument_getMTensor(args[0]); mint n=ld->MTensor_getFlattenedLength(ctT);
  mint *ctm=ld->MTensor_getIntegerData(ctT);
  int mode=(int)MArgument_getInteger(args[1]); int period=(int)MArgument_getInteger(args[2]);
  long iters=(long)MArgument_getInteger(args[3]); int restarts=(int)MArgument_getInteger(args[4]);
  rseed((unsigned long long)MArgument_getInteger(args[5]));
  int *ct=(int*)malloc(sizeof(int)*n), *pl=(int*)malloc(sizeof(int)*n);
  int i; for(i=0;i<n;i++){ int x=(int)ctm[i]; ct[i]= x==9?8:x; }
  int base[25]; { int j=0; for(i=0;i<26;i++) if(i!=9) base[j++]=i; }
  int sq[25], pos[26], bsq[25]; double best=-1e300;
  int r; for(r=0;r<restarts;r++){
    for(i=0;i<25;i++) sq[i]=base[i];
    for(i=24;i>0;i--){ int k=rr26(i+1); int t=sq[i]; sq[i]=sq[k]; sq[k]=t; }
    for(i=0;i<25;i++) pos[sq[i]]=i; pos[9]=pos[8];
    squareDecrypt(ct,n,sq,pos,mode,period,pl);
    double cur=qscore(pl,n);
    long it; for(it=0;it<iters;it++){
      int a=rr26(25), b=rr26(25); if(a==b) continue;
      int t=sq[a]; sq[a]=sq[b]; sq[b]=t; pos[sq[a]]=a; pos[sq[b]]=b; pos[9]=pos[8];
      squareDecrypt(ct,n,sq,pos,mode,period,pl);
      double v=qscore(pl,n); double T=1.0+30.0*(1.0-(double)it/iters);
      if(v>cur || runif()<exp((v-cur)/T)){ cur=v; }
      else { t=sq[a]; sq[a]=sq[b]; sq[b]=t; pos[sq[a]]=a; pos[sq[b]]=b; pos[9]=pos[8]; }
      if(cur>best){ best=cur; memcpy(bsq,sq,sizeof(int)*25); }
    }
  }
  MTensor out; mint dims[1]; dims[0]=26; ld->MTensor_new(MType_Real,1,dims,&out);
  double *od=ld->MTensor_getRealData(out); od[0]=best; for(i=0;i<25;i++) od[1+i]=bsq[i];
  MArgument_setMTensor(res,out);
  free(ct);free(pl); return LIBRARY_NO_ERROR;
}

/* Autokey (Vigenere) hill-climb: keystream is primer(m) then the recovered plaintext.
   plain[i] = (ct[i]-k[i]) mod 26 ; k[i]=primer[i] (i<m) else plain[i-m].
   args: ct, m, iters, restarts, seed.  returns {score, primer0..primer_{m-1}} */
DLLEXPORT int hillClimbAutokey(WolframLibraryData ld, mint argc, MArgument *args, MArgument res){
  MTensor ctT=MArgument_getMTensor(args[0]); mint n=ld->MTensor_getFlattenedLength(ctT);
  mint *ctm=ld->MTensor_getIntegerData(ctT);
  int m=(int)MArgument_getInteger(args[1]); long iters=(long)MArgument_getInteger(args[2]);
  int restarts=(int)MArgument_getInteger(args[3]);
  rseed((unsigned long long)MArgument_getInteger(args[4]));
  int *ct=(int*)malloc(sizeof(int)*n), *pl=(int*)malloc(sizeof(int)*n);
  int i; for(i=0;i<n;i++) ct[i]=(int)ctm[i];
  int *pr=(int*)malloc(sizeof(int)*m), *bp=(int*)malloc(sizeof(int)*m);
  double best=-1e300;
  int r; for(r=0;r<restarts;r++){
    for(i=0;i<m;i++) pr[i]=rr26(26);
    /* decode */
    for(i=0;i<n;i++){ int k=(i<m)?pr[i]:pl[i-m]; pl[i]=(ct[i]-k+26)%26; }
    double cur=qscore(pl,n);
    long it; for(it=0;it<iters;it++){
      int j=rr26(m), old=pr[j], nv=rr26(26); if(nv==old) continue; pr[j]=nv;
      for(i=0;i<n;i++){ int k=(i<m)?pr[i]:pl[i-m]; pl[i]=(ct[i]-k+26)%26; }
      double v=qscore(pl,n); double T=0.02+0.8*(1.0-(double)it/iters); if(T<0.02)T=0.02;
      if(v>cur || runif()<exp((v-cur)/(T*n*0.02))){ cur=v; }
      else { pr[j]=old; for(i=0;i<n;i++){ int k=(i<m)?pr[i]:pl[i-m]; pl[i]=(ct[i]-k+26)%26; } }
      if(cur>best){ best=cur; memcpy(bp,pr,sizeof(int)*m); }
    }
  }
  MTensor out; mint dims[1]; dims[0]=m+1; ld->MTensor_new(MType_Real,1,dims,&out);
  double *od=ld->MTensor_getRealData(out); od[0]=best; for(i=0;i<m;i++) od[1+i]=bp[i];
  MArgument_setMTensor(res,out);
  free(ct);free(pl);free(pr);free(bp); return LIBRARY_NO_ERROR;
}

/* Quagmire I/II/III hill-climb (ciphertext-only, keyword-seedable).
   With mixed alphabet ka (permutation, inv=ka^-1) and per-column offset o=off[j]:
     variant 1 (Quag I) : plain mixed, cipher straight -> p = ka[(c    - o) mod 26]
     variant 2 (Quag II): plain straight, cipher mixed -> p = (inv[c] - o) mod 26
     variant 3 (Quag III, Kryptos): both mixed       -> p = ka[(inv[c]- o) mod 26]
   Offsets are chosen per candidate alphabet by minimum unigram chi-square vs English,
   the alphabet+offsets are annealed on quadgram fitness.  Seeds may be supplied as keyed
   alphabets (each restart r<nSeeds starts from seed row r), else a random permutation.
   args: ct, variant, period, iters, restarts, seed, seedAlphabets (Integer, nSeeds*26 or 0).
   returns Real {score, ka0..ka25, off0..off_{period-1}} (length 27+period). */
static double EFREQ[26];
static void quagOffsets(const int *ct,int n,int period,int variant,const int *inv,const int *ka,int *off){
  int j; for(j=0;j<period;j++){
    int hist[26]; int i; for(i=0;i<26;i++) hist[i]=0;
    for(i=j;i<n;i+=period){ int y = (variant==1)? ct[i] : inv[ct[i]]; hist[y]++; }
    double bestsc=-1e300; int bo=0, o;
    for(o=0;o<26;o++){
      double sc=0.0; int pos;
      for(pos=0;pos<26;pos++){ int cnt=hist[(pos+o)%26]; if(!cnt) continue;
        int letter = (variant==2)? pos : ka[pos];
        double f=EFREQ[letter]; sc += cnt*log(f>1e-9?f:1e-9); }
      if(sc>bestsc){ bestsc=sc; bo=o; }
    }
    off[j]=bo;
  }
}
static inline void quagDecrypt(const int *ct,int n,int period,int variant,const int *inv,const int *ka,const int *off,int *pl){
  int i; for(i=0;i<n;i++){ int o=off[i%period]; int c=ct[i];
    pl[i] = (variant==2) ? ((inv[c]-o)%26+26)%26
          : (variant==1) ? ka[((c-o)%26+26)%26]
                         : ka[((inv[c]-o)%26+26)%26]; }
}
DLLEXPORT int hillClimbQuagmire(WolframLibraryData ld, mint argc, MArgument *args, MArgument res){
  MTensor ctT=MArgument_getMTensor(args[0]); mint n=ld->MTensor_getFlattenedLength(ctT);
  mint *ctm=ld->MTensor_getIntegerData(ctT);
  int variant=(int)MArgument_getInteger(args[1]);
  int period=(int)MArgument_getInteger(args[2]);
  long iters=(long)MArgument_getInteger(args[3]);
  int restarts=(int)MArgument_getInteger(args[4]);
  rseed((unsigned long long)MArgument_getInteger(args[5]));
  MTensor seedT=MArgument_getMTensor(args[6]); mint seedLen=ld->MTensor_getFlattenedLength(seedT);
  mint *seedm = seedLen>0 ? ld->MTensor_getIntegerData(seedT) : NULL;
  int nSeeds = (int)(seedLen/26);
  int *ct=(int*)malloc(sizeof(int)*n), *pl=(int*)malloc(sizeof(int)*n);
  int i; for(i=0;i<n;i++) ct[i]=(int)ctm[i];
  int ka[26],inv[26],bka[26]; int *off=(int*)malloc(sizeof(int)*period);
  int *boff=(int*)malloc(sizeof(int)*period);
  double best=-1e300;
  int r; for(r=0;r<restarts;r++){
    if(r<nSeeds){ for(i=0;i<26;i++) ka[i]=(int)seedm[r*26+i]; }
    else { for(i=0;i<26;i++) ka[i]=i; for(i=25;i>0;i--){ int k=rr26(i+1); int t=ka[i]; ka[i]=ka[k]; ka[k]=t; } }
    for(i=0;i<26;i++) inv[ka[i]]=i;
    quagOffsets(ct,n,period,variant,inv,ka,off);
    quagDecrypt(ct,n,period,variant,inv,ka,off,pl);
    double cur=qscore(pl,n), bloc=cur; int bkaL[26]; int *boffL=(int*)malloc(sizeof(int)*period);
    memcpy(bkaL,ka,sizeof ka); memcpy(boffL,off,sizeof(int)*period);
    long it; for(it=0;it<iters;it++){
      double T=0.03+4.0*(1.0-(double)it/iters);
      int mv=rr26(100), a=0,b=0,jc=0,dl=0;
      if(mv<78){ a=rr26(26); b=rr26(26); if(a==b){ continue; }
        int t=ka[a]; ka[a]=ka[b]; ka[b]=t; inv[ka[a]]=a; inv[ka[b]]=b; }
      else { jc=rr26(period); dl=1+rr26(25); off[jc]=(off[jc]+dl)%26; }
      quagDecrypt(ct,n,period,variant,inv,ka,off,pl);
      double v=qscore(pl,n);
      if(v>cur || runif()<exp((v-cur)/T)){ cur=v;
        if(v>bloc){ bloc=v; memcpy(bkaL,ka,sizeof ka); memcpy(boffL,off,sizeof(int)*period); } }
      else { if(mv<78){ int t=ka[a]; ka[a]=ka[b]; ka[b]=t; inv[ka[a]]=a; inv[ka[b]]=b; }
             else off[jc]=((off[jc]-dl)%26+26)%26; }
    }
    if(bloc>best){ best=bloc; memcpy(bka,bkaL,sizeof bka); memcpy(boff,boffL,sizeof(int)*period); }
    free(boffL);
  }
  MTensor out; mint dims[1]; dims[0]=27+period; ld->MTensor_new(MType_Real,1,dims,&out);
  double *od=ld->MTensor_getRealData(out); od[0]=best;
  for(i=0;i<26;i++) od[1+i]=bka[i]; for(i=0;i<period;i++) od[27+i]=boff[i];
  MArgument_setMTensor(res,out);
  free(ct);free(pl);free(off);free(boff);
  return LIBRARY_NO_ERROR;
}

/* --- unigram peel of a double Quagmire (transposition-invariant) --- */
DLLEXPORT int setEnglishFreq(WolframLibraryData ld, mint argc, MArgument *args, MArgument res){
  MTensor t=MArgument_getMTensor(args[0]); double *d=ld->MTensor_getRealData(t);
  int i; for(i=0;i<26;i++) EFREQ[i]=d[i];
  MArgument_setInteger(res,26); return LIBRARY_NO_ERROR;
}
static inline double unichi(const int *p, int n){
  int cnt[26]; int i; for(i=0;i<26;i++) cnt[i]=0;
  for(i=0;i<n;i++) cnt[p[i]]++;
  double s=0.0; for(i=0;i<26;i++){ double e=n*EFREQ[i]; double d=cnt[i]-e; s+=d*d/e; }
  return s; /* lower = more English */
}
/* hill-climb SA(a)+SB(b) minimizing unigram chi2 of de-Quagmired text.
   eff[i]=(SA[i%a]+SB[i%b])%26 ; plain[i]=(ct[i]-eff[i])%26.
   returns {chi2, SA0..SA{a-1}, SB0..SB{b-1}} */
DLLEXPORT int peelDouble(WolframLibraryData ld, mint argc, MArgument *args, MArgument res){
  MTensor ctT=MArgument_getMTensor(args[0]); mint n=ld->MTensor_getFlattenedLength(ctT);
  mint *ctm=ld->MTensor_getIntegerData(ctT);
  int a=(int)MArgument_getInteger(args[1]), b=(int)MArgument_getInteger(args[2]);
  long iters=(long)MArgument_getInteger(args[3]); int restarts=(int)MArgument_getInteger(args[4]);
  rseed((unsigned long long)MArgument_getInteger(args[5]));
  int *ct=(int*)malloc(sizeof(int)*n), *pl=(int*)malloc(sizeof(int)*n);
  int i; for(i=0;i<n;i++) ct[i]=(int)ctm[i];
  int *SA=(int*)malloc(sizeof(int)*a), *SB=(int*)malloc(sizeof(int)*b);
  int *bA=(int*)malloc(sizeof(int)*a), *bB=(int*)malloc(sizeof(int)*b);
  double best=1e300;
  int r; for(r=0;r<restarts;r++){
    for(i=0;i<a;i++) SA[i]=rr26(26); for(i=0;i<b;i++) SB[i]=rr26(26);
    for(i=0;i<n;i++) pl[i]=(ct[i]-((SA[i%a]+SB[i%b])%26)+26)%26;
    double cur=unichi(pl,n);
    long it; for(it=0;it<iters;it++){
      int inA=(rr26(a+b)<a); int k,old,nv;
      if(inA){ k=rr26(a); old=SA[k]; nv=rr26(26); if(nv==old) continue; SA[k]=nv; }
      else   { k=rr26(b); old=SB[k]; nv=rr26(26); if(nv==old) continue; SB[k]=nv; }
      for(i=0;i<n;i++) pl[i]=(ct[i]-((SA[i%a]+SB[i%b])%26)+26)%26;
      double v=unichi(pl,n);
      double T=0.5+ 40.0*(1.0-(double)it/iters); if(T<0.5)T=0.5;
      if(v<cur || runif()<exp(-(v-cur)/T)){ cur=v; }
      else { if(inA) SA[k]=old; else SB[k]=old; }
      if(cur<best){ best=cur; memcpy(bA,SA,sizeof(int)*a); memcpy(bB,SB,sizeof(int)*b); }
    }
  }
  MTensor out; mint dims[1]; dims[0]=a+b+1; ld->MTensor_new(MType_Real,1,dims,&out);
  double *od=ld->MTensor_getRealData(out);
  od[0]=best; for(i=0;i<a;i++) od[1+i]=bA[i]; for(i=0;i<b;i++) od[1+a+i]=bB[i];
  MArgument_setMTensor(res,out);
  free(ct);free(pl);free(SA);free(SB);free(bA);free(bB);
  return LIBRARY_NO_ERROR;
}

/* polish a composite (double-Quagmire + known transposition) by quadgram of the readable,
   de-transposed plaintext.  perm[i] = source position in the de-Quagmired text that plaintext
   position i draws from.  plain[i] = (ct[perm[i]] - (SA[perm[i]%a]+SB[perm[i]%b])) mod 26.
   Hill-climb SA,SB from the given init.  args: ct, perm, a, b, initSA, initSB, iters, seed.
   returns {score, SA..., SB...} */
DLLEXPORT int polishComposite(WolframLibraryData ld, mint argc, MArgument *args, MArgument res){
  MTensor ctT=MArgument_getMTensor(args[0]); mint n=ld->MTensor_getFlattenedLength(ctT);
  mint *ctm=ld->MTensor_getIntegerData(ctT);
  MTensor pmT=MArgument_getMTensor(args[1]); mint *pm=ld->MTensor_getIntegerData(pmT);
  int a=(int)MArgument_getInteger(args[2]), b=(int)MArgument_getInteger(args[3]);
  mint *iA=ld->MTensor_getIntegerData(MArgument_getMTensor(args[4]));
  mint *iB=ld->MTensor_getIntegerData(MArgument_getMTensor(args[5]));
  long iters=(long)MArgument_getInteger(args[6]);
  rseed((unsigned long long)MArgument_getInteger(args[7]));
  int i;
  int *cp=(int*)malloc(sizeof(int)*n), *ai=(int*)malloc(sizeof(int)*n), *bi=(int*)malloc(sizeof(int)*n), *pl=(int*)malloc(sizeof(int)*n);
  for(i=0;i<n;i++){ int sp=(int)pm[i]; cp[i]=(int)ctm[sp]; ai[i]=sp%a; bi[i]=sp%b; }
  int *SA=(int*)malloc(sizeof(int)*a), *SB=(int*)malloc(sizeof(int)*b), *bA=(int*)malloc(sizeof(int)*a), *bB=(int*)malloc(sizeof(int)*b);
  for(i=0;i<a;i++) SA[i]=(int)iA[i]%26; for(i=0;i<b;i++) SB[i]=(int)iB[i]%26;
  for(i=0;i<n;i++) pl[i]=(cp[i]-((SA[ai[i]]+SB[bi[i]])%26)+26)%26;
  double cur=qscore(pl,n), best=cur; memcpy(bA,SA,sizeof(int)*a); memcpy(bB,SB,sizeof(int)*b);
  long it; for(it=0;it<iters;it++){
    int inA=(rr26(a+b)<a); int k,old,nv;
    if(inA){ k=rr26(a); old=SA[k]; nv=rr26(26); if(nv==old)continue; SA[k]=nv; }
    else   { k=rr26(b); old=SB[k]; nv=rr26(26); if(nv==old)continue; SB[k]=nv; }
    for(i=0;i<n;i++) pl[i]=(cp[i]-((SA[ai[i]]+SB[bi[i]])%26)+26)%26;
    double v=qscore(pl,n);
    double T=0.02+0.8*(1.0-(double)it/iters); if(T<0.02)T=0.02;
    if(v>cur || runif()<exp((v-cur)/(T*n*0.02))){ cur=v; }
    else { if(inA)SA[k]=old; else SB[k]=old; }
    if(cur>best){ best=cur; memcpy(bA,SA,sizeof(int)*a); memcpy(bB,SB,sizeof(int)*b); }
  }
  MTensor out; mint dims[1]; dims[0]=a+b+1; ld->MTensor_new(MType_Real,1,dims,&out);
  double *od=ld->MTensor_getRealData(out); od[0]=best;
  for(i=0;i<a;i++) od[1+i]=bA[i]; for(i=0;i<b;i++) od[1+a+i]=bB[i];
  MArgument_setMTensor(res,out);
  free(cp);free(ai);free(bi);free(pl);free(SA);free(SB);free(bA);free(bB);
  return LIBRARY_NO_ERROR;
}

/* Running-key corpus attack: try every offset of a corpus as the running key of a Vigenere
   (op 0: plain=ct-key, op 1: plain=key-ct Beaufort, op 2: plain=ct+key).  ct and corpus are
   index arrays in the same (caller-chosen) space.  returns {bestScore, bestOffset}. */
DLLEXPORT int runKeyCorpus(WolframLibraryData ld, mint argc, MArgument *args, MArgument res){
  MTensor ctT=MArgument_getMTensor(args[0]); mint n=ld->MTensor_getFlattenedLength(ctT);
  mint *ctm=ld->MTensor_getIntegerData(ctT);
  MTensor coT=MArgument_getMTensor(args[1]); mint clen=ld->MTensor_getFlattenedLength(coT);
  mint *com=ld->MTensor_getIntegerData(coT);
  int op=(int)MArgument_getInteger(args[2]);
  int *ct=(int*)malloc(sizeof(int)*n); int i; for(i=0;i<n;i++) ct[i]=(int)ctm[i];
  int *cor=(int*)malloc(sizeof(int)*clen); for(i=0;i<clen;i++) cor[i]=(int)com[i];
  int *pl=(int*)malloc(sizeof(int)*n);
  double best=-1e300; long bestOff=-1, off;
  for(off=0; off+n<=clen; off++){
    for(i=0;i<n;i++){ int k=cor[off+i];
      pl[i]= op==0 ? (ct[i]-k+26)%26 : (op==1 ? (k-ct[i]+26)%26 : (ct[i]+k)%26); }
    double v=qscore(pl,n);
    if(v>best){ best=v; bestOff=off; }
  }
  MTensor out; mint dims[1]; dims[0]=2; ld->MTensor_new(MType_Real,1,dims,&out);
  double *od=ld->MTensor_getRealData(out); od[0]=best; od[1]=(double)bestOff;
  MArgument_setMTensor(res,out);
  free(ct);free(cor);free(pl); return LIBRARY_NO_ERROR;
}

/* brute all ncol! column orders of a columnar transposition, score de-transposed by quadgram.
   textIdx: Integer tensor (the intermediate, letters as indices). ncol<=10.
   returns Real tensor {bestScore, order0..order_{ncol-1}} (0-based read order). */
static int NC; static int *BestOrd; static double BestSc; static int *TXT; static long NN;
static int *plbuf;
static void perm_rec(int *ord, int depth, int used){
  if(depth==NC){
    /* de-transpose: block k (rows each) is source column ord[k]; read row-major */
    int rows=NN/NC; int k,r;
    for(k=0;k<NC;k++){ int c=ord[k]; for(r=0;r<rows;r++) plbuf[r*NC+c]=TXT[k*rows+r]; }
    double v=qscore(plbuf,NN);
    if(v>BestSc){ BestSc=v; memcpy(BestOrd,ord,sizeof(int)*NC); }
    return;
  }
  int c; for(c=0;c<NC;c++) if(!(used&(1<<c))){ ord[depth]=c; perm_rec(ord,depth+1,used|(1<<c)); }
}
DLLEXPORT int columnarBrute(WolframLibraryData ld, mint argc, MArgument *args, MArgument res){
  MTensor tT=MArgument_getMTensor(args[0]);
  NN=ld->MTensor_getFlattenedLength(tT);
  mint *tm=ld->MTensor_getIntegerData(tT);
  NC=(int)MArgument_getInteger(args[1]);
  TXT=(int*)malloc(sizeof(int)*NN); plbuf=(int*)malloc(sizeof(int)*NN);
  BestOrd=(int*)malloc(sizeof(int)*NC); int *ord=(int*)malloc(sizeof(int)*NC);
  long i; for(i=0;i<NN;i++) TXT[i]=(int)tm[i];
  BestSc=-1e300;
  perm_rec(ord,0,0);
  MTensor out; mint dims[1]; dims[0]=NC+1;
  ld->MTensor_new(MType_Real,1,dims,&out);
  double *od=ld->MTensor_getRealData(out);
  od[0]=BestSc; for(i=0;i<NC;i++) od[i+1]=(double)BestOrd[i];
  MArgument_setMTensor(res,out);
  free(TXT);free(plbuf);free(BestOrd);free(ord);
  return LIBRARY_NO_ERROR;
}
