SEAL/seal2/smallntt.cpp
// Inverse negacyclic NTT using Harvey's butterfly. (See Patrick Longa and Michael Naehrig). 
Inverse negacyclic NTT using Harvey's butterfly. (See Patrick Longa and Michael Naehrig).
Inverse negacyclic NTT using Harvey ' s butterfly . ( See Patrick Longa and Michael Naehrig ) .
void inverse_ntt_negacyclic_harvey_lazy(uint64_t *operand, const SmallNTTTables &tables)
        {
            uint64_t modulus = tables.modulus().value();
            uint64_t two_times_modulus = modulus * 2;
            int n = 1 << tables.coeff_count_power();
            int t = 1;
            for (int m = n; m > 1; m >>= 1)
            {
                int j1 = 0;
                int h = m >> 1;
                if (t >= 4)
                {
                    for (int i = 0; i < h; i++)
                    {
                        int j2 = j1 + t;
                        const uint64_t W = tables.get_from_inv_root_powers_div_two(h + i);
                        const uint64_t Wprime = tables.get_from_scaled_inv_root_powers_div_two(h + i);
                        uint64_t *U = operand + j1;
                        uint64_t *V = U + t;
                        uint64_t currU;
                        uint64_t T;
                        uint64_t H;
                        for (int j = j1; j < j2; j += 4)
                        {
                            T = two_times_modulus - *V + *U;
                            currU = *U + *V - (two_times_modulus & static_cast<uint64_t>(-static_cast<int64_t>((*U << 1) >= T)));
                            *U++ = (currU + (modulus & static_cast<uint64_t>(-static_cast<int64_t>(T & 1)))) >> 1;
                            multiply_uint64_hw64(Wprime, T, &H);
                            *V++ = T * W - H * modulus;
                            T = two_times_modulus - *V + *U;
                            currU = *U + *V - (two_times_modulus & static_cast<uint64_t>(-static_cast<int64_t>((*U << 1) >= T)));
                            *U++ = (currU + (modulus & static_cast<uint64_t>(-static_cast<int64_t>(T & 1)))) >> 1;
                            multiply_uint64_hw64(Wprime, T, &H);
                            *V++ = T * W - H * modulus;
                            T = two_times_modulus - *V + *U;
                            currU = *U + *V - (two_times_modulus & static_cast<uint64_t>(-static_cast<int64_t>((*U << 1) >= T)));
                            *U++ = (currU + (modulus & static_cast<uint64_t>(-static_cast<int64_t>(T & 1)))) >> 1;
                            multiply_uint64_hw64(Wprime, T, &H);
                            *V++ = T * W - H * modulus;
                            T = two_times_modulus - *V + *U;
                            currU = *U + *V - (two_times_modulus & static_cast<uint64_t>(-static_cast<int64_t>((*U << 1) >= T)));
                            *U++ = (currU + (modulus & static_cast<uint64_t>(-static_cast<int64_t>(T & 1)))) >> 1;
                            multiply_uint64_hw64(Wprime, T, &H);
                            *V++ = T * W - H * modulus;
                        }
                        j1 += (t << 1);
                    }
                }
                else
                {
                    for (int i = 0; i < h; i++)
                    {
                        int j2 = j1 + t;
                        const uint64_t W = tables.get_from_inv_root_powers_div_two(h + i);
                        const uint64_t Wprime = tables.get_from_scaled_inv_root_powers_div_two(h + i);
                        uint64_t *U = operand + j1;
                        uint64_t *V = U + t;
                        uint64_t currU;
                        uint64_t T;
                        uint64_t H;
                        for (int j = j1; j < j2; j++)
                        {
                            T = two_times_modulus - *V + *U;
                            currU = *U + *V - (two_times_modulus & static_cast<uint64_t>(-static_cast<int64_t>((*U << 1) >= T)));
                            *U++ = (currU + (modulus & static_cast<uint64_t>(-static_cast<int64_t>(T & 1)))) >> 1;
                            multiply_uint64_hw64(Wprime, T, &H);
                            *V++ = W * T - H * modulus;
                        }
                        j1 += (t << 1);
                    }
                }
                t <<= 1;
            }
        }
