package org.bouncycastle.pqc.crypto.hawk;

import org.bouncycastle.crypto.digests.SHAKEDigest;
import org.bouncycastle.util.Arrays;
import org.bouncycastle.util.Integers;
import org.bouncycastle.util.Longs;
import org.bouncycastle.util.Pack;

class HawkEngine
{
    public static byte[] BITS_LIM00 = new byte[]{0, 0, 0, 0, 0, 0, 0, 0, 9, 9, 10};
    public static byte[] BITS_LIM01 = new byte[]{0, 0, 0, 0, 0, 0, 0, 0, 11, 12, 14};
    public static byte[] BITS_LIMS0 = new byte[]{0, 0, 0, 0, 0, 0, 0, 0, 12, 13, 14};
    public static byte[] BITS_LIMS1 = new byte[]{0, 0, 0, 0, 0, 0, 0, 0, 9, 9, 10};
    public static byte[] BITS_LIM11 = new byte[]{0, 0, 0, 0, 0, 0, 0, 0, 13, 15, 17};
    public static final SmallPrime[] PRIMES = {
        new SmallPrime(2147473409L, 2042615807L, 419348484L, 1790111537L, 786166065L, 20478),
        new SmallPrime(2147389441L, 1862176767L, 1141604340L, 677655126L, 2024968256L, 942807490),
        new SmallPrime(2147387393L, 1472104447L, 554514419L, 563781659L, 1438853699L, 511282737),
        new SmallPrime(2147377153L, 3690881023L, 269819887L, 978644358L, 1971237828L, 1936446844),
        new SmallPrime(2147358721L, 3720222719L, 153618407L, 1882929796L, 289507384L, 264920030),
        new SmallPrime(2147352577L, 2147352575L, 3145700L, 875644459L, 1993867586L, 1197387618),
        new SmallPrime(2147346433L, 498984959L, 154699745L, 1268990641L, 1559885885L, 2112514368),
        new SmallPrime(2147338241L, 2478688255L, 826591197L, 304701980L, 207126964L, 842573420),
        new SmallPrime(2147309569L, 1908234239L, 964657100L, 1953449942L, 309167499L, 75092579),
        new SmallPrime(2147297281L, 1774004223L, 1503608772L, 353848817L, 1726802198L, 2084007226),
        new SmallPrime(2147295233L, 1006444543L, 279363522L, 632955619L, 419515149L, 38880066),
        new SmallPrime(2147239937L, 2881243135L, 1383813014L, 710717957L, 756383674L, 706593520),
        new SmallPrime(2147235841L, 867973119L, 848351122L, 1627458017L, 1867341538L, 1260600213),
        new SmallPrime(2147217409L, 4277923839L, 100122496L, 1700895952L, 1076153660L, 370621817),
        new SmallPrime(2147205121L, 1878769663L, 1111621492L, 44667394L, 1556218865L, 32248737),
        new SmallPrime(2147196929L, 1543217151L, 310009707L, 1302747296L, 1775932980L, 1745280407),
        new SmallPrime(2147178497L, 3518715903L, 1006651223L, 995615578L, 1904695444L, 546877831),
        new SmallPrime(2147100673L, 1505372159L, 520918771L, 32821017L, 801131396L, 1402329446),
        new SmallPrime(2147082241L, 4227457023L, 385646296L, 671813454L, 349022278L, 1235641740),
        new SmallPrime(2147074049L, 1878638591L, 1198259916L, 22548867L, 381743482L, 316764378),
        new SmallPrime(2147051521L, 96036863L, 1908098729L, 1744576193L, 129309952L, 1045517086),
        new SmallPrime(2147043329L, 1538869247L, 440645275L, 848789527L, 259769898L, 74455690),
        new SmallPrime(2147039233L, 2209953791L, 2055370389L, 1117428534L, 1471147502L, 119673623),
        new SmallPrime(2146988033L, 1324904447L, 1363381819L, 2041738204L, 478044407L, 147722658),
        new SmallPrime(2146963457L, 2130186239L, 325123596L, 2068278376L, 875801492L, 1306038870),
        new SmallPrime(2146959361L, 2146959359L, 264240644L, 872271011L, 507971579L, 1990286186),
        new SmallPrime(2146938881L, 1727508479L, 1974074844L, 198886428L, 561977686L, 1269843026),
        new SmallPrime(2146908161L, 1672951807L, 98813339L, 1770745661L, 1393783074L, 1823212183),
        new SmallPrime(2146885633L, 1006034943L, 661929322L, 1068953529L, 641063610L, 596476372),
        new SmallPrime(2146871297L, 834054143L, 1378823498L, 1788896577L, 959521530L, 583621658),
        new SmallPrime(2146846721L, 196495359L, 1834738961L, 827307343L, 2002269931L, 1831804644),
        new SmallPrime(2146834433L, 1572214783L, 655434995L, 1825957452L, 2005158872L, 95719048),
        new SmallPrime(2146818049L, 1505089535L, 963224779L, 1698458806L, 1152339536L, 1292562110),
        new SmallPrime(2146775041L, 2532651007L, 1260858461L, 1897922218L, 1980311657L, 1048890741),
        new SmallPrime(2146756609L, 1840572415L, 1934326828L, 1319016594L, 303348307L, 309425167),
        new SmallPrime(2146744321L, 1001699327L, 1031932938L, 166920536L, 125412453L, 206010953),
        new SmallPrime(2146738177L, 469016575L, 1034034169L, 2049530115L, 240883589L, 2109942428),
        new SmallPrime(2146736129L, 1706334207L, 386311155L, 1983901927L, 399193796L, 1866844306),
        new SmallPrime(2146713601L, 1878278143L, 1917713332L, 795783646L, 2145730902L, 666529419),
        new SmallPrime(2146695169L, 3216242687L, 923528062L, 410549424L, 826761640L, 881662233),
        new SmallPrime(2146656257L, 468934655L, 1316739849L, 819155176L, 1775281069L, 769010182),
        new SmallPrime(2146650113L, 3149088767L, 1357138678L, 1776223808L, 1272363648L, 1075150578),
        new SmallPrime(2146646017L, 598947839L, 836424425L, 861255062L, 862632924L, 1452282057),
        new SmallPrime(2146643969L, 1458778111L, 1699760867L, 1128852454L, 1250440960L, 551264863),
        new SmallPrime(2146603009L, 4008873983L, 258845279L, 1693343820L, 1370713748L, 896492399),
        new SmallPrime(2146572289L, 498210815L, 1553576441L, 1357574155L, 2042997788L, 475091953),
        new SmallPrime(2146547713L, 2343679999L, 731429284L, 874518801L, 933977652L, 2038082698),
        new SmallPrime(2146508801L, 1005658111L, 1986115866L, 994764483L, 1486111048L, 77164992),
        new SmallPrime(2146492417L, 3153125375L, 2074458334L, 1607832505L, 1814950582L, 1064849124),
        new SmallPrime(2146490369L, 3383810047L, 1163389142L, 1860189331L, 766301042L, 1792895956),
        new SmallPrime(2146459649L, 1542479871L, 113653858L, 1774578485L, 436467466L, 509354151),
        new SmallPrime(2146447361L, 1995452415L, 521816115L, 984131366L, 920695297L, 1917290506),
        new SmallPrime(2146441217L, 2108692479L, 106094619L, 1700403904L, 55841604L, 574542256),
        new SmallPrime(2146437121L, 2142242815L, 2143320076L, 1575554072L, 1350403963L, 1873207947),
        new SmallPrime(2146430977L, 2129653759L, 84969459L, 244551713L, 118889817L, 571880191),
        new SmallPrime(2146418689L, 1877983231L, 1167996867L, 1775456464L, 1211626453L, 1434078655),
        new SmallPrime(2146406401L, 1324322815L, 1313757074L, 1622299552L, 60879896L, 1023260741),
        new SmallPrime(2146404353L, 1202685951L, 1813342090L, 1162222118L, 193969494L, 1214240996),
        new SmallPrime(2146379777L, 3383699455L, 1842644774L, 237716594L, 1996303667L, 1647718365),
        new SmallPrime(2146363393L, 1303308287L, 1687447266L, 1850889128L, 545315431L, 1827479741),
        new SmallPrime(2146355201L, 61786111L, 269635263L, 691559024L, 756508884L, 1699936240),
        new SmallPrime(2146336769L, 1072594943L, 654341745L, 1871443110L, 107703592L, 1164751303),
        new SmallPrime(2146312193L, 4226686975L, 1106990599L, 230828908L, 2118815775L, 547343662),
        new SmallPrime(2146293761L, 3652048895L, 1401349558L, 1527866384L, 1939695340L, 934813966),
        new SmallPrime(2146283521L, 653111295L, 594294152L, 1585987504L, 113384110L, 901617467),
        new SmallPrime(2146203649L, 128743423L, 1230019607L, 2070978911L, 1602264262L, 195465499),
        new SmallPrime(2146154497L, 732674047L, 1428919080L, 29335864L, 1593288655L, 2129713526),
        new SmallPrime(2146142209L, 2276165631L, 1819536107L, 649122195L, 996268341L, 667538093),
        new SmallPrime(2146127873L, 2263568383L, 2015437475L, 207440004L, 1904659130L, 2030659201),
        new SmallPrime(2146099201L, 1005248511L, 907637264L, 559331750L, 152292400L, 1381346740),
        new SmallPrime(2146093057L, 1135265791L, 453939696L, 1717051717L, 781899671L, 630273220),
        new SmallPrime(2146091009L, 4025139199L, 1800630758L, 1086157903L, 1771609892L, 1947078850),
        new SmallPrime(2146078721L, 4008349695L, 2001965478L, 1399918158L, 329087331L, 2118472151),
        new SmallPrime(2146060289L, 3416934399L, 279720260L, 1683079797L, 1885781100L, 720846962),
        new SmallPrime(2146048001L, 2645170175L, 1358862595L, 455982466L, 1271329796L, 1757914587),
        new SmallPrime(2146041857L, 2146041855L, 1278996706L, 70834413L, 1910051359L, 2005325121),
        new SmallPrime(2146019329L, 1101637631L, 1427296360L, 223814614L, 804128830L, 1360420507),
        new SmallPrime(2145986561L, 2846435327L, 1292076979L, 1571468773L, 1576376806L, 1922069675),
        new SmallPrime(2145976321L, 2145976319L, 2074996602L, 755073721L, 414238476L, 200778967),
        new SmallPrime(2145964033L, 3605581823L, 689794868L, 1452485437L, 856502235L, 416725139),
        new SmallPrime(2145906689L, 2129129471L, 921247209L, 676947204L, 654865993L, 1812356432),
        new SmallPrime(2145875969L, 3081205759L, 1842525491L, 1326145670L, 309498173L, 2087042424),
        new SmallPrime(2145871873L, 2779211775L, 962993434L, 1019421319L, 1501089909L, 653721111),
        new SmallPrime(2145841153L, 1592193023L, 1869982816L, 214096862L, 1556972726L, 503214451),
        new SmallPrime(2145832961L, 384225279L, 384678957L, 669316727L, 783241415L, 1108169518),
        new SmallPrime(2145816577L, 1860603903L, 1173007304L, 1845637035L, 826920248L, 854681726),
        new SmallPrime(2145785857L, 1571166207L, 669709063L, 1939498551L, 1558887263L, 633181476),
        new SmallPrime(2145755137L, 3689259007L, 1290750531L, 954377268L, 91681545L, 1222833471),
        new SmallPrime(2145742849L, 4008013823L, 223279603L, 943919281L, 916877791L, 1449746232),
        new SmallPrime(2145728513L, 1134901247L, 1222351254L, 2118240515L, 1429489422L, 863342952),
        new SmallPrime(2145699841L, 2745485311L, 1723896025L, 1632678465L, 1030096814L, 1093745233),
        new SmallPrime(2145691649L, 3517229055L, 931453090L, 54415538L, 530452177L, 1965948870),
        new SmallPrime(2145687553L, 1705285631L, 1809739911L, 1581677400L, 1990115608L, 440621047),
        new SmallPrime(2145673217L, 1541693439L, 578267174L, 1396659387L, 1723941977L, 1763556058),
        new SmallPrime(2145630209L, 2879633407L, 328648448L, 40361658L, 1539467270L, 1674729976),
        new SmallPrime(2145595393L, 1457729535L, 255202829L, 1400378563L, 1431207189L, 1703296854),
        new SmallPrime(2145587201L, 518197247L, 2028299732L, 1486218552L, 644522534L, 1270684848),
        new SmallPrime(2145525761L, 4225900543L, 1358930970L, 1145038363L, 878769963L, 1208930401),
        new SmallPrime(2145495041L, 3248596991L, 1641759544L, 431141042L, 1250759371L, 401371792),
        new SmallPrime(2145466369L, 60897279L, 4913761L, 1346487656L, 1003126503L, 1472826275),
        new SmallPrime(2145445889L, 2913003519L, 1950341575L, 209604178L, 2010474299L, 190633762),
        new SmallPrime(2145390593L, 2128613375L, 2111959069L, 1977319511L, 253671817L, 711951677),
        new SmallPrime(2145372161L, 1939851263L, 1176002444L, 230466302L, 1564380182L, 318701388),
        new SmallPrime(2145361921L, 1541382143L, 879247163L, 1458754908L, 789794840L, 1032907577),
        new SmallPrime(2145359873L, 1436522495L, 1349830443L, 160850239L, 842688712L, 977726147),
        new SmallPrime(2145355777L, 1201637375L, 246501130L, 1714048525L, 1460503641L, 167731422),
        new SmallPrime(2145343489L, 295655423L, 2034128553L, 85533909L, 89392725L, 1282860004),
        new SmallPrime(2145325057L, 2665418751L, 1621650965L, 1422485923L, 288923379L, 911614733),
        new SmallPrime(2145318913L, 1872689151L, 659138019L, 507124624L, 1410089651L, 131619480),
        new SmallPrime(2145312769L, 1004462079L, 2144542130L, 601028071L, 309708965L, 1694238751),
        new SmallPrime(2145300481L, 3336482815L, 1732525390L, 1940872306L, 2008521418L, 313780675),
        new SmallPrime(2145282049L, 4120799231L, 1238750391L, 1797682851L, 557030104L, 1782843353),
        new SmallPrime(2145232897L, 27109375L, 1795543839L, 1037744642L, 1557646622L, 900304126),
        new SmallPrime(2145218561L, 4007489535L, 1315715750L, 1070968735L, 736559602L, 1443800839),
        new SmallPrime(2145187841L, 1335687167L, 13023648L, 1309044639L, 1866831222L, 399747323),
        new SmallPrime(2145181697L, 3151814655L, 1518459244L, 1993483355L, 691838010L, 615546001),
        new SmallPrime(2145175553L, 597477375L, 1181348151L, 703506449L, 1514472518L, 488565539),
        new SmallPrime(2145079297L, 2262519807L, 1805182441L, 1948507193L, 198785189L, 638515383),
        new SmallPrime(2145021953L, 2396680191L, 1021023200L, 896525833L, 286005075L, 1441021596),
        new SmallPrime(2145015809L, 2174375935L, 2117792680L, 1601002661L, 347118679L, 1157904101),
        new SmallPrime(2145003521L, 1503275007L, 929358646L, 1668692954L, 37145528L, 641759047),
        new SmallPrime(2144960513L, 1071218687L, 2014663077L, 813829821L, 819596313L, 63874550),
        new SmallPrime(2144944129L, 4023992319L, 1528910090L, 1405891381L, 2038603849L, 1893894911),
        new SmallPrime(2144935937L, 1004085247L, 1020776636L, 1264145057L, 131908564L, 1538042171),
        new SmallPrime(2144894977L, 1071153151L, 117657395L, 1658068046L, 1863791040L, 315971662),
        new SmallPrime(2144888833L, 3583535103L, 1035349752L, 631245709L, 221833619L, 2112295623),
        new SmallPrime(2144880641L, 3952625663L, 2014869161L, 1728523150L, 1503761981L, 2007376047),
        new SmallPrime(2144864257L, 4288153599L, 1298675209L, 544582669L, 1272986021L, 1915075671),
        new SmallPrime(2144827393L, 3080157183L, 1123663006L, 160045786L, 473495242L, 813254440),
        new SmallPrime(2144806913L, 1234642943L, 731069394L, 2014053847L, 1766552348L, 857353667),
        new SmallPrime(2144796673L, 2144796671L, 1796197228L, 54648661L, 347997993L, 2024676365),
        new SmallPrime(2144778241L, 1536604159L, 685152946L, 1271425633L, 321887990L, 1348054165),
        new SmallPrime(2144759809L, 248934399L, 154114551L, 1526637592L, 1903326258L, 1432928978),
        new SmallPrime(2144757761L, 1972791295L, 1931452899L, 132970580L, 1653105390L, 1502698593),
        new SmallPrime(2144729089L, 3751147519L, 319799485L, 856415206L, 2139235838L, 1528560673),
        new SmallPrime(2144727041L, 1054207999L, 456984744L, 648538322L, 1211291222L, 1040592824),
        new SmallPrime(2144696321L, 2543155199L, 117528426L, 396088388L, 1351318511L, 1705928562),
        new SmallPrime(2144667649L, 798296063L, 1911610943L, 338234827L, 1761076141L, 1009985854),
        new SmallPrime(2144573441L, 2509477887L, 1537057360L, 570563298L, 647497611L, 1089716523),
        new SmallPrime(2144555009L, 2261995519L, 1259090311L, 1645617705L, 28737134L, 311642518),
        new SmallPrime(2144550913L, 4023599103L, 852574554L, 787018986L, 2112202798L, 1495324375),
        new SmallPrime(2144536577L, 1335035903L, 1561905341L, 2144070617L, 1829236274L, 1925884510),
        new SmallPrime(2144524289L, 1771231231L, 418502709L, 1096979205L, 1715236225L, 4942642),
        new SmallPrime(2144512001L, 1905436671L, 486401965L, 371227774L, 621999585L, 248381878),
        new SmallPrime(2144468993L, 2144468991L, 1685175757L, 1750820831L, 1117048683L, 375676619),
        new SmallPrime(2144458753L, 428988415L, 481506649L, 208627545L, 1669824636L, 994022373),
        new SmallPrime(2144421889L, 3683731455L, 111545270L, 490157982L, 598291463L, 1506819637),
        new SmallPrime(2144409601L, 1301354495L, 1696205610L, 1773221211L, 1190821653L, 703296234),
        new SmallPrime(2144370689L, 1070628863L, 410368360L, 1852155417L, 1241424587L, 1916372058),
        new SmallPrime(2144348161L, 2039490559L, 2058332258L, 380925864L, 1763931741L, 1080707148),
        new SmallPrime(2144335873L, 2140141567L, 189861841L, 1779816687L, 1155526432L, 1839197694),
        new SmallPrime(2144329729L, 2077220863L, 782112649L, 526684710L, 1471877186L, 752775607),
        new SmallPrime(2144327681L, 2039470079L, 1046835057L, 1247119374L, 1480573420L, 1180445680),
        new SmallPrime(2144309249L, 1322225663L, 655091351L, 1601388861L, 546007835L, 774279999),
        new SmallPrime(2144296961L, 466575359L, 478811653L, 1358954307L, 277575578L, 1153736907),
        new SmallPrime(2144290817L, 4220471295L, 845002172L, 1883197481L, 20750133L, 894854412),
        new SmallPrime(2144243713L, 1859031039L, 134828934L, 37546773L, 2020974628L, 1786100430),
        new SmallPrime(2144237569L, 663648255L, 981840700L, 819447626L, 828884995L, 1801721724),
        new SmallPrime(2144161793L, 1187860479L, 1312055195L, 492584100L, 284208050L, 991575551),
        new SmallPrime(2144155649L, 3280812031L, 1909512015L, 1038118142L, 557855457L, 351830619),
        new SmallPrime(2144137217L, 516747263L, 1231150697L, 1145329724L, 1222551687L, 1031682140),
        new SmallPrime(2144120833L, 2261561343L, 1725803932L, 85962159L, 1987475484L, 1541754604),
        new SmallPrime(2144071681L, 4274778111L, 1126468398L, 956241311L, 2057260329L, 1928592240),
        new SmallPrime(2144065537L, 965466111L, 1878788832L, 1053853819L, 587985232L, 409128944),
        new SmallPrime(2144028673L, 998983679L, 2034483463L, 85221004L, 692082439L, 2038010847),
        new SmallPrime(2144010241L, 2144010239L, 842246168L, 579253174L, 1365083779L, 893394979),
        new SmallPrime(2143952897L, 4224327679L, 754017578L, 344456872L, 574149142L, 951344542),
        new SmallPrime(2143940609L, 3200905215L, 1412359304L, 691093159L, 1240091757L, 209561083),
        new SmallPrime(2143918081L, 1971951615L, 1479162717L, 1879554798L, 1607927459L, 1573571606),
        new SmallPrime(2143899649L, 3335081983L, 1640239719L, 1873544500L, 512920036L, 980684517),
        new SmallPrime(2143891457L, 382283775L, 1157894649L, 1152913815L, 2082762113L, 174332902),
        new SmallPrime(2143885313L, 3448313855L, 77697446L, 1272010615L, 639087423L, 788921708),
        new SmallPrime(2143854593L, 466132991L, 1365712904L, 234096830L, 1624926716L, 267754392),
        new SmallPrime(2143836161L, 3783809023L, 200636173L, 802456521L, 972626486L, 1579347589),
        new SmallPrime(2143762433L, 1669806079L, 1374323479L, 122119316L, 523389252L, 528105857),
        new SmallPrime(2143756289L, 3150389247L, 226179777L, 2014672070L, 977903451L, 205998472),
        new SmallPrime(2143713281L, 2810607615L, 674480231L, 499335673L, 1923111979L, 170087335),
        new SmallPrime(2143690753L, 4224065535L, 1119041321L, 293184621L, 548955895L, 1709245285),
        new SmallPrime(2143670273L, 3066417151L, 966374918L, 1855603277L, 1939385627L, 1369708529),
        new SmallPrime(2143666177L, 1875230719L, 1768658380L, 1592043954L, 626312185L, 1973259130),
        new SmallPrime(2143645697L, 4005916671L, 1369531559L, 1496045094L, 868812782L, 743685730),
        new SmallPrime(2143641601L, 2613403647L, 836342892L, 137205166L, 1106494107L, 1904307099),
        new SmallPrime(2143635457L, 2609203199L, 1360942100L, 611274460L, 1430230424L, 634419595),
        new SmallPrime(2143621121L, 3737456639L, 190915397L, 456284400L, 1315071924L, 1321642891),
        new SmallPrime(2143598593L, 998553599L, 155040255L, 1748188596L, 1338143543L, 1507236419),
        new SmallPrime(2143567873L, 1875132415L, 1411420224L, 961363405L, 100809682L, 4006593),
        new SmallPrime(2143561729L, 964962303L, 1285824486L, 1745824462L, 52145752L, 1898334618),
        new SmallPrime(2143553537L, 1065617407L, 1589845870L, 642609604L, 876094006L, 1160521531),
        new SmallPrime(2143541249L, 3112425471L, 912656057L, 832462820L, 1196335625L, 50513987),
        new SmallPrime(2143531009L, 3871584255L, 917220898L, 733999433L, 1060086084L, 957280110),
        new SmallPrime(2143522817L, 3468922879L, 1098374569L, 2134640637L, 1812047677L, 957540519),
        new SmallPrime(2143506433L, 2260946943L, 933735606L, 453429574L, 1461625836L, 1505133897),
        new SmallPrime(2143488001L, 1333987327L, 1985268644L, 2074073649L, 1479394204L, 1064535019),
        new SmallPrime(2143469569L, 4022517759L, 1477872272L, 1513292205L, 1862896587L, 1004190392),
        new SmallPrime(2143426561L, 495065087L, 186028039L, 1375811092L, 173677829L, 2089274222),
        new SmallPrime(2143383553L, 1858170879L, 887211384L, 2134877926L, 1608809321L, 1661543915),
        new SmallPrime(2143377409L, 2978043903L, 1281348890L, 1288031453L, 483115862L, 1970346133),
        new SmallPrime(2143363073L, 1002512383L, 1950991422L, 1122143235L, 1166724603L, 1788468869),
        new SmallPrime(2143260673L, 1321177087L, 2069683714L, 1805244102L, 1261092720L, 903254409),
        new SmallPrime(2143246337L, 293558271L, 1176536351L, 1950051640L, 309589975L, 1529083832),
        new SmallPrime(2143209473L, 58640383L, 1560740565L, 218139536L, 1255708568L, 1698136215),
        new SmallPrime(2143203329L, 3334385663L, 1971538547L, 1999593400L, 1859827267L, 1956345089),
        new SmallPrime(2143160321L, 2675836927L, 477597631L, 884313283L, 2113130436L, 1457071330),
        new SmallPrime(2143129601L, 2394787839L, 546440653L, 1295503643L, 2080414454L, 1206105984),
        new SmallPrime(2143123457L, 394098687L, 612601193L, 1693130333L, 43008062L, 16697358),
        new SmallPrime(2143100929L, 1002250239L, 591926265L, 1101885974L, 315537419L, 1388418046),
        new SmallPrime(2143092737L, 2143092735L, 1205498739L, 1071117401L, 276869279L, 1816776329),
        new SmallPrime(2143082497L, 2306660351L, 1658985163L, 2126748415L, 1783080325L, 1909268843),
        new SmallPrime(2143062017L, 2004649983L, 806738297L, 1026831688L, 2035800522L, 754902198),
        new SmallPrime(2143051777L, 1539071999L, 1644097744L, 1553535853L, 2052185358L, 1958510946),
        new SmallPrime(2143025153L, 1065089023L, 48408341L, 1525464941L, 1852830390L, 755368639),
        new SmallPrime(2143006721L, 3871059967L, 1620020706L, 2138083553L, 30951415L, 1336294009),
        new SmallPrime(2142996481L, 2273019903L, 576994614L, 1383533L, 1219848252L, 1992320100),
        new SmallPrime(2142976001L, 2742761471L, 1018092510L, 1458681506L, 608641646L, 190719147),
        new SmallPrime(2142965761L, 515575807L, 359250737L, 1354388786L, 2131235811L, 1323949132),
        new SmallPrime(2142916609L, 649744383L, 1208571887L, 1878268755L, 1560639011L, 1431405645),
        new SmallPrime(2142892033L, 3199856639L, 340166218L, 771047236L, 1294387139L, 1350670764),
        new SmallPrime(2142885889L, 2574899199L, 1416994273L, 869163275L, 1050539694L, 503308105),
        new SmallPrime(2142871553L, 3686375423L, 823272682L, 1330177723L, 1618603945L, 680882765),
        new SmallPrime(2142861313L, 3615062015L, 185662521L, 59302870L, 1086768468L, 733960061),
        new SmallPrime(2142830593L, 2142830591L, 1184916005L, 1079082138L, 801265712L, 1392890675),
        new SmallPrime(2142803969L, 2776143871L, 1898723413L, 612877597L, 748398044L, 719565838),
        new SmallPrime(2142785537L, 1723355135L, 1938177810L, 901717203L, 1921085502L, 1027865724),
        new SmallPrime(2142779393L, 4084742143L, 1843695270L, 956271994L, 1027905184L, 437886242),
        new SmallPrime(2142724097L, 465002495L, 1786272468L, 2057509300L, 596569659L, 833374098),
        new SmallPrime(2142707713L, 3149340671L, 1726077360L, 110521970L, 490838647L, 570919333),
        new SmallPrime(2142658561L, 3686162431L, 1631852094L, 261726828L, 931917077L, 2005364883),
        new SmallPrime(2142638081L, 3199602687L, 1965839564L, 408833423L, 1753224304L, 1152041797),
        new SmallPrime(2142564353L, 515174399L, 1080848266L, 1517077457L, 1548531958L, 823650319),
        new SmallPrime(2142533633L, 2272557055L, 392959315L, 729816849L, 921579743L, 129344001),
        new SmallPrime(2142529537L, 359950335L, 17286407L, 1033141274L, 1935370361L, 921394177),
        new SmallPrime(2142527489L, 1538547711L, 2022542562L, 591503450L, 522147359L, 1308604835),
        new SmallPrime(2142502913L, 2142502911L, 862145305L, 1943481144L, 733213664L, 364106514),
        new SmallPrime(2142498817L, 4273205247L, 1497774797L, 362230737L, 2058280749L, 1402886894),
        new SmallPrime(2142416897L, 1186115583L, 1820910794L, 869347984L, 1125890573L, 823863638),
        new SmallPrime(2142363649L, 4222738431L, 1712603348L, 839133619L, 770847044L, 892019817),
        new SmallPrime(2142351361L, 2796662783L, 639467496L, 472998882L, 259948593L, 1009074045),
        new SmallPrime(2142330881L, 2612092927L, 119930462L, 898272391L, 1611927337L, 1395844887),
        new SmallPrime(2142314497L, 1001463807L, 417987874L, 952188691L, 1472816610L, 1525176551),
        new SmallPrime(2142289921L, 4021338111L, 626441030L, 1087932551L, 1833827076L, 1287261065),
        new SmallPrime(2142283777L, 292595711L, 1437184719L, 1990366532L, 718034518L, 2141169066),
        new SmallPrime(2142277633L, 783323135L, 409102935L, 1109091844L, 1748524698L, 1404314290),
        new SmallPrime(2142263297L, 1634752511L, 1332648256L, 1458675808L, 1235027887L, 862728790),
        new SmallPrime(2142208001L, 1068466175L, 1397470467L, 1098539573L, 1097351447L, 1965727963),
        new SmallPrime(2142164993L, 695130111L, 1066940847L, 3394841L, 337813155L, 2063673696),
        new SmallPrime(2142097409L, 3064844287L, 297225318L, 190118328L, 120193342L, 1425544062),
        new SmallPrime(2142087169L, 23963647L, 1242139544L, 508378060L, 2119898775L, 1644552134),
        new SmallPrime(2142078977L, 1735231487L, 34593522L, 662801663L, 1837465297L, 1388457155),
        new SmallPrime(2142074881L, 393050111L, 1775241888L, 546076620L, 753976398L, 1286409850),
        new SmallPrime(2142044161L, 2142044159L, 1993353265L, 20421055L, 1749390716L, 453736187),
        new SmallPrime(2142025729L, 3144464383L, 1054015161L, 216592874L, 174350842L, 2000851460),
        new SmallPrime(2142011393L, 1068269567L, 1022043540L, 154322400L, 1364900586L, 1216790673),
        new SmallPrime(2141974529L, 4272680959L, 878244511L, 450223676L, 840019780L, 936239114),
        new SmallPrime(2141943809L, 4151015423L, 181002276L, 1766661673L, 824377997L, 898611110),
        new SmallPrime(2141933569L, 2259374079L, 921001808L, 138101311L, 1597307312L, 1683735850),
        new SmallPrime(2141931521L, 996886527L, 2026958630L, 786548950L, 491009818L, 1487647832),
        new SmallPrime(2141902849L, 3916093439L, 1774507217L, 1503936548L, 2011562721L, 803375820),
        new SmallPrime(2141890561L, 4050298879L, 1242184656L, 1606622264L, 973689178L, 1298783383),
        new SmallPrime(2141857793L, 1500129279L, 47495456L, 1818547837L, 1933124842L, 1327760978),
        new SmallPrime(2141833217L, 3546925055L, 534546202L, 500926182L, 2118308644L, 1496583315),
        new SmallPrime(2141820929L, 1969854463L, 457737750L, 690473287L, 1631943651L, 664359263),
        new SmallPrime(2141786113L, 1588137983L, 1837204275L, 2031397019L, 1967667817L, 905298889),
        new SmallPrime(2141771777L, 1231607807L, 1585138177L, 54300848L, 314627109L, 927112984),
        new SmallPrime(2141759489L, 2439555071L, 1154859258L, 1383763879L, 1971680935L, 2108097075),
        new SmallPrime(2141749249L, 1068007423L, 2080935967L, 859296785L, 1197038946L, 1880678306),
        new SmallPrime(2141685761L, 3345451007L, 517526213L, 1550564945L, 1632439809L, 954282998),
        new SmallPrime(2141673473L, 2439469055L, 21649850L, 2046161117L, 552392242L, 103443123),
        new SmallPrime(2141669377L, 2070366207L, 126251361L, 1742664644L, 1482481206L, 605966977),
        new SmallPrime(2141655041L, 2661748735L, 484235305L, 1528710429L, 354771509L, 593774164),
        new SmallPrime(2141587457L, 526780415L, 574888543L, 955760611L, 687141330L, 1030240293),
        new SmallPrime(2141583361L, 3748001791L, 1371847173L, 1144908715L, 184363050L, 1463331359),
        new SmallPrime(2141575169L, 1499846655L, 1229054288L, 1538452901L, 324225902L, 1425521960),
        new SmallPrime(2141546497L, 1164273663L, 1768171221L, 430711736L, 1956344257L, 358993721),
        new SmallPrime(2141515777L, 514125823L, 1883057193L, 1188334395L, 254771003L, 712002260),
        new SmallPrime(2141495297L, 463773695L, 1894308447L, 2130610251L, 69814245L, 2121489732),
        new SmallPrime(2141483009L, 3466883071L, 950896971L, 1258381554L, 825951706L, 195111043),
        new SmallPrime(2141458433L, 4272164863L, 566843170L, 64378978L, 1469333580L, 1466451638),
        new SmallPrime(2141360129L, 4003631103L, 510897768L, 2054940459L, 1603364306L, 1896684552),
        new SmallPrime(2141325313L, 1600260095L, 309414728L, 783698702L, 845623609L, 2118606683),
        new SmallPrime(2141317121L, 3714181119L, 41910923L, 1219136678L, 1953984968L, 1546372147),
        new SmallPrime(2141286401L, 1856073727L, 1171249093L, 986447182L, 1530161498L, 1359644991),
        new SmallPrime(2141267969L, 694233087L, 1211133465L, 844862374L, 1144145329L, 372321440),
        new SmallPrime(2141255681L, 2405496831L, 1328934139L, 735916721L, 1296250995L, 1612100237),
        new SmallPrime(2141243393L, 3814770687L, 520479708L, 525351865L, 2026115355L, 538775863),
        new SmallPrime(2141214721L, 1633703935L, 504183101L, 663335303L, 715380266L, 34424233),
        new SmallPrime(2141212673L, 1721782271L, 603205901L, 1433302473L, 1227056696L, 1431443425),
        new SmallPrime(2141202433L, 2036344831L, 1604585501L, 1229299775L, 752118923L, 1881364338),
        new SmallPrime(2141175809L, 1872740351L, 1733516714L, 1807146636L, 363795364L, 1986179172),
        new SmallPrime(2141165569L, 1432328191L, 1490189496L, 1320489732L, 1641816108L, 578038958),
        new SmallPrime(2141073409L, 916336639L, 873759781L, 1556749288L, 651084845L, 1015729028),
        new SmallPrime(2141052929L, 4221427711L, 979186233L, 1762104735L, 715530680L, 622504930),
        new SmallPrime(2141040641L, 647868415L, 949895441L, 2025808022L, 160106133L, 561883724),
        new SmallPrime(2141028353L, 1067286527L, 2135838697L, 916257035L, 1402241650L, 750070220),
        new SmallPrime(2141011969L, 4020060159L, 611765852L, 1170281544L, 1552266705L, 329393240),
        new SmallPrime(2140999681L, 3734835199L, 351284530L, 998193986L, 2031434746L, 1932350134),
        new SmallPrime(2140997633L, 2942109695L, 1853355265L, 1185563841L, 237281087L, 440781243),
        new SmallPrime(2140993537L, 1331492863L, 676775069L, 495308340L, 1540448795L, 750732506),
        new SmallPrime(2140942337L, 4137431039L, 1645148093L, 2102068815L, 2126166725L, 1570735656),
        new SmallPrime(2140925953L, 4271632383L, 759158315L, 871789395L, 712612398L, 1673104634),
        new SmallPrime(2140917761L, 4137406463L, 1126409570L, 487401697L, 790501062L, 2003679138),
        new SmallPrime(2140887041L, 3512424447L, 356555373L, 692136311L, 1873117117L, 2056578551),
        new SmallPrime(2140837889L, 1163565055L, 367306155L, 549695576L, 1539499132L, 1848189620),
        new SmallPrime(2140788737L, 2572802047L, 558045408L, 1566269988L, 1578573646L, 1949949243),
        new SmallPrime(2140766209L, 1067024383L, 1525495467L, 488941581L, 1332090451L, 1459824799),
        new SmallPrime(2140764161L, 3612964863L, 453729911L, 1063914887L, 792300614L, 2072819909),
        new SmallPrime(2140696577L, 1318612991L, 1154564043L, 1249648119L, 1281812338L, 1631278485),
        new SmallPrime(2140684289L, 1872248831L, 2118938259L, 262551938L, 1786968202L, 55549365),
        new SmallPrime(2140653569L, 4082616319L, 1285628803L, 80889460L, 220865202L, 1452802925),
        new SmallPrime(2140594177L, 999743487L, 1382704526L, 1592957652L, 139669975L, 754621876),
        new SmallPrime(2140579841L, 827762687L, 1601414172L, 1043280L, 1283241479L, 868432049),
        new SmallPrime(2140569601L, 1066827775L, 1241785107L, 1615837280L, 1527749306L, 716582032),
        new SmallPrime(2140567553L, 2807461887L, 1271173854L, 308478868L, 1555967545L, 1677706733),
        new SmallPrime(2140557313L, 2794868735L, 1924690389L, 870598748L, 61433129L, 210103227),
        new SmallPrime(2140549121L, 915812351L, 914841856L, 1101232091L, 1853042589L, 1050014014),
        new SmallPrime(2140477441L, 2371164159L, 1750444466L, 1671883855L, 1913885020L, 2062319628),
        new SmallPrime(2140469249L, 3478452223L, 1728225499L, 134343014L, 1041454719L, 2059902640),
        new SmallPrime(2140426241L, 3868479487L, 1380594798L, 1812066246L, 2026797853L, 1212120681),
        new SmallPrime(2140420097L, 3008641023L, 100654027L, 100993678L, 1430665666L, 1692048952),
        new SmallPrime(2140413953L, 2073305087L, 1265120041L, 624789355L, 1980101954L, 1492294345),
        new SmallPrime(2140383233L, 559130623L, 945242106L, 782288940L, 181435266L, 453459249),
        new SmallPrime(2140366849L, 3847448575L, 1028300358L, 75248903L, 312878281L, 294392408),
        new SmallPrime(2140354561L, 592656383L, 368901374L, 2124984708L, 1584238035L, 557749653),
        new SmallPrime(2140348417L, 3146981375L, 495209562L, 1555795152L, 1821951283L, 1580622851),
        new SmallPrime(0L, 0L, 0L, 0L, 0L, 0)
    };

