🚀 go-pugleaf

No nic, programik działa, interpoluje (zapewne nieoptymalnie, ale
cóż...). Kod w ASMie wrzucę, tak dla potomności, oraz może ktoś
skomentuje i poprawi moje wypociny :)

.586
.387

.MODEL FLAT, STDCALL

OPTION CASEMAP:NONE

INCLUDE c:\masm32\include\windows.inc;

.DATA
.CODE

DllEntry PROC ;hInstDLL:HINSTANCE, reason:DWORD, reserved1:DWORD
	mov eax, TRUE
	ret
DllEntry ENDP

;******************************************************************
;****** Procedura obliczajaca wspolczynnik beta dla j-ego elementu
;****** Wg wzoru:
;****** double beta (int j) //obliczanie wspólczynnika beta - wzór 26
;****** {
;******		return ((j* 2 * pi) / (n + 1));
;****** }

Beta PROC n:DWORD, j:DWORD, wynik:PTR
	pushad
	mov ebx, wynik		; adres wyniku

	finit
	fldpi				; PI -> st0
	fadd st(0), st(0)	; 2*PI -> st0
	fild j				; 2*PI -> st1, j -> st0
	fmulp st(1), st(0)	; 2*PI*j -> st0
	fld1				; 1 -> st0
	fild n				; n -> st0
	faddp st(1), st(0)	; n+1 -> st0
	fdivp st(1), st(0)	; 2*PI*j / n+1 -> st0

	fstp QWORD PTR [ebx]; zapisanie wyniku pod adres przekazanego
wskaznika

	popad
	ret					;koniec funkcji
Beta ENDP

;******************************************************************
;****** Procedura obliczajaca pierwsze 2 wyrazy ciagu Un
;****** Wg wzoru:
;****** double U (int k, double beta, Element *tab)  //obliczanie
wspólczynnika U - wzór 28
;****** {
;******		if ( k > n )
;******			return (0.0);
;******		else
;******			return (tab[k].f_x + (2.0) * ( U(k+1, beta, tab) * cosf
(beta) ) - U(k+2, beta, tab));
;****** }

;****** definicja elementu tablicy

;****** struct Element
;****** {
;****** 	double x;
;******		double f_x;
;****** 	double g_x;
;****** 	double blad;
;****** };

wspU PROC n:DWORD, k:DWORD, beta:PTR, tab:PTR, wynik:PTR
LOCAL t1:QWORD, t2:QWORD, t3:QWORD, t4:QWORD, temp:DWORD, tabptr:PTR,
twynik:QWORD
	pushad
	finit

	mov eax, k
	cmp eax, n				; if (k > n)

	ja Zwroc0				; return 0;

	mov ebx, tab			; zaladowanie adresu tablicy
	fld QWORD PTR [ebx+8]	; zaladowanie wartosci f(x)
	fstp t1					; zachowanie wartosci do zmiennej lokalnej, zwolnienie
