Falcon source files (reference implementation)


falcon-vrfy.c

    1 /*
    2  * Falcon signature verification.
    3  *
    4  * ==========================(LICENSE BEGIN)============================
    5  *
    6  * Copyright (c) 2017  Falcon Project
    7  *
    8  * Permission is hereby granted, free of charge, to any person obtaining
    9  * a copy of this software and associated documentation files (the
   10  * "Software"), to deal in the Software without restriction, including
   11  * without limitation the rights to use, copy, modify, merge, publish,
   12  * distribute, sublicense, and/or sell copies of the Software, and to
   13  * permit persons to whom the Software is furnished to do so, subject to
   14  * the following conditions:
   15  *
   16  * The above copyright notice and this permission notice shall be
   17  * included in all copies or substantial portions of the Software.
   18  *
   19  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
   20  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
   21  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
   22  * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
   23  * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
   24  * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
   25  * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
   26  *
   27  * ===========================(LICENSE END)=============================
   28  *
   29  * @author   Thomas Pornin <thomas.pornin@nccgroup.trust>
   30  */
   31 
   32 #include "internal.h"
   33 
   34 /* ===================================================================== */
   35 /*
   36  * Constants for NTT.
   37  *
   38  * Binary case:
   39  *   n = 2^logn  (2 <= n <= 1024)
   40  *   phi = X^n + 1
   41  *   q = 12289
   42  *
   43  * Ternary case:
   44  *   n = 1.5*2^n  (6 <= n <= 768)
   45  *   phi = X^n - X^(n/2) + 1
   46  *   q = 18433
   47  *
   48  * In both cases:
   49  *   q0i = -1/q mod 2^16
   50  *   R = 2^16 mod q
   51  *   R2 = 2^32 mod q
   52  */
   53 
   54 #define Qb     12289
   55 #define Q0Ib   12287
   56 #define Rb      4091
   57 #define R2b    10952
   58 
   59 #define Qt     18433
   60 #define Q0It   18431
   61 #define Rt     10237
   62 #define R2t     4564
   63 
   64 /*
   65  * Inverse of n = 1.5*2^logn modulo q, for the ternary case. Values are
   66  * in Montogmery representation. This is defined only for logn from 2 to 9.
   67  */
   68 static const uint16_t INVNQt[] = {
   69         0, 0, 17067, 17750, 8875, 13654, 6827, 12630, 6315, 12374
   70 };
   71 
   72 /*
   73  * Table for NTT, binary case:
   74  *   GMb[x] = R*(g^rev(x)) mod q
   75  * where g = 7 (it is a 2048-th primitive root of 1 modulo q)
   76  * and rev() is the bit-reversal function over 10 bits.
   77  */
   78 static const uint16_t GMb[] = {
   79          4091,  7888, 11060, 11208,  6960,  4342,  6275,  9759,
   80          1591,  6399,  9477,  5266,   586,  5825,  7538,  9710,
   81          1134,  6407,  1711,   965,  7099,  7674,  3743,  6442,
   82         10414,  8100,  1885,  1688,  1364, 10329, 10164,  9180,
   83         12210,  6240,   997,   117,  4783,  4407,  1549,  7072,
   84          2829,  6458,  4431,  8877,  7144,  2564,  5664,  4042,
   85         12189,   432, 10751,  1237,  7610,  1534,  3983,  7863,
   86          2181,  6308,  8720,  6570,  4843,  1690,    14,  3872,
   87          5569,  9368, 12163,  2019,  7543,  2315,  4673,  7340,
   88          1553,  1156,  8401, 11389,  1020,  2967, 10772,  7045,
   89          3316, 11236,  5285, 11578, 10637, 10086,  9493,  6180,
   90          9277,  6130,  3323,   883, 10469,   489,  1502,  2851,
   91         11061,  9729,  2742, 12241,  4970, 10481, 10078,  1195,
   92           730,  1762,  3854,  2030,  5892, 10922,  9020,  5274,
   93          9179,  3604,  3782, 10206,  3180,  3467,  4668,  2446,
   94          7613,  9386,   834,  7703,  6836,  3403,  5351, 12276,
   95          3580,  1739, 10820,  9787, 10209,  4070, 12250,  8525,
   96         10401,  2749,  7338, 10574,  6040,   943,  9330,  1477,
   97          6865,  9668,  3585,  6633, 12145,  4063,  3684,  7680,
   98          8188,  6902,  3533,  9807,  6090,   727, 10099,  7003,
   99          6945,  1949,  9731, 10559,  6057,   378,  7871,  8763,
  100          8901,  9229,  8846,  4551,  9589, 11664,  7630,  8821,
  101          5680,  4956,  6251,  8388, 10156,  8723,  2341,  3159,
  102          1467,  5460,  8553,  7783,  2649,  2320,  9036,  6188,
  103           737,  3698,  4699,  5753,  9046,  3687,    16,   914,
  104          5186, 10531,  4552,  1964,  3509,  8436,  7516,  5381,
  105         10733,  3281,  7037,  1060,  2895,  7156,  8887,  5357,
  106          6409,  8197,  2962,  6375,  5064,  6634,  5625,   278,
  107           932, 10229,  8927,  7642,   351,  9298,   237,  5858,
  108          7692,  3146, 12126,  7586,  2053, 11285,  3802,  5204,
  109          4602,  1748, 11300,   340,  3711,  4614,   300, 10993,
  110          5070, 10049, 11616, 12247,  7421, 10707,  5746,  5654,
  111          3835,  5553,  1224,  8476,  9237,  3845,   250, 11209,
  112          4225,  6326,  9680, 12254,  4136,  2778,   692,  8808,
  113          6410,  6718, 10105, 10418,  3759,  7356, 11361,  8433,
  114          6437,  3652,  6342,  8978,  5391,  2272,  6476,  7416,
  115          8418, 10824, 11986,  5733,   876,  7030,  2167,  2436,
  116          3442,  9217,  8206,  4858,  5964,  2746,  7178,  1434,
  117          7389,  8879, 10661, 11457,  4220,  1432, 10832,  4328,
  118          8557,  1867,  9454,  2416,  3816,  9076,   686,  5393,
  119          2523,  4339,  6115,   619,   937,  2834,  7775,  3279,
  120          2363,  7488,  6112,  5056,   824, 10204, 11690,  1113,
  121          2727,  9848,   896,  2028,  5075,  2654, 10464,  7884,
  122         12169,  5434,  3070,  6400,  9132, 11672, 12153,  4520,
  123          1273,  9739, 11468,  9937, 10039,  9720,  2262,  9399,
  124         11192,   315,  4511,  1158,  6061,  6751, 11865,   357,
  125          7367,  4550,   983,  8534,  8352, 10126,  7530,  9253,
  126          4367,  5221,  3999,  8777,  3161,  6990,  4130, 11652,
  127          3374, 11477,  1753,   292,  8681,  2806, 10378, 12188,
  128          5800, 11811,  3181,  1988,  1024,  9340,  2477, 10928,
  129          4582,  6750,  3619,  5503,  5233,  2463,  8470,  7650,
  130          7964,  6395,  1071,  1272,  3474, 11045,  3291, 11344,
  131          8502,  9478,  9837,  1253,  1857,  6233,  4720, 11561,
  132          6034,  9817,  3339,  1797,  2879,  6242,  5200,  2114,
  133          7962,  9353, 11363,  5475,  6084,  9601,  4108,  7323,
  134         10438,  9471,  1271,   408,  6911,  3079,   360,  8276,
  135         11535,  9156,  9049, 11539,   850,  8617,   784,  7919,
  136          8334, 12170,  1846, 10213, 12184,  7827, 11903,  5600,
  137          9779,  1012,   721,  2784,  6676,  6552,  5348,  4424,
  138          6816,  8405,  9959,  5150,  2356,  5552,  5267,  1333,
  139          8801,  9661,  7308,  5788,  4910,   909, 11613,  4395,
  140          8238,  6686,  4302,  3044,  2285, 12249,  1963,  9216,
  141          4296, 11918,   695,  4371,  9793,  4884,  2411, 10230,
  142          2650,   841,  3890, 10231,  7248,  8505, 11196,  6688,
  143          4059,  6060,  3686,  4722, 11853,  5816,  7058,  6868,
  144         11137,  7926,  4894, 12284,  4102,  3908,  3610,  6525,
  145          7938,  7982, 11977,  6755,   537,  4562,  1623,  8227,
  146         11453,  7544,   906, 11816,  9548, 10858,  9703,  2815,
  147         11736,  6813,  6979,   819,  8903,  6271, 10843,   348,
  148          7514,  8339,  6439,   694,   852,  5659,  2781,  3716,
  149         11589,  3024,  1523,  8659,  4114, 10738,  3303,  5885,
  150          2978,  7289, 11884,  9123,  9323, 11830,    98,  2526,
  151          2116,  4131, 11407,  1844,  3645,  3916,  8133,  2224,
  152         10871,  8092,  9651,  5989,  7140,  8480,  1670,   159,
  153         10923,  4918,   128,  7312,   725,  9157,  5006,  6393,
  154          3494,  6043, 10972,  6181, 11838,  3423, 10514,  7668,
  155          3693,  6658,  6905, 11953, 10212, 11922,  9101,  8365,
  156          5110,    45,  2400,  1921,  4377,  2720,  1695,    51,
  157          2808,   650,  1896,  9997,  9971, 11980,  8098,  4833,
  158          4135,  4257,  5838,  4765, 10985, 11532,   590, 12198,
  159           482, 12173,  2006,  7064, 10018,  3912, 12016, 10519,
  160         11362,  6954,  2210,   284,  5413,  6601,  3865, 10339,
  161         11188,  6231,   517,  9564, 11281,  3863,  1210,  4604,
  162          8160, 11447,   153,  7204,  5763,  5089,  9248, 12154,
  163         11748,  1354,  6672,   179,  5532,  2646,  5941, 12185,
  164           862,  3158,   477,  7279,  5678,  7914,  4254,   302,
  165          2893, 10114,  6890,  9560,  9647, 11905,  4098,  9824,
  166         10269,  1353, 10715,  5325,  6254,  3951,  1807,  6449,
  167          5159,  1308,  8315,  3404,  1877,  1231,   112,  6398,
  168         11724, 12272,  7286,  1459, 12274,  9896,  3456,   800,
  169          1397, 10678,   103,  7420,  7976,   936,   764,   632,
  170          7996,  8223,  8445,  7758, 10870,  9571,  2508,  1946,
  171          6524, 10158,  1044,  4338,  2457,  3641,  1659,  4139,
  172          4688,  9733, 11148,  3946,  2082,  5261,  2036, 11850,
  173          7636, 12236,  5366,  2380,  1399,  7720,  2100,  3217,
  174         10912,  8898,  7578, 11995,  2791,  1215,  3355,  2711,
  175          2267,  2004,  8568, 10176,  3214,  2337,  1750,  4729,
  176          4997,  7415,  6315, 12044,  4374,  7157,  4844,   211,
  177          8003, 10159,  9290, 11481,  1735,  2336,  5793,  9875,
  178          8192,   986,  7527,  1401,   870,  3615,  8465,  2756,
  179          9770,  2034, 10168,  3264,  6132,    54,  2880,  4763,
  180         11805,  3074,  8286,  9428,  4881,  6933,  1090, 10038,
  181          2567,   708,   893,  6465,  4962, 10024,  2090,  5718,
  182         10743,   780,  4733,  4623,  2134,  2087,  4802,   884,
  183          5372,  5795,  5938,  4333,  6559,  7549,  5269, 10664,
  184          4252,  3260,  5917, 10814,  5768,  9983,  8096,  7791,
  185          6800,  7491,  6272,  1907, 10947,  6289, 11803,  6032,
  186         11449,  1171,  9201,  7933,  2479,  7970, 11337,  7062,
  187          8911,  6728,  6542,  8114,  8828,  6595,  3545,  4348,
  188          4610,  2205,  6999,  8106,  5560, 10390,  9321,  2499,
  189          2413,  7272,  6881, 10582,  9308,  9437,  3554,  3326,
  190          5991, 11969,  3415, 12283,  9838, 12063,  4332,  7830,
  191         11329,  6605, 12271,  2044, 11611,  7353, 11201, 11582,
  192          3733,  8943,  9978,  1627,  7168,  3935,  5050,  2762,
  193          7496, 10383,   755,  1654, 12053,  4952, 10134,  4394,
  194          6592,  7898,  7497,  8904, 12029,  3581, 10748,  5674,
  195         10358,  4901,  7414,  8771,   710,  6764,  8462,  7193,
  196          5371,  7274, 11084,   290,  7864,  6827, 11822,  2509,
  197          6578,  4026,  5807,  1458,  5721,  5762,  4178,  2105,
  198         11621,  4852,  8897,  2856, 11510,  9264,  2520,  8776,
  199          7011,  2647,  1898,  7039,  5950, 11163,  5488,  6277,
  200          9182, 11456,   633, 10046, 11554,  5633,  9587,  2333,
  201          7008,  7084,  5047,  7199,  9865,  8997,   569,  6390,
  202         10845,  9679,  8268, 11472,  4203,  1997,     2,  9331,
  203           162,  6182,  2000,  3649,  9792,  6363,  7557,  6187,
  204          8510,  9935,  5536,  9019,  3706, 12009,  1452,  3067,
  205          5494,  9692,  4865,  6019,  7106,  9610,  4588, 10165,
  206          6261,  5887,  2652, 10172,  1580, 10379,  4638,  9949
  207 };
  208 
  209 /*
  210  * Table for inverse NTT, binary case:
  211  *   iGMb[x] = R*((1/g)^rev(x)) mod q
  212  * Since g = 7, 1/g = 8778 mod 12289.
  213  */
  214 static const uint16_t iGMb[] = {
  215          4091,  4401,  1081,  1229,  2530,  6014,  7947,  5329,
  216          2579,  4751,  6464, 11703,  7023,  2812,  5890, 10698,
  217          3109,  2125,  1960, 10925, 10601, 10404,  4189,  1875,
  218          5847,  8546,  4615,  5190, 11324, 10578,  5882, 11155,
  219          8417, 12275, 10599,  7446,  5719,  3569,  5981, 10108,
  220          4426,  8306, 10755,  4679, 11052,  1538, 11857,   100,
  221          8247,  6625,  9725,  5145,  3412,  7858,  5831,  9460,
  222          5217, 10740,  7882,  7506, 12172, 11292,  6049,    79,
  223            13,  6938,  8886,  5453,  4586, 11455,  2903,  4676,
  224          9843,  7621,  8822,  9109,  2083,  8507,  8685,  3110,
  225          7015,  3269,  1367,  6397, 10259,  8435, 10527, 11559,
  226         11094,  2211,  1808,  7319,    48,  9547,  2560,  1228,
  227          9438, 10787, 11800,  1820, 11406,  8966,  6159,  3012,
  228          6109,  2796,  2203,  1652,   711,  7004,  1053,  8973,
  229          5244,  1517,  9322, 11269,   900,  3888, 11133, 10736,
  230          4949,  7616,  9974,  4746, 10270,   126,  2921,  6720,
  231          6635,  6543,  1582,  4868,    42,   673,  2240,  7219,
  232          1296, 11989,  7675,  8578, 11949,   989, 10541,  7687,
  233          7085,  8487,  1004, 10236,  4703,   163,  9143,  4597,
  234          6431, 12052,  2991, 11938,  4647,  3362,  2060, 11357,
  235         12011,  6664,  5655,  7225,  5914,  9327,  4092,  5880,
  236          6932,  3402,  5133,  9394, 11229,  5252,  9008,  1556,
  237          6908,  4773,  3853,  8780, 10325,  7737,  1758,  7103,
  238         11375, 12273,  8602,  3243,  6536,  7590,  8591, 11552,
  239          6101,  3253,  9969,  9640,  4506,  3736,  6829, 10822,
  240          9130,  9948,  3566,  2133,  3901,  6038,  7333,  6609,
  241          3468,  4659,   625,  2700,  7738,  3443,  3060,  3388,
  242          3526,  4418, 11911,  6232,  1730,  2558, 10340,  5344,
  243          5286,  2190, 11562,  6199,  2482,  8756,  5387,  4101,
  244          4609,  8605,  8226,   144,  5656,  8704,  2621,  5424,
  245         10812,  2959, 11346,  6249,  1715,  4951,  9540,  1888,
  246          3764,    39,  8219,  2080,  2502,  1469, 10550,  8709,
  247          5601,  1093,  3784,  5041,  2058,  8399, 11448,  9639,
  248          2059,  9878,  7405,  2496,  7918, 11594,   371,  7993,
  249          3073, 10326,    40, 10004,  9245,  7987,  5603,  4051,
  250          7894,   676, 11380,  7379,  6501,  4981,  2628,  3488,
  251         10956,  7022,  6737,  9933,  7139,  2330,  3884,  5473,
  252          7865,  6941,  5737,  5613,  9505, 11568, 11277,  2510,
  253          6689,   386,  4462,   105,  2076, 10443,   119,  3955,
  254          4370, 11505,  3672, 11439,   750,  3240,  3133,   754,
  255          4013, 11929,  9210,  5378, 11881, 11018,  2818,  1851,
  256          4966,  8181,  2688,  6205,  6814,   926,  2936,  4327,
  257         10175,  7089,  6047,  9410, 10492,  8950,  2472,  6255,
  258           728,  7569,  6056, 10432, 11036,  2452,  2811,  3787,
  259           945,  8998,  1244,  8815, 11017, 11218,  5894,  4325,
  260          4639,  3819,  9826,  7056,  6786,  8670,  5539,  7707,
  261          1361,  9812,  2949, 11265, 10301,  9108,   478,  6489,
  262           101,  1911,  9483,  3608, 11997, 10536,   812,  8915,
  263           637,  8159,  5299,  9128,  3512,  8290,  7068,  7922,
  264          3036,  4759,  2163,  3937,  3755, 11306,  7739,  4922,
  265         11932,   424,  5538,  6228, 11131,  7778, 11974,  1097,
  266          2890, 10027,  2569,  2250,  2352,   821,  2550, 11016,
  267          7769,   136,   617,  3157,  5889,  9219,  6855,   120,
  268          4405,  1825,  9635,  7214, 10261, 11393,  2441,  9562,
  269         11176,   599,  2085, 11465,  7233,  6177,  4801,  9926,
  270          9010,  4514,  9455, 11352, 11670,  6174,  7950,  9766,
  271          6896, 11603,  3213,  8473,  9873,  2835, 10422,  3732,
  272          7961,  1457, 10857,  8069,   832,  1628,  3410,  4900,
  273         10855,  5111,  9543,  6325,  7431,  4083,  3072,  8847,
  274          9853, 10122,  5259, 11413,  6556,   303,  1465,  3871,
  275          4873,  5813, 10017,  6898,  3311,  5947,  8637,  5852,
  276          3856,   928,  4933,  8530,  1871,  2184,  5571,  5879,
  277          3481, 11597,  9511,  8153,    35,  2609,  5963,  8064,
  278          1080, 12039,  8444,  3052,  3813, 11065,  6736,  8454,
  279          2340,  7651,  1910, 10709,  2117,  9637,  6402,  6028,
  280          2124,  7701,  2679,  5183,  6270,  7424,  2597,  6795,
  281          9222, 10837,   280,  8583,  3270,  6753,  2354,  3779,
  282          6102,  4732,  5926,  2497,  8640, 10289,  6107, 12127,
  283          2958, 12287, 10292,  8086,   817,  4021,  2610,  1444,
  284          5899, 11720,  3292,  2424,  5090,  7242,  5205,  5281,
  285          9956,  2702,  6656,   735,  2243, 11656,   833,  3107,
  286          6012,  6801,  1126,  6339,  5250, 10391,  9642,  5278,
  287          3513,  9769,  3025,   779,  9433,  3392,  7437,   668,
  288         10184,  8111,  6527,  6568, 10831,  6482,  8263,  5711,
  289          9780,   467,  5462,  4425, 11999,  1205,  5015,  6918,
  290          5096,  3827,  5525, 11579,  3518,  4875,  7388,  1931,
  291          6615,  1541,  8708,   260,  3385,  4792,  4391,  5697,
  292          7895,  2155,  7337,   236, 10635, 11534,  1906,  4793,
  293          9527,  7239,  8354,  5121, 10662,  2311,  3346,  8556,
  294           707,  1088,  4936,   678, 10245,    18,  5684,   960,
  295          4459,  7957,   226,  2451,     6,  8874,   320,  6298,
  296          8963,  8735,  2852,  2981,  1707,  5408,  5017,  9876,
  297          9790,  2968,  1899,  6729,  4183,  5290, 10084,  7679,
  298          7941,  8744,  5694,  3461,  4175,  5747,  5561,  3378,
  299          5227,   952,  4319,  9810,  4356,  3088, 11118,   840,
  300          6257,   486,  6000,  1342, 10382,  6017,  4798,  5489,
  301          4498,  4193,  2306,  6521,  1475,  6372,  9029,  8037,
  302          1625,  7020,  4740,  5730,  7956,  6351,  6494,  6917,
  303         11405,  7487, 10202, 10155,  7666,  7556, 11509,  1546,
  304          6571, 10199,  2265,  7327,  5824, 11396, 11581,  9722,
  305          2251, 11199,  5356,  7408,  2861,  4003,  9215,   484,
  306          7526,  9409, 12235,  6157,  9025,  2121, 10255,  2519,
  307          9533,  3824,  8674, 11419, 10888,  4762, 11303,  4097,
  308          2414,  6496,  9953, 10554,   808,  2999,  2130,  4286,
  309         12078,  7445,  5132,  7915,   245,  5974,  4874,  7292,
  310          7560, 10539,  9952,  9075,  2113,  3721, 10285, 10022,
  311          9578,  8934, 11074,  9498,   294,  4711,  3391,  1377,
  312          9072, 10189,  4569, 10890,  9909,  6923,    53,  4653,
  313           439, 10253,  7028, 10207,  8343,  1141,  2556,  7601,
  314          8150, 10630,  8648,  9832,  7951, 11245,  2131,  5765,
  315         10343,  9781,  2718,  1419,  4531,  3844,  4066,  4293,
  316         11657, 11525, 11353,  4313,  4869, 12186,  1611, 10892,
  317         11489,  8833,  2393,    15, 10830,  5003,    17,   565,
  318          5891, 12177, 11058, 10412,  8885,  3974, 10981,  7130,
  319          5840, 10482,  8338,  6035,  6964,  1574, 10936,  2020,
  320          2465,  8191,   384,  2642,  2729,  5399,  2175,  9396,
  321         11987,  8035,  4375,  6611,  5010, 11812,  9131, 11427,
  322           104,  6348,  9643,  6757, 12110,  5617, 10935,   541,
  323           135,  3041,  7200,  6526,  5085, 12136,   842,  4129,
  324          7685, 11079,  8426,  1008,  2725, 11772,  6058,  1101,
  325          1950,  8424,  5688,  6876, 12005, 10079,  5335,   927,
  326          1770,   273,  8377,  2271,  5225, 10283,   116, 11807,
  327            91, 11699,   757,  1304,  7524,  6451,  8032,  8154,
  328          7456,  4191,   309,  2318,  2292, 10393, 11639,  9481,
  329         12238, 10594,  9569,  7912, 10368,  9889, 12244,  7179,
  330          3924,  3188,   367,  2077,   336,  5384,  5631,  8596,
  331          4621,  1775,  8866,   451,  6108,  1317,  6246,  8795,
  332          5896,  7283,  3132, 11564,  4977, 12161,  7371,  1366,
  333         12130, 10619,  3809,  5149,  6300,  2638,  4197,  1418,
  334         10065,  4156,  8373,  8644, 10445,   882,  8158, 10173,
  335          9763, 12191,   459,  2966,  3166,   405,  5000,  9311,
  336          6404,  8986,  1551,  8175,  3630, 10766,  9265,   700,
  337          8573,  9508,  6630, 11437, 11595,  5850,  3950,  4775,
  338         11941,  1446,  6018,  3386, 11470,  5310,  5476,   553,
  339          9474,  2586,  1431,  2741,   473, 11383,  4745,   836,
  340          4062, 10666,  7727, 11752,  5534,   312,  4307,  4351,
  341          5764,  8679,  8381,  8187,     5,  7395,  4363,  1152,
  342          5421,  5231,  6473,   436,  7567,  8603,  6229,  8230
  343 };
  344 
  345 /*
  346  * Tables for NTT and inverse NTT, ternary case.
  347  *
  348  *   GMt_square: for degree halving and doubling (inverse: iGMt_square)
  349  *   GMt_cubic: for ternary split/merge (inverse: iGMt_cubic)
  350  */
  351 
  352 static const uint16_t GMt_square[] = {
  353          9358,  9358,  8086, 11703,  5774, 14509, 12722,  9851,
  354         16750, 12828,  7852,   806,  1458, 10770,  9265, 12609,
  355          8775,  1328, 14458, 11372,  1705,  1823,  7105,  6894,
  356          9026,    72, 11666,  7057,  4201,  8427, 18263, 14143,
  357         11872,  6834,  4390,  7775, 13188, 11852, 17658,  7550,
  358          8671,  4125, 12457, 11838, 14110,  5843,  1013, 16889,
  359         15905,  5600, 12331, 18417,  5191,  4134,  9307, 10416,
  360           141, 17654, 14588, 12484,   374,  9438, 14640,  1869,
  361         17998, 16130, 13431, 13647,  6690,  6180, 12094,   509,
  362         12223, 13523,  4231,  1594,  5883,  7501, 10569, 12987,
  363         17487, 15162,  2004,   694,  7557,  9626,  4139,  9031,
  364         10013, 13052,  3184,  2280,  6819,   761, 10145,  8793,
  365         15397,  5792,   430,  6514,  4105,  8173, 14998, 17409,
  366          9415, 15310,  5503, 14176,  9024,  5443, 11982,  6357,
  367          4076,  3104,  1147,  7259,   876,  6926,  9056, 11672,
  368          8610, 11260,  3662,  8921, 16955,  6074, 12328, 17257,
  369         18117,   700, 13062, 18431,  2953,  5125, 12684,  1302,
  370          6930,  6815, 11040, 10777,  4655,  5788,  1830,  7146,
  371          5330,  8726,  5778,  3767, 14007, 15171, 17287, 17705,
  372         16342,  2532, 17017,  5470, 15632, 10638,   166, 15032,
  373          3832, 13211,  2833, 14024, 12256,  7850, 17450, 13144,
  374         14661,  9989,  6120,  6976,  3983,  4010, 15841, 11575,
  375          8164, 10848,   398,   285, 12373, 16224, 17397, 17228,
  376          7857, 15028, 12038,  3433,  9467,  4695, 15720, 13943,
  377          4443,  3691, 16893,  6678,  1588, 11882,  7158,  2810,
  378         12578,  9470,  3440, 15246, 14407, 10085,  9386, 10241,
  379          9408,  6459,  6609, 11726, 13581, 16348, 10863, 16069,
  380         14175,  6399,  9176,  2773,  7008,   109, 17149,  1211,
  381         15887, 17073, 15175, 12117, 16909,   576,  1163,  1157,
  382          4969, 10459,  7517,  6448,  9389, 11401,  9611,  5076,
  383          2811, 17806, 16687,  6901, 13339,  2651, 12233,  5101,
  384         14069, 14567,  7491,  2539,  2282,  9878,  8104,  6081
  385 };
  386 
  387 static const uint16_t GMt_cubic[] = {
  388             0,     0,  5863, 13355,  3451, 10413, 17691,  5912,
  389          1638, 13729,   651,  6638,  5184, 14889,  7732, 13716,
  390          6193, 12464,  2225,  4481,  6221, 17110, 16888,  3019,
  391         18432,  3784, 17251, 11902,  2776,  2426,   158, 10417,
  392          2894, 16739, 10603,  6889,  3044,  2129,  3573,  9590,
  393          9271, 14968,  9120, 14929, 14605, 15247,  9822, 12913,
  394          5467, 13131, 10444,   256, 12400,  8818,  2565,  8231,
  395          1966,  7588,  1254, 10578, 16985,  4631,  2733, 17674,
  396          8301, 17281,  5426,  2378, 16107,  9043, 15618, 16119,
  397         12790,  7698,  2720, 11567, 15351, 12632,  6810,   294,
  398          3399,  4418, 17657,  5537,  2072, 12010, 15948,  2410,
  399         16169, 14064, 15170, 15515, 17644, 17863,  7485,  8281,
  400         14259, 15768,  6376,  2013, 11100,  6407, 14337, 15544,
  401         11571, 12144, 18069, 13334,  7623,  2213, 15082, 16713,
  402          5319,  1740,  1405, 10617, 17722, 17639,  7516,  1575,
  403          3339, 10262,  2036,   770,  2735, 10106,  6995,   708,
  404         10542, 16517, 18369,  2547,  7012, 10112, 11767,  7800,
  405         16536,  7811,  6572, 16102, 12667, 12305,  4798,   873,
  406          2005,  7476, 10486,  7225,   886,  2182, 15004, 16937,
  407          4939,  1886, 13070, 17292,  3488, 17869, 12257, 15373,
  408          7292,  1273, 10933, 11613, 15275,  5288,  9143,  1629,
  409         11564,  1766,  9795,  4483,  8622,   762, 16188, 15900,
  410         14917, 14351,  9946,  4522,  9359, 13770,  2538, 18234,
  411          6214,  6732,  8614, 12601,  3224,  3030, 13570,  5458,
  412         11553,  6524, 15226,  6374,  2292,  9015, 17926,  1456,
  413          6036, 16696,   981, 11362, 18094, 10899,  4828, 16384,
  414          2832, 11718, 11051,  7493,  9259,  5077, 13369, 10289,
  415          4117, 15590, 18415, 12813, 18101,  2844, 13102,  6802,
  416          4802,  4170, 17033,  7329, 15140,     4, 15470,  4728,
  417         11498, 11881,  5515, 15829,  7508, 13414,  8183,  2968,
  418            81,  6857,  3577, 12887, 14773,  6257,  5635,  4141,
  419          4355, 18215,  4803,   386,  2568, 15312, 12364, 16011,
  420          9971,  2087,  7035, 15245,  6870, 12883,  9820,  2048,
  421          9503,  3431,  6849,   182, 15728,  5405, 10032, 10892,
  422          7427,  6557,  4606,  8514,  9175,  9572,  6246, 14675,
  423         12678,  7547, 17800, 17415, 12902,  7849,  6073,  5719,
  424         15262, 17614, 12210,  8891, 10155,  6285,  3327,   371,
  425         11109,  9217,  6542,   591, 18258, 17045, 14346, 18354,
  426         12065,  4581, 12121, 13873,   321,  1914, 10762, 13522,
  427          8759, 16911, 12225,  7430, 16576,  3915, 16986,   847,
  428         11952,  8214,  7586, 13190,   648, 17990, 10183, 10931,
  429          7690,  6747,  2111, 11898, 16407, 16689,  1558,  3088,
  430          4854, 10165,  4765, 15147, 18252,  2883,  7254, 16034,
  431          1550, 14927,  7233,  3333, 10522,    32, 13162,   958,
  432         18150,  1758, 15721, 13460, 11422,  4537,  7848, 17164,
  433           259, 15326, 11210, 14126, 18336, 16821, 14377, 11648,
  434         13534, 12651, 15777,  4319, 14503, 14122, 18289, 10339,
  435          4223,  1579, 14676,  4645,   340,  3750, 14787,  8580,
  436         10662,  4829, 12745, 12081,  5686, 13920, 11240, 11204,
  437         13430,   661,  3447,  7116,  8279,  8364, 16288,  6160,
  438          5385, 10058,  5685, 17704,   403,  4987, 15521, 14507,
  439          3474, 15546, 14142, 16104, 15068, 14390,  4098, 13754,
  440         12138,  4844,  6242, 11378,   436,  9146, 17661,  8834,
  441          4719,  4881, 11092, 18246,  5919, 17032, 10151,  2988,
  442         10093,  1264,  3775,   975, 18425, 11839,  8977,  3051,
  443         13104, 17667,  5208, 16238, 10038,  6621, 12497, 10430,
  444         12967,  1518,  9171,  6275,  3257,  7189, 15710, 18218,
  445         10604,  3105, 17921,  1943,   797,  7164,  1971,  7101,
  446         11938,  5891,  9471, 13921,  2646, 15088, 12395,  9305,
  447         16040,  4509, 10156,  2501,  7088, 17456,  9434,  6465,
  448          2304,   473, 13677,  6096,   347, 14128,  4628, 17431,
  449          3037, 10184, 13732,   739, 11602,  5438, 17845, 13032,
  450          9597, 16395,  7359,  5807, 12846, 16990, 13613,  8643,
  451          8738,  4210,  5836, 17743,  1140, 17995,  1871, 16841
  452 };
  453 
  454 static const uint16_t iGMt_square[] = {
  455          3318,   879,  6730, 10347,  8582,  5711,  3924, 12659,
  456          5824,  9168,  7663, 16975, 17627, 10581,  5605,  1683,
  457          4290,   170, 10006, 14232, 11376,  6767, 18361,  9407,
  458         11539, 11328, 16610, 16728,  7061,  3975, 17105,  9658,
  459         16564,  3793,  8995, 18059,  5949,  3845,   779, 18292,
  460          8017,  9126, 14299, 13242,    16,  6102, 12833,  2528,
  461          1544, 17420, 12590,  4323,  6595,  5976, 14308,  9762,
  462         10883,   775,  6581,  5245, 10658, 14043, 11599,  6561,
  463          1176,  6105, 12359,  1478,  9512, 14771,  7173,  9823,
  464          6761,  9377, 11507, 17557, 11174, 17286, 15329, 14357,
  465         12076,  6451, 12990,  9409,  4257, 12930,  3123,  9018,
  466          1024,  3435, 10260, 14328, 11919, 18003, 12641,  3036,
  467          9640,  8288, 17672, 11614, 16153, 15249,  5381,  8420,
  468          9402, 14294,  8807, 10876, 17739, 16429,  3271,   946,
  469          5446,  7864, 10932, 12550, 16839, 14202,  4910,  6210,
  470         17924,  6339, 12253, 11743,  4786,  5002,  2303,   435,
  471         12352, 10329,  8555, 16151, 15894, 10942,  3866,  4364,
  472         13332,  6200, 15782,  5094, 11532,  1746,   627, 15622,
  473         13357,  8822,  7032,  9044, 11985, 10916,  7974, 13464,
  474         17276, 17270, 17857,  1524,  6316,  3258,  1360,  2546,
  475         17222,  1284, 18324, 11425, 15660,  9257, 12034,  4258,
  476          2364,  7570,  2085,  4852,  6707, 11824, 11974,  9025,
  477          8192,  9047,  8348,  4026,  3187, 14993,  8963,  5855,
  478         15623, 11275,  6551, 16845, 11755,  1540, 14742, 13990,
  479          4490,  2713, 13738,  8966, 15000,  6395,  3405, 10576,
  480          1205,  1036,  2209,  6060, 18148, 18035,  7585, 10269,
  481          6858,  2592, 14423, 14450, 11457, 12313,  8444,  3772,
  482          5289,   983, 10583,  6177,  4409, 15600,  5222, 14601,
  483          3401, 18267,  7795,  2801, 12963,  1416, 15901,  2091,
  484           728,  1146,  3262,  4426, 14666, 12655,  9707, 13103,
  485         11287, 16603, 12645, 13778,  7656,  7393, 11618, 11503,
  486         17131,  5749, 13308, 15480,     2,  5371, 17733,   316
  487 };
  488 
  489 static const uint16_t iGMt_cubic[] = {
  490             0,     0, 10467, 10693, 11779, 12521, 11471,  8020,
  491         12449,  4717,  8728,  3544, 12446, 11795,  6342,  4704,
  492          8174,  8016,   350, 16007,  5349,  6531, 14648, 14649,
  493         13869, 15414,  7544,  1323, 16177, 13952, 12162,  5969,
  494          3492,   759, 12354, 13802,  9109,  7855, 12811, 10845,
  495         12767, 10202,  3582,  9615, 10188, 18177, 10769,  5302,
  496         15342,  5520, 17791,  3186, 12624,  3504, 12736,  3465,
  497         12416,  8843,   915, 16304,  3714, 11544,  4588,  1694,
  498          6287, 17725, 11062,  8327,  1266, 17663, 11510,  8171,
  499          5941, 16858,    83,   794,  9221,  7816,  3579, 16693,
  500         16802,  1720,  5410, 16220,  4735,  5099, 17860,  6289,
  501         17226,  2889,  4693, 12026,  4363, 16420, 16924,  2665,
  502         17637, 10152, 18214,   570, 18088,  2918,  2105,  4369,
  503         13538, 16023,  8495,  6423, 12120, 12896, 17414, 14015,
  504          6516, 18139,  2719,  5801,  9586,  6866,  5092, 10735,
  505         17932,  2314,  7064,  9390,  3048, 16055,  9453,  1152,
  506         14786,  2422,  5689,  3121,  4417, 18047,  4573,   218,
  507          1494, 14292,  8516, 12176,  9123,  5546, 11657, 11576,
  508          5215, 15465, 12527,  5019,  8119,  2604, 18050,  6552,
  509         10742, 13705, 15136, 18429,  9704, 11104,   632, 14263,
  510          6300, 11631, 15257, 15589,  5602,  5620,  6960,  2843,
  511          3080,  8144,  4182, 13356,  3558, 10940,  9547,  6715,
  512          6877,  2049,  7195,  7534,  8052,  7071,  7773,  1737,
  513         16470, 16977, 11710,  9418,  8852, 12059,  5029, 11909,
  514          8112, 12975,   194, 15403, 14446,  5832, 17915, 11701,
  515          2737,   199, 14022,  4663,  5424, 13911,   566,  4082,
  516           288,  2533,  7860, 17671,  5312, 13950,  9798, 16667,
  517          7514, 16804,  9987, 13145, 17753,  6820,  6019, 17160,
  518         15317,  3060,  4052,   564, 14211,  1141,  3053, 16547,
  519         16500,  1496, 17137, 16251,  3261, 11208, 12962, 10957,
  520          3925, 17560,   362,  6128,  8903,  2331,  8725, 10622,
  521          3967, 10633, 15333,  8321, 15822, 15886, 12458,  1916,
  522          3463,  1592,  1578,   438,  6526,   690,  4528, 14223,
  523          4970,  9790, 14289,  1443,  1552, 12626, 11635,  2038,
  524          4813,  5401,  6164, 12995, 12993, 17694, 11286,  8249,
  525          5630,  1002,  4652,  4305,  7581, 12337,  1831, 17960,
  526          2969, 11968,  8065,   977,  7655, 15932, 11531, 13924,
  527          3090,  9128,  5991,  3345, 13983,  4512,  6047, 12542,
  528         13303, 11332, 12066, 11269, 15978, 16490,  7499, 15328,
  529         15925,   215, 14501, 11244,  2896, 12158, 11449, 16915,
  530          2067,  8003,  3417, 11812,  7403,  2195, 13870,   766,
  531          5926, 15382,  6586,  6594,  2800, 17458,  8829, 17169,
  532          7163, 15445,  7320,  1401, 11279,   187, 18271, 13552,
  533          8827,  9599,  9723,  9287, 13297,  7055,  7294, 13589,
  534          8777,  4679,   678,  4043, 16471,  2329,  6361,  2887,
  535          1014,  3926, 13849, 13446,  6414,   729, 13760,  8375,
  536         10128, 12273, 18348, 10069, 14764, 11317, 12769, 17772,
  537            36,  7229, 10199,  4513,   664,  6352,  5833, 13604,
  538          6207,  9853, 15023, 14683, 10031, 13788,  2644, 16854,
  539          7950,  8094,   381,  4311, 11458, 14114,   883,  5782,
  540          2729,  6785,  1515,  1612, 15517,  4307,  3366,  3107,
  541          9117,  1269,  6885, 13896,  2261,  4973, 16392, 16675,
  542         12204, 17475, 10490, 18401,  3900, 15100,  5056,  3506,
  543          9653,  2399, 15369, 15550,  8051,  3286, 13122,  8268,
  544         16903, 15345, 18151,  1744,  8646,  6535,   943, 11686,
  545         17685,  7502,  1091,   443, 12829,  5243,  3738, 10219,
  546         16139, 17586, 12661, 14518,  4795, 11003, 10281,  1522,
  547         15673,  4911, 16840, 16519, 16681,  4560,  7484, 13852,
  548         14425,    79,  1213,  1388,  5951, 17842,  1892,  9216,
  549          2956, 18062,  3870, 12148,  3319,  9542, 16081,   819,
  550           354, 12714,  5053, 10584,   385,  1018,  5131, 10886,
  551         10004,  3758, 18036,  8861, 14525,  9919,   870, 11876,
  552         17573,  7541, 10323, 13028,  6667, 18251,  6072, 15002,
  553          7772, 16385, 12420,  5550, 10223,  3188,  7884, 16346
  554 };
  555 
  556 /*
  557  * Reduce a small signed integer modulo q. The source integer MUST
  558  * be between -q/2 and +q/2.
  559  */
  560 static inline uint32_t
  561 mq_conv_small(int x, uint32_t q)
  562 {
  563         /*
  564          * If x < 0, the cast to uint32_t will set the high bit to 1.
  565          */
  566         uint32_t y;
  567 
  568         y = (uint32_t)x;
  569         y += q & -(y >> 31);
  570         return y;
  571 }
  572 
  573 /*
  574  * Addition modulo q. Operands must be in the 0..q-1 range.
  575  */
  576 static inline uint32_t
  577 mq_add(uint32_t x, uint32_t y, uint32_t q)
  578 {
  579         /*
  580          * We compute x + y - q. If the result is negative, then the
  581          * high bit will be set, and 'd >> 31' will be equal to 1;
  582          * thus '-(d >> 31)' will be an all-one pattern. Otherwise,
  583          * it will be an all-zero pattern. In other words, this
  584          * implements a conditional addition of q.
  585          */
  586         uint32_t d;
  587 
  588         d = x + y - q;
  589         d += q & -(d >> 31);
  590         return d;
  591 }
  592 
  593 /*
  594  * Subtraction modulo q. Operands must be in the 0..q-1 range.
  595  */
  596 static inline uint32_t
  597 mq_sub(uint32_t x, uint32_t y, uint32_t q)
  598 {
  599         /*
  600          * As in mq_add(), we use a conditional addition to ensure the
  601          * result is in the 0..q-1 range.
  602          */
  603         uint32_t d;
  604 
  605         d = x - y;
  606         d += q & -(d >> 31);
  607         return d;
  608 }
  609 
  610 /*
  611  * Division by 2 modulo q. Operand must be in the 0..q-1 range.
  612  */
  613 static inline uint32_t
  614 mq_rshift1(uint32_t x, uint32_t q)
  615 {
  616         x += q & -(x & 1);
  617         return (x >> 1);
  618 }
  619 
  620 /*
  621  * Montgomery multiplication modulo q. If we set R = 2^16 mod q, then
  622  * this function computes: x * y / R mod q
  623  * Operands must be in the 0..q-1 range.
  624  */
  625 static inline uint32_t
  626 mq_montymul(uint32_t x, uint32_t y, uint32_t q, uint32_t q0i)
  627 {
  628         uint32_t z, w;
  629 
  630         /*
  631          * We compute x*y + k*q with a value of k chosen so that the 16
  632          * low bits of the result are 0. We can then shift the value.
  633          * After the shift, result may still be larger than q, but it
  634          * will be lower than 2*q, so a conditional subtraction works.
  635          */
  636 
  637         z = x * y;
  638         w = ((z * q0i) & 0xFFFF) * q;
  639 
  640         /*
  641          * When adding z and w, the result will have its low 16 bits
  642          * equal to 0. Since x, y and z are lower than q, the sum will
  643          * be no more than (2^15 - 1) * q + (q - 1)^2, which will
  644          * fit on 29 bits.
  645          */
  646         z = (z + w) >> 16;
  647 
  648         /*
  649          * After the shift, analysis shows that the value will be less
  650          * than 2q. We do a subtraction then conditional subtraction to
  651          * ensure the result is in the expected range.
  652          */
  653         z -= q;
  654         z += q & -(z >> 31);
  655         return z;
  656 }
  657 
  658 /*
  659  * Montgomery squaring (computes (x^2)/R).
  660  */
  661 static inline uint32_t
  662 mq_montysqr(uint32_t x, uint32_t q, uint32_t q0i)
  663 {
  664         return mq_montymul(x, x, q, q0i);
  665 }
  666 
  667 /*
  668  * Divide x by y modulo q = 12289.
  669  */
  670 static inline uint32_t
  671 mq_div_12289(uint32_t x, uint32_t y)
  672 {
  673         /*
  674          * We invert y by computing y^(q-2) mod q.
  675          *
  676          * We use the following addition chain for exponent e = 12287:
  677          *
  678          *   e0 = 1
  679          *   e1 = 2 * e0 = 2
  680          *   e2 = e1 + e0 = 3
  681          *   e3 = e2 + e1 = 5
  682          *   e4 = 2 * e3 = 10
  683          *   e5 = 2 * e4 = 20
  684          *   e6 = 2 * e5 = 40
  685          *   e7 = 2 * e6 = 80
  686          *   e8 = 2 * e7 = 160
  687          *   e9 = e8 + e2 = 163
  688          *   e10 = e9 + e8 = 323
  689          *   e11 = 2 * e10 = 646
  690          *   e12 = 2 * e11 = 1292
  691          *   e13 = e12 + e9 = 1455
  692          *   e14 = 2 * e13 = 2910
  693          *   e15 = 2 * e14 = 5820
  694          *   e16 = e15 + e10 = 6143
  695          *   e17 = 2 * e16 = 12286
  696          *   e18 = e17 + e0 = 12287
  697          *
  698          * Additions on exponents are converted to Montgomery
  699          * multiplications. We define all intermediate results as so
  700          * many local variables, and let the C compiler work out which
  701          * must be kept around.
  702          */
  703         uint32_t y0, y1, y2, y3, y4, y5, y6, y7, y8, y9;
  704         uint32_t y10, y11, y12, y13, y14, y15, y16, y17, y18;
  705 
  706         y0 = mq_montymul(y, R2b, Qb, Q0Ib);
  707         y1 = mq_montysqr(y0, Qb, Q0Ib);
  708         y2 = mq_montymul(y1, y0, Qb, Q0Ib);
  709         y3 = mq_montymul(y2, y1, Qb, Q0Ib);
  710         y4 = mq_montysqr(y3, Qb, Q0Ib);
  711         y5 = mq_montysqr(y4, Qb, Q0Ib);
  712         y6 = mq_montysqr(y5, Qb, Q0Ib);
  713         y7 = mq_montysqr(y6, Qb, Q0Ib);
  714         y8 = mq_montysqr(y7, Qb, Q0Ib);
  715         y9 = mq_montymul(y8, y2, Qb, Q0Ib);
  716         y10 = mq_montymul(y9, y8, Qb, Q0Ib);
  717         y11 = mq_montysqr(y10, Qb, Q0Ib);
  718         y12 = mq_montysqr(y11, Qb, Q0Ib);
  719         y13 = mq_montymul(y12, y9, Qb, Q0Ib);
  720         y14 = mq_montysqr(y13, Qb, Q0Ib);
  721         y15 = mq_montysqr(y14, Qb, Q0Ib);
  722         y16 = mq_montymul(y15, y10, Qb, Q0Ib);
  723         y17 = mq_montysqr(y16, Qb, Q0Ib);
  724         y18 = mq_montymul(y17, y0, Qb, Q0Ib);
  725 
  726         /*
  727          * Final multiplication with x, which is not in Montgomery
  728          * representation, computes the correct division result.
  729          */
  730         return mq_montymul(y18, x, Qb, Q0Ib);
  731 }
  732 
  733 /*
  734  * Divide x by y modulo q = 18433.
  735  */
  736 static inline uint32_t
  737 mq_div_18433(uint32_t x, uint32_t y)
  738 {
  739         /*
  740          * We invert y by computing y^(q-2) mod q.
  741          *
  742          * We use the following addition chain for exponent e = 18431:
  743          *
  744          *   e0             = 1
  745          *   e1  = 2 * e0   = 2
  746          *   e2  = e1 + e0  = 3
  747          *   e3  = 2 * e2   = 6
  748          *   e4  = e3 + e0  = 7
  749          *   e5  = 2 * e4   = 14
  750          *   e6  = 2 * e5   = 28
  751          *   e7  = e6 + e4  = 35
  752          *   e8  = e7 + e6  = 63
  753          *   e9  = 2 * e8   = 126
  754          *   e10 = 2 * e9   = 252
  755          *   e11 = e10 + e7 = 287
  756          *   e12 = 2 * e11  = 574
  757          *   e13 = 2 * e12  = 1148
  758          *   e14 = 2 * e13  = 2296
  759          *   e15 = 2 * e14  = 4592
  760          *   e16 = 2 * e15  = 9184
  761          *   e17 = 2 * e16  = 18368
  762          *   e18 = e17 + e8 = 18431
  763          *
  764          * Additions on exponents are converted to Montgomery
  765          * multiplications. We define all intermediate results as so
  766          * many local variables, and let the C compiler work out which
  767          * must be kept around.
  768          */
  769         uint32_t y0, y1, y2, y3, y4, y5, y6, y7, y8, y9;
  770         uint32_t y10, y11, y12, y13, y14, y15, y16, y17, y18;
  771 
  772         y0  = mq_montymul(y, R2t, Qt, Q0It);          
  773         y1  = mq_montysqr(y0, Qt, Q0It);
  774         y2  = mq_montymul(y1, y0, Qt, Q0It);
  775         y3  = mq_montysqr(y2, Qt, Q0It);
  776         y4  = mq_montymul(y3, y0, Qt, Q0It);
  777         y5  = mq_montysqr(y4, Qt, Q0It);
  778         y6  = mq_montysqr(y5, Qt, Q0It);
  779         y7  = mq_montymul(y6, y4, Qt, Q0It);
  780         y8  = mq_montymul(y7, y6, Qt, Q0It);
  781         y9  = mq_montysqr(y8, Qt, Q0It);
  782         y10 = mq_montysqr(y9, Qt, Q0It);
  783         y11 = mq_montymul(y10, y7, Qt, Q0It);
  784         y12 = mq_montysqr(y11, Qt, Q0It);
  785         y13 = mq_montysqr(y12, Qt, Q0It);
  786         y14 = mq_montysqr(y13, Qt, Q0It);
  787         y15 = mq_montysqr(y14, Qt, Q0It);
  788         y16 = mq_montysqr(y15, Qt, Q0It);
  789         y17 = mq_montysqr(y16, Qt, Q0It);
  790         y18 = mq_montymul(y17, y8, Qt, Q0It);
  791 
  792         /*
  793          * Final multiplication with x, which is not in Montgomery
  794          * representation, computes the correct division result.
  795          */
  796         return mq_montymul(y18, x, Qt, Q0It);
  797 }
  798 
  799 /*
  800  * Compute NTT on a ring element, binary case.
  801  */
  802 static void
  803 mq_NTT_binary(uint16_t *a, unsigned logn)
  804 {
  805         size_t n, t, m;
  806 
  807         n = (size_t)1 << logn;
  808         t = n;
  809         for (m = 1; m < n; m <<= 1) {
  810                 size_t ht, i, j1;
  811 
  812                 ht = t >> 1;
  813                 for (i = 0, j1 = 0; i < m; i ++, j1 += t) {
  814                         size_t j, j2;
  815                         uint32_t s;
  816 
  817                         s = GMb[m + i];
  818                         j2 = j1 + ht;
  819                         for (j = j1; j < j2; j ++) {
  820                                 uint32_t u, v;
  821 
  822                                 u = a[j];
  823                                 v = mq_montymul(a[j + ht], s, Qb, Q0Ib);
  824                                 a[j] = (uint16_t)mq_add(u, v, Qb);
  825                                 a[j + ht] = (uint16_t)mq_sub(u, v, Qb);
  826                         }
  827                 }
  828                 t = ht;
  829         }
  830 }
  831 
  832 /*
  833  * Compute the inverse NTT on a ring element, binary case.
  834  */
  835 static void
  836 mq_iNTT_binary(uint16_t *a, unsigned logn)
  837 {
  838         size_t n, t, m;
  839         uint32_t ni;
  840 
  841         n = (size_t)1 << logn;
  842         t = 1;
  843         m = n;
  844         while (m > 1) {
  845                 size_t hm, dt, i, j1;
  846 
  847                 hm = m >> 1;
  848                 dt = t << 1;
  849                 for (i = 0, j1 = 0; i < hm; i ++, j1 += dt) {
  850                         size_t j, j2;
  851                         uint32_t s;
  852 
  853                         j2 = j1 + t;
  854                         s = iGMb[hm + i];
  855                         for (j = j1; j < j2; j ++) {
  856                                 uint32_t u, v, w;
  857 
  858                                 u = a[j];
  859                                 v = a[j + t];
  860                                 a[j] = (uint16_t)mq_add(u, v, Qb);
  861                                 w = mq_sub(u, v, Qb);
  862                                 a[j + t] = (uint16_t)
  863                                         mq_montymul(w, s, Qb, Q0Ib);
  864                         }
  865                 }
  866                 t = dt;
  867                 m = hm;
  868         }
  869 
  870         /*
  871          * To complete the inverse NTT, we must now divide all values by
  872          * n (the vector size). We thus need the inverse of n, i.e. we
  873          * need to divide 1 by 2 logn times. But we also want it in
  874          * Montgomery representation, i.e. we also want to multiply it
  875          * by R = 2^16. In the common case, this should be a simple right
  876          * shift. The loop below is generic and works also in corner cases;
  877          * its computation time is negligible.
  878          */
  879         ni = Rb;
  880         for (m = n; m > 1; m >>= 1) {
  881                 ni = mq_rshift1(ni, Qb);
  882         }
  883         for (m = 0; m < n; m ++) {
  884                 a[m] = (uint16_t)mq_montymul(a[m], ni, Qb, Q0Ib);
  885         }
  886 }
  887 
  888 /*
  889  * Compute NTT on a ring element, ternary case.
  890  */
  891 static void
  892 mq_NTT_ternary(uint16_t *a, unsigned logn)
  893 {
  894         size_t n, hn, u, v, t, m;
  895         uint32_t r, w;
  896 
  897         n = (size_t)3 << (logn - 1);
  898         hn = n >> 1;
  899 
  900         /*
  901          * Modulo X^2-X+1.
  902          */
  903         r = GMt_square[1];
  904         for (u = 0; u < hn; u ++) {
  905                 uint32_t a0, a1, b;
  906 
  907                 a0 = a[u];
  908                 a1 = a[u + hn];
  909                 b = mq_montymul(a1, r, Qt, Q0It);
  910                 a[u] = mq_add(a0, b, Qt);
  911                 a[u + hn] = mq_sub(mq_add(a0, a1, Qt), b, Qt);
  912         }
  913 
  914         /*
  915          * Intermediate steps for degree doubling.
  916          */
  917         t = hn;
  918         for (m = 2; t > 3; m <<= 1) {
  919                 size_t ht, u1, v1;
  920 
  921                 ht = t >> 1;
  922                 for (u1 = 0, v1 = 0; u1 < m; u1 ++, v1 += t) {
  923                         size_t v2;
  924                         uint32_t s;
  925 
  926                         s = GMt_square[m + u1];
  927                         v2 = v1 + ht;
  928                         for (v = v1; v < v2; v ++) {
  929                                 uint32_t a0, a1;
  930 
  931                                 a0 = a[v];
  932                                 a1 = a[v + ht];
  933                                 a1 = mq_montymul(a1, s, Qt, Q0It);
  934                                 a[v] = mq_add(a0, a1, Qt);
  935                                 a[v + ht] = mq_sub(a0, a1, Qt);
  936                         }
  937                 }
  938                 t = ht;
  939         }
  940 
  941         /*
  942          * Degree tripling.
  943          */
  944         w = mq_montymul(GMt_square[1], GMt_square[1], Qt, Q0It);
  945         for (u = 0, v = (size_t)1 << (logn - 1); u < n; u += 3, v ++) {
  946                 uint32_t fA, fB, fC, x, x2;
  947                 uint32_t fB0, fB1, fB2, fC0, fC1, fC2;
  948 
  949                 fA = a[u + 0];
  950                 fB = a[u + 1];
  951                 fC = a[u + 2];
  952                 x = GMt_cubic[v];
  953                 x2 = mq_montysqr(x, Qt, Q0It);
  954                 fB0 = mq_montymul(fB, x, Qt, Q0It);
  955                 fB1 = mq_montymul(fB0, w, Qt, Q0It);
  956                 fB2 = mq_montymul(fB1, w, Qt, Q0It);
  957                 fC0 = mq_montymul(fC, x2, Qt, Q0It);
  958                 fC1 = mq_montymul(fC0, w, Qt, Q0It);
  959                 fC2 = mq_montymul(fC1, w, Qt, Q0It);
  960                 a[u + 0] = mq_add(fA, mq_add(fB0, fC0, Qt), Qt);
  961                 a[u + 1] = mq_add(fA, mq_add(fB1, fC2, Qt), Qt);
  962                 a[u + 2] = mq_add(fA, mq_add(fB2, fC1, Qt), Qt);
  963         }
  964 }
  965 
  966 /*
  967  * Compute the inverse NTT on a ring element, ternary case.
  968  */
  969 static void
  970 mq_iNTT_ternary(uint16_t *a, unsigned logn)
  971 {
  972         size_t n, hn, u, v, t, m;
  973         uint32_t r, w, ni;
  974 
  975         n = (size_t)3 << (logn - 1);
  976         hn = n >> 1;
  977 
  978         /*
  979          * Dividing degree by 3.
  980          */
  981         w = mq_montymul(iGMt_square[1], iGMt_square[1], Qt, Q0It);
  982         for (u = 0, v = (size_t)1 << (logn - 1); u < n; u += 3, v ++) {
  983                 uint32_t f0, f1, f2, x, x2;
  984                 uint32_t f11, f12, f21, f22;
  985 
  986                 f0 = a[u + 0];
  987                 f1 = a[u + 1];
  988                 f2 = a[u + 2];
  989                 x = iGMt_cubic[v];
  990                 x2 = mq_montysqr(x, Qt, Q0It);
  991                 f11 = mq_montymul(f1, w, Qt, Q0It);
  992                 f12 = mq_montymul(f11, w, Qt, Q0It);
  993                 f21 = mq_montymul(f2, w, Qt, Q0It);
  994                 f22 = mq_montymul(f21, w, Qt, Q0It);
  995                 a[u + 0] = mq_add(f0, mq_add(f1, f2, Qt), Qt);
  996                 a[u + 1] = mq_montymul(x,
  997                         mq_add(f0, mq_add(f11, f22, Qt), Qt), Qt, Q0It);
  998                 a[u + 2] = mq_montymul(x2,
  999                         mq_add(f0, mq_add(f12, f21, Qt), Qt), Qt, Q0It);
 1000         }
 1001 
 1002         /*
 1003          * Intermediate steps for degree halving.
 1004          */
 1005         t = 6;
 1006         for (m = (size_t)1 << (logn - 2); t < n; m >>= 1) {
 1007                 size_t ht, u1, v1;
 1008 
 1009                 ht = t >> 1;
 1010                 for (u1 = 0, v1 = 0; u1 < m; u1 ++, v1 += t) {
 1011                         size_t v2;
 1012                         uint32_t s;
 1013 
 1014                         s = iGMt_square[m + u1];
 1015                         v2 = v1 + ht;
 1016                         for (v = v1; v < v2; v ++) {
 1017                                 uint32_t a0, a1;
 1018 
 1019                                 a0 = a[v];
 1020                                 a1 = a[v + ht];
 1021                                 a[v] = mq_add(a0, a1, Qt);
 1022                                 a[v + ht] = mq_montymul(
 1023                                         mq_sub(a0, a1, Qt), s, Qt, Q0It);
 1024                         }
 1025                 }
 1026                 t <<= 1;
 1027         }
 1028 
 1029         /*
 1030          * Modulo X^2-X+1.
 1031          */
 1032         r = iGMt_square[0];
 1033         for (u = 0; u < hn; u ++) {
 1034                 uint32_t a0, a1, b;
 1035 
 1036                 a0 = a[u];
 1037                 a1 = a[u + hn];
 1038                 b = mq_montymul(r, mq_sub(a0, a1, Qt), Qt, Q0It);
 1039                 a[u] = mq_sub(mq_add(a0, a1, Qt), b, Qt);
 1040                 a[u + hn] = mq_add(b, b, Qt);
 1041         }
 1042 
 1043         /*
 1044          * Corrective factor: all values have been (cumulatively)
 1045          * multiplied by n. Inverses of n modulo q have been precomputed.
 1046          */
 1047         ni = INVNQt[logn];
 1048         for (u = 0; u < n; u ++) {
 1049                 a[u] = mq_montymul(a[u], ni, Qt, Q0It);
 1050         }
 1051 }
 1052 
 1053 static void
 1054 mq_NTT(uint16_t *a, unsigned logn, int ternary)
 1055 {
 1056         if (ternary) {
 1057                 mq_NTT_ternary(a, logn);
 1058         } else {
 1059                 mq_NTT_binary(a, logn);
 1060         }
 1061 }
 1062 
 1063 static void
 1064 mq_iNTT(uint16_t *a, unsigned logn, int ternary)
 1065 {
 1066         if (ternary) {
 1067                 mq_iNTT_ternary(a, logn);
 1068         } else {
 1069                 mq_iNTT_binary(a, logn);
 1070         }
 1071 }
 1072 
 1073 /*
 1074  * Convert a polynomial (mod q) to Montgomery representation.
 1075  */
 1076 static void
 1077 mq_poly_tomonty(uint16_t *f, unsigned logn, int ternary)
 1078 {
 1079         size_t u, n;
 1080 
 1081         if (ternary) {
 1082                 n = (size_t)3 << (logn - 1);
 1083                 for (u = 0; u < n; u ++) {
 1084                         f[u] = (uint16_t)mq_montymul(f[u], R2t, Qt, Q0It);
 1085                 }
 1086         } else {
 1087                 n = (size_t)1 << logn;
 1088                 for (u = 0; u < n; u ++) {
 1089                         f[u] = (uint16_t)mq_montymul(f[u], R2b, Qb, Q0Ib);
 1090                 }
 1091         }
 1092 }
 1093 
 1094 /*
 1095  * Multiply two polynomials together (NTT representation, and using
 1096  * a Montgomery multiplication). Result f*g is written over f.
 1097  */
 1098 static void
 1099 mq_poly_montymul_ntt(uint16_t *f, const uint16_t *g, unsigned logn, int ternary)
 1100 {
 1101         size_t u, n;
 1102 
 1103         if (ternary) {
 1104                 n = (size_t)3 << (logn - 1);
 1105                 for (u = 0; u < n; u ++) {
 1106                         f[u] = (uint16_t)mq_montymul(f[u], g[u], Qt, Q0It);
 1107                 }
 1108         } else {
 1109                 n = (size_t)1 << logn;
 1110                 for (u = 0; u < n; u ++) {
 1111                         f[u] = (uint16_t)mq_montymul(f[u], g[u], Qb, Q0Ib);
 1112                 }
 1113         }
 1114 }
 1115 
 1116 /*
 1117  * Subtract polynomial g from polynomial f.
 1118  */
 1119 static void
 1120 mq_poly_sub(uint16_t *f, const uint16_t *g, unsigned logn, int ternary)
 1121 {
 1122         size_t u, n;
 1123 
 1124         if (ternary) {
 1125                 n = (size_t)3 << (logn - 1);
 1126                 for (u = 0; u < n; u ++) {
 1127                         f[u] = (uint16_t)mq_sub(f[u], g[u], Qt);
 1128                 }
 1129         } else {
 1130                 n = (size_t)1 << logn;
 1131                 for (u = 0; u < n; u ++) {
 1132                         f[u] = (uint16_t)mq_sub(f[u], g[u], Qb);
 1133                 }
 1134         }
 1135 }
 1136 
 1137 /* ===================================================================== */
 1138 
 1139 /*
 1140  * Falcon verification context: it contains the decoded public key,
 1141  * and the running hash context.
 1142  */
 1143 struct falcon_vrfy_ {
 1144         shake_context sc;
 1145         uint16_t h[1024];
 1146         unsigned logn;
 1147         int ternary;
 1148 };
 1149 
 1150 /* see falcon.h */
 1151 falcon_vrfy *
 1152 falcon_vrfy_new(void)
 1153 {
 1154         falcon_vrfy *fv;
 1155 
 1156         fv = malloc(sizeof *fv);
 1157         if (fv == NULL) {
 1158                 return NULL;
 1159         }
 1160         fv->logn = 0;
 1161         fv->ternary = 0;
 1162         return fv;
 1163 }
 1164 
 1165 /* see falcon.h */
 1166 void
 1167 falcon_vrfy_free(falcon_vrfy *fv)
 1168 {
 1169         if (fv != NULL) {
 1170                 free(fv);
 1171         }
 1172 }
 1173 
 1174 /* see falcon.h */
 1175 int
 1176 falcon_vrfy_set_public_key(falcon_vrfy *fv,
 1177         const void *pkey, size_t len)
 1178 {
 1179         const unsigned char *buf;
 1180         unsigned fb;
 1181 
 1182         buf = pkey;
 1183 
 1184         /*
 1185          * Read first byte: it defines the modulus and degree:
 1186          *   t 000 dddd
 1187          *   ^ ^^^ ^^^^
 1188          *   |  |    |
 1189          *   |  |    +----- degree log, over four bits (1 to 10 only)
 1190          *   |  |
 1191          *   |  |
 1192          *   |  |
 1193          *   |  +---------- reserved, must be 0
 1194          *   |
 1195          *   +------------- 1 for ternary, 0 for binary
 1196          */
 1197         if (len <= 1) {
 1198                 goto bad_pkey;
 1199         }
 1200         fb = *buf ++;
 1201         len --;
 1202         fv->logn = fb & 0x0F;
 1203         if ((fb >> 7) != 0) {
 1204                 fv->ternary = 1;
 1205                 if (fv->logn < 2 || fv->logn > 9) {
 1206                         goto bad_pkey;
 1207                 }
 1208         } else {
 1209                 fv->ternary = 0;
 1210                 if (fv->logn < 1 || fv->logn > 10) {
 1211                         goto bad_pkey;
 1212                 }
 1213         }
 1214         if (((fb >> 4) & 0x07) != 0) {
 1215                 goto bad_pkey;
 1216         }
 1217 
 1218         /*
 1219          * Decode public vector.
 1220          */
 1221         if (fv->ternary) {
 1222                 if (falcon_decode_18433(fv->h, fv->logn, buf, len) != len) {
 1223                         goto bad_pkey;
 1224                 }
 1225         } else {
 1226                 if (falcon_decode_12289(fv->h, fv->logn, buf, len) != len) {
 1227                         goto bad_pkey;
 1228                 }
 1229         }
 1230 
 1231         /*
 1232          * We apply NTT and Montgomery representation immediately, since
 1233          * it will be used for signature verification.
 1234          * (Conceptually this could be part of the public key encoding,
 1235          * but it would require making both NTT and Montgomery
 1236          * representations part of the specification.)
 1237          */
 1238         mq_NTT(fv->h, fv->logn, fv->ternary);
 1239         mq_poly_tomonty(fv->h, fv->logn, fv->ternary);
 1240         return 1;
 1241 
 1242 bad_pkey:
 1243         fv->logn = 0;
 1244         return 0;
 1245 }
 1246 
 1247 /* see falcon.h */
 1248 void
 1249 falcon_vrfy_start(falcon_vrfy *fv, const void *r, size_t rlen)
 1250 {
 1251         /*
 1252          * We use SHAKE-256, which has an internal capacity of 512 bits.
 1253          */
 1254         shake_init(&fv->sc, 512);
 1255         shake_inject(&fv->sc, r, rlen);
 1256 }
 1257 
 1258 /* see falcon.h */
 1259 void
 1260 falcon_vrfy_update(falcon_vrfy *fv, const void *data, size_t len)
 1261 {
 1262         shake_inject(&fv->sc, data, len);
 1263 }
 1264 
 1265 /* see internal.h */
 1266 int
 1267 falcon_vrfy_verify_raw(const uint16_t *c0, const int16_t *s2,
 1268         const uint16_t *h, unsigned logn, int ternary)
 1269 {
 1270         uint16_t x[1024];
 1271         size_t u, n;
 1272         uint32_t q;
 1273 
 1274         if (ternary) {
 1275                 n = (size_t)3 << (logn - 1);
 1276                 q = Qt;
 1277         } else {
 1278                 n = (size_t)1 << logn;
 1279                 q = Qb;
 1280         }
 1281 
 1282         /*
 1283          * Reduce s2 elements modulo q ([0..q-1] range).
 1284          */
 1285         for (u = 0; u < n; u ++) {
 1286                 uint32_t w;
 1287 
 1288                 w = (uint32_t)s2[u];
 1289                 w += q & -(w >> 31);
 1290                 x[u] = (uint16_t)w;
 1291         }
 1292 
 1293         /*
 1294          * Compute s1 = s2*h - c0 mod phi mod q (in x[]).
 1295          */
 1296         mq_NTT(x, logn, ternary);
 1297         mq_poly_montymul_ntt(x, h, logn, ternary);
 1298         mq_iNTT(x, logn, ternary);
 1299         mq_poly_sub(x, c0, logn, ternary);
 1300 
 1301         /*
 1302          * Normalize s1 elements into the [-q/2..q/2] range.
 1303          */
 1304         for (u = 0; u < n; u ++) {
 1305                 int32_t w;
 1306 
 1307                 w = (int32_t)x[u];
 1308                 w -= (int32_t)(q & -(((q >> 1) - (uint32_t)w) >> 31));
 1309                 ((int16_t *)x)[u] = (int16_t)w;
 1310         }
 1311 
 1312         /*
 1313          * Signature is valid if and only if the aggregate (s1,s2) vector
 1314          * is short enough.
 1315          */
 1316         return falcon_is_short((int16_t *)x, s2, logn, ternary);
 1317 }
 1318 
 1319 /* see falcon.h */
 1320 int
 1321 falcon_vrfy_verify(falcon_vrfy *fv, const void *sig, size_t len)
 1322 {
 1323         const unsigned char *sig_buf;
 1324         unsigned q;
 1325         int fb;
 1326         uint16_t c0[1024];
 1327         int16_t s2[1024];
 1328 
 1329         /*
 1330          * Public key must have been set.
 1331          */
 1332         if (fv->logn == 0) {
 1333                 return -2;
 1334         }
 1335 
 1336         /*
 1337          * Signature cannot be too short.
 1338          */
 1339         if (len <= 2) {
 1340                 return -1;
 1341         }
 1342 
 1343         /*
 1344          * Read first byte: it defines the modulus and degree:
 1345          *   t cc 0 dddd
 1346          *   ^ ^^ ^ ^^^^
 1347          *   |  | |   |
 1348          *   |  | |   +----- degree log, over four bits (1 to 10 only)
 1349          *   |  | |
 1350          *   |  | +--------- reserved, must be 0
 1351          *   |  |
 1352          *   |  +----------- compression type
 1353          *   |
 1354          *   +-------------- 1 for ternary, 0 for binary
 1355          *
 1356          * Compression is:
 1357          *   00   uncompressed, 16 bits per integer (signed)
 1358          *   01   compressed with static codes
 1359          * Other compression values are reserved.
 1360          */
 1361         sig_buf = sig;
 1362         fb = *sig_buf ++;
 1363         len --;
 1364 
 1365         /* Check reserved bit. */
 1366         if ((fb & 0x10) != 0) {
 1367                 return -1;
 1368         }
 1369 
 1370         /* Check degree. */
 1371         if ((fb & 0x0F) != fv->logn) {
 1372                 return -1;
 1373         }
 1374         if ((fb >> 7) != fv->ternary) {
 1375                 return -1;
 1376         }
 1377         q = fv->ternary ? Qt : Qb;
 1378 
 1379         /* Decode value. */
 1380         if (falcon_decode_small(s2, fv->logn,
 1381                 (fb >> 5) & 0x03, q, sig_buf, len) != len)
 1382         {
 1383                 return -1;
 1384         }
 1385 
 1386         /*
 1387          * Hash the message into a vector, then verify the signature.
 1388          */
 1389         shake_flip(&fv->sc);
 1390         falcon_hash_to_point(&fv->sc, q, c0, fv->logn);
 1391         return (int)falcon_vrfy_verify_raw(
 1392                 c0, s2, fv->h, fv->logn, fv->ternary);
 1393 }
 1394 
 1395 /* see internal.h */
 1396 int
 1397 falcon_compute_public(uint16_t *h,
 1398         const int16_t *f, const int16_t *g, unsigned logn, int ternary)
 1399 {
 1400         size_t u, n;
 1401         uint16_t t[1024];
 1402         uint32_t q;
 1403 
 1404         if (ternary) {
 1405                 n = (size_t)3 << (logn - 1);
 1406                 q = Qt;
 1407         } else {
 1408                 n = (size_t)1 << logn;
 1409                 q = Qb;
 1410         }
 1411         for (u = 0; u < n; u ++) {
 1412                 t[u] = mq_conv_small(f[u], q);
 1413                 h[u] = mq_conv_small(g[u], q);
 1414         }
 1415         mq_NTT(h, logn, ternary);
 1416         mq_NTT(t, logn, ternary);
 1417         for (u = 0; u < n; u ++) {
 1418                 if (t[u] == 0) {
 1419                         return 0;
 1420                 }
 1421                 h[u] = ternary
 1422                         ? mq_div_18433(h[u], t[u])
 1423                         : mq_div_12289(h[u], t[u]);
 1424         }
 1425         mq_iNTT(h, logn, ternary);
 1426         return 1;
 1427 }
 1428 
 1429 /* see internal.h */
 1430 int
 1431 falcon_complete_private(int16_t *G,
 1432         const int16_t *f, const int16_t *g, const int16_t *F,
 1433         unsigned logn, int ternary)
 1434 {
 1435         size_t u, n;
 1436         uint16_t t1[1024], t2[1024];
 1437         uint32_t q;
 1438 
 1439         if (ternary) {
 1440                 n = (size_t)3 << (logn - 1);
 1441                 q = Qt;
 1442         } else {
 1443                 n = (size_t)1 << logn;
 1444                 q = Qb;
 1445         }
 1446         for (u = 0; u < n; u ++) {
 1447                 t1[u] = mq_conv_small(g[u], q);
 1448                 t2[u] = mq_conv_small(F[u], q);
 1449         }
 1450         mq_NTT(t1, logn, ternary);
 1451         mq_NTT(t2, logn, ternary);
 1452         mq_poly_tomonty(t1, logn, ternary);
 1453         mq_poly_montymul_ntt(t1, t2, logn, ternary);
 1454         for (u = 0; u < n; u ++) {
 1455                 t2[u] = mq_conv_small(f[u], q);
 1456         }
 1457         mq_NTT(t2, logn, ternary);
 1458         for (u = 0; u < n; u ++) {
 1459                 if (t2[u] == 0) {
 1460                         return 0;
 1461                 }
 1462                 t1[u] = ternary
 1463                         ? mq_div_18433(t1[u], t2[u])
 1464                         : mq_div_12289(t1[u], t2[u]);
 1465         }
 1466         mq_iNTT(t1, logn, ternary);
 1467         for (u = 0; u < n; u ++) {
 1468                 uint32_t w;
 1469 
 1470                 w = t1[u];
 1471                 G[u] = (int32_t)w - (int32_t)(q & ~-((w - (q >> 1)) >> 31));
 1472         }
 1473         return 1;
 1474 }