include/flint/fft.h

0001 /*
0002     Copyright (C) 2009, 2011 William Hart
0003
0004     This file is part of FLINT.
0005
0006     FLINT is free software: you can redistribute it and/or modify it under
0007     the terms of the GNU Lesser General Public License (LGPL) as published
0008     by the Free Software Foundation; either version 3 of the License, or
0009     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
0010 */
0011
0012 #ifndef FFT_H
0013 #define FFT_H
0014
0015 #ifdef FFT_INLINES_C
0016 #define FFT_INLINE
0017 #else
0018 #define FFT_INLINE static inline
0019 #endif
0020
0021 #include <gmp.h>
0022 #include "flint.h"
0023
0024 #ifdef __cplusplus
0025 extern "C" {
0026 #endif
0027
0028 /* defined in mpn_extras.h */
0029 mp_limb_t flint_mpn_sumdiff_n(mp_ptr s, mp_ptr d, mp_srcptr x, mp_srcptr y, mp_size_t n);
0030
0031 #define fft_sumdiff(t, u, r, s, n) \
0032    (n == 0 ? 0 : flint_mpn_sumdiff_n(t, u, r, s, n))
0033
0034 #define SWAP_PTRS(xx, yy) FLINT_SWAP(mp_ptr, xx, yy)
0035
0036 /* used for generating random values mod p in test code */
0037 /* NOTE: One needs to include mpn_extras.h for this macro. */
0038 #define random_fermat(nn, state, limbs) \
0039    do { \
0040       if (n_randint(state, 10) == 0) { \
0041          flint_mpn_zero(nn, limbs); \
0042          nn[limbs] = 1; \
0043       } else { \
0044          if (n_randint(state, 2) == 0) \
0045             flint_mpn_rrandom(nn, state, limbs); \
0046          else \
0047             flint_mpn_urandomb(nn, state, limbs * FLINT_BITS); \
0048          nn[limbs] = n_randint(state, 1024); \
0049       } \
0050       if (n_randint(state, 2)) \
0051          nn[limbs] = -nn[limbs]; \
0052    } while (0)
0053
0054 FFT_INLINE
0055 void mpn_addmod_2expp1_1(mp_limb_t * r, mp_size_t limbs, mp_limb_signed_t c)
0056 {
0057    mp_limb_t sum = r[0] + c;
0058
0059    /* check if adding c would cause a carry to propagate */
0060    if ((mp_limb_signed_t)(sum ^ r[0]) >= 0)
0061       r[0] = sum;
0062    else
0063    {
0064       if (c >= 0) mpn_add_1(r, r, limbs + 1, c);
0065       else mpn_sub_1(r, r, limbs + 1, -c);
0066    }
0067 }
0068
0069 void fft_combine_limbs(mp_limb_t * res, mp_limb_t ** poly, slong length,
0070             mp_size_t coeff_limbs, mp_size_t output_limbs, mp_size_t total_limbs);
0071
0072 void fft_combine_bits(mp_limb_t * res, mp_limb_t ** poly, slong length,
0073                  flint_bitcnt_t bits, mp_size_t output_limbs, mp_size_t total_limbs);
0074
0075 mp_size_t fft_split_limbs(mp_limb_t ** poly, mp_srcptr limbs,
0076             mp_size_t total_limbs, mp_size_t coeff_limbs, mp_size_t output_limbs);
0077
0078 mp_size_t fft_split_bits(mp_limb_t ** poly, mp_srcptr limbs,
0079                  mp_size_t total_limbs, flint_bitcnt_t bits, mp_size_t output_limbs);
0080
0081 void fermat_to_mpz(mpz_t m, mp_limb_t * i, mp_size_t limbs);
0082
0083 void mpn_normmod_2expp1(mp_limb_t * t, mp_size_t limbs);
0084
0085 void butterfly_lshB(mp_limb_t * t, mp_limb_t * u, mp_limb_t * i1,
0086                        mp_limb_t * i2, mp_size_t limbs, mp_size_t x, mp_size_t y);
0087
0088 void butterfly_rshB(mp_limb_t * t, mp_limb_t * u, mp_limb_t * i1,
0089                        mp_limb_t * i2, mp_size_t limbs, mp_size_t x, mp_size_t y);
0090
0091 void mpn_negmod_2expp1(mp_limb_t* z, const mp_limb_t* a, mp_size_t limbs);
0092
0093 void mpn_mul_2expmod_2expp1(mp_limb_t * t,
0094                                   mp_limb_t * i1, mp_size_t limbs, flint_bitcnt_t d);
0095
0096 void mpn_div_2expmod_2expp1(mp_limb_t * t,
0097                                   mp_limb_t * i1, mp_size_t limbs, flint_bitcnt_t d);
0098
0099 void fft_adjust(mp_limb_t * r, mp_limb_t * i1,
0100                                      mp_size_t i, mp_size_t limbs, flint_bitcnt_t w);
0101
0102 void fft_butterfly(mp_limb_t * s, mp_limb_t * t, mp_limb_t * i1,
0103                      mp_limb_t * i2, mp_size_t i, mp_size_t limbs, flint_bitcnt_t w);
0104
0105 void ifft_butterfly(mp_limb_t * s, mp_limb_t * t, mp_limb_t * i1,
0106                      mp_limb_t * i2, mp_size_t i, mp_size_t limbs, flint_bitcnt_t w);
0107
0108 void fft_radix2(mp_limb_t ** ii,
0109                     mp_size_t n, flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2);
0110
0111 void fft_truncate1(mp_limb_t ** ii, mp_size_t n, flint_bitcnt_t w,
0112                                mp_limb_t ** t1, mp_limb_t ** t2, mp_size_t trunc);
0113
0114 void fft_truncate(mp_limb_t ** ii,  mp_size_t n, flint_bitcnt_t w,
0115                                mp_limb_t ** t1, mp_limb_t ** t2, mp_size_t trunc);
0116
0117 void ifft_radix2(mp_limb_t ** ii, mp_size_t n,
0118                                  flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2);
0119
0120 void ifft_truncate1(mp_limb_t ** ii, mp_size_t n, flint_bitcnt_t w,
0121                                mp_limb_t ** t1, mp_limb_t ** t2, mp_size_t trunc);
0122
0123 void ifft_truncate(mp_limb_t ** ii, mp_size_t n, flint_bitcnt_t w,
0124                                mp_limb_t ** t1, mp_limb_t ** t2, mp_size_t trunc);
0125
0126 void fft_butterfly_sqrt2(mp_limb_t * s, mp_limb_t * t,
0127                          mp_limb_t * i1, mp_limb_t * i2, mp_size_t i,
0128                                 mp_size_t limbs, flint_bitcnt_t w, mp_limb_t * temp);
0129
0130 void ifft_butterfly_sqrt2(mp_limb_t * s, mp_limb_t * t, mp_limb_t * i1,
0131    mp_limb_t * i2, mp_size_t i, mp_size_t limbs, flint_bitcnt_t w, mp_limb_t * temp);
0132
0133 void fft_adjust_sqrt2(mp_limb_t * r, mp_limb_t * i1,
0134                    mp_size_t i, mp_size_t limbs, flint_bitcnt_t w, mp_limb_t * temp);
0135
0136 void fft_truncate_sqrt2(mp_limb_t ** ii, mp_size_t n, flint_bitcnt_t w,
0137             mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** temp, mp_size_t trunc);
0138
0139 void ifft_truncate_sqrt2(mp_limb_t ** ii, mp_size_t n, flint_bitcnt_t w,
0140             mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** temp, mp_size_t trunc);
0141
0142 void mul_truncate_sqrt2(mp_ptr r1, mp_srcptr i1, mp_size_t n1,
0143                         mp_srcptr i2, mp_size_t n2, flint_bitcnt_t depth, flint_bitcnt_t w);
0144
0145 void fft_butterfly_twiddle(mp_limb_t * u, mp_limb_t * v,
0146    mp_limb_t * s, mp_limb_t * t, mp_size_t limbs, flint_bitcnt_t b1, flint_bitcnt_t b2);
0147
0148 void ifft_butterfly_twiddle(mp_limb_t * u, mp_limb_t * v,
0149    mp_limb_t * s, mp_limb_t * t, mp_size_t limbs, flint_bitcnt_t b1, flint_bitcnt_t b2);
0150
0151 void fft_radix2_twiddle(mp_limb_t ** ii, mp_size_t is,
0152       mp_size_t n, flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2,
0153                             mp_size_t ws, mp_size_t r, mp_size_t c, mp_size_t rs);
0154
0155 void ifft_radix2_twiddle(mp_limb_t ** ii, mp_size_t is,
0156         mp_size_t n, flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2,
0157                             mp_size_t ws, mp_size_t r, mp_size_t c, mp_size_t rs);
0158
0159 void fft_truncate1_twiddle(mp_limb_t ** ii, mp_size_t is,
0160         mp_size_t n, flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2,
0161            mp_size_t ws, mp_size_t r, mp_size_t c, mp_size_t rs, mp_size_t trunc);
0162
0163 void ifft_truncate1_twiddle(mp_limb_t ** ii, mp_size_t is,
0164         mp_size_t n, flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2,
0165            mp_size_t ws, mp_size_t r, mp_size_t c, mp_size_t rs, mp_size_t trunc);
0166
0167 void fft_mfa_truncate_sqrt2(mp_limb_t ** ii, mp_size_t n,
0168                        flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2,
0169                                 mp_limb_t ** temp, mp_size_t n1, mp_size_t trunc);
0170
0171 void ifft_mfa_truncate_sqrt2(mp_limb_t ** ii, mp_size_t n,
0172                       flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2,
0173                                 mp_limb_t ** temp, mp_size_t n1, mp_size_t trunc);
0174
0175 void mul_mfa_truncate_sqrt2(mp_ptr r1, mp_srcptr i1, mp_size_t n1,
0176                         mp_srcptr i2, mp_size_t n2, flint_bitcnt_t depth, flint_bitcnt_t w);
0177
0178 void fft_mfa_truncate_sqrt2_outer(mp_limb_t ** ii, mp_size_t n,
0179                       flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2,
0180                                 mp_limb_t ** temp, mp_size_t n1, mp_size_t trunc);
0181
0182 void fft_mfa_truncate_sqrt2_inner(mp_limb_t ** ii, mp_limb_t ** jj,
0183             mp_size_t n, flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2,
0184                 mp_limb_t ** FLINT_UNUSED(temp), mp_size_t n1, mp_size_t trunc, mp_limb_t ** tt);
0185
0186 void ifft_mfa_truncate_sqrt2_outer(mp_limb_t ** ii, mp_size_t n,
0187                         flint_bitcnt_t w, mp_limb_t ** t1, mp_limb_t ** t2,
0188                                 mp_limb_t ** temp, mp_size_t n1, mp_size_t trunc);
0189
0190 void fft_negacyclic(mp_limb_t ** ii, mp_size_t n, flint_bitcnt_t w,
0191                              mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** temp);
0192
0193 void ifft_negacyclic(mp_limb_t ** ii, mp_size_t n, flint_bitcnt_t w,
0194                              mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** temp);
0195
0196 void fft_naive_convolution_1(mp_limb_t * r, mp_limb_t * ii,
0197                                                      mp_limb_t * jj, mp_size_t m);
0198
0199 void _fft_mulmod_2expp1(mp_limb_t * r1, mp_limb_t * i1, mp_limb_t * i2,
0200                              mp_size_t r_limbs, flint_bitcnt_t depth, flint_bitcnt_t w);
0201
0202 slong fft_adjust_limbs(mp_size_t limbs);
0203
0204 void fft_mulmod_2expp1(mp_limb_t * r, mp_limb_t * i1, mp_limb_t * i2,
0205                                         mp_size_t n, mp_size_t w, mp_limb_t * tt);
0206
0207 void flint_mpn_mul_fft_main(mp_ptr r1, mp_srcptr i1, mp_size_t n1,
0208                         mp_srcptr i2, mp_size_t n2);
0209
0210 void fft_convolution_basic(mp_limb_t ** ii, mp_limb_t ** jj,
0211              slong depth, slong limbs, slong trunc, mp_limb_t ** t1,
0212                             mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t ** tt);
0213
0214 void fft_convolution(mp_limb_t ** ii, mp_limb_t ** jj, slong depth,
0215                                     slong limbs, slong trunc, mp_limb_t ** t1,
0216                             mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t ** tt);
0217
0218 /***** FFT Precaching *****/
0219
0220 void fft_precache(mp_limb_t ** jj, slong depth, slong limbs,
0221                slong trunc, mp_limb_t ** t1, mp_limb_t ** t2, mp_limb_t ** s1);
0222
0223 void fft_convolution_precache(mp_limb_t ** ii, mp_limb_t ** jj,
0224                slong depth, slong limbs, slong trunc, mp_limb_t ** t1,
0225                         mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t ** tt);
0226
0227 #ifdef __cplusplus
0228 }
0229 #endif
0230
0231 #endif