stosu

	mov ebx, beta			; zaladowanie adresu zmiennej beta (liczonej przez
funkcje "Beta(j, n)"
	fld QWORD PTR [ebx]		; zaladowanie wartosci spod wskaznika beta
	fcos					; cosinus (beta)
	fadd st(0), st(0)		; 2 * cos (beta)
	fstp t2					; zachowanie powyzszej wartosci do lokalnej zmiennej

	mov eax, k				; zachowanie zmiennej k
	mov temp, eax
	inc temp				; temp = k+1
	mov esi, tab			; zachowanie adresu tablicy
	mov tabptr, esi
	add tabptr, 32			; przesuniecie sie do nastepnego elementu tablicy

	lea edi, twynik			; wyluskanie adresu zmiennej twynik

	INVOKE wspU, n, temp, beta, tabptr, edi				; wywolanie funkcji wspU
(n, k+1, beta, tab[k+1])

	lea edi, twynik			; zaladowanie adresu wyniku
	fld QWORD PTR [edi]		; odczytanie wyniku zwroconego przez wspU
	fstp t3					; zapisanie wyniku

	mov eax, k				; zachowanie zmiennej k
	mov temp, eax
	inc temp				; temp = k+1
	inc temp				; temp = k+2
	mov esi, tab			; zachowanie adresu tablicy
	mov tabptr, esi
	add tabptr, 64			; przesuniecie sie w tablicy o 2 elementy

	lea edi, twynik			; wyluskanie adresu zmiennej twynik

	INVOKE wspU, n, temp, beta, tabptr, edi				; wywolanie funkcji wspU
(n, k+2, beta, tab[k+2])

	lea edi, twynik			; zaladowanie adresu wyniku
	fld QWORD PTR [edi]		; odczytanie wyniku zwroconego przez wspU
	fstp t4					; zapisanie wyniku

	fld t3					; st1 = U(k+1)
	fld t2					; st0 = 2 * cos(beta)

	fmulp st(1), st(0)		; st0 = 2 * cos(beta) * U(k+1)
	fld t4					; st0 = U(k+2)
	fsubp st(1), st(0)		; st0 = 2 * cos(beta) * U(k+1) - U(k+1)

	mov ebx, tab
	fld QWORD PTR [ebx+8]	; st0 = f(x)
	faddp st(1), st(0)		; st0 = f(x) + 2 * cos(beta) * U(k+1)

	mov esi, wynik			; wczytanie adresu zmiennej wynikowej
	fstp QWORD PTR [esi]	; zapisanie wyniku

	jmp Koniec

Zwroc0:
	fldz					; ladujemy 0 do st0

	mov esi, wynik			; wczytanie adresu zmiennej wynikowej
	fstp QWORD PTR [esi]	; zapisanie zera pod adres wynikowy

Koniec:
	popad
	ret ;koniec funkcji

wspU ENDP


;******************************************************************
;****** Procedura obliczajaca pierwsze 2 wyrazy ciagu Un
;****** Wg wzoru:
;****** double U (int k, double beta, Element *tab)  //obliczanie
wspólczynnika U - wzór 28
;****** {
;******		if ( k > n )
;******			return (0.0);
;******		else
;******			return (tab[k].f_x + (2.0) * ( U(k+1, beta, tab) * cosf
(beta) ) - U(k+2, beta, tab));
;****** }

;****** definicja elementu tablicy

;****** struct Element
;****** {
;****** 	double x;
;******		double f_x;
;****** 	double g_x;
;****** 	double blad;
;****** };

wspU_2 PROC n:DWORD, k:DWORD, beta:QWORD, tab:PTR, wynik:PTR
LOCAL Uk1:QWORD, Uk2:QWORD, tabT:PTR
	pushad
	finit

	mov eax, k
	cmp eax, n				; if (k > n)

	ja Zwroc0				; return 0

	mov ebx, k
	inc ebx					; ebx = k+1
	mov ecx, tab
	add ecx, 32
	mov tabT, ecx			; tabT = tab + 32 (przejscie do kolejnego elementu
tablicy)

	INVOKE wspU_2, n, ebx, beta, tabT, wynik

	mov ebx, wynik
	fld QWORD PTR [ebx]
	fstp Uk1				; zapisanie wyniku do zmiennej lokalnej

	mov ebx, k
	add ebx, 2				; ebx = k+2
	mov ecx, tab
	add ecx, 64
	mov tabT, ecx			; tabT = tab + 64 (przejscie o 2 elementy tablicy)

	INVOKE wspU_2, n, ebx, beta, tabT, wynik

	mov ebx, wynik
	fld QWORD PTR [ebx]
	fstp Uk2				; zapisanie wyniku do zmiennej lokalnej

	fld beta
	fcos					; cosinus (beta)
	fadd st(0), st(0)		; 2 * cos (beta)

	fld Uk1					; zaladowanie wartosci U(k+1)
	fmulp st(1), st(0)		; 2 * cos (beta) * U(k+1) -> st0

	fld Uk2
	fsubp st(1), st(0)		; 2 * cos (beta) * U(k+1) - U(k+2) -> st0

	mov esi, tab
	fld QWORD PTR [esi+8]
	faddp st(1), st(0)		; f(xk) + 2 * cos (beta) * U(k+1) - U(k+2) -> st0

	mov esi, wynik
	fstp QWORD PTR [esi]	; zapisanie wyniku

	jmp Koniec

Zwroc0:
	mov esi, wynik
	fldz
	fstp QWORD PTR [esi]	; zwrocenie zera

Koniec:
	popad
	ret;

wspU_2 ENDP


;******************************************************************
;****** Procedura liczaca wspolczynnik A
;****** double wspA (int j, Element *tab) //obliczanie wspolczynnika A
- wzór 27
;****** {
;****** 	double beta_j = beta (j);
;****** 	return ( (2.0/(n+1)) * (tab[0].f_x + U(1, beta_j, tab) * cosf
(beta_j) - U(2, beta_j, tab)) );
;****** }
wspA PROC n:DWORD, j:DWORD, tab:PTR, wynik:PTR
LOCAL Tbeta:QWORD, U2:QWORD, temp:DWORD
	pushad
	finit

	lea ebx, Tbeta			; wyluskanie adresu zmiennej Tbeta (bedzie w niej
przechowywany wynik funkcji Beta())
	INVOKE Beta, n, j, ebx

	mov esi, tab
	add esi, 64

	INVOKE wspU_2, n, 2, Tbeta, esi, wynik

	mov esi, wynik
	fld QWORD PTR [esi]		; zapamietanie wyniku
	fstp U2

	mov esi, tab
	add esi, 32

	INVOKE wspU_2, n, 1, Tbeta, esi, wynik

	mov esi, wynik
	fld QWORD PTR [esi]		; U(1) -> st0

	fld Tbeta
	fcos					; cos(beta) -> st0

	fmulp st(1), st(0)		; cos(beta) * U1 -> st0
	fld U2					; U(2) -> st0
	fsubp st(1), st(0)		; cos(beta) * U1 - U2 -> st0

	mov esi, tab
	fld QWORD PTR [esi+8]	; f(x0) -> st0
	faddp st(1), st(0)		; f(x0) + cos(beta) * U1 - U2 -> st(0)
	fadd st(0), st(0)		; 2 * st0

	mov ebx, n
	inc ebx					; n+1 -> ebx
	mov temp, ebx

	fild temp				; n+1 -> st0

	fdivp st(1), st(0)		; [2 * (x0) + cos(beta) * U1 - U2 -> st(0)] / (n
+1) -> st0

	mov esi, wynik
	fstp QWORD PTR [esi]	; zwrocenie wyniku

	popad

	ret
wspA ENDP

;******************************************************************
;****** Procedura liczaca wspolczynnik B
;****** double B (int j, Element* tab) //obliczanie współczynnika B -
wzór 27
;****** {
;****** 	double beta_j = beta (j);
;****** 	return (2 * ( U(n, 1, beta_j, tab) * sinf(beta_j) ) ) / (n +
1) ;
;****** }
wspB PROC n:DWORD, j:DWORD, tab:PTR, wynik:PTR
LOCAL Tbeta:QWORD, temp:DWORD
	pushad
	finit

	lea esi, Tbeta			; wyluskanie adresu zmiennej Tbeta (bedzie w niej
przechowywany wynik funkcji Beta())
	INVOKE Beta, n, j, esi

	mov esi, tab
	add esi, 32

	INVOKE wspU_2, n, 1, Tbeta, esi, wynik

	mov esi, wynik
	fld QWORD PTR [esi]		; U(1) -> st0

	fld Tbeta
	fsin					; sin(beta) -> st0

	fmulp st(1), st(0)		; U(1) * sin(beta) -> st0
	fadd st(0), st(0)		; 2 * U(1) * sin(beta) -> st0

	mov ebx, n
	inc ebx					; n+1 -> ebx
	mov temp, ebx

	fild temp				; n+1 -> st0

	fdivp st(1), st(0)		; [2 * sin(beta) * U1] / (n+1) -> st0

	mov esi, wynik
	fstp QWORD PTR [esi]	; zwrocenie wyniku

	popad
	ret
wspB ENDP

;******************************************************************
;****** Procedura obliczajaca wartosc funkcji interpolujacej w danym
punkcie.
;****** wg wzoru:
;****** t(x) =	A(0)/2 +													// lokalnie jako wA02
;******			[ delta * A(m+1) * cos ((m+1)*x) ] / 2 +					// lokalnie
jako wA
;******			SUMA (od j=1 do m) [ A(j) * cos (j*x) + B(j) * sin (j*x) ]
;****** gdzie:
;******			dla n parzystych:		delta = 0, m =	n	/ 2
;******			dla n nieparzystych:	delta = 1, m =(n-1)	/ 2
;******
;****** funkcja bedzie wywolywana osobno dla kazdego punktu
;****** interpolacji, stad 2 wskazniki na tablice tab[]
;****** jeden wskaznik na poczatek tablicy oraz jeden wskaznik
;****** na aktualnie liczony element

Interpoluj PROC n:DWORD, tab:PTR, AKT_tab:PTR
LOCAL m:DWORD, j:DWORD, wA:QWORD, wA02:QWORD, wAsuma:QWORD,
wBsuma:QWORD, SUMA:QWORD
	pushad
	finit

	mov eax, n				; zaladowanie n
	shl eax, 31				; wyluskanie najmlodszego bitu (1 - nieparzysta)

	cmp eax, 0				; sprawdzenie parzystosci liczby
	je Parzysta

Nieparzysta:				; tutaj bedziemy liczyc: [ A(m+1) * cos ((m+1)*x) ] /
2
	mov ebx, n
	dec ebx					; n - 1
	shr ebx, 1				; (n - 1) / 2
	mov m, ebx				; m = (n - 1) / 2

	lea esi, wA
	mov ebx, m
	inc ebx
	mov j, ebx				; j = m+1

	INVOKE wspA, n, j, tab, esi		; obliczamy A(m+1), wynik -> wA

	fld wA					; A(m+1) -> st0
	fild j					; m+1 -> st0
	fld QWORD PTR AKT_tab	; zaladowanie wartosci x dla aktualnie
liczonego wezla
	fmulp st(1), st(0)		; [ A(m+1) * x ] -> st0
	fcos					; cos ((m+1)*x) -> st0
	fmulp st(1), st(0)		; [ A(m+1) * cos ((m+1)*x) ] -> st0
	mov j, 2
	fild j					; zaladowanie 2 do st0
	fdivp st(1), st(0)		; [ A(m+1) * cos ((m+1)*x) ] / 2 -> st0
	fstp wA					; zapisanie powyzszego wyniku do wA

	jmp M_obliczone

Parzysta:					; ustawienie odpowiedniego m oraz wA = 0
	mov ebx, n
	shr ebx, 1				; n / 2
	mov m, ebx				; m = n / 2
	fldz
	fstp wA

M_obliczone:				; obliczenie A(0) / 2 i petla liczaca sume

	lea esi, wA02			; zaladowanie adresu zmiennej do esi

	INVOKE wspA, n, 0, tab, esi		; obliczamy A(0), wynik -> wA02

	fld wA02				; A(0) -> st0
	mov j, 2
	fild j					; zaladowanie 2 do st0
	fdivp st(1), st(0)		; A(0) / 2 -> st0
	fstp wA02				; zaladowanie st0 do zmiennej wA02

;* petla obliczajaca sume
;** SUMA (od j=1 do m) [ A(j) * cos (j*x) + B(j) * sin (j*x) ]

	mov ecx, 1								; for (j = 1...)
	fldz
	fstp SUMA								; SUMA = 0
Petla:
		mov j, ecx
		cmp ecx, m
		ja Endfor							; for (j = 1; j < m...)

		lea esi, wAsuma
		INVOKE wspA, n, j, tab, esi		; wAsuma = A(j)

		lea esi, wBsuma
		INVOKE wspB, n, j, tab, esi		; wBsuma = B(j)

		mov esi, AKT_tab
		fld QWORD PTR [esi]					; tab[j].x -> st0
		fild j								; j -> st0
		fmulp st(1), st(0)					; j * x -> st0
		fcos								; cos (j * x) -> st0
		fld wAsuma
		fmulp st(1), st(0)					; A(j) * cos (j * x) -> st0

		mov esi, AKT_tab
		fld QWORD PTR [esi]					; tab[j].x -> st0
		fild j								; j -> st0
		fmulp st(1), st(0)					; j * x -> st0
		fsin								; sin (j * x) -> st0
		fld wBsuma
		fmulp st(1), st(0)					; B(j) * cos (j * x) -> st0

		faddp st(1), st(0)					; [ A(j) * cos (j*x) + B(j) * sin (j*x) ] ->
st0

		fld SUMA
		faddp st(1), st(0)
		fstp SUMA							; SUMA += [ A(j) * cos (j*x) + B(j) * sin (j*x) ]

		inc ecx								; for (j = 1, j < m, j++)
		jmp Petla
Endfor:

	fld SUMA
	fld wA02
	fld wA
	faddp st(1), st(0)
	faddp st(1), st(0)						; SUMA += A(0)/2 + [ delta * A(m+1) * cos ((m
+1)*x) ] / 2

	mov esi, AKT_tab						; zaladowanie wskaznika na aktualny element
tablicy
	add esi, 16								; przesuniecie sie na element g_x struktury
	fstp QWORD PTR [esi]					; zapisanie wartosci funkcji interpolacyjnej
dla danego elementu

Liczenie_bledu:
	mov esi, AKT_tab
	fld QWORD PTR [esi+8]					; f(x) -> st0 __zaladowanie wartosci
oryginalnej funkcji
	fld QWORD PTR [esi+16]					; g(x) -> st0 __zaladowanie wartosci
funkcji interpolacyjnej

	fsubp st(1), st(0)						; f(x) - g(x) -> st(0)

	fldz
	fcomip st(0), st(1)						; porownanie st0 i st1
	jbe bladdodatni

bladujemny:
	fldz
	fld1
	fsubp st(1), st(0)						; -1 -> st0
	fmulp st(1), st(0)						; zmiana znaku bledu

bladdodatni:
	fstp QWORD PTR [esi+24]					; zapisanie wartosci bledu
Koniec:
	popad
	ret
Interpoluj ENDP

END DllEntry