Модул 2*Pi с помощта на SSE/SSE2

Все още съм доста нов в използването на SSE и се опитвам да внедря модул от 2*Pi за входове с двойна точност от порядъка 1e8 (резултатът от който ще бъде въведен в някои векторизирани тригонометрични изчисления).

Сегашният ми опит за кода се основава на идеята, че mod(x, 2*Pi) = x - floor(x/(2*Pi))*2*Pi и изглежда така:

#define _PD_CONST(Name, Val)                                            \
static const double _pd_##Name[2] __attribute__((aligned(16))) = { Val, Val }  

_PD_CONST(2Pi, 6.283185307179586);  /* = 2*pi  */  
_PD_CONST(recip_2Pi, 0.159154943091895); /* = 1/(2*pi)  */

void vec_mod_2pi(const double * vec, int Size, double * modAns)
{
    __m128d sse_a, sse_b, sse_c;
    int i;
    int k = 0;
    double t = 0;

    unsigned int initial_mode;
    initial_mode = _MM_GET_ROUNDING_MODE();

    _MM_SET_ROUNDING_MODE(_MM_ROUND_DOWN);

    for (i = 0; i < Size; i += 2)
    {
        sse_a = _mm_loadu_pd(vec+i);
        sse_b = _mm_mul_pd( _mm_cvtepi32_pd( _mm_cvtpd_epi32( _mm_mul_pd(sse_a, *(__m128d*)_pd_recip_2Pi) ) ), *(__m128d*)_pd_2Pi);
        sse_c = _mm_sub_pd(sse_a, sse_b);
        _mm_storeu_pd(modAns+i,sse_c);
    }

    k = i-2;
    for (i = 0; i < Size%2; i++)
    {
        t = (double)((int)(vec[k+i] * 0.159154943091895)) * 6.283185307179586;
        modAns[k+i] = vec[k+i] - t;
    }

    _MM_SET_ROUNDING_MODE(initial_mode);
}

За съжаление, това в момента връща много NaN с няколко отговора от 1.128e119 също (донякъде извън диапазона от 0 -> 2*Pi, към който се стремях!). Подозирам, че там, където греша, е преобразуването от двойно към int-към-double, което се опитвам да използвам, за да направя floor.

Може ли някой да ми предложи къде съм сбъркал и как да го подобря?

P.S. съжалявам за формата на този код, за първи път публикувам въпрос тук и не мога да го накарам да ми даде празни редове в кодовия блок, за да го направя четим.


c optimization simd sse
person BigA    schedule 20.07.2012    source източник
comment
Може би трябва да го направя малко по-ясно - не съм много загрижен за точността на резултата тук (всъщност, ако се върна с точност с единична точност, това би било добре). Това, от което съм най-загрижен, е а) да го накарам да работи в SSE, за да мога да се отърва от ненужно хранилище/зареждане, което имам в момента между два други SSE блока и б) да го накарам да работи бързо, ако е възможно.   -  person BigA    schedule 22.07.2012
comment
Първо, използвайте _mm_set1_pd(6.283185307179586) като нормален човек; компилаторът вече е достатъчно интелигентен, за да превърне това в еквивалент на static const double[2], но с елиминиране на дублиране (като за дублиращи се низови литерали). Второ, можете да използвате _mm_cvttpd_epi32 (обърнете внимание на допълнителното t за отрязване), за да конвертирате в цяло число с отрязване вместо текущия режим на закръгляване. Ако вашите числа са положителни, това е същото като floor. Ако имате SSE4.1, можете да използвате _mm_floor_pd (roundpd), за да закръглите директно до цяло число, вместо да конвертирате в int32_t със знак и обратно.   -  person Peter Cordes    schedule 29.04.2018
comment
Не е ясно обаче защо бихте получили NaN. Cvt извън диапазона -2^31..2^31-1 ще ви даде 0x80000000 (т.е. -2^31, най-отрицателното цяло число), което Intel използва като целочислена неопределена стойност. Конвертирането обратно в double трябва да успее. Това може да са вашите 1.128e119 резултати? Очаква се да получите катастрофални грешки при анулиране, особено след като вашето 2Pi * recip_2Pi вероятно не е точно 1.0, защото грешките при закръгляването ще бъдат различни. (Всъщност 1.0- 2Pi*recip_2Pi с gcc -O0 е само ~2,1E-15, така че това е няколко порядъка по-малко от вашите входове.)   -  person Peter Cordes    schedule 29.04.2018

Отговори (2)


arrow_upward
7
arrow_downward

Ако искате някаква точност, простият алгоритъм е ужасно лош. За точен алгоритъм за намаляване на обхвата вижте напр. Ng et al., НАМАЛЯВАНЕ НА АРГУМЕНТИ ЗА ОГРОМНИ АРГУМЕНТИ: Добро до последния момент (вече достъпно чрез Wayback Machine: 2012-12-24)

person janneb    schedule 20.07.2012
comment
Трябва да прегледам този документ отново, за да се уверя, че съм го разбрал, но този метод не изисква ли повече битове, отколкото са налични в XMM регистрите? Този метод действително ли е работещ в SSE? - person BigA; 20.07.2012
comment
@BigA: Наистина, не можете да поставите всички битове в една XMM reg наведнъж. Ще трябва да направите нещо по-сложно. - person janneb; 20.07.2012
comment
Предимство на простия алгоритъм fmod 2pi, което не бива да се пренебрегва, е, че въведените по този начин грешки при закръгляване често ще противодействат на по-ранни грешки при закръгляване при изчисляването на ъгъла. Например резултатът от израза Math.Sin(x*(2*Math.Pi)) като цяло ще бъде по-точен, ако Sin изпълни редукция на аргумент mod Math.Pi, отколкото ако извърши редукция mod π [за някои стойности на x редукция mod π ще работи малко по-добре, но за някои работи много по-зле]. - person supercat; 04.06.2014

arrow_upward
1
arrow_downward

За големи аргументи алгоритъмът на Hayne-Panek обикновено е използвани. Документът на Hayne-Panek обаче е доста труден за четене и предлагам да погледнете Глава 11 в Наръчника за аритметика с плаваща запетая за по-достъпно обяснение.

person Marat Dukhan    schedule 20.07.2012