Быстрая реализация реализации двоичного возведения в степень в OpenCL

Я пытался разработать быструю реализацию двоичного возведения в степень в OpenCL. Моя текущая реализация очень похожа на ту, что в этой книге о пи.

// Returns 16^n mod ak
inline double expm (long n, double ak)
{
    double r = 16.0;
    long nt;

    if (ak == 1) return 0.;
    if (n == 0) return 1;
    if (n == 1) return fmod(16.0, ak);

    for (nt=1; nt <= n; nt <<=1);

    nt >>= 2;

    do
    {
        r = fmod(r*r, ak);
        if ((n & nt) != 0)
            r = fmod(16.0*r, ak);
        nt >>= 1;
    } while (nt != 0);
    return r;
}

Есть ли место для улучшения? Прямо сейчас моя программа тратит подавляющее большинство своего времени на эту функцию.

Источник

2 ответа

Решение

Моя первая мысль состоит в том, чтобы векторизовать его, для потенциальной скорости ~1,6x. При этом используется 5 умножений на цикл по сравнению с 2 умножениями в оригинале, но примерно с четвертью числа циклов для достаточно большого N. Преобразование всех doubleс longс, и обменять fmodс для %Может обеспечить некоторую скорость в зависимости от используемого графического процессора и т. д.

inline double expm(long n, double ak) {

    double4 r = (1.0, 1.0, 1.0, 1.0);
    long4 ns = n & (0x1111111111111111, 0x2222222222222222, 0x4444444444444444,
            0x8888888888888888);
    long nt;

    if(ak == 1) return 0.;

    for(nt=15; nt<n; nt<<=4); //This can probably be vectorized somehow as well.

    do {
        double4 tmp = r*r;
        tmp = tmp*tmp;
        tmp = tmp*tmp;
        r = fmod(tmp*tmp, ak); //Raise it to the 16th power, 
                                       //same as multiplying the exponent 
                                       //(of the result) by 16, same as
                                       //bitshifting the exponent to the right 4 bits.

        r = select(fmod(r*(16.0,256.0,65536.0, 4294967296.0), ak), r, (ns & nt) - 1);
        nt >>= 4;
    } while(nt != 0); //Process n four bits at a time.

    return fmod(r.x*r.y*r.z*r.w, ak); //And then combine all of them.
}

Изменить: я уверен, что это работает сейчас.

Источник
  • Цикл для извлечения nt = log2(n); можно заменить на
    if (n & 1) ...; n >>= 1;
    в цикле do-while.
  • Учитывая, что изначально r = 16;, fmod(r*r, ak) против fmod(16*r,ak) можно легко отложить, чтобы вычислить по модулю только каждую N-ю итерацию или около того - развертывание цикла?
  • И почему фмод?
Источник
Другие вопросы по тегам