Использование% с SSE2? - PullRequest
Новая
       125

Использование% с SSE2?

0 голосов
/ 02 января 2019

Вот код, который я пытаюсь преобразовать в SSE2: 

double *pA = a;
double *pB = b[voiceIndex];
double *pC = c[voiceIndex];
double *left = audioLeft;
double *right = audioRight;
double phase = 0.0;
double bp0 = mNoteFrequency * mHostPitch;

for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex++) {
    // some other code (that will use phase)

    phase += std::clamp(mRadiansPerSample * (bp0 * pB[sampleIndex] + pC[sampleIndex]), 0.0, PI);

    while (phase >= TWOPI) { phase -= TWOPI; }
}

Вот что я достиг: 

double *pA = a;
double *pB = b[voiceIndex];
double *pC = c[voiceIndex];
double *left = audioLeft;
double *right = audioRight;
double phase = 0.0;
double bp0 = mNoteFrequency * mHostPitch;

__m128d v_boundLower = _mm_set1_pd(0.0);
__m128d v_boundUpper = _mm_set1_pd(PI);
__m128d v_bp0 = _mm_set1_pd(bp0);
__m128d v_radiansPerSample = _mm_set1_pd(mRadiansPerSample);

__m128d v_phase = _mm_set1_pd(phase);
__m128d v_pB = _mm_load_pd(pB);
__m128d v_pC = _mm_load_pd(pC);
__m128d v_result = _mm_mul_pd(v_bp0, v_pB);
v_result = _mm_add_pd(v_result, v_pC);
v_result = _mm_mul_pd(v_result, v_radiansPerSample);
v_result = _mm_max_pd(v_result, v_boundLower);
v_result = _mm_min_pd(v_result, v_boundUpper);

for (int sampleIndex = 0; sampleIndex < roundintup8(blockSize); sampleIndex += 8, pB += 8, pC += 8) {
    // some other code (that will use v_phase)

    v_phase = _mm_add_pd(v_phase, v_result);

    v_pB = _mm_load_pd(pB + 2);
    v_pC = _mm_load_pd(pC + 2);
    v_result = _mm_mul_pd(v_bp0, v_pB);
    v_result = _mm_add_pd(v_result, v_pC);
    v_result = _mm_mul_pd(v_result, v_radiansPerSample);
    v_result = _mm_max_pd(v_result, v_boundLower);
    v_result = _mm_min_pd(v_result, v_boundUpper);
    v_phase = _mm_add_pd(v_phase, v_result);

    v_pB = _mm_load_pd(pB + 4);
    v_pC = _mm_load_pd(pC + 4);
    v_result = _mm_mul_pd(v_bp0, v_pB);
    v_result = _mm_add_pd(v_result, v_pC);
    v_result = _mm_mul_pd(v_result, v_radiansPerSample);
    v_result = _mm_max_pd(v_result, v_boundLower);
    v_result = _mm_min_pd(v_result, v_boundUpper);
    v_phase = _mm_add_pd(v_phase, v_result);

    v_pB = _mm_load_pd(pB + 6);
    v_pC = _mm_load_pd(pC + 6);
    v_result = _mm_mul_pd(v_bp0, v_pB);
    v_result = _mm_add_pd(v_result, v_pC);
    v_result = _mm_mul_pd(v_result, v_radiansPerSample);
    v_result = _mm_max_pd(v_result, v_boundLower);
    v_result = _mm_min_pd(v_result, v_boundUpper);
    v_phase = _mm_add_pd(v_phase, v_result);

    v_pB = _mm_load_pd(pB + 8);
    v_pC = _mm_load_pd(pC + 8);
    v_result = _mm_mul_pd(v_bp0, v_pB);
    v_result = _mm_add_pd(v_result, v_pC);
    v_result = _mm_mul_pd(v_result, v_radiansPerSample);
    v_result = _mm_max_pd(v_result, v_boundLower);
    v_result = _mm_min_pd(v_result, v_boundUpper);

    // ... fmod?
}

Но я не совсем уверен, как заменитьwhile (phase >= TWOPI) { phase -= TWOPI; } (что по сути является классическим fmod в C ++).

Есть какие-нибудь причудливые свойства?Не могу найти ни одного в этом списке .Подразделение + какая-то ракетная бит-сдвиг?

1 Ответ

0 голосов
/ 02 января 2019

Как говорится в комментариях, похоже, что в этом вы можете сделать это просто замаскированным вычитанием со сравнением + andpd.Это работает до тех пор, пока вы никогда не сможете отойти от желаемого диапазона более чем на одно вычитание.

Как 

const __m128d v2pi = _mm_set1_pd(TWOPI);


__m128d needs_range_reduction = _mm_cmpge_pd(vphase, v2pi);
__m128d offset = _mm_and_pd(needs_range_reduction, v2pi);  // 0.0 or 2*Pi
vphase = _mm_sub_pd(vphase, offset);

Для реализации фактического (медленного) fmod не слишком заботясь о последних нескольких битах значения, вы должны сделать integer_quotient = floor(x/y) (или, может быть, rint(x/y) или ceil), а затем x - y * integer_quotient.floor / rint / ceil дешевы с SSE4.1 _mm_round_pd или _mm_floor_pd().Это даст вам остаток, который может быть отрицательным, как при целочисленном делении.

Я уверен, что существуют численные методы, которые лучше избегают ошибок округления, прежде чем катастрофическая отмена вычитает два соседних числа.Если вы заботитесь о точности, проверьте.(Использование double векторов, когда вас не очень заботит точность, довольно глупо; можно также использовать float и получать вдвое больше работы на каждый вектор).Если входное значение намного больше, чем модуль, неизбежна потеря точности, и минимизация ошибки округления во временной области, вероятно, очень важна.Но в противном случае точность будет проблемой только в том случае, если вы не заботитесь об относительной ошибке в результатах, очень близких к нулю, когда x является почти точным кратным y.(Результат, близкий к нулю, только немногие нижние биты значения и значения оставлены для точности.)

Без SSE4.1 существуют хитрости, такие как сложение и вычитание достаточно большого числа.Преобразование в целое число и обратно еще хуже для pd, потому что инструкция упакованного преобразования также декодирует некоторые случайные операции.Не говоря уже о том, что 32-разрядное целое число не покрывает весь диапазон double, но вам не хватает точности уменьшения диапазона, если ваш ввод был таким огромным.

Если у вас FMA , вы можете избежать ошибки округления в y * integer_quotient части умножения и подпункта._mm_fmsub_pd.

...