📄 emuoperalib.s
字号:
incl %ecx cmpw $0x7fff,%cx je add19 add15: orw %cx,%cx #renormalize if needed jle add28 ret add20: sbbl (%edi),%eax #subtract mantissas sbbl 4(%edi),%edx jnc add24 xorl $0x80000000,%ecx notl %edx notl %eax negb %bl sbbl $-1,%eax sbbl $-1,%edx add24: jns add25 #weed out not-normal, possible all-zero orw %cx,%cx jle add19 ret add25: je add27 #here needs normalizing, or is zero add26: decw %cx addb %bl,%bl adcl %eax,%eax adcl %edx,%edx jns add26 orw %cx,%cx jle add19 ret add28: orl %edx,%edx #this may be 0+0 jnz add19 movb %bl,%dl orl %eax,%edx jnz add22 btsl $16,%ecx ret add27: orb %bl,%dl #here may be all zero orl %eax,%edx jz addzer add22: xorl %edx,%edx add19: call renorm #normalize ret addzer: movl $0x00010000,%ecx #result went all the way to zero movb h_ctrl+1,%bh andb $0x0c,%bh cmpb $0x04,%bh jnz add18 btsl $31,%ecx add18: ret # multiply registers to stack# input: cx:edx:eax = EP number# result: registers = ST * number .globl epmul,epmul1mulzer: shrl $16,%ecx movb $1,%cl shll $16,%ecx xorl %edx,%edx jmp mulin2 mulinf: orl $0x00027fff,%ecx #return INF movl $0x80000000,%edx mulin2: xorl %eax,%eax xorb %bl,%bl jmp mul21 mulw: call retnan #weed out NaN arguments movl %ecx,%esi #weed out special cases addl 8(%edi),%esi btrl $31,%esi cmpl $0x37fff,%esi jz stdnan testb $1,%bl jnz mulzer cmpw $0x7fff,%si jnc mulinf call tnormal #normalize arguments if needed xchgl %eax,(%edi) xchgl %edx,4(%edi) xchgl %ecx,8(%edi) call tnormal jmp epmul1 epmul: shldl $16,%ecx,%ebx #weed out special arguments orb 10(%edi),%bl jnz mulw epmul1: pushl %ebp #save exponent, sign, FP stack pointer pushl %ecx movl %eax,%esi #save A0123 movl %edx,%ecx mull (%edi) # A23 * B23 xchgl %eax,%esi movl %edx,%ebx mull 4(%edi) # A23 * B01 addl %eax,%ebx adcl $0,%edx movl %edx,%ebp movl (%edi),%eax # A01 * B23 mull %ecx addl %eax,%ebx adcl %edx,%ebp movl $0,%eax adcl $0,%eax xchgl %ecx,%eax mull 4(%edi) # A01 * B01 addl %ebp,%eax adcl %ecx,%edx popl %ecx #exponent js mul42 #possibly normalize addl %ebx,%ebx adcl %eax,%eax adcl %edx,%edx decl %ecx mul42: addl $-1,%esi #set rounding bits adcl %ebx,%ebx jz mul43 movb $0x80,%bl mul43: rcrb $1,%bl popl %ebp #restore base pointer subw $tBIAS-2,%cx #calculate new exponent addw 8(%edi),%cx jo mul44 decw %cx jle mul45 mul21: roll $8,%ecx #put in sign xorb 11(%edi),%cl rorl $8,%ecx ret #returnmul44: stc #here re-normalize jmp mul46 mul45: clc mul46: call renorm jmp mul21 # divide registers to stack# input: cx:edx:eax = EP number# result: registers = ST / number .globl epdiv,epdivr,epdiv1divw61: orb $zerod,h_stat call exproc jmp mulinf divw: call retnan #weed out NaN arguments cmpw $0x7fff,%cx #weed out A = INF jnz divw3 cmpw %cx,8(%edi) jz stdnan jmp mulinf divw3: cmpw $0x7fff,8(%edi) #weed out B = INF jz mulzer btl $16,%ecx #if A = 0 and B = 0 return NaN jnc divw7 testb $1,10(%edi) jnz stdnan jmp mulzer divw7: testb $1,10(%edi) #if B = 0 "zero divide" jnz divw61 call tnormal #normalize arguments xchgl %eax,(%edi) xchgl %edx,4(%edi) xchgl %ecx,8(%edi) call tnormal xchgl %eax,(%edi) xchgl %edx,4(%edi) xchgl %ecx,8(%edi) jmp epdiv1 epdivr: xchgl (%edi),%eax xchgl 4(%edi),%edx xchgl 8(%edi),%ecx epdiv: shldl $16,%ecx,%ebx #weed out special arguments orb 10(%edi),%bl jnz divw epdiv1: pushl %ebp #save exponent, sign pushl %ecx movl $1,%ebp #first bit = 0 cmpl 4(%edi),%edx #jump if x < y jc div50 subl 0(%edi),%eax #get first bit sbbl 4(%edi),%edx jnc div49 movl %eax,%edx #prevent overflow xorl %eax,%eax movl %eax,%ecx jmp div31 div49: decl %ebp #first bit = 1div50: divl 4(%edi) # q = A0123 / B01 -> ebx movl %eax,%ecx # r -> esi movl %edx,%esi mull 0(%edi) # x * B23 -> edx:eax xchgl %edx,%esi #subtract that from qA23 negl %eax sbbl %esi,%edx jnc div41 div31: decl %ecx addl 0(%edi),%eax adcl 4(%edi),%edx jnc div31 div41: cmpl 4(%edi),%edx #prevent overflow jnz div51 movl %eax,%esi xorl %eax,%eax movl %eax,%ebx subl 0(%edi),%esi jmp div32 div51: divl 4(%edi) # x = A0123 / B01 -> bx movl %eax,%ebx movl %edx,%esi mull 0(%edi) #then q * B23 -> edx:eax negl %eax #subtract that from qA23 sbbl %edx,%esi jnc div42 div32: decl %ebx addl 0(%edi),%eax adcl 4(%edi),%esi jnc div32 div42: movl %esi,%edx #set rounding bit, first bit flag orl %eax,%edx xchgl %ebx,%edx jz div54 movb $0x80,%bl addl %eax,%eax adcl %esi,%esi jc div43 subl 0(%edi),%eax sbbl 4(%edi),%esi cmc div43: rcrb $1,%bl div54: movl %edx,%eax #load result to edx:eax movl %ecx,%edx popl %ecx #restore registers cmpl $1,%ebp #handle first bit jnc div57 rcrl $1,%edx rcrl $1,%eax rcrb $1,%bl incw %cx div57: popl %ebp subw 8(%edi),%cx #calculate new exponent addw $tBIAS-0,%cx jo div23 decw %cx jle div24 div21: roll $8,%ecx #put in sign xorb 11(%edi),%cl rorl $8,%ecx ret div23: stc #renormalize if needed jmp div26 div24: clc div26: call renorm jmp div21 # partial remainder, Intel version# input: cx:edx:eax = EP number# result: registers = ST % number .globl epremres88: cmpw $-1,%cx #here exponent of A < exponent of B jnz rem83 # return A unless A > B/2 pushl %eax # in which case return A-B pushl %edx shll $1,(%edi) rcll $1,4(%edi) btrl $31,%edx jmp res17 rem87: btl $16,%ebx #branch if IEEE remainder jc res88 rem83: jmp rem16 remw: andb $0x0fb,h_stat+1 #clear C2 call retnan #remove NaNs and invalids andb $0x0b8,h_stat+1 #clear condition codes btl $16,%ecx #if B = 0 or A = INF return NaN jc stdnan cmpw $0x7fff,%si jz stdnan call tnormal #normalize arguments xchgl %eax,(%edi) xchgl %edx,4(%edi) xchgl %ecx,%esi cmpw $0x7fff,%si jz rem17 btl $16,%ecx jc rem19 call tnormal jmp rem3 eprem: movl 8(%edi),%esi shldl $16,%ecx,%ebx #weed out special arguments orb 10(%edi),%bl jnz remw andb $0x0b8,h_stat+1 #clear condition codes xchgl (%edi),%eax #swap registers xchgl 4(%edi),%edx xchgl %esi,%ecx rem3: movb $0,%bl #initialize quotient subw %si,%cx #get the shift count js rem87 cmpw $64,%cx #maximum shift count 64 jl rem4 addw %cx,%si orb $32,%cl andw $63,%cx subw %cx,%si orb $0x04,h_stat+1 rem4: movw %si,8(%edi) pushl %ebp #save registers movl (%edi),%ebp movl 4(%edi),%esi rem5: subl %ebp,%eax #divide using shift-subtract sbbl %esi,%edx jnc rem9 addl %ebp,%eax adcl %esi,%edx rem9: cmc rem11: adcb %bl,%bl decb %cl js rem15 addl %eax,%eax adcl %edx,%edx jnc rem5 subl %ebp,%eax sbbl %esi,%edx jmp rem11 rem15: popl %ebp #restore registers testb $0x04,h_stat+1 #if C2 is set (incomplete) jnz rem16 # leave Q zero testl $0x10000,%ebx #branch if IEEE remainder jnz res18 rem41: andl $7,%ebx #massage quotient for the required form rorw $1,%bx # bits are discontinuous and reversed shrb $4,%bh rorw $1,%bx shrb $2,%bh shlb $7,%bl shll $1,%ebx orb %bh,h_stat+1 rem16: movw 8(%edi),%cx #normalizerem17: xorb %bl,%bl call renorm rem19: ret res18: pushl %eax #save remainder pushl %edx res17: addl %eax,%eax #if next quotent bit would be 0, we are adcl %edx,%edx # done setc %bh subl (%edi),%eax sbbl 4(%edi),%edx sbbb $0,%bh jc res42 movl %ebx,%esi #next bit is one, if all the rest are 0, andl $1,%esi # and Q is even, we are likewise done orl %edx,%esi orl %eax,%esi jz res42 incl %ebx #here we round up the quotient popl %edx popl %eax subl (%edi),%eax sbbl 4(%edi),%edx notl %edx negl %eax sbbl $-1,%edx btcl $31,%ecx jmp rem41 res42: popl %edx popl %eax jmp rem41 # square-root# input: ecx:edx:eax = EP number# result: ST = sqrt(number) .globl epsqrtsqr89: ret sqrw: xorl %ebx,%ebx btl $16,%ecx #zero returns as such jc sqr89 orw %cx,%cx #check for unsupported jz sqrw5 orl %edx,%edx jns stdnan cmpw $0x7fff,%cx #check for NaN jnz sqrw5 movl %edx,%ebx addl %ebx,%ebx orl %eax,%ebx jnz nanret sqrw5: orl %ecx,%ecx #negative is illegal js stdnan cmpw $0x7fff,%cx #+INF returns as such jz sqr89 call tnormal #normalize movl %ecx,8(%edi) jmp sqr2 epsqrt: testl $0x80030000,%ecx #jump if special value jnz sqrw sqr2: pushl %ebp #save registers pushl %edi xorl %ebx,%ebx #save x, possibly shifted right by 1 movl %eax,%esi andb $1,%cl shrdl %cl,%esi,%ebx shrdl %cl,%edx,%esi pushl %esi pushl %ebx incl %ecx #load high 30 or 31 bits -> esi:edi shrdl %cl,%edx,%eax shrl %cl,%edx movl %edx,%esi movl %eax,%edi testb $1,%cl #get first estimate with a*x+b movl $0x4d12cfe4,%ebx jz sqr3 movl $0x6cff9300,%ebx sqr3: shldl $17,%esi,%eax mulw %bx shrl $16,%ebx addw %dx,%bx movl %esi,%eax #16-bit iteration y = (x/y + y)/2 shldl $16,%esi,%edx divw %bx shrw $1,%bx addw %ax,%bx sbbw $0,%bx shll $16,%ebx movl %esi,%edx #32-bit iteration y = (x/y + y)/2 movl %edi,%eax divl %ebx shrl $1,%ebx addl %eax,%ebx sbbl $0,%ebx movl %esi,%edx #64-bit iteration y = (x/y + y)/2 movl %edi,%eax divl %ebx movl %eax,%ecx xorl %eax,%eax divl %ebx xorl %edx,%edx shrl $1,%ebx rcrl $1,%edx addl %edx,%eax adcl %ebx,%ecx movl %eax,%ebx #mantissa of y -> ecx:ebx mull %eax #calculate y^2 -> (none):si:edi:ebp movl %eax,%ebp movl %edx,%edi movl %ecx,%eax mull %eax movl %eax,%esi movl %ebx,%eax mull %ecx addl %eax,%edi adcl %edx,%esi addl %eax,%edi adcl %edx,%esi movl %ebx,%eax #mantissa of y -> ecx:eax popl %ebx #now y^2 - x into si:edi:ebp subl %ebx,%edi popl %ebx sbbl %ebx,%esi js sqr37 jz sqr35 sqr31: subl $1,%eax #try out next smallest y sbbl $0,%ecx stc #new residue = old - 2y + 1 sbbl %eax,%ebp sbbl %ecx,%edi sbbl $0,%esi subl %eax,%ebp sbbl %ecx,%edi sbbl $0,%esi jns sqr31 #more if we are still above real sqrt(x) cmc sqr32: sbbb %bl,%bl #1/2 bit here, never exactly half orb $0x40,%bl sqr33: movl %ecx,%edx popl %edi #restore registers popl %ebp movl 8(%edi),%ecx #calculate new exponent shrw $1,%cx adcw $0x1fff,%cx ret sqr37: stc #try next largest x adcl %eax,%ebp adcl %ecx,%edi adcl $0,%esi addl %eax,%ebp adcl %ecx,%edi adcl $0,%esi jns sqr32 #quit if differential > 0 addl $1,%eax adcl $0,%ecx jmp sqr37 sqr35: movl %edi,%ebx #see if we have an exact answer orl %ebp,%ebx jz sqr33 jmp sqr31 #__lib endp#emuseg ends# end⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -