ntt64.c

//  ntt64.c
//  Copyright (c) 2023 Raccoon Signature Team. See LICENSE.

//  === 64-bit Number Theoretic Transform

#ifndef POLYR_Q32

#include <stddef.h>
#include <stdbool.h>

#include "polyr.h"
#include "mont64.h"

//  === Roots of unity constants

/*
    (g is the smallest full generator for both q1 and q2)

    n   = 512
    q   = 549824583172097
    r   = 2^64 % q
    g   = Mod(15, q)
    h   = g^(znorder(g)/(2*n))

    bitrev(n,i) = \
      b = binary(n + (i%n));\
      sum(i = 2, length(b), 2^(i-2) * b[i])

    Generating the Montgomery ("r") scaled constants racc_w_64[511]:

    w = vector(511,i,lift(r*h^bitrev(n,i)))

    (   the function polyr_fntt() evaluates the input, intepreted
        as a polynomial with the 0-degree coefficient first,
        at the roots of unity h^(2*bitrev(n,i)+1) for i= 0, 1, 2,.. )

    e = vector(512,i,lift(h^(2 *bitrev(n,i-1)+1)));
    evalv(v,x) = sum(i=0,511,x^i*v[i+1])
    fntt(v) = vector(512,i,lift(evalv(v, Mod(e[i],q))))

*/

