分享一下最近的成果0.0

复制内容到剪贴板

代码:

template <typename type, size_t Rank, typename Allocator>
std::complex<type> DFT(const multiarraycase<type, Rank, Allocator> &buffer, size_t idx)
{
        typedef multiarraycase<type, Rank, Allocator>::base_array base_array;
        size_t d = idx / buffer.size(0, 2);
        type factor = -(type)2 * (type)M_PI * (type)d / (type)buffer.size(0);
               
        std::complex<type> result(0, 0);
        __if_exists(base_array::type)
        {
                for(size_t i = 0; i < buffer.size(0); ++i)
                        result += DFT(buffer, idx % buffer.size(0, 2)) * polar<type>(1, (type)i * factor);
        }
        __if_not_exists(base_array::type)
        {
                for(size_t i = 0; i < buffer.size(0); ++i)
                        result += polar<type>(buffer, (type)i * factor);
        }
        return result;
}

复制内容到剪贴板

代码:

template <typename type, typename arg_type, size_t Rank>
multiarray<std::complex<type>, Rank> &_FourierTransform(const multiarray<arg_type, Rank> &buffer, const bool inverse,
        multiarray<std::complex<type>, Rank> &result)
{
        size_t size = buffer.size(0);
        result._base::reshape(buffer.shape().begin(),
                (multiarray<arg_type, Rank>::pointer)nullptr, (multiarray<arg_type, Rank>::pointer)nullptr);
        parallel_for(size, [&, inverse](size_t idx) -> void
        {
                multiarray<std::complex<type>, Rank - 1> row;
                FourierTransform<type, arg_type, Rank - 1>(buffer[idx], inverse, row);
                std::copy(row.begin(), row.end(), result[idx].begin());
        });
        return result;
}
template <typename type, typename arg_type>
multiarray<std::complex<type>> &CooleyTukey(const multiarray<arg_type> &buffer, const size_t radix, const bool inverse,
        multiarray<std::complex<type>> &result)
{
        const size_t N = buffer.size(), index[2][2] = {{radix, N / radix}, {N / radix, radix}};
        const std::complex<type> TwiddleFactor = std::polar<type>(1, (type)(inverse ? -1 : 1) * (type)2 * (type)M_PI / (type)N);

        multiarray<arg_type, 2> temp(index[1], buffer.begin(), buffer.end());
        Transpose(temp, temp)._base::reshape(index[0], (multiarray<arg_type, 2>::pointer)nullptr, (multiarray<arg_type, 2>::pointer)nullptr);

        multiarray<std::complex<type>, 2> reorder;
        Transpose(_FourierTransform<type, arg_type>(temp, inverse, reorder), reorder);
        for(size_t i = 0; i < reorder.size(0); ++i)
                for(size_t j = 0; j < reorder.size(1); ++j)
                        reorder[j] *= std::pow<type>(TwiddleFactor, (type)i * (type)j);
        Transpose(_FourierTransform<type, std::complex<type>>(reorder, inverse, reorder), reorder);
               
        result.reshape(N);
        std::copy(reorder.begin(), reorder.end(), result.begin());
        return result;
}
template <typename type, typename arg_type, size_t Rank>
multiarray<std::complex<type>, Rank> &FourierTransform(const multiarray<arg_type, Rank> &buffer, const bool inverse,
        multiarray<std::complex<type>, Rank> &result)
{
        typedef multiarray<type, Rank>::base_array base_array;
        __if_exists(base_array::type)
        {
                const size_t index[2] = {buffer.size(0), buffer.size(0, 2)};
                multiarray<std::complex<type>, 2> reorder(index[0], index[1]);
                       
                multiarray<std::complex<type>, Rank> rows;
                _FourierTransform<type, arg_type, Rank>(buffer, inverse, rows);
                std::copy(rows.begin(), rows.end(), reorder.begin());
                Transpose(reorder, reorder);
                _FourierTransform<type, std::complex<type>>(reorder, inverse, reorder);
                Transpose(reorder, reorder);

                for(size_t i = 0; i < Rank; ++i)
                        result.reshape(i + 1, buffer.size(i));
                std::copy(reorder.begin(), reorder.end(), result.begin());
        }
        __if_not_exists(base_array::type)
        {
                const size_t N = buffer.size();
                const std::complex<type> TwiddleFactor = std::polar<type>(1, (type)(inverse ? -1 : 1) * (type)2 * (type)M_PI / (type)N);
                switch(N)
                {
                case 0:
                        break;
                case 1:
                        result.reshape(N);
                        result[0] = buffer[0];
                        break;
                case 2:
                        {
                                multiarray<arg_type> temp(buffer);
                                result.reshape(N);
                                result[0] = temp[0] + temp[1];
                                result[1] = temp[0] + temp[1] * TwiddleFactor;
                        }
                        break;
                default:
                        CooleyTukey<type, arg_type>(buffer, 2, inverse, result);   //这里只要把2给改了就可以变成mixed radix...至於合数分解这方面问题还未有时间处理
                }
                if(inverse)
                {
                        for(auto &value : result)
                                value /= (type)N;
                }
        }
        return result;
}

复制内容到剪贴板

代码:

void parallel_for(const size_t count, const std::function<void(size_t)> &functor)
{
        if(count > MAX_PARALLEL_NUM)
        {
                std::thread pThread[MAX_PARALLEL_NUM];
                const size_t blocksize = count % MAX_PARALLEL_NUM ? count / MAX_PARALLEL_NUM + 1 : count / MAX_PARALLEL_NUM;
                for(size_t i = 0; i < MAX_PARALLEL_NUM; ++i)
                {
                        const size_t start = blocksize * i, end = count < start + blocksize ? count : start + blocksize;
                        pThread = std::thread([](const std::function<void(size_t)> &functor, const size_t start, const size_t end) -> void
                        {
                                for(size_t i = start; i < end; ++i)
                                        functor(i);
                        }, std::ref(functor), start, end);
                }
                for(size_t i = 0; i < MAX_PARALLEL_NUM; ++i)
                        pThread.join();
        } else
                for(size_t i = 0; i < count; ++i)
                        functor(i);
}