    static final short[] REV10 = new short[]{
        0, 512, 256, 768, 128, 640, 384, 896, 64, 576, 320, 832,
        192, 704, 448, 960, 32, 544, 288, 800, 160, 672, 416, 928,
        96, 608, 352, 864, 224, 736, 480, 992, 16, 528, 272, 784,
        144, 656, 400, 912, 80, 592, 336, 848, 208, 720, 464, 976,
        48, 560, 304, 816, 176, 688, 432, 944, 112, 624, 368, 880,
        240, 752, 496, 1008, 8, 520, 264, 776, 136, 648, 392, 904,
        72, 584, 328, 840, 200, 712, 456, 968, 40, 552, 296, 808,
        168, 680, 424, 936, 104, 616, 360, 872, 232, 744, 488, 1000,
        24, 536, 280, 792, 152, 664, 408, 920, 88, 600, 344, 856,
        216, 728, 472, 984, 56, 568, 312, 824, 184, 696, 440, 952,
        120, 632, 376, 888, 248, 760, 504, 1016, 4, 516, 260, 772,
        132, 644, 388, 900, 68, 580, 324, 836, 196, 708, 452, 964,
        36, 548, 292, 804, 164, 676, 420, 932, 100, 612, 356, 868,
        228, 740, 484, 996, 20, 532, 276, 788, 148, 660, 404, 916,
        84, 596, 340, 852, 212, 724, 468, 980, 52, 564, 308, 820,
        180, 692, 436, 948, 116, 628, 372, 884, 244, 756, 500, 1012,
        12, 524, 268, 780, 140, 652, 396, 908, 76, 588, 332, 844,
        204, 716, 460, 972, 44, 556, 300, 812, 172, 684, 428, 940,
        108, 620, 364, 876, 236, 748, 492, 1004, 28, 540, 284, 796,
        156, 668, 412, 924, 92, 604, 348, 860, 220, 732, 476, 988,
        60, 572, 316, 828, 188, 700, 444, 956, 124, 636, 380, 892,
        252, 764, 508, 1020, 2, 514, 258, 770, 130, 642, 386, 898,
        66, 578, 322, 834, 194, 706, 450, 962, 34, 546, 290, 802,
        162, 674, 418, 930, 98, 610, 354, 866, 226, 738, 482, 994,
        18, 530, 274, 786, 146, 658, 402, 914, 82, 594, 338, 850,
        210, 722, 466, 978, 50, 562, 306, 818, 178, 690, 434, 946,
        114, 626, 370, 882, 242, 754, 498, 1010, 10, 522, 266, 778,
        138, 650, 394, 906, 74, 586, 330, 842, 202, 714, 458, 970,
        42, 554, 298, 810, 170, 682, 426, 938, 106, 618, 362, 874,
        234, 746, 490, 1002, 26, 538, 282, 794, 154, 666, 410, 922,
        90, 602, 346, 858, 218, 730, 474, 986, 58, 570, 314, 826,
        186, 698, 442, 954, 122, 634, 378, 890, 250, 762, 506, 1018,
        6, 518, 262, 774, 134, 646, 390, 902, 70, 582, 326, 838,
        198, 710, 454, 966, 38, 550, 294, 806, 166, 678, 422, 934,
        102, 614, 358, 870, 230, 742, 486, 998, 22, 534, 278, 790,
        150, 662, 406, 918, 86, 598, 342, 854, 214, 726, 470, 982,
        54, 566, 310, 822, 182, 694, 438, 950, 118, 630, 374, 886,
        246, 758, 502, 1014, 14, 526, 270, 782, 142, 654, 398, 910,
        78, 590, 334, 846, 206, 718, 462, 974, 46, 558, 302, 814,
        174, 686, 430, 942, 110, 622, 366, 878, 238, 750, 494, 1006,
        30, 542, 286, 798, 158, 670, 414, 926, 94, 606, 350, 862,
        222, 734, 478, 990, 62, 574, 318, 830, 190, 702, 446, 958,
        126, 638, 382, 894, 254, 766, 510, 1022, 1, 513, 257, 769,
        129, 641, 385, 897, 65, 577, 321, 833, 193, 705, 449, 961,
        33, 545, 289, 801, 161, 673, 417, 929, 97, 609, 353, 865,
        225, 737, 481, 993, 17, 529, 273, 785, 145, 657, 401, 913,
        81, 593, 337, 849, 209, 721, 465, 977, 49, 561, 305, 817,
        177, 689, 433, 945, 113, 625, 369, 881, 241, 753, 497, 1009,
        9, 521, 265, 777, 137, 649, 393, 905, 73, 585, 329, 841,
        201, 713, 457, 969, 41, 553, 297, 809, 169, 681, 425, 937,
        105, 617, 361, 873, 233, 745, 489, 1001, 25, 537, 281, 793,
        153, 665, 409, 921, 89, 601, 345, 857, 217, 729, 473, 985,
        57, 569, 313, 825, 185, 697, 441, 953, 121, 633, 377, 889,
        249, 761, 505, 1017, 5, 517, 261, 773, 133, 645, 389, 901,
        69, 581, 325, 837, 197, 709, 453, 965, 37, 549, 293, 805,
        165, 677, 421, 933, 101, 613, 357, 869, 229, 741, 485, 997,
        21, 533, 277, 789, 149, 661, 405, 917, 85, 597, 341, 853,
        213, 725, 469, 981, 53, 565, 309, 821, 181, 693, 437, 949,
        117, 629, 373, 885, 245, 757, 501, 1013, 13, 525, 269, 781,
        141, 653, 397, 909, 77, 589, 333, 845, 205, 717, 461, 973,
        45, 557, 301, 813, 173, 685, 429, 941, 109, 621, 365, 877,
        237, 749, 493, 1005, 29, 541, 285, 797, 157, 669, 413, 925,
        93, 605, 349, 861, 221, 733, 477, 989, 61, 573, 317, 829,
        189, 701, 445, 957, 125, 637, 381, 893, 253, 765, 509, 1021,
        3, 515, 259, 771, 131, 643, 387, 899, 67, 579, 323, 835,
        195, 707, 451, 963, 35, 547, 291, 803, 163, 675, 419, 931,
        99, 611, 355, 867, 227, 739, 483, 995, 19, 531, 275, 787,
        147, 659, 403, 915, 83, 595, 339, 851, 211, 723, 467, 979,
        51, 563, 307, 819, 179, 691, 435, 947, 115, 627, 371, 883,
        243, 755, 499, 1011, 11, 523, 267, 779, 139, 651, 395, 907,
        75, 587, 331, 843, 203, 715, 459, 971, 43, 555, 299, 811,
        171, 683, 427, 939, 107, 619, 363, 875, 235, 747, 491, 1003,
        27, 539, 283, 795, 155, 667, 411, 923, 91, 603, 347, 859,
        219, 731, 475, 987, 59, 571, 315, 827, 187, 699, 443, 955,
        123, 635, 379, 891, 251, 763, 507, 1019, 7, 519, 263, 775,
        135, 647, 391, 903, 71, 583, 327, 839, 199, 711, 455, 967,
        39, 551, 295, 807, 167, 679, 423, 935, 103, 615, 359, 871,
        231, 743, 487, 999, 23, 535, 279, 791, 151, 663, 407, 919,
        87, 599, 343, 855, 215, 727, 471, 983, 55, 567, 311, 823,
        183, 695, 439, 951, 119, 631, 375, 887, 247, 759, 503, 1015,
        15, 527, 271, 783, 143, 655, 399, 911, 79, 591, 335, 847,
        207, 719, 463, 975, 47, 559, 303, 815, 175, 687, 431, 943,
        111, 623, 367, 879, 239, 751, 495, 1007, 31, 543, 287, 799,
        159, 671, 415, 927, 95, 607, 351, 863, 223, 735, 479, 991,
        63, 575, 319, 831, 191, 703, 447, 959, 127, 639, 383, 895,
        255, 767, 511, 1023
    };


