[ Branimir Maksimovic @ 30.08.2017. 04:46 ] @
Uradio sam ovo malo pa reko da podelim ako nekom treba, koristio sam cephes math lib da pretabam u avx verziju. Inace ne handlujem greske za nan i inf,ali to sad pa i nije bitno.
asembler je fasm za Linux. jedino je sin avx2 pa to mozete da vratite na avx1 lako gledajuci ostale f-je.

Code:

format elf64
define AVX2 1

macro polevl p1,p2,p3 {
        local .L0
        ; %1 = ymm value,%2 = addr coef array , %3 = N
        mov rbx,p2
        lea rcx,[rbx+8*p3]
        vbroadcastsd ymm15,[rbx]
.L0:
        add rbx,8
if AVX2
        vbroadcastsd ymm14,[rbx]
        vfmadd213pd ymm15,p1,ymm14
else
        vmulpd ymm15,ymm15,p1
        vbroadcastsd ymm14,[rbx]
        vaddpd ymm15,ymm15,ymm14
end if
        cmp rbx,rcx
        jl .L0
}

macro p1evl p1,p2,p3 {
        local .L0
        ; %1 = ymm value,%2 = addr coef array , %3 = N
        mov rbx,p2
        lea rcx,[rbx+8*(p3-1)]
        vbroadcastsd ymm15,[rbx]
        vaddpd ymm15,ymm15,p1
.L0:
        add rbx,8
if AVX2
        vbroadcastsd ymm14,[rbx]
        vfmadd213pd ymm15,p1,ymm14
else
        vmulpd ymm15,ymm15,p1
        vbroadcastsd ymm14,[rbx]
        vaddpd ymm15,ymm15,ymm14
end if
        cmp rbx,rcx
        jl .L0
}

section '.text' executable

public tangent
public sine
public cosine
public sin_pd
public cos_pd
public tan_pd
public x87sincos
public sincos_pd

tangent:
        sub rsp,8
        vmovsd [rsp],xmm0
        fld qword[rsp]
        fptan
        fstp st0
        fstp qword[rsp]
        vmovsd xmm0,[rsp]
        add rsp,8
        ret

sine:
        sub rsp,8
        vmovsd [rsp],xmm0
        fld qword[rsp]
        fsin
        fstp qword[rsp]
        vmovsd xmm0,[rsp]
        add rsp,8
        ret
cosine:
        sub rsp,8
        vmovsd [rsp],xmm0
        fld qword[rsp]
        fcos
        fstp qword[rsp]
        vmovsd xmm0,[rsp]
        add rsp,8
        ret
x87sincos:
        sub rsp,8
        vmovsd [rsp],xmm0
        fld qword[rsp]
        fsincos
        fstp qword[rsi]
        fstp qword[rdi]
        add rsp,8
        ret
; ymm0 -> x
sin_pd:
        push rbx
        vmovapd ymm1,ymm0
        vmovapd ymm2,[sign_mask]
        vandnpd ymm0,ymm2,ymm0 ; make positive x
        vandpd ymm1,ymm1,[sign_mask] ; save sign bit

        vbroadcastsd ymm2,[O4PI]
        vmulpd ymm3,ymm0,ymm2 ; 
        vroundpd ymm3,ymm3,3 ; truncate y

        vmulpd ymm2,ymm3,[OD16]
        vroundpd ymm2,ymm2,3

;       vmulpd ymm2,ymm2,[S16]
;       vsubpd ymm2,ymm3,ymm2

        vfnmadd132pd ymm2,ymm3,[S16]

        vcvttpd2dq xmm2,ymm2 ; j

        vpand xmm4,xmm2,[mask_1]
        vpaddd xmm2,xmm2,xmm4 ; j += 1
        vcvtdq2pd ymm4,xmm4
        vaddpd ymm3,ymm3,ymm4 ; y += 1.0

        vpand xmm4,xmm2,[mask_4]
        vpslld xmm4,xmm4,29 ; move mask to highest position
if AVX2 ; just example, too lazy to repeat for other functions ;)
        vpmovzxdq ymm4,xmm4
        vpsllq ymm4,ymm4,32
else
        vpmovzxdq xmm5,xmm4
        vpsllq xmm5,xmm5,32
        vpsrldq xmm4,xmm4,8
        vpmovzxdq xmm6,xmm4
        vpsllq xmm6,xmm6,32
        vmovapd xmm4,xmm5
        vinsertf128 ymm4,ymm4,xmm6,1
