flib.h

FMUL in assembly. Enjoy!

#pragma once
#include <stdio.h>

int asmMUL(double src1, double src2, double* dest)	// 0 -succes; 1- overflow ; 2 -underflow
{
	double rez=0;
	int valret=-1;
	__asm
	{	
		mov	edx,dword ptr [src1+4]	; edx <- xu
		mov	eax,dword ptr [src1]	; eax <- xl
		mov	edi,dword ptr [src2+4]	; edi <- yu
		mov	esi,dword ptr [src2]	; esi <- yl

	; Step 1: exponent (and sign-bit) of product -> ecx

		mov	ecx,edx		; exponent ...
		and	ecx,7ff00000h	; ... of x ...
		jz	ret0		; (return 0 if x == 0)
		shr	ecx,20		; in cx
		mov	ebx,edi		; exponent ...
		and	ebx,7ff00000h	; ... of y ...
		jz	ret0		; (return 0 if y == 0)
		shr	ebx,20		; in bx
		add	cx,bx
		sub	cx,3feh	; exponent of product (exponent offset = 0x3fe)
		jmp	continue

	ret0:
		jmp	retz

	continue:
		mov	ebx,edi		; yu
		xor	ebx,edx		; 'xor'ed with xu
		and	ebx,80000000h	; sign-bit of product ...
		or	ecx,ebx		; added to ecx

	; Step 2: product of mantissae -> (edx,eax)

	; 2a: (xu,xl) -> (edx,eax)

		shld	edx,eax,11	; extract mantissa of x ...
		shl	eax,11
		or	edx,80000000h	; ... and restore leading bit

	; 2b: (yu,yl) -> (edi,esi)

		shld	edi,esi,11	; extract mantissa of y ...
		shl	esi,11
		or	edi,80000000h	; ... and restore leading bit

	; 2c: xl * yu (upper dword) -> esi

		mov	ebx,edx		; save xu in ebx
		mul	edi		; xl * yu
		mov	eax,esi		; eax <- yl
		mov	esi,edx		; save upper word of xl  yu in esi

	; 2c: yl * xu (upper dword) -> ebx

		mul	ebx		; yl * xu
		mov	eax,ebx		; eax <- xu
		mov	ebx,edx		; save upper dword of yl  xu in ebx

	; 2d: xu * yu

		mul	edi		; xu * yu -> (edx,eax)

	; 2e: add the parts and store in (edx,eax)

		add	eax,esi
		adc	edx,0
		add	eax,ebx
		adc	edx,0

	; Step 3: normalise product

		test	edx,80000000h	; is mantissa too small? (can only be by 1 bit)
		jnz	mul4
		shld	edx,eax,1	; yes; multiply mantissa ...
		shl	eax,1		; ... by 2 ...
		dec	cx		; ... and decrement exponent

	; Step 4: shift and round
	mul4:
		mov	esi,ebx		; save xu for sign-bit
		shrd	eax,edx,11	; shift mantissa of result ...
		jc	roundup		; (if last bit shifted out is 1, round up)
		shr	edx,11		; ... to proper position
		jmp	mul5

	roundup:
		shr	edx,11		; complet shift to proper position
		add	eax,1	; add rounding bit
		adc	edx,0
		test	edx,00200000h	; has mantissa overflowed (very unlikely)
		jz	mul5
		inc	cx		; if so increment exponent

	; Step 5: put it all together
	mul5:
		mov	ebx,ecx		; extract ...
		and	ebx,80000000h	; ... sign-bit
		cmp	cx,7ffh	; is exponent too big ?
		jge	overflow
		cmp	cx,0
		jle	underflow
		shl	ecx,20

		and	edx,000fffffh	; remove leading bit
		or	edx,ecx		; exponent added ...
		or	edx,ebx		; ... and sign-bit

	return:
		mov	dword ptr [rez+4],edx	; push upper dword of product ...
		mov	dword ptr [rez],eax	; ... and lower dword
		mov [valret],0
		jmp end
	underflow:
		mov [valret],2
		jmp end
	retz:
		xor	edx,edx		; clear upper dword
		xor	eax,eax		; ... and lower dword
		jmp	return

	overflow:
		mov [valret],1

	end:

	}
	*dest=rez;
	return valret;
}

int flib(int tip, double* src1, double* src2, double* dest)
{
	int ret=-1;
	switch (tip)
	{
	case 3:
		ret=asmMUL(*src1,*src2,dest);
		if (ret==-1)
		{
			perror("Floating point error:Overflow");
			return 1;
		}
		else if (ret==-2)
		{
			perror("Floating point error:Underflow");
			return 2;
		}
		else
			return 0;
		break;
	default:
		return -1;
	}
}