package org.bouncycastle.pqc.crypto.falcon;

import org.bouncycastle.crypto.digests.SHAKEDigest;
import org.bouncycastle.util.Pack;

class FalconKeyGen
{
    //    FPREngine fpr;
    //FalconSmallPrimeList primes;
    //FalconFFT fft;
//    FalconCodec codec;
//    FalconVrfy vrfy;

    FalconKeyGen()
    {
//        this.fpr = new FPREngine();
        //this.primes = new FalconSmallPrimeList();
        //this.fft = new FalconFFT();
//        this.codec = new FalconCodec();
//        this.vrfy = new FalconVrfy();
    }

    private static int mkn(int logn)
    {
        return 1 << logn;
    }

    /*
     * Reduce a small signed integer modulo a small prime. The source
     * value x MUST be such that -p < x < p.
     */
    private static int modp_set(int x, int p)
    {
        int w = x;
        w += p & -(w >>> 31);
        return w;
    }

    /*
     * Normalize a modular integer around 0.
     */
    private static int modp_norm(int x, int p)
    {
        return (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.
     */
    private static int modp_ninv31(int p)
    {
        int y = 2 - p;
        y *= 2 - p * y;
        y *= 2 - p * y;
        y *= 2 - p * y;
        y *= 2 - p * y;
        return 0x7FFFFFFF & -y;
    }

    /*
     * Compute R = 2^31 mod p.
     */
    private static int modp_R(int p)
    {
        /*
         * Since 2^30 < p < 2^31, we know that 2^31 mod p is simply
         * 2^31 - p.
         */
        return (1 << 31) - p;
    }

    /*
     * Addition modulo p.
     */
    private static int modp_add(int a, int b, int p)
    {
        int d = a + b - p;
        d += p & -(d >>> 31);
        return d;
    }

    /*
     * Subtraction modulo p.
     */
    private static int modp_sub(int a, int b, int p)
    {
        int d = a - b;
        d += p & -(d >>> 31);
        return d;
    }

    /*
     * Montgomery multiplication modulo p. The 'p0i' value is -1/p mod 2^31.
     * It is required that p is an odd integer.
     */
    private static int modp_montymul(int a, int b, int p, int p0i)
    {
        long z, w;
        int d;

        z = toUnsignedLong(a) * toUnsignedLong(b);
        w = ((z * p0i) & 0x7FFFFFFFL) * p;
        d = (int)((z + w) >>> 31) - p;
        d += p & -(d >>> 31);
        return d;
    }

    /*
     * Compute R2 = 2^62 mod p.
     */
    private static int modp_R2(int p, int p0i)
    {
        int 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.
     */
    private static int modp_Rx(int x, int p, int p0i, int R2)
    {
        int i;
        int 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; (1 << i) <= x; i++)
        {
            if ((x & (1 << 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
     */
    private static int modp_div(int a, int b, int p, int p0i, int R)
    {
        int z, e;
        int i;

        e = p - 2;
        z = R;
        for (i = 30; i >= 0; i--)
        {
            int z2;

            z = modp_montymul(z, z, p, p0i);
            z2 = modp_montymul(z, b, p, p0i);
            z ^= (z ^ z2) & -(int)((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.
     */
    private static final short[] 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.
     */
    private static void modp_mkgm2(int[] srcgm, int gm, int[] srcigm, int igm, int logn,
                                   int g, int p, int p0i)
    {
        int u, n;
        int k;
        int ig, x1, x2, R2;

        n = mkn(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++)
        {
            int v;

            v = REV10[u << k];
            srcgm[gm + v] = x1;
            srcigm[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]...
     */
    private static void modp_NTT2_ext(int[] srca, int a, int stride, int[] srcgm, int gm, int logn,
                                      int p, int p0i)
    {
        int t, m, n;

        if (logn == 0)
        {
            return;
        }
        n = mkn(logn);
        t = n;
        for (m = 1; m < n; m <<= 1)
        {
            int ht, u, v1;

            ht = t >> 1;
            for (u = 0, v1 = 0; u < m; u++, v1 += t)
            {
                int s;
                int v;
                int r1, r2;

                s = srcgm[gm + m + u];
                r1 = a + v1 * stride;
                r2 = r1 + ht * stride;
                for (v = 0; v < ht; v++, r1 += stride, r2 += stride)
                {
                    int x, y;

                    x = srca[r1];
                    y = modp_montymul(srca[r2], s, p, p0i);
                    srca[r1] = modp_add(x, y, p);
                    srca[r2] = modp_sub(x, y, p);
                }
            }
            t = ht;
        }
    }

    /*
     * Compute the inverse NTT over a polynomial (binary case).
     */
    private static void modp_iNTT2_ext(int[] srca, int a, int stride, int[] srcigm, int igm, int logn,
                                       int p, int p0i)
    {
        int t, m, n, k;
        int ni;
        int r;

        if (logn == 0)
        {
            return;
        }
        n = mkn(logn);
        t = 1;
        for (m = n; m > 1; m >>= 1)
        {
            int hm, dt, u, v1;

            hm = m >> 1;
            dt = t << 1;
            for (u = 0, v1 = 0; u < hm; u++, v1 += dt)
            {
                int s;
                int v;
                int r1, r2;

                s = srcigm[igm + hm + u];
                r1 = a + v1 * stride;
                r2 = r1 + t * stride;
                for (v = 0; v < t; v++, r1 += stride, r2 += stride)
                {
                    int x, y;

                    x = srca[r1];
                    y = srca[r2];
                    srca[r1] = modp_add(x, y, p);
                    srca[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 = 1 << (31 - logn);
        for (k = 0, r = a; k < n; k++, r += stride)
        {
            srca[r] = modp_montymul(srca[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)
    private static void modp_NTT2(int[] srca, int a, int[] srcgm, int gm, int logn, int p, int p0i)
    {
        modp_NTT2_ext(srca, a, 1, srcgm, gm, logn, p, p0i);
    }

    //    #define modp_iNTT2(a, igm, logn, p, p0i) modp_iNTT2_ext(a, 1, igm, logn, p, p0i)
    private static void modp_iNTT2(int[] srca, int a, int[] srcigm, int igm, int logn, int p, int p0i)
    {
        modp_iNTT2_ext(srca, a, 1, srcigm, 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().
     */
    private static void modp_poly_rec_res(int[] srcf, int f, int logn,
                                          int p, int p0i, int R2)
    {
        int hn, u;

        hn = 1 << (logn - 1);
        for (u = 0; u < hn; u++)
        {
            int w0, w1;

            w0 = srcf[f + (u << 1)];
            w1 = srcf[f + (u << 1) + 1];
            srcf[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.
     */
    private static void zint_sub(int[] srca, int a, int[] srcb, int b, int len,
                                 int ctl)
    {
        int u;
        int cc, m;
        int aw, w, au;
        cc = 0;
        m = -ctl;
        for (u = 0; u < len; u++)
        {
            au = a + u;
            aw = srca[au];
            w = aw - srcb[b + u] - cc;
            cc = w >>> 31;
            aw ^= ((w & 0x7FFFFFFF) ^ aw) & m;
            srca[au] = aw;
        }
    }

    /*
     * Mutiply the provided big integer m with a small value x.
     * This function assumes that x < 2^31. The carry word is returned.
     */
    private static int zint_mul_small(int[] srcm, int m, int mlen, int x)
    {
        int u;
        int cc;

        cc = 0;
        for (u = 0; u < mlen; u++)
        {
            long z;

            z = toUnsignedLong(srcm[m + u]) * toUnsignedLong(x) + cc;
            srcm[m + u] = (int)z & 0x7FFFFFFF;
            cc = (int)(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
     */
    private static int zint_mod_small_unsigned(int[] srcd, int d, int dlen,
                                               int p, int p0i, int R2)
    {
        /*
         * Algorithm: we inject words one by one, starting with the high
         * word. Each step is:
         *  - multiply x by 2^31
         *  - add new word
         */
        int x = 0;
        int u = dlen;
        while (u-- > 0)
        {
            int w;

            x = modp_montymul(x, R2, p, p0i);
            w = srcd[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.
     */
    private static int zint_mod_small_signed(int[] srcd, int d, int dlen,
                                             int p, int p0i, int R2, int Rx)
    {
        int z;

        if (dlen == 0)
        {
            return 0;
        }
        z = zint_mod_small_unsigned(srcd, d, dlen, p, p0i, R2);
        z = modp_sub(z, Rx & -(srcd[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.
     */
    private static void zint_add_mul_small(int[] srcx, int x,
                                           int[] srcy, int y, int len, int s)
    {
        int u;
        int cc;

        cc = 0;
        for (u = 0; u < len; u++)
        {
            int xw, yw;
            long z;

            xw = srcx[x + u];
            yw = srcy[y + u];
            z = toUnsignedLong(yw) * toUnsignedLong(s) + toUnsignedLong(xw) + toUnsignedLong(cc);
            srcx[x + u] = (int)z & 0x7FFFFFFF;
            cc = (int)(z >>> 31);
        }
        srcx[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.
     */
    private static void zint_norm_zero(int[] srcx, int x, int[] srcp, int p, int len)
    {
        int u;
        int 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)
        {
            int wx, wp, cc;

            /*
             * Get the two words to compare in wx and wp (both over
             * 31 bits exactly).
             */
            wx = srcx[x + u];
            wp = (srcp[p + u] >>> 1) | (bb << 30);
            bb = srcp[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(srcx, x, srcp, 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.
     */
    private static void zint_rebuild_CRT(int[] srcxx, int xx, int xlen, int xstride,
                                         int num, int normalize_signed,
                                         int[] srctmp, int tmp)
    {
        int u;
        int x;

        srctmp[tmp] = FalconSmallPrimeList.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.
             */
            int p, p0i, s, R2;
            int v;

            p = FalconSmallPrimeList.PRIMES[u].p;
            s = FalconSmallPrimeList.PRIMES[u].s;
            p0i = modp_ninv31(p);
            R2 = modp_R2(p, p0i);

            for (v = 0, x = xx; v < num; v++, x += xstride)
            {
                int xp, xq, xr;
                /*
                 * xp = the integer x modulo the prime p for this
                 *      iteration
                 * xq = (x mod q) mod p
                 */
                xp = srcxx[x + u];
                xq = zint_mod_small_unsigned(srcxx, 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(srcxx, x, srctmp, tmp, u, xr);
            }

            /*
             * Update product of primes in tmp[].
             */
            srctmp[tmp + u] = zint_mul_small(srctmp, tmp, u, p);
        }

        /*
         * Normalize the reconstructed values around 0.
         */
        if (normalize_signed != 0)
        {
            for (u = 0, x = xx; u < num; u++, x += xstride)
            {
                zint_norm_zero(srcxx, x, srctmp, 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.
     */
    private static void zint_negate(int[] srca, int a, int len, int ctl)
    {
        int u;
        int 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++)
        {
            int aw;

            aw = srca[a + u];
            aw = (aw ^ m) + cc;
            srca[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.
     */
    private static int zint_co_reduce(int[] srca, int a, int[] srcb, int b, int len,
                                      long xa, long xb, long ya, long yb)
    {
        int u;
        long cca, ccb;
        int nega, negb;

        cca = 0;
        ccb = 0;
        for (u = 0; u < len; u++)
        {
            int wa, wb;
            long za, zb;

            wa = srca[a + u];
            wb = srcb[b + u];
            za = wa * xa + wb * xb + cca;
            zb = wa * ya + wb * yb + ccb;
            if (u > 0)
            {
                srca[a + u - 1] = (int)za & 0x7FFFFFFF;
                srcb[b + u - 1] = (int)zb & 0x7FFFFFFF;
            }
//            cca = *(int64_t *)&za >> 31;
            cca = za >> 31;
//            ccb = *(int64_t *)&zb >> 31;
            ccb = zb >> 31;
        }
        srca[a + len - 1] = (int)cca;
        srcb[b + len - 1] = (int)ccb;

        nega = (int)(cca >>> 63);
        negb = (int)(ccb >>> 63);
        zint_negate(srca, a, len, nega);
        zint_negate(srcb, 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.
     */
    private static void zint_finish_mod(int[] srca, int a, int len, int[] srcm, int m, int neg)
    {
        int u;
        int 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 = (srca[a + u] - srcm[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++)
        {
            int aw, mw;

            aw = srca[a + u];
            mw = (srcm[m + u] ^ xm) & ym;
            aw = aw - mw - cc;
            srca[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.
     */
    private static void zint_co_reduce_mod(int[] srca, int a, int[] srcb, int b, int[] srcm, int m, int len,
                                           int m0i, long xa, long xb, long ya, long yb)
    {
        int u;
        long cca, ccb;
        int fa, fb;

        /*
         * These are actually four combined Montgomery multiplications.
         */
        cca = 0;
        ccb = 0;
        fa = ((srca[a] * (int)xa + srcb[b] * (int)xb) * m0i) & 0x7FFFFFFF;
        fb = ((srca[a] * (int)ya + srcb[b] * (int)yb) * m0i) & 0x7FFFFFFF;
        for (u = 0; u < len; u++)
        {
            int wa, wb;
            long za, zb;

            wa = srca[a + u];
            wb = srcb[b + u];
            za = wa * xa + wb * xb
                + srcm[m + u] * toUnsignedLong(fa) + cca;
            zb = wa * ya + wb * yb
                + srcm[m + u] * toUnsignedLong(fb) + ccb;
            if (u > 0)
            {
                srca[a + u - 1] = (int)za & 0x7FFFFFFF;
                srcb[b + u - 1] = (int)zb & 0x7FFFFFFF;
            }
            cca = za >> 31;
            ccb = zb >> 31;
        }
        srca[a + len - 1] = (int)cca;
        srcb[b + len - 1] = (int)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(srca, a, len, srcm, m, (int)(cca >>> 63));
        zint_finish_mod(srcb, b, len, srcm, m, (int)(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.
     */
    private static int zint_bezout(int[] srcu, int u, int[] srcv, int v,
                                   int[] srcx, int x, int[] srcy, int y,
                                   int len, int[] srctmp, int 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.
         */
        int u0, u1, v0, v1, a, b;
        int x0i, y0i;
        int num, rc;
        int 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(srcx[x]);
        y0i = modp_ninv31(srcy[y]);

        /*
         * 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);
        System.arraycopy(srcx, x, srctmp, a, len);
        // memcpy(b, y, len * sizeof *y);
        System.arraycopy(srcy, y, srctmp, b, len);
        // u0[0] = 1;
        srcu[u0] = 1;
        // memset(u0 + 1, 0, (len - 1) * sizeof *u0);
        // memset(v0, 0, len * sizeof *v0);
        srcv[v0] = 0;
        for (int i = 1; i < len; i++)
        {
            srcu[u0 + i] = 0;
            srcv[v0 + i] = 0;
        }
        // memcpy(u1, y, len * sizeof *u1);
        System.arraycopy(srcy, y, srctmp, u1, len);
        // memcpy(v1, x, len * sizeof *v1);
        System.arraycopy(srcx, x, srctmp, v1, len);
        // v1[0] --;
        srctmp[v1]--;
        /*
         * 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 * len + 30; num >= 30; num -= 30)
        {
            int c0, c1;
            int a0, a1, b0, b1;
            long a_hi, b_hi;
            int a_lo, b_lo;
            long pa, pb, qa, qb;
            int i;
            int 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 = -1;
            c1 = -1;
            a0 = 0;
            a1 = 0;
            b0 = 0;
            b1 = 0;
            j = len;
            while (j-- > 0)
            {
                int aw, bw;

                aw = srctmp[a + j];
                bw = srctmp[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) - 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 = (toUnsignedLong(a0) << 31) + toUnsignedLong(a1);
            b_hi = (toUnsignedLong(b0) << 31) + toUnsignedLong(b1);
            a_lo = srctmp[a];
            b_lo = srctmp[b];

            /*
             * 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.
                 */
                int rt, oa, ob, cAB, cBA, cA;
                long rz;

                /*
                 * rt = 1 if a_hi > b_hi, 0 otherwise.
                 */
                rz = b_hi - a_hi;
                rt = (int)((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 & -toUnsignedLong(cAB);
                pa -= qa & -(long)cAB;
                pb -= qb & -(long)cAB;
                b_lo -= a_lo & -cBA;
                b_hi -= a_hi & -toUnsignedLong(cBA);
                qa -= pa & -(long)cBA;
                qb -= pb & -(long)cBA;

                /*
                 * Shifting.
                 */
                a_lo += a_lo & (cA - 1);
                pa += pa & ((long)cA - 1);
                pb += pb & ((long)cA - 1);
                a_hi ^= (a_hi ^ (a_hi >> 1)) & -toUnsignedLong(cA);
                b_lo += b_lo & -cA;
                qa += qa & -(long)cA;
                qb += qb & -(long)cA;
                b_hi ^= (b_hi ^ (b_hi >> 1)) & (toUnsignedLong(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(srctmp, a, srctmp, b, len, pa, pb, qa, qb);
            pa -= (pa + pa) & -(long)(r & 1);
            pb -= (pb + pb) & -(long)(r & 1);
            qa -= (qa + qa) & -(long)(r >>> 1);
            qb -= (qb + qb) & -(long)(r >>> 1);
            zint_co_reduce_mod(srcu, u0, srctmp, u1, srcy, y, len, y0i, pa, pb, qa, qb);
            zint_co_reduce_mod(srcv, v0, srctmp, v1, srcx, 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 = srctmp[a] ^ 1;
        for (j = 1; j < len; j++)
        {
            rc |= srctmp[a + j];
        }
        return ((1 - ((rc | -rc) >>> 31)) & srcx[x] & srcy[y]);
    }

    /*
     * 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.
     */
    private static void zint_add_scaled_mul_small(int[] srcx, int x, int xlen,
                                                  int[] srcy, int y, int ylen, int k,
                                                  int sch, int scl)
    {
        int u;
        int ysign, tw;
        int cc;

        if (ylen == 0)
        {
            return;
        }

        ysign = -(srcy[y + ylen - 1] >>> 30) >>> 1;
        tw = 0;
        cc = 0;
        for (u = sch; u < xlen; u++)
        {
            int v;
            int wy, wys, ccu;
            long z;

            /*
             * Get the next word of y (scaled).
             */
            v = u - sch;
            wy = v < ylen ? srcy[y + v] : ysign;
            wys = ((wy << scl) & 0x7FFFFFFF) | tw;
            tw = wy >>> (31 - scl);

            /*
             * The expression below does not overflow.
             */
            z = (toUnsignedLong(wys) * (long)k + toUnsignedLong(srcx[x + u]) + cc);
            srcx[x + u] = (int)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 = (int)(z >>> 31);
            cc = 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.
     */
    private static void zint_sub_scaled(int[] srcx, int x, int xlen,
                                        int[] srcy, int y, int ylen, int sch, int scl)
    {
        int u;
        int ysign, tw;
        int cc;

        if (ylen == 0)
        {
            return;
        }

        ysign = -(srcy[y + ylen - 1] >>> 30) >>> 1;
        tw = 0;
        cc = 0;
        for (u = sch; u < xlen; u++)
        {
            int v;
            int w, wy, wys;

            /*
             * Get the next word of y (scaled).
             */
            v = u - sch;
            wy = v < ylen ? srcy[y + v] : ysign;
            wys = ((wy << scl) & 0x7FFFFFFF) | tw;
            tw = wy >>> (31 - scl);

            w = srcx[x + u] - wys - cc;
            srcx[x + u] = w & 0x7FFFFFFF;
            cc = w >>> 31;
        }
    }

    /*
     * Convert a one-word signed big integer into a signed value.
     */
    private static int zint_one_to_plain(int[] srcx, int x)
    {
        int w;

        w = srcx[x];
        w |= (w & 0x40000000) << 1;
        return 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.
     */
    private static void poly_big_to_fp(double[] srcd, int doff, int[] srcf, int f, int flen, int fstride,
                                       int logn)
    {
        int n, u;

        n = mkn(logn);
        if (flen == 0)
        {
            for (u = 0; u < n; u++)
            {
                srcd[doff + u] = FPREngine.fpr_zero;
            }
            return;
        }
        for (u = 0; u < n; u++, f += fstride)
        {
            int v;
            int neg, cc, xm;
            double x, fsc;

            /*
             * Get sign of the integer; if it is negative, then we
             * will load its absolute value instead, and negate the
             * result.
             */
            neg = -(srcf[f + flen - 1] >>> 30);
            xm = neg >>> 1;
            cc = neg & 1;
            x = FPREngine.fpr_zero;
            fsc = FPREngine.fpr_one;
            for (v = 0; v < flen; v++, fsc *= FPREngine.fpr_ptwo31)
            {
                int w = (srcf[f + v] ^ xm) + cc;
                cc = w >>> 31;
                w &= 0x7FFFFFFF;
                w -= (w << 1) & neg;
                x += w * fsc;
            }
            srcd[doff + 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.
     */
    private static int poly_big_to_small(byte[] srcd, int d, int[] srcs, int s, int lim, int logn)
    {
        int n, u;

        n = mkn(logn);
        for (u = 0; u < n; u++)
        {
            int z;

            z = zint_one_to_plain(srcs, s + u);
            if (z < -lim || z > lim)
            {
                return 0;
            }
            srcd[d + u] = (byte)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.
     */
    private static void poly_sub_scaled(int[] srcF, int F, int Flen, int Fstride,
                                        int[] srcf, int f, int flen, int fstride,
                                        int[] srck, int sch, int scl, int logn)
    {
        int n, u;

        n = mkn(logn);
        for (u = 0; u < n; u++)
        {
            int kf;
            int v;
            int x;
            int y;

            kf = -srck[u];
            x = F + u * Fstride;
            y = f;
            for (v = 0; v < n; v++)
            {
                zint_add_scaled_mul_small(srcF,
                    x, Flen, srcf, 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.
     */
    private static void poly_sub_scaled_ntt(int[] srcF, int F, int Flen, int Fstride,
                                            int[] srcf, int f, int flen, int fstride,
                                            int[] srck, int sch, int scl, int logn,
                                            int[] srctmp, int tmp)
    {
        int gm, igm, fk, t1, x;
        int y;
        int n, u, tlen;
//        FalconSmallPrime[] primes;

        n = mkn(logn);
        tlen = flen + 1;
        gm = tmp;
        igm = gm + mkn(logn);
        fk = igm + mkn(logn);
        t1 = fk + n * tlen;

//        primes = this.primes.PRIMES;

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

            p = FalconSmallPrimeList.PRIMES[u].p;
            p0i = modp_ninv31(p);
            R2 = modp_R2(p, p0i);
            Rx = modp_Rx(flen, p, p0i, R2);
            modp_mkgm2(srctmp, gm, srctmp, igm, logn, FalconSmallPrimeList.PRIMES[u].g, p, p0i);

            for (v = 0; v < n; v++)
            {
                srctmp[t1 + v] = modp_set(srck[v], p);
            }
            modp_NTT2(srctmp, t1, srctmp, gm, logn, p, p0i);
            for (v = 0, y = f, x = fk + u;
                 v < n; v++, y += fstride, x += tlen)
            {
                srctmp[x] = zint_mod_small_signed(srcf, y, flen, p, p0i, R2, Rx);
            }
            modp_NTT2_ext(srctmp, fk + u, tlen, srctmp, gm, logn, p, p0i);
            for (v = 0, x = fk + u; v < n; v++, x += tlen)
            {
                srctmp[x] = modp_montymul(
                    modp_montymul(srctmp[t1 + v], srctmp[x], p, p0i), R2, p, p0i);
            }
            modp_iNTT2_ext(srctmp, fk + u, tlen, srctmp, igm, logn, p, p0i);
        }

        /*
         * Rebuild k*f.
         */
        zint_rebuild_CRT(srctmp, fk, tlen, tlen, n, 1, srctmp, t1);

        /*
         * Subtract k*f, scaled, from F.
         */
        for (u = 0, x = F, y = fk; u < n; u++, x += Fstride, y += tlen)
        {
            zint_sub_scaled(srcF, x, Flen, srctmp, y, tlen, sch, scl);
        }
    }

    /* ==================================================================== */

    /*
     * 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.
     *
     * The caller supplies an 8-byte scratch buffer; mkgauss reads f and g
     * coefficient-by-coefficient (2 calls per Gaussian sample, ~16K calls
     * per Falcon-1024 keygen attempt) so a per-call allocation here would
     * be the dominant byte[] allocator in the keygen path.
     */
    private static long get_rng_u64(SHAKEDigest rng, byte[] tmp8)
    {
        rng.doOutput(tmp8, 0, 8);
        return Pack.littleEndianToLong(tmp8, 0);
    }


    /*
     * 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.
     */
    private static final long[] gauss_1024_12289 = {
        1283868770400643928L, 6416574995475331444L, 4078260278032692663L,
        2353523259288686585L, 1227179971273316331L, 575931623374121527L,
        242543240509105209L, 91437049221049666L, 30799446349977173L,
        9255276791179340L, 2478152334826140L, 590642893610164L,
        125206034929641L, 23590435911403L, 3948334035941L,
        586753615614L, 77391054539L, 9056793210L,
        940121950L, 86539696L, 7062824L,
        510971L, 32764L, 1862L,
        94L, 4L, 0L
    };

    /*
     * 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.
     */
    private static int mkgauss(SHAKEDigest rng, int logn, byte[] rngBuf8)
    {
        int u, g;
        int val;

        g = 1 << (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).
             */
            long r;
            int 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, rngBuf8);
            neg = (int)(r >>> 63);
            r &= ~(1L << 63);
            f = (int)((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, rngBuf8);
            r &= ~(1L << 63);
            for (k = 1; k < gauss_1024_12289.length; k++)
            {
                int t;

                t = (int)((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 += 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.
     */

    private static final int[] MAX_BL_SMALL = {
        1, 1, 2, 2, 4, 7, 14, 27, 53, 106, 209
    };

    private static final int[] 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.
     */
    private static final int[] bitlength_avg = {
        4,
        11,
        24,
        50,
        102,
        202,
        401,
        794,
        1577,
        3138,
        6308
    };
    private static final int[] bitlength_std = {
        0,
        1,
        1,
        1,
        1,
        2,
        4,
        5,
        8,
        13,
        25
    };

    /*
     * Minimal recursion depth at which we rebuild intermediate values
     * when reconstructing f and g.
     */
    private static final int 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.
     */
    private static int poly_small_sqnorm(byte[] srcf, int logn)
    {
        int n, u;
        int s, ng;

        n = mkn(logn);
        s = 0;
        ng = 0;
        for (u = 0; u < n; u++)
        {
            int z;

            z = srcf[u];
            s += (z * z);
            ng |= s;
        }
        return s | -(ng >>> 31);
    }

    /*
     * Convert a small vector to floating point.
     */
    private static void poly_small_to_fp(double[] srcx, int x, byte[] srcf, int logn)
    {
        int n, u;

        n = mkn(logn);
        for (u = 0; u < n; u++)
        {
            srcx[x + u] = srcf[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.
     */
    private static void make_fg_step(int[] srcdata, int data, int logn, int depth,
                                     int in_ntt, int out_ntt)
    {
        int n, hn, u;
        int slen, tlen;
        int fd, gd, fs, gs, gm, igm, t1;
        //FalconSmallPrime[] primes;

        n = 1 << logn;
        hn = n >> 1;
        slen = MAX_BL_SMALL[depth];
        tlen = MAX_BL_SMALL[depth + 1];
        //primes = FalconSmallPrimeList.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);
        System.arraycopy(srcdata, data, srcdata, fs, 2 * n * slen);

        /*
         * First slen words: we use the input values directly, and apply
         * inverse NTT as we go.
         */
        for (u = 0; u < slen; u++)
        {
            int p, p0i, R2;
            int v;
            int x;

            p = FalconSmallPrimeList.PRIMES[u].p;
            p0i = modp_ninv31(p);
            R2 = modp_R2(p, p0i);
            modp_mkgm2(srcdata, gm, srcdata, igm, logn, FalconSmallPrimeList.PRIMES[u].g, p, p0i);

            for (v = 0, x = fs + u; v < n; v++, x += slen)
            {
                srcdata[t1 + v] = srcdata[x];
            }
            if (in_ntt == 0)
            {
                modp_NTT2(srcdata, t1, srcdata, gm, logn, p, p0i);
            }
            for (v = 0, x = fd + u; v < hn; v++, x += tlen)
            {
                int w0, w1;

                w0 = srcdata[t1 + (v << 1)];
                w1 = srcdata[t1 + (v << 1) + 1];
                srcdata[x] = modp_montymul(
                    modp_montymul(w0, w1, p, p0i), R2, p, p0i);
            }
            if (in_ntt != 0)
            {
                modp_iNTT2_ext(srcdata, fs + u, slen, srcdata, igm, logn, p, p0i);
            }

            for (v = 0, x = gs + u; v < n; v++, x += slen)
            {
                srcdata[t1 + v] = srcdata[x];
            }
            if (in_ntt == 0)
            {
                modp_NTT2(srcdata, t1, srcdata, gm, logn, p, p0i);
            }
            for (v = 0, x = gd + u; v < hn; v++, x += tlen)
            {
                int w0, w1;

                w0 = srcdata[t1 + (v << 1)];
                w1 = srcdata[t1 + (v << 1) + 1];
                srcdata[x] = modp_montymul(
                    modp_montymul(w0, w1, p, p0i), R2, p, p0i);
            }
            if (in_ntt != 0)
            {
                modp_iNTT2_ext(srcdata, gs + u, slen, srcdata, igm, logn, p, p0i);
            }

            if (out_ntt == 0)
            {
                modp_iNTT2_ext(srcdata, fd + u, tlen, srcdata, igm, logn - 1, p, p0i);
                modp_iNTT2_ext(srcdata, gd + u, tlen, srcdata, 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(srcdata, fs, slen, slen, n, 1, srcdata, gm);
        zint_rebuild_CRT(srcdata, gs, slen, slen, n, 1, srcdata, gm);

        /*
         * Remaining words: use modular reductions to extract the values.
         */
        for (u = slen; u < tlen; u++)
        {
            int p, p0i, R2, Rx;
            int v;
            int x;

            p = FalconSmallPrimeList.PRIMES[u].p;
            p0i = modp_ninv31(p);
            R2 = modp_R2(p, p0i);
            Rx = modp_Rx(slen, p, p0i, R2);
            modp_mkgm2(srcdata, gm, srcdata, igm, logn, FalconSmallPrimeList.PRIMES[u].g, p, p0i);
            for (v = 0, x = fs; v < n; v++, x += slen)
            {
                srcdata[t1 + v] = zint_mod_small_signed(srcdata, x, slen, p, p0i, R2, Rx);
            }
            modp_NTT2(srcdata, t1, srcdata, gm, logn, p, p0i);
            for (v = 0, x = fd + u; v < hn; v++, x += tlen)
            {
                int w0, w1;

                w0 = srcdata[t1 + (v << 1)];
                w1 = srcdata[t1 + (v << 1) + 1];
                srcdata[x] = modp_montymul(
                    modp_montymul(w0, w1, p, p0i), R2, p, p0i);
            }
            for (v = 0, x = gs; v < n; v++, x += slen)
            {
                srcdata[t1 + v] = zint_mod_small_signed(srcdata, x, slen, p, p0i, R2, Rx);
            }
            modp_NTT2(srcdata, t1, srcdata, gm, logn, p, p0i);
            for (v = 0, x = gd + u; v < hn; v++, x += tlen)
            {
                int w0, w1;

                w0 = srcdata[t1 + (v << 1)];
                w1 = srcdata[t1 + (v << 1) + 1];
                srcdata[x] = modp_montymul(
                    modp_montymul(w0, w1, p, p0i), R2, p, p0i);
            }

            if (out_ntt == 0)
            {
                modp_iNTT2_ext(srcdata, fd + u, tlen, srcdata, igm, logn - 1, p, p0i);
                modp_iNTT2_ext(srcdata, gd + u, tlen, srcdata, igm, logn - 1, p, p0i);
            }
        }
    }

    /* d 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).
     */
    private static void make_fg(int[] srcdata, int data, byte[] srcf, byte[] srcg,
                                int logn, int depth, int out_ntt)
    {
        int n, u;
        int ft, gt, p0;
        int d;


        n = mkn(logn);
        ft = data;
        gt = ft + n;

        p0 = FalconSmallPrimeList.PRIMES[0].p;
        for (u = 0; u < n; u++)
        {
            srcdata[ft + u] = modp_set(srcf[u], p0);
            srcdata[gt + u] = modp_set(srcg[u], p0);
        }

        if (depth == 0 && out_ntt != 0)
        {
            int gm, igm;
            int p, p0i;

            p = FalconSmallPrimeList.PRIMES[0].p;
            p0i = modp_ninv31(p);
            gm = gt + n;
            igm = gm + n;
            modp_mkgm2(srcdata, gm, srcdata, igm, logn, FalconSmallPrimeList.PRIMES[0].g, p, p0i);
            modp_NTT2(srcdata, ft, srcdata, gm, logn, p, p0i);
            modp_NTT2(srcdata, gt, srcdata, gm, logn, p, p0i);
            return;
        }

        for (d = 0; d < depth; d++)
        {
            make_fg_step(srcdata, data, logn - d, d,
                d != 0 ? 1 : 0, ((d + 1) < depth || out_ntt != 0) ? 1 : 0);
        }
    }

    /*
     * 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.
     */
    private static int solve_NTRU_deepest(int logn_top,
                                          byte[] srcf, byte[] srcg, int[] srctmp)
    {
        int len;
        int Fp, Gp, fp, gp, t1, q;
        //FalconSmallPrime[] primes;

        len = MAX_BL_SMALL[logn_top];
        //primes = this.primes.PRIMES;

        Fp = 0;
        Gp = Fp + len;
        fp = Gp + len;
        gp = fp + len;
        t1 = gp + len;

        make_fg(srctmp, fp, srcf, srcg, 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(srctmp, fp, len, len, 2, 0, srctmp, 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(srctmp, Gp, srctmp, Fp, srctmp, fp, srctmp, gp, len, srctmp, t1) == 0)
        {
            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(srctmp, Fp, len, q) != 0
            || zint_mul_small(srctmp, 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.
     */
    private static int solve_NTRU_intermediate(int logn_top,
                                               byte[] srcf, byte[] srcg, int depth, int[] srctmp,
                                               double[] fprBuf, int[] kBuf)
    {
        /*
         * 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.
         */
        int logn;
        int n, hn, slen, dlen, llen, rlen, FGlen, u;
        int Fd, Gd, Ft, Gt, ft, gt, t1;
        // rt1..rt5 are offsets into the caller-supplied fprBuf scratch,
        // so the recursion does not re-allocate these doubles per call.
        // rt1..rt4 are size n; rt5 is size n/2.
        int rt1, rt2, rt3, rt4, rt5;
        int scale_fg, minbl_fg, maxbl_fg, maxbl_FG, scale_k;
        int x, y;
        // k is the caller-supplied kBuf scratch (size n).
        int[] k;
        //FalconSmallPrime[] primes;

        logn = logn_top - depth;
        n = 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 = this.primes.PRIMES;

        /*
         * Fd and Gd are the F and G from the deeper level.
         */
        Fd = 0;
        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(srctmp, ft, srcf, srcg, logn_top, depth, 1);

        /*
         * Move the newly computed f and g to make room for our candidate
         * F and G (unreduced).
         */
        Ft = 0;
        Gt = Ft + n * llen;
        t1 = Gt + n * llen;
//        memmove(t1, ft, 2 * n * slen * sizeof *ft);
        int tmp = n * slen;
        System.arraycopy(srctmp, ft, srctmp, t1, tmp + tmp);
        ft = t1;
        gt = ft + tmp;
        t1 = gt + tmp;

        /*
         * Move Fd and Gd _after_ f and g.
         */
//        memmove(t1, Fd, 2 * hn * dlen * sizeof *Fd);
        tmp = hn * dlen;
        System.arraycopy(srctmp, Fd, srctmp, t1, tmp + tmp);
        Fd = t1;
        Gd = Fd + tmp;

        /*
         * 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++)
        {
            int p, p0i, R2, Rx;
            int v;
            int xs, ys, xd, yd;

            p = FalconSmallPrimeList.PRIMES[u].p;
            p0i = modp_ninv31(p);
            R2 = modp_R2(p, p0i);
            Rx = modp_Rx(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)
            {
                srctmp[xd] = zint_mod_small_signed(srctmp, xs, dlen, p, p0i, R2, Rx);
                srctmp[yd] = zint_mod_small_signed(srctmp, 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++)
        {
            int p, p0i, R2;
            int gm, igm, fx, gx, Fp, Gp;
            int v;

            /*
             * All computations are done modulo p.
             */
            p = FalconSmallPrimeList.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(srctmp, ft, slen, slen, n, 1, srctmp, t1);
                zint_rebuild_CRT(srctmp, gt, slen, slen, n, 1, srctmp, t1);
            }

            gm = t1;
            igm = gm + n;
            fx = igm + n;
            gx = fx + n;

            modp_mkgm2(srctmp, gm, srctmp, igm, logn, FalconSmallPrimeList.PRIMES[u].g, p, p0i);

            if (u < slen)
            {
                for (v = 0, x = ft + u, y = gt + u;
                     v < n; v++, x += slen, y += slen)
                {
                    srctmp[fx + v] = srctmp[x];
                    srctmp[gx + v] = srctmp[y];
                }
                modp_iNTT2_ext(srctmp, ft + u, slen, srctmp, igm, logn, p, p0i);
                modp_iNTT2_ext(srctmp, gt + u, slen, srctmp, igm, logn, p, p0i);
            }
            else
            {
                int Rx;

                Rx = modp_Rx(slen, p, p0i, R2);
                for (v = 0, x = ft, y = gt;
                     v < n; v++, x += slen, y += slen)
                {
                    srctmp[fx + v] = zint_mod_small_signed(srctmp, x, slen,
                        p, p0i, R2, Rx);
                    srctmp[gx + v] = zint_mod_small_signed(srctmp, y, slen,
                        p, p0i, R2, Rx);
                }
                modp_NTT2(srctmp, fx, srctmp, gm, logn, p, p0i);
                modp_NTT2(srctmp, gx, srctmp, 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)
            {
                srctmp[Fp + v] = srctmp[x];
                srctmp[Gp + v] = srctmp[y];
            }
            modp_NTT2(srctmp, Fp, srctmp, gm, logn - 1, p, p0i);
            modp_NTT2(srctmp, Gp, srctmp, 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))
            {
                int ftA, ftB, gtA, gtB;
                int mFp, mGp;

                ftA = srctmp[fx + (v << 1)];
                ftB = srctmp[fx + (v << 1) + 1];
                gtA = srctmp[gx + (v << 1)];
                gtB = srctmp[gx + (v << 1) + 1];
                mFp = modp_montymul(srctmp[Fp + v], R2, p, p0i);
                mGp = modp_montymul(srctmp[Gp + v], R2, p, p0i);
                srctmp[x] = modp_montymul(gtB, mFp, p, p0i);
                srctmp[x + llen] = modp_montymul(gtA, mFp, p, p0i);
                srctmp[y] = modp_montymul(ftB, mGp, p, p0i);
                srctmp[y + llen] = modp_montymul(ftA, mGp, p, p0i);
            }
            modp_iNTT2_ext(srctmp, Ft + u, llen, srctmp, igm, logn, p, p0i);
            modp_iNTT2_ext(srctmp, Gt + u, llen, srctmp, igm, logn, p, p0i);
        }

        /*
         * Rebuild F and G with the CRT.
         */
        zint_rebuild_CRT(srctmp, Ft, llen, llen, n, 1, srctmp, t1);
        zint_rebuild_CRT(srctmp, Gt, llen, llen, n, 1, srctmp, 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).
         */
        // Carve the 5 FFT regions out of the persistent fprBuf scratch.
        // rt1..rt4 each occupy n doubles; rt5 occupies n/2 doubles.
        rt1 = 0;
        rt2 = rt1 + n;
        rt3 = rt2 + n;
        rt4 = rt3 + n;
        rt5 = rt4 + n;
        k = kBuf;

        /*
         * 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 = Math.min(slen, 10);
        poly_big_to_fp(fprBuf, rt3, srctmp, ft + slen - rlen, rlen, slen, logn);
        poly_big_to_fp(fprBuf, rt4, srctmp, gt + slen - rlen, rlen, slen, logn);

        /*
         * Values in rt3 and rt4 are downscaled by 2^(scale_fg).
         */
        scale_fg = 31 * (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_avg[depth] - 6 * bitlength_std[depth];
        maxbl_fg = bitlength_avg[depth] + 6 * bitlength_std[depth];

        /*
         * Compute 1/(f*adj(f)+g*adj(g)) in rt5. We also keep adj(f)
         * and adj(g) in rt3 and rt4, respectively.
         */
        FalconFFT.FFT(fprBuf, rt3, logn);
        FalconFFT.FFT(fprBuf, rt4, logn);
        FalconFFT.poly_invnorm2_fft(fprBuf, rt5, fprBuf, rt3, fprBuf, rt4, logn);
        FalconFFT.poly_adj_fft(fprBuf, rt3, logn);
        FalconFFT.poly_adj_fft(fprBuf, 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 * 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;
            int scl, sch;
            double pdc, pt;

            /*
             * Convert current F and G into floating-point. We apply
             * scaling if the current length is more than 10 words.
             */
            rlen = Math.min(FGlen, 10);
            scale_FG = 31 * (FGlen - rlen);
            poly_big_to_fp(fprBuf, rt1, srctmp, Ft + FGlen - rlen, rlen, llen, logn);
            poly_big_to_fp(fprBuf, rt2, srctmp, Gt + FGlen - rlen, rlen, llen, logn);

            /*
             * Compute (F*adj(f)+G*adj(g))/(f*adj(f)+g*adj(g)) in rt2.
             */
            FalconFFT.FFT(fprBuf, rt1, logn);
            FalconFFT.FFT(fprBuf, rt2, logn);
            FalconFFT.poly_mul_fft(fprBuf, rt1, fprBuf, rt3, logn);
            FalconFFT.poly_mul_fft(fprBuf, rt2, fprBuf, rt4, logn);
            FalconFFT.poly_add(fprBuf, rt2, fprBuf, rt1, logn);
            FalconFFT.poly_mul_autoadj_fft(fprBuf, rt2, fprBuf, rt5, logn);
            FalconFFT.iFFT(fprBuf, 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 = FPREngine.fpr_two;
            }
            else
            {
                pt = FPREngine.fpr_onehalf;
            }
            pdc = FPREngine.fpr_one;
            while (dc != 0)
            {
                if ((dc & 1) != 0)
                {
                    pdc = pdc * pt;
                }
                dc >>= 1;
                pt *= pt;
            }

            for (u = 0; u < n; u++)
            {
                double xv = fprBuf[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 (FPREngine.fpr_mtwo31m1 >= xv
                    || xv >= FPREngine.fpr_ptwo31m1)
                {
                    return 0;
                }
                k[u] = (int)FPREngine.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 = (scale_k / 31);
            scl = (scale_k % 31);
            if (depth <= DEPTH_INT_FG)
            {
                poly_sub_scaled_ntt(srctmp, Ft, FGlen, llen, srctmp, ft, slen, slen,
                    k, sch, scl, logn, srctmp, t1);
                poly_sub_scaled_ntt(srctmp, Gt, FGlen, llen, srctmp, gt, slen, slen,
                    k, sch, scl, logn, srctmp, t1);
            }
            else
            {
                poly_sub_scaled(srctmp, Ft, FGlen, llen, srctmp, ft, slen, slen,
                    k, sch, scl, logn);
                poly_sub_scaled(srctmp, Gt, FGlen, llen, srctmp, 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 (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)
            {
                int v;
                int sw;

                sw = -(srctmp[Ft + FGlen - 1] >>> 30) >>> 1;
                for (v = FGlen; v < slen; v++)
                {
                    srctmp[Ft + v] = sw;
                }
                sw = -(srctmp[Gt + FGlen - 1] >>> 30) >>> 1;
                for (v = FGlen; v < slen; v++)
                {
                    srctmp[Gt + v] = sw;
                }
            }
        }

        /*
         * Compress encoding of all values to 'slen' words (this is the
         * expected output format).
         */
        for (u = 0, x = 0, y = 0;
             u < (n << 1); u++, x += slen, y += llen)
        {
//            memmove(x, y, slen * sizeof *y);
            System.arraycopy(srctmp, y, srctmp, x, slen);
        }
        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.
     */
    private static int solve_NTRU_binary_depth1(int logn_top,
                                                byte[] srcf, byte[] srcg, int[] srctmp,
                                                double[] fprBuf)
    {
        /*
         * 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.
         */
        int depth, logn;
        int n_top, n, hn, slen, dlen, llen, u;
        int Fd, Gd, Ft, Gt, ft, gt, t1;
        // rt1..rt6 are offsets into the caller-supplied fprBuf scratch.
        // rt1..rt5 are size n; rt6 is size n/2.
        int rt1, rt2, rt3, rt4, rt5, rt6;
        int x, y;

        depth = 1;
        n_top = 1 << logn_top;
        logn = logn_top - depth;
        n = 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 = 0;
        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++)
        {
            int p, p0i, R2, Rx;
            int v;
            int xs, ys, xd, yd;

            p = FalconSmallPrimeList.PRIMES[u].p;
            p0i = modp_ninv31(p);
            R2 = modp_R2(p, p0i);
            Rx = modp_Rx(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)
            {
                srctmp[xd] = zint_mod_small_signed(srctmp, xs, dlen, p, p0i, R2, Rx);
                srctmp[yd] = zint_mod_small_signed(srctmp, 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));
        System.arraycopy(srctmp, Ft, srctmp, 0, llen * n);
        Ft = 0;
//        memmove(Ft + llen * n, Gt, llen * n * sizeof(uint32_t));
        System.arraycopy(srctmp, Gt, srctmp, Ft + llen * n, llen * n);
        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++)
        {
            int p, p0i, R2;
            int gm, igm, fx, gx, Fp, Gp;
            int e;
            int v;

            /*
             * All computations are done modulo p.
             */
            p = FalconSmallPrimeList.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(srctmp, gm, srctmp, igm, logn_top, FalconSmallPrimeList.PRIMES[u].g, p, p0i);

            /*
             * Set ft and gt to f and g modulo p, respectively.
             */
            for (v = 0; v < n_top; v++)
            {
                srctmp[fx + v] = modp_set(srcf[v], p);
                srctmp[gx + v] = modp_set(srcg[v], p);
            }

            /*
             * Convert to NTT and compute our f and g.
             */
            modp_NTT2(srctmp, fx, srctmp, gm, logn_top, p, p0i);
            modp_NTT2(srctmp, gx, srctmp, gm, logn_top, p, p0i);
            for (e = logn_top; e > logn; e--)
            {
                modp_poly_rec_res(srctmp, fx, e, p, p0i, R2);
                modp_poly_rec_res(srctmp, 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);
            System.arraycopy(srctmp, igm, srctmp, gm + n, n);
            igm = gm + n;
//                memmove(igm + n, fx, n * sizeof *ft);
            System.arraycopy(srctmp, fx, srctmp, igm + n, n);
            fx = igm + n;
//                memmove(fx + n, gx, n * sizeof *gt);
            System.arraycopy(srctmp, gx, srctmp, fx + n, n);
            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)
            {
                srctmp[Fp + v] = srctmp[x];
                srctmp[Gp + v] = srctmp[y];
            }
            modp_NTT2(srctmp, Fp, srctmp, gm, logn - 1, p, p0i);
            modp_NTT2(srctmp, Gp, srctmp, 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))
            {
                int ftA, ftB, gtA, gtB;
                int mFp, mGp;

                ftA = srctmp[fx + (v << 1)];
                ftB = srctmp[fx + (v << 1) + 1];
                gtA = srctmp[gx + (v << 1)];
                gtB = srctmp[gx + (v << 1) + 1];
                mFp = modp_montymul(srctmp[Fp + v], R2, p, p0i);
                mGp = modp_montymul(srctmp[Gp + v], R2, p, p0i);
                srctmp[x] = modp_montymul(gtB, mFp, p, p0i);
                srctmp[x + llen] = modp_montymul(gtA, mFp, p, p0i);
                srctmp[y] = modp_montymul(ftB, mGp, p, p0i);
                srctmp[y + llen] = modp_montymul(ftA, mGp, p, p0i);
            }
            modp_iNTT2_ext(srctmp, Ft + u, llen, srctmp, igm, logn, p, p0i);
            modp_iNTT2_ext(srctmp, Gt + u, llen, srctmp, igm, logn, p, p0i);

            /*
             * Also save ft and gt (only up to size slen).
             */
            if (u < slen)
            {
                modp_iNTT2(srctmp, fx, srctmp, igm, logn, p, p0i);
                modp_iNTT2(srctmp, gx, srctmp, igm, logn, p, p0i);
                for (v = 0, x = ft + u, y = gt + u;
                     v < n; v++, x += slen, y += slen)
                {
                    srctmp[x] = srctmp[fx + v];
                    srctmp[y] = srctmp[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(srctmp, Ft, llen, llen, n << 1, 1, srctmp, t1);
        zint_rebuild_CRT(srctmp, ft, slen, slen, n << 1, 1, srctmp, 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.
         */

        /*
         * Carve rt1..rt6 from the persistent fprBuf scratch.
         * rt1..rt5 each size n; rt6 size n/2.
         */
        rt1 = 0;
        rt2 = rt1 + n;
        rt3 = rt2 + n;
        rt4 = rt3 + n;
        rt5 = rt4 + n;
        rt6 = rt5 + n;

        /*
         * Convert F and G into floating point (rt1 and rt2).
         */
        poly_big_to_fp(fprBuf, rt1, srctmp, Ft, llen, llen, logn);
        poly_big_to_fp(fprBuf, rt2, srctmp, 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);
        System.arraycopy(srctmp, ft, srctmp, 0, 2 * slen * n);
        ft = 0;
        gt = ft + slen * n;

        /*
         * Convert f and g into floating point (rt3 and rt4).
         */
        poly_big_to_fp(fprBuf, rt3, srctmp, ft, slen, slen, logn);
        poly_big_to_fp(fprBuf, rt4, srctmp, gt, slen, slen, logn);

        /*
         * We now have:
         *   rt1 = F
         *   rt2 = G
         *   rt3 = f
         *   rt4 = g
         * in that order in RAM. We convert all of them to FFT.
         */
        FalconFFT.FFT(fprBuf, rt1, logn);
        FalconFFT.FFT(fprBuf, rt2, logn);
        FalconFFT.FFT(fprBuf, rt3, logn);
        FalconFFT.FFT(fprBuf, rt4, logn);

        /*
         * Compute:
         *   rt5 = F*adj(f) + G*adj(g)
         *   rt6 = 1 / (f*adj(f) + g*adj(g))
         * (Note that rt6 is half-length.)
         */
        FalconFFT.poly_add_muladj_fft(fprBuf, rt5, fprBuf, rt1, fprBuf, rt2, fprBuf, rt3, fprBuf, rt4, logn);
        FalconFFT.poly_invnorm2_fft(fprBuf, rt6, fprBuf, rt3, fprBuf, rt4, logn);

        /*
         * Compute:
         *   rt5 = (F*adj(f)+G*adj(g)) / (f*adj(f)+g*adj(g))
         */
        FalconFFT.poly_mul_autoadj_fft(fprBuf, rt5, fprBuf, 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.
         */
        FalconFFT.iFFT(fprBuf, rt5, logn);
        for (u = 0; u < n; u++)
        {
            double z;

            z = fprBuf[rt5 + u];
            if (z >= FPREngine.fpr_ptwo63m1 || FPREngine.fpr_mtwo63m1 >= z)
            {
                return 0;
            }
            fprBuf[rt5 + u] = FPREngine.fpr_rint(z);
        }
        FalconFFT.FFT(fprBuf, rt5, logn);

        /*
         * Subtract k*f from F, and k*g from G.
         */
        FalconFFT.poly_mul_fft(fprBuf, rt3, fprBuf, rt5, logn);
        FalconFFT.poly_mul_fft(fprBuf, rt4, fprBuf, rt5, logn);
        FalconFFT.poly_sub(fprBuf, rt1, fprBuf, rt3, logn);
        FalconFFT.poly_sub(fprBuf, rt2, fprBuf, rt4, logn);
        FalconFFT.iFFT(fprBuf, rt1, logn);
        FalconFFT.iFFT(fprBuf, rt2, logn);

        /*
         * Convert back F and G to integers, and return.
         */
        Gt = Ft + n;
        for (u = 0; u < n; u++)
        {
            srctmp[Ft + u] = (int)FPREngine.fpr_rint(fprBuf[rt1 + u]);
            srctmp[Gt + u] = (int)FPREngine.fpr_rint(fprBuf[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.
     */
    private static int solve_NTRU_binary_depth0(int logn,
                                                byte[] srcf, byte[] srcg, int[] srctmp,
                                                double[] fprBuf)
    {
        int n, hn, u;
        int p, p0i, R2;
        int Fp, Gp, t1, t2, t3, t4, t5;
        int gm, igm, ft, gt;
        int rt1, rt2, rt3;

        n = 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 = FalconSmallPrimeList.PRIMES[0].p;
        p0i = modp_ninv31(p);
        R2 = modp_R2(p, p0i);

        Fp = 0;
        Gp = Fp + hn;
        ft = Gp + hn;
        gt = ft + n;
        gm = gt + n;
        igm = gm + n;

        modp_mkgm2(srctmp, gm, srctmp, igm, logn, FalconSmallPrimeList.PRIMES[0].g, p, p0i);

        /*
         * Convert F' anf G' in NTT representation.
         */
        for (u = 0; u < hn; u++)
        {
            srctmp[Fp + u] = modp_set(zint_one_to_plain(srctmp, Fp + u), p);
            srctmp[Gp + u] = modp_set(zint_one_to_plain(srctmp, Gp + u), p);
        }
        modp_NTT2(srctmp, Fp, srctmp, gm, logn - 1, p, p0i);
        modp_NTT2(srctmp, Gp, srctmp, gm, logn - 1, p, p0i);

        /*
         * Load f and g and convert them to NTT representation.
         */
        for (u = 0; u < n; u++)
        {
            srctmp[ft + u] = modp_set(srcf[u], p);
            srctmp[gt + u] = modp_set(srcg[u], p);
        }
        modp_NTT2(srctmp, ft, srctmp, gm, logn, p, p0i);
        modp_NTT2(srctmp, gt, srctmp, gm, logn, p, p0i);

        /*
         * Build the unreduced F,G in ft and gt.
         */
        for (u = 0; u < n; u += 2)
        {
            int ftA, ftB, gtA, gtB;
            int mFp, mGp;

            ftA = srctmp[ft + u];
            ftB = srctmp[ft + u + 1];
            gtA = srctmp[gt + u];
            gtB = srctmp[gt + u + 1];
            mFp = modp_montymul(srctmp[Fp + (u >> 1)], R2, p, p0i);
            mGp = modp_montymul(srctmp[Gp + (u >> 1)], R2, p, p0i);
            srctmp[ft + u] = modp_montymul(gtB, mFp, p, p0i);
            srctmp[ft + u + 1] = modp_montymul(gtA, mFp, p, p0i);
            srctmp[gt + u] = modp_montymul(ftB, mGp, p, p0i);
            srctmp[gt + u + 1] = modp_montymul(ftA, mGp, p, p0i);
        }
        modp_iNTT2(srctmp, ft, srctmp, igm, logn, p, p0i);
        modp_iNTT2(srctmp, gt, srctmp, igm, logn, p, p0i);

        Gp = Fp + n;
        t1 = Gp + n;
//        memmove(Fp, ft, 2 * n * sizeof *ft);
        System.arraycopy(srctmp, ft, srctmp, Fp, 2 * n);

        /*
         * 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(srctmp, t1, srctmp, t2, logn, FalconSmallPrimeList.PRIMES[0].g, p, p0i);

        /*
         * Convert F and G to NTT.
         */
        modp_NTT2(srctmp, Fp, srctmp, t1, logn, p, p0i);
        modp_NTT2(srctmp, Gp, srctmp, t1, logn, p, p0i);

        /*
         * Load f and adj(f) in t4 and t5, and convert them to NTT
         * representation.
         */
        srctmp[t4] = srctmp[t5] = modp_set(srcf[0], p);
        for (u = 1; u < n; u++)
        {
            srctmp[t4 + u] = modp_set(srcf[u], p);
            srctmp[t5 + n - u] = modp_set(-srcf[u], p);
        }
        modp_NTT2(srctmp, t4, srctmp, t1, logn, p, p0i);
        modp_NTT2(srctmp, t5, srctmp, t1, logn, p, p0i);

        /*
         * Compute F*adj(f) in t2, and f*adj(f) in t3.
         */
        for (u = 0; u < n; u++)
        {
            int w;

            w = modp_montymul(srctmp[t5 + u], R2, p, p0i);
            srctmp[t2 + u] = modp_montymul(w, srctmp[Fp + u], p, p0i);
            srctmp[t3 + u] = modp_montymul(w, srctmp[t4 + u], p, p0i);
        }

        /*
         * Load g and adj(g) in t4 and t5, and convert them to NTT
         * representation.
         */
        srctmp[t4] = srctmp[t5] = modp_set(srcg[0], p);
        for (u = 1; u < n; u++)
        {
            srctmp[t4 + u] = modp_set(srcg[u], p);
            srctmp[t5 + n - u] = modp_set(-srcg[u], p);
        }
        modp_NTT2(srctmp, t4, srctmp, t1, logn, p, p0i);
        modp_NTT2(srctmp, t5, srctmp, t1, logn, p, p0i);

        /*
         * Add G*adj(g) to t2, and g*adj(g) to t3.
         */
        for (u = 0; u < n; u++)
        {
            int w;

            w = modp_montymul(srctmp[t5 + u], R2, p, p0i);
            srctmp[t2 + u] = modp_add(srctmp[t2 + u],
                modp_montymul(w, srctmp[Gp + u], p, p0i), p);
            srctmp[t3 + u] = modp_add(srctmp[t3 + u],
                modp_montymul(w, srctmp[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(srctmp, t1, srctmp, t4, logn, FalconSmallPrimeList.PRIMES[0].g, p, p0i);
        modp_iNTT2(srctmp, t2, srctmp, t4, logn, p, p0i);
        modp_iNTT2(srctmp, t3, srctmp, t4, logn, p, p0i); // TODO fix binary_depth0 -> t1 value is wrong
        for (u = 0; u < n; u++)
        {
            srctmp[t1 + u] = modp_norm(srctmp[t2 + u], p);
            srctmp[t2 + u] = modp_norm(srctmp[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.
         */
        // rt1..rt3 are offsets into the caller-supplied fprBuf scratch.
        // fprBuf is sized in keygen() to at least 3*n (in fact 5*n) doubles.
        rt1 = 0;
        rt2 = rt1 + n;
        rt3 = rt2 + n;
        for (u = 0; u < n; u++)
        {
            fprBuf[rt3 + u] = srctmp[t2 + u];
        }
        FalconFFT.FFT(fprBuf, rt3, logn);
        System.arraycopy(fprBuf, rt3, fprBuf, rt2, hn);

        /*
         * Convert F*adj(f)+G*adj(g) in FFT representation.
         */
        rt3 = rt2 + hn;
        for (u = 0; u < n; u++)
        {
            fprBuf[rt3 + u] = srctmp[t1 + u];
        }
        FalconFFT.FFT(fprBuf, rt3, logn);

        /*
         * Compute (F*adj(f)+G*adj(g))/(f*adj(f)+g*adj(g)) and get
         * its rounded normal representation in t1.
         */
        FalconFFT.poly_div_autoadj_fft(fprBuf, rt3, fprBuf, rt2, logn);
        FalconFFT.iFFT(fprBuf, rt3, logn);
        for (u = 0; u < n; u++)
        {
            srctmp[t1 + u] = modp_set((int)FPREngine.fpr_rint(fprBuf[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(srctmp, t2, srctmp, t3, logn, FalconSmallPrimeList.PRIMES[0].g, p, p0i);
        for (u = 0; u < n; u++)
        {
            srctmp[t4 + u] = modp_set(srcf[u], p);
            srctmp[t5 + u] = modp_set(srcg[u], p);
        }
        modp_NTT2(srctmp, t1, srctmp, t2, logn, p, p0i);
        modp_NTT2(srctmp, t4, srctmp, t2, logn, p, p0i);
        modp_NTT2(srctmp, t5, srctmp, t2, logn, p, p0i);
        for (u = 0; u < n; u++)
        {
            int kw;

            kw = modp_montymul(srctmp[t1 + u], R2, p, p0i);
            srctmp[Fp + u] = modp_sub(srctmp[Fp + u],
                modp_montymul(kw, srctmp[t4 + u], p, p0i), p);
            srctmp[Gp + u] = modp_sub(srctmp[Gp + u],
                modp_montymul(kw, srctmp[t5 + u], p, p0i), p);
        }
        modp_iNTT2(srctmp, Fp, srctmp, t3, logn, p, p0i);
        modp_iNTT2(srctmp, Gp, srctmp, t3, logn, p, p0i);
        for (u = 0; u < n; u++)
        {
            srctmp[Fp + u] = modp_norm(srctmp[Fp + u], p);
            srctmp[Gp + u] = modp_norm(srctmp[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.
     */
    private static int solve_NTRU(int logn, byte[] srcF, //byte[] srcG,
                                  byte[] srcf, byte[] srcg, int lim, int[] srctmp,
                                  double[] fprBuf, int[] kBuf)
    {
        int n, u;
        int ft, gt, Ft, Gt, gm;
        int p, p0i, r;
        int G = 0;
        n = mkn(logn);

        if (solve_NTRU_deepest(logn, srcf, srcg, srctmp) == 0)
        {
            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)
        {
            int depth;

            depth = logn;
            while (depth-- > 0)
            {
                if (solve_NTRU_intermediate(logn, srcf, srcg, depth, srctmp, fprBuf, kBuf) == 0)
                {
                    return 0;
                }
            }
        }
        else
        {
            int depth = logn;
            while (depth-- > 2)
            {
                if (solve_NTRU_intermediate(logn, srcf, srcg, depth, srctmp, fprBuf, kBuf) == 0)
                {
                    return 0;
                }
            }
            if (solve_NTRU_binary_depth1(logn, srcf, srcg, srctmp, fprBuf) == 0)
            {
                return 0;
            }
            if (solve_NTRU_binary_depth0(logn, srcf, srcg, srctmp, fprBuf) == 0)
            {
                return 0;
            }
        }

        /*
         * If no buffer has been provided for G, use a temporary one.
         */
//        if (srcG == null)
//        {
//            G = 0;
//            srcG = new byte[n];
//        }
        byte[] srcG = new byte[n];

        /*
         * Final F and G are in fk->tmp, one word per coefficient
         * (signed value over 31 bits).
         */
        if (poly_big_to_small(srcF, 0, srctmp, 0, lim, logn) == 0
            || poly_big_to_small(srcG, G, srctmp, n, lim, logn) == 0)
        {
            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 = 0;
        ft = Gt + n;
        gt = ft + n;
        Ft = gt + n;
        gm = Ft + n;

        p = FalconSmallPrimeList.PRIMES[0].p;
        p0i = modp_ninv31(p);
        modp_mkgm2(srctmp, gm, srctmp, 0, logn, FalconSmallPrimeList.PRIMES[0].g, p, p0i);
        for (u = 0; u < n; u++)
        {
            srctmp[Gt + u] = modp_set(srcG[G + u], p);
        }
        for (u = 0; u < n; u++)
        {
            srctmp[ft + u] = modp_set(srcf[u], p);
            srctmp[gt + u] = modp_set(srcg[u], p);
            srctmp[Ft + u] = modp_set(srcF[u], p);
        }
        modp_NTT2(srctmp, ft, srctmp, gm, logn, p, p0i);
        modp_NTT2(srctmp, gt, srctmp, gm, logn, p, p0i);
        modp_NTT2(srctmp, Ft, srctmp, gm, logn, p, p0i);
        modp_NTT2(srctmp, Gt, srctmp, gm, logn, p, p0i);
        r = modp_montymul(12289, 1, p, p0i);
        for (u = 0; u < n; u++)
        {
            int z;

            z = modp_sub(modp_montymul(srctmp[ft + u], srctmp[Gt + u], p, p0i),
                modp_montymul(srctmp[gt + u], srctmp[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.
     */
    private static void poly_small_mkgauss(SHAKEDigest rng, byte[] srcf, int logn, byte[] rngBuf8)
    {
        int n, u;
        int mod2;

        n = mkn(logn);
        mod2 = 0;
        for (u = 0; u < n; u++)
        {
            int s;

            for (; ; )
            {
                s = mkgauss(rng, logn, rngBuf8);

                /*
                 * 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)
                {
                    continue;
                }

                /*
                 * 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 ^ (s & 1)) == 0)
                    {
                        continue;
                    }
                }
                else
                {
                    mod2 ^= (s & 1);
                }
                srcf[u] = (byte)s;
                break;
            }
        }
    }

    /* see falcon.h */
    static void keygen(SHAKEDigest rc,
                       byte[] srcf, byte[] srcg, byte[] srcF, short[] srch,
                       int logn)
    {
        /*
         * 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.
         */
        int n, u;
        int[] itmp;
        //byte[] btmp;
        short[] stmp;
        double[] ftmp;
        int h2, tmp2;
        //SHAKE256 rc;

        n = mkn(logn);
        //rc = rng;

        /*
         * Allocate the per-keygen scratch buffers once, outside the
         * retry loop. Falcon expects ~4 retries on average per keygen
         * (norm check, b-norm check, NTRU-solve failure can each force
         * a restart), so a per-retry alloc was 3-4x redundant. None of
         * these buffers carry state across retries; every retry
         * overwrites them from scratch.
         */
        ftmp = new double[3 * n];
        stmp = new short[2 * n];
        itmp = logn > 2 ? new int[28 * n] : new int[28 * n * 3];
        byte[] rngBuf8 = new byte[8];
        /*
         * Persistent FFT/k scratch threaded through solve_NTRU* so the
         * recursive solver does not re-allocate rt1..rt5/rt6/k per call.
         * Sizing rules (max-by-callee, with n meaning the local degree):
         *   - solve_NTRU_intermediate: 5 regions of size n_depth, used
         *     as rt1..rt4 (size n) + rt5 (size n/2) + k stored in kBuf;
         *     needs 4.5 * n_depth doubles. n_depth <= top_n/4 for the
         *     realistic logn >= 3 path, <= top_n for logn <= 2.
         *   - solve_NTRU_binary_depth1: rt1..rt5 (size n=top_n/2) + rt6
         *     (size n/2); needs 5.5 * top_n/2 = 2.75 * top_n.
         *   - solve_NTRU_binary_depth0: rt1, rt2, rt3 in a single
         *     double[3*n] (n = top_n); needs 3 * top_n.
         * Pick 5*n to cover all of the above with headroom.
         */
        double[] fprBuf = new double[5 * n];
        int[] kBuf = new int[n];

        /*
         * 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 (; ; )
        {
            int rt1, rt2, rt3;
            double bnorm;
            int 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, srcf, logn, rngBuf8);
            poly_small_mkgauss(rc, srcg, logn, rngBuf8);

            /*
             * 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 << (FalconCodec.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 (srcf[u] >= lim || srcf[u] <= -lim
                    || srcg[u] >= lim || srcg[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(srcf, logn);
            normg = poly_small_sqnorm(srcg, logn);
            norm = (normf + normg) | -((normf | normg) >>> 31);
            if ((norm & 0xffffffffL) >= 16823L)
            {
                continue;
            }

            /*
             * We compute the orthogonalized vector norm.
             */
            rt1 = 0;
            rt2 = rt1 + n;
            rt3 = rt2 + n;
            poly_small_to_fp(ftmp, rt1, srcf, logn);
            poly_small_to_fp(ftmp, rt2, srcg, logn);
            FalconFFT.FFT(ftmp, rt1, logn);
            FalconFFT.FFT(ftmp, rt2, logn);
            FalconFFT.poly_invnorm2_fft(ftmp, rt3, ftmp, rt1, ftmp, rt2, logn);
            FalconFFT.poly_adj_fft(ftmp, rt1, logn);
            FalconFFT.poly_adj_fft(ftmp, rt2, logn);
            FalconFFT.poly_mulconst(ftmp, rt1, FPREngine.fpr_q, logn);
            FalconFFT.poly_mulconst(ftmp, rt2, FPREngine.fpr_q, logn);
            FalconFFT.poly_mul_autoadj_fft(ftmp, rt1, ftmp, rt3, logn);
            FalconFFT.poly_mul_autoadj_fft(ftmp, rt2, ftmp, rt3, logn);
            FalconFFT.iFFT(ftmp, rt1, logn);
            FalconFFT.iFFT(ftmp, rt2, logn);
            bnorm = FPREngine.fpr_zero;
            for (u = 0; u < n; u++)
            {
                bnorm += ftmp[rt1 + u] * ftmp[rt1 + u] + ftmp[rt2 + u] * ftmp[rt2 + u];
            }
            //if (!FPREngine.fpr_lt(bnorm, FPREngine.fpr_bnorm_max))
            if (bnorm >= FPREngine.fpr_bnorm_max)
            {
                continue;
            }

            /*
             * Compute public key h = g/f mod X^N+1 mod q. If this
             * fails, we must restart.
             */
            if (srch == null)
            {
                h2 = 0;
                srch = stmp;
                tmp2 = h2 + n;
            }
            else
            {
                h2 = 0;
                tmp2 = 0;
            }
            if (FalconVrfy.compute_public(srch, h2, srcf, srcg, logn, stmp, tmp2) == 0)
            {
                continue;
            }

            /*
             * Solve the NTRU equation to get F and G.
             */
            lim = (1 << (FalconCodec.max_FG_bits[logn] - 1)) - 1;
            if (solve_NTRU(logn, srcF, srcf, srcg, lim, itmp, fprBuf, kBuf) == 0)
            {
                continue;
            }

            /*
             * Key pair is generated.
             */
            break;
        }
    }

    private static long toUnsignedLong(int x)
    {
        return x & 0xffffffffL;
    }
}