end if
        vxorpd ymm1,ymm1,ymm4 ; invert sign

        vpand xmm4,xmm2,[mask_3]
        vpcmpeqd xmm4,xmm4,[mask_0] 
if AVX2
        vpmovsxdq ymm4,xmm4
else
        vpmovsxdq xmm5,xmm4
        vpsrldq xmm4,xmm4,8
        vpmovsxdq xmm6,xmm4
        vmovapd xmm4,xmm5
        vinsertf128 ymm4,ymm4,xmm6,1 ; selection mask 
end if

; Extended precision modular arithmetic
;       vmulpd ymm5,ymm3,[DP1]
;       vmulpd ymm6,ymm3,[DP2]
;       vmulpd ymm7,ymm3,[DP3]
;       vsubpd ymm0,ymm0,ymm5
;       vsubpd ymm0,ymm0,ymm6
;       vsubpd ymm0,ymm0,ymm7

        vfnmadd231pd ymm0,ymm3,[DP1]
        vfnmadd231pd ymm0,ymm3,[DP2]
        vfnmadd231pd ymm0,ymm3,[DP3]

        vmulpd ymm5,ymm0,ymm0 ; x^2

        ; first
        polevl ymm5,sincof,5
        vmulpd ymm15,ymm15,ymm5
;       vmulpd ymm15,ymm15,ymm0
;       vaddpd ymm6,ymm15,ymm0 ; y1

        vfmadd132pd ymm15,ymm0,ymm0
        vmovapd ymm6,ymm15

        ; second
        polevl ymm5,coscof,5
        vmulpd ymm15,ymm15,ymm5
        vmulpd ymm15,ymm15,ymm5
;       vmulpd ymm7,ymm5,[OD2]
;       vsubpd ymm7,ymm15,ymm7

        vfnmadd231pd ymm15,ymm5,[OD2]
        vmovapd ymm7,ymm15

        vaddpd ymm7,ymm7,[one]; y2

        ; combine
        vandpd ymm6,ymm4,ymm6
        vandnpd ymm7,ymm4,ymm7
        vaddpd ymm0,ymm6,ymm7

        vxorpd ymm0,ymm0,ymm1
        pop rbx
        ret
; ymm0 -> x
cos_pd:
        push rbx
        vmovapd ymm1,ymm0
        vmovapd ymm2,[sign_mask]
        vandnpd ymm0,ymm2,ymm0 ; make positive x
        vandpd ymm1,ymm1,[sign_mask] ; save sign bit

        vbroadcastsd ymm2,[O4PI]
        vmulpd ymm3,ymm0,ymm2 ; 
        vroundpd ymm3,ymm3,3 ; truncate y

        vmulpd ymm2,ymm3,[OD16]
        vroundpd ymm2,ymm2,3
        vmulpd ymm2,ymm2,[S16]
        vsubpd ymm2,ymm3,ymm2
        vcvttpd2dq xmm2,ymm2 ; j

        vpand xmm4,xmm2,[mask_1]
        vpaddd xmm2,xmm2,xmm4 ; j += 1
        vcvtdq2pd ymm4,xmm4
        vaddpd ymm3,ymm3,ymm4 ; y += 1.0

        vpand xmm4,xmm2,[mask_4]
        vpslld xmm4,xmm4,29 ; move mask to highest position
        vpmovzxdq xmm5,xmm4
        vpsllq xmm5,xmm5,32
        vpsrldq xmm4,xmm4,8
        vpmovzxdq xmm6,xmm4
        vpsllq xmm6,xmm6,32
        vmovapd xmm4,xmm5
        vinsertf128 ymm4,ymm4,xmm6,1
        vxorpd ymm1,ymm1,ymm4 ; invert sign

        vpand xmm4,xmm2,[mask_3]
        vpcmpeqd xmm4,xmm4,[mask_0] 
        vpmovsxdq xmm5,xmm4
        vpsrldq xmm4,xmm4,8
        vpmovsxdq xmm6,xmm4
        vmovapd xmm4,xmm5
        vinsertf128 ymm4,ymm4,xmm6,1 ; selection mask 

        ; Extended precision modular arithmetic
        vmulpd ymm5,ymm3,[DP1]
        vmulpd ymm6,ymm3,[DP2]
        vmulpd ymm7,ymm3,[DP3]
        vsubpd ymm0,ymm0,ymm5
        vsubpd ymm0,ymm0,ymm6
        vsubpd ymm0,ymm0,ymm7

        vmulpd ymm5,ymm0,ymm0 ; x^2

        ; first
        polevl ymm5,sincof,5
        vmulpd ymm15,ymm15,ymm5
        vmulpd ymm15,ymm15,ymm0
        vaddpd ymm6,ymm15,ymm0 ; y1

        ; second
        polevl ymm5,coscof,5
        vmulpd ymm15,ymm15,ymm5
        vmulpd ymm15,ymm15,ymm5
        vmulpd ymm7,ymm5,[OD2]
        vsubpd ymm7,ymm15,ymm7
        vaddpd ymm7,ymm7,[one]; y2

        ; combine
        vandnpd ymm6,ymm4,ymm6
        vandpd ymm7,ymm4,ymm7
        vaddpd ymm0,ymm6,ymm7

        vxorpd ymm0,ymm0,ymm1
        pop rbx
        ret