    // Calculate parity of coefficients
    public static int parity(int logn, byte[] f)
    {
        int n = 1 << logn;
        int pp = 0;
        for (int u = 0; u < n; u++)
        {
            pp += f[u];
        }
        return pp & 1;
    }

    // Convert integer to fixed-point
    public static long fxrOf(int j)
    {
        return ((long)j << 32);
    }

    // Fixed-point reciprocal (inverse)
    public static long fxrInv(long y)
    {
        // Check for division by zero
        if (y == 0)
        {
            throw new ArithmeticException("Division by zero in fxrInv");
        }
        return innerFxrDiv(1L << 32, y);
    }

    public static long innerFxrDiv(long x, long y)
    {
        // Handle division by zero
        if (y == 0)
        {
            throw new ArithmeticException("Division by zero");
        }

        /*
         * Get absolute values and signs. From now on, we can suppose
         * that x and y fit on 63 bits (we ignore edge conditions).
         */
        long sx = x >>> 63;  // Use unsigned shift for sign extraction
        x = (x ^ -sx) + sx;

        long sy = y >>> 63;  // Use unsigned shift for sign extraction
        y = (y ^ -sy) + sy;

        /*
         * Do a bit by bit division, assuming that the quotient fits.
         * The numerator starts at x*2^31, and is shifted one bit a time.
         */
        long q = 0;
        long num = x >>> 31;  // Use unsigned shift

        // First loop: bits 63 to 33
        for (int i = 63; i >= 33; i--)
        {
            long diff = num - y;
            long b = 1 - (diff >>> 63);  // Use unsigned shift for sign bit
            q |= b << i;
            num -= y & -b;
            num <<= 1;
            num |= (x >>> (i - 33)) & 1;
        }

        // Second loop: bits 32 to 0
        for (int i = 32; i >= 0; i--)
        {
            long diff = num - y;
            long b = 1 - (diff >>> 63);  // Use unsigned shift for sign bit
            q |= b << i;
            num -= y & -b;
            num <<= 1;
        }

        /*
         * Rounding: if the remainder is at least y/2 (scaled), we add
         * 2^(-32) to the quotient.
         */
        long diff = num - y;
        long b = 1 - (diff >>> 63);  // Use unsigned shift for sign bit
        q += b;

        /*
         * Sign management: if the original x and y had different signs,
         * then we must negate the quotient.
         */
        sx ^= sy;
        q = (q ^ -sx) + sx;

        return q;
    }

    // Convert f and g to RNS+NTT form using precomputed primes
    public static void makeFgZero(int logn, byte[] f, int fOff, byte[] g, int gOff, int[] tmp, int tmpOff)
    {
        int n = 1 << logn;
        SmallPrime prime = PRIMES[0];  // Use first prime

        // Split tmp buffer into three segments
        int ft = tmpOff;
        int gt = ft + n;
        int gm = gt + n;

        polyMpSetSmall(logn, tmp, ft, f, fOff, (int)prime.p);
        polyMpSetSmall(logn, tmp, gt, g, gOff, (int)prime.p);
        mpMkgm(logn, tmp, gm, (int)prime.g, (int)prime.p, (int)prime.p0i);
        mpNTT(logn, tmp, ft, tmp, gm, (int)prime.p, (int)prime.p0i);
        mpNTT(logn, tmp, gt, tmp, gm, (int)prime.p, (int)prime.p0i);
    }

    // Convert polynomial coefficients to residues mod p
    private static void polyMpSetSmall(int logn, int[] d, int dOff, byte[] f, int fOff, int p)
    {
        int n = 1 << logn;
        for (int u = 0; u < n; u++)
        {
            d[u + dOff] = mp_set(f[fOff + u], p);
        }
    }

    // Precompute NTT roots (gm array)
    public static void mpMkgm(int logn, int[] gm, int gmOff, int g, int p, int p0i)
    {
        int n = 1 << logn;

        // Precompute g^(2^(10-logn))
        for (int j = logn; j < 10; j++)
        {
            g = mpMontyMul(g, g, p, p0i);
        }

        int k = 10 - logn;
        int x1 = mpR(p);
        int u = 0;
        while (true)
        {
            int v = REV10[u << k];
            gm[v + gmOff] = x1;
            u++;
            if (u >= n)
            {
                break;
            }
            x1 = mpMontyMul(x1, g, p, p0i);
        }
    }

    // Compute R = 2^32 mod p
    private static int mpR(int p)
    {
        return -(p << 1);  // Equivalent to 2^32 - 2p
    }

    // Number Theoretic Transform (in-place)
    public static void mpNTT(int logn, int[] a, int aOff, int[] gm, int gmOff, int p, int p0i)
    {
        if (logn == 0)
        {
            return;
        }

        int t = 1 << logn;
        for (int lm = 0; lm < logn; lm++)
        {
            int m = 1 << lm;
            int ht = t >>> 1;
            int v0 = 0;
            for (int u = 0; u < m; u++)
            {
                int s = gm[gmOff + u + m];
                for (int v = 0; v < ht; v++)
                {
                    int k1 = v0 + v;
                    int k2 = k1 + ht;
                    int x1 = a[aOff + k1];
                    int x2 = mpMontyMul(a[aOff + k2], s, p, p0i);
                    a[aOff + k1] = mpAdd(x1, x2, p);
                    a[aOff + k2] = mpSub(x1, x2, p);
                }
                v0 += t;
            }
            t = ht;
        }
    }

    // Modular arithmetic functions
    static int tbmask(int x)
    {
        return x >> 31;
    }

    public static int mp_set(int v, int p)
    {
        int w = v;
        return w + (p & tbmask(w));
    }

    public static int mpAdd(int a, int b, int p)
    {
        int d = a + b - p;
        return d + (p & tbmask(d));
    }

    private static int mpSub(int a, int b, int p)
    {
        int d = a - b;
        return d + (p & tbmask(d));
    }

    public static int mpMontyMul(int a, int b, int p, int p0i)
    {
        long z = (long)a * (long)b;
        long w = (z * p0i) & 0xFFFFFFFFL;
        long temp = z + w * (long)p;
        int d = (int)(temp >>> 32) - p;
        return d + (p & tbmask(d));
    }

    // Additional helper function for normalized values
    public static int mpNorm(int x, int p)
    {
        int halfP = p >>> 1;
        int diff = halfP - x;//(((long)(halfP - x)) & 0xFFFFFFFFL);
        int mask = tbmask(diff);
        return x - (p & mask);
    }

    // Rebuild integer from CRT representation
    public static void zintRebuildCRT(int[] xx, int xOff, int xlen, int n, int numSets,
                                      int normalizeSigned, int[] tmp, int tmpOff, SmallPrime[] primes)
    {
        int uu = 0;
        tmp[tmpOff] = (int)primes[0].p;

        for (int u = 1; u < xlen; u++)
        {
            SmallPrime prime = primes[u];
            int p = (int)prime.p;
            int p0i = (int)prime.p0i;
            int R2 = (int)prime.R2;
            int s = (int)prime.s;
            uu += n;
            int kk = 0;

            for (int k = 0; k < numSets; k++)
            {
                for (int v = 0; v < n; v++)
                {
                    int xp = xx[xOff + kk + v + uu];
                    int xq = zintModSmallUnsigned(xx, xOff + kk + v, u, n, p, p0i, R2);

                    int xr = mpMontyMul(s, mpSub(xp, xq, p), p, p0i);
                    zintAddMulSmall(xx, xOff + kk + v, u, n, tmp, tmpOff, xr);
                }
                kk += n * xlen;
            }
            tmp[tmpOff + u] = zintMulSmall(tmp, tmpOff, u, p);
        }

        if (normalizeSigned != 0)
        {
            int kk = 0;
            for (int k = 0; k < numSets; k++)
            {
                for (int v = 0; v < n; v++)
                {
                    zintNormZero(xx, xOff + kk + v, xlen, n, tmp, tmpOff);
                }
                kk += n * xlen;
            }
        }
    }

    private static int mpHalf(int a, int p)
    {
        return (a + (p & -(a & 1))) >>> 1;
    }

    // Add multiple of small integer to big integer
    private static void zintAddMulSmall(int[] x, int offset, int len,
                                        int stride, int[] y, int yOff, int s)
    {
        long cc = 0;
        int index = offset;

        for (int u = 0; u < len; u++)
        {
            long xw = x[index];
            long yw = y[yOff + u];
            long z = yw * s + xw + cc;

            x[index] = (int)z & 0x7FFFFFFF;
            cc = z >>> 31;
            index += stride;
        }
        x[index] = (int)cc;
    }

    // Multiply big integer by small integer
    private static int zintMulSmall(int[] x, int xOff, int len, int s)
    {
        long cc = 0;
        for (int u = 0; u < len; u++)
        {
            long z = (long)x[xOff + u] * s + cc;
            x[xOff + u] = (int)z & 0x7FFFFFFF;
            cc = z >>> 31;
        }
        return (int)cc;
    }

    public static void zintNormZero(int[] x, int offset, int len, int stride, int[] p, int pOff)
    {
        /*
         * Compare x with p/2. We use the shifted version of p, and p
         * is odd, so we really compare with (p-1)/2; we want to perform
         * the subtraction if and only if x > (p-1)/2.
         */
        int r = 0;
        int bb = 0;

        // Start from the most significant word and move backwards
        int currentOffset = offset + (len - 1) * stride;
        for (int u = len - 1; u >= 0; u--)
        {
            /*
             * Get the two words to compare in wx and wp (both over
             * 31 bits exactly).
             */
            int wx = x[currentOffset];
            int wp = (p[pOff + u] >>> 1) | (bb << 30);  // Use unsigned shift
            bb = p[pOff + u] & 1;

            /*
             * We set cc to -1, 0 or 1, depending on whether wp is
             * lower than, equal to, or greater than wx.
             */
            int cc = wp - wx;
            cc = ((-cc) >>> 31) | -(cc >>> 31);  // Use unsigned shifts

            /*
             * If r != 0 then it is either 1 or -1, and we keep its
             * value. Otherwise, if r = 0, then we replace it with cc.
             */
            r |= cc & ((r & 1) - 1);

            // Move to next (less significant) word
            currentOffset -= stride;
        }

        /*
         * At this point, r = -1, 0 or 1, depending on whether (p-1)/2
         * is lower than, equal to, or greater than x. We thus want to
         * do the subtraction only if r = -1.
         */
        int cc = 0;
        int m = tbmask(r);

        // Reset to least significant word
        currentOffset = offset;
        for (int j = 0; j < len; j++)
        {
            int xw = x[currentOffset];
            int w = xw - p[pOff + j] - cc;
            cc = w >>> 31;  // Use unsigned shift for carry
            xw ^= ((w & 0x7FFFFFFF) ^ xw) & m;
            x[currentOffset] = xw;
            currentOffset += stride;
        }
    }

    // Fixed-point comparison
    public static boolean fxrLt(long x, long y)
    {
        return Longs.compareUnsigned(x, y) < 0;
    }

    // Main recursive step for (f,g) computation
    public static void makeFgStep(HawkParameters prof, int lognTop, int depth,
                                  int[] tmp, int tmpOff)
    {
        int logn = lognTop - depth;
        int n = 1 << logn;
        int hn = n >>> 1;
        int slen = prof.maxBlSmall[depth];
        int tlen = prof.maxBlSmall[depth + 1];

        /*
         * Layout:
         *   fd    output f' (hn*tlen)
         *   gd    output g' (hn*tlen)
         *   fs    source (n*slen)
         *   gs    source (n*slen)
         *   t1    NTT support (n)
         *   t2    extra (max(n, slen - n))
         */
        int fdStart = 0;
        int gdStart = fdStart + hn * tlen;
        int fsStart = gdStart + hn * tlen;
        int gsStart = fsStart + n * slen;
        int t1Start = gsStart + n * slen;
        int t2Start = t1Start + n;

        // Move source polynomials to their segments
        System.arraycopy(tmp, tmpOff, tmp, tmpOff + fsStart, 2 * n * slen);

        /*
         * First slen words: we use the input values directly, and apply
         * inverse NTT as we go, so that we get the sources in RNS (non-NTT).
         */
        for (int u = 0; u < slen; u++)
        {
            SmallPrime prime = PRIMES[u];
            int p = (int)prime.p;
            int p0i = (int)prime.p0i;
            int R2 = (int)prime.R2;

            int fsOffset = fsStart + u * n;
            int gsOffset = gsStart + u * n;
            int fdOffset = fdStart + u * hn;
            int gdOffset = gdStart + u * hn;

            // Compute f' and g' for this prime
            for (int v = 0; v < hn; v++)
            {
                int idx1 = 2 * v;
                int idx2 = 2 * v + 1;

                int fProd = mpMontyMul(tmp[tmpOff + fsOffset + idx1], tmp[tmpOff + fsOffset + idx2], (int)p, (int)p0i);
                tmp[tmpOff + fdOffset + v] = mpMontyMul(fProd, (int)R2, (int)p, (int)p0i);

                int gProd = mpMontyMul(tmp[tmpOff + gsOffset + idx1], tmp[tmpOff + gsOffset + idx2], (int)p, (int)p0i);
                tmp[tmpOff + gdOffset + v] = mpMontyMul(gProd, (int)R2, (int)p, (int)p0i);
            }

            // Apply inverse NTT to source polynomials
            mpMkigm(logn, tmp, tmpOff + t1Start, (int)prime.ig, p, p0i);
            mpINTT(logn, tmp, tmpOff + fsOffset, tmp, tmpOff + t1Start, p, p0i);
            mpINTT(logn, tmp, tmpOff + gsOffset, tmp, tmpOff + t1Start, p, p0i);
        }

        /*
         * Now that fs and gs are in RNS, rebuild their plain integer
         * coefficients.
         */
        zintRebuildCRT(tmp, tmpOff + fsStart, slen, n, 2, 1, tmp, tmpOff + t1Start, PRIMES);

        // Process remaining primes
        for (int u = slen; u < tlen; u++)
        {
            SmallPrime prime = PRIMES[u];
            int p = (int)prime.p;
            int p0i = (int)prime.p0i;
            int R2 = (int)prime.R2;
            int Rx = mpRx31(slen, p, p0i, R2);

            int fdOffset = fdStart + u * hn;
            int gdOffset = gdStart + u * hn;
            int[] gm = new int[n];
            mpMkgm(logn, gm, 0, (int)prime.g, p, p0i);
            // Process f polynomial
            for (int v = 0; v < n; v++)
            {
                tmp[tmpOff + t2Start + v] = zintModSmallSigned(tmp, tmpOff + fsStart + v, slen, n, p, p0i, R2, Rx);
            }

            mpNTT(logn, tmp, tmpOff + t2Start, gm, 0, p, p0i);
            for (int v = 0; v < hn; v++)
            {
                int prod = mpMontyMul(tmp[tmpOff + t2Start + 2 * v], tmp[tmpOff + t2Start + 2 * v + 1], p, p0i);
                tmp[tmpOff + fdOffset + v] = mpMontyMul(prod, R2, p, p0i);
            }

            // Process g polynomial
            for (int v = 0; v < n; v++)
            {
                tmp[tmpOff + t2Start + v] = zintModSmallSigned(tmp, tmpOff + gsStart + v, slen, n, p, p0i, R2, Rx);
            }
            mpNTT(logn, tmp, tmpOff + t2Start, gm, 0, p, p0i);
            for (int v = 0; v < hn; v++)
            {
                int prod = mpMontyMul(tmp[tmpOff + t2Start + 2 * v], tmp[tmpOff + t2Start + 2 * v + 1], p, p0i);
                tmp[tmpOff + gdOffset + v] = mpMontyMul(prod, R2, p, p0i);
            }
        }

        // Move results to front of buffer
        //System.arraycopy(tmp, tmpOff + fdStart, tmp, tmpOff, 2 * hn * tlen);
    }

    // Compute inverse NTT roots
    public static void mpMkigm(int logn, int[] igm, int igmOff, int ig, int p, int p0i)
    {
        int n = 1 << logn;


        // Precompute ig^(2^(10-logn))
        for (int j = logn; j < 10; j++)
        {
            ig = mpMontyMul(ig, ig, p, p0i);
        }

        int k = 10 - logn;
        int x2 = mpHR(p);
        int u = 0;
        while (true)
        {
            int v = REV10[u << k];
            igm[v + igmOff] = x2;
            //System.out.println(u + " " + v + " " + x2);
            u++;
            if (u >= n)
            {
                break;
            }
            x2 = mpMontyMul(x2, ig, p, p0i);
        }
    }

    // Compute R/2 = 2^31 mod p
    public static int mpHR(int p)
    {
        return (1 << 31) - p;
    }

    // Inverse NTT with offset support
    public static void mpINTT(int logn, int[] a, int offset,
                              int[] igm, int igmOff, int p, int p0i)
    {
//        if (logn == 0)
//        {
//            return;
//        }

        int n = 1 << logn;
        int t = 1;
        for (int lm = 0; lm < logn; lm++)
        {
            int hm = 1 << (logn - 1 - lm);
            int dt = t << 1;
            int v0 = 0;
            for (int u = 0; u < hm; u++)
            {
                int s = igm[igmOff + u + hm];
                for (int v = 0; v < t; v++)
                {
                    int idx1 = offset + v0 + v;
                    int idx2 = idx1 + t;

                    int x1 = a[idx1];
                    int x2 = a[idx2];

                    a[idx1] = mpHalf(mpAdd(x1, x2, p), p);
                    a[idx2] = mpMontyMul(mpSub(x1, x2, p), s, p, p0i);
                }
                v0 += dt;
            }
            t = dt;
        }
    }

    // Compute R = 2^(31*len) mod p
    private static int mpRx31(int e, int p, int p0i, int R2)
    {
        /* x <- 2^63 mod p = Montgomery representation of 2^31 */
        int x = mpHalf(R2, p);
        int d = 1;

        while (true)
        {
            if ((e & 1) != 0)
            {
                d = mpMontyMul(d, x, p, p0i);
            }
            e >>>= 1;  // Use unsigned shift
            if (e == 0)
            {
                return d;
            }
            x = mpMontyMul(x, x, p, p0i);
        }
    }

    // Compute signed big integer modulo p
    public static int zintModSmallSigned(int[] d, int offset, int len,
                                         int stride, int p, int p0i,
                                         int R2, int Rx)
    {
        if (len == 0)
        {
            return 0;
        }

        // First compute the unsigned modulus
        int z = zintModSmallUnsigned(d, offset, len, stride, p, p0i, R2);

        // Check the sign bit (bit 30) of the most significant word
        int signBit = d[offset + (len - 1) * stride] >>> 30;  // Use unsigned shift
        int adjustment = Rx & (-signBit);  // Rx if sign bit set, 0 otherwise

        // Subtract adjustment to handle negative numbers
        z = mpSub(z, adjustment, p);

        return z;
    }

    // Compute unsigned big integer modulo p
    private static int zintModSmallUnsigned(int[] d, int offset, int len,
                                            int stride, int p, int p0i,
                                            int R2)
    {
        /*
         * Algorithm: we inject words one by one, starting with the high
         * word. Each step is:
         *  - multiply x by 2^31
         *  - add new word
         */
        int x = 0;
        int z = mpHalf(R2, p);

        // Start from the highest word (end of the array)
        int index = offset + len * stride;

        for (int u = len; u > 0; u--)
        {
            index -= stride;
            int w = d[index] - p;
            w += p & tbmask(w);
            x = mpMontyMul(x, z, p, p0i);
            x = mpAdd(x, w, p);
        }

        return x;
    }

    // Main function to regenerate f and g polynomials
    public static void hawkRegenFg(int logn, byte[] f, byte[] g, byte[] seed)
    {
        switch (logn)
        {
        case 8:
            regenFg8(f, g, seed);
            break;
        case 9:
            regenFg9(f, g, seed);
            break;
        case 10:
            regenFg10(f, g, seed);
            break;
        default:
            throw new IllegalArgumentException("Unsupported logn value: " + logn);
        }
    }

    // Helper function to convert little-endian bytes to long
    private static long dec64le(byte[] buf)
    {
        long value = 0;
        for (int i = 0; i < 8; i++)
        {
            value |= (long)(buf[i] & 0xFF) << (8 * i);
        }
        return value;
    }

