text: #include "inner.h"
  
  /*
   * Falcon key pair generation.
   *
   * ==========================(LICENSE BEGIN)============================
   *
   * Copyright (c) 2017-2019  Falcon Project
   *
   * Permission is hereby granted, free of charge, to any person obtaining
   * a copy of this software and associated documentation files (the
   * "Software"), to deal in the Software without restriction, including
   * without limitation the rights to use, copy, modify, merge, publish,
   * distribute, sublicense, and/or sell copies of the Software, and to
   * permit persons to whom the Software is furnished to do so, subject to
   * the following conditions:
   *
   * The above copyright notice and this permission notice shall be
   * included in all copies or substantial portions of the Software.
   *
   * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
   * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
   * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
   * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
   * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
   * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
   * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
   *
   * ===========================(LICENSE END)=============================
   *
   * @author   Thomas Pornin <thomas.pornin@nccgroup.com>
   */
  
  
  #define MKN(logn)   ((size_t)1 << (logn))
  
  /* ==================================================================== */
  /*
   * Modular arithmetics.
   *
   * We implement a few functions for computing modulo a small integer p.
   *
   * All functions require that 2^30 < p < 2^31. Moreover, operands must
   * be in the 0..p-1 range.
   *
   * Modular addition and subtraction work for all such p.
   *
   * Montgomery multiplication requires that p is odd, and must be provided
   * with an additional value p0i = -1/p mod 2^31. See below for some basics
   * on Montgomery multiplication.
   *
   * Division computes an inverse modulo p by an exponentiation (with
   * exponent p-2): this works only if p is prime. Multiplication
   * requirements also apply, i.e. p must be odd and p0i must be provided.
   *
   * The NTT and inverse NTT need all of the above, and also that
   * p = 1 mod 2048.
   *
   * -----------------------------------------------------------------------
   *
   * We use Montgomery representation with 31-bit values:
   *
   *   Let R = 2^31 mod p. When 2^30 < p < 2^31, R = 2^31 - p.
   *   Montgomery representation of an integer x modulo p is x*R mod p.
   *
   *   Montgomery multiplication computes (x*y)/R mod p for
   *   operands x and y. Therefore:
   *
   *    - if operands are x*R and y*R (Montgomery representations of x and
   *      y), then Montgomery multiplication computes (x*R*y*R)/R = (x*y)*R
   *      mod p, which is the Montgomery representation of the product x*y;
   *
   *    - if operands are x*R and y (or x and y*R), then Montgomery
   *      multiplication returns x*y mod p: mixed-representation
   *      multiplications yield results in normal representation.
   *
   * To convert to Montgomery representation, we multiply by R, which is done
   * by Montgomery-multiplying by R^2. Stand-alone conversion back from
   * Montgomery representation is Montgomery-multiplication by 1.
   */
  
  /*
   * Precomputed small primes. Each element contains the following:
   *
   *  p   The prime itself.
   *
   *  g   A primitive root of phi = X^N+1 (in field Z_p).
   *
   *  s   The inverse of the product of all previous primes in the array,
   *      computed modulo p and in Montgomery representation.
   *
   * All primes are such that p = 1 mod 2048, and are lower than 2^31. They
   * are listed in decreasing order.
   */
  
  typedef struct {
      uint32_t p;
      uint32_t g;
      uint32_t s;
  } small_prime;
  
  static const small_prime PRIMES[] = {
      { 2147473409,  383167813,      10239 },
      { 2147389441,  211808905,  471403745 },
      { 2147387393,   37672282, 1329335065 },
      { 2147377153, 1977035326,  968223422 },
      { 2147358721, 1067163706,  132460015 },
      { 2147352577, 1606082042,  598693809 },
      { 2147346433, 2033915641, 1056257184 },
      { 2147338241, 1653770625,  421286710 },
      { 2147309569,  631200819, 1111201074 },
      { 2147297281, 2038364663, 1042003613 },
      { 2147295233, 1962540515,   19440033 },
      { 2147239937, 2100082663,  353296760 },
      { 2147235841, 1991153006, 1703918027 },
      { 2147217409,  516405114, 1258919613 },
      { 2147205121,  409347988, 1089726929 },
      { 2147196929,  927788991, 1946238668 },
      { 2147178497, 1136922411, 1347028164 },
      { 2147100673,  868626236,  701164723 },
      { 2147082241, 1897279176,  617820870 },
      { 2147074049, 1888819123,  158382189 },
      { 2147051521,   25006327,  522758543 },
      { 2147043329,  327546255,   37227845 },
      { 2147039233,  766324424, 1133356428 },
      { 2146988033, 1862817362,   73861329 },
      { 2146963457,  404622040,  653019435 },
      { 2146959361, 1936581214,  995143093 },
      { 2146938881, 1559770096,  634921513 },
      { 2146908161,  422623708, 1985060172 },
      { 2146885633, 1751189170,  298238186 },
      { 2146871297,  578919515,  291810829 },
      { 2146846721, 1114060353,  915902322 },
      { 2146834433, 2069565474,   47859524 },
      { 2146818049, 1552824584,  646281055 },
      { 2146775041, 1906267847, 1597832891 },
      { 2146756609, 1847414714, 1228090888 },
      { 2146744321, 1818792070, 1176377637 },
      { 2146738177, 1118066398, 1054971214 },
      { 2146736129,   52057278,  933422153 },
      { 2146713601,  592259376, 1406621510 },
      { 2146695169,  263161877, 1514178701 },
      { 2146656257,  685363115,  384505091 },
      { 2146650113,  927727032,  537575289 },
      { 2146646017,   52575506, 1799464037 },
      { 2146643969, 1276803876, 1348954416 },
      { 2146603009,  814028633, 1521547704 },
      { 2146572289, 1846678872, 1310832121 },
      { 2146547713,  919368090, 1019041349 },
      { 2146508801,  671847612,   38582496 },
      { 2146492417,  283911680,  532424562 },
      { 2146490369, 1780044827,  896447978 },
      { 2146459649,  327980850, 1327906900 },
      { 2146447361, 1310561493,  958645253 },
      { 2146441217,  412148926,  287271128 },
      { 2146437121,  293186449, 2009822534 },
      { 2146430977,  179034356, 1359155584 },
      { 2146418689, 1517345488, 1790248672 },
      { 2146406401, 1615820390, 1584833571 },
      { 2146404353,  826651445,  607120498 },
      { 2146379777,    3816988, 1897049071 },
      { 2146363393, 1221409784, 1986921567 },
      { 2146355201, 1388081168,  849968120 },
      { 2146336769, 1803473237, 1655544036 },
      { 2146312193, 1023484977,  273671831 },
      { 2146293761, 1074591448,  467406983 },
      { 2146283521,  831604668, 1523950494 },
      { 2146203649,  712865423, 1170834574 },
      { 2146154497, 1764991362, 1064856763 },
      { 2146142209,  627386213, 1406840151 },
      { 2146127873, 1638674429, 2088393537 },
      { 2146099201, 1516001018,  690673370 },
      { 2146093057, 1294931393,  315136610 },
      { 2146091009, 1942399533,  973539425 },
      { 2146078721, 1843461814, 2132275436 },
      { 2146060289, 1098740778,  360423481 },
      { 2146048001, 1617213232, 1951981294 },
      { 2146041857, 1805783169, 2075683489 },
      { 2146019329,  272027909, 1753219918 },
      { 2145986561, 1206530344, 2034028118 },
      { 2145976321, 1243769360, 1173377644 },
      { 2145964033,  887200839, 1281344586 },
      { 2145906689, 1651026455,  906178216 },
      { 2145875969, 1673238256, 1043521212 },
      { 2145871873, 1226591210, 1399796492 },
      { 2145841153, 1465353397, 1324527802 },
      { 2145832961, 1150638905,  554084759 },
      { 2145816577,  221601706,  427340863 },
      { 2145785857,  608896761,  316590738 },
      { 2145755137, 1712054942, 1684294304 },
      { 2145742849, 1302302867,  724873116 },
      { 2145728513,  516717693,  431671476 },
      { 2145699841,  524575579, 1619722537 },
      { 2145691649, 1925625239,  982974435 },
      { 2145687553,  463795662, 1293154300 },
      { 2145673217,  771716636,  881778029 },
      { 2145630209, 1509556977,  837364988 },
      { 2145595393,  229091856,  851648427 },
      { 2145587201, 1796903241,  635342424 },
      { 2145525761,  715310882, 1677228081 },
      { 2145495041, 1040930522,  200685896 },
      { 2145466369,  949804237, 1809146322 },
      { 2145445889, 1673903706,   95316881 },
      { 2145390593,  806941852, 1428671135 },
      { 2145372161, 1402525292,  159350694 },
      { 2145361921, 2124760298, 1589134749 },
      { 2145359873, 1217503067, 1561543010 },
      { 2145355777,  338341402,   83865711 },
      { 2145343489, 1381532164,  641430002 },
      { 2145325057, 1883895478, 1528469895 },
      { 2145318913, 1335370424,   65809740 },
      { 2145312769, 2000008042, 1919775760 },
      { 2145300481,  961450962, 1229540578 },
      { 2145282049,  910466767, 1964062701 },
      { 2145232897,  816527501,  450152063 },
      { 2145218561, 1435128058, 1794509700 },
      { 2145187841,   33505311, 1272467582 },
      { 2145181697,  269767433, 1380363849 },
      { 2145175553,   56386299, 1316870546 },
      { 2145079297, 2106880293, 1391797340 },
      { 2145021953, 1347906152,  720510798 },
      { 2145015809,  206769262, 1651459955 },
      { 2145003521, 1885513236, 1393381284 },
      { 2144960513, 1810381315,   31937275 },
      { 2144944129, 1306487838, 2019419520 },
      { 2144935937,   37304730, 1841489054 },
      { 2144894977, 1601434616,  157985831 },
      { 2144888833,   98749330, 2128592228 },
      { 2144880641, 1772327002, 2076128344 },
      { 2144864257, 1404514762, 2029969964 },
      { 2144827393,  801236594,  406627220 },
      { 2144806913,  349217443, 1501080290 },
      { 2144796673, 1542656776, 2084736519 },
      { 2144778241, 1210734884, 1746416203 },
      { 2144759809, 1146598851,  716464489 },
      { 2144757761,  286328400, 1823728177 },
      { 2144729089, 1347555695, 1836644881 },
      { 2144727041, 1795703790,  520296412 },
      { 2144696321, 1302475157,  852964281 },
      { 2144667649, 1075877614,  504992927 },
      { 2144573441,  198765808, 1617144982 },
      { 2144555009,  321528767,  155821259 },
      { 2144550913,  814139516, 1819937644 },
      { 2144536577,  571143206,  962942255 },
      { 2144524289, 1746733766,    2471321 },
      { 2144512001, 1821415077,  124190939 },
      { 2144468993,  917871546, 1260072806 },
      { 2144458753,  378417981, 1569240563 },
      { 2144421889,  175229668, 1825620763 },
      { 2144409601, 1699216963,  351648117 },
      { 2144370689, 1071885991,  958186029 },
      { 2144348161, 1763151227,  540353574 },
      { 2144335873, 1060214804,  919598847 },
      { 2144329729,  663515846, 1448552668 },
      { 2144327681, 1057776305,  590222840 },
      { 2144309249, 1705149168, 1459294624 },
      { 2144296961,  325823721, 1649016934 },
      { 2144290817,  738775789,  447427206 },
      { 2144243713,  962347618,  893050215 },
      { 2144237569, 1655257077,  900860862 },
      { 2144161793,  242206694, 1567868672 },
      { 2144155649,  769415308, 1247993134 },
      { 2144137217,  320492023,  515841070 },
      { 2144120833, 1639388522,  770877302 },
      { 2144071681, 1761785233,  964296120 },
      { 2144065537,  419817825,  204564472 },
      { 2144028673,  666050597, 2091019760 },
      { 2144010241, 1413657615, 1518702610 },
      { 2143952897, 1238327946,  475672271 },
      { 2143940609,  307063413, 1176750846 },
      { 2143918081, 2062905559,  786785803 },
      { 2143899649, 1338112849, 1562292083 },
      { 2143891457,   68149545,   87166451 },
      { 2143885313,  921750778,  394460854 },
      { 2143854593,  719766593,  133877196 },
      { 2143836161, 1149399850, 1861591875 },
      { 2143762433, 1848739366, 1335934145 },
      { 2143756289, 1326674710,  102999236 },
      { 2143713281,  808061791, 1156900308 },
      { 2143690753,  388399459, 1926468019 },
      { 2143670273, 1427891374, 1756689401 },
      { 2143666177, 1912173949,  986629565 },
      { 2143645697, 2041160111,  371842865 },
      { 2143641601, 1279906897, 2023974350 },
      { 2143635457,  720473174, 1389027526 },
      { 2143621121, 1298309455, 1732632006 },
      { 2143598593, 1548762216, 1825417506 },
      { 2143567873,  620475784, 1073787233 },
      { 2143561729, 1932954575,  949167309 },
      { 2143553537,  354315656, 1652037534 },
      { 2143541249,  577424288, 1097027618 },
      { 2143531009,  357862822,  478640055 },
      { 2143522817, 2017706025, 1550531668 },
      { 2143506433, 2078127419, 1824320165 },
      { 2143488001,  613475285, 1604011510 },
      { 2143469569, 1466594987,  502095196 },
      { 2143426561, 1115430331, 1044637111 },
      { 2143383553,    9778045, 1902463734 },
      { 2143377409, 1557401276, 2056861771 },
      { 2143363073,  652036455, 1965915971 },
      { 2143260673, 1464581171, 1523257541 },
      { 2143246337, 1876119649,  764541916 },
      { 2143209473, 1614992673, 1920672844 },
      { 2143203329,  981052047, 2049774209 },
      { 2143160321, 1847355533,  728535665 },
      { 2143129601,  965558457,  603052992 },
      { 2143123457, 2140817191,    8348679 },
      { 2143100929, 1547263683,  694209023 },
      { 2143092737,  643459066, 1979934533 },
      { 2143082497,  188603778, 2026175670 },
      { 2143062017, 1657329695,  377451099 },
      { 2143051777,  114967950,  979255473 },
      { 2143025153, 1698431342, 1449196896 },
      { 2143006721, 1862741675, 1739650365 },
      { 2142996481,  756660457,  996160050 },
      { 2142976001,  927864010, 1166847574 },
      { 2142965761,  905070557,  661974566 },
      { 2142916609,   40932754, 1787161127 },
      { 2142892033, 1987985648,  675335382 },
      { 2142885889,  797497211, 1323096997 },
      { 2142871553, 2068025830, 1411877159 },
      { 2142861313, 1217177090, 1438410687 },
      { 2142830593,  409906375, 1767860634 },
      { 2142803969, 1197788993,  359782919 },
      { 2142785537,  643817365,  513932862 },
      { 2142779393, 1717046338,  218943121 },
      { 2142724097,   89336830,  416687049 },
      { 2142707713,    5944581, 1356813523 },
      { 2142658561,  887942135, 2074011722 },
      { 2142638081,  151851972, 1647339939 },
      { 2142564353, 1691505537, 1483107336 },
      { 2142533633, 1989920200, 1135938817 },
      { 2142529537,  959263126, 1531961857 },
      { 2142527489,  453251129, 1725566162 },
      { 2142502913, 1536028102,  182053257 },
      { 2142498817,  570138730,  701443447 },
      { 2142416897,  326965800,  411931819 },
      { 2142363649, 1675665410, 1517191733 },
      { 2142351361,  968529566, 1575712703 },
      { 2142330881, 1384953238, 1769087884 },
      { 2142314497, 1977173242, 1833745524 },
      { 2142289921,   95082313, 1714775493 },
      { 2142283777,  109377615, 1070584533 },
      { 2142277633,   16960510,  702157145 },
      { 2142263297,  553850819,  431364395 },
      { 2142208001,  241466367, 2053967982 },
      { 2142164993, 1795661326, 1031836848 },
      { 2142097409, 1212530046,  712772031 },
      { 2142087169, 1763869720,  822276067 },
      { 2142078977,  644065713, 1765268066 },
      { 2142074881,  112671944,  643204925 },
      { 2142044161, 1387785471, 1297890174 },
      { 2142025729,  783885537, 1000425730 },
      { 2142011393,  905662232, 1679401033 },
      { 2141974529,  799788433,  468119557 },
      { 2141943809, 1932544124,  449305555 },
      { 2141933569, 1527403256,  841867925 },
      { 2141931521, 1247076451,  743823916 },
      { 2141902849, 1199660531,  401687910 },
      { 2141890561,  150132350, 1720336972 },
      { 2141857793, 1287438162,  663880489 },
      { 2141833217,  618017731, 1819208266 },
      { 2141820929,  999578638, 1403090096 },
      { 2141786113,   81834325, 1523542501 },
      { 2141771777,  120001928,  463556492 },
      { 2141759489,  122455485, 2124928282 },
      { 2141749249,  141986041,  940339153 },
      { 2141685761,  889088734,  477141499 },
      { 2141673473,  324212681, 1122558298 },
      { 2141669377, 1175806187, 1373818177 },
      { 2141655041, 1113654822,  296887082 },
      { 2141587457,  991103258, 1585913875 },
      { 2141583361, 1401451409, 1802457360 },
      { 2141575169, 1571977166,  712760980 },
      { 2141546497, 1107849376, 1250270109 },
      { 2141515777,  196544219,  356001130 },
      { 2141495297, 1733571506, 1060744866 },
      { 2141483009,  321552363, 1168297026 },
      { 2141458433,  505818251,  733225819 },
      { 2141360129, 1026840098,  948342276 },
      { 2141325313,  945133744, 2129965998 },
      { 2141317121, 1871100260, 1843844634 },
      { 2141286401, 1790639498, 1750465696 },
      { 2141267969, 1376858592,  186160720 },
      { 2141255681, 2129698296, 1876677959 },
      { 2141243393, 2138900688, 1340009628 },
      { 2141214721, 1933049835, 1087819477 },
      { 2141212673, 1898664939, 1786328049 },
      { 2141202433,  990234828,  940682169 },
      { 2141175809, 1406392421,  993089586 },
      { 2141165569, 1263518371,  289019479 },
      { 2141073409, 1485624211,  507864514 },
      { 2141052929, 1885134788,  311252465 },
      { 2141040641, 1285021247,  280941862 },
      { 2141028353, 1527610374,  375035110 },
      { 2141011969, 1400626168,  164696620 },
      { 2140999681,  632959608,  966175067 },
      { 2140997633, 2045628978, 1290889438 },
      { 2140993537, 1412755491,  375366253 },
      { 2140942337,  719477232,  785367828 },
      { 2140925953,   45224252,  836552317 },
      { 2140917761, 1157376588, 1001839569 },
      { 2140887041,  278480752, 2098732796 },
      { 2140837889, 1663139953,  924094810 },
      { 2140788737,  802501511, 2045368990 },
      { 2140766209, 1820083885, 1800295504 },
      { 2140764161, 1169561905, 2106792035 },
      { 2140696577,  127781498, 1885987531 },
      { 2140684289,   16014477, 1098116827 },
      { 2140653569,  665960598, 1796728247 },
      { 2140594177, 1043085491,  377310938 },
      { 2140579841, 1732838211, 1504505945 },
      { 2140569601,  302071939,  358291016 },
      { 2140567553,  192393733, 1909137143 },
      { 2140557313,  406595731, 1175330270 },
      { 2140549121, 1748850918,  525007007 },
      { 2140477441,  499436566, 1031159814 },
      { 2140469249, 1886004401, 1029951320 },
      { 2140426241, 1483168100, 1676273461 },
      { 2140420097, 1779917297,  846024476 },
      { 2140413953,  522948893, 1816354149 },
      { 2140383233, 1931364473, 1296921241 },
      { 2140366849, 1917356555,  147196204 },
      { 2140354561,   16466177, 1349052107 },
      { 2140348417, 1875366972, 1860485634 },
      { 2140323841,  456498717, 1790256483 },
      { 2140321793, 1629493973,  150031888 },
      { 2140315649, 1904063898,  395510935 },
      { 2140280833, 1784104328,  831417909 },
      { 2140250113,  256087139,  697349101 },
      { 2140229633,  388553070,  243875754 },
      { 2140223489,  747459608, 1396270850 },
      { 2140200961,  507423743, 1895572209 },
      { 2140162049,  580106016, 2045297469 },
      { 2140149761,  712426444,  785217995 },
      { 2140137473, 1441607584,  536866543 },
      { 2140119041,  346538902, 1740434653 },
      { 2140090369,  282642885,   21051094 },
      { 2140076033, 1407456228,  319910029 },
      { 2140047361, 1619330500, 1488632070 },
      { 2140041217, 2089408064, 2012026134 },
      { 2140008449, 1705524800, 1613440760 },
      { 2139924481, 1846208233, 1280649481 },
      { 2139906049,  989438755, 1185646076 },
      { 2139867137, 1522314850,  372783595 },
      { 2139842561, 1681587377,  216848235 },
      { 2139826177, 2066284988, 1784999464 },
      { 2139824129,  480888214, 1513323027 },
      { 2139789313,  847937200,  858192859 },
      { 2139783169, 1642000434, 1583261448 },
      { 2139770881,  940699589,  179702100 },
      { 2139768833,  315623242,  964612676 },
      { 2139666433,  331649203,  764666914 },
      { 2139641857, 2118730799, 1313764644 },
      { 2139635713,  519149027,  519212449 },
      { 2139598849, 1526413634, 1769667104 },
      { 2139574273,  551148610,  820739925 },
      { 2139568129, 1386800242,  472447405 },
      { 2139549697,  813760130, 1412328531 },
      { 2139537409, 1615286260, 1609362979 },
      { 2139475969, 1352559299, 1696720421 },
      { 2139455489, 1048691649, 1584935400 },
      { 2139432961,  836025845,  950121150 },
      { 2139424769, 1558281165, 1635486858 },
      { 2139406337, 1728402143, 1674423301 },
      { 2139396097, 1727715782, 1483470544 },
      { 2139383809, 1092853491, 1741699084 },
      { 2139369473,  690776899, 1242798709 },
      { 2139351041, 1768782380, 2120712049 },
      { 2139334657, 1739968247, 1427249225 },
      { 2139332609, 1547189119,  623011170 },
      { 2139310081, 1346827917, 1605466350 },
      { 2139303937,  369317948,  828392831 },
      { 2139301889, 1560417239, 1788073219 },
      { 2139283457, 1303121623,  595079358 },
      { 2139248641, 1354555286,  573424177 },
      { 2139240449,   60974056,  885781403 },
      { 2139222017,  355573421, 1221054839 },
      { 2139215873,  566477826, 1724006500 },
      { 2139150337,  871437673, 1609133294 },
      { 2139144193, 1478130914, 1137491905 },
      { 2139117569, 1854880922,  964728507 },
      { 2139076609,  202405335,  756508944 },
      { 2139062273, 1399715741,  884826059 },
      { 2139045889, 1051045798, 1202295476 },
      { 2139033601, 1707715206,  632234634 },
      { 2139006977, 2035853139,  231626690 },
      { 2138951681,  183867876,  838350879 },
      { 2138945537, 1403254661,  404460202 },
      { 2138920961,  310865011, 1282911681 },
      { 2138910721, 1328496553,  103472415 },
      { 2138904577,   78831681,  993513549 },
      { 2138902529, 1319697451, 1055904361 },
      { 2138816513,  384338872, 1706202469 },
      { 2138810369, 1084868275,  405677177 },
      { 2138787841,  401181788, 1964773901 },
      { 2138775553, 1850532988, 1247087473 },
      { 2138767361,  874261901, 1576073565 },
      { 2138757121, 1187474742,  993541415 },
      { 2138748929, 1782458888, 1043206483 },
      { 2138744833, 1221500487,  800141243 },
      { 2138738689,  413465368, 1450660558 },
      { 2138695681,  739045140,  342611472 },
      { 2138658817, 1355845756,  672674190 },
      { 2138644481,  608379162, 1538874380 },
      { 2138632193, 1444914034,  686911254 },
      { 2138607617,  484707818, 1435142134 },
      { 2138591233,  539460669, 1290458549 },
      { 2138572801, 2093538990, 2011138646 },
      { 2138552321, 1149786988, 1076414907 },
      { 2138546177,  840688206, 2108985273 },
      { 2138533889,  209669619,  198172413 },
      { 2138523649, 1975879426, 1277003968 },
      { 2138490881, 1351891144, 1976858109 },
      { 2138460161, 1817321013, 1979278293 },
      { 2138429441, 1950077177,  203441928 },
      { 2138400769,  908970113,  628395069 },
      { 2138398721,  219890864,  758486760 },
      { 2138376193, 1306654379,  977554090 },
      { 2138351617,  298822498, 2004708503 },
      { 2138337281,  441457816, 1049002108 },
      { 2138320897, 1517731724, 1442269609 },
      { 2138290177, 1355911197, 1647139103 },
      { 2138234881,  531313247, 1746591962 },
      { 2138214401, 1899410930,  781416444 },
      { 2138202113, 1813477173, 1622508515 },
      { 2138191873, 1086458299, 1025408615 },
      { 2138183681, 1998800427,  827063290 },
      { 2138173441, 1921308898,  749670117 },
      { 2138103809, 1620902804, 2126787647 },
      { 2138099713,  828647069, 1892961817 },
      { 2138085377,  179405355, 1525506535 },
      { 2138060801,  615683235, 1259580138 },
      { 2138044417, 2030277840, 1731266562 },
      { 2138042369, 2087222316, 1627902259 },
      { 2138032129,  126388712, 1108640984 },
      { 2138011649,  715026550, 1017980050 },
      { 2137993217, 1693714349, 1351778704 },
      { 2137888769, 1289762259, 1053090405 },
      { 2137853953,  199991890, 1254192789 },
      { 2137833473,  941421685,  896995556 },
      { 2137817089,  750416446, 1251031181 },
      { 2137792513,  798075119,  368077456 },
      { 2137786369,  878543495, 1035375025 },
      { 2137767937,    9351178, 1156563902 },
      { 2137755649, 1382297614, 1686559583 },
      { 2137724929, 1345472850, 1681096331 },
      { 2137704449,  834666929,  630551727 },
      { 2137673729, 1646165729, 1892091571 },
      { 2137620481,  778943821,   48456461 },
      { 2137618433, 1730837875, 1713336725 },
      { 2137581569,  805610339, 1378891359 },
      { 2137538561,  204342388, 1950165220 },
      { 2137526273, 1947629754, 1500789441 },
      { 2137516033,  719902645, 1499525372 },
      { 2137491457,  230451261,  556382829 },
      { 2137440257,  979573541,  412760291 },
      { 2137374721,  927841248, 1954137185 },
      { 2137362433, 1243778559,  861024672 },
      { 2137313281, 1341338501,  980638386 },
      { 2137311233,  937415182, 1793212117 },
      { 2137255937,  795331324, 1410253405 },
      { 2137243649,  150756339, 1966999887 },
      { 2137182209,  163346914, 1939301431 },
      { 2137171969, 1952552395,  758913141 },
      { 2137159681,  570788721,  218668666 },
      { 2137147393, 1896656810, 2045670345 },
      { 2137141249,  358493842,  518199643 },
      { 2137139201, 1505023029,  674695848 },
      { 2137133057,   27911103,  830956306 },
      { 2137122817,  439771337, 1555268614 },
      { 2137116673,  790988579, 1871449599 },
      { 2137110529,  432109234,  811805080 },
      { 2137102337, 1357900653, 1184997641 },
      { 2137098241,  515119035, 1715693095 },
      { 2137090049,  408575203, 2085660657 },
      { 2137085953, 2097793407, 1349626963 },
      { 2137055233, 1556739954, 1449960883 },
      { 2137030657, 1545758650, 1369303716 },
      { 2136987649,  332602570,  103875114 },
      { 2136969217, 1499989506, 1662964115 },
      { 2136924161,  857040753,    4738842 },
      { 2136895489, 1948872712,  570436091 },
      { 2136893441,   58969960, 1568349634 },
      { 2136887297, 2127193379,  273612548 },
      { 2136850433,  111208983, 1181257116 },
      { 2136809473, 1627275942, 1680317971 },
      { 2136764417, 1574888217,   14011331 },
      { 2136741889,   14011055, 1129154251 },
      { 2136727553,   35862563, 1838555253 },
      { 2136721409,  310235666, 1363928244 },
      { 2136698881, 1612429202, 1560383828 },
      { 2136649729, 1138540131,  800014364 },
      { 2136606721,  602323503, 1433096652 },
      { 2136563713,  182209265, 1919611038 },
      { 2136555521,  324156477,  165591039 },
      { 2136549377,  195513113,  217165345 },
      { 2136526849, 1050768046,  939647887 },
      { 2136508417, 1886286237, 1619926572 },
      { 2136477697,  609647664,   35065157 },
      { 2136471553,  679352216, 1452259468 },
      { 2136457217,  128630031,  824816521 },
      { 2136422401,   19787464, 1526049830 },
      { 2136420353,  698316836, 1530623527 },
      { 2136371201, 1651862373, 1804812805 },
      { 2136334337,  326596005,  336977082 },
      { 2136322049,   63253370, 1904972151 },
      { 2136297473,  312176076,  172182411 },
      { 2136248321,  381261841,  369032670 },
      { 2136242177,  358688773, 1640007994 },
      { 2136229889,  512677188,   75585225 },
      { 2136219649, 2095003250, 1970086149 },
      { 2136207361, 1909650722,  537760675 },
      { 2136176641, 1334616195, 1533487619 },
      { 2136158209, 2096285632, 1793285210 },
      { 2136143873, 1897347517,  293843959 },
      { 2136133633,  923586222, 1022655978 },
      { 2136096769, 1464868191, 1515074410 },
      { 2136094721, 2020679520, 2061636104 },
      { 2136076289,  290798503, 1814726809 },
      { 2136041473,  156415894, 1250757633 },
      { 2135996417,  297459940, 1132158924 },
      { 2135955457,  538755304, 1688831340 },
      { 0, 0, 0 }
  };
  
  /*
   * Reduce a small signed integer modulo a small prime. The source
   * value x MUST be such that -p < x < p.
   */
  static inline uint32_t
  modp_set(int32_t x, uint32_t p) {
      uint32_t w;
  
      w = (uint32_t)x;
      w += p & -(w >> 31);
      return w;
  }
  
  /*
   * Normalize a modular integer around 0.
   */
  static inline int32_t
  modp_norm(uint32_t x, uint32_t p) {
      return (int32_t)(x - (p & (((x - ((p + 1) >> 1)) >> 31) - 1)));
  }
  
  /*
   * Compute -1/p mod 2^31. This works for all odd integers p that fit
   * on 31 bits.
   */
  static uint32_t
  modp_ninv31(uint32_t p) {
      uint32_t y;
  
      y = 2 - p;
      y *= 2 - p * y;
      y *= 2 - p * y;
      y *= 2 - p * y;
      y *= 2 - p * y;
      return (uint32_t)0x7FFFFFFF & -y;
  }
  
  /*
   * Compute R = 2^31 mod p.
   */
  static inline uint32_t
  modp_R(uint32_t p) {
      /*
       * Since 2^30 < p < 2^31, we know that 2^31 mod p is simply
       * 2^31 - p.
       */
      return ((uint32_t)1 << 31) - p;
  }
  
  /*
   * Addition modulo p.
   */
  static inline uint32_t
  modp_add(uint32_t a, uint32_t b, uint32_t p) {
      uint32_t d;
  
      d = a + b - p;
      d += p & -(d >> 31);
      return d;
  }
  
  /*
   * Subtraction modulo p.
   */
  static inline uint32_t
  modp_sub(uint32_t a, uint32_t b, uint32_t p) {
      uint32_t d;
  
      d = a - b;
      d += p & -(d >> 31);
      return d;
  }
  
  /*
   * Halving modulo p.
   */
  /* unused
  static inline uint32_t
  modp_half(uint32_t a, uint32_t p)
  {
      a += p & -(a & 1);
      return a >> 1;
  }
  */
  
  /*
   * Montgomery multiplication modulo p. The 'p0i' value is -1/p mod 2^31.
   * It is required that p is an odd integer.
   */
  static inline uint32_t
  modp_montymul(uint32_t a, uint32_t b, uint32_t p, uint32_t p0i) {
      uint64_t z, w;
      uint32_t d;
  
      z = (uint64_t)a * (uint64_t)b;
      w = ((z * p0i) & (uint64_t)0x7FFFFFFF) * p;
      d = (uint32_t)((z + w) >> 31) - p;
      d += p & -(d >> 31);
      return d;
  }
  
  /*
   * Compute R2 = 2^62 mod p.
   */
  static uint32_t
  modp_R2(uint32_t p, uint32_t p0i) {
      uint32_t z;
  
      /*
       * Compute z = 2^31 mod p (this is the value 1 in Montgomery
       * representation), then double it with an addition.
       */
      z = modp_R(p);
      z = modp_add(z, z, p);
  
      /*
       * Square it five times to obtain 2^32 in Montgomery representation
       * (i.e. 2^63 mod p).
       */
      z = modp_montymul(z, z, p, p0i);
      z = modp_montymul(z, z, p, p0i);
      z = modp_montymul(z, z, p, p0i);
      z = modp_montymul(z, z, p, p0i);
      z = modp_montymul(z, z, p, p0i);
  
      /*
       * Halve the value mod p to get 2^62.
       */
      z = (z + (p & -(z & 1))) >> 1;
      return z;
  }
  
  /*
   * Compute 2^(31*x) modulo p. This works for integers x up to 2^11.
   * p must be prime such that 2^30 < p < 2^31; p0i must be equal to
   * -1/p mod 2^31; R2 must be equal to 2^62 mod p.
   */
  static inline uint32_t
  modp_Rx(unsigned x, uint32_t p, uint32_t p0i, uint32_t R2) {
      int i;
      uint32_t r, z;
  
      /*
       * 2^(31*x) = (2^31)*(2^(31*(x-1))); i.e. we want the Montgomery
       * representation of (2^31)^e mod p, where e = x-1.
       * R2 is 2^31 in Montgomery representation.
       */
      x --;
      r = R2;
      z = modp_R(p);
      for (i = 0; (1U << i) <= x; i ++) {
          if ((x & (1U << i)) != 0) {
              z = modp_montymul(z, r, p, p0i);
          }
          r = modp_montymul(r, r, p, p0i);
      }
      return z;
  }
  
  /*
   * Division modulo p. If the divisor (b) is 0, then 0 is returned.
   * This function computes proper results only when p is prime.
   * Parameters:
   *   a     dividend
   *   b     divisor
   *   p     odd prime modulus
   *   p0i   -1/p mod 2^31
   *   R     2^31 mod R
   */
  static uint32_t
  modp_div(uint32_t a, uint32_t b, uint32_t p, uint32_t p0i, uint32_t R) {
      uint32_t z, e;
      int i;
  
      e = p - 2;
      z = R;
      for (i = 30; i >= 0; i --) {
          uint32_t z2;
  
          z = modp_montymul(z, z, p, p0i);
          z2 = modp_montymul(z, b, p, p0i);
          z ^= (z ^ z2) & -(uint32_t)((e >> i) & 1);
      }
  
      /*
       * The loop above just assumed that b was in Montgomery
       * representation, i.e. really contained b*R; under that
       * assumption, it returns 1/b in Montgomery representation,
       * which is R/b. But we gave it b in normal representation,
       * so the loop really returned R/(b/R) = R^2/b.
       *
       * We want a/b, so we need one Montgomery multiplication with a,
       * which also remove one of the R factors, and another such
       * multiplication to remove the second R factor.
       */
      z = modp_montymul(z, 1, p, p0i);
      return modp_montymul(a, z, p, p0i);
  }
  
  /*
   * Bit-reversal index table.
   */
  static const uint16_t REV10[] = {
      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
  };
  
  /*
   * Compute the roots for NTT and inverse NTT (binary case). Input
   * parameter g is a primitive 2048-th root of 1 modulo p (i.e. g^1024 =
   * -1 mod p). This fills gm[] and igm[] with powers of g and 1/g:
   *   gm[rev(i)] = g^i mod p
   *   igm[rev(i)] = (1/g)^i mod p
   * where rev() is the "bit reversal" function over 10 bits. It fills
   * the arrays only up to N = 2^logn values.
   *
   * The values stored in gm[] and igm[] are in Montgomery representation.
   *
   * p must be a prime such that p = 1 mod 2048.
   */
  static void
  modp_mkgm2(uint32_t *gm, uint32_t *igm, unsigned logn,
             uint32_t g, uint32_t p, uint32_t p0i) {
      size_t u, n;
      unsigned k;
      uint32_t ig, x1, x2, R2;
  
      n = (size_t)1 << logn;
  
      /*
       * We want g such that g^(2N) = 1 mod p, but the provided
       * generator has order 2048. We must square it a few times.
       */
      R2 = modp_R2(p, p0i);
      g = modp_montymul(g, R2, p, p0i);
      for (k = logn; k < 10; k ++) {
          g = modp_montymul(g, g, p, p0i);
      }
  
      ig = modp_div(R2, g, p, p0i, modp_R(p));
      k = 10 - logn;
      x1 = x2 = modp_R(p);
      for (u = 0; u < n; u ++) {
          size_t v;
  
          v = REV10[u << k];
          gm[v] = x1;
          igm[v] = x2;
          x1 = modp_montymul(x1, g, p, p0i);
          x2 = modp_montymul(x2, ig, p, p0i);
      }
  }
  
  /*
   * Compute the NTT over a polynomial (binary case). Polynomial elements
   * are a[0], a[stride], a[2 * stride]...
   */
  static void
  modp_NTT2_ext(uint32_t *a, size_t stride, const uint32_t *gm, unsigned logn,
                uint32_t p, uint32_t p0i) {
      size_t t, m, n;
  
      if (logn == 0) {
          return;
      }
      n = (size_t)1 << logn;
      t = n;
      for (m = 1; m < n; m <<= 1) {
          size_t ht, u, v1;
  
          ht = t >> 1;
          for (u = 0, v1 = 0; u < m; u ++, v1 += t) {
              uint32_t s;
              size_t v;
              uint32_t *r1, *r2;
  
              s = gm[m + u];
              r1 = a + v1 * stride;
              r2 = r1 + ht * stride;
              for (v = 0; v < ht; v ++, r1 += stride, r2 += stride) {
                  uint32_t x, y;
  
                  x = *r1;
                  y = modp_montymul(*r2, s, p, p0i);
                  *r1 = modp_add(x, y, p);
                  *r2 = modp_sub(x, y, p);
              }
          }
          t = ht;
      }
  }
  
  /*
   * Compute the inverse NTT over a polynomial (binary case).
   */
  static void
  modp_iNTT2_ext(uint32_t *a, size_t stride, const uint32_t *igm, unsigned logn,
                 uint32_t p, uint32_t p0i) {
      size_t t, m, n, k;
      uint32_t ni;
      uint32_t *r;
  
      if (logn == 0) {
          return;
      }
      n = (size_t)1 << logn;
      t = 1;
      for (m = n; m > 1; m >>= 1) {
          size_t hm, dt, u, v1;
  
          hm = m >> 1;
          dt = t << 1;
          for (u = 0, v1 = 0; u < hm; u ++, v1 += dt) {
              uint32_t s;
              size_t v;
              uint32_t *r1, *r2;
  
              s = igm[hm + u];
              r1 = a + v1 * stride;
              r2 = r1 + t * stride;
              for (v = 0; v < t; v ++, r1 += stride, r2 += stride) {
                  uint32_t x, y;
  
                  x = *r1;
                  y = *r2;
                  *r1 = modp_add(x, y, p);
                  *r2 = modp_montymul(
                            modp_sub(x, y, p), s, p, p0i);;
              }
          }
          t = dt;
      }
  
      /*
       * We need 1/n in Montgomery representation, i.e. R/n. Since
       * 1 <= logn <= 10, R/n is an integer; morever, R/n <= 2^30 < p,
       * thus a simple shift will do.
       */
      ni = (uint32_t)1 << (31 - logn);
      for (k = 0, r = a; k < n; k ++, r += stride) {
          *r = modp_montymul(*r, ni, p, p0i);
      }
  }
  
  /*
   * Simplified macros for NTT and iNTT (binary case) when the elements
   * are consecutive in RAM.
   */
  #define modp_NTT2(a, gm, logn, p, p0i)   modp_NTT2_ext(a, 1, gm, logn, p, p0i)
  #define modp_iNTT2(a, igm, logn, p, p0i) modp_iNTT2_ext(a, 1, igm, logn, p, p0i)
  
  /*
   * Given polynomial f in NTT representation modulo p, compute f' of degree
   * less than N/2 such that f' = f0^2 - X*f1^2, where f0 and f1 are
   * polynomials of degree less than N/2 such that f = f0(X^2) + X*f1(X^2).
   *
   * The new polynomial is written "in place" over the first N/2 elements
   * of f.
   *
   * If applied logn times successively on a given polynomial, the resulting
   * degree-0 polynomial is the resultant of f and X^N+1 modulo p.
   *
   * This function applies only to the binary case; it is invoked from
   * solve_NTRU_binary_depth1().
   */
  static void
  modp_poly_rec_res(uint32_t *f, unsigned logn,
                    uint32_t p, uint32_t p0i, uint32_t R2) {
      size_t hn, u;
  
      hn = (size_t)1 << (logn - 1);
      for (u = 0; u < hn; u ++) {
          uint32_t w0, w1;
  
          w0 = f[(u << 1) + 0];
          w1 = f[(u << 1) + 1];
          f[u] = modp_montymul(modp_montymul(w0, w1, p, p0i), R2, p, p0i);
      }
  }
  
  /* ==================================================================== */
  /*
   * Custom bignum implementation.
   *
   * This is a very reduced set of functionalities. We need to do the
   * following operations:
   *
   *  - Rebuild the resultant and the polynomial coefficients from their
   *    values modulo small primes (of length 31 bits each).
   *
   *  - Compute an extended GCD between the two computed resultants.
   *
   *  - Extract top bits and add scaled values during the successive steps
   *    of Babai rounding.
   *
   * When rebuilding values using CRT, we must also recompute the product
   * of the small prime factors. We always do it one small factor at a
   * time, so the "complicated" operations can be done modulo the small
   * prime with the modp_* functions. CRT coefficients (inverses) are
   * precomputed.
   *
   * All values are positive until the last step: when the polynomial
   * coefficients have been rebuilt, we normalize them around 0. But then,
   * only additions and subtractions on the upper few bits are needed
   * afterwards.
   *
   * We keep big integers as arrays of 31-bit words (in uint32_t values);
   * the top bit of each uint32_t is kept equal to 0. Using 31-bit words
   * makes it easier to keep track of carries. When negative values are
   * used, two's complement is used.
   */
  
  /*
   * Subtract integer b from integer a. Both integers are supposed to have
   * the same size. The carry (0 or 1) is returned. Source arrays a and b
   * MUST be distinct.
   *
   * The operation is performed as described above if ctr = 1. If
   * ctl = 0, the value a[] is unmodified, but all memory accesses are
   * still performed, and the carry is computed and returned.
   */
  static uint32_t
  zint_sub(uint32_t *a, const uint32_t *b, size_t len,
           uint32_t ctl) {
      size_t u;
      uint32_t cc, m;
  
      cc = 0;
      m = -ctl;
      for (u = 0; u < len; u ++) {
          uint32_t aw, w;
  
          aw = a[u];
          w = aw - b[u] - cc;
          cc = w >> 31;
          aw ^= ((w & 0x7FFFFFFF) ^ aw) & m;
          a[u] = aw;
      }
      return cc;
  }
  
  /*
   * Mutiply the provided big integer m with a small value x.
   * This function assumes that x < 2^31. The carry word is returned.
   */
  static uint32_t
  zint_mul_small(uint32_t *m, size_t mlen, uint32_t x) {
      size_t u;
      uint32_t cc;
  
      cc = 0;
      for (u = 0; u < mlen; u ++) {
          uint64_t z;
  
          z = (uint64_t)m[u] * (uint64_t)x + cc;
          m[u] = (uint32_t)z & 0x7FFFFFFF;
          cc = (uint32_t)(z >> 31);
      }
      return cc;
  }
  
  /*
   * Reduce a big integer d modulo a small integer p.
   * Rules:
   *  d is unsigned
   *  p is prime
   *  2^30 < p < 2^31
   *  p0i = -(1/p) mod 2^31
   *  R2 = 2^62 mod p
   */
  static uint32_t
  zint_mod_small_unsigned(const uint32_t *d, size_t dlen,
                          uint32_t p, uint32_t p0i, uint32_t R2) {
      uint32_t x;
      size_t u;
  
      /*
       * Algorithm: we inject words one by one, starting with the high
       * word. Each step is:
       *  - multiply x by 2^31
       *  - add new word
       */
      x = 0;
      u = dlen;
      while (u -- > 0) {
          uint32_t w;
  
          x = modp_montymul(x, R2, p, p0i);
          w = d[u] - p;
          w += p & -(w >> 31);
          x = modp_add(x, w, p);
      }
      return x;
  }
  
  /*
   * Similar to zint_mod_small_unsigned(), except that d may be signed.
   * Extra parameter is Rx = 2^(31*dlen) mod p.
   */
  static uint32_t
  zint_mod_small_signed(const uint32_t *d, size_t dlen,
                        uint32_t p, uint32_t p0i, uint32_t R2, uint32_t Rx) {
      uint32_t z;
  
      if (dlen == 0) {
          return 0;
      }
      z = zint_mod_small_unsigned(d, dlen, p, p0i, R2);
      z = modp_sub(z, Rx & -(d[dlen - 1] >> 30), p);
      return z;
  }
  
  /*
   * Add y*s to x. x and y initially have length 'len' words; the new x
   * has length 'len+1' words. 's' must fit on 31 bits. x[] and y[] must
   * not overlap.
   */
  static void
  zint_add_mul_small(uint32_t *x,
                     const uint32_t *y, size_t len, uint32_t s) {
      size_t u;
      uint32_t cc;
  
      cc = 0;
      for (u = 0; u < len; u ++) {
          uint32_t xw, yw;
          uint64_t z;
  
          xw = x[u];
          yw = y[u];
          z = (uint64_t)yw * (uint64_t)s + (uint64_t)xw + (uint64_t)cc;
          x[u] = (uint32_t)z & 0x7FFFFFFF;
          cc = (uint32_t)(z >> 31);
      }
      x[len] = cc;
  }
  
  /*
   * Normalize a modular integer around 0: if x > p/2, then x is replaced
   * with x - p (signed encoding with two's complement); otherwise, x is
   * untouched. The two integers x and p are encoded over the same length.
   */
  static void
  zint_norm_zero(uint32_t *x, const uint32_t *p, size_t len) {
      size_t u;
      uint32_t r, bb;
  
      /*
       * 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.
       */
      r = 0;
      bb = 0;
      u = len;
      while (u -- > 0) {
          uint32_t wx, wp, cc;
  
          /*
           * Get the two words to compare in wx and wp (both over
           * 31 bits exactly).
           */
          wx = x[u];
          wp = (p[u] >> 1) | (bb << 30);
          bb = p[u] & 1;
  
          /*
           * We set cc to -1, 0 or 1, depending on whether wp is
           * lower than, equal to, or greater than wx.
           */
          cc = wp - wx;
          cc = ((-cc) >> 31) | -(cc >> 31);
  
          /*
           * 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);
      }
  
      /*
       * 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.
       */
      zint_sub(x, p, len, r >> 31);
  }
  
  /*
   * Rebuild integers from their RNS representation. There are 'num'
   * integers, and each consists in 'xlen' words. 'xx' points at that
   * first word of the first integer; subsequent integers are accessed
   * by adding 'xstride' repeatedly.
   *
   * The words of an integer are the RNS representation of that integer,
   * using the provided 'primes' are moduli. This function replaces
   * each integer with its multi-word value (little-endian order).
   *
   * If "normalize_signed" is non-zero, then the returned value is
   * normalized to the -m/2..m/2 interval (where m is the product of all
   * small prime moduli); two's complement is used for negative values.
   */
  static void
  zint_rebuild_CRT(uint32_t *xx, size_t xlen, size_t xstride,
                   size_t num, const small_prime *primes, int normalize_signed,
                   uint32_t *tmp) {
      size_t u;
      uint32_t *x;
  
      tmp[0] = primes[0].p;
      for (u = 1; u < xlen; u ++) {
          /*
           * At the entry of each loop iteration:
           *  - the first u words of each array have been
           *    reassembled;
           *  - the first u words of tmp[] contains the
           * product of the prime moduli processed so far.
           *
           * We call 'q' the product of all previous primes.
           */
          uint32_t p, p0i, s, R2;
          size_t v;
  
          p = primes[u].p;
          s = primes[u].s;
          p0i = modp_ninv31(p);
          R2 = modp_R2(p, p0i);
  
          for (v = 0, x = xx; v < num; v ++, x += xstride) {
              uint32_t xp, xq, xr;
              /*
               * xp = the integer x modulo the prime p for this
               *      iteration
               * xq = (x mod q) mod p
               */
              xp = x[u];
              xq = zint_mod_small_unsigned(x, u, p, p0i, R2);
  
              /*
               * New value is (x mod q) + q * (s * (xp - xq) mod p)
               */
              xr = modp_montymul(s, modp_sub(xp, xq, p), p, p0i);
              zint_add_mul_small(x, tmp, u, xr);
          }
  
          /*
           * Update product of primes in tmp[].
           */
          tmp[u] = zint_mul_small(tmp, u, p);
      }
  
      /*
       * Normalize the reconstructed values around 0.
       */
      if (normalize_signed) {
          for (u = 0, x = xx; u < num; u ++, x += xstride) {
              zint_norm_zero(x, tmp, xlen);
          }
      }
  }
  
  /*
   * Negate a big integer conditionally: value a is replaced with -a if
   * and only if ctl = 1. Control value ctl must be 0 or 1.
   */
  static void
  zint_negate(uint32_t *a, size_t len, uint32_t ctl) {
      size_t u;
      uint32_t cc, m;
  
      /*
       * If ctl = 1 then we flip the bits of a by XORing with
       * 0x7FFFFFFF, and we add 1 to the value. If ctl = 0 then we XOR
       * with 0 and add 0, which leaves the value unchanged.
       */
      cc = ctl;
      m = -ctl >> 1;
      for (u = 0; u < len; u ++) {
          uint32_t aw;
  
          aw = a[u];
          aw = (aw ^ m) + cc;
          a[u] = aw & 0x7FFFFFFF;
          cc = aw >> 31;
      }
  }
  
  /*
   * Replace a with (a*xa+b*xb)/(2^31) and b with (a*ya+b*yb)/(2^31).
   * The low bits are dropped (the caller should compute the coefficients
   * such that these dropped bits are all zeros). If either or both
   * yields a negative value, then the value is negated.
   *
   * Returned value is:
   *  0  both values were positive
   *  1  new a had to be negated
   *  2  new b had to be negated
   *  3  both new a and new b had to be negated
   *
   * Coefficients xa, xb, ya and yb may use the full signed 32-bit range.
   */
  static uint32_t
  zint_co_reduce(uint32_t *a, uint32_t *b, size_t len,
                 int64_t xa, int64_t xb, int64_t ya, int64_t yb) {
      size_t u;
      int64_t cca, ccb;
      uint32_t nega, negb;
  
      cca = 0;
      ccb = 0;
      for (u = 0; u < len; u ++) {
          uint32_t wa, wb;
          uint64_t za, zb;
  
          wa = a[u];
          wb = b[u];
          za = wa * (uint64_t)xa + wb * (uint64_t)xb + (uint64_t)cca;
          zb = wa * (uint64_t)ya + wb * (uint64_t)yb + (uint64_t)ccb;
          if (u > 0) {
              a[u - 1] = (uint32_t)za & 0x7FFFFFFF;
              b[u - 1] = (uint32_t)zb & 0x7FFFFFFF;
          }
          cca = *(int64_t *)&za >> 31;
          ccb = *(int64_t *)&zb >> 31;
      }
      a[len - 1] = (uint32_t)cca;
      b[len - 1] = (uint32_t)ccb;
  
      nega = (uint32_t)((uint64_t)cca >> 63);
      negb = (uint32_t)((uint64_t)ccb >> 63);
      zint_negate(a, len, nega);
      zint_negate(b, len, negb);
      return nega | (negb << 1);
  }
  
  /*
   * Finish modular reduction. Rules on input parameters:
   *
   *   if neg = 1, then -m <= a < 0
   *   if neg = 0, then 0 <= a < 2*m
   *
   * If neg = 0, then the top word of a[] is allowed to use 32 bits.
   *
   * Modulus m must be odd.
   */
  static void
  zint_finish_mod(uint32_t *a, size_t len, const uint32_t *m, uint32_t neg) {
      size_t u;
      uint32_t cc, xm, ym;
  
      /*
       * First pass: compare a (assumed nonnegative) with m. Note that
       * if the top word uses 32 bits, subtracting m must yield a
       * value less than 2^31 since a < 2*m.
       */
      cc = 0;
      for (u = 0; u < len; u ++) {
          cc = (a[u] - m[u] - cc) >> 31;
      }
  
      /*
       * If neg = 1 then we must add m (regardless of cc)
       * If neg = 0 and cc = 0 then we must subtract m
       * If neg = 0 and cc = 1 then we must do nothing
       *
       * In the loop below, we conditionally subtract either m or -m
       * from a. Word xm is a word of m (if neg = 0) or -m (if neg = 1);
       * but if neg = 0 and cc = 1, then ym = 0 and it forces mw to 0.
       */
      xm = -neg >> 1;
      ym = -(neg | (1 - cc));
      cc = neg;
      for (u = 0; u < len; u ++) {
          uint32_t aw, mw;
  
          aw = a[u];
          mw = (m[u] ^ xm) & ym;
          aw = aw - mw - cc;
          a[u] = aw & 0x7FFFFFFF;
          cc = aw >> 31;
      }
  }
  
  /*
   * Replace a with (a*xa+b*xb)/(2^31) mod m, and b with
   * (a*ya+b*yb)/(2^31) mod m. Modulus m must be odd; m0i = -1/m[0] mod 2^31.
   */
  static void
  zint_co_reduce_mod(uint32_t *a, uint32_t *b, const uint32_t *m, size_t len,
                     uint32_t m0i, int64_t xa, int64_t xb, int64_t ya, int64_t yb) {
      size_t u;
      int64_t cca, ccb;
      uint32_t fa, fb;
  
      /*
       * These are actually four combined Montgomery multiplications.
       */
      cca = 0;
      ccb = 0;
      fa = ((a[0] * (uint32_t)xa + b[0] * (uint32_t)xb) * m0i) & 0x7FFFFFFF;
      fb = ((a[0] * (uint32_t)ya + b[0] * (uint32_t)yb) * m0i) & 0x7FFFFFFF;
      for (u = 0; u < len; u ++) {
          uint32_t wa, wb;
          uint64_t za, zb;
  
          wa = a[u];
          wb = b[u];
          za = wa * (uint64_t)xa + wb * (uint64_t)xb
               + m[u] * (uint64_t)fa + (uint64_t)cca;
          zb = wa * (uint64_t)ya + wb * (uint64_t)yb
               + m[u] * (uint64_t)fb + (uint64_t)ccb;
          if (u > 0) {
              a[u - 1] = (uint32_t)za & 0x7FFFFFFF;
              b[u - 1] = (uint32_t)zb & 0x7FFFFFFF;
          }
          cca = *(int64_t *)&za >> 31;
          ccb = *(int64_t *)&zb >> 31;
      }
      a[len - 1] = (uint32_t)cca;
      b[len - 1] = (uint32_t)ccb;
  
      /*
       * At this point:
       *   -m <= a < 2*m
       *   -m <= b < 2*m
       * (this is a case of Montgomery reduction)
       * The top words of 'a' and 'b' may have a 32-th bit set.
       * We want to add or subtract the modulus, as required.
       */
      zint_finish_mod(a, len, m, (uint32_t)((uint64_t)cca >> 63));
      zint_finish_mod(b, len, m, (uint32_t)((uint64_t)ccb >> 63));
  }
  
  /*
   * Compute a GCD between two positive big integers x and y. The two
   * integers must be odd. Returned value is 1 if the GCD is 1, 0
   * otherwise. When 1 is returned, arrays u and v are filled with values
   * such that:
   *   0 <= u <= y
   *   0 <= v <= x
   *   x*u - y*v = 1
   * x[] and y[] are unmodified. Both input values must have the same
   * encoded length. Temporary array must be large enough to accommodate 4
   * extra values of that length. Arrays u, v and tmp may not overlap with
   * each other, or with either x or y.
   */
  static int
  zint_bezout(uint32_t *u, uint32_t *v,
              const uint32_t *x, const uint32_t *y,
              size_t len, uint32_t *tmp) {
      /*
       * Algorithm is an extended binary GCD. We maintain 6 values
       * a, b, u0, u1, v0 and v1 with 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
       *
       * Initial values are:
       *
       *  a = x   u0 = 1   v0 = 0
       *  b = y   u1 = y   v1 = x-1
       *
       * Each iteration reduces either a or b, and maintains the
       * invariants. Algorithm stops when a = b, at which point their
       * common value is GCD(a,b) and (u0,v0) (or (u1,v1)) contains
       * the values (u,v) we want to return.
       *
       * The formal definition of the algorithm is a sequence of steps:
       *
       *  - If a is even, then:
       *        a <- a/2
       *        u0 <- u0/2 mod y
       *        v0 <- v0/2 mod x
       *
       *  - Otherwise, if b is even, then:
       *        b <- b/2
       *        u1 <- u1/2 mod y
       *        v1 <- v1/2 mod x
       *
       *  - Otherwise, if a > b, then:
       *        a <- (a-b)/2
       *        u0 <- (u0-u1)/2 mod y
       *        v0 <- (v0-v1)/2 mod x
       *
       *  - Otherwise:
       *        b <- (b-a)/2
       *        u1 <- (u1-u0)/2 mod y
       *        v1 <- (v1-v0)/2 mod y
       *
       * We can show that the operations above preserve the invariants:
       *
       *  - If a is even, then u0 and v0 are either both even or both
       *    odd (since a = x*u0 - y*v0, and x and y are both odd).
       *    If u0 and v0 are both even, then (u0,v0) <- (u0/2,v0/2).
       *    Otherwise, (u0,v0) <- ((u0+y)/2,(v0+x)/2). Either way,
       *    the a = x*u0 - y*v0 invariant is preserved.
       *
       *  - The same holds for the case where b is even.
       *
       *  - If a and b are odd, and a > b, then:
       *
       *      a-b = x*(u0-u1) - y*(v0-v1)
       *
       *    In that situation, if u0 < u1, then x*(u0-u1) < 0, but
       *    a-b > 0; therefore, it must be that v0 < v1, and the
       *    first part of the update is: (u0,v0) <- (u0-u1+y,v0-v1+x),
       *    which preserves the invariants. Otherwise, if u0 > u1,
       *    then u0-u1 >= 1, thus x*(u0-u1) >= x. But a <= x and
       *    b >= 0, hence a-b <= x. It follows that, in that case,
       *    v0-v1 >= 0. The first part of the update is then:
       *    (u0,v0) <- (u0-u1,v0-v1), which again preserves the
       *    invariants.
       *
       *    Either way, once the subtraction is done, the new value of
       *    a, which is the difference of two odd values, is even,
       *    and the remaining of this step is a subcase of the
       *    first algorithm case (i.e. when a is even).
       *
       *  - If a and b are odd, and b > a, then the a similar
       *    argument holds.
       *
       * The values a and b start at x and y, respectively. Since x
       * and y are odd, their GCD is odd, and it is easily seen that
       * all steps conserve the GCD (GCD(a-b,b) = GCD(a, b);
       * GCD(a/2,b) = GCD(a,b) if GCD(a,b) is odd). Moreover, either a
       * or b is reduced by at least one bit at each iteration, so
       * the algorithm necessarily converges on the case a = b, at
       * which point the common value is the GCD.
       *
       * In the algorithm expressed above, when a = b, the fourth case
       * applies, and sets b = 0. Since a contains the GCD of x and y,
       * which are both odd, a must be odd, and subsequent iterations
       * (if any) will simply divide b by 2 repeatedly, which has no
       * consequence. Thus, the algorithm can run for more iterations
       * than necessary; the final GCD will be in a, and the (u,v)
       * coefficients will be (u0,v0).
       *
       *
       * The presentation above is bit-by-bit. It can be sped up by
       * noticing that all decisions are taken based on the low bits
       * and high bits of a and b. We can extract the two top words
       * and low word of each of a and b, and compute reduction
       * parameters pa, pb, qa and qb such that the new values for
       * a and b are:
       *    a' = (a*pa + b*pb) / (2^31)
       *    b' = (a*qa + b*qb) / (2^31)
       * the two divisions being exact. The coefficients are obtained
       * just from the extracted words, and may be slightly off, requiring
       * an optional correction: if a' < 0, then we replace pa with -pa
       * and pb with -pb. Each such step will reduce the total length
       * (sum of lengths of a and b) by at least 30 bits at each
       * iteration.
       */
      uint32_t *u0, *u1, *v0, *v1, *a, *b;
      uint32_t x0i, y0i;
      uint32_t num, rc;
      size_t j;
  
      if (len == 0) {
          return 0;
      }
  
      /*
       * u0 and v0 are the u and v result buffers; the four other
       * values (u1, v1, a and b) are taken from tmp[].
       */
      u0 = u;
      v0 = v;
      u1 = tmp;
      v1 = u1 + len;
      a = v1 + len;
      b = a + len;
  
      /*
       * We'll need the Montgomery reduction coefficients.
       */
      x0i = modp_ninv31(x[0]);
      y0i = modp_ninv31(y[0]);
  
      /*
       * Initialize a, b, u0, u1, v0 and v1.
       *  a = x   u0 = 1   v0 = 0
       *  b = y   u1 = y   v1 = x-1
       * Note that x is odd, so computing x-1 is easy.
       */
      memcpy(a, x, len * sizeof * x);
      memcpy(b, y, len * sizeof * y);
      u0[0] = 1;
      memset(u0 + 1, 0, (len - 1) * sizeof * u0);
      memset(v0, 0, len * sizeof * v0);
      memcpy(u1, y, len * sizeof * u1);
      memcpy(v1, x, len * sizeof * v1);
      v1[0] --;
  
      /*
       * Each input operand may be as large as 31*len bits, and we
       * reduce the total length by at least 30 bits at each iteration.
       */
      for (num = 62 * (uint32_t)len + 30; num >= 30; num -= 30) {
          uint32_t c0, c1;
          uint32_t a0, a1, b0, b1;
          uint64_t a_hi, b_hi;
          uint32_t a_lo, b_lo;
          int64_t pa, pb, qa, qb;
          int i;
          uint32_t r;
  
          /*
           * Extract the top words of a and b. If j is the highest
           * index >= 1 such that a[j] != 0 or b[j] != 0, then we
           * want (a[j] << 31) + a[j-1] and (b[j] << 31) + b[j-1].
           * If a and b are down to one word each, then we use
           * a[0] and b[0].
           */
          c0 = (uint32_t) -1;
          c1 = (uint32_t) -1;
          a0 = 0;
          a1 = 0;
          b0 = 0;
          b1 = 0;
          j = len;
          while (j -- > 0) {
              uint32_t aw, bw;
  
              aw = a[j];
              bw = b[j];
              a0 ^= (a0 ^ aw) & c0;
              a1 ^= (a1 ^ aw) & c1;
              b0 ^= (b0 ^ bw) & c0;
              b1 ^= (b1 ^ bw) & c1;
              c1 = c0;
              c0 &= (((aw | bw) + 0x7FFFFFFF) >> 31) - (uint32_t)1;
          }
  
          /*
           * If c1 = 0, then we grabbed two words for a and b.
           * If c1 != 0 but c0 = 0, then we grabbed one word. It
           * is not possible that c1 != 0 and c0 != 0, because that
           * would mean that both integers are zero.
           */
          a1 |= a0 & c1;
          a0 &= ~c1;
          b1 |= b0 & c1;
          b0 &= ~c1;
          a_hi = ((uint64_t)a0 << 31) + a1;
          b_hi = ((uint64_t)b0 << 31) + b1;
          a_lo = a[0];
          b_lo = b[0];
  
          /*
           * Compute reduction factors:
           *
           *   a' = a*pa + b*pb
           *   b' = a*qa + b*qb
           *
           * such that a' and b' are both multiple of 2^31, but are
           * only marginally larger than a and b.
           */
          pa = 1;
          pb = 0;
          qa = 0;
          qb = 1;
          for (i = 0; i < 31; i ++) {
              /*
               * At each iteration:
               *
               *   a <- (a-b)/2 if: a is odd, b is odd, a_hi > b_hi
               *   b <- (b-a)/2 if: a is odd, b is odd, a_hi <= b_hi
               *   a <- a/2 if: a is even
               *   b <- b/2 if: a is odd, b is even
               *
               * We multiply a_lo and b_lo by 2 at each
               * iteration, thus a division by 2 really is a
               * non-multiplication by 2.
               */
              uint32_t rt, oa, ob, cAB, cBA, cA;
              uint64_t rz;
  
              /*
               * rt = 1 if a_hi > b_hi, 0 otherwise.
               */
              rz = b_hi - a_hi;
              rt = (uint32_t)((rz ^ ((a_hi ^ b_hi)
                                     & (a_hi ^ rz))) >> 63);
  
              /*
               * cAB = 1 if b must be subtracted from a
               * cBA = 1 if a must be subtracted from b
               * cA = 1 if a must be divided by 2
               *
               * Rules:
               *
               *   cAB and cBA cannot both be 1.
               *   If a is not divided by 2, b is.
               */
              oa = (a_lo >> i) & 1;
              ob = (b_lo >> i) & 1;
              cAB = oa & ob & rt;
              cBA = oa & ob & ~rt;
              cA = cAB | (oa ^ 1);
  
              /*
               * Conditional subtractions.
               */
              a_lo -= b_lo & -cAB;
              a_hi -= b_hi & -(uint64_t)cAB;
              pa -= qa & -(int64_t)cAB;
              pb -= qb & -(int64_t)cAB;
              b_lo -= a_lo & -cBA;
              b_hi -= a_hi & -(uint64_t)cBA;
              qa -= pa & -(int64_t)cBA;
              qb -= pb & -(int64_t)cBA;
  
              /*
               * Shifting.
               */
              a_lo += a_lo & (cA - 1);
              pa += pa & ((int64_t)cA - 1);
              pb += pb & ((int64_t)cA - 1);
              a_hi ^= (a_hi ^ (a_hi >> 1)) & -(uint64_t)cA;
              b_lo += b_lo & -cA;
              qa += qa & -(int64_t)cA;
              qb += qb & -(int64_t)cA;
              b_hi ^= (b_hi ^ (b_hi >> 1)) & ((uint64_t)cA - 1);
          }
  
          /*
           * Apply the computed parameters to our values. We
           * may have to correct pa and pb depending on the
           * returned value of zint_co_reduce() (when a and/or b
           * had to be negated).
           */
          r = zint_co_reduce(a, b, len, pa, pb, qa, qb);
          pa -= (pa + pa) & -(int64_t)(r & 1);
          pb -= (pb + pb) & -(int64_t)(r & 1);
          qa -= (qa + qa) & -(int64_t)(r >> 1);
          qb -= (qb + qb) & -(int64_t)(r >> 1);
          zint_co_reduce_mod(u0, u1, y, len, y0i, pa, pb, qa, qb);
          zint_co_reduce_mod(v0, v1, x, len, x0i, pa, pb, qa, qb);
      }
  
      /*
       * At that point, array a[] should contain the GCD, and the
       * results (u,v) should already be set. We check that the GCD
       * is indeed 1. We also check that the two operands x and y
       * are odd.
       */
      rc = a[0] ^ 1;
      for (j = 1; j < len; j ++) {
          rc |= a[j];
      }
      return (int)((1 - ((rc | -rc) >> 31)) & x[0] & y[0]);
  }
  
  /*
   * Add k*y*2^sc to x. The result is assumed to fit in the array of
   * size xlen (truncation is applied if necessary).
   * Scale factor 'sc' is provided as sch and scl, such that:
   *   sch = sc / 31
   *   scl = sc % 31
   * xlen MUST NOT be lower than ylen.
   *
   * x[] and y[] are both signed integers, using two's complement for
   * negative values.
   */
  static void
  zint_add_scaled_mul_small(uint32_t *x, size_t xlen,
                            const uint32_t *y, size_t ylen, int32_t k,
                            uint32_t sch, uint32_t scl) {
      size_t u;
      uint32_t ysign, tw;
      int32_t cc;
  
      if (ylen == 0) {
          return;
      }
  
      ysign = -(y[ylen - 1] >> 30) >> 1;
      tw = 0;
      cc = 0;
      for (u = sch; u < xlen; u ++) {
          size_t v;
          uint32_t wy, wys, ccu;
          uint64_t z;
  
          /*
           * Get the next word of y (scaled).
           */
          v = u - sch;
          if (v < ylen) {
              wy = y[v];
          } else {
              wy = ysign;
          }
          wys = ((wy << scl) & 0x7FFFFFFF) | tw;
          tw = wy >> (31 - scl);
  
          /*
           * The expression below does not overflow.
           */
          z = (uint64_t)((int64_t)wys * (int64_t)k + (int64_t)x[u] + cc);
          x[u] = (uint32_t)z & 0x7FFFFFFF;
  
          /*
           * Right-shifting the signed value z would yield
           * implementation-defined results (arithmetic shift is
           * not guaranteed). However, we can cast to unsigned,
           * and get the next carry as an unsigned word. We can
           * then convert it back to signed by using the guaranteed
           * fact that 'int32_t' uses two's complement with no
           * trap representation or padding bit, and with a layout
           * compatible with that of 'uint32_t'.
           */
          ccu = (uint32_t)(z >> 31);
          cc = *(int32_t *)&ccu;
      }
  }
  
  /*
   * Subtract y*2^sc from x. The result is assumed to fit in the array of
   * size xlen (truncation is applied if necessary).
   * Scale factor 'sc' is provided as sch and scl, such that:
   *   sch = sc / 31
   *   scl = sc % 31
   * xlen MUST NOT be lower than ylen.
   *
   * x[] and y[] are both signed integers, using two's complement for
   * negative values.
   */
  static void
  zint_sub_scaled(uint32_t *x, size_t xlen,
                  const uint32_t *y, size_t ylen, uint32_t sch, uint32_t scl) {
      size_t u;
      uint32_t ysign, tw;
      uint32_t cc;
  
      if (ylen == 0) {
          return;
      }
  
      ysign = -(y[ylen - 1] >> 30) >> 1;
      tw = 0;
      cc = 0;
      for (u = sch; u < xlen; u ++) {
          size_t v;
          uint32_t w, wy, wys;
  
          /*
           * Get the next word of y (scaled).
           */
          v = u - sch;
          if (v < ylen) {
              wy = y[v];
          } else {
              wy = ysign;
          }
          wys = ((wy << scl) & 0x7FFFFFFF) | tw;
          tw = wy >> (31 - scl);
  
          w = x[u] - wys - cc;
          x[u] = w & 0x7FFFFFFF;
          cc = w >> 31;
      }
  }
  
  /*
   * Convert a one-word signed big integer into a signed value.
   */
  static inline int32_t
  zint_one_to_plain(const uint32_t *x) {
      uint32_t w;
  
      w = x[0];
      w |= (w & 0x40000000) << 1;
      return *(int32_t *)&w;
  }
  
  /* ==================================================================== */
  
  /*
   * Convert a polynomial to floating-point values.
   *
   * Each coefficient has length flen words, and starts fstride words after
   * the previous.
   *
   * IEEE-754 binary64 values can represent values in a finite range,
   * roughly 2^(-1023) to 2^(+1023); thus, if coefficients are too large,
   * they should be "trimmed" by pointing not to the lowest word of each,
   * but upper.
   */
  static void
  poly_big_to_fp(fpr *d, const uint32_t *f, size_t flen, size_t fstride,
                 unsigned logn) {
      size_t n, u;
  
      n = MKN(logn);
      if (flen == 0) {
          for (u = 0; u < n; u ++) {
              d[u] = fpr_zero;
          }
          return;
      }
      for (u = 0; u < n; u ++, f += fstride) {
          size_t v;
          uint32_t neg, cc, xm;
          fpr x, fsc;
  
          /*
           * Get sign of the integer; if it is negative, then we
           * will load its absolute value instead, and negate the
           * result.
           */
          neg = -(f[flen - 1] >> 30);
          xm = neg >> 1;
          cc = neg & 1;
          x = fpr_zero;
          fsc = fpr_one;
          for (v = 0; v < flen; v ++, fsc = fpr_mul(fsc, fpr_ptwo31)) {
              uint32_t w;
  
              w = (f[v] ^ xm) + cc;
              cc = w >> 31;
              w &= 0x7FFFFFFF;
              w -= (w << 1) & neg;
              x = fpr_add(x, fpr_mul(fpr_of(*(int32_t *)&w), fsc));
          }
          d[u] = x;
      }
  }
  
  /*
   * Convert a polynomial to small integers. Source values are supposed
   * to be one-word integers, signed over 31 bits. Returned value is 0
   * if any of the coefficients exceeds the provided limit (in absolute
   * value), or 1 on success.
   *
   * This is not constant-time; this is not a problem here, because on
   * any failure, the NTRU-solving process will be deemed to have failed
   * and the (f,g) polynomials will be discarded.
   */
  static int
  poly_big_to_small(int8_t *d, const uint32_t *s, int lim, unsigned logn) {
      size_t n, u;
  
      n = MKN(logn);
      for (u = 0; u < n; u ++) {
          int32_t z;
  
          z = zint_one_to_plain(s + u);
          if (z < -lim || z > lim) {
              return 0;
          }
          d[u] = (int8_t)z;
      }
      return 1;
  }
  
  /*
   * Subtract k*f from F, where F, f and k are polynomials modulo X^N+1.
   * Coefficients of polynomial k are small integers (signed values in the
   * -2^31..2^31 range) scaled by 2^sc. Value sc is provided as sch = sc / 31
   * and scl = sc % 31.
   *
   * This function implements the basic quadratic multiplication algorithm,
   * which is efficient in space (no extra buffer needed) but slow at
   * high degree.
   */
  static void
  poly_sub_scaled(uint32_t *F, size_t Flen, size_t Fstride,
                  const uint32_t *f, size_t flen, size_t fstride,
                  const int32_t *k, uint32_t sch, uint32_t scl, unsigned logn) {
      size_t n, u;
  
      n = MKN(logn);
      for (u = 0; u < n; u ++) {
          int32_t kf;
          size_t v;
          uint32_t *x;
          const uint32_t *y;
  
          kf = -k[u];
          x = F + u * Fstride;
          y = f;
          for (v = 0; v < n; v ++) {
              zint_add_scaled_mul_small(
                  x, Flen, y, flen, kf, sch, scl);
              if (u + v == n - 1) {
                  x = F;
                  kf = -kf;
              } else {
                  x += Fstride;
              }
              y += fstride;
          }
      }
  }
  
  /*
   * Subtract k*f from F. Coefficients of polynomial k are small integers
   * (signed values in the -2^31..2^31 range) scaled by 2^sc. This function
   * assumes that the degree is large, and integers relatively small.
   * The value sc is provided as sch = sc / 31 and scl = sc % 31.
   */
  static void
  poly_sub_scaled_ntt(uint32_t *F, size_t Flen, size_t Fstride,
                      const uint32_t *f, size_t flen, size_t fstride,
                      const int32_t *k, uint32_t sch, uint32_t scl, unsigned logn,
                      uint32_t *tmp) {
      uint32_t *gm, *igm, *fk, *t1, *x;
      const uint32_t *y;
      size_t n, u, tlen;
      const small_prime *primes;
  
      n = MKN(logn);
      tlen = flen + 1;
      gm = tmp;
      igm = gm + MKN(logn);
      fk = igm + MKN(logn);
      t1 = fk + n * tlen;
  
      primes = PRIMES;
  
      /*
       * Compute k*f in fk[], in RNS notation.
       */
      for (u = 0; u < tlen; u ++) {
          uint32_t p, p0i, R2, Rx;
          size_t v;
  
          p = primes[u].p;
          p0i = modp_ninv31(p);
          R2 = modp_R2(p, p0i);
          Rx = modp_Rx((unsigned)flen, p, p0i, R2);
          modp_mkgm2(gm, igm, logn, primes[u].g, p, p0i);
  
          for (v = 0; v < n; v ++) {
              t1[v] = modp_set(k[v], p);
          }
          modp_NTT2(t1, gm, logn, p, p0i);
          for (v = 0, y = f, x = fk + u;
                  v < n; v ++, y += fstride, x += tlen) {
              *x = zint_mod_small_signed(y, flen, p, p0i, R2, Rx);
          }
          modp_NTT2_ext(fk + u, tlen, gm, logn, p, p0i);
          for (v = 0, x = fk + u; v < n; v ++, x += tlen) {
              *x = modp_montymul(
                       modp_montymul(t1[v], *x, p, p0i), R2, p, p0i);
          }
          modp_iNTT2_ext(fk + u, tlen, igm, logn, p, p0i);
      }
  
      /*
       * Rebuild k*f.
       */
      zint_rebuild_CRT(fk, tlen, tlen, n, primes, 1, t1);
  
      /*
       * Subtract k*f, scaled, from F.
       */
      for (u = 0, x = F, y = fk; u < n; u ++, x += Fstride, y += tlen) {
          zint_sub_scaled(x, Flen, y, tlen, sch, scl);
      }
  }
  
  /* ==================================================================== */
  
  
  #define RNG_CONTEXT   inner_shake256_context
  
  /*
   * Get a random 8-byte integer from a SHAKE-based RNG. This function
   * ensures consistent interpretation of the SHAKE output so that
   * the same values will be obtained over different platforms, in case
   * a known seed is used.
   */
  static inline uint64_t
  get_rng_u64(inner_shake256_context *rng) {
      /*
       * We enforce little-endian representation.
       */
  
      uint8_t tmp[8];
  
      inner_shake256_extract(rng, tmp, sizeof tmp);
      return (uint64_t)tmp[0]
             | ((uint64_t)tmp[1] << 8)
             | ((uint64_t)tmp[2] << 16)
             | ((uint64_t)tmp[3] << 24)
             | ((uint64_t)tmp[4] << 32)
             | ((uint64_t)tmp[5] << 40)
             | ((uint64_t)tmp[6] << 48)
             | ((uint64_t)tmp[7] << 56);
  }
  
  /*
   * Table below incarnates a discrete Gaussian distribution:
   *    D(x) = exp(-(x^2)/(2*sigma^2))
   * where sigma = 1.17*sqrt(q/(2*N)), q = 12289, and N = 1024.
   * Element 0 of the table is P(x = 0).
   * For k > 0, element k is P(x >= k+1 | x > 0).
   * Probabilities are scaled up by 2^63.
   */
  static const uint64_t gauss_1024_12289[] = {
      1283868770400643928u,  6416574995475331444u,  4078260278032692663u,
      2353523259288686585u,  1227179971273316331u,   575931623374121527u,
      242543240509105209u,    91437049221049666u,    30799446349977173u,
      9255276791179340u,     2478152334826140u,      590642893610164u,
      125206034929641u,       23590435911403u,        3948334035941u,
      586753615614u,          77391054539u,           9056793210u,
      940121950u,             86539696u,              7062824u,
      510971u,                32764u,                 1862u,
      94u,                    4u,                    0u
  };
  
  /*
   * Generate a random value with a Gaussian distribution centered on 0.
   * The RNG must be ready for extraction (already flipped).
   *
   * Distribution has standard deviation 1.17*sqrt(q/(2*N)). The
   * precomputed table is for N = 1024. Since the sum of two independent
   * values of standard deviation sigma has standard deviation
   * sigma*sqrt(2), then we can just generate more values and add them
   * together for lower dimensions.
   */
  static int
  mkgauss(RNG_CONTEXT *rng, unsigned logn) {
      unsigned u, g;
      int val;
  
      g = 1U << (10 - logn);
      val = 0;
      for (u = 0; u < g; u ++) {
          /*
           * Each iteration generates one value with the
           * Gaussian distribution for N = 1024.
           *
           * We use two random 64-bit values. First value
           * decides on whether the generated value is 0, and,
           * if not, the sign of the value. Second random 64-bit
           * word is used to generate the non-zero value.
           *
           * For constant-time code we have to read the complete
           * table. This has negligible cost, compared with the
           * remainder of the keygen process (solving the NTRU
           * equation).
           */
          uint64_t r;
          uint32_t f, v, k, neg;
  
          /*
           * First value:
           *  - flag 'neg' is randomly selected to be 0 or 1.
           *  - flag 'f' is set to 1 if the generated value is zero,
           *    or set to 0 otherwise.
           */
          r = get_rng_u64(rng);
          neg = (uint32_t)(r >> 63);
          r &= ~((uint64_t)1 << 63);
          f = (uint32_t)((r - gauss_1024_12289[0]) >> 63);
  
          /*
           * We produce a new random 63-bit integer r, and go over
           * the array, starting at index 1. We store in v the
           * index of the first array element which is not greater
           * than r, unless the flag f was already 1.
           */
          v = 0;
          r = get_rng_u64(rng);
          r &= ~((uint64_t)1 << 63);
          for (k = 1; k < (uint32_t)((sizeof gauss_1024_12289)
                                     / (sizeof gauss_1024_12289[0])); k ++) {
              uint32_t t;
  
              t = (uint32_t)((r - gauss_1024_12289[k]) >> 63) ^ 1;
              v |= k & -(t & (f ^ 1));
              f |= t;
          }
  
          /*
           * We apply the sign ('neg' flag). If the value is zero,
           * the sign has no effect.
           */
          v = (v ^ -neg) + neg;
  
          /*
           * Generated value is added to val.
           */
          val += *(int32_t *)&v;
      }
      return val;
  }
  
  /*
   * The MAX_BL_SMALL[] and MAX_BL_LARGE[] contain the lengths, in 31-bit
   * words, of intermediate values in the computation:
   *
   *   MAX_BL_SMALL[depth]: length for the input f and g at that depth
   *   MAX_BL_LARGE[depth]: length for the unreduced F and G at that depth
   *
   * Rules:
   *
   *  - Within an array, values grow.
   *
   *  - The 'SMALL' array must have an entry for maximum depth, corresponding
   *    to the size of values used in the binary GCD. There is no such value
   *    for the 'LARGE' array (the binary GCD yields already reduced
   *    coefficients).
   *
   *  - MAX_BL_LARGE[depth] >= MAX_BL_SMALL[depth + 1].
   *
   *  - Values must be large enough to handle the common cases, with some
   *    margins.
   *
   *  - Values must not be "too large" either because we will convert some
   *    integers into floating-point values by considering the top 10 words,
   *    i.e. 310 bits; hence, for values of length more than 10 words, we
   *    should take care to have the length centered on the expected size.
   *
   * The following average lengths, in bits, have been measured on thousands
   * of random keys (fg = max length of the absolute value of coefficients
   * of f and g at that depth; FG = idem for the unreduced F and G; for the
   * maximum depth, F and G are the output of binary GCD, multiplied by q;
   * for each value, the average and standard deviation are provided).
   *
   * Binary case:
   *    depth: 10    fg: 6307.52 (24.48)    FG: 6319.66 (24.51)
   *    depth:  9    fg: 3138.35 (12.25)    FG: 9403.29 (27.55)
   *    depth:  8    fg: 1576.87 ( 7.49)    FG: 4703.30 (14.77)
   *    depth:  7    fg:  794.17 ( 4.98)    FG: 2361.84 ( 9.31)
   *    depth:  6    fg:  400.67 ( 3.10)    FG: 1188.68 ( 6.04)
   *    depth:  5    fg:  202.22 ( 1.87)    FG:  599.81 ( 3.87)
   *    depth:  4    fg:  101.62 ( 1.02)    FG:  303.49 ( 2.38)
   *    depth:  3    fg:   50.37 ( 0.53)    FG:  153.65 ( 1.39)
   *    depth:  2    fg:   24.07 ( 0.25)    FG:   78.20 ( 0.73)
   *    depth:  1    fg:   10.99 ( 0.08)    FG:   39.82 ( 0.41)
   *    depth:  0    fg:    4.00 ( 0.00)    FG:   19.61 ( 0.49)
   *
   * Integers are actually represented either in binary notation over
   * 31-bit words (signed, using two's complement), or in RNS, modulo
   * many small primes. These small primes are close to, but slightly
   * lower than, 2^31. Use of RNS loses less than two bits, even for
   * the largest values.
   *
   * IMPORTANT: if these values are modified, then the temporary buffer
   * sizes (FALCON_KEYGEN_TEMP_*, in inner.h) must be recomputed
   * accordingly.
   */
  
  static const size_t MAX_BL_SMALL[] = {
      1, 1, 2, 2, 4, 7, 14, 27, 53, 106, 209
  };
  
  static const size_t MAX_BL_LARGE[] = {
      2, 2, 5, 7, 12, 21, 40, 78, 157, 308
  };
  
  /*
   * Average and standard deviation for the maximum size (in bits) of
   * coefficients of (f,g), depending on depth. These values are used
   * to compute bounds for Babai's reduction.
   */
  static const struct {
      int avg;
      int std;
  } BITLENGTH[] = {
      {    4,  0 },
      {   11,  1 },
      {   24,  1 },
      {   50,  1 },
      {  102,  1 },
      {  202,  2 },
      {  401,  4 },
      {  794,  5 },
      { 1577,  8 },
      { 3138, 13 },
      { 6308, 25 }
  };
  
  /*
   * Minimal recursion depth at which we rebuild intermediate values
   * when reconstructing f and g.
   */
  #define DEPTH_INT_FG   4
  
  /*
   * Compute squared norm of a short vector. Returned value is saturated to
   * 2^32-1 if it is not lower than 2^31.
   */
  static uint32_t
  poly_small_sqnorm(const int8_t *f, unsigned logn) {
      size_t n, u;
      uint32_t s, ng;
  
      n = MKN(logn);
      s = 0;
      ng = 0;
      for (u = 0; u < n; u ++) {
          int32_t z;
  
          z = f[u];
          s += (uint32_t)(z * z);
          ng |= s;
      }
      return s | -(ng >> 31);
  }
  
  /*
   * Align (upwards) the provided 'data' pointer with regards to 'base'
   * so that the offset is a multiple of the size of 'fpr'.
   */
  static fpr *
  align_fpr(void *base, void *data) {
      uint8_t *cb, *cd;
      size_t k, km;
  
      cb = base;
      cd = data;
      k = (size_t)(cd - cb);
      km = k % sizeof(fpr);
      if (km) {
          k += (sizeof(fpr)) - km;
      }
      return (fpr *)(cb + k);
  }
  
  /*
   * Align (upwards) the provided 'data' pointer with regards to 'base'
   * so that the offset is a multiple of the size of 'uint32_t'.
   */
  static uint32_t *
  align_u32(void *base, void *data) {
      uint8_t *cb, *cd;
      size_t k, km;
  
      cb = base;
      cd = data;
      k = (size_t)(cd - cb);
      km = k % sizeof(uint32_t);
      if (km) {
          k += (sizeof(uint32_t)) - km;
      }
      return (uint32_t *)(cb + k);
  }
  
  /*
   * Convert a small vector to floating point.
   */
  static void
  poly_small_to_fp(fpr *x, const int8_t *f, unsigned logn) {
      size_t n, u;
  
      n = MKN(logn);
      for (u = 0; u < n; u ++) {
          x[u] = fpr_of(f[u]);
      }
  }
  
  /*
   * Input: f,g of degree N = 2^logn; 'depth' is used only to get their
   * individual length.
   *
   * Output: f',g' of degree N/2, with the length for 'depth+1'.
   *
   * Values are in RNS; input and/or output may also be in NTT.
   */
  static void
  make_fg_step(uint32_t *data, unsigned logn, unsigned depth,
               int in_ntt, int out_ntt) {
      size_t n, hn, u;
      size_t slen, tlen;
      uint32_t *fd, *gd, *fs, *gs, *gm, *igm, *t1;
      const small_prime *primes;
  
      n = (size_t)1 << logn;
      hn = n >> 1;
      slen = MAX_BL_SMALL[depth];
      tlen = MAX_BL_SMALL[depth + 1];
      primes = PRIMES;
  
      /*
       * Prepare room for the result.
       */
      fd = data;
      gd = fd + hn * tlen;
      fs = gd + hn * tlen;
      gs = fs + n * slen;
      gm = gs + n * slen;
      igm = gm + n;
      t1 = igm + n;
      memmove(fs, data, 2 * n * slen * sizeof * data);
  
      /*
       * First slen words: we use the input values directly, and apply
       * inverse NTT as we go.
       */
      for (u = 0; u < slen; u ++) {
          uint32_t p, p0i, R2;
          size_t v;
          uint32_t *x;
  
          p = primes[u].p;
          p0i = modp_ninv31(p);
          R2 = modp_R2(p, p0i);
          modp_mkgm2(gm, igm, logn, primes[u].g, p, p0i);
  
          for (v = 0, x = fs + u; v < n; v ++, x += slen) {
              t1[v] = *x;
          }
          if (!in_ntt) {
              modp_NTT2(t1, gm, logn, p, p0i);
          }
          for (v = 0, x = fd + u; v < hn; v ++, x += tlen) {
              uint32_t w0, w1;
  
              w0 = t1[(v << 1) + 0];
              w1 = t1[(v << 1) + 1];
              *x = modp_montymul(
                       modp_montymul(w0, w1, p, p0i), R2, p, p0i);
          }
          if (in_ntt) {
              modp_iNTT2_ext(fs + u, slen, igm, logn, p, p0i);
          }
  
          for (v = 0, x = gs + u; v < n; v ++, x += slen) {
              t1[v] = *x;
          }
          if (!in_ntt) {
              modp_NTT2(t1, gm, logn, p, p0i);
          }
          for (v = 0, x = gd + u; v < hn; v ++, x += tlen) {
              uint32_t w0, w1;
  
              w0 = t1[(v << 1) + 0];
              w1 = t1[(v << 1) + 1];
              *x = modp_montymul(
                       modp_montymul(w0, w1, p, p0i), R2, p, p0i);
          }
          if (in_ntt) {
              modp_iNTT2_ext(gs + u, slen, igm, logn, p, p0i);
          }
  
          if (!out_ntt) {
              modp_iNTT2_ext(fd + u, tlen, igm, logn - 1, p, p0i);
              modp_iNTT2_ext(gd + u, tlen, igm, logn - 1, p, p0i);
          }
      }
  
      /*
       * Since the fs and gs words have been de-NTTized, we can use the
       * CRT to rebuild the values.
       */
      zint_rebuild_CRT(fs, slen, slen, n, primes, 1, gm);
      zint_rebuild_CRT(gs, slen, slen, n, primes, 1, gm);
  
      /*
       * Remaining words: use modular reductions to extract the values.
       */
      for (u = slen; u < tlen; u ++) {
          uint32_t p, p0i, R2, Rx;
          size_t v;
          uint32_t *x;
  
          p = primes[u].p;
          p0i = modp_ninv31(p);
          R2 = modp_R2(p, p0i);
          Rx = modp_Rx((unsigned)slen, p, p0i, R2);
          modp_mkgm2(gm, igm, logn, primes[u].g, p, p0i);
          for (v = 0, x = fs; v < n; v ++, x += slen) {
              t1[v] = zint_mod_small_signed(x, slen, p, p0i, R2, Rx);
          }
          modp_NTT2(t1, gm, logn, p, p0i);
          for (v = 0, x = fd + u; v < hn; v ++, x += tlen) {
              uint32_t w0, w1;
  
              w0 = t1[(v << 1) + 0];
              w1 = t1[(v << 1) + 1];
              *x = modp_montymul(
                       modp_montymul(w0, w1, p, p0i), R2, p, p0i);
          }
          for (v = 0, x = gs; v < n; v ++, x += slen) {
              t1[v] = zint_mod_small_signed(x, slen, p, p0i, R2, Rx);
          }
          modp_NTT2(t1, gm, logn, p, p0i);
          for (v = 0, x = gd + u; v < hn; v ++, x += tlen) {
              uint32_t w0, w1;
  
              w0 = t1[(v << 1) + 0];
              w1 = t1[(v << 1) + 1];
              *x = modp_montymul(
                       modp_montymul(w0, w1, p, p0i), R2, p, p0i);
          }
  
          if (!out_ntt) {
              modp_iNTT2_ext(fd + u, tlen, igm, logn - 1, p, p0i);
              modp_iNTT2_ext(gd + u, tlen, igm, logn - 1, p, p0i);
          }
      }
  }
  
  /*
   * Compute f and g at a specific depth, in RNS notation.
   *
   * Returned values are stored in the data[] array, at slen words per integer.
   *
   * Conditions:
   *   0 <= depth <= logn
   *
   * Space use in data[]: enough room for any two successive values (f', g',
   * f and g).
   */
  static void
  make_fg(uint32_t *data, const int8_t *f, const int8_t *g,
          unsigned logn, unsigned depth, int out_ntt) {
      size_t n, u;
      uint32_t *ft, *gt, p0;
      unsigned d;
      const small_prime *primes;
  
      n = MKN(logn);
      ft = data;
      gt = ft + n;
      primes = PRIMES;
      p0 = primes[0].p;
      for (u = 0; u < n; u ++) {
          ft[u] = modp_set(f[u], p0);
          gt[u] = modp_set(g[u], p0);
      }
  
      if (depth == 0 && out_ntt) {
          uint32_t *gm, *igm;
          uint32_t p, p0i;
  
          p = primes[0].p;
          p0i = modp_ninv31(p);
          gm = gt + n;
          igm = gm + MKN(logn);
          modp_mkgm2(gm, igm, logn, primes[0].g, p, p0i);
          modp_NTT2(ft, gm, logn, p, p0i);
          modp_NTT2(gt, gm, logn, p, p0i);
          return;
      }
  
      if (depth == 0) {
          return;
      }
      if (depth == 1) {
          make_fg_step(data, logn, 0, 0, out_ntt);
          return;
      }
      make_fg_step(data, logn, 0, 0, 1);
      for (d = 1; d + 1 < depth; d ++) {
          make_fg_step(data, logn - d, d, 1, 1);
      }
      make_fg_step(data, logn - depth + 1, depth - 1, 1, out_ntt);
  }
  
  /*
   * Solving the NTRU equation, deepest level: compute the resultants of
   * f and g with X^N+1, and use binary GCD. The F and G values are
   * returned in tmp[].
   *
   * Returned value: 1 on success, 0 on error.
   */
  static int
  solve_NTRU_deepest(unsigned logn_top,
                     const int8_t *f, const int8_t *g, uint32_t *tmp) {
      size_t len;
      uint32_t *Fp, *Gp, *fp, *gp, *t1, q;
      const small_prime *primes;
  
      len = MAX_BL_SMALL[logn_top];
      primes = PRIMES;
  
      Fp = tmp;
      Gp = Fp + len;
      fp = Gp + len;
      gp = fp + len;
      t1 = gp + len;
  
      make_fg(fp, f, g, logn_top, logn_top, 0);
  
      /*
       * We use the CRT to rebuild the resultants as big integers.
       * There are two such big integers. The resultants are always
       * nonnegative.
       */
      zint_rebuild_CRT(fp, len, len, 2, primes, 0, t1);
  
      /*
       * Apply the binary GCD. The zint_bezout() function works only
       * if both inputs are odd.
       *
       * We can test on the result and return 0 because that would
       * imply failure of the NTRU solving equation, and the (f,g)
       * values will be abandoned in that case.
       */
      if (!zint_bezout(Gp, Fp, fp, gp, len, t1)) {
          return 0;
      }
  
      /*
       * Multiply the two values by the target value q. Values must
       * fit in the destination arrays.
       * We can again test on the returned words: a non-zero output
       * of zint_mul_small() means that we exceeded our array
       * capacity, and that implies failure and rejection of (f,g).
       */
      q = 12289;
      if (zint_mul_small(Fp, len, q) != 0
              || zint_mul_small(Gp, len, q) != 0) {
          return 0;
      }
  
      return 1;
  }
  
  /*
   * Solving the NTRU equation, intermediate level. Upon entry, the F and G
   * from the previous level should be in the tmp[] array.
   * This function MAY be invoked for the top-level (in which case depth = 0).
   *
   * Returned value: 1 on success, 0 on error.
   */
  static int
  solve_NTRU_intermediate(unsigned logn_top,
                          const int8_t *f, const int8_t *g, unsigned depth, uint32_t *tmp) {
      /*
       * In this function, 'logn' is the log2 of the degree for
       * this step. If N = 2^logn, then:
       *  - the F and G values already in fk->tmp (from the deeper
       *    levels) have degree N/2;
       *  - this function should return F and G of degree N.
       */
      unsigned logn;
      size_t n, hn, slen, dlen, llen, rlen, FGlen, u;
      uint32_t *Fd, *Gd, *Ft, *Gt, *ft, *gt, *t1;
      fpr *rt1, *rt2, *rt3, *rt4, *rt5;
      int scale_fg, minbl_fg, maxbl_fg, maxbl_FG, scale_k;
      uint32_t *x, *y;
      int32_t *k;
      const small_prime *primes;
  
      logn = logn_top - depth;
      n = (size_t)1 << logn;
      hn = n >> 1;
  
      /*
       * slen = size for our input f and g; also size of the reduced
       *        F and G we return (degree N)
       *
       * dlen = size of the F and G obtained from the deeper level
       *        (degree N/2 or N/3)
       *
       * llen = size for intermediary F and G before reduction (degree N)
       *
       * We build our non-reduced F and G as two independent halves each,
       * of degree N/2 (F = F0 + X*F1, G = G0 + X*G1).
       */
      slen = MAX_BL_SMALL[depth];
      dlen = MAX_BL_SMALL[depth + 1];
      llen = MAX_BL_LARGE[depth];
      primes = PRIMES;
  
      /*
       * Fd and Gd are the F and G from the deeper level.
       */
      Fd = tmp;
      Gd = Fd + dlen * hn;
  
      /*
       * Compute the input f and g for this level. Note that we get f
       * and g in RNS + NTT representation.
       */
      ft = Gd + dlen * hn;
      make_fg(ft, f, g, logn_top, depth, 1);
  
      /*
       * Move the newly computed f and g to make room for our candidate
       * F and G (unreduced).
       */
      Ft = tmp;
      Gt = Ft + n * llen;
      t1 = Gt + n * llen;
      memmove(t1, ft, 2 * n * slen * sizeof * ft);
      ft = t1;
      gt = ft + slen * n;
      t1 = gt + slen * n;
  
      /*
       * Move Fd and Gd _after_ f and g.
       */
      memmove(t1, Fd, 2 * hn * dlen * sizeof * Fd);
      Fd = t1;
      Gd = Fd + hn * dlen;
  
      /*
       * We reduce Fd and Gd modulo all the small primes we will need,
       * and store the values in Ft and Gt (only n/2 values in each).
       */
      for (u = 0; u < llen; u ++) {
          uint32_t p, p0i, R2, Rx;
          size_t v;
          uint32_t *xs, *ys, *xd, *yd;
  
          p = primes[u].p;
          p0i = modp_ninv31(p);
          R2 = modp_R2(p, p0i);
          Rx = modp_Rx((unsigned)dlen, p, p0i, R2);
          for (v = 0, xs = Fd, ys = Gd, xd = Ft + u, yd = Gt + u;
                  v < hn;
                  v ++, xs += dlen, ys += dlen, xd += llen, yd += llen) {
              *xd = zint_mod_small_signed(xs, dlen, p, p0i, R2, Rx);
              *yd = zint_mod_small_signed(ys, dlen, p, p0i, R2, Rx);
          }
      }
  
      /*
       * We do not need Fd and Gd after that point.
       */
  
      /*
       * Compute our F and G modulo sufficiently many small primes.
       */
      for (u = 0; u < llen; u ++) {
          uint32_t p, p0i, R2;
          uint32_t *gm, *igm, *fx, *gx, *Fp, *Gp;
          size_t v;
  
          /*
           * All computations are done modulo p.
           */
          p = primes[u].p;
          p0i = modp_ninv31(p);
          R2 = modp_R2(p, p0i);
  
          /*
           * If we processed slen words, then f and g have been
           * de-NTTized, and are in RNS; we can rebuild them.
           */
          if (u == slen) {
              zint_rebuild_CRT(ft, slen, slen, n, primes, 1, t1);
              zint_rebuild_CRT(gt, slen, slen, n, primes, 1, t1);
          }
  
          gm = t1;
          igm = gm + n;
          fx = igm + n;
          gx = fx + n;
  
          modp_mkgm2(gm, igm, logn, primes[u].g, p, p0i);
  
          if (u < slen) {
              for (v = 0, x = ft + u, y = gt + u;
                      v < n; v ++, x += slen, y += slen) {
                  fx[v] = *x;
                  gx[v] = *y;
              }
              modp_iNTT2_ext(ft + u, slen, igm, logn, p, p0i);
              modp_iNTT2_ext(gt + u, slen, igm, logn, p, p0i);
          } else {
              uint32_t Rx;
  
              Rx = modp_Rx((unsigned)slen, p, p0i, R2);
              for (v = 0, x = ft, y = gt;
                      v < n; v ++, x += slen, y += slen) {
                  fx[v] = zint_mod_small_signed(x, slen,
                                                p, p0i, R2, Rx);
                  gx[v] = zint_mod_small_signed(y, slen,
                                                p, p0i, R2, Rx);
              }
              modp_NTT2(fx, gm, logn, p, p0i);
              modp_NTT2(gx, gm, logn, p, p0i);
          }
  
          /*
           * Get F' and G' modulo p and in NTT representation
           * (they have degree n/2). These values were computed in
           * a previous step, and stored in Ft and Gt.
           */
          Fp = gx + n;
          Gp = Fp + hn;
          for (v = 0, x = Ft + u, y = Gt + u;
                  v < hn; v ++, x += llen, y += llen) {
              Fp[v] = *x;
              Gp[v] = *y;
          }
          modp_NTT2(Fp, gm, logn - 1, p, p0i);
          modp_NTT2(Gp, gm, logn - 1, p, p0i);
  
          /*
           * Compute our F and G modulo p.
           *
           * General case:
           *
           *   we divide degree by d = 2 or 3
           *   f'(x^d) = N(f)(x^d) = f * adj(f)
           *   g'(x^d) = N(g)(x^d) = g * adj(g)
           *   f'*G' - g'*F' = q
           *   F = F'(x^d) * adj(g)
           *   G = G'(x^d) * adj(f)
           *
           * We compute things in the NTT. We group roots of phi
           * such that all roots x in a group share the same x^d.
           * If the roots in a group are x_1, x_2... x_d, then:
           *
           *   N(f)(x_1^d) = f(x_1)*f(x_2)*...*f(x_d)
           *
           * Thus, we have:
           *
           *   G(x_1) = f(x_2)*f(x_3)*...*f(x_d)*G'(x_1^d)
           *   G(x_2) = f(x_1)*f(x_3)*...*f(x_d)*G'(x_1^d)
           *   ...
           *   G(x_d) = f(x_1)*f(x_2)*...*f(x_{d-1})*G'(x_1^d)
           *
           * In all cases, we can thus compute F and G in NTT
           * representation by a few simple multiplications.
           * Moreover, in our chosen NTT representation, roots
           * from the same group are consecutive in RAM.
           */
          for (v = 0, x = Ft + u, y = Gt + u; v < hn;
                  v ++, x += (llen << 1), y += (llen << 1)) {
              uint32_t ftA, ftB, gtA, gtB;
              uint32_t mFp, mGp;
  
              ftA = fx[(v << 1) + 0];
              ftB = fx[(v << 1) + 1];
              gtA = gx[(v << 1) + 0];
              gtB = gx[(v << 1) + 1];
              mFp = modp_montymul(Fp[v], R2, p, p0i);
              mGp = modp_montymul(Gp[v], R2, p, p0i);
              x[0] = modp_montymul(gtB, mFp, p, p0i);
              x[llen] = modp_montymul(gtA, mFp, p, p0i);
              y[0] = modp_montymul(ftB, mGp, p, p0i);
              y[llen] = modp_montymul(ftA, mGp, p, p0i);
          }
          modp_iNTT2_ext(Ft + u, llen, igm, logn, p, p0i);
          modp_iNTT2_ext(Gt + u, llen, igm, logn, p, p0i);
      }
  
      /*
       * Rebuild F and G with the CRT.
       */
      zint_rebuild_CRT(Ft, llen, llen, n, primes, 1, t1);
      zint_rebuild_CRT(Gt, llen, llen, n, primes, 1, t1);
  
      /*
       * At that point, Ft, Gt, ft and gt are consecutive in RAM (in that
       * order).
       */
  
      /*
       * Apply Babai reduction to bring back F and G to size slen.
       *
       * We use the FFT to compute successive approximations of the
       * reduction coefficient. We first isolate the top bits of
       * the coefficients of f and g, and convert them to floating
       * point; with the FFT, we compute adj(f), adj(g), and
       * 1/(f*adj(f)+g*adj(g)).
       *
       * Then, we repeatedly apply the following:
       *
       *   - Get the top bits of the coefficients of F and G into
       *     floating point, and use the FFT to compute:
       *        (F*adj(f)+G*adj(g))/(f*adj(f)+g*adj(g))
       *
       *   - Convert back that value into normal representation, and
       *     round it to the nearest integers, yielding a polynomial k.
       *     Proper scaling is applied to f, g, F and G so that the
       *     coefficients fit on 32 bits (signed).
       *
       *   - Subtract k*f from F and k*g from G.
       *
       * Under normal conditions, this process reduces the size of F
       * and G by some bits at each iteration. For constant-time
       * operation, we do not want to measure the actual length of
       * F and G; instead, we do the following:
       *
       *   - f and g are converted to floating-point, with some scaling
       *     if necessary to keep values in the representable range.
       *
       *   - For each iteration, we _assume_ a maximum size for F and G,
       *     and use the values at that size. If we overreach, then
       *     we get zeros, which is harmless: the resulting coefficients
       *     of k will be 0 and the value won't be reduced.
       *
       *   - We conservatively assume that F and G will be reduced by
       *     at least 25 bits at each iteration.
       *
       * Even when reaching the bottom of the reduction, reduction
       * coefficient will remain low. If it goes out-of-range, then
       * something wrong occurred and the whole NTRU solving fails.
       */
  
      /*
       * Memory layout:
       *  - We need to compute and keep adj(f), adj(g), and
       *    1/(f*adj(f)+g*adj(g)) (sizes N, N and N/2 fp numbers,
       *    respectively).
       *  - At each iteration we need two extra fp buffer (N fp values),
       *    and produce a k (N 32-bit words). k will be shared with one
       *    of the fp buffers.
       *  - To compute k*f and k*g efficiently (with the NTT), we need
       *    some extra room; we reuse the space of the temporary buffers.
       *
       * Arrays of 'fpr' are obtained from the temporary array itself.
       * We ensure that the base is at a properly aligned offset (the
       * source array tmp[] is supposed to be already aligned).
       */
  
      rt3 = align_fpr(tmp, t1);
      rt4 = rt3 + n;
      rt5 = rt4 + n;
      rt1 = rt5 + (n >> 1);
      k = (int32_t *)align_u32(tmp, rt1);
      rt2 = align_fpr(tmp, k + n);
      if (rt2 < (rt1 + n)) {
          rt2 = rt1 + n;
      }
      t1 = (uint32_t *)k + n;
  
      /*
       * Get f and g into rt3 and rt4 as floating-point approximations.
       *
       * We need to "scale down" the floating-point representation of
       * coefficients when they are too big. We want to keep the value
       * below 2^310 or so. Thus, when values are larger than 10 words,
       * we consider only the top 10 words. Array lengths have been
       * computed so that average maximum length will fall in the
       * middle or the upper half of these top 10 words.
       */
      rlen = slen;
      if (rlen > 10) {
          rlen = 10;
      }
      poly_big_to_fp(rt3, ft + slen - rlen, rlen, slen, logn);
      poly_big_to_fp(rt4, gt + slen - rlen, rlen, slen, logn);
  
      /*
       * Values in rt3 and rt4 are downscaled by 2^(scale_fg).
       */
      scale_fg = 31 * (int)(slen - rlen);
  
      /*
       * Estimated boundaries for the maximum size (in bits) of the
       * coefficients of (f,g). We use the measured average, and
       * allow for a deviation of at most six times the standard
       * deviation.
       */
      minbl_fg = BITLENGTH[depth].avg - 6 * BITLENGTH[depth].std;
      maxbl_fg = BITLENGTH[depth].avg + 6 * BITLENGTH[depth].std;
  
      /*
       * Compute 1/(f*adj(f)+g*adj(g)) in rt5. We also keep adj(f)
       * and adj(g) in rt3 and rt4, respectively.
       */
      FFT(rt3, logn);
      FFT(rt4, logn);
      poly_invnorm2_fft(rt5, rt3, rt4, logn);
      poly_adj_fft(rt3, logn);
      poly_adj_fft(rt4, logn);
  
      /*
       * Reduce F and G repeatedly.
       *
       * The expected maximum bit length of coefficients of F and G
       * is kept in maxbl_FG, with the corresponding word length in
       * FGlen.
       */
      FGlen = llen;
      maxbl_FG = 31 * (int)llen;
  
      /*
       * Each reduction operation computes the reduction polynomial
       * "k". We need that polynomial to have coefficients that fit
       * on 32-bit signed integers, with some scaling; thus, we use
       * a descending sequence of scaling values, down to zero.
       *
       * The size of the coefficients of k is (roughly) the difference
       * between the size of the coefficients of (F,G) and the size
       * of the coefficients of (f,g). Thus, the maximum size of the
       * coefficients of k is, at the start, maxbl_FG - minbl_fg;
       * this is our starting scale value for k.
       *
       * We need to estimate the size of (F,G) during the execution of
       * the algorithm; we are allowed some overestimation but not too
       * much (poly_big_to_fp() uses a 310-bit window). Generally
       * speaking, after applying a reduction with k scaled to
       * scale_k, the size of (F,G) will be size(f,g) + scale_k + dd,
       * where 'dd' is a few bits to account for the fact that the
       * reduction is never perfect (intuitively, dd is on the order
       * of sqrt(N), so at most 5 bits; we here allow for 10 extra
       * bits).
       *
       * The size of (f,g) is not known exactly, but maxbl_fg is an
       * upper bound.
       */
      scale_k = maxbl_FG - minbl_fg;
  
      for (;;) {
          int scale_FG, dc, new_maxbl_FG;
          uint32_t scl, sch;
          fpr pdc, pt;
  
          /*
           * Convert current F and G into floating-point. We apply
           * scaling if the current length is more than 10 words.
           */
          rlen = FGlen;
          if (rlen > 10) {
              rlen = 10;
          }
          scale_FG = 31 * (int)(FGlen - rlen);
          poly_big_to_fp(rt1, Ft + FGlen - rlen, rlen, llen, logn);
          poly_big_to_fp(rt2, Gt + FGlen - rlen, rlen, llen, logn);
  
          /*
           * Compute (F*adj(f)+G*adj(g))/(f*adj(f)+g*adj(g)) in rt2.
           */
          FFT(rt1, logn);
          FFT(rt2, logn);
          poly_mul_fft(rt1, rt3, logn);
          poly_mul_fft(rt2, rt4, logn);
          poly_add(rt2, rt1, logn);
          poly_mul_autoadj_fft(rt2, rt5, logn);
          iFFT(rt2, logn);
  
          /*
           * (f,g) are scaled by 'scale_fg', meaning that the
           * numbers in rt3/rt4 should be multiplied by 2^(scale_fg)
           * to have their true mathematical value.
           *
           * (F,G) are similarly scaled by 'scale_FG'. Therefore,
           * the value we computed in rt2 is scaled by
           * 'scale_FG-scale_fg'.
           *
           * We want that value to be scaled by 'scale_k', hence we
           * apply a corrective scaling. After scaling, the values
           * should fit in -2^31-1..+2^31-1.
           */
          dc = scale_k - scale_FG + scale_fg;
  
          /*
           * We will need to multiply values by 2^(-dc). The value
           * 'dc' is not secret, so we can compute 2^(-dc) with a
           * non-constant-time process.
           * (We could use ldexp(), but we prefer to avoid any
           * dependency on libm. When using FP emulation, we could
           * use our fpr_ldexp(), which is constant-time.)
           */
          if (dc < 0) {
              dc = -dc;
              pt = fpr_two;
          } else {
              pt = fpr_onehalf;
          }
          pdc = fpr_one;
          while (dc != 0) {
              if ((dc & 1) != 0) {
                  pdc = fpr_mul(pdc, pt);
              }
              dc >>= 1;
              pt = fpr_sqr(pt);
          }
  
          for (u = 0; u < n; u ++) {
              fpr xv;
  
              xv = fpr_mul(rt2[u], pdc);
  
              /*
               * Sometimes the values can be out-of-bounds if
               * the algorithm fails; we must not call
               * fpr_rint() (and cast to int32_t) if the value
               * is not in-bounds. Note that the test does not
               * break constant-time discipline, since any
               * failure here implies that we discard the current
               * secret key (f,g).
               */
              if (!fpr_lt(fpr_mtwo31m1, xv)
                      || !fpr_lt(xv, fpr_ptwo31m1)) {
                  return 0;
              }
              k[u] = (int32_t)fpr_rint(xv);
          }
  
          /*
           * Values in k[] are integers. They really are scaled
           * down by maxbl_FG - minbl_fg bits.
           *
           * If we are at low depth, then we use the NTT to
           * compute k*f and k*g.
           */
          sch = (uint32_t)(scale_k / 31);
          scl = (uint32_t)(scale_k % 31);
          if (depth <= DEPTH_INT_FG) {
              poly_sub_scaled_ntt(Ft, FGlen, llen, ft, slen, slen,
                                  k, sch, scl, logn, t1);
              poly_sub_scaled_ntt(Gt, FGlen, llen, gt, slen, slen,
                                  k, sch, scl, logn, t1);
          } else {
              poly_sub_scaled(Ft, FGlen, llen, ft, slen, slen,
                              k, sch, scl, logn);
              poly_sub_scaled(Gt, FGlen, llen, gt, slen, slen,
                              k, sch, scl, logn);
          }
  
          /*
           * We compute the new maximum size of (F,G), assuming that
           * (f,g) has _maximal_ length (i.e. that reduction is
           * "late" instead of "early". We also adjust FGlen
           * accordingly.
           */
          new_maxbl_FG = scale_k + maxbl_fg + 10;
          if (new_maxbl_FG < maxbl_FG) {
              maxbl_FG = new_maxbl_FG;
              if ((int)FGlen * 31 >= maxbl_FG + 31) {
                  FGlen --;
              }
          }
  
          /*
           * We suppose that scaling down achieves a reduction by
           * at least 25 bits per iteration. We stop when we have
           * done the loop with an unscaled k.
           */
          if (scale_k <= 0) {
              break;
          }
          scale_k -= 25;
          if (scale_k < 0) {
              scale_k = 0;
          }
      }
  
      /*
       * If (F,G) length was lowered below 'slen', then we must take
       * care to re-extend the sign.
       */
      if (FGlen < slen) {
          for (u = 0; u < n; u ++, Ft += llen, Gt += llen) {
              size_t v;
              uint32_t sw;
  
              sw = -(Ft[FGlen - 1] >> 30) >> 1;
              for (v = FGlen; v < slen; v ++) {
                  Ft[v] = sw;
              }
              sw = -(Gt[FGlen - 1] >> 30) >> 1;
              for (v = FGlen; v < slen; v ++) {
                  Gt[v] = sw;
              }
          }
      }
  
      /*
       * Compress encoding of all values to 'slen' words (this is the
       * expected output format).
       */
      for (u = 0, x = tmp, y = tmp;
              u < (n << 1); u ++, x += slen, y += llen) {
          memmove(x, y, slen * sizeof * y);
      }
      return 1;
  }
  
  /*
   * Solving the NTRU equation, binary case, depth = 1. Upon entry, the
   * F and G from the previous level should be in the tmp[] array.
   *
   * Returned value: 1 on success, 0 on error.
   */
  static int
  solve_NTRU_binary_depth1(unsigned logn_top,
                           const int8_t *f, const int8_t *g, uint32_t *tmp) {
      /*
       * The first half of this function is a copy of the corresponding
       * part in solve_NTRU_intermediate(), for the reconstruction of
       * the unreduced F and G. The second half (Babai reduction) is
       * done differently, because the unreduced F and G fit in 53 bits
       * of precision, allowing a much simpler process with lower RAM
       * usage.
       */
      unsigned depth, logn;
      size_t n_top, n, hn, slen, dlen, llen, u;
      uint32_t *Fd, *Gd, *Ft, *Gt, *ft, *gt, *t1;
      fpr *rt1, *rt2, *rt3, *rt4, *rt5, *rt6;
      uint32_t *x, *y;
  
      depth = 1;
      n_top = (size_t)1 << logn_top;
      logn = logn_top - depth;
      n = (size_t)1 << logn;
      hn = n >> 1;
  
      /*
       * Equations are:
       *
       *   f' = f0^2 - X^2*f1^2
       *   g' = g0^2 - X^2*g1^2
       *   F' and G' are a solution to f'G' - g'F' = q (from deeper levels)
       *   F = F'*(g0 - X*g1)
       *   G = G'*(f0 - X*f1)
       *
       * f0, f1, g0, g1, f', g', F' and G' are all "compressed" to
       * degree N/2 (their odd-indexed coefficients are all zero).
       */
  
      /*
       * slen = size for our input f and g; also size of the reduced
       *        F and G we return (degree N)
       *
       * dlen = size of the F and G obtained from the deeper level
       *        (degree N/2)
       *
       * llen = size for intermediary F and G before reduction (degree N)
       *
       * We build our non-reduced F and G as two independent halves each,
       * of degree N/2 (F = F0 + X*F1, G = G0 + X*G1).
       */
      slen = MAX_BL_SMALL[depth];
      dlen = MAX_BL_SMALL[depth + 1];
      llen = MAX_BL_LARGE[depth];
  
      /*
       * Fd and Gd are the F and G from the deeper level. Ft and Gt
       * are the destination arrays for the unreduced F and G.
       */
      Fd = tmp;
      Gd = Fd + dlen * hn;
      Ft = Gd + dlen * hn;
      Gt = Ft + llen * n;
  
      /*
       * We reduce Fd and Gd modulo all the small primes we will need,
       * and store the values in Ft and Gt.
       */
      for (u = 0; u < llen; u ++) {
          uint32_t p, p0i, R2, Rx;
          size_t v;
          uint32_t *xs, *ys, *xd, *yd;
  
          p = PRIMES[u].p;
          p0i = modp_ninv31(p);
          R2 = modp_R2(p, p0i);
          Rx = modp_Rx((unsigned)dlen, p, p0i, R2);
          for (v = 0, xs = Fd, ys = Gd, xd = Ft + u, yd = Gt + u;
                  v < hn;
                  v ++, xs += dlen, ys += dlen, xd += llen, yd += llen) {
              *xd = zint_mod_small_signed(xs, dlen, p, p0i, R2, Rx);
              *yd = zint_mod_small_signed(ys, dlen, p, p0i, R2, Rx);
          }
      }
  
      /*
       * Now Fd and Gd are not needed anymore; we can squeeze them out.
       */
      memmove(tmp, Ft, llen * n * sizeof(uint32_t));
      Ft = tmp;
      memmove(Ft + llen * n, Gt, llen * n * sizeof(uint32_t));
      Gt = Ft + llen * n;
      ft = Gt + llen * n;
      gt = ft + slen * n;
  
      t1 = gt + slen * n;
  
      /*
       * Compute our F and G modulo sufficiently many small primes.
       */
      for (u = 0; u < llen; u ++) {
          uint32_t p, p0i, R2;
          uint32_t *gm, *igm, *fx, *gx, *Fp, *Gp;
          unsigned e;
          size_t v;
  
          /*
           * All computations are done modulo p.
           */
          p = PRIMES[u].p;
          p0i = modp_ninv31(p);
          R2 = modp_R2(p, p0i);
  
          /*
           * We recompute things from the source f and g, of full
           * degree. However, we will need only the n first elements
           * of the inverse NTT table (igm); the call to modp_mkgm()
           * below will fill n_top elements in igm[] (thus overflowing
           * into fx[]) but later code will overwrite these extra
           * elements.
           */
          gm = t1;
          igm = gm + n_top;
          fx = igm + n;
          gx = fx + n_top;
          modp_mkgm2(gm, igm, logn_top, PRIMES[u].g, p, p0i);
  
          /*
           * Set ft and gt to f and g modulo p, respectively.
           */
          for (v = 0; v < n_top; v ++) {
              fx[v] = modp_set(f[v], p);
              gx[v] = modp_set(g[v], p);
          }
  
          /*
           * Convert to NTT and compute our f and g.
           */
          modp_NTT2(fx, gm, logn_top, p, p0i);
          modp_NTT2(gx, gm, logn_top, p, p0i);
          for (e = logn_top; e > logn; e --) {
              modp_poly_rec_res(fx, e, p, p0i, R2);
              modp_poly_rec_res(gx, e, p, p0i, R2);
          }
  
          /*
           * From that point onward, we only need tables for
           * degree n, so we can save some space.
           */
          if (depth > 0) { /* always true */
              memmove(gm + n, igm, n * sizeof * igm);
              igm = gm + n;
              memmove(igm + n, fx, n * sizeof * ft);
              fx = igm + n;
              memmove(fx + n, gx, n * sizeof * gt);
              gx = fx + n;
          }
  
          /*
           * Get F' and G' modulo p and in NTT representation
           * (they have degree n/2). These values were computed
           * in a previous step, and stored in Ft and Gt.
           */
          Fp = gx + n;
          Gp = Fp + hn;
          for (v = 0, x = Ft + u, y = Gt + u;
                  v < hn; v ++, x += llen, y += llen) {
              Fp[v] = *x;
              Gp[v] = *y;
          }
          modp_NTT2(Fp, gm, logn - 1, p, p0i);
          modp_NTT2(Gp, gm, logn - 1, p, p0i);
  
          /*
           * Compute our F and G modulo p.
           *
           * Equations are:
           *
           *   f'(x^2) = N(f)(x^2) = f * adj(f)
           *   g'(x^2) = N(g)(x^2) = g * adj(g)
           *
           *   f'*G' - g'*F' = q
           *
           *   F = F'(x^2) * adj(g)
           *   G = G'(x^2) * adj(f)
           *
           * The NTT representation of f is f(w) for all w which
           * are roots of phi. In the binary case, as well as in
           * the ternary case for all depth except the deepest,
           * these roots can be grouped in pairs (w,-w), and we
           * then have:
           *
           *   f(w) = adj(f)(-w)
           *   f(-w) = adj(f)(w)
           *
           * and w^2 is then a root for phi at the half-degree.
           *
           * At the deepest level in the ternary case, this still
           * holds, in the following sense: the roots of x^2-x+1
           * are (w,-w^2) (for w^3 = -1, and w != -1), and we
           * have:
           *
           *   f(w) = adj(f)(-w^2)
           *   f(-w^2) = adj(f)(w)
           *
           * In all case, we can thus compute F and G in NTT
           * representation by a few simple multiplications.
           * Moreover, the two roots for each pair are consecutive
           * in our bit-reversal encoding.
           */
          for (v = 0, x = Ft + u, y = Gt + u;
                  v < hn; v ++, x += (llen << 1), y += (llen << 1)) {
              uint32_t ftA, ftB, gtA, gtB;
              uint32_t mFp, mGp;
  
              ftA = fx[(v << 1) + 0];
              ftB = fx[(v << 1) + 1];
              gtA = gx[(v << 1) + 0];
              gtB = gx[(v << 1) + 1];
              mFp = modp_montymul(Fp[v], R2, p, p0i);
              mGp = modp_montymul(Gp[v], R2, p, p0i);
              x[0] = modp_montymul(gtB, mFp, p, p0i);
              x[llen] = modp_montymul(gtA, mFp, p, p0i);
              y[0] = modp_montymul(ftB, mGp, p, p0i);
              y[llen] = modp_montymul(ftA, mGp, p, p0i);
          }
          modp_iNTT2_ext(Ft + u, llen, igm, logn, p, p0i);
          modp_iNTT2_ext(Gt + u, llen, igm, logn, p, p0i);
  
          /*
           * Also save ft and gt (only up to size slen).
           */
          if (u < slen) {
              modp_iNTT2(fx, igm, logn, p, p0i);
              modp_iNTT2(gx, igm, logn, p, p0i);
              for (v = 0, x = ft + u, y = gt + u;
                      v < n; v ++, x += slen, y += slen) {
                  *x = fx[v];
                  *y = gx[v];
              }
          }
      }
  
      /*
       * Rebuild f, g, F and G with the CRT. Note that the elements of F
       * and G are consecutive, and thus can be rebuilt in a single
       * loop; similarly, the elements of f and g are consecutive.
       */
      zint_rebuild_CRT(Ft, llen, llen, n << 1, PRIMES, 1, t1);
      zint_rebuild_CRT(ft, slen, slen, n << 1, PRIMES, 1, t1);
  
      /*
       * Here starts the Babai reduction, specialized for depth = 1.
       *
       * Candidates F and G (from Ft and Gt), and base f and g (ft and gt),
       * are converted to floating point. There is no scaling, and a
       * single pass is sufficient.
       */
  
      /*
       * Convert F and G into floating point (rt1 and rt2).
       */
      rt1 = align_fpr(tmp, gt + slen * n);
      rt2 = rt1 + n;
      poly_big_to_fp(rt1, Ft, llen, llen, logn);
      poly_big_to_fp(rt2, Gt, llen, llen, logn);
  
      /*
       * Integer representation of F and G is no longer needed, we
       * can remove it.
       */
      memmove(tmp, ft, 2 * slen * n * sizeof * ft);
      ft = tmp;
      gt = ft + slen * n;
      rt3 = align_fpr(tmp, gt + slen * n);
      memmove(rt3, rt1, 2 * n * sizeof * rt1);
      rt1 = rt3;
      rt2 = rt1 + n;
      rt3 = rt2 + n;
      rt4 = rt3 + n;
  
      /*
       * Convert f and g into floating point (rt3 and rt4).
       */
      poly_big_to_fp(rt3, ft, slen, slen, logn);
      poly_big_to_fp(rt4, gt, slen, slen, logn);
  
      /*
       * Remove unneeded ft and gt.
       */
      memmove(tmp, rt1, 4 * n * sizeof * rt1);
      rt1 = (fpr *)tmp;
      rt2 = rt1 + n;
      rt3 = rt2 + n;
      rt4 = rt3 + n;
  
      /*
       * We now have:
       *   rt1 = F
       *   rt2 = G
       *   rt3 = f
       *   rt4 = g
       * in that order in RAM. We convert all of them to FFT.
       */
      FFT(rt1, logn);
      FFT(rt2, logn);
      FFT(rt3, logn);
      FFT(rt4, logn);
  
      /*
       * Compute:
       *   rt5 = F*adj(f) + G*adj(g)
       *   rt6 = 1 / (f*adj(f) + g*adj(g))
       * (Note that rt6 is half-length.)
       */
      rt5 = rt4 + n;
      rt6 = rt5 + n;
      poly_add_muladj_fft(rt5, rt1, rt2, rt3, rt4, logn);
      poly_invnorm2_fft(rt6, rt3, rt4, logn);
  
      /*
       * Compute:
       *   rt5 = (F*adj(f)+G*adj(g)) / (f*adj(f)+g*adj(g))
       */
      poly_mul_autoadj_fft(rt5, rt6, logn);
  
      /*
       * Compute k as the rounded version of rt5. Check that none of
       * the values is larger than 2^63-1 (in absolute value)
       * because that would make the fpr_rint() do something undefined;
       * note that any out-of-bounds value here implies a failure and
       * (f,g) will be discarded, so we can make a simple test.
       */
      iFFT(rt5, logn);
      for (u = 0; u < n; u ++) {
          fpr z;
  
          z = rt5[u];
          if (!fpr_lt(z, fpr_ptwo63m1) || !fpr_lt(fpr_mtwo63m1, z)) {
              return 0;
          }
          rt5[u] = fpr_of(fpr_rint(z));
      }
      FFT(rt5, logn);
  
      /*
       * Subtract k*f from F, and k*g from G.
       */
      poly_mul_fft(rt3, rt5, logn);
      poly_mul_fft(rt4, rt5, logn);
      poly_sub(rt1, rt3, logn);
      poly_sub(rt2, rt4, logn);
      iFFT(rt1, logn);
      iFFT(rt2, logn);
  
      /*
       * Convert back F and G to integers, and return.
       */
      Ft = tmp;
      Gt = Ft + n;
      rt3 = align_fpr(tmp, Gt + n);
      memmove(rt3, rt1, 2 * n * sizeof * rt1);
      rt1 = rt3;
      rt2 = rt1 + n;
      for (u = 0; u < n; u ++) {
          Ft[u] = (uint32_t)fpr_rint(rt1[u]);
          Gt[u] = (uint32_t)fpr_rint(rt2[u]);
      }
  
      return 1;
  }
  
  /*
   * Solving the NTRU equation, top level. Upon entry, the F and G
   * from the previous level should be in the tmp[] array.
   *
   * Returned value: 1 on success, 0 on error.
   */
  static int
  solve_NTRU_binary_depth0(unsigned logn,
                           const int8_t *f, const int8_t *g, uint32_t *tmp) {
      size_t n, hn, u;
      uint32_t p, p0i, R2;
      uint32_t *Fp, *Gp, *t1, *t2, *t3, *t4, *t5;
      uint32_t *gm, *igm, *ft, *gt;
      fpr *rt2, *rt3;
  
      n = (size_t)1 << logn;
      hn = n >> 1;
  
      /*
       * Equations are:
       *
       *   f' = f0^2 - X^2*f1^2
       *   g' = g0^2 - X^2*g1^2
       *   F' and G' are a solution to f'G' - g'F' = q (from deeper levels)
       *   F = F'*(g0 - X*g1)
       *   G = G'*(f0 - X*f1)
       *
       * f0, f1, g0, g1, f', g', F' and G' are all "compressed" to
       * degree N/2 (their odd-indexed coefficients are all zero).
       *
       * Everything should fit in 31-bit integers, hence we can just use
       * the first small prime p = 2147473409.
       */
      p = PRIMES[0].p;
      p0i = modp_ninv31(p);
      R2 = modp_R2(p, p0i);
  
      Fp = tmp;
      Gp = Fp + hn;
      ft = Gp + hn;
      gt = ft + n;
      gm = gt + n;
      igm = gm + n;
  
      modp_mkgm2(gm, igm, logn, PRIMES[0].g, p, p0i);
  
      /*
       * Convert F' anf G' in NTT representation.
       */
      for (u = 0; u < hn; u ++) {
          Fp[u] = modp_set(zint_one_to_plain(Fp + u), p);
          Gp[u] = modp_set(zint_one_to_plain(Gp + u), p);
      }
      modp_NTT2(Fp, gm, logn - 1, p, p0i);
      modp_NTT2(Gp, gm, logn - 1, p, p0i);
  
      /*
       * Load f and g and convert them to NTT representation.
       */
      for (u = 0; u < n; u ++) {
          ft[u] = modp_set(f[u], p);
          gt[u] = modp_set(g[u], p);
      }
      modp_NTT2(ft, gm, logn, p, p0i);
      modp_NTT2(gt, gm, logn, p, p0i);
  
      /*
       * Build the unreduced F,G in ft and gt.
       */
      for (u = 0; u < n; u += 2) {
          uint32_t ftA, ftB, gtA, gtB;
          uint32_t mFp, mGp;
  
          ftA = ft[u + 0];
          ftB = ft[u + 1];
          gtA = gt[u + 0];
          gtB = gt[u + 1];
          mFp = modp_montymul(Fp[u >> 1], R2, p, p0i);
          mGp = modp_montymul(Gp[u >> 1], R2, p, p0i);
          ft[u + 0] = modp_montymul(gtB, mFp, p, p0i);
          ft[u + 1] = modp_montymul(gtA, mFp, p, p0i);
          gt[u + 0] = modp_montymul(ftB, mGp, p, p0i);
          gt[u + 1] = modp_montymul(ftA, mGp, p, p0i);
      }
      modp_iNTT2(ft, igm, logn, p, p0i);
      modp_iNTT2(gt, igm, logn, p, p0i);
  
      Gp = Fp + n;
      t1 = Gp + n;
      memmove(Fp, ft, 2 * n * sizeof * ft);
  
      /*
       * We now need to apply the Babai reduction. At that point,
       * we have F and G in two n-word arrays.
       *
       * We can compute F*adj(f)+G*adj(g) and f*adj(f)+g*adj(g)
       * modulo p, using the NTT. We still move memory around in
       * order to save RAM.
       */
      t2 = t1 + n;
      t3 = t2 + n;
      t4 = t3 + n;
      t5 = t4 + n;
  
      /*
       * Compute the NTT tables in t1 and t2. We do not keep t2
       * (we'll recompute it later on).
       */
      modp_mkgm2(t1, t2, logn, PRIMES[0].g, p, p0i);
  
      /*
       * Convert F and G to NTT.
       */
      modp_NTT2(Fp, t1, logn, p, p0i);
      modp_NTT2(Gp, t1, logn, p, p0i);
  
      /*
       * Load f and adj(f) in t4 and t5, and convert them to NTT
       * representation.
       */
      t4[0] = t5[0] = modp_set(f[0], p);
      for (u = 1; u < n; u ++) {
          t4[u] = modp_set(f[u], p);
          t5[n - u] = modp_set(-f[u], p);
      }
      modp_NTT2(t4, t1, logn, p, p0i);
      modp_NTT2(t5, t1, logn, p, p0i);
  
      /*
       * Compute F*adj(f) in t2, and f*adj(f) in t3.
       */
      for (u = 0; u < n; u ++) {
          uint32_t w;
  
          w = modp_montymul(t5[u], R2, p, p0i);
          t2[u] = modp_montymul(w, Fp[u], p, p0i);
          t3[u] = modp_montymul(w, t4[u], p, p0i);
      }
  
      /*
       * Load g and adj(g) in t4 and t5, and convert them to NTT
       * representation.
       */
      t4[0] = t5[0] = modp_set(g[0], p);
      for (u = 1; u < n; u ++) {
          t4[u] = modp_set(g[u], p);
          t5[n - u] = modp_set(-g[u], p);
      }
      modp_NTT2(t4, t1, logn, p, p0i);
      modp_NTT2(t5, t1, logn, p, p0i);
  
      /*
       * Add G*adj(g) to t2, and g*adj(g) to t3.
       */
      for (u = 0; u < n; u ++) {
          uint32_t w;
  
          w = modp_montymul(t5[u], R2, p, p0i);
          t2[u] = modp_add(t2[u],
                           modp_montymul(w, Gp[u], p, p0i), p);
          t3[u] = modp_add(t3[u],
                           modp_montymul(w, t4[u], p, p0i), p);
      }
  
      /*
       * Convert back t2 and t3 to normal representation (normalized
       * around 0), and then
       * move them to t1 and t2. We first need to recompute the
       * inverse table for NTT.
       */
      modp_mkgm2(t1, t4, logn, PRIMES[0].g, p, p0i);
      modp_iNTT2(t2, t4, logn, p, p0i);
      modp_iNTT2(t3, t4, logn, p, p0i);
      for (u = 0; u < n; u ++) {
          t1[u] = (uint32_t)modp_norm(t2[u], p);
          t2[u] = (uint32_t)modp_norm(t3[u], p);
      }
  
      /*
       * At that point, array contents are:
       *
       *   F (NTT representation) (Fp)
       *   G (NTT representation) (Gp)
       *   F*adj(f)+G*adj(g) (t1)
       *   f*adj(f)+g*adj(g) (t2)
       *
       * We want to divide t1 by t2. The result is not integral; it
       * must be rounded. We thus need to use the FFT.
       */
  
      /*
       * Get f*adj(f)+g*adj(g) in FFT representation. Since this
       * polynomial is auto-adjoint, all its coordinates in FFT
       * representation are actually real, so we can truncate off
       * the imaginary parts.
       */
      rt3 = align_fpr(tmp, t3);
      for (u = 0; u < n; u ++) {
          rt3[u] = fpr_of(((int32_t *)t2)[u]);
      }
      FFT(rt3, logn);
      rt2 = align_fpr(tmp, t2);
      memmove(rt2, rt3, hn * sizeof * rt3);
  
      /*
       * Convert F*adj(f)+G*adj(g) in FFT representation.
       */
      rt3 = rt2 + hn;
      for (u = 0; u < n; u ++) {
          rt3[u] = fpr_of(((int32_t *)t1)[u]);
      }
      FFT(rt3, logn);
  
      /*
       * Compute (F*adj(f)+G*adj(g))/(f*adj(f)+g*adj(g)) and get
       * its rounded normal representation in t1.
       */
      poly_div_autoadj_fft(rt3, rt2, logn);
      iFFT(rt3, logn);
      for (u = 0; u < n; u ++) {
          t1[u] = modp_set((int32_t)fpr_rint(rt3[u]), p);
      }
  
      /*
       * RAM contents are now:
       *
       *   F (NTT representation) (Fp)
       *   G (NTT representation) (Gp)
       *   k (t1)
       *
       * We want to compute F-k*f, and G-k*g.
       */
      t2 = t1 + n;
      t3 = t2 + n;
      t4 = t3 + n;
      t5 = t4 + n;
      modp_mkgm2(t2, t3, logn, PRIMES[0].g, p, p0i);
      for (u = 0; u < n; u ++) {
          t4[u] = modp_set(f[u], p);
          t5[u] = modp_set(g[u], p);
      }
      modp_NTT2(t1, t2, logn, p, p0i);
      modp_NTT2(t4, t2, logn, p, p0i);
      modp_NTT2(t5, t2, logn, p, p0i);
      for (u = 0; u < n; u ++) {
          uint32_t kw;
  
          kw = modp_montymul(t1[u], R2, p, p0i);
          Fp[u] = modp_sub(Fp[u],
                           modp_montymul(kw, t4[u], p, p0i), p);
          Gp[u] = modp_sub(Gp[u],
                           modp_montymul(kw, t5[u], p, p0i), p);
      }
      modp_iNTT2(Fp, t3, logn, p, p0i);
      modp_iNTT2(Gp, t3, logn, p, p0i);
      for (u = 0; u < n; u ++) {
          Fp[u] = (uint32_t)modp_norm(Fp[u], p);
          Gp[u] = (uint32_t)modp_norm(Gp[u], p);
      }
  
      return 1;
  }
  
  /*
   * Solve the NTRU equation. Returned value is 1 on success, 0 on error.
   * G can be NULL, in which case that value is computed but not returned.
   * If any of the coefficients of F and G exceeds lim (in absolute value),
   * then 0 is returned.
   */
  static int
  solve_NTRU(unsigned logn, int8_t *F, int8_t *G,
             const int8_t *f, const int8_t *g, int lim, uint32_t *tmp) {
      size_t n, u;
      uint32_t *ft, *gt, *Ft, *Gt, *gm;
      uint32_t p, p0i, r;
      const small_prime *primes;
  
      n = MKN(logn);
  
      if (!solve_NTRU_deepest(logn, f, g, tmp)) {
          return 0;
      }
  
      /*
       * For logn <= 2, we need to use solve_NTRU_intermediate()
       * directly, because coefficients are a bit too large and
       * do not fit the hypotheses in solve_NTRU_binary_depth0().
       */
      if (logn <= 2) {
          unsigned depth;
  
          depth = logn;
          while (depth -- > 0) {
              if (!solve_NTRU_intermediate(logn, f, g, depth, tmp)) {
                  return 0;
              }
          }
      } else {
          unsigned depth;
  
          depth = logn;
          while (depth -- > 2) {
              if (!solve_NTRU_intermediate(logn, f, g, depth, tmp)) {
                  return 0;
              }
          }
          if (!solve_NTRU_binary_depth1(logn, f, g, tmp)) {
              return 0;
          }
          if (!solve_NTRU_binary_depth0(logn, f, g, tmp)) {
              return 0;
          }
      }
  
      /*
       * If no buffer has been provided for G, use a temporary one.
       */
      if (G == NULL) {
          G = (int8_t *)(tmp + 2 * n);
      }
  
      /*
       * Final F and G are in fk->tmp, one word per coefficient
       * (signed value over 31 bits).
       */
      if (!poly_big_to_small(F, tmp, lim, logn)
              || !poly_big_to_small(G, tmp + n, lim, logn)) {
          return 0;
      }
  
      /*
       * Verify that the NTRU equation is fulfilled. Since all elements
       * have short lengths, verifying modulo a small prime p works, and
       * allows using the NTT.
       *
       * We put Gt[] first in tmp[], and process it first, so that it does
       * not overlap with G[] in case we allocated it ourselves.
       */
      Gt = tmp;
      ft = Gt + n;
      gt = ft + n;
      Ft = gt + n;
      gm = Ft + n;
  
      primes = PRIMES;
      p = primes[0].p;
      p0i = modp_ninv31(p);
      modp_mkgm2(gm, tmp, logn, primes[0].g, p, p0i);
      for (u = 0; u < n; u ++) {
          Gt[u] = modp_set(G[u], p);
      }
      for (u = 0; u < n; u ++) {
          ft[u] = modp_set(f[u], p);
          gt[u] = modp_set(g[u], p);
          Ft[u] = modp_set(F[u], p);
      }
      modp_NTT2(ft, gm, logn, p, p0i);
      modp_NTT2(gt, gm, logn, p, p0i);
      modp_NTT2(Ft, gm, logn, p, p0i);
      modp_NTT2(Gt, gm, logn, p, p0i);
      r = modp_montymul(12289, 1, p, p0i);
      for (u = 0; u < n; u ++) {
          uint32_t z;
  
          z = modp_sub(modp_montymul(ft[u], Gt[u], p, p0i),
                       modp_montymul(gt[u], Ft[u], p, p0i), p);
          if (z != r) {
              return 0;
          }
      }
  
      return 1;
  }
  
  /*
   * Generate a random polynomial with a Gaussian distribution. This function
   * also makes sure that the resultant of the polynomial with phi is odd.
   */
  static void
  poly_small_mkgauss(RNG_CONTEXT *rng, int8_t *f, unsigned logn) {
      size_t n, u;
      unsigned mod2;
  
      n = MKN(logn);
      mod2 = 0;
      for (u = 0; u < n; u ++) {
          int s;
  
  restart:
          s = mkgauss(rng, logn);
  
          /*
           * We need the coefficient to fit within -127..+127;
           * realistically, this is always the case except for
           * the very low degrees (N = 2 or 4), for which there
           * is no real security anyway.
           */
          if (s < -127 || s > 127) {
              goto restart;
          }
  
          /*
           * We need the sum of all coefficients to be 1; otherwise,
           * the resultant of the polynomial with X^N+1 will be even,
           * and the binary GCD will fail.
           */
          if (u == n - 1) {
              if ((mod2 ^ (unsigned)(s & 1)) == 0) {
                  goto restart;
              }
          } else {
              mod2 ^= (unsigned)(s & 1);
          }
          f[u] = (int8_t)s;
      }
  }
  
  /* see falcon.h */
  void
  keygen(inner_shake256_context *rng,
                                 int8_t *f, int8_t *g, int8_t *F, int8_t *G, uint16_t *h,
                                 unsigned logn, uint8_t *tmp) {
      /*
       * Algorithm is the following:
       *
       *  - Generate f and g with the Gaussian distribution.
       *
       *  - If either Res(f,phi) or Res(g,phi) is even, try again.
       *
       *  - If ||(f,g)|| is too large, try again.
       *
       *  - If ||B~_{f,g}|| is too large, try again.
       *
       *  - If f is not invertible mod phi mod q, try again.
       *
       *  - Compute h = g/f mod phi mod q.
       *
       *  - Solve the NTRU equation fG - gF = q; if the solving fails,
       *    try again. Usual failure condition is when Res(f,phi)
       *    and Res(g,phi) are not prime to each other.
       */
      size_t n, u;
      uint16_t *h2, *tmp2;
      RNG_CONTEXT *rc;
  
      n = MKN(logn);
      rc = rng;
  
      /*
       * We need to generate f and g randomly, until we find values
       * such that the norm of (g,-f), and of the orthogonalized
       * vector, are satisfying. The orthogonalized vector is:
       *   (q*adj(f)/(f*adj(f)+g*adj(g)), q*adj(g)/(f*adj(f)+g*adj(g)))
       * (it is actually the (N+1)-th row of the Gram-Schmidt basis).
       *
       * In the binary case, coefficients of f and g are generated
       * independently of each other, with a discrete Gaussian
       * distribution of standard deviation 1.17*sqrt(q/(2*N)). Then,
       * the two vectors have expected norm 1.17*sqrt(q), which is
       * also our acceptance bound: we require both vectors to be no
       * larger than that (this will be satisfied about 1/4th of the
       * time, thus we expect sampling new (f,g) about 4 times for that
       * step).
       *
       * We require that Res(f,phi) and Res(g,phi) are both odd (the
       * NTRU equation solver requires it).
       */
      for (;;) {
          fpr *rt1, *rt2, *rt3;
          fpr bnorm;
          uint32_t normf, normg, norm;
          int lim;
  
          /*
           * The poly_small_mkgauss() function makes sure
           * that the sum of coefficients is 1 modulo 2
           * (i.e. the resultant of the polynomial with phi
           * will be odd).
           */
          poly_small_mkgauss(rc, f, logn);
          poly_small_mkgauss(rc, g, logn);
  
          /*
           * Verify that all coefficients are within the bounds
           * defined in max_fg_bits. This is the case with
           * overwhelming probability; this guarantees that the
           * key will be encodable with FALCON_COMP_TRIM.
           */
          lim = 1 << (max_fg_bits[logn] - 1);
          for (u = 0; u < n; u ++) {
              /*
               * We can use non-CT tests since on any failure
               * we will discard f and g.
               */
              if (f[u] >= lim || f[u] <= -lim
                      || g[u] >= lim || g[u] <= -lim) {
                  lim = -1;
                  break;
              }
          }
          if (lim < 0) {
              continue;
          }
  
          /*
           * Bound is 1.17*sqrt(q). We compute the squared
           * norms. With q = 12289, the squared bound is:
           *   (1.17^2)* 12289 = 16822.4121
           * Since f and g are integral, the squared norm
           * of (g,-f) is an integer.
           */
          normf = poly_small_sqnorm(f, logn);
          normg = poly_small_sqnorm(g, logn);
          norm = (normf + normg) | -((normf | normg) >> 31);
          if (norm >= 16823) {
              continue;
          }
  
          /*
           * We compute the orthogonalized vector norm.
           */
          rt1 = (fpr *)tmp;
          rt2 = rt1 + n;
          rt3 = rt2 + n;
          poly_small_to_fp(rt1, f, logn);
          poly_small_to_fp(rt2, g, logn);
          FFT(rt1, logn);
          FFT(rt2, logn);
          poly_invnorm2_fft(rt3, rt1, rt2, logn);
          poly_adj_fft(rt1, logn);
          poly_adj_fft(rt2, logn);
          poly_mulconst(rt1, fpr_q, logn);
          poly_mulconst(rt2, fpr_q, logn);
          poly_mul_autoadj_fft(rt1, rt3, logn);
          poly_mul_autoadj_fft(rt2, rt3, logn);
          iFFT(rt1, logn);
          iFFT(rt2, logn);
          bnorm = fpr_zero;
          for (u = 0; u < n; u ++) {
              bnorm = fpr_add(bnorm, fpr_sqr(rt1[u]));
              bnorm = fpr_add(bnorm, fpr_sqr(rt2[u]));
          }
          if (!fpr_lt(bnorm, fpr_bnorm_max)) {
              continue;
          }
  
          /*
           * Compute public key h = g/f mod X^N+1 mod q. If this
           * fails, we must restart.
           */
          if (h == NULL) {
              h2 = (uint16_t *)tmp;
              tmp2 = h2 + n;
          } else {
              h2 = h;
              tmp2 = (uint16_t *)tmp;
          }
          if (!compute_public(h2, f, g, logn, (uint8_t *)tmp2)) {
              continue;
          }
  
          /*
           * Solve the NTRU equation to get F and G.
           */
          lim = (1 << (max_FG_bits[logn] - 1)) - 1;
          if (!solve_NTRU(logn, F, G, f, g, lim, (uint32_t *)tmp)) {
              continue;
          }
  
          /*
           * Key pair is generated.
           */
          break;
      }
  }
  