; ymm0 -> x
sincos_pd:
        push rbx
        vmovapd ymm1,ymm0
        vmovapd ymm2,[sign_mask]
        vandnpd ymm0,ymm2,ymm0 ; make positive x
        vandpd ymm1,ymm1,[sign_mask] ; save sign bit

        vbroadcastsd ymm2,[O4PI]
        vmulpd ymm3,ymm0,ymm2 ; 
        vroundpd ymm3,ymm3,3 ; truncate y

        vmulpd ymm2,ymm3,[OD16]
        vroundpd ymm2,ymm2,3
        vmulpd ymm2,ymm2,[S16]
        vsubpd ymm2,ymm3,ymm2
        vcvttpd2dq xmm2,ymm2 ; j

        vpand xmm4,xmm2,[mask_1]
        vpaddd xmm2,xmm2,xmm4 ; j += 1
        vcvtdq2pd ymm4,xmm4
        vaddpd ymm3,ymm3,ymm4 ; y += 1.0

        vpand xmm4,xmm2,[mask_4]
        vpslld xmm4,xmm4,29 ; move mask to highest position
        vpmovzxdq xmm5,xmm4
        vpsllq xmm5,xmm5,32
        vpsrldq xmm4,xmm4,8
        vpmovzxdq xmm6,xmm4
        vpsllq xmm6,xmm6,32
        vmovapd xmm4,xmm5
        vinsertf128 ymm4,ymm4,xmm6,1
        vxorpd ymm1,ymm1,ymm4 ; invert sign

        vpand xmm4,xmm2,[mask_3]
        vpcmpeqd xmm4,xmm4,[mask_0] 
        vpmovsxdq xmm5,xmm4
        vpsrldq xmm4,xmm4,8
        vpmovsxdq xmm6,xmm4
        vmovapd xmm4,xmm5
        vinsertf128 ymm4,ymm4,xmm6,1 ; selection mask 

        ; Extended precision modular arithmetic
        vmulpd ymm5,ymm3,[DP1]
        vmulpd ymm6,ymm3,[DP2]
        vmulpd ymm7,ymm3,[DP3]
        vsubpd ymm0,ymm0,ymm5
        vsubpd ymm0,ymm0,ymm6
        vsubpd ymm0,ymm0,ymm7

        vmulpd ymm5,ymm0,ymm0 ; x^2

        ; first
        polevl ymm5,sincof,5
        vmulpd ymm15,ymm15,ymm5
        vmulpd ymm15,ymm15,ymm0
        vaddpd ymm6,ymm15,ymm0 ; y1

        ; second
        polevl ymm5,coscof,5
        vmulpd ymm15,ymm15,ymm5
        vmulpd ymm15,ymm15,ymm5
        vmulpd ymm7,ymm5,[OD2]
        vsubpd ymm7,ymm15,ymm7
        vaddpd ymm7,ymm7,[one]; y2

        ; combine sin
        vandpd ymm8,ymm4,ymm6
        vandnpd ymm9,ymm4,ymm7
        vaddpd ymm0,ymm8,ymm9

        vxorpd ymm0,ymm0,ymm1
        vmovupd [rdi],ymm0

        ; combine cos
        vandnpd ymm6,ymm4,ymm6
        vandpd ymm7,ymm4,ymm7
        vaddpd ymm0,ymm6,ymm7

        vxorpd ymm0,ymm0,ymm1
        vmovupd [rsi],ymm0

        pop rbx
        ret