static const int64_t racc_w_64[511] = {
    209408931609897, 131076814237138, 46162320619533,  192435386452264,
    45147673384707,  199893044697986, 73227418396945,  37721780595870,
    417006390137132, 530441435114090, 293293967316741, 77537632964119,
    397683601422803, 355158546321898, 97000472436250,  9633177676322,
    2540248838235,   133027841664207, 542469485928921, 398979303753866,
    324032137440331, 461282595569708, 249216434398905, 373384824028144,
    171404720899221, 125987723682393, 449035531792827, 248577544161206,
    263822611259530, 348988546829164, 102450415821103, 489329478724880,
    83289879148447,  80268611727648,  290625880469644, 47471514278669,
    71942871131238,  415089145165827, 9564037332679,   413384901302967,
    138722501121331, 24875604016876,  333873791079803, 243806517656296,
    31220627722436,  43247531123334,  341022237839572, 96550402747764,
    199496477777928, 156648627347877, 138092088059953, 46217145081118,
    540932350368814, 6879053396002,   262308844368444, 457468571783904,
    178532388281334, 57709566308463,  338919316083985, 45258448649820,
    432375600852479, 349893754650677, 123781869489683, 252323176356550,
    434426191548097, 167323232543988, 86827205937764,  199107214819245,
    70706228320110,  536721620900782, 52182436153238,  441064991745544,
    548454424775682, 405285370495304, 245424562571178, 527428102845750,
    472240708273764, 31949557679847,  323267515331585, 432916134513901,
    291358237599102, 136523991411120, 413460116276708, 282297447295031,
    256650373660298, 335928320528141, 385799155722515, 224008605394786,
    157070826751996, 262845589447882, 222518606426787, 472936778670735,
    167462952897800, 292721318786275, 385665278640984, 316375138619932,
    109804615019434, 121813705958567, 234390170722256, 64753722389307,
    226177994089190, 493983720478949, 202598179185893, 280748684684490,
    120557505069973, 1279402686367,   403132393364157, 59443425320688,
    267570283561684, 464436526155180, 33995302249681,  477102774398853,
    379159891585345, 150320113205964, 10766296740086,  498302386366140,
    140715050694842, 491229317852731, 130015534262650, 431859439565299,
    76566911957484,  299902094060553, 338039289793685, 246617738944365,
    394781847683933, 435691047042537, 280197530993299, 453150527671934,
    448509752488338, 203343744250240, 480783394616494, 127427685221182,
    3226137114466,   437580988001252, 175912055457621, 485397369594776,
    132259143016852, 496972583725889, 303726932107212, 169964826931121,
    364182014812580, 106010716015933, 400949249209847, 489213600710291,
    327920368362372, 317118730496425, 525646850710563, 505281666416937,
    549230012764282, 162881645841977, 421040665308142, 152586264481105,
    55240039421504,  242380069509750, 256949899110475, 28900446367198,
    400737996131246, 441816438107175, 65990781490661,  411905908218796,
    361666241755396, 413080517160721, 61872464009638,  206077790377089,
    35691738614920,  513523432625931, 543778954751551, 463135208359438,
    175146450709671, 88814037929636,  491016955282,    169763039968971,
    32069143859017,  419747025907333, 524456666745855, 481637469887022,
    283485578235052, 84293750227657,  84739787618337,  252499783197300,
    45490563370658,  251116276790200, 127494951710650, 221718508081815,
    377257286321911, 426957694002132, 273483748114935, 230515989870991,
    155839402922807, 42173322106903,  530779529119950, 252254089305824,
    437439347249343, 374715973207311, 377438202807575, 454635050866964,
    445690073958189, 482836966305777, 501474227211310, 51779491964020,
    161427998519961, 294675711448947, 513247978393719, 222573935667997,
    279016404493033, 465305230276051, 467077293582288, 241030454841721,
    300799038375642, 395426361885876, 54852749495588,  120213016122640,
    323399198250458, 431940263453913, 11655351490709,  251217189677313,
    204706088503911, 147690595465162, 283271661992647, 449442014274831,
    60668939127614,  396216382941742, 325654597720256, 429691845100243,
    416902256748047, 324455300753385, 170761109981736, 8431291022912,
    426721780597060, 23484981744328,  480538655455442, 359050325679881,
    444800518849271, 12567023223460,  205934588245436, 124808406729275,
    284147441592356, 367019776661797, 316050528223979, 375957039674907,
    487128768089616, 336320304813994, 466132500606615, 251423492013773,
    511526162511020, 447353315292259, 270101388153399, 501365804951712,
    137117230574403, 55696042166782,  545465267107733, 146285008599712,
    112179228517447, 76938745909585,  208464480769227, 353442478876430,
    277625379208541, 500668743679268, 298714720962488, 410459671628176,
    507448961634440, 386662356104673, 21248630989155,  223249706673423,
    276090669689405, 526730732188275, 502083139900483, 339349306057810,
    396617212295211, 129831209728404, 418031859449479, 362497905066509,
    242139454400181, 64512067787428,  104096906620971, 233248286118624,
    395625175838205, 220706181998868, 178011780558488, 297608875847049,
    56300677758822,  411686396018384, 446433637531621, 424902650816766,
    471232381855058, 11480274684976,  66255131154300,  456867515230800,
    31323128470944,  89259333705787,  282664902181958, 132101047255159,
    2130660490928,   225139101731979, 378253454469362, 50580845857958,
    455093890927570, 251297637847462, 362518083378031, 458466212144218,
    282201701100770, 520171521062184, 197112012111216, 252110537867508,
    463941921169400, 326169317056317, 450335195193715, 481145328536034,
    516238661511508, 18998775653199,  488725981008283, 516554533210258,
    388613415295435, 122565080659048, 336609986081208, 547754322522410,
    457289401879026, 490743670357629, 65607711716482,  524470253535877,
    113617727931557, 31436313021226,  390501739822830, 87378501258482,
    185748637869801, 142101801786678, 254094604641235, 258433476747023,
    76556381512981,  343564511406371, 300217343872725, 194710678590173,
    138930253490310, 202013898798094, 502861168844488, 418158273978293,
    8183197180961,   1973852184081,   140150542983303, 116938241314859,
    451104337372646, 218983032891607, 469846812616670, 187021050421659,
    250923558679262, 151968380866200, 242228879552541, 327416975973725,
    379749068719938, 121732084852046, 329016831648138, 509991452494934,
    421512523538528, 204777905639040, 384880513489892, 448572637110610,
    307570741537882, 448099811371701, 239142294833278, 488619825633974,
    483476975795421, 58557446422668,  311914828802751, 49272608253955,
    412364941486067, 134985800185735, 296829380860897, 358392482407352,
    445820860845508, 356999441512402, 449599463217424, 459671140884604,
    372482767126226, 123795254301554, 69033636211214,  357768739580106,
    18676019291022,  306908722114695, 280253957569726, 449896578152105,
    40087191175247,  98505378694646,  479469403306358, 300957805660779,
    467044652391770, 404723813618415, 469246343748684, 233040423876983,
    106261189626398, 376100717535891, 366526811751098, 508654291449351,
    150397450393545, 346844257806196, 543404937757387, 131277337751109,
    48877741985483,  15378174789751,  417317734189149, 89420920236901,
    344448593618386, 81573808672651,  409451398424702, 354095171726010,
    195238242691812, 311389272519781, 540748375704259, 269421380536186,
    365004010201246, 249268186002513, 85706110288659,  529292973770802,
    241027961524662, 228408410982028, 96779619954579,  207610196644409,
    273464461369757, 184559744695575, 319071856053126, 132715189879535,
    323465490848965, 76821774525391,  389301739358667, 198372696828630,
    403965795506306, 77764285044059,  500861755393268, 279807245617648,
    41962022638291,  537827067611816, 37167445102122,  519660805424029,
    183143257441005, 111200915763639, 236682148757055, 198414470464350,
    8092241188896,   136915629127150, 271916721719633, 17932091883836,
    311285773713158, 62447473878422,  377193201318740, 182327859434527,
    151207547371060, 520722673915319, 106118117435405, 81679115335636,
    522180437629763, 498926532800398, 419410125985911, 485404179032773,
    132197727217027, 14129928545491,  118871705078881, 65273163517139,
    359884895202291, 425775782986874, 123621257503745, 539830154081926,
    172916720253805, 325434057902356, 7689206981481,   548159825584473,
    6377613068405,   28391976233608,  173208810046062, 181736852155961,
    310324961443086, 398266946742291, 451883604750610, 417422724956278,
    490272312147488, 160972484954135, 490304981099835, 339838900206265,
    289390782209915, 204400439950910, 10276200581659,  270712993723384,
    452348841055289, 453402064308017, 321596800408166, 390805527875731,
    266703067060790, 338373242294334, 86246683807021,  9882114944603,
    102810217166605, 199447235925556, 4536260255266,   513983304167486,
    372038633526341, 81267445502953,  189618291911104, 97805238992534,
    225727196151879, 62449389214119,  219226695843158, 282848852690504,
    540337360809628, 324633569793958, 545290044295941, 186191092042498,
    161675321641364, 491504876037951, 413884745894778, 323627212742799,
    118642261844448, 403988630153751, 470675365723139};