    // Regeneration for logn=8
    private static void regenFg8(byte[] f, byte[] g, byte[] seed)
    {
        int seedLen = 16;
        for (int j = 0; j < 4; j++)
        {
            SHAKEDigest shake = new SHAKEDigest(256);
            shake.update(seed, 0, seedLen);
            byte jx = (byte)j;
            shake.update(jx);

            // In BouncyCastle, we don't need explicit flip - just extract
            for (int u = 0; u < 512; u += 64)
            {
                byte[] qb = new byte[8];
                shake.doOutput(qb, 0, 8);
                long q = dec64le(qb);

                // Bit manipulation to convert to coefficients
                q = (q & 0x5555555555555555L) + ((q >>> 1) & 0x5555555555555555L);
                q = (q & 0x3333333333333333L) + ((q >>> 2) & 0x3333333333333333L);

                byte[] vv = new byte[16];
                for (int i = 0; i < 16; i++)
                {
                    int val = (int)(q & 0x0F) - 2;
                    vv[i] = (byte)val;
                    q >>>= 4;
                }

                if (u < 256)
                {
                    int destPos = u + (j << 4);
                    System.arraycopy(vv, 0, f, destPos, 16);
                }
                else
                {
                    int destPos = (u - 256) + (j << 4);
                    System.arraycopy(vv, 0, g, destPos, 16);
                }
            }
        }
    }

    // Regeneration for logn=9
    private static void regenFg9(byte[] f, byte[] g, byte[] seed)
    {
        int seedLen = 24;
        for (int j = 0; j < 4; j++)
        {
            SHAKEDigest shake = new SHAKEDigest(256);
            shake.update(seed, 0, seedLen);
            byte jx = (byte)j;
            shake.update(jx);

            for (int u = 0; u < 1024; u += 32)
            {
                byte[] qb = new byte[8];
                shake.doOutput(qb, 0, 8);
                long q = dec64le(qb);

                // Additional bit manipulation step
                q = (q & 0x5555555555555555L) + ((q >>> 1) & 0x5555555555555555L);
                q = (q & 0x3333333333333333L) + ((q >>> 2) & 0x3333333333333333L);
                q = (q & 0x0F0F0F0F0F0F0F0FL) + ((q >>> 4) & 0x0F0F0F0F0F0F0F0FL);

                byte[] vv = new byte[8];
                for (int i = 0; i < 8; i++)
                {
                    int val = (int)(q & 0xFF) - 4;
                    vv[i] = (byte)val;
                    q >>>= 8;
                }

                if (u < 512)
                {
                    int destPos = u + (j << 3);
                    System.arraycopy(vv, 0, f, destPos, 8);
                }
                else
                {
                    int destPos = (u - 512) + (j << 3);
                    System.arraycopy(vv, 0, g, destPos, 8);
                }
            }
        }
    }

    // Regeneration for logn=10
    private static void regenFg10(byte[] f, byte[] g, byte[] seed)
    {
        int seedLen = 40;
        for (int j = 0; j < 4; j++)
        {
            SHAKEDigest shake = new SHAKEDigest(256);
            shake.update(seed, 0, seedLen);
            byte jx = (byte)j;
            shake.update(jx);

            for (int u = 0; u < 2048; u += 16)
            {
                byte[] qb = new byte[8];
                shake.doOutput(qb, 0, 8);
                long q = dec64le(qb);

                // Additional bit manipulation steps
                q = (q & 0x5555555555555555L) + ((q >>> 1) & 0x5555555555555555L);
                q = (q & 0x3333333333333333L) + ((q >>> 2) & 0x3333333333333333L);
                q = (q & 0x0F0F0F0F0F0F0F0FL) + ((q >>> 4) & 0x0F0F0F0F0F0F0F0FL);
                q = (q & 0x00FF00FF00FF00FFL) + ((q >>> 8) & 0x00FF00FF00FF00FFL);

                byte[] vv = new byte[4];
                for (int i = 0; i < 4; i++)
                {
                    int val = (int)(q & 0xFFFF) - 8;
                    vv[i] = (byte)val;
                    q >>>= 16;
                }

                if (u < 1024)
                {
                    int destPos = u + (j << 2);
                    System.arraycopy(vv, 0, f, destPos, 4);
                }
                else
                {
                    int destPos = (u - 1024) + (j << 2);
                    System.arraycopy(vv, 0, g, destPos, 4);
                }
            }
        }
    }

    // Compute the squared norm of a polynomial
    public static int polySqNorm(int logn, byte[] f, int fOff)
    {
        int n = 1 << logn;
        int s = 0;
        for (int u = 0; u < n; u++)
        {
            int x = f[u + fOff];
            s += x * x;
        }
        return s;
    }

    // Fast Fourier Transform for fixed-point complex numbers
    public static void vectFFT(int logn, long[] f, int fOff)
    {
        int hn = 1 << (logn - 1);
        int t = hn;

        for (int lm = 1; lm < logn; lm++)
        {
            int m = 1 << lm;
            int ht = t >>> 1;
            int j0 = 0;
            int hm = m >>> 1;

            for (int i = 0; i < hm; i++)
            {
                // Get twiddle factor from precomputed table
                Fxc s = Utils.GM_TAB[m + i];
                long sRe = s.re;
                long sIm = s.im;

                for (int j = j0; j < j0 + ht; j++)
                {
                    // Extract complex numbers
                    long xRe = f[fOff + j];
                    long xIm = f[fOff + j + hn];
                    long yRe = f[fOff + j + ht];
                    long yIm = f[fOff + j + ht + hn];

                    // y * s (complex multiplication) using Karatsuba-like trick from fxcMul
                    long z0 = fxrMulInline(sRe, yRe);
                    long z1 = fxrMulInline(sIm, yIm);
                    long z2 = fxrMulInline(sRe + sIm, yRe + yIm);
                    long yMulRe = z0 - z1;
                    long yMulIm = z2 - (z0 + z1);

                    // Butterfly: (x + y*s, x - y*s)
                    f[fOff + j]           = xRe + yMulRe;
                    f[fOff + j + hn]      = xIm + yMulIm;
                    f[fOff + j + ht]      = xRe - yMulRe;
                    f[fOff + j + ht + hn] = xIm - yMulIm;
                }
                j0 += t;
            }
            t = ht;
        }
    }

    // Inline fixed-point multiplication: returns the upper 64 bits of x*y (Q32.32 * Q32.32 >> 32).
    private static long fxrMulInline(long x, long y)
    {
        int xh = (int)(x >> 32);
        int yh = (int)(y >> 32);
        long xl = x & 0xFFFFFFFFL;
        long yl = y & 0xFFFFFFFFL;
        long z0 = (xl * yl) >>> 32;
        long z1 = xl * yh;
        long z2 = yl * xh;
        long z3 = ((long)xh * yh) << 32;
        return z0 + z1 + z2 + z3;
    }

    // Inverse Fast Fourier Transform
    public static void vectIFFT(int logn, long[] f, int fOff)
    {
        int hn = 1 << (logn - 1);
        int ht = 1;

        for (int lm = logn - 1; lm > 0; lm--)
        {
            int m = 1 << lm;
            int t = ht << 1;
            int j0 = 0;
            int hm = m >>> 1;

            for (int i = 0; i < hm; i++)
            {
                Fxc gm = Utils.GM_TAB[m + i];
                long sRe = gm.re;
                long sIm = -gm.im;  // conjugate

                for (int j = j0; j < j0 + ht; j++)
                {
                    long xRe = f[j + fOff];
                    long xIm = f[j + fOff + hn];
                    long yRe = f[j + fOff + ht];
                    long yIm = f[j + fOff + ht + hn];

                    // z1 = (x + y) / 2 (matches fxrDiv2e(.,1): (v+1)>>1)
                    long z1Re = (xRe + yRe + 1) >> 1;
                    long z1Im = (xIm + yIm + 1) >> 1;
                    f[j + fOff]      = z1Re;
                    f[j + fOff + hn] = z1Im;

                    // z2 = s * ((x - y) / 2)
                    long dRe = (xRe - yRe + 1) >> 1;
                    long dIm = (xIm - yIm + 1) >> 1;
                    long w0 = fxrMulInline(sRe, dRe);
                    long w1 = fxrMulInline(sIm, dIm);
                    long w2 = fxrMulInline(sRe + sIm, dRe + dIm);
                    f[j + fOff + ht]      = w0 - w1;
                    f[j + fOff + ht + hn] = w2 - (w0 + w1);
                }
                j0 += t;
            }
            ht = t;
        }
    }

    // Error codes
    public static final int SOLVE_OK = 0;
    public static final int SOLVE_ERR_GCD = -1;
    public static final int SOLVE_ERR_REDUCE = -2;
    public static final int SOLVE_ERR_LIMIT = -3;


    public static void make_fg_intermediate(HawkParameters prof, int logn_top,
                                            byte[] f, int fOff, byte[] g, int gOff, int depth, int[] tmp, int tmpOff)
    {
        makeFgZero(logn_top, f, fOff, g, gOff, tmp, tmpOff);
        for (int d = 0; d < depth; d++)
        {
            makeFgStep(prof, logn_top, d, tmp, tmpOff);
        }
    }

    public static int make_fg_deepest(HawkParameters prof, int logn_top,
                                      byte[] f, int fOff, byte[] g, int gOff, int[] tmp, int tmpOff, int sav_off)
    {
        makeFgZero(logn_top, f, fOff, g, gOff, tmp, tmpOff);
        int r = 1;

        int n = 1 << logn_top;
        int b = 0;
        for (int u = 0; u < n; u++)
        {
            b |= tmp[tmpOff + u] - 1;
        }
        r = 1 - (b >>> 31);

        for (int d = 0; d < logn_top; d++)
        {
            makeFgStep(prof, logn_top, d, tmp, tmpOff);

            int d2 = d + 1;
            if (d2 < logn_top && d2 >= prof.minSaveFg[logn_top])
            {
                int slen = prof.maxBlSmall[d2];
                int fglen = slen << (logn_top + 1 - d2);
                sav_off -= fglen;
                System.arraycopy(tmp, tmpOff, tmp, sav_off, fglen);
            }
        }
        return r;
    }

    public static int zint_bezout(int[] tmp, int u, int v, int x, int y, int len, int tmpOff)
    {
        if (len == 0)
        {
            return 0;
        }
        /*
         * Algorithm is basically the optimized binary GCD as described in:
         *    https://eprint.iacr.org/2020/972
         * The paper shows that with registers of size 2*k bits, one can
         * do k-1 inner iterations and get a reduction by k-1 bits. In
         * fact, it also works with registers of 2*k-1 bits (though not
         * 2*k-2; the "upper half" of the approximation must have at
         * least one extra bit). Here, we want to perform 31 inner
         * iterations (since that maps well to Montgomery reduction with
         * our 31-bit words) so we must use 63-bit approximations.
         *
         * We also slightly expand the original algorithm by maintaining
         * four coefficients (u0, u1, v0 and v1) instead of the two
         * coefficients (u, v), because we want a full Bezout relation,
         * not just a modular inverse.
         *
         * We set up integers u0, v0, u1, v1, a and b. Throughout the
         * algorithm, they maintain the following invariants:
         *   a = x*u0 - y*v0
         *   b = x*u1 - y*v1
         *   0 <= a <= x
         *   0 <= b <= y
         *   0 <= u0 < y
         *   0 <= v0 < x
         *   0 <= u1 <= y
         *   0 <= v1 < x
         */
        int u0 = tmpOff;
        int v0 = u0 + len;
        int u1 = u;
        int v1 = v;
        int a = v0 + len;
        int b = a + len;

        /*
         * We'll need the Montgomery reduction coefficients.
         */
        int x0i = mp_ninv31(tmp[x]);
        int y0i = mp_ninv31(tmp[y]);

        /*
         * Initial values:
         *   a = x   u0 = 1   v0 = 0
         *   b = y   u1 = y   v1 = x - 1
         * Note that x is odd, so computing x-1 is easy.
         */
        System.arraycopy(tmp, x, tmp, a, len);
        System.arraycopy(tmp, y, tmp, b, len);
        tmp[u0] = 1;
        Arrays.fill(tmp, u0 + 1, u0 + len, 0);
        Arrays.fill(tmp, v0, v0 + len, 0);
        System.arraycopy(tmp, y, tmp, u1, len);
        System.arraycopy(tmp, x, tmp, v1, len);
        tmp[v1]--;

        /*
         * Each input operand may be as large as 31*len bits, and we
         * reduce the total length by at least 31 bits at each iteration.
         */
        for (int num = 62 * len + 31; num >= 30; num -= 31)
        {
            /*
             * Extract the top 32 bits of a and b: if j is such that:
             *   2^(j-1) <= max(a,b) < 2^j
             * then we want:
             *   xa = (2^31)*floor(a / 2^(j-32)) + (a mod 2^31)
             *   xb = (2^31)*floor(a / 2^(j-32)) + (b mod 2^31)
             * (if j < 63 then xa = a and xb = b).
             */
            int c0 = 0xFFFFFFFF;
            int c1 = 0xFFFFFFFF;
            int cp = 0xFFFFFFFF;
            int a0 = 0;
            int a1 = 0;
            int b0 = 0;
            int b1 = 0;
            int j = len;
            while (j-- > 0)
            {
                int aw = tmp[a + j];
                int bw = tmp[b + j];
                a1 ^= c1 & (a1 ^ aw);
                a0 ^= c0 & (a0 ^ aw);
                b1 ^= c1 & (b1 ^ bw);
                b0 ^= c0 & (b0 ^ bw);
                cp = c0;
                c0 = c1;
                int cond = ((aw | bw) + 0x7FFFFFFF) >>> 31;
                c1 &= cond - 1;
            }

            /*
             * Possible situations:
             *   cp = 0, c0 = 0, c1 = 0
             *     j >= 63, top words of a and b are in a0:a1 and b0:b1
             *     (a1 and b1 are highest, a1|b1 != 0)
             *
             *   cp = -1, c0 = 0, c1 = 0
             *     32 <= j <= 62, a0:a1 and b0:b1 contain a and b, exactly
             *
             *   cp = -1, c0 = -1, c1 = 0
             *     j <= 31, a0 and a1 both contain a, b0 and b1 contain b
             *
             * When j >= 63, we align the top words to ensure that we get
             * the full 32 bits. We also take care to always call
             * lzcnt() with a non-zero operand.
             */
            int s = lzcnt_nonzero(a1 | b1 | ((cp & c0) >>> 1));
            int ha = (a1 << s) | (a0 >>> (31 - s));
            int hb = (b1 << s) | (b0 >>> (31 - s));

            /*
             * If j <= 62, then we instead use the non-aligned bits.
             */
            ha ^= (cp & (ha ^ a1));
            hb ^= (cp & (hb ^ b1));

            /*
             * If j <= 31, then all of the above was bad, and we simply
             * clear the upper bits.
             */
            ha &= ~c0;
            hb &= ~c0;

            /*
             * Assemble the approximate values xa and xb (63 bits each).
             */
            long xa = ((ha & 0xFFFFFFFFL) << 31) | (tmp[a] & 0x7FFFFFFFL);
            long xb = ((hb & 0xFFFFFFFFL) << 31) | (tmp[b] & 0x7FFFFFFFL);

            /*
             * Compute reduction factors:
             *   a' = a*pa + b*pb
             *   b' = a*qa + b*qb
             * such that a' and b' are both multiples of 2^31, but are
             * only marginally larger than a and b.
             * Each coefficient is in the -(2^31-1)..+2^31 range. To keep
             * them on 32-bit values, we compute pa+(2^31-1)[] and so on.
             */
            long fg0 = 1;
            long fg1 = 1L << 32;
            for (int i = 0; i < 31; i++)
            {
                long a_odd = -(xa & 1);
                long dx = xa - xb;
                dx = dx >> 63;
                long swap = a_odd & dx;
                long t1 = swap & (xa ^ xb);
                xa ^= t1;
                xb ^= t1;
                long t2 = swap & (fg0 ^ fg1);
                fg0 ^= t2;
                fg1 ^= t2;
                xa -= a_odd & xb;
                fg0 -= a_odd & fg1;
                xa >>= 1;
                fg1 <<= 1;
            }

            fg0 += 0x7FFFFFFF7FFFFFFFL;
            fg1 += 0x7FFFFFFF7FFFFFFFL;
            long f0 = (fg0 & 0xFFFFFFFFL) - 0x7FFFFFFFL;
            long g0 = (fg0 >>> 32) - 0x7FFFFFFFL;
            long f1 = (fg1 & 0xFFFFFFFFL) - 0x7FFFFFFFL;
            long g1 = (fg1 >>> 32) - 0x7FFFFFFFL;

            /*
             * Apply the update factors.
             */
            int negab = zint_co_reduce(tmp, a, b, len, f0, g0, f1, g1);
            f0 -= (f0 + f0) & -(negab & 1);
            g0 -= (g0 + g0) & -(negab & 1);
            f1 -= (f1 + f1) & -((negab >>> 1) & 1);
            g1 -= (g1 + g1) & -((negab >>> 1) & 1);
            zint_co_reduce_mod(tmp, u0, u1, y, len, y0i, f0, g0, f1, g1);
            zint_co_reduce_mod(tmp, v0, v1, x, len, x0i, f0, g0, f1, g1);
        }

        /*
         * b contains GCD(x,y), provided that x and y were indeed odd.
         * Result is correct if the GCD is 1.
         */
        int r = tmp[b] ^ 1;
        for (int j = 1; j < len; j++)
        {
            r |= tmp[b + j];
        }
        r |= (tmp[x] & tmp[y] & 1) ^ 1;
        return 1 - (r >>> 31);
    }

    private static int mp_ninv31(int x)
    {
        int y = 2 - x;
        y *= 2 - x * y;
        y *= 2 - x * y;
        y *= 2 - x * y;
        y *= 2 - x * y;
        return -y & 0x7FFFFFFF;
    }

    private static int lzcnt_nonzero(int x)
    {
        return Integers.numberOfLeadingZeros(x);
    }

    private static int zint_co_reduce(int[] tmp, int a, int b, int len,
                                      long xa, long xb, long ya, long yb)
    {
        long cca = 0;
        long ccb = 0;
        for (int u = 0; u < len; u++)
        {
            int wa = tmp[a + u];
            int wb = tmp[b + u];
            long za = wa * xa + wb * xb + cca;
            long zb = wa * ya + wb * yb + ccb;
            if (u > 0)
            {
                tmp[a + u - 1] = (int)za & 0x7FFFFFFF;
                tmp[b + u - 1] = (int)zb & 0x7FFFFFFF;
            }
            cca = za >> 31;
            ccb = zb >> 31;
        }
        tmp[a + len - 1] = (int)cca & 0x7FFFFFFF;
        tmp[b + len - 1] = (int)ccb & 0x7FFFFFFF;

        int nega = (int)(cca >>> 63);
        int negb = (int)(ccb >>> 63);
        zint_negate(tmp, a, len, nega);
        zint_negate(tmp, b, len, negb);
        return nega | (negb << 1);
    }

    private static void zint_co_reduce_mod(int[] tmp, int a, int b, int m, int len,
                                           int m0i, long xa, long xb, long ya, long yb)
    {
        long cca = 0;
        long ccb = 0;
        int fa = (int)((tmp[a] * xa + tmp[b] * xb) * m0i) & 0x7FFFFFFF;
        int fb = (int)((tmp[a] * ya + tmp[b] * yb) * m0i) & 0x7FFFFFFF;
        for (int u = 0; u < len; u++)
        {
            int wa = tmp[a + u];
            int wb = tmp[b + u];
            int mw = tmp[m + u];
            long za = wa * xa + wb * xb + mw * (long)fa + cca;
            long zb = wa * ya + wb * yb + mw * (long)fb + ccb;
            if (u > 0)
            {
                tmp[a + u - 1] = (int)za & 0x7FFFFFFF;
                tmp[b + u - 1] = (int)zb & 0x7FFFFFFF;
            }
            cca = za >> 31;
            ccb = zb >> 31;
        }
        tmp[a + len - 1] = (int)cca;
        tmp[b + len - 1] = (int)ccb;

        zint_finish_mod(tmp, a, m, len, (int)(cca >>> 63));
        zint_finish_mod(tmp, b, m, len, (int)(ccb >>> 63));
    }

    private static void zint_finish_mod(int[] tmp, int a, int m, int len, int neg)
    {
        int cc = 0;
        for (int u = 0; u < len; u++)
        {
            cc = (tmp[a + u] - tmp[m + u] - cc) >>> 31;
        }

        int xm = -neg >>> 1;
        int ym = -(neg | (1 - cc));
        cc = neg;
        for (int u = 0; u < len; u++)
        {
            int mw = (tmp[m + u] ^ xm) & ym;
            int aw = tmp[a + u] - mw - cc;
            tmp[a + u] = aw & 0x7FFFFFFF;
            cc = aw >>> 31;
        }
    }

    public static void zint_negate(int[] a, int aOff, int len, int ctl)
    {
        int cc = ctl;
        int m = (-ctl) >>> 1; // This generates 0x7FFFFFFF if ctl=1, 0 otherwise
        for (int u = 0; u < len; u++)
        {
            int aw = a[aOff + u];
            aw = (aw ^ m) + cc;
            a[aOff + u] = aw & 0x7FFFFFFF;
            cc = aw >>> 31; // Unsigned shift to get carry
        }
    }

    public static int solve_NTRU_deepest(HawkParameters prof, int logn_top,
                                         byte[] f, int fOff, byte[] g, int gOff, int[] tmp, int tmpOff)
    {
        // Call make_fg_deepest and check return value
        int sav_off = 6 << logn_top;
        /*
         * Get (f,g) at the deepest level (i.e. Res(f,X^n+1) and Res(g,X^n+1)).
         * Obtained (f,g) are in RNS+NTT (since degree n = 1, this is
         * equivalent to RNS).
         */
        int r = make_fg_deepest(prof, logn_top, f, fOff, g, gOff, tmp, tmpOff, sav_off + tmpOff);
        if (r == 0)
        {
            return SOLVE_ERR_GCD;
        }

        int len = prof.maxBlSmall[logn_top];
        /*
         * Reorganize memory:
         *    Fp   output F (len)
         *    Gp   output G (len)
         *    fp   Res(f,X^n+1) (len)
         *    gp   Res(g,X^n+1) (len)
         *    t1   rest of temporary
         */
        int Fp = tmpOff;
        int Gp = Fp + len;
        int fp = Gp + len;
        int gp = fp + len;
        int t1 = gp + len;
        System.arraycopy(tmp, tmpOff, tmp, fp, 2 * len);

        // Rebuild CRT for the resultants
        zintRebuildCRT(tmp, fp, len, 1, 2, 0, tmp, t1, PRIMES);

        // Compute Bezout coefficients
        if (zint_bezout(tmp, Gp, Fp, fp, gp, len, t1) == 0)
        {
            return SOLVE_ERR_GCD;
        }

        // Multiply by q if needed
        if (prof.q != 1)
        {
            if (zint_mul_small(tmp, Fp, len, prof.q) != 0 ||
                zint_mul_small(tmp, Gp, len, prof.q) != 0)
            {
                return SOLVE_ERR_REDUCE;
            }
        }

        return SOLVE_OK;
    }

