|
|
subttl emfsqrt.asm - FSQRT instruction page;*******************************************************************************;emfsqrt.asm - FSQRT instruction; by Tim Paterson;; Microsoft Confidential;; Copyright (c) Microsoft Corporation 1991; All Rights Reserved;;Inputs:; edi = [CURstk];;Revision History:;; [] 09/05/91 TP Initial 32-bit version.;;*******************************************************************************
;A linear approximation of the square root function is used to get the;intial value for Newton-Raphson iteration. This approximation gives;nearly 5-bit accuracy over the required input interval, [1,4). The;equation for the linear approximation of y = sqrt(x) is y = mx + b,;where m is the slope (named SQRT_COEF) and b is the y-intercept (named;SQRT_INTERCEPT).;;(The values for m and b were computed with Excel Solver in two passes: ;the first pass computed them full precision, minimizing absolute error;;the second computed only b after m was rounded to an 8-bit value.);;The resulting values have the following maximum error:;;inp. value --> 1 2.18972 3.82505;----------------------------------------------------------------;abs. err., full prec. 0.04544 -0.03233 0.04423;abs. err., truncated 0.04544 -0.04546 0.04423;;The three input values shown represent the left end point, the maximum ;error (derivative of absolute error == 0), and the right end point. ;The right end point is not 4 because the approximation reaches 2.000;at the value given--we abandon the linear approximation at that point;and use that same value for all greater input values. This linear;approximation is computed with 8-bit operations, so truncations can;add a negative error. This increases maximum error only when it is;already negative, as shown in the table.;;Each iteration of Newton-Raphson approximation more than doubles the;number of bits of accuracy. Suppose the current guess is A, and it has;an absolute error of e (i.e., A+e or A-e is the root). Then the absolute;error after the next iteration is e^2/2A. This error is always positive.;However, the divide instruction truncates, which introduces an error;that is always negative. Sometimes a constant or rounding bit is added;to balance the positive and negative errors. The maximum possible error ;is given in comments below for each iteration. (Note that when we compute ;the error from e^2/2A, A could be in the range 1 to 2--we use 1 to get;max error.) Remember that the binary point is to the RIGHT of the MSB;when looking at these error numbers.
;SQRT_INTERCEPT is used when the binary point is to the right of the MSB.;Multiplying it by 64K would put the binary point to the left of the MSB,;so it must be divided by two to be aligned.SQRT_INTERCEPT equ 23185 ; 0.70755 * 65536 / 2
;SQRT_COEF would have the binary point to the left of the MSB if multiplied;by 256. However, this would leave it with a leading zero, so we multiply;it by two more to normalize it.SQRT_COEF equ 173 ; 0.33789 * 256 * 2
SqrtSpcl: cmp al,bTAG_DEN jz SqrtDen cmp al,bTAG_INF jnz SpclDestNotDen;Have infinity or ah,ah ;Is it negative? js ReturnIndefiniteSqrtRet: ret
MaxStartRoot:;The first iteration is calculated as (ax / bh) * 100H + bx. The first ;trial root in bx should be 10000H (which is too big). But it's very;easy to calculate (ax / 100H) * 100H + 10000H = ax. mov bx,ax cmp ax,-1 ;Would subsequent DIV overflow? jb FirstTrialRoot;The reduced argument is so close to 4.0 that the 16-bit DIV instruction;used in the next iteration would overflow. If the argument is 4-A ;then a guess of 2.0 is in error by approximately A/4. [This is not;an upper bound. The error is a little by more than this by an;addition with the magnitude of A^2. This is an insignificant amount;when A is small.] This means that the first guess of 2.0 is quite;accurate, and we'll use it to bypass some of the iteration steps. ;This will eliminate the DIV overflow by skipping the DIV.;;One iteration is performed by: (Arg/Guess + Guess)/2. When Guess = 2,;this becomes (Arg/2 + 2)/2 = Arg/4 + 1. We get Arg/2 just by assuming;the binary point is one bit further left; then a single right shift is;needed to get Arg/4. By shifting in a 1 bit on the left, we account for;adding 1 at the same time. [Note that if Arg = 4 - A, then Arg/4 + 1; = (4 - A)/4 + 1 = 1 - A/4 + 1 = 2 - A/4. In other words, we just;subtract out exactly what we estimate our error to be, A/4.];;Since the upper 16 bits are 0FFFFH, A <= 2^-14, so error <= 2^-16 =; +0.00001526, -0. mov ebx,esi ;Return root in ebx sar ebx,1 ;Trial root = arg/2 cmp esi,ebx ;Will 32-bit division overflow? jb StartThirdIteration ;No, our 32-bit guess is good;Argument is really, really close to 4.0: with an initial trial root of;2.0, max absolute error is 2^-32 = +2.328E-10, -0. One trivial;iteration will get us 65-bit accuracy, max abs. error = +2.71E-20, -0. mov ebx,esi mov eax,ecx ;65-bit root*2 in ebx:eax (MSB implied) shl ecx,2 ;ecx = low half*4 jmp RoundRoot
SqrtDen: mov EMSEG:[CURerr],Denormal test EMSEG:[CWmask],Denormal ;Is denormal exception masked? jnz SqrtRet ;If not, quit
;******EM_ENTRY eFSQRTeFSQRT:;****** mov eax,EMSEG:[edi].ExpSgn cmp al,bTAG_ZERO jz SqrtRet ja SqrtSpcl or ah,ah js ReturnIndefinite mov esi,EMSEG:[edi].lManHi mov ecx,EMSEG:[edi].lManLo sar EMSEG:[edi].wExp,1 ;Divide exponent by two mov edi,0 ;Extend mantissa jc RootAligned ;If odd exponent, leave it normalized shrd edi,ecx,1 shrd ecx,esi,1 shr esi,1 ;Denormalize, extending into ediRootAligned:;esi:ecx:edi has mantissa, 2 MSBs are left of binary point. Range is [1,4). shld eax,esi,16 ;Get high word of mantissa movzx ebx,ah ;High byte to bl;UNDONE: MASM 6 bug!!;UNDONE: SQRT_COEF (=0AEH) get sign extended!! mov dx,SQRT_COEF ;UNDONE imul bx,dx ;UNDONE;UNDONE imul bx,SQRT_COEF ;Product in bx;Multiply by SQRT_COEF causes binary point to shift left 1 bit. add bx,SQRT_INTERCEPT ;5-bit approx. square root in bh jc MaxStartRoot;Max absolute error is +/- 0.04546 div bh ;See how close we are add bh,al ;quotient + divisor (always sets CY)FirstTrialRoot:;Avoid RCR because it takes 9 clocks on 386. Use SHRD (3 clocks) instead. mov dl,1 ;Need bit set shrd bx,dx,1 ;(quotient + divisor)/2;bx has 9-bit approx. square root, normalized;Max absolute error is +0.001033, -0.003906 movzx eax,si shld edx,esi,16 ;dx:ax has high half mantissa div bx ;Test our approximation add ebx,eax ;quotient + divisor shl ebx,15 ;Normalize (quotient + divisor)/2;ebx has 17-bit approx. square root, normalized;Max absolute error is +0.000007629, -0.00001526;Add adjustment factor to center the error range at +/-0.00001144 or bh,20H ;Add in 0.000003815StartThirdIteration: mov edx,esi mov eax,ecx div ebx ;Test approximation stc ;Set bit for rounding (= 2.328E-10) adc ebx,eax ;quotient + divisor + round bit;Avoid RCR because it takes 9 clocks on 386. Use SHRD (3 clocks) instead. mov dl,1 ;Need bit set shrd ebx,edx,1 ;(quotient + divisor)/2, rounded;ebx has 32-bit approx. square root, normalized;Max absolute error is +2.983E-10, -2.328E-10 mov edx,esi ;Last time we need high half mov eax,ecx shld ecx,edi,2 ;ecx = low half*4, w/extension back in div ebx ;Test approximation xchg edi,eax ;Save 1st quotient, get extension mov esi,eax or esi,edx ;Any remainder? jz HaveRoot ;Result is ebx:esi div ebx ;edi:eax is 64-bit quotient add ebx,edi ;quotient + divisor (always sets CY)RoundRoot: mov esi,eax ;Save low half root*2
;We have 65-bit root*2 in ebx:esi (eax==esi) (MSB is implied one).;Max absolute error is +4.450E-20, -5.421E-20. This maximum error ;corresponds to just less than +/- 1 in the last (65th) bit. ; ;We have to determine if this error is positive or negative so;we can tell if we rounded up or down (and set the status bit;accordingly). This is done by squaring the root and comparing the;that result with the input.;;Squaring the sample root requires summing partial products:; lo*lo + lo*hi + hi*lo + hi*hi. lo*hi == hi*lo, so only one multiply;is needed there. The low half of lo*lo isn't relevant, we know it;is non-zero. Only the low few bits of hi*hi are needed, so we can use;an 8-bit multiply there. Since the MSB is implied, we need to add in;two 1*lo products (shifted up 64 bits). We only need bits 64 - 71 of;the 130-bit product (the action happens near bit 65). What we're ;squaring is root*2, so the result is square*4. ecx already has arg*4.
mul eax ;Low partial product of square mov edi,edx ;Only high half counts mov eax,ebx mul esi ;Middle partial product of square add eax,eax ;There are two of these adc edx,edx add edi,eax adc edx,0 ;edx:edi = lo*lo + lo*hi + hi*lo add edx,esi ;lo*implied msb add edx,esi ;lo*implied msb again mov al,bl mul al ;hi*hi - only low 8 bits are valid add al,dl ;Bits 64 - 71 of product or al,1 ;Account for sticky bits 0 - 63 sub cl,al ;Compare product with argument;Sign flag set if product is larger. In this case, subtract 1 from root. add cl,cl ;Set CY if sign is setSubOneFromRoot: sbb esi,0 ;Reduce root if product was too big sbb ebx,0ShiftRoot:;ebx:esi = root*2;Absolute error is in the range (0, -5.421E-20). This is equivalent to;less than +1, -0 in last bit. Thus LSB is correct rounding bit as ;long as we set a sticky bit below it.;;Now divide root*2 by 2, preserving LSB as rounding bit and filling;eax with 1's as sticky bits.;;Avoid RCR because it takes 9 clocks on 386. Use SHRD (3 clocks) instead. mov eax,-1 shrd eax,esi,1 ;Move round bit to MSB of eax shrd esi,ebx,1 shrd ebx,eax,1 ;Shift 1 into MSB of ebxStoreRoot: mov edi,EMSEG:[CURstk] mov EMSEG:[Result],edi mov ecx,EMSEG:[edi].ExpSgn;mantissa in ebx:esi:eax, exponent in high ebx, sign in bh bit 7 jmp EMSEG:[RoundMode]
HaveRoot:;esi = eax = edx = 0 cmp edi,ebx ;Does quotient == divisor? jz StoreRoot ;If so, we're done;Quotient != divisor, so answer is not exact. Since remainder is zero,;the division was exact. The only error in the result is e^2/2A, which;is always positive. We need the error to be only negative so that;the rounding routine can properly tell if it rounded up. add ebx,edi ;quotient + divisor (always sets CY) jmp SubOneFromRoot ;Reduce root to ensure negative error
|
