;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Pour archive : ça utilise une série de Taylor. Mais on va utiliser le sinus hard-codé en vrai.

section .data align=16
piX2 dq 6.283185307179586, 6.283185307179586
pi dq 3.141592653589793, 3.141592653589793
borneA dq -5.49778714378, -3.92699081699, -2.35619449019, -0.78539816339, 0.78539816339, 0.
a dq -6.283185307179586, -4.71238898038469, -3.141592653589793, -1.5707963267949, 0., 1.5707963267949, 3.141592653589793
  dq 4.71238898038469
sinA dq 0., 1., 0., -1., 0., 1., 0., -1.
cosA dq 1., 0., -1., 0., 1., 0., -1., 0.
diviseur dq 1., 1.
tailleDouble dq 8., 8.
section .text align=16	;Oui, une section .text. Admettons 

sinDVold:
	movapd xmm9, oword [piX2];On met dans un registre du processeur à virgule flottante 2 double : pi * 2
	movapd xmm10, [tailleDouble]
	xor r8, r8		;Une MàZ d'un registre qui va se souvenir de la parité du nombre d'éléments à sommer. Ah ben oui, on fait 2 sommes d'un coup, d'accord, mais si c'est impair ?
	mov rcx, rdi	;Nombre d'éléments à sommer mis dans RCX parce que c'est le registre compteur
	shr rcx, 1		;qu'on divise par 2, parce qu'on va les sommer par 2
	jnc .somme	;Et si la retenue (carry flag) est posée, la taille est impaire
	inc r8			;On se souvient que c'est impair pour plus tard
	test cl, cl;
	jz .dernier
.somme:
	movapd xmm0, oword [rsi];On met dans un registre du processeur à virgule flottante 2 double
	;On divise par 2 * pi puis multiplie la partie entière du quotient par 2 * pi
	movapd xmm1, xmm0
	divpd xmm1, xmm9
	cvtpd2dq xmm2, xmm1
	cvtdq2pd xmm1, xmm2
	mulpd xmm1, xmm9
	subpd xmm0, xmm1
	pxor xmm1, xmm1;On va descendre de 2*Pi tout ce qui est supérieur à 0
	cmppd xmm1, xmm0, 1
	andpd xmm1, xmm9;AND entre xmm1 et pi * 2
	subpd xmm0, xmm1
	;Déjà mieux : ]-2 * pi, 0]
	xor rax, rax
	xor rdx, rdx
	xor r10, r10
	xor r11, r11
	mov r9, 5
.chercheA:
	add r10, rdx
	add r11, rax
	movlpd xmm1, [r10+borneA]		; a
	movhpd xmm1, [r11+borneA]
	cmppd xmm1, xmm0, 1
;	andpd xmm1, xmm10;AND entre xmm1 et 8
;	cvtpd2dq xmm2, xmm1 
;	movq rax, xmm2
	movmskpd rax, xmm1
	test rax, rax
	jz .chargement
	mov edx, eax
	and edx, 1
	shl edx, 3
	and eax, 2
	shl eax, 2
;	shr rax, 32
	dec r9
	jnz .chercheA
	;boucler
.chargement:
	movlpd xmm3, [r10+a]		; a
	movhpd xmm3, [r11+a]
	movlpd xmm1, [r10+sinA]		;On stocke sin(a)
	movhpd xmm1, [r11+sinA]		;On stocke sin(a)
	movlpd xmm2, [r10+cosA]		;cos(a)
	movhpd xmm2, [r11+cosA]		;cos(a)
	movapd xmm6, [diviseur]		; 1
	movapd xmm7, xmm6		; 1
	movapd xmm8, xmm6		; 1
	subpd xmm0, xmm3		;xmm0 = x - a
	xor rdx, rdx
	inc rdx
	movapd xmm5, xmm0
	movapd xmm3, xmm1
.taylor:
	movapd xmm4, xmm2		;
	movapd xmm2, xmm1		;On échange cos(a) et sin(a)
	movapd xmm1, xmm4
	mulpd xmm4, xmm5		;xmm4 = cos(a)(x-a)
	divpd xmm4, xmm6		;xmm4 = cos(a)(x-a)/1
	test rdx, 2				;Comment dire : 2 tous les 4 : soustraction, les autres : addition
	jnz .sub
	addpd xmm3, xmm4		;xmm3 = sin(a) + cos(a)(x-a)
.taylor2:
	inc rdx
	mulpd xmm5, xmm0		;xmm5 = (x-a)*(x-a)
	addpd xmm7, xmm8
	mulpd xmm6, xmm7
	cmp rdx, 16
	jb .taylor

.stockage:
	movntpd oword [rsi], xmm3;On stocke le résultat dans l'accumulateur
	add rsi, 16		;On passe à la paire de double suivante
	dec rcx
	jnz .somme
;	loop .somme		;Et on boucle sur RCX

.dernier:
	;~ cmp r8, 0		;Si la taille des vecteurs est paire
	;~ jz .retour		;Merci, c'est fini. Sinon, c'est impair.
;	movsd xmm0, qword [rax]	;On met dans un registre du processeur à virgule flottante le dernier double
;	rsqrtsd xmm0, xmm0	;On fait une seule somme
;	movsd qword [rax], xmm0	;Et donc une seule sauvegarde. Et fin.
	
.retour:
	ret
.sub:
	subpd xmm3, xmm4		;xmm3 = sin(a) + cos(a)(x-a) - sin(a)(x-a)*(x-a)/2
	jmp .taylor2
