include/flint/acb_dft.h

0001 /*
0002     Copyright (C) 2016 Pascal Molin
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 ACB_DFT_H
0013 #define ACB_DFT_H
0014
0015 #include "limb_types.h"
0016 #include "acb.h"
0017
0018 #ifdef __cplusplus
0019 extern "C" {
0020 #endif
0021
0022 void _acb_dft_naive(acb_ptr w, acb_srcptr v, slong dv, acb_srcptr z, slong dz, slong len, slong prec);
0023 void acb_dft_naive(acb_ptr w, acb_srcptr v, slong len, slong prec);
0024 void acb_dft_crt(acb_ptr w, acb_srcptr v, slong len, slong prec);
0025 void acb_dft_cyc(acb_ptr w, acb_srcptr v, slong len, slong prec);
0026 void acb_dft_inverse_cyc(acb_ptr w, acb_srcptr v, slong len, slong prec);
0027 void acb_dft_rad2_inplace(acb_ptr v, int e, slong prec);
0028 void acb_dft_rad2(acb_ptr w, acb_srcptr v, int e, slong prec);
0029 void acb_dft_bluestein(acb_ptr w, acb_srcptr v, slong len, slong prec);
0030 void acb_dft_prod(acb_ptr w, acb_srcptr v, slong * cyc, slong num, slong prec);
0031
0032 void acb_dft_rad2_inplace_threaded(acb_ptr v, int e, slong prec);
0033
0034 void acb_dft_convol_naive(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec);
0035 void acb_dft_convol_dft(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec);
0036 void acb_dft_convol_rad2(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec);
0037 void acb_dft_convol_mullow(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec);
0038 void acb_dft_convol(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec);
0039
0040 #define CRT_MAX 15
0041 typedef struct
0042 {
0043     int num;
0044     nmod_t n;
0045     slong m[CRT_MAX];
0046     ulong M[CRT_MAX];
0047     ulong vM[CRT_MAX];
0048 }
0049 crt_struct;
0050
0051 typedef crt_struct crt_t[1];
0052
0053 void crt_init(crt_t c, ulong n);
0054
0055 void crt_decomp(acb_ptr y, acb_srcptr x, slong dx, const crt_t c, ulong len);
0056 void crt_recomp(acb_ptr y, acb_srcptr x, const crt_t c, ulong len);
0057
0058 typedef struct acb_dft_step_struct acb_dft_step_struct;
0059 typedef acb_dft_step_struct * acb_dft_step_ptr;
0060
0061 typedef struct
0062 {
0063     slong n;
0064     acb_ptr z;
0065     int zclear;
0066     slong num;
0067     acb_dft_step_ptr cyc;
0068 }
0069 acb_dft_cyc_struct;
0070
0071 typedef acb_dft_cyc_struct acb_dft_cyc_t[1];
0072
0073 typedef struct
0074 {
0075     int e;
0076     slong n; /* = 1 << e */
0077     slong dv;
0078     slong nz; /* = n but could be bigger */
0079     acb_ptr z;
0080 }
0081 acb_dft_rad2_struct;
0082
0083 typedef acb_dft_rad2_struct acb_dft_rad2_t[1];
0084
0085 typedef struct
0086 {
0087     slong n;
0088     slong dv;
0089     acb_ptr z; /* z[k] = e(k^2/2n) */
0090     acb_ptr g; /* g[k] = dft( z ) */
0091     acb_dft_rad2_t rad2;
0092 }
0093 acb_dft_bluestein_struct;
0094
0095 typedef acb_dft_bluestein_struct acb_dft_bluestein_t[1];
0096
0097 typedef struct
0098 {
0099     slong n;
0100     slong num;
0101     acb_dft_step_ptr cyc;
0102 }
0103 acb_dft_prod_struct;
0104
0105 typedef acb_dft_prod_struct acb_dft_prod_t[1];
0106
0107 typedef struct
0108 {
0109     slong n;
0110     crt_t c;
0111     slong dv;
0112     /* then a product */
0113     acb_dft_step_ptr cyc;
0114 }
0115 acb_dft_crt_struct;
0116
0117 typedef acb_dft_crt_struct acb_dft_crt_t[1];
0118
0119 typedef struct
0120 {
0121     slong n;
0122     slong dv;
0123     int zclear;
0124     acb_ptr z;
0125     slong dz;
0126 }
0127 acb_dft_naive_struct;
0128
0129 typedef acb_dft_naive_struct acb_dft_naive_t[1];
0130
0131 typedef struct
0132 {
0133     slong n;
0134     int type;
0135     union
0136     {
0137         acb_dft_rad2_t rad2;
0138         acb_dft_cyc_t cyc;
0139         acb_dft_prod_t prod;
0140         acb_dft_crt_t crt;
0141         acb_dft_naive_t naive;
0142         acb_dft_bluestein_t bluestein;
0143     } t;
0144 }
0145 acb_dft_pre_struct;
0146
0147 typedef acb_dft_pre_struct acb_dft_pre_t[1];
0148
0149 /* covers both product and cyclic case */
0150 struct
0151 acb_dft_step_struct
0152 {
0153     /* [G:H] */
0154     slong m;
0155     /* card H */
0156     slong M;
0157     slong dv; /* = M for prod, also = M if cyc is reordered */
0158     /* pointer on some roots of unity, if needed */
0159     acb_srcptr z;
0160     /* index of mM in z */
0161     slong dz;
0162     /* dft to call on H */
0163     acb_dft_pre_t pre;
0164     /* dft to call on G/H ? */
0165 };
0166
0167 /*typedef acb_dft_pre_struct acb_dft_pre_t[1];*/
0168
0169 #define DFT_VERB 0
0170
0171 enum
0172 {
0173     DFT_NAIVE, DFT_CYC, DFT_PROD, DFT_CRT , DFT_RAD2 , DFT_CONV
0174 };
0175
0176 void acb_dft_step(acb_ptr w, acb_srcptr v, acb_dft_step_ptr cyc, slong num, slong prec);
0177
0178 void acb_dft_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, slong prec);
0179 void acb_dft_inverse_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, slong prec);
0180 void acb_dft_naive_precomp(acb_ptr w, acb_srcptr v, const acb_dft_naive_t pol, slong prec);
0181 void acb_dft_cyc_precomp(acb_ptr w, acb_srcptr v, const acb_dft_cyc_t cyc, slong prec);
0182
0183 void acb_dft_rad2_precomp_inplace(acb_ptr v, const acb_dft_rad2_t rad2, slong prec);
0184 void acb_dft_rad2_precomp(acb_ptr w, acb_srcptr v, const acb_dft_rad2_t rad2, slong prec);
0185 void acb_dft_crt_precomp(acb_ptr w, acb_srcptr v, const acb_dft_crt_t crt, slong prec);
0186 void acb_dft_prod_precomp(acb_ptr w, acb_srcptr v, const acb_dft_prod_t prod, slong prec);
0187 void acb_dft_bluestein_precomp(acb_ptr w, acb_srcptr v, const acb_dft_bluestein_t t, slong prec);
0188
0189 void acb_dft_rad2_precomp_inplace_threaded(acb_ptr v, const acb_dft_rad2_t rad2, slong prec);
0190
0191 void acb_dft_inverse_rad2_precomp_inplace(acb_ptr v, const acb_dft_rad2_t rad2, slong prec);
0192 void acb_dft_inverse_rad2_precomp(acb_ptr w, acb_srcptr v, const acb_dft_rad2_t rad2, slong prec);
0193 void acb_dft_convol_rad2_precomp(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, const acb_dft_rad2_t rad2, slong prec);
0194
0195 void _acb_dft_precomp_init(acb_dft_pre_t pre, slong dv, acb_ptr z, slong dz, slong len, slong prec);
0196 void acb_dft_precomp_init(acb_dft_pre_t pre, slong len, slong prec);
0197 void acb_dft_precomp_clear(acb_dft_pre_t pre);
0198
0199 void acb_dft(acb_ptr w, acb_srcptr v, slong len, slong prec);
0200 void acb_dft_inverse(acb_ptr w, acb_srcptr v, slong len, slong prec);
0201
0202 acb_dft_step_ptr _acb_dft_steps_prod(slong * m, slong num, slong prec);
0203
0204 FLINT_FORCE_INLINE void
0205 acb_dft_prod_init(acb_dft_prod_t t, slong * cyc, slong num, slong prec)
0206 {
0207     t->num = num;
0208     t->cyc = _acb_dft_steps_prod(cyc, num, prec);
0209 }
0210
0211 void acb_dft_prod_clear(acb_dft_prod_t t);
0212
0213 void _acb_dft_cyc_init_z_fac(acb_dft_cyc_t t, n_factor_t fac, slong dv, acb_ptr z, slong dz, slong len, slong prec);
0214 void _acb_dft_cyc_init(acb_dft_cyc_t t, slong dv, slong len, slong prec);
0215
0216 FLINT_FORCE_INLINE void
0217 acb_dft_cyc_init(acb_dft_cyc_t t, slong len, slong prec)
0218 {
0219     _acb_dft_cyc_init(t, 1, len, prec);
0220 }
0221
0222 void acb_dft_cyc_clear(acb_dft_cyc_t t);
0223
0224 void _acb_dft_naive_init(acb_dft_naive_t pol, slong dv, acb_ptr z, slong dz, slong len, slong prec);
0225
0226 FLINT_FORCE_INLINE void
0227 acb_dft_naive_init(acb_dft_naive_t pol, slong len, slong prec)
0228 {
0229     _acb_dft_naive_init(pol, 1, NULL, 0, len, prec);
0230 }
0231
0232 FLINT_FORCE_INLINE void
0233 acb_dft_naive_clear(acb_dft_naive_t pol)
0234 {
0235     if (pol->zclear)
0236         _acb_vec_clear(pol->z, pol->n);
0237 }
0238
0239 void _acb_dft_rad2_init(acb_dft_rad2_t t, slong dv, int e, slong prec);
0240
0241 FLINT_FORCE_INLINE void
0242 acb_dft_rad2_init(acb_dft_rad2_t t, int e, slong prec)
0243 {
0244     _acb_dft_rad2_init(t, 1, e, prec);
0245 }
0246
0247 FLINT_FORCE_INLINE void
0248 acb_dft_rad2_clear(acb_dft_rad2_t t)
0249 {
0250     _acb_vec_clear(t->z, t->nz);
0251 }
0252
0253 void _acb_dft_bluestein_init(acb_dft_bluestein_t t, slong dv, slong n, slong prec);
0254
0255 FLINT_FORCE_INLINE void
0256 acb_dft_bluestein_init(acb_dft_bluestein_t t, slong n, slong prec)
0257 {
0258     _acb_dft_bluestein_init(t, 1, n, prec);
0259 }
0260
0261 FLINT_FORCE_INLINE void
0262 acb_dft_bluestein_clear(acb_dft_bluestein_t t)
0263 {
0264     if (t->n != 0)
0265     {
0266         _acb_vec_clear(t->z, t->n);
0267         _acb_vec_clear(t->g, t->rad2->n);
0268         acb_dft_rad2_clear(t->rad2);
0269     }
0270 }
0271
0272 void _acb_dft_crt_init(acb_dft_crt_t crt, slong dv, slong len, slong prec);
0273 void acb_dft_crt_init(acb_dft_crt_t crt, slong len, slong prec);
0274 void acb_dft_crt_clear(acb_dft_crt_t crt);
0275
0276 /* utils, could be moved elsewhere */
0277
0278 FLINT_FORCE_INLINE void
0279 acb_swap_ri(acb_t x)
0280 {
0281     arb_swap(acb_realref(x), acb_imagref(x));
0282 }
0283
0284 FLINT_FORCE_INLINE void
0285 acb_vec_swap_ri(acb_ptr v, slong len)
0286 {
0287     slong k;
0288     for (k = 0; k < len; k++)
0289         acb_swap_ri(v + k);
0290 }
0291
0292 FLINT_FORCE_INLINE void
0293 _acb_vec_kronecker_mul(acb_ptr z, acb_srcptr x, acb_srcptr y, slong len, slong prec)
0294 {
0295     slong k;
0296     for (k = 0; k < len; k++)
0297         acb_mul(z + k, x + k, y + k, prec);
0298 }
0299
0300 FLINT_FORCE_INLINE void
0301 _acb_vec_kronecker_mul_step(acb_ptr z, acb_srcptr x, acb_srcptr y, slong step, slong len, slong prec)
0302 {
0303     slong k;
0304     for (k = 0; k < len; k++)
0305     {
0306         acb_mul(z + k, x + k, y + k * step, prec);
0307     }
0308 }
0309
0310 #ifdef __cplusplus
0311 }
0312 #endif
0313
0314 #endif