//  Forward NTT (negacyclic -- evaluate polynomial at factors of x^n+1).

void polyr_fntt(int64_t *v)
{
    size_t i, j, k;
    int64_t x, y, z;
    int64_t *p0, *p1, *p2;

    const int64_t *w = racc_w_64;

    for (k = 1, j = RACC_N >> 1; j > 0; k <<= 1, j >>= 1) {

        p0 = v;
        for (i = 0; i < k; i++) {
            z = *w++;
            p1 = p0 + j;
            p2 = p1 + j;

            while (p1 < p2) {
                x = *p0;
                y = *p1;
                y = mont64_mulq(y, z);
                *p0++ = mont64_add(x, y);
                *p1++ = mont64_sub(x, y);
            }
            p0 = p2;
        }
    }
}

//  Reverse NTT (negacyclic -- x^n+1), normalize by 1/(n*r).

void polyr_intt(int64_t *v)
{
    size_t i, j, k;
    int64_t x, y, z;
    int64_t *p0, *p1, *p2;

    const int64_t *w = &racc_w_64[RACC_N - 2];

    for (j = 1, k = RACC_N >> 1; k > 0; j <<= 1, k >>= 1) {

        p0 = v;

        for (i = 0; i < k; i++) {
            z = *w--;
            p1 = p0 + j;
            p2 = p1 + j;

            while (p1 < p2) {
                x = *p0;
                y = *p1;
                *p0++ = mont64_add(x, y);
                y = mont64_sub(y, x);
                *p1++ = mont64_mulq(y, z);
            }
            p0 = p2;
        }
    }

    //  normalization
    polyr_ntt_smul(v, v, MONT_NI);
}

//  Scalar multiplication, Montgomery reduction.

void polyr_ntt_smul(int64_t *r, const int64_t *a, int64_t c)
{
    size_t i;

    for (i = 0; i < RACC_N; i++) {
        r[i] = mont64_cadd(mont64_mulq(a[i], c), RACC_Q);
    }
}

//  Coefficient multiply:  r = a * b,  Montgomery reduction.

void polyr_ntt_cmul(int64_t *r, const int64_t *a, const int64_t *b)
{
    size_t i;

    for (i = 0; i < RACC_N; i++) {
        r[i] = mont64_cadd(mont64_mulq(a[i], b[i]), RACC_Q);
    }
}

//  Coefficient multiply and add:  r = a * b + c, Montgomery reduction.

void polyr_ntt_mula(int64_t *r, const int64_t *a, const int64_t *b,
                    const int64_t *c)
{
    size_t i;

    for (i = 0; i < RACC_N; i++) {
        r[i] = mont64_csub(mont64_cadd(mont64_mulq(a[i], b[i]), RACC_Q) + c[i],
                           RACC_Q);
    }
}

//  POLYR_Q32
#endif