tan_pd:
        push rbx
        vmovapd ymm1,ymm0
        vmovapd ymm2,[sign_mask]
        vandnpd ymm0,ymm2,ymm0 ; make positive x
        vandpd ymm1,ymm1,[sign_mask] ; save sign bit

        vbroadcastsd ymm2,[O4PI]
        vmulpd ymm3,ymm0,ymm2 ; 
        vroundpd ymm3,ymm3,3 ; truncate y

        vmulpd ymm2,ymm3,[OD16]
        vroundpd ymm2,ymm2,3
        vmulpd ymm2,ymm2,[S16]
        vsubpd ymm2,ymm3,ymm2
        vcvttpd2dq xmm2,ymm2 ; j

        vpand xmm4,xmm2,[mask_1]
        vpaddd xmm2,xmm2,xmm4 ; j += 1
        vcvtdq2pd ymm4,xmm4
        vaddpd ymm3,ymm3,ymm4 ; y += 1.0

        vpand xmm4,xmm2,[mask_2]
        vpcmpeqd xmm4,xmm4,[mask_0] 
        vpmovsxdq xmm5,xmm4
        vpsrldq xmm4,xmm4,8
        vpmovsxdq xmm6,xmm4
        vmovapd xmm4,xmm5
        vinsertf128 ymm4,ymm4,xmm6,1 ; selection mask 2

        ; Extended precision modular arithmetic
        vmulpd ymm5,ymm3,[DP1]
        vmulpd ymm6,ymm3,[DP2]
        vmulpd ymm7,ymm3,[DP3]
        vsubpd ymm0,ymm0,ymm5
        vsubpd ymm0,ymm0,ymm6
        vsubpd ymm0,ymm0,ymm7

        vmulpd ymm5,ymm0,ymm0 ; x^2

        vcmpnlepd ymm6,ymm5,[prec] ; selection mask 1

        ;calculate polynom
        polevl ymm5,P,2
        vmovapd ymm13,ymm15
        p1evl ymm5,Q,4
        vdivpd ymm13,ymm13,ymm15
        vmulpd ymm13,ymm13,ymm5
        vmulpd ymm13,ymm13,ymm0
        vaddpd ymm13,ymm13,ymm0

        vandpd ymm13,ymm6,ymm13
        vandnpd ymm0,ymm6,ymm0 ; select according to mask 1
        vaddpd ymm0,ymm13,ymm0

        vmovapd ymm6,[mone]
        vdivpd ymm7,ymm6,ymm0

        vandpd ymm0,ymm4,ymm0
        vandnpd ymm7,ymm4,ymm7 ; select according to mask 2
        vaddpd ymm0,ymm0,ymm7

        vxorpd ymm0,ymm0,ymm1 ; invert sign
        pop rbx
        ret

section '.data' writeable align 32
sincof dq \
        1.58962301576546568060E-10,\
        -2.50507477628578072866E-8,\
        2.75573136213857245213E-6,\
        -1.98412698295895385996E-4,\
        8.33333333332211858878E-3,\
        -1.66666666666666307295E-1
coscof dq \
        -1.13585365213876817300E-11,\
        2.08757008419747316778E-9,\
        -2.75573141792967388112E-7,\
        2.48015872888517045348E-5,\
        -1.38888888888730564116E-3,\
        4.16666666666665929218E-2
P dq \
-1.30936939181383777646E4, \
 1.15351664838587416140E6, \
 -1.79565251976484877988E7 
Q dq \
 1.36812963470692954678E4, \
 -1.32089234440210967447E6,\
  2.50083801823357915839E7,\
  -5.38695755929454629881E7 
O4PI dq 1.273239544735162
align 32
DP1: times 4 dq   7.85398125648498535156E-1;
DP2: times 4 dq   3.77489470793079817668E-8;
DP3: times 4 dq   2.69515142907905952645E-15;
sign_mask: times 4 dq 0x8000000000000000
prec: times 4 dq 1.0e-14
one: times 4 dq 1.0
mone: times 4 dq -1.0
OD16: times 4 dq 0.0625
S16: times 4 dq 16.0
OD2: times 4 dq 0.5
mask_0: times 4 dd 0
mask_1: times 4 dd 1
mask_2: times 4 dd 2
mask_3: times 4 dd 3
mask_4: times 4 dd 4


