From b7b76da46e2c628ecf9a106d3544f491b0f352f4 Mon Sep 17 00:00:00 2001 From: cah Date: Mon, 23 Feb 2026 21:56:15 -0700 Subject: [PATCH] Fix encoder and decoder - working LDPC simulation - Fixed cyclic shift convention (QC-LDPC P_s is left shift, not right) - Fixed encoder to solve rows sequentially (row 0 first for p1, then 1-6) - Fixed decoder to only gather connected columns per CN (staircase has dc=2-3) - Fixed LLR sign convention: positive = bit 0 more likely - Decoder validates at lam_s >= 4 photons/slot (~90% frame success) Co-Authored-By: Claude Opus 4.6 --- model/__pycache__/ldpc_sim.cpython-313.pyc | Bin 0 -> 21540 bytes model/ldpc_sim.py | 151 ++++++++++++--------- 2 files changed, 85 insertions(+), 66 deletions(-) create mode 100644 model/__pycache__/ldpc_sim.cpython-313.pyc diff --git a/model/__pycache__/ldpc_sim.cpython-313.pyc b/model/__pycache__/ldpc_sim.cpython-313.pyc new file mode 100644 index 0000000000000000000000000000000000000000..94e7d624b36921957d579ad0631aa713063fac22 GIT binary patch literal 21540 zcmcJ1dvqJudFKq?Z$7}cNNOmN7C{kwhHn5K@y~3lK?pbN+LFC zm^L{D?bL#7HwGPT1ZBM&Dp5B~>fI`*=~iyLjeDBz0SzUh5uP5m>vK~7v(@Wt*7;+1 zf8QMp08+G)&33QEJ9qB=?l<@Oz3zAMz;3s22)(WUCh-1Ij{9%)q6Suimj^HL9QTKu z$Vr^Yi<)P6NyDC6Nz0x(NynagNiS+8gQ$~?qFyqcF`O~E^&BUe&zL038Dr+nDjFo4 zXgp&RO>TU_ZQ?j3CE3^JHgS?eG)V<#Gy{U#q}+20*UE|Jejc%f#a0&ESZrspgT)1i zjli)bS65Amh5Iy8u~<~aNlvjCaf#?eTq>3zE<0n}$%~~`+{?UFE|wvuLM%sId8PpQ z6)e9>+$>hI)F!bCW!0h}Zes5>Vx3scQng|YN;aPn#M;h5j`N|vXEtx)c4ZnNaO0YB zjcbc<6@zR!dF1qA;fQ}EH0qaxL;iEV^MR1$n-F@DF(J5xLxHgS=tbX1SQzk+`6YjF z#4n76Bw^s$lfrX;DKHio@r47Spv7`xa%#dq=?{iwAv887`@_OyAn2B7CWTQ#9|(>M z7XsmP!h~-Yct!t3P%u$wiD8YMIayLRl_-7LJ))v=>l zc#*}!KG`o!`odD+qOj-UURR67a`>Dt7(|OFsokBz=}E}K*QY27hph!azrk`_}B$mdm;XeJ} z2omkVjUt!nQ!4>!shsL;1jON${v1dJe6tFEolomW!?RQV5#(e84J?hugY%494P^@kXn59zE!+?d1q(wo z2Q2K=w?Jd(*(I88h93lsB0|jqS-@may=J0F39J>IFu5O%YzvPjAn?Yg@vJF zH8ohgQ$v%l!UnY-t*GflxevYQt>yB+UZU1gq4D*0VAka@do)8?D5AxKhux@@mxsPn zj*1}&3O>l^m7|VFM-U$UPR%wvdpkJJ$_ZTWXyzGFo*u5AbLk?7k6xVe1xGHLVFwrS%cN z6q3_?Z<_B*>jJ^>-n3rw1;_nqL$CKxw|F#d=u@6r6s8S5$_xLZ%aG=K(!7-BN78!v zTwpAm<|Aq2NNB?^xBeXR9dXwu(0;yc$W?HJwBN))ekY--%sIx^v#8TEUo zq|nPC8m){LWoTC`qj5@>LY~sj;hDh1sF!-|J?B*@qtt}&bmC8b6VV*^iOF)w^s4Et zy0`jcWmlxD7Zxroo3_n$uiA^>+Wq?3D8E`%_SX5Y*F?KlOUkbp;?wc)f-P0j5Iwq5 zP!gMdZT6$Wicd@8)-CuG)7NI->%JCRuHBWY+8r~lG`eE^TfWN; z@w&^d`18r)`cz@V13gz(`@qT-Z@#DH9GmVIb2j^2&jU1z8_C}*G%{>Qy@SLbazmLl zjt%?_S@2}nkQ?=a8V_Uy#_tLN&kfL`UKdmjSkL8l;J^mEZVNZ4Y2rkU4tPZ1 zZQ^uXt43`BR-HD|d3bhYA}}IIf$?);H=AR@C%ovM@CV1i;}kkB49^PNSlKpVDkKMp zPs%N9o%Ol&l7Qx=wWyKSjfAFVr3!ktj)a2XsQw`ST-vl2n4~ADfR|jjBu&tS)EY151%-fn7%Ff16c1CG7 zM)6k8`FWAA!W@Yj#M+_k+~np#9bG5)oKDuXa1|f1M^GQsc{E!v!$VjvN)UC_uw9w<2Mr!A$_;tt`uuW_UNjJW%JoJx zehv#|VR|61|o+LDL#a zZyv04g9m0>Z^>&znTkF4L>nxku@mi8S4Lx|a1{Hb9Fn-lGC=%%T~PP9gX_?!X+Q5V z*P}hGO}a5PWg6pL7GE!h-y-^_XCS!-AQ(czRfJXl#i@x9_AI2v@HxM5VxU`)!@htt zLc*sUmS#r6Gm;;oXHHxVo+l)sV5gw)f+6Yb(8yV#OOS<0zb_~Y^31S%cp?NI6p)1x z;DX#33=5~dvgN`#NX3*Ad{a{s0SL%oVHrAdQ5X-L_sgAvWMIX4(`XyQDRiU&d$gf@uP-hX|rt?i3Vf#6t3*si>`f3az*9qnlo1or`fB^g-z z7hiI>3r%Vrq^K4|Tf!wNxTLkXBnh}MXc^Q}v$UPC58**QZewyWeMHTn4=GT}w;EZ- z#z-OL)95yV7Tg#_bpSJSA}J&3q8h;5G@DmrwppNIM-4{&)oILD zSwk;r-#GF2;0Pkn%j|01SVjajqD{m?j^{OHX%cu3A8GE!0-066t1O=lmqx>jXIrjo z6~3aQ-ZGZs(xmx|62vy_14Jr6G19V|^-iU$nXJLKfxl%x7OQVUKG2dW9j0ou?Z9IS ziY{Gv^+Ig=N=JP98@uD%QUzO*rY(>1jO-X3LcQO13{<+qMxaU=0h-PL_M9a9*>ZBr zUu=qiqw4b1OZ<>cN>l+!FnUv}& zyMeMzBx_y>Nuy~!lZs_hGO()Ephy|kys?sw1EL>)@<|A0bKJvP&Q=heeeJ-N{*<*g zsjE#|Ywznc`mOg(8vT}s7Oth^z1=s?rc6(M=jZ3IUYFPiwcXuz=E%qVnwipt{UiK1&I@19@uEskFMLaK7-a?#G4{Bq&G zsQx2oMZEM~Lt=W-cfBETaam|fIoqRr)HmPon(3aIbJRX`aD~N}e&^NSiO%UbiX)j4@A55U+AIvrIl$a{$bIeemd1mK@aQ)!tAz^W?WIE2!l#iW(@mhM@vA z$D>toLvZRQ!ATK_f+mRf+RnG2mSPDHR^)v*E54~XI1XY|rk`%CY84}?qhtUa9Nz}e zUQj@fDT5O|7j>EZE?iY*Ep_!taZ ze&5JBsG(spF5%JD-YJBkM#__curM_<0bv3h4Z;T0JTQI5U?qbpO!y{2-?p>oaHcL; za4__1YYhvH!@iMM#wA$2MukvtVivq_z(3-L5Fx?x1zVV${ZON7DjN+zs0j?upcN1v zPR)d;X2PAq$&&-rme6?P>BiYFw7Y~>ff8*l7qd38&PfVtnDk~RheH#YerGiAM$(gB zX$^+lEa8GRFwD$gG#;mwQo`FMw8PwpacFcmuIakVl?T|i29V%}T^5DrmSR#FO&(Vk z!V9-tW=2;P2B>EOR@#7}9}P{WjZ=z+kdfH5acm;w3-8{^=r44#&^YQI1)j8?wMBUh zr}-cXeUoXsGEH7(ERi<>Wj`Q2i5fEP7?p@-bhU!a$GjuT%8)#O`3eDl9}%OZW?Qu9 zHTzsQ2pvdCVwcNjW-to`-DD8FK>1bxC2G`0P;eJiC^?-qZ%>|14vrPph2HY_(j zm8v_qTyt=_^kDSxYsWy-3X1?Yvr@S^rhU`!ab<0+`%S|=C+8@;TLx+m>Xp>hJj&up zXkj$Uz-W}ndw5ZkHMX-9l)p8(O;7}BPty(s6#@b!Y%CZJ>S*_>J9#7JJ`ckNbWE)Y zilu&t8z5G_E`V}~IqLXpzOwB%VG6p)9)Tb%QIO>PX$!%z_WiQXV{_oWb$T{=!34uW>D~r(u@iH zZqRad2FbX1ZqkpFg_f*chmCx{xA#O}T8o2_v>xg1L6~%SUz#6Qgi~r2NQE-^X4x0^ z0=suwdKvkDhClfUM2a=1`!&N`&RGAOW%C6v*;yT5n|W*GN<*UP&B>ItKB=pJbik;i zY{06J$RDslew{+mn%6z(hTu*(?s;4vdDcRRz2|^h3LsjeO0b>)5qd(C^_Nsm93DNf zmc;twnv``@Qn%^RCey}cP5zLPk-YiLZhhGKpin)h&DtIyRD?Pp++-})qUIEAO|U}1 z0%YX`kQR`cwKY%zfr-AuLQdpBEeI2=mOkG~l@oYCuFysZELjp0H0GsEA|EslWkn{e zCiH^>gqf|VXBm8vD5b~9S_|^wF$PUa9xOzbpt%=JECXXA7_y%ALLJD!n0gzO_O)_B z>*svu0KQ`z{2OKoMQxAi6pydibLRshKe$Y!SzX~l`{1a@M(uKAEQ1b@qmi17;KLYs z>`fdn*nwsO!&yPI?Y9C>tKd7|MW6fGyv{`ElQZAjZ0|R0?wf1;WG0F$&uy1B+ z6n4*i!70S#5aDH#7KmO#P8d8jfD=Xm%8w9&Wk$kMg{Z8AQ=&wBZt)upCQ?@(8%ZjNhm@F zU8jn6BqTvFB6E;%-ZwGhClL%(B$$sKz59*OE4n6;DaIsRMkd%m0UVYk3_2%~%%doe zqey{KOilTNqprfVG1H5*LFsK;GdhyyC&3Lcm~vW&0g=-tdhuQiNq7$g+gZFLZK5=t zSV+$jP7}j~7N_}_LDGEVP(8;X52151$>1~u`9Yy-ML^<(is~H+dX?Fgz69V8@F)LI zL`?Xx*rU(CRy}vht?fI;M0?ElX2X1cy#GUc?RDea(WI^Iqw1}obnhIGO@Hlh zjDP*u!tuEyNo)OTL2+!?Yms>O{O`>5Kt!??%^g{_7S0`6DJY!lAvtKr>o1V}Q!u|L z_Uvm1R-EP0BZw~V!K1w5)t(i5(fsqV@E^RyGG9G%&(1l@?-oKC%&%Tn{^sQT=weA! zyXq*4)y2fiXHt$WQSFMOEIz$ZvFxaeYX82hGAhrH#${xi3NIDR7sRzcGgaTS5<&-O zD^k9L1z9Cy4xCpr;x7Q;IXQ?YepHhXe3}(&yE;Lw7?C{ z$Z*@xd*YG`mMm_t4s7Wekn*X&&oRWjR==Jx_h(i_uprl#+_!RV>e+{u6nY9h1tc2f zAQ)pS@;LHx&^}KQK?REk&*6phI9OaranPA}=&SP-vl6GuR0jV}f!qejUZ_RX!`?CF znQPx~Jys?8D<#+Byk0&-z5JXn*|RGEb01pu*mJ$Ye1HkE`6va`B-1Smo>CQeuxw+j zJ+v?Fp0W(K5)~^Ofs9zu+nX7IjpsS|FFMcIRsZ+%?8&y9&2hOpkL;~t2s_ed?vBeC;Y>JwKY6`zhslobUaqlJe-VKPs~`LHnc9Q@Md4r|yrcVC4{*W?Ar5 zW#?m)N7IK;(DgXi|XN0j-aJx%U6IEn|J`L0@^#rk%PCL@s?+_}VCCer%~dTUwc)f(`87 zfzVh^W}LXzd3rIoE#~e8v8)jqlH1RG$E5YXV{|-~le9=Bo$eG=?Q>Mf*hxrHhK$vx z9;^E#$u~=qSVnVJlrU8oQ=~BHUEV<0j{_#8V5yrF=#6=2k#s1)9w*_P19_9M#pC2l zfMF<@v&!ROj{H8>Noz`07vL2`iyvhker!WR5mx&x+CEJn{ z8vCAg6|l1(MMPvL2<(sn`yz9yKnpx=_T%f!C5ps~Wa4noh`znVCpiazifdq3;-n{W#$qlchuSriEaM^9yyB(po8WK{=yvVDy5_ zE(%NAge*w}ZaN?IhC`4i)0&Y$T8jeaAR!UCNt@XTsDIQ8e;X`6l3vrs5L{#E8z#u0 zsChn^HY3Hf203j(Qc*JGv@z(vz#z^3i?D?T<-iEjZqquN@pQqu&dWqD@;RFPPqMsybyI!f zrJrrudD$4#tX6EgaxBr2*t=ZaoT_lg^sAM%SEj#N5$j&5+mi5HtGwa6RuDU}Qd{?K z_rmz1dAYVdb{uaD=MwV5JzO?#+%0Fwu85}rJi^?)^qtJ%p9j4n7OJg z_w-!#u7|~sZ&21pY0v(f1MkhIs=I$MyHuPQUO2N@cg6et*_b(5-Mw1VvN)Tn*%z~{ zR&7pbzUhh`T6v;vsrXtj))ODdd{WD8s=sH1?-CA&svG`ayZJ;yT9{aDUkIVkLgO{t zVm}SW>bAC}%A3bi+m6JJ$B!?U)URx7i5FoN@)^*Rrb7c+p26+xzn6>p=7ZP=` zB3PVC=ey(GAKGfJT_A<4b+ydBIGQT!m^(!_yuGjPjU9;>U#W}rq^vdZFD=zAZ|+zc zy!oZ1wJ)jb`|N&z=L(+q?7oS!Z~5#)>%PxEwDf;=-^SIR;y=4@;VO1~cE5z;f2{8Q z?EV(6a3>Yo+GMb9?yIFIOLc!#ezIQozOMY_X3P7VY$vM<-`{4Wbc+%Bf6-KY(q;Nf z&OotHa@fZH_M($}Eq~Kiadho!>*mfwSg|ug9BPfvRQ*s zo>bqkPO3L}3_00|+%gAu6Uol$c1(UZdh~fY5ZXLOf`XChpyHn6hHMnn^@36>GE7!- z(m+l>OR@xH9Xnw7N8BxWm_oIY1BVqxu%6q%*plO1GeOug1M(Pou_&IG;OOwY%s9U> z8#hebHE5c0oXz+U<46MK)gf}3(`KT&iYeQm@HFN8no)X#f|xN1g#lcmsxd_Jqen2( z?@}Q}otP2%50HRG+wPbKA=&@B$k=B2Oa}-?7uDm7kzzQ7=nVM1a&j0PxxA8kAwv-uT|Qd9~)Yt%X!D) zO`m(qLD6B5Xei)BW6mQO*iGt})?hbzc;@Am*N`~F6^v&GQYsIfDIBfh{f23#nei=MG`cM&Y%f@X8c+3XU$Vk@@>Cw1J;%|4_Q@2 z9Z`-O*py4FDbXpGiI)-TfD7n9Y@g%c0$Z!xWfGSLTiiMoT zHdr*YHJe)mLlR&q%wQb`i-XQ!2^2`Hr-Ts%YhL^76}bEu3RJPj-Z`E-qzjgMOrBC! z*W|I}`oVH6o>I|C$N9{{{RRk+SyW)ny^xC-al(kQ1Anw;@KqXfSoWR4vcWglkzx7Z z57@IJx85)}Wi;DR*vi?nOsrtfa1*ond zsH|mrkit!ohn7b#fCwQvV)qHmPh8bhAe$@a#hJ7BwYm=-^Ls~!y{rmU5}0DM735ZV zBE-&OMh^*%?QKmsnQm`W&d1rgxS+TPbi(&Nd*Gk>7Tqji>O#;vCiy1GeK1F%Rp%f0 z0ieBfgM^(-GI&*P?^G`Cpgb#^%a&hRdNTg=okGOiOdWj%wEtJe!?53X3XN&LIbw7n z8`oNXrETdL`<0Gv?L9!frEZns+ptRu5f(fC>t84kw+V=nzF; zA(Y09aty6XRd#VgjjA|az?Wp{y99cJA{sm9KDb^}%Dzi3LRl=$1ict>x_*_pS|!~= z0Xd$@{|gb04N5q3b&6cuS9Rv-)3NEdFI>J53&dw{6o0SsdgUG0fo0c$o6~=O@rM_G zqW{w`rCg_yji;dw-PR3Ih2w8^$2wxZD|}3jZ+jQ|$hEDBkwoQU{f+6xk)&nEZCwYz zk~SexnlL6Ujkk3!s`g@Bm$cQy`{Nf9a(oi*@}8xgizjdLHydv0ZdTvwzZFRiJeNH7 ze9|&_TlWIhX#C2lj~#_`$5%_XL`^G&#j(MJ=JMI)!mV&wgP>u-Aw<-sr0A7tYhU+9 zkA4Dyec$DMcbxUh&iX|EB7dVj<#b2M6Ru)wVrQy+dz8*Qs<$uhy?HdX>CpVi=!w|$ zYKf4jOqFbp9z|L6A^ujw2NkL91IxC?=<(PK34LOE#a0$?Shj7BJ`ISHvbT4?0oTpM zvkR>$XA_P@P+GritBpR5)RnsVQ%DpzW4=GQxKdc1fG=I$8@78c=m~e5xu%^F|LaLS>4R|_hmYt1J(YV8m^|%C9(yk7c>cC&aK%(`$uaL(sSsA2WpR1QbWg9Xu>6wKnyvSa@?7cGyL~)o zbI$cWDlHXZaS+FQ@bAIHHkOu%m*??{tbvbl_nqK)hMed;kQLT}GZQv)g%cqz=Y)2Z z=i;CuD?bf}Rj!JI379gKF!<%1lL>L?sdw^>=~<}+@e5d9)rg? zkTuAvn8ANDNssBVrIyE*njc$geQc@iaiw_@l3nGk9*3u32#Sdk@@8Q}L+#*?^Q82c4mFc`0CXtll4Ivmoa+52%GI6j z7R89+SoevO@NB27h(RFo?TVNMT=k`!t3pHru5T0?2}4)J)(n9rW8Uma^Db8rs zYY{TzXpH;i^E|bX4y~m>r6|j^;PuDu1IvtYNk5`O8WQP0BT5@!-t%4vOKr+a7(PF; zT#T12bIK+ZE>(&vPN3ZSd5)YM;nY2l?*jeDhX%biNs za-}O(xqYQ(%hfXrXTI^$O7+&OeG7f5>gLtTnt0$Fn^&qfuWW9>&AuzcSI6ENOFVzg zn_{>6?v>jLEK%*tnRCD#fI|mlG-PyEvdDGs@Ww#8s@+KEHJGR^=`02kl8m?;TU$JvW@M@3!U9 zbqrrNKdg>*F8}y>O|by+NTvh55fot(_jFwP zb>!)4Ls@NFtlk*c#(0ohy$&^A#hQIOnEbN*VnBElW+1Cc=gNLfjfNbgni_gSO{pU% z*5Vf*gk2e~1-A~8k6FMy?K+iLX~-JXq?QactEt=*oULBmB5tL%FS37XUGT)L6d~3R zwW@V;-fI}$29I{Fm1i8L!U3X`>q2%$KUg!PPpJLD*K%uX$GXaR-NK36#75DTm8JAU z>ayep`!||riLjTU+ie2SpYbu^qOR9YgzbuxspfPpSFpG(tMgeX#>F>_HVG_r%wC zqU3NcM@^|%Hm)b`BKhod%63zk8xd>u9$13%t(Lh}xe>)Xs(wHB6nBYRa$}NRuco}- z?ah^=NBcacwR>pjxQcD?_to792cz8Rbmm4H^`G>ZA2%Xvaf|zLRPj2r8cW`O+o)&1 z1uY*sk!|-o*k>EiV|6rUNVoaX&@(E0&fy2nHNDx~tR}#8n|}@#LblG$#~~Fi_XIWj zr+h6SSe!R%UAQ@X0BiCoPa#I_pr-(Fx2G86lV?ZnjVY4&I9JCtb1Kyug~T?d5xC79 z4rAd;()W6ob~87|KL&>H`%?wRb*L{=^Vpx$aq(&Es9(mJGJb?66pYx{{3=Y^0&fJ# z@5Rkxf21<28PE^yyI&u;9BFK&^YQJ^WlDRr_kQkh8sTGw2oajfPshQOFzX75HeT@-J}A1N-Qx_ z-|M@$fvGSLT}PWMVsyJ`%kxqQhSP{m{UspsN{;NO-}Djb7k)nfhmA5GGh~om`zuO;U(v7nTu0+FBY*cP zmsU}l-lGwUNOP&X?m5?HoBvMQkJJEtaOVjZDHV4A~Q=cQ-0Gr zT48A`G{vlbNP2AELP^pq%|cp(Uz^cQ`NHSYx|eZJ2!_>w6eQcL$?F}%<#g>ndL?Uw^mi0}g`%~}w)9_-mo_Nl zEj#g&)va?0MazbZdLStM5+I8_lK_R>4{`VZ~($@bIB_WrqJD>};^UHP)E zJYKM5TdCQ5SIg-u?s2@n^W$uZeQ{!?rXyRjD_3H>(Q#9^Qgb+4c9dURsQJ7?^9}uy zyizlet@&(T%}n9ZY~k^|LesnQ(xJ62wHDp6Rxewtah-bZ$R|3(@1OYciQn)0a$jsG znSyWYYCbkv?-U+vHg--nLR! zkeEaO>voNVTEm0$!dSLtAb9v9(`!DZ*SMzP-HRHDnt`)pz zN(uXJ?z~l&JaH!FJPR9yy(rcgKYZIRB)V5^g_lmedLlNySoEQ7`!DfJK6PZ{HC=JW zpZ`X6+>w}A+Pl=0baegHbl_vN{Zj3vT56`_#@9g#8ynRC{=b7K6eaXVDiJnAf($Vr$Q!Bb{b|f7;ZksyDG+|`s z|GF#P_zj$M?@r$+x+X7nBm#@!_qvxll8!yMO?yAFZH8+V%ryYF#M|ladfb(Ga&h{5 zGuLOnd*LSk=f)ozZ}C60+-m>vp7-~B(Efhs^4{L%*1iwqWQ(Xw2K`dPhNd3iv|sVX z#e_C-B(Xc`Xu56M{z=aML!Z3T9^aa*aW9_x-o*8Z?*^0ZuB79@ZPQZ^VDxE<9$mFL zE}fV^aj9>VMqyep9NhmyG!*a>oLzZO3lSq}q_}jvi0w z6X#cPsr0!wo{Jw}IF%}GhU0^+2p&+LL`VFY#21%5OULhW{C*gQ8so1lRxa&cYP_rE zo%{a>$J@K^!?mI5mrr3G{LC)=>X#LngiHeAUpvZ9Z`b}!N7?BvZF0Nibf+=dW$CvW ziI&lZ_Ih!B>>~g^3vTptJf)veL^bA}bO1(D7$$AivA`@j4&CG`1k48UZL z9*ecl_srYl{IaQXPQPNYM#p0P^XKN>NmB*h3G8k)XR*(@g?`i`0ae{-7`E_-}I0(@m=@Jw0y&Tn})AOgSWXK za*ce+!{_)TJb#M6-=`_yiyzkS p1 = P_{h01}^-1 * P_{h00} * info + Row r: P_{hr0}*info + P_{hrr}*p[r] + P_{hr,r+1}*p[r+1] = 0 + -> p[r+1] = P_{hr,r+1}^-1 * (P_{hr0}*info + P_{hrr}*p[r]) + + P_s * x = np.roll(x, -s) (left circular shift) + P_s^-1*y = np.roll(y, +s) (right circular shift) """ - # info_bits goes in columns 0..K-1 (first base column = info) - # Parity bits in columns K..N-1 + assert len(info_bits) == K - # We need to solve: H[:,K:] * p = H[:,:K] * info (mod 2) - H_p = H[:, K:].copy() # M x (N-K) = 224 x 224 - H_i = H[:, :K].copy() # M x K = 224 x 32 + def apply_p(x, s): + """Apply circulant P_s: result[i] = x[(i+s)%Z].""" + return np.roll(x, -int(s)) - syndrome = H_i @ info_bits % 2 # M-vector + def inv_p(y, s): + """Apply P_s inverse: P_{-s}.""" + return np.roll(y, int(s)) - # Gaussian elimination on H_p to solve for parity - n_parity = N - K # 224 - assert H_p.shape == (M, n_parity) + # Parity blocks: p[c] is a Z-length binary vector for base column c (1..7) + p = [np.zeros(Z, dtype=np.int8) for _ in range(N_BASE)] - # Augmented matrix [H_p | syndrome] - aug = np.hstack([H_p, syndrome.reshape(-1, 1)]).astype(np.int8) + # Step 1: Solve row 0 for p[1] + # P_{H[0][0]} * info + P_{H[0][1]} * p1 = 0 + # p1 = P_{H[0][1]}^-1 * P_{H[0][0]} * info + accum = apply_p(info_bits, H_BASE[0, 0]) + p[1] = inv_p(accum, H_BASE[0, 1]) - # Forward elimination - pivot_row = 0 - for col in range(n_parity): - # Find pivot - found = False - for row in range(pivot_row, M): - if aug[row, col] == 1: - aug[[pivot_row, row]] = aug[[row, pivot_row]] - found = True - break - if not found: - continue # skip this column (rank deficient) + # Step 2: Solve rows 1-6 sequentially for p[2]..p[7] + for r in range(1, M_BASE): + # P_{H[r][0]}*info + P_{H[r][r]}*p[r] + P_{H[r][r+1]}*p[r+1] = 0 + accum = apply_p(info_bits, H_BASE[r, 0]) + accum = (accum + apply_p(p[r], H_BASE[r, r])) % 2 + p[r + 1] = inv_p(accum, H_BASE[r, r + 1]) - # Eliminate - for row in range(M): - if row != pivot_row and aug[row, col] == 1: - aug[row] = (aug[row] + aug[pivot_row]) % 2 - pivot_row += 1 - - parity = aug[:n_parity, -1] # solution - codeword = np.concatenate([info_bits, parity]) + # Assemble codeword: [info | p1 | p2 | ... | p7] + codeword = np.concatenate([info_bits] + [p[c] for c in range(1, N_BASE)]) # Verify check = H @ codeword % 2 @@ -140,18 +155,22 @@ def poisson_channel(codeword, lam_s, lam_b): # Compute exact LLR for each observation # P(y|1) = (lam_s+lam_b)^y * exp(-(lam_s+lam_b)) / y! # P(y|0) = lam_b^y * exp(-lam_b) / y! - # LLR = y * log((lam_s+lam_b)/lam_b) - lam_s + # + # Convention: LLR = log(P(y|0) / P(y|1)) (positive = bit 0 more likely) + # This matches decoder: positive belief -> hard decision 0, negative -> 1 + # + # LLR = lam_s - y * log((lam_s + lam_b) / lam_b) llr = np.zeros(n, dtype=np.float64) for i in range(n): y = photon_counts[i] if lam_b > 0: - llr[i] = y * np.log((lam_s + lam_b) / lam_b) - lam_s + llr[i] = lam_s - y * np.log((lam_s + lam_b) / lam_b) else: # No background: click = definitely bit 1, no click = definitely bit 0 if y > 0: - llr[i] = 100.0 # strong positive (bit=1) + llr[i] = -100.0 # strong negative (bit=1 likely) else: - llr[i] = -lam_s # no photons, likely bit=0 + llr[i] = lam_s # no photons, likely bit=0 return llr, photon_counts @@ -246,37 +265,37 @@ def decode_layered_min_sum(llr_q, max_iter=30, early_term=True): for iteration in range(max_iter): # Process each base matrix row (layer) for row in range(M_BASE): + # Find which columns are connected in this row + connected_cols = [c for c in range(N_BASE) if H_BASE[row, c] >= 0] + dc = len(connected_cols) + # Step 1: Compute VN->CN messages by subtracting old CN->VN - vn_to_cn = [[0]*Z for _ in range(N_BASE)] - for col in range(N_BASE): + # vn_to_cn[col_pos][z] where col_pos indexes into connected_cols + vn_to_cn = [[0]*Z for _ in range(dc)] + for ci, col in enumerate(connected_cols): shift = int(H_BASE[row, col]) - if shift < 0: - continue for z in range(Z): shifted_z = (z + shift) % Z bit_idx = col * Z + shifted_z old_msg = msg[row][col][z] - vn_to_cn[col][z] = sat_sub_q(beliefs[bit_idx], old_msg) + vn_to_cn[ci][z] = sat_sub_q(beliefs[bit_idx], old_msg) - # Step 2: CN min-sum update - cn_to_vn = [[0]*Z for _ in range(N_BASE)] + # Step 2: CN min-sum update (only over connected columns) + cn_to_vn = [[0]*Z for _ in range(dc)] for z in range(Z): - # Gather messages from all columns for this check node - cn_inputs = [vn_to_cn[col][z] for col in range(N_BASE)] + cn_inputs = [vn_to_cn[ci][z] for ci in range(dc)] cn_outputs = min_sum_cn_update(cn_inputs) - for col in range(N_BASE): - cn_to_vn[col][z] = cn_outputs[col] + for ci in range(dc): + cn_to_vn[ci][z] = cn_outputs[ci] # Step 3: Update beliefs and store new messages - for col in range(N_BASE): + for ci, col in enumerate(connected_cols): shift = int(H_BASE[row, col]) - if shift < 0: - continue for z in range(Z): shifted_z = (z + shift) % Z bit_idx = col * Z + shifted_z - new_msg = cn_to_vn[col][z] - extrinsic = vn_to_cn[col][z] + new_msg = cn_to_vn[ci][z] + extrinsic = vn_to_cn[ci][z] beliefs[bit_idx] = sat_add_q(extrinsic, new_msg) msg[row][col][z] = new_msg