    // Helper method for multiplication
    public static int zint_mul_small(int[] a, int aOff, int len, int x)
    {
        long cc = 0;
        for (int u = 0; u < len; u++)
        {
            long z = (long)a[aOff + u] * (long)x + cc;
            a[aOff + u] = (int)z & 0x7FFFFFFF;
            cc = z >>> 31;
        }
        return (int)cc;
    }

    private static final int MIN_LOGN_FGNTT = 4;

    /*
     * Solving the NTRU equation, intermediate level.
     * Input is (F,G) from one level deeper (half-degree), in plain
     * representation, at the start of tmp[]; output is (F,G) from this
     * level, written at the start of tmp[].
     *
     * Returned value: 0 on success, a negative error code otherwise.
     */
    public static int solve_NTRU_intermediate(HawkParameters prof, int logn_top,
                                              byte[] f, int fOff, byte[] g, int gOff, int depth, int[] tmp, int tmpOff)
    {
        /*
         * MAX SIZE:
         *    input: 2 * hn * max_bl_small[depth + 1]
         */
        int logn = logn_top - depth;
        int n = 1 << logn;
        int hn = n >>> 1;

        /*
         * slen   size for (f,g) at this level (also output (F,G))
         * llen   size for unreduced (F,G) at this level
         * dlen   size for input (F,G) from deeper level
         * Note: we always have llen >= dlen (constraint enforced in profiles)
         */
        int slen = prof.maxBlSmall[depth];
        int llen = prof.maxBlLarge[depth];
        int dlen = prof.maxBlSmall[depth + 1];

        /*
         * Fd   F from deeper level (dlen*hn)
         * Gd   G from deeper level (dlen*hn)
         * ft   f from this level (slen*n)
         * gt   g from this level (slen*n)
         */
        int Fd = tmpOff;
        int Gd = Fd + dlen * hn;
        int fgt = Gd + dlen * hn;

        // Get (f,g) for this level (in RNS+NTT).
        if (depth < prof.minSaveFg[logn_top])
        {
            make_fg_intermediate(prof, logn_top, f, fOff, g, gOff, depth, tmp, tmpOff + fgt);
        }
        else
        {
            int sav_fg = tmpOff + (6 << logn_top);
            for (int d = prof.minSaveFg[logn_top]; d <= depth; d++)
            {
                sav_fg -= prof.maxBlSmall[d] << (logn_top + 1 - d);
            }
            System.arraycopy(tmp, tmpOff + sav_fg, tmp, fgt, 2 * slen * n);
        }

        /* obsolete
        kg_stats_set_max_size(logn_top, depth,
            2 * dlen * hn + 3 * ((size_t)1 << logn_top));
        */

        /*
         * Move buffers so that we have room for the unreduced (F,G) at
         * this level.
         *   Ft   F from this level (unreduced) (llen*n)
         *   Gt   G from this level (unreduced) (llen*n)
         *   ft   f from this level (slen*n)
         *   gt   g from this level (slen*n)
         *   Fd   F from deeper level (dlen*hn)
         *   Gd   G from deeper level (dlen*hn)
         */
        int Ft = tmpOff;
        int Gt = Ft + llen * n;
        int ft = Gt + llen * n;
        int gt = ft + slen * n;
        Fd = gt + slen * n;
        Gd = Fd + dlen * hn;
        int t1 = Gd + dlen * hn;

        // Move data to new positions
        System.arraycopy(tmp, tmpOff + fgt, tmp, tmpOff + ft, 2 * n * slen);
        System.arraycopy(tmp, tmpOff, tmp, tmpOff + Fd, 2 * hn * dlen);

        /* obsolete
        kg_stats_set_max_size(logn_top, depth,
            2 * llen * n + 2 * slen * n + 2 * dlen * hn);
        */

        /*
         * Convert Fd and Gd to RNS, with output temporarily stored
         * in (Ft, Gt). Fd and Gd have degree hn only; we store the
         * values for each modulus p in the _last_ hn slots of the
         * n-word line for that modulus.
         */
        for (int u = 0; u < llen; u++)
        {
            SmallPrime prime = PRIMES[u];
            int p = (int)prime.p;
            int p0i = (int)prime.p0i;
            int R2 = (int)prime.R2;
            int Rx = mpRx31(dlen, p, p0i, R2);

            int xt = Ft + u * n + hn;
            int yt = Gt + u * n + hn;

            for (int v = 0; v < hn; v++)
            {
                tmp[tmpOff + xt + v] = zintModSmallSigned(tmp, Fd + v, dlen, hn, p, p0i, R2, Rx);
                tmp[tmpOff + yt + v] = zintModSmallSigned(tmp, Gd + v, dlen, hn, p, p0i, R2, Rx);
            }
        }

        /*
         * Fd and Gd are no longer needed.
         */
        t1 = Fd;

        /*
         * Compute (F,G) (unreduced) modulo sufficiently many small primes.
         * We also un-NTT (f,g) as we go; when slen primes have been
         * processed, we obtain (f,g) in RNS, and we apply the CRT to
         * get (f,g) in plain representation.
         */
        for (int u = 0; u < llen; u++)
        {
            /*
             * If we have processed exactly slen primes, then (f,g)
             * are in RNS, and we can rebuild them.
             */
            if (u == slen)
            {
                zintRebuildCRT(tmp, ft, slen, n, 2, 1, tmp, t1, PRIMES);
            }

            SmallPrime prime = PRIMES[u];
            int p = (int)prime.p;
            int p0i = (int)prime.p0i;
            int R2 = (int)prime.R2;
            /*
             * Memory layout: we keep Ft, Gt, ft and gt; we append:
             *   gm    NTT support (n)
             *   igm   iNTT support (n)
             *   fx    temporary f mod p (NTT) (n)
             *   gx    temporary g mod p (NTT) (n)
             */
            int gm = t1;
            int igm = gm + n;
            int fx = igm + n;
            int gx = fx + n;


            mpMkgmigm(logn, tmp, tmpOff + gm, tmp, tmpOff + igm, (int)prime.g, (int)prime.ig, p, p0i);

            if (u < slen)
            {
                System.arraycopy(tmp, tmpOff + ft + u * n, tmp, tmpOff + fx, n);
                System.arraycopy(tmp, tmpOff + gt + u * n, tmp, tmpOff + gx, n);
                mpINTT(logn, tmp, tmpOff + ft + u * n, tmp, tmpOff + igm, p, p0i);
                mpINTT(logn, tmp, tmpOff + gt + u * n, tmp, tmpOff + igm, p, p0i);
            }
            else
            {
                int Rx = mpRx31(slen, p, p0i, R2);
                for (int v = 0; v < n; v++)
                {
                    tmp[tmpOff + fx + v] = zintModSmallSigned(tmp, tmpOff + ft + v, slen, n, p, p0i, R2, Rx);
                    tmp[tmpOff + gx + v] = zintModSmallSigned(tmp, tmpOff + gt + v, slen, n, p, p0i, R2, Rx);
                }
                mpNTT(logn, tmp, tmpOff + fx, tmp, tmpOff + gm, p, p0i);
                mpNTT(logn, tmp, tmpOff + gx, tmp, tmpOff + gm, p, p0i);
            }

            /*
             * We have (F,G) from deeper level in Ft and Gt, in
             * RNS. We apply the NTT modulo p.
             */
            int Fe = Ft + u * n;
            int Ge = Gt + u * n;
            mpNTT(logn - 1, tmp, tmpOff + Fe + hn, tmp, tmpOff + gm, p, p0i);
            mpNTT(logn - 1, tmp, tmpOff + Ge + hn, tmp, tmpOff + gm, p, p0i);

            /*
             * Compute F and G (unreduced) modulo p.
             */
            for (int v = 0; v < hn; v++)
            {
                int fa = tmp[tmpOff + fx + (v << 1)];
                int fb = tmp[tmpOff + fx + (v << 1) + 1];
                int ga = tmp[tmpOff + gx + (v << 1)];
                int gb = tmp[tmpOff + gx + (v << 1) + 1];
                int mFp = mpMontyMul(tmp[tmpOff + Fe + v + hn], R2, p, p0i);
                int mGp = mpMontyMul(tmp[tmpOff + Ge + v + hn], R2, p, p0i);
                tmp[tmpOff + Fe + (v << 1)] = mpMontyMul(gb, mFp, p, p0i);
                tmp[tmpOff + Fe + (v << 1) + 1] = mpMontyMul(ga, mFp, p, p0i);
                tmp[tmpOff + Ge + (v << 1)] = mpMontyMul(fb, mGp, p, p0i);
                tmp[tmpOff + Ge + (v << 1) + 1] = mpMontyMul(fa, mGp, p, p0i);
            }

            /*
             * We want the new (F,G) in RNS only (no NTT).
             */
            mpINTT(logn, tmp, tmpOff + Fe, tmp, tmpOff + igm, p, p0i);
            mpINTT(logn, tmp, tmpOff + Ge, tmp, tmpOff + igm, p, p0i);
        }

        /*
         * Edge case: if slen == llen, then we have not rebuilt (f,g)
         * into plain representation yet, so we do it now.
         */
        if (slen == llen)
        {
            zintRebuildCRT(tmp, ft, slen, n, 2, 1, tmp, t1, PRIMES);
        }

        /*
         * We now have the unreduced (F,G) in RNS. We rebuild their
         * plain representation.
         */
        zintRebuildCRT(tmp, Ft, llen, n, 2, 1, tmp, t1, PRIMES);

        /*
         * We now reduce these (F,G) with Babai's nearest plane
         * algorithm. The reduction conceptually goes as follows:
         *   k <- round((F*adj(f) + G*adj(g))/(f*adj(f) + g*adj(g)))
         *   (F, G) <- (F - k*f, G - k*g)
         * We use fixed-point approximations of (f,g) and (F, G) to get
         * a value k as a small polynomial with scaling; we then apply
         * k on the full-width polynomial. Each iteration "shaves" a
         * a few bits off F and G.
         *
         * We apply the process sufficiently many times to reduce (F, G)
         * to the size of (f, g) with a reasonable probability of success.
         * Since we want full constant-time processing, the number of
         * iterations and the accessed slots work on some assumptions on
         * the sizes of values (sizes have been measured over many samples,
         * and a margin of 5 times the standard deviation).
         */

        /*
         * If depth is at least 2, and we will use the NTT to subtract
         * k*(f,g) from (F,G), then we will need to convert (f,g) to NTT over
         * slen+1 words, which requires an extra word to ft and gt.
         */
        boolean use_sub_ntt = (depth > 1 && logn >= MIN_LOGN_FGNTT);
        if (use_sub_ntt)
        {
            System.arraycopy(tmp, tmpOff + gt, tmp, tmpOff + gt + n, n * slen);
            gt += n;
            t1 += 2 * n;
        }

        /*
         * New layout:
         *   Ft    F from this level (unreduced) (llen*n)
         *   Gt    G from this level (unreduced) (llen*n)
         *   ft    f from this level (slen*n) (+n if use_sub_ntt)
         *   gt    g from this level (slen*n) (+n if use_sub_ntt)
         *   rt3   (n fxr = 2*n)
         *   rt4   (n fxr = 2*n)
         *   rt1   (hn fxr = n)
         */
        int rt3 = t1 / 2;
        int rt4 = rt3 + n;
        int rt1 = rt4 + n;
        int rlen = Math.min(prof.wordWin[depth], slen);
        int blen = slen - rlen;
        int ftb = ft + blen * n;
        int gtb = gt + blen * n;
        int scale_fg = 31 * blen;
        int scale_FG = 31 * llen;

        /*
         * Convert f and g into fixed-point approximations, in rt3 and rt4,
         * respectively. They are scaled down by 2^(scale_fg + scale_x).
         * scale_fg is public (it depends only on the recursion depth), but
         * scale_x comes from a measurement on the actual values of (f,g) and
         * is thus secret.
         *
         * The value scale_x is adjusted so that the largest coefficient is
         * close to, but lower than, some limit t (in absolute value). The
         * limit t is chosen so that f*adj(f) + g*adj(g) does not overflow,
         * i.e. all coefficients must remain below 2^31.
         *
         * Let n be the degree; we know that n <= 2^10. The squared norm
         * of a polynomial is the sum of the squared norms of the
         * coefficients, with the squared norm of a complex number being
         * the product of that number with its complex conjugate. If all
         * coefficients of f are less than t (in absolute value), then
         * the squared norm of f is less than n*t^2. The squared norm of
         * FFT(f) (f in FFT representation) is exactly n times the
         * squared norm of f, so this leads to n^2*t^2 as a maximum
         * bound. adj(f) has the same norm as f. This implies that each
         * complex coefficient of FFT(f) has a maximum squared norm of
         * n^2*t^2 (with a maximally imbalanced polynomial with all
         * coefficient but one being zero). The computation of f*adj(f)
         * exactly is, in FFT representation, the product of each
         * coefficient with its conjugate; thus, the coefficients of
         * f*adj(f), in FFT representation, are at most n^2*t^2.
         *
         * Since we want the coefficients of f*adj(f)+g*adj(g) not to exceed
         * 2^31, we need n^2*t^2 <= 2^30, i.e. n*t <= 2^15. We can adjust t
         * accordingly (called scale_t in the code below). We also need to
         * take care that t must not exceed scale_x. Approximation of f and
         * g are extracted with scale scale_fg + scale_x - scale_t, and
         * later fixed by dividing them by 2^scale_t.
         */
        int scale_xf = poly_max_bitlength(logn, tmp, tmpOff + ftb, rlen);
        int scale_xg = poly_max_bitlength(logn, tmp, tmpOff + gtb, rlen);
        int scale_x = scale_xf;
        scale_x ^= (scale_xf ^ scale_xg) & tbmask(scale_xf - scale_xg);
        int scale_t = 15 - logn;
        scale_t ^= (scale_t ^ scale_x) & tbmask(scale_x - scale_t);
        int scdiff = scale_x - scale_t;

        /*
         * rt is the 64-bit fixed-point view that the reference C aliases over
         * tmp (fxr *rt = (fxr *)tmp). We keep it as a persistent scratch array:
         * poly_big_to_fixed() writes it directly from the int[] data and the
         * fxr_round() loop reads it back into int[] tmp[k_offset], so the two
         * domains never need a bulk re-sync (the FFT scratch region and the
         * int F/G/f/g region are disjoint in the buffer layout).
         */
        long[] rt = new long[(tmp.length - tmpOff) >>> 1];

        poly_big_to_fixed(logn, rt, rt3, tmp, tmpOff + ftb, rlen, scdiff);
        poly_big_to_fixed(logn, rt, rt4, tmp, tmpOff + gtb, rlen, scdiff);

        /*
         * Compute adj(f)/(f*adj(f) + g*adj(g)) into rt3 (FFT).
         * Compute adj(g)/(f*adj(f) + g*adj(g)) into rt4 (FFT).
         */

        vectFFT(logn, rt, rt3);
        vectFFT(logn, rt, rt4);
        vect_norm_fft(logn, rt, rt1, rt, rt3, rt, rt4);
        vect_mul2e(logn, rt, rt3, scale_t);
        vect_mul2e(logn, rt, rt4, scale_t);

        for (int u = 0; u < hn; u++)
        {
            rt[rt3 + u] = fxr_div(rt[rt3 + u], rt[rt1 + u]);
            rt[rt3 + u + hn] = fxr_div(fxr_neg(rt[rt3 + u + hn]), rt[rt1 + u]);
            rt[rt4 + u] = fxr_div(rt[rt4 + u], rt[rt1 + u]);
            rt[rt4 + u + hn] = fxr_div(fxr_neg(rt[rt4 + u + hn]), rt[rt1 + u]);
        }

        /*
         * New layout:
         *   Ft    F from this level (unreduced) (llen*n)
         *   Gt    G from this level (unreduced) (llen*n)
         *   ft    f from this level (slen*n) (+n if use_sub_ntt)
         *   gt    g from this level (slen*n) (+n if use_sub_ntt)
         *   rt3   (n fxr = 2*n)
         *   rt4   (n fxr = 2*n)
         *   rt1   (n fxr = 2*n)     |   k    (n)
         *   rt2   (n fxr = 2*n)     |   t2   (3*n)
         * Exception: at depth == 1, we omit ft and gt:
         *   Ft    F from this level (unreduced) (llen*n)
         *   Gt    G from this level (unreduced) (llen*n)
         *   rt3   (n fxr = 2*n)
         *   rt4   (n fxr = 2*n)
         *   rt1   (n fxr = 2*n)     |   k    (n)
         *   rt2   (n fxr = 2*n)     |   t2   (3*n)
         */
        if (depth == 1)
        {
            t1 = ft;
            int nrt3 = t1 / 2;
            System.arraycopy(rt, rt3, rt, nrt3, 2 * n);
            rt3 = nrt3;
            rt4 = rt3 + n;
            rt1 = rt4 + n;
        }

// Prepare k array
        int k_offset = rt1 * 2;
        int t2 = k_offset + n;
        int rt2 = t2 / 2;
        if (rt2 < rt1 + n)
        {
            rt2 = rt1 + n;
        }

        /*
         * If we are going to use poly_sub_scaled_ntt(), then we convert
         * f and g to the NTT representation. Since poly_sub_scaled_ntt()
         * itself will use more than n*(slen+2) words in t2[], we can do
         * the same here.
         */
        if (use_sub_ntt)
        {
            int gm_offset = t2;
            int tn_offset = gm_offset + n;

            for (int u = 0; u <= slen; u++)
            {
                SmallPrime prime = PRIMES[u];
                int p = (int)prime.p;
                int p0i = (int)prime.p0i;
                int R2 = (int)prime.R2;
                int Rx = mpRx31(slen, p, p0i, R2);

                mpMkgm(logn, tmp, tmpOff + gm_offset, (int)prime.g, p, p0i);

                for (int v = 0; v < n; v++)
                {
                    tmp[tn_offset + v] = zintModSmallSigned(tmp, ft + v, slen, n, p, p0i, R2, Rx);
                }
                mpNTT(logn, tmp, tmpOff + tn_offset, tmp, tmpOff + gm_offset, p, p0i);
                tn_offset += n;
            }

            System.arraycopy(tmp, tmpOff + gm_offset + n, tmp, tmpOff + ft, (slen + 1) * n);
            tn_offset = gm_offset + n;

            for (int u = 0; u <= slen; u++)
            {
                SmallPrime prime = PRIMES[u];
                int p = (int)prime.p;
                int p0i = (int)prime.p0i;
                int R2 = (int)prime.R2;
                int Rx = mpRx31(slen, p, p0i, R2);

                mpMkgm(logn, tmp, tmpOff + gm_offset, (int)prime.g, p, p0i);

                for (int v = 0; v < n; v++)
                {
                    tmp[tmpOff + tn_offset + v] = zintModSmallSigned(tmp, tmpOff + gt + v, slen, n, p, p0i, R2, Rx);
                }
                mpNTT(logn, tmp, tmpOff + tn_offset, tmp, tmpOff + gm_offset, p, p0i);
                tn_offset += n;
            }

            System.arraycopy(tmp, tmpOff + gm_offset + n, tmp, tmpOff + gt, (slen + 1) * n);
        }

        /*
         * Reduce F and G repeatedly.
         */
        int FGlen = llen;
        int loop = 0;
        while (true)
        {
            loop++;
            /*
             * Convert the current F and G into fixed-point. We want
             * to apply scaling scale_FG + scale_x.
             */
            int[] divRem = new int[2];
            DIVREM31(divRem, scale_FG);
            int tlen = divRem[0];
            int toff = divRem[1];

            poly_big_to_fixed(logn, rt, rt1, tmp, tmpOff + Ft + tlen * n, FGlen - tlen, scale_x + toff);
            poly_big_to_fixed(logn, rt, rt2, tmp, tmpOff + Gt + tlen * n, FGlen - tlen, scale_x + toff);

            // Compute k
            vectFFT(logn, rt, rt1);
            vectFFT(logn, rt, rt2);
            vect_mul_fft(logn, rt, rt1, rt, rt3);
            vect_mul_fft(logn, rt, rt2, rt, rt4);
            vect_add(logn, rt, rt2, rt, rt1);
            vectIFFT(logn, rt, rt2);
            for (int u = 0; u < n; u++)
            {
                tmp[tmpOff + k_offset + u] = fxr_round(rt[rt2 + u]);
            }

            // Apply reduction
            int scale_k = scale_FG - scale_fg;
            if (depth == 1)
            {
                poly_sub_kfg_scaled_depth1(logn_top, tmp, tmpOff + Ft, tmp, tmpOff + Gt, FGlen, tmp, tmpOff + k_offset, scale_k, f, g, tmp, tmpOff + t2);
            }
            else if (use_sub_ntt)
            {
                poly_sub_scaled_ntt(logn, tmp, tmpOff + Ft, FGlen, tmp, tmpOff + ft, slen, tmp, tmpOff + k_offset, scale_k, tmp, tmpOff + t2);
                poly_sub_scaled_ntt(logn, tmp, tmpOff + Gt, FGlen, tmp, tmpOff + gt, slen, tmp, tmpOff + k_offset, scale_k, tmp, tmpOff + t2);
            }
            else
            {
                poly_sub_scaled(logn, tmp, tmpOff + Ft, FGlen, tmp, tmpOff + ft, slen, tmp, tmpOff + k_offset, scale_k);
                poly_sub_scaled(logn, tmp, tmpOff + Gt, FGlen, tmp, tmpOff + gt, slen, tmp, tmpOff + k_offset, scale_k);
            }
            // Check if reduction is complete
            if (scale_FG <= scale_fg)
            {
                break;
            }
            if (scale_FG <= (scale_fg + prof.reduceBits))
            {
                scale_FG = scale_fg;
            }
            else
            {
                scale_FG -= prof.reduceBits;
            }
            while (FGlen > slen && 31 * (FGlen - slen) > scale_FG - scale_fg + 30)
            {
                FGlen--;
            }
        }
// Move G to final position
        System.arraycopy(tmp, tmpOff + Gt, tmp, tmpOff + slen * n, slen * n);
        int Gt_final = slen * n;

// Final verification (skip if depth == 1)
        if (depth == 1)
        {
            return SOLVE_OK;
        }

        t2 = t1 + n;
        int t3 = t2 + n;
        int t4 = t3 + n;
        SmallPrime prime0 = PRIMES[0];
        int p = (int)prime0.p;
        int p0i = (int)prime0.p0i;
        int R2 = (int)prime0.R2;
        int Rx = mpRx31(slen, p, p0i, R2);

        mpMkgm(logn, tmp, tmpOff + t4, (int)prime0.g, p, p0i);

        if (use_sub_ntt)
        {
            t1 = ft;
            for (int u = 0; u < n; u++)
            {
                tmp[tmpOff + t2 + u] = zintModSmallSigned(tmp, tmpOff + Gt_final + u, slen, n, p, p0i, R2, Rx);
            }
            mpNTT(logn, tmp, t2, tmp, t4, p, p0i);
        }
        else
        {
            for (int u = 0; u < n; u++)
            {
                tmp[tmpOff + t1 + u] = zintModSmallSigned(tmp, tmpOff + ft + u, slen, n, p, p0i, R2, Rx);
                tmp[tmpOff + t2 + u] = zintModSmallSigned(tmp, tmpOff + Gt_final + u, slen, n, p, p0i, R2, Rx);
            }
            mpNTT(logn, tmp, tmpOff + t1, tmp, tmpOff + t4, p, p0i);
            mpNTT(logn, tmp, tmpOff + t2, tmp, tmpOff + t4, p, p0i);
        }

        for (int u = 0; u < n; u++)
        {
            tmp[tmpOff + t3 + u] = mpMontyMul(tmp[tmpOff + t1 + u], tmp[tmpOff + t2 + u], p, p0i);
        }

        if (use_sub_ntt)
        {
            t1 = gt;
            for (int u = 0; u < n; u++)
            {
                tmp[tmpOff + t2 + u] = zintModSmallSigned(tmp, tmpOff + Ft + u, slen, n, p, p0i, R2, Rx);
            }
            mpNTT(logn, tmp, tmpOff + t2, tmp, tmpOff + t4, p, p0i);
        }
        else
        {
            for (int u = 0; u < n; u++)
            {
                tmp[tmpOff + t1 + u] = zintModSmallSigned(tmp, tmpOff + gt + u, slen, n, p, p0i, R2, Rx);
                tmp[tmpOff + t2 + u] = zintModSmallSigned(tmp, tmpOff + Ft + u, slen, n, p, p0i, R2, Rx);
            }
            mpNTT(logn, tmp, tmpOff + t1, tmp, tmpOff + t4, p, p0i);
            mpNTT(logn, tmp, tmpOff + t2, tmp, tmpOff + t4, p, p0i);
        }

        int rv = mpMontyMul(prof.q, 1, p, p0i);
        for (int u = 0; u < n; u++)
        {
            int x = mpMontyMul(tmp[tmpOff + t1 + u], tmp[tmpOff + t2 + u], p, p0i);
            if (mpSub(tmp[tmpOff + t3 + u], x, p) != rv)
            {
                return SOLVE_ERR_REDUCE;
            }
        }

        return SOLVE_OK;
    }