Evo kako se koristi:

Code:

#include <math.h>
#include <time.h>
#include <stdio.h>
#include <immintrin.h>
extern "C" double tangent(double x);
extern "C" double sine(double x);
extern "C" double cosine(double x);
extern "C" void x87sincos(double x,double* rsin,double* rcos);
extern "C" double md_sin(double x);
extern "C" double md_cos(double x);
extern "C" double md_tan(double x);
extern "C" __m256d sin_pd(__m256d x);
extern "C" __m256d cos_pd(__m256d x);
extern "C" void sincos_pd(__m256d x,__m256d* rsin,__m256d* rcos);
extern "C" __m256d tan_pd(__m256d x);
typedef double(*fp_t)(double);
typedef void(*fpb_t)(double,double*,double*);
typedef __m256d(*fpv_t)(__m256d);
typedef void(*fpvb_t)(__m256d,__m256d*,__m256d*);

#define NITER 10000000
void test(fp_t f, const char* str)
{
  volatile double x = -1000.0, y;
  clock_t begin = clock();
  for(int i = 0;i<NITER;++i)
  {
    y = f(x);
    x+=0.001;
  }
  clock_t end = clock();
  double diff = double(end-begin)/CLOCKS_PER_SEC;
  printf("%.15f %s time %f\n",y,str,diff);
}
void test(fpb_t f, const char* str)
{
  double x = -1000.0, s,c;
  clock_t begin = clock();
  for(int i = 0;i<NITER;++i)
  {
    f(x,&s,&c);
    x+=0.001;
  }
  clock_t end = clock();
  double diff = double(end-begin)/CLOCKS_PER_SEC;
  printf("%.15f %.15f %s time %f\n",s,c,str,diff);
}
void test(fpv_t f, const char* str)
{
  volatile __m256d x = _mm256_set1_pd(-1000.0),one = _mm256_set1_pd(0.001),y;
  double res[4];
  clock_t begin = clock();
  for(int i = 0;i<NITER;++i)
  {
    y = f(x);
    x = _mm256_add_pd(x,one);
  }
  clock_t end = clock();
  double diff = double(end-begin)/CLOCKS_PER_SEC;
  _mm256_storeu_pd(res,y);
  printf("%.15f %s time %f\n",res[2],str,diff);
}
void test(fpvb_t f, const char* str)
{
  __m256d x = _mm256_set1_pd(-1000.0),one = _mm256_set1_pd(0.001);
  double ress[4],resc[4];
  clock_t begin = clock();
  for(int i = 0;i<NITER;++i)
  {
    f(x,(__m256d*)ress,(__m256d*)resc);
    x = _mm256_add_pd(x,one);
  }
  clock_t end = clock();
  double diff = double(end-begin)/CLOCKS_PER_SEC;
  printf("%.15f %.15f %s time %f\n",
  ress[0],resc[3],str,diff);
}

int main()
{
  test(tangent,"x87 tangent");
  test(tangent,"x87 tangent");
  test(tan,"clib tan");
  test(md_tan,"cephes tan");
  test(tan_pd,"tan_pd");
  test(sine,"x87 sine");
  test(sin,"clib sin");
  test(md_sin,"cephes sin");
  test(sin_pd,"sin_pd");
  test(cosine,"x87 cosine");
  test(cos,"clib cos");
  test(md_cos,"cephes cos");
  test(cos_pd,"cos_pd");
  test(x87sincos,"x87 sincos");
  test(sincos,"clib sincos");
  test(sincos_pd,"sincos_pd");
}


[ Living Light @ 30.08.2017. 05:17 ] @
Branimire, Matematicaru, ----Svaka Ti Cast !!!

Ja ne "hendlujem" sa taj matiš, aLi, ali...

...kad citam tvoje postove, dodje mi DRAGO,

da je neko tako lepo i duboko u tim Integralima, Diferencialima, Limesima i jos +100Q x e na X !!!

EEee jesi MAHER ! ...svaka ti dala,