void inverse_ntt_negacyclic_harvey_lazy ( uint64_t * operand , const SmallNTTTables & tables ) { uint64_t modulus = tables . modulus ( ) . value ( ) ; uint64_t two_times_modulus = modulus * 2 ; int n = 1 << tables . coeff_count_power ( ) ; int t = 1 ; for ( int m = n ; m > 1 ; m >>= 1 ) { int j1 = 0 ; int h = m >> 1 ; if ( t >= 4 ) { for ( int i = 0 ; i < h ; i ++ ) { int j2 = j1 + t ; const uint64_t W = tables . get_from_inv_root_powers_div_two ( h + i ) ; const uint64_t Wprime = tables . get_from_scaled_inv_root_powers_div_two ( h + i ) ; uint64_t * U = operand + j1 ; uint64_t * V = U + t ; uint64_t currU ; uint64_t T ; uint64_t H ; for ( int j = j1 ; j < j2 ; j += 4 ) { T = two_times_modulus - * V + * U ; currU = * U + * V - ( two_times_modulus & static_cast < uint64_t > ( - static_cast < int64_t > ( ( * U << 1 ) >= T ) ) ) ; * U ++ = ( currU + ( modulus & static_cast < uint64_t > ( - static_cast < int64_t > ( T & 1 ) ) ) ) >> 1 ; multiply_uint64_hw64 ( Wprime , T , & H ) ; * V ++ = T * W - H * modulus ; T = two_times_modulus - * V + * U ; currU = * U + * V - ( two_times_modulus & static_cast < uint64_t > ( - static_cast < int64_t > ( ( * U << 1 ) >= T ) ) ) ; * U ++ = ( currU + ( modulus & static_cast < uint64_t > ( - static_cast < int64_t > ( T & 1 ) ) ) ) >> 1 ; multiply_uint64_hw64 ( Wprime , T , & H ) ; * V ++ = T * W - H * modulus ; T = two_times_modulus - * V + * U ; currU = * U + * V - ( two_times_modulus & static_cast < uint64_t > ( - static_cast < int64_t > ( ( * U << 1 ) >= T ) ) ) ; * U ++ = ( currU + ( modulus & static_cast < uint64_t > ( - static_cast < int64_t > ( T & 1 ) ) ) ) >> 1 ; multiply_uint64_hw64 ( Wprime , T , & H ) ; * V ++ = T * W - H * modulus ; T = two_times_modulus - * V + * U ; currU = * U + * V - ( two_times_modulus & static_cast < uint64_t > ( - static_cast < int64_t > ( ( * U << 1 ) >= T ) ) ) ; * U ++ = ( currU + ( modulus & static_cast < uint64_t > ( - static_cast < int64_t > ( T & 1 ) ) ) ) >> 1 ; multiply_uint64_hw64 ( Wprime , T , & H ) ; * V ++ = T * W - H * modulus ; } j1 += ( t << 1 ) ; } } else { for ( int i = 0 ; i < h ; i ++ ) { int j2 = j1 + t ; const uint64_t W = tables . get_from_inv_root_powers_div_two ( h + i ) ; const uint64_t Wprime = tables . get_from_scaled_inv_root_powers_div_two ( h + i ) ; uint64_t * U = operand + j1 ; uint64_t * V = U + t ; uint64_t currU ; uint64_t T ; uint64_t H ; for ( int j = j1 ; j < j2 ; j ++ ) { T = two_times_modulus - * V + * U ; currU = * U + * V - ( two_times_modulus & static_cast < uint64_t > ( - static_cast < int64_t > ( ( * U << 1 ) >= T ) ) ) ; * U ++ = ( currU + ( modulus & static_cast < uint64_t > ( - static_cast < int64_t > ( T & 1 ) ) ) ) >> 1 ; multiply_uint64_hw64 ( Wprime , T , & H ) ; * V ++ = W * T - H * modulus ; } j1 += ( t << 1 ) ; } } t <<= 1 ; } }
Inverse negacyclic NTT using Harvey's butterfly.
Inverse negacyclic NTT using Harvey ' s butterfly .
// return the bit-reversed order of NTT.  // Need the powers of phi^{-1} in bit-reversed order // Need the powers of  phi^{-1} in bit-reversed order // U = x[i], V = x[i+m] // Compute U - V + 2q // Cleverly check whether currU + currV >= two_times_modulus // Need to make it so that div2_uint_mod takes values that are > q.  //div2_uint_mod(U, modulusptr, coeff_uint64_count, U);  // We use also the fact that parity of currU is same as parity of T. // Since our modulus is always so small that currU + masked_modulus < 2^64, // we never need to worry about wrapping around when adding masked_modulus. //uint64_t masked_modulus = modulus & static_cast<uint64_t>(-static_cast<int64_t>(T & 1)); //uint64_t carry = add_uint64(currU, masked_modulus, 0, &currU); //currU += modulus & static_cast<uint64_t>(-static_cast<int64_t>(T & 1)); // effectively, the next two multiply perform multiply modulo beta = 2**wordsize. 