    public static void mpMkgmigm(int logn, int[] gm, int gmOffset, int[] igm, int igmOffset,
                                 int g, int ig, int p, int p0i)
    {
        int n = 1 << logn;
        for (int j = logn; j < 10; j++)
        {
            g = mpMontyMul(g, g, p, p0i);
            ig = mpMontyMul(ig, ig, p, p0i);
        }

        int k = 10 - logn;
        int x1 = mpR(p);
        int x2 = mpHR(p);
        int u = 0;
        while (true)
        {
            int v = REV10[u << k];
            gm[gmOffset + v] = x1;
            igm[igmOffset + v] = x2;
            u++;
            if (u >= n)
            {
                break;
            }
            x1 = mpMontyMul(x1, g, p, p0i);
            x2 = mpMontyMul(x2, ig, p, p0i);
        }
    }

    private static long inner_fxr_div(long x, long y)
    {
        long sx = x >>> 63;
        x = (x ^ -sx) + sx;
        long sy = y >>> 63;
        y = (y ^ -sy) + sy;

        long q = 0;
        long num = x >>> 31;
        for (int i = 63; i >= 33; i--)
        {
            long b = 1 - ((num - y) >>> 63);
            q |= b << i;
            num -= y & -b;
            num <<= 1;
            num |= (x >>> (i - 33)) & 1;
        }
        for (int i = 32; i >= 0; i--)
        {
            long b = 1 - ((num - y) >>> 63);
            q |= b << i;
            num -= y & -b;
            num <<= 1;
        }

        long b = 1 - ((num - y) >>> 63);
        q += b;

        sx ^= sy;
        q = (q ^ -sx) + sx;
        return q;
    }

    // Convert big integer polynomial to fixed-point representation
    public static void poly_big_to_fixed(int logn, long[] d, int dOff, int[] f, int fOff, int len, int sc)
    {
        int n = 1 << logn;

        if (len == 0)
        {
            Arrays.fill(d, dOff, dOff + n, 0);
            return;
        }

        /*
         * We split the bit length into sch and scl such that:
         *   sc = 31*sch + scl
         * We also want scl in the 1..31 range, not 0..30. It may happen
         * that sch becomes -1, which will "wrap around" (harmlessly).
         */
        int sch, scl;
        int[] tmp = new int[2];
        DIVREM31(tmp, sc);
        sch = tmp[0];
        scl = tmp[1];

        // Adjust scl to be in range [1, 31]
        int z = (scl - 1) >>> 31;  // -1 if scl==0, 0 otherwise
        sch -= z;
        scl |= 31 & -z;  // If scl was 0, set it to 31

        int t0 = (sch - 1) & 0xFFFFFF;
        int t1 = sch & 0xFFFFFF;
        int t2 = (sch + 1) & 0xFFFFFF;

        for (int u = 0; u < n; u++)
        {
            int w0 = 0, w1 = 0, w2 = 0;

            // Extract words at positions t0, t1, t2
            for (int v = 0; v < len; v++)
            {
                int w = f[fOff + u + (v << logn)];  // f[v * n]

                // Constant-time selection based on word index
                int mask0 = ((v ^ t0) - 1) >>> 31;  // -1 if v == t0, 0 otherwise
                int mask1 = ((v ^ t1) - 1) >>> 31;
                int mask2 = ((v ^ t2) - 1) >>> 31;

                w0 |= w & -mask0;
                w1 |= w & -mask1;
                w2 |= w & -mask2;
            }

            /*
             * If there were not enough words for the requested
             * scaling, then we must supply copies with the proper
             * sign.
             */
            int lastWordIndex = (len - 1) << logn;
            int lastWord = f[fOff + u + lastWordIndex];
            int ws = -(lastWord >>> 30) >>> 1;  // Sign extension word

            // Conditionally apply sign extension
            int mask0 = -((len - sch) >>> 31);      // -1 if len < sch
            int mask1 = -((len - sch - 1) >>> 31);  // -1 if len < sch+1
            int mask2 = -((len - sch - 2) >>> 31);  // -1 if len < sch+2

            w0 |= ws & mask0;
            w1 |= ws & mask1;
            w2 |= ws & mask2;

            /*
             * Sign-extend w2 to ensure the sign bit is properly
             * set in the fixed-point value.
             */
            w2 |= (w2 & 0x40000000L) << 1;

            /*
             * Assemble the 64-bit value with the shifts.
             * Since the shift count (scl) is guaranteed to be in 1..31,
             * we do not have special cases to handle.
             */
            int xl = (w0 >>> (scl - 1)) | (w1 << (32 - scl));
            int xh = (w1 >>> scl) | (w2 << (31 - scl));

            // Combine into 64-bit fixed-point value
            long value = ((long)xl & 0xFFFFFFFFL) | (((long)xh & 0xFFFFFFFFL) << 32);
            d[dOff + u] = value;

            // Move to next coefficient
            // In the original code: f++ moves to next coefficient
            // Since we're using array indexing, we don't need to modify f
        }
    }

// Additional helper methods would be implemented similarly

    private static void DIVREM31(int[] result, int x)
    {
        int q = (int)((x * 67651L) >>> 21);
        int r = x - 31 * q;
        result[0] = q;
        result[1] = r;
    }

    public static int poly_max_bitlength(int logn, int[] f, int fOff, int flen)
    {
        if (flen == 0)
        {
            return 0;
        }

        int n = 1 << logn;
        int t = 0;
        int tk = 0;

        for (int u = 0; u < n; u++)
        {
            // Extend sign bit into a 31-bit mask
            int m = -(f[fOff + u + ((flen - 1) << logn)] >>> 30) & 0x7FFFFFFF;

            // Get top non-zero sign-adjusted word, with index
            int c = 0;
            int ck = 0;
            for (int v = 0; v < flen; v++)
            {
                int w = f[fOff + u + (v << logn)] ^ m; // sign-adjusted word
                int nz = ((w - 1) >>> 31) - 1; // 0xFFFFFFFF if w != 0, else 0
                c ^= nz & (c ^ w);
                ck ^= nz & (ck ^ v);
            }

            // If ck > tk, or tk == ck but c > t, then (c,ck) must
            // replace (t,tk) as current candidate for top word/index
            int rr = tbmask((tk - ck) | (((tk ^ ck) - 1) & (t - c)));
            t ^= rr & (t ^ c);
            tk ^= rr & (tk ^ ck);
        }

        // Get bit length of the top word (which has been sign-adjusted)
        // and return the result
        return 31 * tk + 32 - Integers.numberOfLeadingZeros(t);
    }

    static void fromLongArrayToByte32Array(int[] out, int outOff, long[] in)
    {
        for (int i = 0; i != in.length; i++)
        {
            out[outOff + 2 * i] = (int)in[i];
            out[outOff + 2 * i + 1] = (int)(in[i] >>> 32);
        }
    }

    static void fromByte32ArrayToShortArray(short[] out, int outOff, int[] in, int inOff, int len)
    {
        for (int i = 0; i != len; i++)
        {
            out[outOff + 2 * i] = (short)in[i + inOff];
            out[outOff + 2 * i + 1] = (short)(in[i + inOff] >>> 16);
        }
    }

    public static long fxr_sqr(long x)
    {
        long xl = x & 0xFFFFFFFFL;          // Extract low 32 bits (unsigned)
        int xh = (int)(x >>> 32);            // Extract high 32 bits (signed)

        long z0 = (xl * xl) >>> 32;         // Unsigned multiplication and shift
        long z1 = xl * (long)xh;            // Signed multiplication
        long z3 = ((long)xh * (long)xh) << 32;  // Signed multiplication and shift

        return z0 + (z1 << 1) + z3;         // Combine results
    }

    public static void vect_mul2e(int logn, long[] a, int aOff, int e)
    {
        int n = 1 << logn;
        for (int u = 0; u < n; u++)
        {
            a[u + aOff] = fxr_mul2e(a[u + aOff], e);
        }
    }

    private static long fxr_mul2e(long x, int n)
    {
        return x << n;
    }

    public static void vect_norm_fft(int logn, long[] d, int dPos, long[] a, int aPos, long[] b, int bPos)
    {
        int hn = 1 << (logn - 1);
        for (int u = 0; u < hn; u++)
        {
            // Calculate squared norms for complex numbers in a and b
            long aRealSqr = fxr_sqr(a[aPos + u]);
            long aImagSqr = fxr_sqr(a[aPos + u + hn]);
            long bRealSqr = fxr_sqr(b[bPos + u]);
            long bImagSqr = fxr_sqr(b[bPos + u + hn]);

            // Sum the squared components
            long aNorm = fxr_add(aRealSqr, aImagSqr);
            long bNorm = fxr_add(bRealSqr, bImagSqr);

            // Store the total norm
            d[dPos + u] = fxr_add(aNorm, bNorm);
        }
    }

    private static long fxr_add(long x, long y)
    {
        return x + y;
    }

    private static long fxr_div(long x, long y)
    {
        return inner_fxr_div(x, y);
    }

    private static long fxr_neg(long x)
    {
        return -x;
    }

    public static void vect_mul_fft(int logn, long[] a, int aPos, long[] b, int bPos)
    {
        int hn = 1 << (logn - 1);
        for (int u = 0; u < hn; u++)
        {
            long xRe = a[aPos + u];
            long xIm = a[aPos + u + hn];
            long yRe = b[bPos + u];
            long yIm = b[bPos + u + hn];
            long z0 = fxrMulInline(xRe, yRe);
            long z1 = fxrMulInline(xIm, yIm);
            long z2 = fxrMulInline(xRe + xIm, yRe + yIm);
            a[aPos + u]      = z0 - z1;
            a[aPos + u + hn] = z2 - (z0 + z1);
        }
    }

    public static void vect_add(int logn, long[] a, int aOff, long[] b, int bOff)
    {
        int n = 1 << logn;
        for (int u = 0; u < n; u++)
        {
            a[u + aOff] += b[u + bOff];
        }
    }

    public static int fxr_round(long x)
    {
        x += 0x80000000L;
        return (int)(x >>> 32);
    }

    public static void poly_sub_kfg_scaled_depth1(int logn_top, int[] F, int fOff, int[] G, int gOff, int FGlen,
                                                  int[] k, int kOff, int sc, byte[] f, byte[] g, int[] tmp, int tmpOff)
    {
        int logn = logn_top - 1;
        int n = 1 << logn;
        int hn = n >>> 1;
        int gm = tmpOff;
        int t1 = gm + n;
        int t2 = t1 + n;

        /*
         * Step 1: convert F and G to RNS. Since FGlen is equal to 1 or 2,
         * we do it with some specialized code. We assume that the RNS
         * representation does not lose information (i.e. each signed
         * coefficient is lower than (p0*p1)/2, with FGlen = 2 and the two
         * prime moduli are p0 and p1).
         */
        if (FGlen == 1)
        {
            int p = (int)PRIMES[0].p;
            for (int u = 0; u < n; u++)
            {
                int xf = F[u + fOff];
                int xg = G[u + gOff];
                xf |= (xf & 0x40000000) << 1;
                xg |= (xg & 0x40000000) << 1;
                F[u + fOff] = mp_set(xf, p);
                G[u + gOff] = mp_set(xg, p);
            }
        }
        else
        {
            int p0 = (int)PRIMES[0].p;
            int p0_0i = (int)PRIMES[0].p0i;
            int z0 = mpHalf((int)PRIMES[0].R2, p0);
            int p1 = (int)PRIMES[1].p;
            int p1_0i = (int)PRIMES[1].p0i;
            int z1 = mpHalf((int)PRIMES[1].R2, p1);

            for (int u = 0; u < n; u++)
            {
                int xl = F[u + fOff];
                int xh = F[u + fOff + n] | ((F[u + fOff + n] & 0x40000000) << 1);
                int yl0 = xl - (p0 & ~tbmask(xl - p0));
                int yh0 = mp_set(xh, p0);
                int r0 = mpAdd(yl0, mpMontyMul(yh0, z0, p0, p0_0i), p0);
                int yl1 = xl - (p1 & ~tbmask(xl - p1));
                int yh1 = mp_set(xh, p1);
                int r1 = mpAdd(yl1, mpMontyMul(yh1, z1, p1, p1_0i), p1);
                F[u + fOff] = r0;
                F[u + fOff + n] = r1;

                xl = G[u + gOff];
                xh = G[u + gOff + n] | ((G[u + gOff + n] & 0x40000000) << 1);
                yl0 = xl - (p0 & ~tbmask(xl - p0));
                yh0 = mp_set(xh, p0);
                r0 = mpAdd(yl0, mpMontyMul(yh0, z0, p0, p0_0i), p0);
                yl1 = xl - (p1 & ~tbmask(xl - p1));
                yh1 = mp_set(xh, p1);
                r1 = mpAdd(yl1, mpMontyMul(yh1, z1, p1, p1_0i), p1);
                G[u + gOff] = r0;
                G[u + gOff + n] = r1;
            }
        }

        /*
         * Step 2: for FGlen small primes, convert F and G to RNS+NTT,
         * and subtract (2^sc)*(ft,gt). The (ft,gt) polynomials are computed
         * in RNS+NTT dynamically.
         */
        for (int u = 0; u < FGlen; u++)
        {
            SmallPrime prime = PRIMES[u];
            int p = (int)prime.p;
            int p0i = (int)prime.p0i;
            int R2 = (int)prime.R2;
            int R3 = mpMontyMul(R2, R2, p, p0i);
            mpMkgm(logn, tmp, gm, (int)prime.g, p, p0i);

            /*
             * k <- (2^sc)*k (and into NTT).
             */
            int scv = mpMontyMul(1 << (sc & 31), R2, p, p0i);
            for (int m = sc >>> 5; m > 0; m--)
            {
                scv = mpMontyMul(scv, R2, p, p0i);
            }
            for (int v = 0; v < n; v++)
            {
                int x = mp_set(k[v + kOff], p);
                k[v + kOff] = mpMontyMul(scv, x, p, p0i);
            }
            mpNTT(logn, k, kOff, tmp, gm, p, p0i);

            /*
             * Convert F and G to NTT.
             */
            int Fu = fOff + (u << logn);
            int Gu = gOff + (u << logn);
            mpNTT(logn, F, Fu, tmp, gm, p, p0i);
            mpNTT(logn, G, Gu, tmp, gm, p, p0i);

            /*
             * Given the top-level f, we obtain ft = N(f) (the f at
             * depth 1) with:
             *    f = f_e(X^2) + X*f_o(X^2)
             * with f_e and f_o being modulo X^n+1. Then:
             *    N(f) = f_e^2 - X*f_o^2
             * The NTT representation of X is obtained from the gm[] tab:
             *    NTT(X)[2*j + 0] = gm[j + n/2]
             *    NTT(X)[2*j + 1] = -NTT(X)[2*j + 0]
             * Note that the values in gm[] are in Montgomery
             * representation.
             */
            for (int v = 0; v < n; v++)
            {
                tmp[t1 + v] = mp_set(f[(v << 1)], p);
                tmp[t2 + v] = mp_set(f[(v << 1) + 1], p);
            }
            mpNTT(logn, tmp, t1, tmp, gm, p, p0i);
            mpNTT(logn, tmp, t2, tmp, gm, p, p0i);

            for (int v = 0; v < hn; v++)
            {
                int xe0 = tmp[t1 + (v << 1)];
                int xe1 = tmp[t1 + (v << 1) + 1];
                int xo0 = tmp[t2 + (v << 1)];
                int xo1 = tmp[t2 + (v << 1) + 1];
                int xv0 = tmp[gm + hn + v];
                int xv1 = p - xv0;
                xe0 = mpMontyMul(xe0, xe0, p, p0i);
                xe1 = mpMontyMul(xe1, xe1, p, p0i);
                xo0 = mpMontyMul(xo0, xo0, p, p0i);
                xo1 = mpMontyMul(xo1, xo1, p, p0i);
                int xf0 = mpSub(xe0, mpMontyMul(xo0, xv0, p, p0i), p);
                int xf1 = mpSub(xe1, mpMontyMul(xo1, xv1, p, p0i), p);

                int xkf0 = mpMontyMul(mpMontyMul(xf0, k[kOff + (v << 1)], p, p0i), R3, p, p0i);
                int xkf1 = mpMontyMul(mpMontyMul(xf1, k[kOff + (v << 1) + 1], p, p0i), R3, p, p0i);
                F[Fu + (v << 1)] = mpSub(F[Fu + (v << 1)], xkf0, p);
                F[Fu + (v << 1) + 1] = mpSub(F[Fu + (v << 1) + 1], xkf1, p);
            }

            /*
             * Same treatment for G and gt.
             */
            for (int v = 0; v < n; v++)
            {
                tmp[t1 + v] = mp_set(g[(v << 1)], p);
                tmp[t2 + v] = mp_set(g[(v << 1) + 1], p);
            }
            mpNTT(logn, tmp, t1, tmp, gm, p, p0i);
            mpNTT(logn, tmp, t2, tmp, gm, p, p0i);

            for (int v = 0; v < hn; v++)
            {
                int xe0 = tmp[t1 + (v << 1)];
                int xe1 = tmp[t1 + (v << 1) + 1];
                int xo0 = tmp[t2 + (v << 1)];
                int xo1 = tmp[t2 + (v << 1) + 1];
                int xv0 = tmp[gm + hn + v];
                int xv1 = p - xv0;
                xe0 = mpMontyMul(xe0, xe0, p, p0i);
                xe1 = mpMontyMul(xe1, xe1, p, p0i);
                xo0 = mpMontyMul(xo0, xo0, p, p0i);
                xo1 = mpMontyMul(xo1, xo1, p, p0i);
                int xg0 = mpSub(xe0, mpMontyMul(xo0, xv0, p, p0i), p);
                int xg1 = mpSub(xe1, mpMontyMul(xo1, xv1, p, p0i), p);

                int xkg0 = mpMontyMul(mpMontyMul(xg0, k[kOff + (v << 1)], p, p0i), R3, p, p0i);
                int xkg1 = mpMontyMul(mpMontyMul(xg1, k[kOff + (v << 1) + 1], p, p0i), R3, p, p0i);
                G[Gu + (v << 1)] = mpSub(G[Gu + (v << 1)], xkg0, p);
                G[Gu + (v << 1) + 1] = mpSub(G[Gu + (v << 1) + 1], xkg1, p);
            }

            /*
             * Convert back F and G to RNS.
             */
            mpMkigm(logn, tmp, t1, (int)prime.ig, p, p0i);
            mpINTT(logn, F, Fu, tmp, t1, p, p0i);
            mpINTT(logn, G, Gu, tmp, t1, p, p0i);

            /*
             * We replaced k (plain 32-bit) with (2^sc)*k (NTT). We must
             * put it back to its initial value if there should be another
             * iteration.
             */
            if ((u + 1) < FGlen)
            {
                mpINTT(logn, k, kOff, tmp, t1, p, p0i);
                scv = 1 << (-sc & 31);
                for (int m = sc >>> 5; m > 0; m--)
                {
                    scv = mpMontyMul(scv, 1, p, p0i);
                }
                for (int v = 0; v < n; v++)
                {
                    k[v + kOff] = mpNorm(mpMontyMul(scv, k[v + kOff], p, p0i), p);
                }
            }
        }

        /*
         * Output F and G are in RNS (non-NTT), but we want plain integers.
         */
        if (FGlen == 1)
        {
            int p = (int)PRIMES[0].p;
            for (int u = 0; u < n; u++)
            {
                F[u + fOff] = mpNorm(F[u + fOff], p) & 0x7FFFFFFF;
                G[u + gOff] = mpNorm(G[u + gOff], p) & 0x7FFFFFFF;
            }
        }
        else
        {
            int p0 = (int)PRIMES[0].p;
            int p1 = (int)PRIMES[1].p;
            int p1_0i = (int)PRIMES[1].p0i;
            int s = (int)PRIMES[1].s;
            long pp = (long)p0 * (long)p1;
            long hpp = pp >>> 1;

            for (int u = 0; u < n; u++)
            {
                /*
                 * Apply CRT with two primes on the coefficient of F.
                 */
                int x0 = F[u + fOff];
                int x1 = F[u + fOff + n];
                int x0m1 = x0 - (p1 & ~tbmask(x0 - p1));
                int y = mpMontyMul(mpSub(x1, x0m1, p1), s, p1, p1_0i);
                long z = (long)x0 + (long)p0 * (long)y;
                z -= pp & -((hpp - z) >>> 63);
                F[u + fOff] = (int)z & 0x7FFFFFFF;
                F[u + fOff + n] = (int)(z >>> 31) & 0x7FFFFFFF;
            }

            for (int u = 0; u < n; u++)
            {
                /*
                 * Apply CRT with two primes on the coefficient of G.
                 */
                int x0 = G[u + gOff];
                int x1 = G[u + gOff + n];
                int x0m1 = x0 - (p1 & ~tbmask(x0 - p1));
                int y = mpMontyMul(mpSub(x1, x0m1, p1), s, p1, p1_0i);
                long z = (long)x0 + (long)p0 * (long)y;
                z -= pp & -((hpp - z) >>> 63);
                G[u + gOff] = (int)z & 0x7FFFFFFF;
                G[u + gOff + n] = (int)(z >>> 31) & 0x7FFFFFFF;
            }
        }
    }

