<commit_before>/// Implementation of the Number-Theoretical Transform in GF(P) Galois field
#include <iostream>
#include <algorithm>
#include <stdint.h>
#include <string.h>

#if defined(_M_X64) || defined(_M_AMD64) || defined(__x86_64__)
#define MY_CPU_AMD64
#endif

#if defined(MY_CPU_AMD64) || defined(_M_IA64)
#define MY_CPU_64BIT
#endif


/***********************************************************************************************************************
*** GF(P) **************************************************************************************************************
************************************************************************************************************************/

typedef uint32_t ElementT;        // data items type, 32-bit unsigned integer for GF(P) computations with P>65536
typedef uint64_t DoubleElementT;  // twice wider type to hold intermediate results


template <typename T, T P>
ElementT GF_Add (ElementT X, ElementT Y)
{
    ElementT res = X + Y;
    return res - ((res>=P)+(res<X))*P;   // (res>=P || res<X)? res-P : res;
}

template <typename T, T P>
ElementT GF_Sub (ElementT X, ElementT Y)
{
    ElementT res = X - Y;
    return res + (res>X)*P;   // res<=X? res : res+P
}

#if __GNUC__ && defined(MY_CPU_64BIT)
// Alternative GF_Mul64 implementation for GCC

#include <inttypes.h>
typedef unsigned __int128 FourElement;    // 4x wider type to hold intermediate MUL results

template <ElementT P>
ElementT GF_Mul64 (ElementT X, ElementT Y)
{
    // See chapter "16.9 Division" in the http://www.agner.org/optimize/optimizing_assembly.pdf
    int BITS = 95;  // for P in 2^31..2^32, in general bits(DoubleElementT) + rounded_down_logb(P)
    DoubleElementT extra = ((FourElement(2) << BITS) / P) & 1;    // 0;
    DoubleElementT invP  =  (FourElement(1) << BITS) / P + extra; // 9225624384409992840ULL;

    DoubleElementT res = DoubleElementT(X)*Y;
    res -= DoubleElementT(((res+1-extra)*FourElement(invP)) >> BITS) * P;
    return ElementT(res);
}

#elif _MSC_VER
// Alternative GF_Mul64 implementation made with MSVC intrinsics

template <ElementT P>
ElementT GF_Mul64 (ElementT X, ElementT Y)
{
    DoubleElementT estInvP = ((DoubleElementT(1)<<63) / P) << 1;<commit_msg>NTT: new GF_Mul64 implementation made compatible with any P value, including 65537<commit_after>/// Implementation of the Number-Theoretical Transform in GF(P) Galois field
#include <iostream>
#include <algorithm>
#include <stdint.h>
#include <string.h>

#if defined(_M_X64) || defined(_M_AMD64) || defined(__x86_64__)
#define MY_CPU_AMD64
#endif

#if defined(MY_CPU_AMD64) || defined(_M_IA64)
#define MY_CPU_64BIT
#endif


/***********************************************************************************************************************
*** GF(P) **************************************************************************************************************
************************************************************************************************************************/

typedef uint32_t ElementT;        // data items type, 32-bit unsigned integer for GF(P) computations with P>65536
typedef uint64_t DoubleElementT;  // twice wider type to hold intermediate results


template <typename T, T P>
ElementT GF_Add (ElementT X, ElementT Y)
{
    ElementT res = X + Y;
    return res - ((res>=P)+(res<X))*P;   // (res>=P || res<X)? res-P : res;
}

template <typename T, T P>
ElementT GF_Sub (ElementT X, ElementT Y)
{
    ElementT res = X - Y;
    return res + (res>X)*P;   // res<=X? res : res+P
}

#if __GNUC__ && defined(MY_CPU_64BIT)
// Alternative GF_Mul64 implementation for GCC

#include <inttypes.h>
typedef unsigned __int128 FourElement;    // 4x wider type to hold intermediate MUL results

// The binary logarithm, rounded down
template <typename T>
static constexpr int trunc_log2 (T x)
{return x<=1? 0 : 1+trunc_log2(x/2);}

template <ElementT P>
ElementT GF_Mul64 (ElementT X, ElementT Y)
{
    // See chapter "16.9 Division" in the http://www.agner.org/optimize/optimizing_assembly.pdf
    constexpr int BITS = trunc_log2(P) + 8*sizeof(DoubleElementT);
    constexpr FourElement    invP2 = (FourElement(2) << BITS) / P;  // double the invP value
    constexpr DoubleElementT invP  = (invP2+1) / 2;                 // rounded invP value
    constexpr DoubleElementT extra = 1 - (invP2 & 1);               // 1 if invP was rounded down, 0 otherwise

    DoubleElementT res = DoubleElementT(X)*Y;
    res -= DoubleElementT(((res+extra)*FourElement(invP)) >> BITS) * P;
    return ElementT(res);
}

#elif _MSC_VER
// Alternative GF_Mul64 implementation made with MSVC intrinsics

template <ElementT P>
ElementT GF_Mul64 (ElementT X, ElementT Y)
{
    DoubleElementT estInvP = ((DoubleElementT(1)<<63) / P) << 1;
{'fraction_non_alphanumeric': 0.2524617640896711, 'fraction_numerical': 0.03205531112507857, 'mean_word_length': 4.7173652694610775, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '838488', 'n_tokens_mistral': 1605, 'n_tokens_neox': 1436, 'n_words': 428}