    public static void poly_sub_scaled_ntt(int logn, int[] F, int FOff, int Flen,
                                           int[] f, int fOff, int flen,
                                           int[] k, int kOff, int sc, int[] tmp, int tmpOff)
    {
        int n = 1 << logn;
        int tlen = flen + 1;
        int gm = tmpOff;
        int igm = gm + n;
        int fk = igm + n;
        int t1 = fk + (tlen << logn);

        // DIVREM31 implementation
        int sch = (int)(sc / 31);
        int scl = sc % 31;

        // Compute k*f in fk[], in RNS notation
        for (int u = 0; u < tlen; u++)
        {
            SmallPrime prime = PRIMES[u];
            int p = (int)prime.p;
            int p0i = (int)prime.p0i;
            int R2 = (int)prime.R2;

            mpMkgmigm(logn, tmp, gm, tmp, igm, (int)prime.g, (int)prime.ig, p, p0i);

            for (int v = 0; v < n; v++)
            {
                tmp[t1 + v] = mp_set(k[kOff + v], p);
            }
            mpNTT(logn, tmp, t1, tmp, gm, p, p0i);

            int fs = fOff + (u << logn);
            int ff = fk + (u << logn);

            for (int v = 0; v < n; v++)
            {
                int term = mpMontyMul(tmp[t1 + v], f[fs + v], p, p0i);
                tmp[ff + v] = mpMontyMul(term, R2, p, p0i);
            }
            mpINTT(logn, tmp, ff, tmp, igm, p, p0i);
        }

        // Rebuild k*f
        zintRebuildCRT(tmp, fk, tlen, n, 1, 1, tmp, t1, PRIMES);

        // Subtract k*f, scaled, from F
        for (int u = 0; u < n; u++)
        {
            zint_sub_scaled(F, FOff + u, Flen, tmp, fk + u, tlen, n, sch, scl);
        }
    }

    public static void poly_sub_scaled(int logn,
                                       int[] F, int F_offset, int Flen,
                                       int[] f, int f_offset, int flen,
                                       int[] k, int k_offset, int sc)
    {
        if (flen == 0)
        {
            return;
        }
        int[] divRem = divRem31(sc);
        int sch = divRem[0];
        int scl = divRem[1];
        if (sch >= Flen)
        {
            return;
        }
        F_offset += sch << logn;
        Flen -= sch;
        switch (logn)
        {
        case 1:
        {
            int t0 = 0;
            int t1 = 0;
            int signf0 = -(f[f_offset + (flen << 1) - 2] >>> 30) >>> 1;
            int signf1 = -(f[f_offset + (flen << 1) - 1] >>> 30) >>> 1;
            int k0 = k[k_offset + 0];
            int k1 = k[k_offset + 1];
            long cc0 = 0;
            long cc1 = 0;
            for (int u = 0; u < Flen; u++)
            {
                int f0, f1;
                if (u < flen)
                {
                    f0 = f[f_offset + (u << 1) + 0];
                    f1 = f[f_offset + (u << 1) + 1];
                }
                else
                {
                    f0 = signf0;
                    f1 = signf1;
                }
                int fs0 = ((f0 << scl) & 0x7FFFFFFF) | t0;
                int fs1 = ((f1 << scl) & 0x7FFFFFFF) | t1;
                t0 = f0 >> (31 - scl);
                t1 = f1 >> (31 - scl);

                int F0 = F[F_offset + (u << 1) + 0];
                int F1 = F[F_offset + (u << 1) + 1];
                long z0 = (long)F0 + cc0
                    - (long)fs0 * (long)k0
                    + (long)fs1 * (long)k1;
                long z1 = (long)F1 + cc1
                    - (long)fs0 * (long)k1
                    - (long)fs1 * (long)k0;
                F[F_offset + (u << 1) + 0] = (int)z0 & 0x7FFFFFFF;
                F[F_offset + (u << 1) + 1] = (int)z1 & 0x7FFFFFFF;
                cc0 = z0 >> 31;
                cc1 = z1 >> 31;
            }
            return;
        }

        case 2:
        {
            int t0 = 0;
            int t1 = 0;
            int t2 = 0;
            int t3 = 0;
            int signf0 = -(f[f_offset + (flen << 2) - 4] >>> 30) >>> 1;
            int signf1 = -(f[f_offset + (flen << 2) - 3] >>> 30) >>> 1;
            int signf2 = -(f[f_offset + (flen << 2) - 2] >>> 30) >>> 1;
            int signf3 = -(f[f_offset + (flen << 2) - 1] >>> 30) >>> 1;
            int k0 = k[k_offset + 0];
            int k1 = k[k_offset + 1];
            int k2 = k[k_offset + 2];
            int k3 = k[k_offset + 3];
            long cc0 = 0;
            long cc1 = 0;
            long cc2 = 0;
            long cc3 = 0;
            for (int u = 0; u < Flen; u++)
            {
                int f0, f1, f2, f3;
                if (u < flen)
                {
                    f0 = f[f_offset + (u << 2) + 0];
                    f1 = f[f_offset + (u << 2) + 1];
                    f2 = f[f_offset + (u << 2) + 2];
                    f3 = f[f_offset + (u << 2) + 3];
                }
                else
                {
                    f0 = signf0;
                    f1 = signf1;
                    f2 = signf2;
                    f3 = signf3;
                }
                int fs0 = ((f0 << scl) & 0x7FFFFFFF) | t0;
                int fs1 = ((f1 << scl) & 0x7FFFFFFF) | t1;
                int fs2 = ((f2 << scl) & 0x7FFFFFFF) | t2;
                int fs3 = ((f3 << scl) & 0x7FFFFFFF) | t3;
                t0 = f0 >>> (31 - scl);
                t1 = f1 >>> (31 - scl);
                t2 = f2 >>> (31 - scl);
                t3 = f3 >>> (31 - scl);

                int F0 = F[F_offset + (u << 2) + 0];
                int F1 = F[F_offset + (u << 2) + 1];
                int F2 = F[F_offset + (u << 2) + 2];
                int F3 = F[F_offset + (u << 2) + 3];
                long z0 = (long)F0 + cc0
                    - (long)fs0 * (long)k0
                    + (long)fs1 * (long)k3
                    + (long)fs2 * (long)k2
                    + (long)fs3 * (long)k1;
                long z1 = (long)F1 + cc1
                    - (long)fs0 * (long)k1
                    - (long)fs1 * (long)k0
                    + (long)fs2 * (long)k3
                    + (long)fs3 * (long)k2;
                long z2 = (long)F2 + cc2
                    - (long)fs0 * (long)k2
                    - (long)fs1 * (long)k1
                    - (long)fs2 * (long)k0
                    + (long)fs3 * (long)k3;
                long z3 = (long)F3 + cc3
                    - (long)fs0 * (long)k3
                    - (long)fs1 * (long)k2
                    - (long)fs2 * (long)k1
                    - (long)fs3 * (long)k0;
                F[F_offset + (u << 2) + 0] = (int)z0 & 0x7FFFFFFF;
                F[F_offset + (u << 2) + 1] = (int)z1 & 0x7FFFFFFF;
                F[F_offset + (u << 2) + 2] = (int)z2 & 0x7FFFFFFF;
                F[F_offset + (u << 2) + 3] = (int)z3 & 0x7FFFFFFF;
                cc0 = z0 >> 31;
                cc1 = z1 >> 31;
                cc2 = z2 >> 31;
                cc3 = z3 >> 31;
            }
            return;
        }

        case 3:
        {
            int t0 = 0;
            int t1 = 0;
            int t2 = 0;
            int t3 = 0;
            int t4 = 0;
            int t5 = 0;
            int t6 = 0;
            int t7 = 0;
            int signf0 = -(f[f_offset + (flen << 3) - 8] >>> 30) >>> 1;
            int signf1 = -(f[f_offset + (flen << 3) - 7] >>> 30) >>> 1;
            int signf2 = -(f[f_offset + (flen << 3) - 6] >>> 30) >>> 1;
            int signf3 = -(f[f_offset + (flen << 3) - 5] >>> 30) >>> 1;
            int signf4 = -(f[f_offset + (flen << 3) - 4] >>> 30) >>> 1;
            int signf5 = -(f[f_offset + (flen << 3) - 3] >>> 30) >>> 1;
            int signf6 = -(f[f_offset + (flen << 3) - 2] >>> 30) >>> 1;
            int signf7 = -(f[f_offset + (flen << 3) - 1] >>> 30) >>> 1;
            int k0 = k[k_offset + 0];
            int k1 = k[k_offset + 1];
            int k2 = k[k_offset + 2];
            int k3 = k[k_offset + 3];
            int k4 = k[k_offset + 4];
            int k5 = k[k_offset + 5];
            int k6 = k[k_offset + 6];
            int k7 = k[k_offset + 7];
            long cc0 = 0;
            long cc1 = 0;
            long cc2 = 0;
            long cc3 = 0;
            long cc4 = 0;
            long cc5 = 0;
            long cc6 = 0;
            long cc7 = 0;
            for (int u = 0; u < Flen; u++)
            {
                int f0, f1, f2, f3, f4, f5, f6, f7;
                if (u < flen)
                {
                    f0 = f[f_offset + (u << 3) + 0];
                    f1 = f[f_offset + (u << 3) + 1];
                    f2 = f[f_offset + (u << 3) + 2];
                    f3 = f[f_offset + (u << 3) + 3];
                    f4 = f[f_offset + (u << 3) + 4];
                    f5 = f[f_offset + (u << 3) + 5];
                    f6 = f[f_offset + (u << 3) + 6];
                    f7 = f[f_offset + (u << 3) + 7];
                }
                else
                {
                    f0 = signf0;
                    f1 = signf1;
                    f2 = signf2;
                    f3 = signf3;
                    f4 = signf4;
                    f5 = signf5;
                    f6 = signf6;
                    f7 = signf7;
                }
                int fs0 = ((f0 << scl) & 0x7FFFFFFF) | t0;
                int fs1 = ((f1 << scl) & 0x7FFFFFFF) | t1;
                int fs2 = ((f2 << scl) & 0x7FFFFFFF) | t2;
                int fs3 = ((f3 << scl) & 0x7FFFFFFF) | t3;
                int fs4 = ((f4 << scl) & 0x7FFFFFFF) | t4;
                int fs5 = ((f5 << scl) & 0x7FFFFFFF) | t5;
                int fs6 = ((f6 << scl) & 0x7FFFFFFF) | t6;
                int fs7 = ((f7 << scl) & 0x7FFFFFFF) | t7;
                t0 = f0 >>> (31 - scl);
                t1 = f1 >>> (31 - scl);
                t2 = f2 >>> (31 - scl);
                t3 = f3 >>> (31 - scl);
                t4 = f4 >>> (31 - scl);
                t5 = f5 >>> (31 - scl);
                t6 = f6 >>> (31 - scl);
                t7 = f7 >>> (31 - scl);

                int F0 = F[F_offset + (u << 3) + 0];
                int F1 = F[F_offset + (u << 3) + 1];
                int F2 = F[F_offset + (u << 3) + 2];
                int F3 = F[F_offset + (u << 3) + 3];
                int F4 = F[F_offset + (u << 3) + 4];
                int F5 = F[F_offset + (u << 3) + 5];
                int F6 = F[F_offset + (u << 3) + 6];
                int F7 = F[F_offset + (u << 3) + 7];
                long z0 = (long)F0 + cc0
                    - (long)fs0 * (long)k0
                    + (long)fs1 * (long)k7
                    + (long)fs2 * (long)k6
                    + (long)fs3 * (long)k5
                    + (long)fs4 * (long)k4
                    + (long)fs5 * (long)k3
                    + (long)fs6 * (long)k2
                    + (long)fs7 * (long)k1;
                long z1 = (long)F1 + cc1
                    - (long)fs0 * (long)k1
                    - (long)fs1 * (long)k0
                    + (long)fs2 * (long)k7
                    + (long)fs3 * (long)k6
                    + (long)fs4 * (long)k5
                    + (long)fs5 * (long)k4
                    + (long)fs6 * (long)k3
                    + (long)fs7 * (long)k2;
                long z2 = (long)F2 + cc2
                    - (long)fs0 * (long)k2
                    - (long)fs1 * (long)k1
                    - (long)fs2 * (long)k0
                    + (long)fs3 * (long)k7
                    + (long)fs4 * (long)k6
                    + (long)fs5 * (long)k5
                    + (long)fs6 * (long)k4
                    + (long)fs7 * (long)k3;
                long z3 = (long)F3 + cc3
                    - (long)fs0 * (long)k3
                    - (long)fs1 * (long)k2
                    - (long)fs2 * (long)k1
                    - (long)fs3 * (long)k0
                    + (long)fs4 * (long)k7
                    + (long)fs5 * (long)k6
                    + (long)fs6 * (long)k5
                    + (long)fs7 * (long)k4;
                long z4 = (long)F4 + cc4
                    - (long)fs0 * (long)k4
                    - (long)fs1 * (long)k3
                    - (long)fs2 * (long)k2
                    - (long)fs3 * (long)k1
                    - (long)fs4 * (long)k0
                    + (long)fs5 * (long)k7
                    + (long)fs6 * (long)k6
                    + (long)fs7 * (long)k5;
                long z5 = (long)F5 + cc5
                    - (long)fs0 * (long)k5
                    - (long)fs1 * (long)k4
                    - (long)fs2 * (long)k3
                    - (long)fs3 * (long)k2
                    - (long)fs4 * (long)k1
                    - (long)fs5 * (long)k0
                    + (long)fs6 * (long)k7
                    + (long)fs7 * (long)k6;
                long z6 = (long)F6 + cc6
                    - (long)fs0 * (long)k6
                    - (long)fs1 * (long)k5
                    - (long)fs2 * (long)k4
                    - (long)fs3 * (long)k3
                    - (long)fs4 * (long)k2
                    - (long)fs5 * (long)k1
                    - (long)fs6 * (long)k0
                    + (long)fs7 * (long)k7;
                long z7 = (long)F7 + cc7
                    - (long)fs0 * (long)k7
                    - (long)fs1 * (long)k6
                    - (long)fs2 * (long)k5
                    - (long)fs3 * (long)k4
                    - (long)fs4 * (long)k3
                    - (long)fs5 * (long)k2
                    - (long)fs6 * (long)k1
                    - (long)fs7 * (long)k0;
                F[F_offset + (u << 3) + 0] = (int)z0 & 0x7FFFFFFF;
                F[F_offset + (u << 3) + 1] = (int)z1 & 0x7FFFFFFF;
                F[F_offset + (u << 3) + 2] = (int)z2 & 0x7FFFFFFF;
                F[F_offset + (u << 3) + 3] = (int)z3 & 0x7FFFFFFF;
                F[F_offset + (u << 3) + 4] = (int)z4 & 0x7FFFFFFF;
                F[F_offset + (u << 3) + 5] = (int)z5 & 0x7FFFFFFF;
                F[F_offset + (u << 3) + 6] = (int)z6 & 0x7FFFFFFF;
                F[F_offset + (u << 3) + 7] = (int)z7 & 0x7FFFFFFF;
                cc0 = z0 >> 31;
                cc1 = z1 >> 31;
                cc2 = z2 >> 31;
                cc3 = z3 >> 31;
                cc4 = z4 >> 31;
                cc5 = z5 >> 31;
                cc6 = z6 >> 31;
                cc7 = z7 >> 31;
            }
            return;
        }
        }
        int n = 1 << logn;
        for (int u = 0; u < n; u++)
        {
            int kf = -k[k_offset + u];
            int x_offset = F_offset + u;
            for (int v = 0; v < n; v++)
            {
                zint_add_scaled_mul_small(
                    F, x_offset, Flen,
                    f, f_offset + v, flen, n,
                    kf, 0, scl);
                if (u + v == n - 1)
                {
                    x_offset = F_offset;
                    kf = -kf;
                }
                else
                {
                    x_offset++;
                }
            }
        }
    }

    public static void zint_sub_scaled(int[] x, int xOffset, int xlen,
                                       int[] y, int yOffset, int ylen, int stride,
                                       int sch, int scl)
    {
        if (ylen == 0)
        {
            return;
        }

        int ysign = -(y[yOffset + stride * (ylen - 1)] >>> 30) >>> 1;
        int tw = 0;
        int cc = 0;
        xOffset += sch * stride;

        for (int u = sch; u < xlen; u++)
        {
            // Get the next word of (2^sc)*y
            int wy;
            if (ylen > 0)
            {
                wy = y[yOffset];
                yOffset += stride;
                ylen--;
            }
            else
            {
                wy = ysign;
            }
            int wys = ((wy << scl) & 0x7FFFFFFF) | tw;
            tw = wy >>> (31 - scl);

            int w = x[xOffset] - wys - cc;
            x[xOffset] = w & 0x7FFFFFFF;
            cc = w >>> 31;
            xOffset += stride;
        }
    }

    /**
     * Compute q = x / 31 and r = x % 31 for an unsigned integer x. This
     * method is constant-time and works for values x up to 63487 (inclusive).
     */
    public static int[] divRem31(int x)
    {
        int divrem31_x = x;
        int divrem31_q = (int)((divrem31_x * 67651L) >>> 21);
        int q = divrem31_q;
        int r = divrem31_x - 31 * divrem31_q;
        return new int[]{q, r};
    }

    public static void zint_add_scaled_mul_small(int[] x, int xOffset, int xlen,
                                                 int[] y, int yOffset, int ylen, int stride,
                                                 int k, int sch, int scl)
    {
        if (ylen == 0)
        {
            return;
        }

        int ysign = -(y[yOffset + stride * (ylen - 1)] >>> 30) >>> 1;
        int tw = 0;
        int cc = 0;
        xOffset += sch * stride;

        for (int u = sch; u < xlen; u++)
        {
            // Get the next word of (2^sc)*y
            int wy;
            if (ylen > 0)
            {
                wy = y[yOffset];
                yOffset += stride;
                ylen--;
            }
            else
            {
                wy = ysign;
            }
            int wys = ((wy << scl) & 0x7FFFFFFF) | tw;
            tw = wy >>> (31 - scl);

            // The expression below does not overflow due to 64-bit arithmetic
            long z = (long)wys * (long)k + (long)x[xOffset] + (long)cc;
            x[xOffset] = (int)z & 0x7FFFFFFF;
            xOffset += stride;

            // New carry word is a signed right-shift of z
            // In Java, we can directly cast the long to int for the signed shift
            cc = (int)(z >>> 31);
        }
    }

    /**
     * Solving the NTRU equation, top recursion level. This is a specialized
     * variant for solve_NTRU_intermediate() with depth == 0, for lower RAM
     * usage and faster operation.
     * <p>
     * Returned value: 0 on success, a negative error code otherwise.
     */
    public static int solve_NTRU_depth0(HawkParameters prof, int logn,
                                        byte[] f, int fOff, byte[] g, int gOff,
                                        int[] tmp, int tmpOffset)
    {
        int n = 1 << logn;
        int hn = n >>> 1;

        /*
         * At depth 0, all values fit on 30 bits, so we work with a
         * single modulus p.
         */
        SmallPrime prime = PRIMES[0];
        int p = (int)prime.p;
        int p0i = (int)prime.p0i;
        int R2 = (int)prime.R2;

        /*
         * Buffer layout:
         *   Fd   F from upper level (hn)
         *   Gd   G from upper level (hn)
         *   ft   f (n)
         *   gt   g (n)
         *   gm   helper for NTT
         */
        int Fd = tmpOffset;
        int Gd = Fd + hn;
        int ft = Gd + hn;
        int gt = ft + n;
        int gm = gt + n;

        /*
         * Load f and g, and convert to RNS+NTT.
         */
        mpMkgm(logn, tmp, gm, (int)prime.g, p, p0i);
        polyMpSetSmall(logn, tmp, ft, f, fOff, p);
        polyMpSetSmall(logn, tmp, gt, g, gOff, p);
        mpNTT(logn, tmp, ft, tmp, gm, p, p0i);
        mpNTT(logn, tmp, gt, tmp, gm, p, p0i);

        // Convert Fd and Gd to RNS+NTT
        polyMpSet(logn - 1, tmp, Fd, p);
        polyMpSet(logn - 1, tmp, Gd, p);
        mpNTT(logn - 1, tmp, Fd, tmp, gm, p, p0i);
        mpNTT(logn - 1, tmp, Gd, tmp, gm, p, p0i);

        // Build the unreduced (F,G) into ft and gt
        for (int u = 0; u < hn; u++)
        {
            int fa = tmp[ft + (u << 1) + 0];
            int fb = tmp[ft + (u << 1) + 1];
            int ga = tmp[gt + (u << 1) + 0];
            int gb = tmp[gt + (u << 1) + 1];
            int mFd = mpMontyMul(tmp[Fd + u], R2, p, p0i);
            int mGd = mpMontyMul(tmp[Gd + u], R2, p, p0i);
            tmp[ft + (u << 1) + 0] = mpMontyMul(gb, mFd, p, p0i);
            tmp[ft + (u << 1) + 1] = mpMontyMul(ga, mFd, p, p0i);
            tmp[gt + (u << 1) + 0] = mpMontyMul(fb, mGd, p, p0i);
            tmp[gt + (u << 1) + 1] = mpMontyMul(fa, mGd, p, p0i);
        }

        /*
         * Reorganize buffers:
         *   Fp   unreduced F (RNS+NTT) (n)
         *   Gp   unreduced G (RNS+NTT) (n)
         *   t1   free (n)
         *   t2   NTT support (gm) (n)
         *   t3   free (n)
         *   t4   free (n)
         */
        int Fp = tmpOffset;
        int Gp = Fp + n;
        int t1 = Gp + n;
        int t2 = t1 + n;  // alias on gm
        int t3 = t2 + n;
        int t4 = t3 + n;

        // Copy ft and gt to Fp and Gp
        System.arraycopy(tmp, ft, tmp, Fp, 2 * n);

        // Working modulo p (using the NTT), compute:
        // t1 <- F*adj(f) + G*adj(g)
        // t2 <- f*adj(f) + g*adj(g)

        // t4 <- f (RNS+NTT)
        polyMpSetSmall(logn, tmp, t4, f, fOff, p);
        mpNTT(logn, tmp, t4, tmp, gm, p, p0i);

        // t1 <- F*adj(f) (RNS+NTT)
        // t3 <- f*adj(f) (RNS+NTT)
        for (int u = 0; u < n; u++)
        {
            int w = mpMontyMul(tmp[t4 + (n - 1) - u], R2, p, p0i);
            tmp[t1 + u] = mpMontyMul(w, tmp[Fp + u], p, p0i);
            tmp[t3 + u] = mpMontyMul(w, tmp[t4 + u], p, p0i);
        }

        // t4 <- g (RNS+NTT)
        polyMpSetSmall(logn, tmp, t4, g, gOff, p);
        mpNTT(logn, tmp, t4, tmp, gm, p, p0i);

        // t1 <- t1 + G*adj(g)
        // t3 <- t3 + g*adj(g)
        for (int u = 0; u < n; u++)
        {
            int w = mpMontyMul(tmp[t4 + (n - 1) - u], R2, p, p0i);
            tmp[t1 + u] = mpAdd(tmp[t1 + u], mpMontyMul(w, tmp[Gp + u], p, p0i), p);
            tmp[t3 + u] = mpAdd(tmp[t3 + u], mpMontyMul(w, tmp[t4 + u], p, p0i), p);
        }

        /*
         * Convert back F*adj(f) + G*adj(g) and f*adj(f) + g*adj(g) to
         * plain representation, and move f*adj(f) + g*adj(g) to t2.
         */
        mpMkigm(logn, tmp, t4, (int)prime.ig, p, p0i);
        mpINTT(logn, tmp, t1, tmp, t4, p, p0i);
        mpINTT(logn, tmp, t3, tmp, t4, p, p0i);
        for (int u = 0; u < n; u++)
        {
            /*
             * NOTE: no truncature to 31 bits.
             */
            tmp[t1 + u] = mpNorm(tmp[t1 + u], p);
            tmp[t2 + u] = mpNorm(tmp[t3 + u], p);
        }

        /*
         * Buffer contents:
         *   Fp   unreduced F (RNS+NTT) (n)
         *   Gp   unreduced G (RNS+NTT) (n)
         *   t1   F*adj(f) + G*adj(g) (plain, 32-bit) (n)
         *   t2   f*adj(f) + g*adj(g) (plain, 32-bit) (n)
         */

        /*
         * We need to divide t1 by t2, and round the result. We convert
         * them to FFT representation, downscaled by 2^10 (to avoid overflows).
         * We first convert f*adj(f) + g*adj(g), which is auto-adjoint;
         * thus, its FFT representation only has half-size.
         */
        // Persistent fixed-point scratch (the C aliases this over tmp); the int
        // data and FFT scratch occupy disjoint regions, so no bulk sync is needed.
        long[] rt = new long[(tmp.length - tmpOffset) >>> 1];
        int rt3 = t3 >>> 1;
        for (int u = 0; u < n; u++)
        {
            rt[rt3 + u] = (long)tmp[t2 + u] << 22;
        }
        vectFFT(logn, rt, rt3);
        int rt2 = t2 >>> 1;

        System.arraycopy(rt, rt3, rt, rt2, hn);

        /*
         * Buffer contents:
         *   Fp    unreduced F (RNS+NTT) (n)
         *   Gp    unreduced G (RNS+NTT) (n)
         *   t1    F*adj(f) + G*adj(g) (plain, 32-bit) (n)
         *   rt2   f*adj(f) + g*adj(g) (FFT, auto-ajdoint) (hn fxr values = n)
         *   rt3   free (n fxr values = 2*n)
         */

        /*
         * Convert F*adj(f) + G*adj(g) to FFT (scaled by 2^10) (into rt3).
         */
        for (int u = 0; u < n; u++)
        {
            long x = (long)tmp[t1 + u] << 22;
            rt[rt3 + u] = x;
        }
        vectFFT(logn, rt, rt3);

        /*
         * Divide F*adj(f) + G*adj(g) by f*adj(f) + g*adj(g) and round
         * the result into t1, with conversion to RNS.
         */
        vectDivAutoAdjFFT(logn, rt, rt3, rt, rt2);
        vectIFFT(logn, rt, rt3);
        for (int u = 0; u < n; u++)
        {
            tmp[t1 + u] = mp_set(fxr_round(rt[rt3 + u]), p);
        }

        /*
         * Buffer contents:
         *   Fp    unreduced F (RNS+NTT) (n)
         *   Gp    unreduced G (RNS+NTT) (n)
         *   t1    k (RNS) (n)
         *   t2    free (n)
         *   t3    free (n)
         *   t4    free (n)
         */

        /*
         * Convert k to RNS+NTT+Montgomery.
         */
        mpMkgm(logn, tmp, t4, (int)prime.g, p, p0i);
        mpNTT(logn, tmp, t1, tmp, t4, p, p0i);
        for (int u = 0; u < n; u++)
        {
            tmp[t1 + u] = mpMontyMul(tmp[t1 + u], R2, p, p0i);
        }

        /*
         * Subtract k*f from F and k*g from G.
         * We also compute f*G - g*F (in RNS+NTT) to check that the solution
         * is correct.
         */
        for (int u = 0; u < n; u++)
        {
            tmp[t2 + u] = mp_set(f[fOff + u], p);
            tmp[t3 + u] = mp_set(g[gOff + u], p);
        }
        mpNTT(logn, tmp, t2, tmp, t4, p, p0i);
        mpNTT(logn, tmp, t3, tmp, t4, p, p0i);

        int rv = mpMontyMul(prof.q, 1, p, p0i);
        for (int u = 0; u < n; u++)
        {
            tmp[Fp + u] = mpSub(tmp[Fp + u], mpMontyMul(tmp[t1 + u], tmp[t2 + u], p, p0i), p);
            tmp[Gp + u] = mpSub(tmp[Gp + u], mpMontyMul(tmp[t1 + u], tmp[t3 + u], p, p0i), p);
            int x = mpSub(
                mpMontyMul(tmp[t2 + u], tmp[Gp + u], p, p0i),
                mpMontyMul(tmp[t3 + u], tmp[Fp + u], p, p0i), p);
            if (x != rv)
            {
                return SOLVE_ERR_REDUCE;
            }
        }

        // Convert back F and G into normal representation
        mpMkigm(logn, tmp, t4, (int)prime.ig, p, p0i);
        mpINTT(logn, tmp, Fp, tmp, t4, p, p0i);
        mpINTT(logn, tmp, Gp, tmp, t4, p, p0i);
        polyMpNorm(logn, tmp, Fp, p);
        polyMpNorm(logn, tmp, Gp, p);

        return SOLVE_OK;
    }

    public static void polyMpSet(int logn, int[] f, int fOffset, int p)
    {
        int n = 1 << logn;
        for (int u = 0; u < n; u++)
        {
            int x = f[fOffset + u];
            // Sign extension: if bit 30 is set, set bit 31 to make it negative in two's complement
            x |= (x & 0x40000000) << 1;
            // Convert to signed integer and then set modulo p
            f[fOffset + u] = mp_set(x, p);
        }
    }

    public static void vectDivAutoAdjFFT(int logn, long[] a, int aOffset, long[] b, int bOffset)
    {
        int hn = 1 << (logn - 1);
        for (int u = 0; u < hn; u++)
        {
            a[aOffset + u] = fxr_div(a[aOffset + u], b[bOffset + u]);
            a[aOffset + u + hn] = fxr_div(a[aOffset + u + hn], b[bOffset + u]);
        }
    }

    public static void polyMpNorm(int logn, int[] f, int fOffset, int p)
    {
        int n = 1 << logn;
        for (int u = 0; u < n; u++)
        {
            f[fOffset + u] = mpNorm(f[fOffset + u], p) & 0x7FFFFFFF;
        }
    }

    public static int solve_NTRU(HawkParameters prof, int logn,
                                 byte[] f, int fOff, byte[] g, int gOff, int[] tmp, int tmpOffset)
    {
        int n = 1 << logn;

        // Solve at the deepest level first
        int err = solve_NTRU_deepest(prof, logn, f, fOff, g, gOff, tmp, tmpOffset);
        if (err != SOLVE_OK)
        {
            return err;
        }

        // Solve intermediate levels recursively
        int depth = logn;
        while (depth-- > 1)
        {
            err = solve_NTRU_intermediate(prof, logn, f, fOff, g, gOff, depth, tmp, tmpOffset);
            if (err != SOLVE_OK)
            {
                return err;
            }
        }

        // Solve at depth 0
        err = solve_NTRU_depth0(prof, logn, f, fOff, g, gOff, tmp, tmpOffset);
        if (err != SOLVE_OK)
        {
            return err;
        }

        /*
         * F and G are at the start of tmp[] (plain, 31 bits per value).
         * We need to convert them to 8-bit representation, and check
         * that they are within the expected range.
         */

        // Allocate space for F and G as 8-bit values
        byte[] tmpBytes = Pack.intToLittleEndian(tmp);
        int F = n * 2 * 4;
        int G = F + n;

        int lim = prof.coeffFgLimit[logn];

        // Convert from 31-bit big integer representation to 8-bit small integers
        if (!polyBigToSmall(logn, tmpBytes, F, tmp, tmpOffset, lim))
        {
            return SOLVE_ERR_LIMIT;
        }
        if (!polyBigToSmall(logn, tmpBytes, G, tmp, tmpOffset + n, lim))
        {
            return SOLVE_ERR_LIMIT;
        }
        Pack.littleEndianToInt(tmpBytes, 0, tmp, 0, tmp.length);
        System.arraycopy(tmp, F >>> 2, tmp, tmpOffset, n >>> 1); // 2*n / 4
        return SOLVE_OK;
    }

    public static boolean polyBigToSmall(int logn, byte[] d, int dOffset,
                                         int[] s, int sOffset, int lim)
    {
        int n = 1 << logn;
        for (int u = 0; u < n; u++)
        {
            int x = s[sOffset + u];

            // Sign extension: if bit 30 is set, set bit 31 to make it negative in two's complement
            x |= (x & 0x40000000) << 1;

            // Convert to signed integer (Java integers are already signed)
            int z = x;

            // Check if the value is within the allowed range [-lim, lim]
            if (z < -lim || z > lim)
            {
                return false;
            }

            // Convert to 8-bit (truncation is safe since we checked the range)
            d[dOffset + u] = (byte)z;
        }
        return true;
    }

    private static long dec64le(byte[] buf, int offset)
    {
        // Convert 8 little-endian bytes to a Java long (unsigned 64-bit value)
        long value = 0;
        for (int i = 0; i < 8; i++)
        {
            value |= ((long)(buf[offset + i] & 0xFF)) << (8 * i);
        }
        return value;
    }

    private static void processChunk(SHAKEDigest shake, byte[] f, byte[] g,
                                     int j, int uLimit, int fSize,
                                     int shift, int chunkSize, int subtract)
    {
        byte[] qb = new byte[8];
        byte[] vv = new byte[chunkSize];

        for (int u = 0; u < uLimit; u += chunkSize * 4)
        { // Note: adjusted step
            shake.doOutput(qb, 0, 8);
            long q = dec64le(qb, 0);

            // Bit manipulation to convert random bytes to coefficients
            q = (q & 0x5555555555555555L) + ((q >>> 1) & 0x5555555555555555L);
            q = (q & 0x3333333333333333L) + ((q >>> 2) & 0x3333333333333333L);

            // Additional steps for higher dimensions
            if (chunkSize == 4)
            {
                q = (q & 0x0F0F0F0F0F0F0F0FL) + ((q >>> 4) & 0x0F0F0F0F0F0F0F0FL);
                q = (q & 0x00FF00FF00FF00FFL) + ((q >>> 8) & 0x00FF00FF00FF00FFL);
            }
            else if (chunkSize == 8)
            {
                q = (q & 0x0F0F0F0F0F0F0F0FL) + ((q >>> 4) & 0x0F0F0F0F0F0F0F0FL);
            }

            // Extract coefficients
            for (int i = 0; i < chunkSize; i++)
            {
                int val;
                if (chunkSize == 16)
                {
                    val = (int)(q & 0x0F) - subtract;
                    q >>= 4;
                }
                else if (chunkSize == 8)
                {
                    val = (int)(q & 0xFF) - subtract;
                    q >>= 8;
                }
                else
                { // chunkSize == 4
                    val = (int)(q & 0xFFFF) - subtract;
                    q >>= 16;
                }
                vv[i] = (byte)val;
            }

            // Copy to f or g array
            if (u < fSize)
            {
                System.arraycopy(vv, 0, f, u + (j << shift), chunkSize);
            }
            else
            {
                System.arraycopy(vv, 0, g, (u - fSize) + (j << shift), chunkSize);
            }
        }
    }

    public static void regen_fg_8(byte[] f, byte[] g, byte[] seed)
    {
        int seedLen = 16;
        for (int j = 0; j < 4; j++)
        {
            SHAKEDigest shake = new SHAKEDigest(256);
            shake.update(seed, 0, seedLen);
            shake.update((byte)j);
            // No explicit flip needed in Bouncy Castle - proceed directly to extraction
            processChunk(shake, f, g, j, 512, 256, 4, 16, 2);
        }
    }

    public static void regen_fg_9(byte[] f, byte[] g, byte[] seed)
    {
        int seedLen = 24;
        for (int j = 0; j < 4; j++)
        {
            SHAKEDigest shake = new SHAKEDigest(256);
            shake.update(seed, 0, seedLen);
            shake.update((byte)j);
            processChunk(shake, f, g, j, 1024, 512, 3, 8, 4);
        }
    }

    public static void regen_fg_10(byte[] f, byte[] g, byte[] seed)
    {
        int seedLen = 40;
        for (int j = 0; j < 4; j++)
        {
            SHAKEDigest shake = new SHAKEDigest(256);
            shake.update(seed, 0, seedLen);
            shake.update((byte)j);
            processChunk(shake, f, g, j, 2048, 1024, 2, 4, 8);
        }
    }

    public static void Hawk_regen_fg(int logn, byte[] f, int fOff, byte[] g, int gOff, byte[] seed)
    {
        switch (logn)
        {
        case 8:
            regen_fg_8(f, g, seed);
            break;
        case 9:
            regen_fg_9(f, g, seed);
            break;
        case 10:
            regen_fg_10(f, g, seed);
            break;
        default:
            throw new IllegalArgumentException("Unsupported logn: " + logn);
        }
    }

    public static int parity(int logn, byte[] f, int fOffset)
    {
        int n = 1 << logn;
        int pp = 0;
        for (int u = 0; u < n; u++)
        {
            // Treat the byte as unsigned by using bitwise AND with 0xFF
            pp += f[fOffset + u] & 0xFF;
        }
        return pp & 1;
    }

    public static boolean make_q001(int logn, int lim00, int lim01, int lim11,
                                    byte[] f, int fOffset, byte[] g, int gOffset,
                                    byte[] F, int FOffset, byte[] G, int GOffset,
                                    int[] tmp, int tmpOffset)
    {
        int n = 1 << logn;
        int hn = n >>> 1;

        SmallPrime prime = PRIMES[0];
        int p = (int)prime.p;
        int p0i = (int)prime.p0i;
        int R2 = (int)prime.R2;

        // Allocate temporary arrays within the tmp buffer
        int t1 = tmpOffset;
        int t2 = t1 + n;
        int t3 = t2 + n;
        int t4 = t3 + n;
        int t5 = t4 + n;

        // Generate NTT roots and convert polynomials to modular representation
        mpMkgm(logn, tmp, t1, (int)prime.g, p, p0i);
        polyMpSetSmall(logn, tmp, t2, f, fOffset, p);
        polyMpSetSmall(logn, tmp, t3, g, gOffset, p);
        polyMpSetSmall(logn, tmp, t4, F, FOffset, p);
        polyMpSetSmall(logn, tmp, t5, G, GOffset, p);

        // Apply NTT to all polynomials
        mpNTT(logn, tmp, t2, tmp, t1, p, p0i);
        mpNTT(logn, tmp, t3, tmp, t1, p, p0i);
        mpNTT(logn, tmp, t4, tmp, t1, p, p0i);
        mpNTT(logn, tmp, t5, tmp, t1, p, p0i);

        // Compute q00, q01, q11 in frequency domain
        for (int u = 0; u < hn; u++)
        {
            int xf = tmp[t2 + u];
            int xfa = tmp[t2 + (n - 1) - u];
            int xg = tmp[t3 + u];
            int xga = tmp[t3 + (n - 1) - u];
            int xF = tmp[t4 + u];
            int xFa = tmp[t4 + (n - 1) - u];
            int xG = tmp[t5 + u];
            int xGa = tmp[t5 + (n - 1) - u];

            // Compute q00 = f*adj(f) + g*adj(g)
            int term1 = mpMontyMul(xf, xfa, p, p0i);
            int term2 = mpMontyMul(xg, xga, p, p0i);
            int xq00 = mpMontyMul(R2, mpAdd(term1, term2, p), p, p0i);

            // Compute q11 = F*adj(F) + G*adj(G)
            term1 = mpMontyMul(xF, xFa, p, p0i);
            term2 = mpMontyMul(xG, xGa, p, p0i);
            int xq11 = mpMontyMul(R2, mpAdd(term1, term2, p), p, p0i);

            // Compute q01 = F*adj(f) + G*adj(g)
            int xq01_0 = mpMontyMul(R2, mpAdd(
                mpMontyMul(xF, xfa, p, p0i),
                mpMontyMul(xG, xga, p, p0i), p), p, p0i);

            // Compute adj(F)*f + adj(G)*g (symmetric to above)
            int xq01_1 = mpMontyMul(R2, mpAdd(
                mpMontyMul(xFa, xf, p, p0i),
                mpMontyMul(xGa, xg, p, p0i), p), p, p0i);

            // Check if q00 is invertible
            if (xq00 == 0)
            {
                return false;
            }

            // Store results (exploiting symmetry for auto-adjoint polynomials)
            tmp[t3 + u] = xq00;
            tmp[t3 + (n - 1) - u] = xq00;  // q00 is auto-adjoint
            tmp[t4 + u] = xq01_0;
            tmp[t4 + (n - 1) - u] = xq01_1; // q01 needs both parts
            tmp[t5 + u] = xq11;
            tmp[t5 + (n - 1) - u] = xq11;   // q11 is auto-adjoint
        }

        // Convert back to normal representation
        mpMkigm(logn, tmp, t1, (int)prime.ig, p, p0i);
        mpINTT(logn, tmp, t3, tmp, t1, p, p0i);
        mpINTT(logn, tmp, t4, tmp, t1, p, p0i);
        mpINTT(logn, tmp, t5, tmp, t1, p, p0i);

        // Store results in output arrays (q00, q01, q11).
        // q00, q01 are 16-bit, stored packed two-per-int starting at tmp[tmpOffset];
        // q11 is 32-bit, stored at tmp[tmpOffset + n]. These output regions are
        // disjoint from the source regions t3/t4/t5 (at tmpOffset + 2n..5n), so we
        // can pack each coefficient in place rather than round-tripping the whole
        // buffer through a short[] view per iteration (the C code aliases tmp).
        int q00ShortIdx = tmpOffset << 1;        // short index of q00[0]
        int q01ShortIdx = q00ShortIdx + n;       // short index of q01[0]
        int q11Offset = tmpOffset + n;           // int index of q11[0]
        // Check coefficient bounds and store results
        for (int u = 0; u < n; u++)
        {
            int xq00 = mpNorm(tmp[t3 + u], p);
            int xq01 = mpNorm(tmp[t4 + u], p);
            int xq11 = mpNorm(tmp[t5 + u], p);

            // Check coefficient bounds
            if (u == 0)
            {
                // Special bounds for constant term of q00
                if (xq00 < -32768 || xq00 > 32767)
                {
                    return false;
                }
            }
            else
            {
                // Regular bounds for other coefficients
                if (xq00 <= -lim00 || xq00 >= lim00)
                {
                    return false;
                }
                if (xq11 <= -lim11 || xq11 >= lim11)
                {
                    return false;
                }
            }

            // Check q01 bounds for all coefficients
            if (xq01 <= -lim01 || xq01 >= lim01)
            {
                return false;
            }

            // Store results (convert to appropriate types)
            setPackedShort(tmp, q00ShortIdx + u, (short)xq00);
            setPackedShort(tmp, q01ShortIdx + u, (short)xq01);
            tmp[q11Offset + u] = xq11;          // Store as int32_t
        }

        return true;
    }

    /*
     * Store a 16-bit value at logical short index 'shortIdx' into the int[] buffer,
     * using little-endian two-shorts-per-int packing: short 2k -> low half of int k,
     * short 2k+1 -> high half. This matches the short view the caller reads back.
     */
    private static void setPackedShort(int[] buf, int shortIdx, short val)
    {
        int i = shortIdx >>> 1;
        if ((shortIdx & 1) == 0)
        {
            buf[i] = (buf[i] & 0xFFFF0000) | (val & 0xFFFF);
        }
        else
        {
            buf[i] = (buf[i] & 0x0000FFFF) | (val << 16);
        }
    }
}
