Annotation of coherent/b/lib/libc/crt/i386/daddsub.s, revision 1.1

1.1     ! root        1: //////////
        !             2: / libc/crt/i386/daddsub.s
        !             3: / i386 C runtime library.
        !             4: / IEEE software floating point support.
        !             5: //////////
        !             6: 
        !             7: //////////
        !             8: / double _dadd(double d)
        !             9: / Return d + %edx:eax in %edx:%eax.
        !            10: /
        !            11: / double _dsub(double d)
        !            12: / Return d - %edx:eax in %edx:%eax.
        !            13: /
        !            14: / double _drsub(double d)
        !            15: / Return %edx:eax - d in %edx:%eax.
        !            16: /
        !            17: / Addition is easy: align the mantissa of the smaller operand with the larger,
        !            18: / possibly rounding up; add the aligned mantissas; shift right 1 bit
        !            19: / if necessary to normalize the result mantissa.
        !            20: /
        !            21: / Subtraction is a bit trickier: normalizing the result is nontrivial,
        !            22: / and we must be careful not to lose precision.  A simple example:
        !            23: /      A =     0x10 0000 0000 0001             53 bits
        !            24: /      B =      0x8 0000 0000 0000 . 8         53 bits
        !            25: /      A-B =    0x8 0000 0000 0000 . 8         53 bits
        !            26: / If we justify the mantissa of B with A before subtracting, we round it to
        !            27: /      B' =     0x8 0000 0000 0001             52 bits
        !            28: / and the result of the subtraction is
        !            29: /      A-B' =   0x8 0000 0000 0000             52 bits
        !            30: / so we have lost a bit of precision.
        !            31: / Shifting the larger operand left one bit eliminates this problem:
        !            32: /      A' =    0x10 0000 0000 0001 . 0         54 bits (1 after '.')
        !            33: /      A'-B =   0x8 0000 0000 0000 . 8         53 bits
        !            34: / which is the correct result.  Shifting once is enough, because subtraction
        !            35: / can clear at most the first result bit when operands are not aligned.
        !            36: //////////
        !            37: 
        !            38: d      .equ    8
        !            39: EXPMASK        .equ    0x7FF00000
        !            40: MANMASK        .equ    0x000FFFFF
        !            41: SGNMASK        .equ    0x80000000
        !            42: HIDDEN .equ    0x00100000
        !            43: DMBITS .equ    53
        !            44: 
        !            45:        .globl  _dadd
        !            46:        .globl  _dsub
        !            47:        .globl  _drsub
        !            48: 
        !            49: _drsub:
        !            50:        xchgl   %edx, 8(%esp)
        !            51:        xchgl   %eax, 4(%esp)           / exchange arg order
        !            52: /      jmp     _dsub                   / and fall through to subtract
        !            53: 
        !            54: _dsub:
        !            55:        xorl    %edx, $SGNMASK          / flip rhs sign
        !            56: /      jmp     _dadd                   / and fall through to add
        !            57: 
        !            58: _dadd:
        !            59:        push    %ebp
        !            60:        mov     %ebp, %esp
        !            61:        push    %esi
        !            62:        push    %edi
        !            63:        push    %ebx
        !            64:        push    %ecx
        !            65: 
        !            66:        movl    %esi, d+4(%ebp)
        !            67:        movl    %edi, d(%ebp)           / d to ESI:EDI
        !            68:                                        / now done with EBP as index register
        !            69: 
        !            70:        / Check for special cases +-0.0, +-infinity, NaN on each side.
        !            71:        / Ignore denormals.
        !            72:        movl    %ebx, %esi
        !            73:        andl    %ebx, $EXPMASK
        !            74:        jz      ?done                   / lhs 0.0, return rhs; ignore denormal
        !            75:        cmpl    %ebx, $EXPMASK
        !            76:        jz      ?retlhs                 / lhs +-infinity or NaN, return it
        !            77:        movl    %ecx, %edx
        !            78:        andl    %ecx, $EXPMASK
        !            79:        jz      ?retlhs                 / rhs 0.0, return lhs; ignore denormal
        !            80:        cmpl    %ecx, $EXPMASK
        !            81:        jz      ?done                   / rhs +-infinity or NaN, return it
        !            82: 
        !            83:        / Force arg with larger exponent to EDX:EAX.
        !            84:        shrl    %ebx, $20               / lhs biased exp in EBX
        !            85:        shrl    %ecx, $20               / rhs biased exp in ECX
        !            86:        cmpl    %ecx, %ebx
        !            87:        jge     ?expdiff
        !            88:        xchgl   %esi, %edx
        !            89:        xchgl   %edi, %eax              / force |EDX:EAX| >= |ESI:EDI| (maybe)
        !            90:        xchgl   %ebx, %ecx              / and exchange exponents accordingly
        !            91: 
        !            92:        / Compute the exponent difference in preparation for shifting ESI:EDI.
        !            93:        / We know that exponent(EDX:EAX) >= exponent(ESI:EDI).
        !            94: ?expdiff:
        !            95:        movl    %ebp, %ecx              / save EDX:EAX exponent in EBP
        !            96:        subl    %ecx, %ebx              / nonnegative exponent difference to ECX
        !            97:        cmpl    %ecx, $DMBITS+1
        !            98:        jg      ?done                   / ESI:EDI insignificant, return EDX:EAX as is
        !            99: 
        !           100:        movl    %ebx, %edx
        !           101:        xorl    %ebx, %esi              / sign bit 0 iff arg signs match
        !           102:        pushfl                          / save as add/sub flag
        !           103:        movl    %ebx, %edx
        !           104:        andl    %ebx, $SGNMASK          / save result sign bit in EBX
        !           105:        andl    %edx, $MANMASK
        !           106:        orl     %edx, $HIDDEN           / extract mantissa, restore hidden bit
        !           107:        andl    %esi, $MANMASK
        !           108:        orl     %esi, $HIDDEN           / extract mantissa, restore hidden bit
        !           109: 
        !           110:        / Align mantissas by shifting ESI:EDI to the right position.
        !           111:        jecxz   ?addsub                 / skip shift if exponents match
        !           112:        popfl                           / recover add/sub flag
        !           113:        pushfl
        !           114:        jns     ?align                  / addition, just adjust
        !           115: 
        !           116:        / For subtraction with differing exponents, shift args left one bit
        !           117:        / to avoid loss of precision; see comment at top.
        !           118:        shld    %esi, %edi, $1
        !           119:        shll    %edi, $1                / shift ESI:EDI left 1 bit
        !           120:        shld    %edx, %eax, $1
        !           121:        shll    %eax, $1                / shift EDX:EAX left 1 bit
        !           122:        decl    %ebp                    / adjust result exponent
        !           123: 
        !           124: ?align:
        !           125:        cmpl    %ecx, $32
        !           126:        jl      ?adj
        !           127:        movl    %edi, %esi
        !           128:        subl    %esi, %esi              / shift ESI:EDI right by 32 bits
        !           129: ?adj:
        !           130:        shrd    %edi, %esi, %cl
        !           131:        pushfl
        !           132:        shrl    %esi, %cl               / shift ESI:EDI right by CL mod 32
        !           133:        popfl                           / carry indicates rounding
        !           134:        adcl    %edi, $0                / round up
        !           135:        adcl    %esi, $0
        !           136:        
        !           137: ?addsub:
        !           138:        popfl                           / restore add/sub flag
        !           139:        js      ?sub
        !           140:        addl    %eax, %edi              / add lo mantissas
        !           141:        adcl    %edx, %esi              / and hi mantissas
        !           142:        cmpl    %edx, $HIDDEN<<1        / check for carry past hidden bit
        !           143:        jl      ?pack                   / done if none
        !           144: 
        !           145:        / Shift result mantissa right one bit.
        !           146: ?rshift:
        !           147:        incl    %ebp                    / bump result exponent
        !           148:        cmpl    %ebp, $EXPMASK
        !           149:        je      ?inf                    / exponent overflow, return +-infinity
        !           150:        shrd    %eax, %edx, $1
        !           151:        pushfl
        !           152:        shrl    %edx, $1                / shift result right one bit
        !           153:        popfl
        !           154:        adcl    %eax, $0
        !           155:        adcl    %edx, $0                / round up, carry to hidden bit impossible
        !           156: 
        !           157:        / Pack result mantissa from EDX:EAX, exponent from EBP, sign from EBX.
        !           158: ?pack:
        !           159:        andl    %edx, $MANMASK          / mask off hidden bit
        !           160:        orl     %edx, %ebx              / pack with sign
        !           161:        shll    %ebp, $20               / position exponent
        !           162:        orl     %edx, %ebp              / pack with exponent
        !           163:        
        !           164:        / Return the value from EDX:EAX.
        !           165: ?done:
        !           166:        pop     %ecx
        !           167:        pop     %ebx
        !           168:        pop     %edi
        !           169:        pop     %esi
        !           170:        pop     %ebp
        !           171:        ret
        !           172: 
        !           173:        / Return the value from ESI:EDI.
        !           174: ?retlhs:
        !           175:        xchgl   %esi, %edx
        !           176:        xchgl   %edi, %eax
        !           177:        jmp     ?done
        !           178: 
        !           179:        / Return +-infinity.
        !           180: ?inf:
        !           181:        subl    %eax, %eax
        !           182:        subl    %edx, %edx              / zero mantissa
        !           183:        jmp     ?pack
        !           184: 
        !           185:        / Subtract.
        !           186:        / We want the result mantissa to be nonnegative.
        !           187:        / The exponent check above assures that exp(EDX:EAX) >= exp(ESI:EDI),
        !           188:        / but when the exponents are equal we must compare mantissas.
        !           189: ?sub:
        !           190:        jcxz    ?argtest                / compare mantissas if exponents equal
        !           191: 
        !           192: ?dosub:
        !           193:        subl    %eax, %edi              / subtract lo mantissas
        !           194:        sbbl    %edx, %esi              / and hi mantissas
        !           195: 
        !           196:        / Normalize the difference.
        !           197:        / The code would be simpler if it used iterated left shift by 1
        !           198:        / to find hidden bit, but the reverse bit scan used here is faster.
        !           199:        / Bit 21 can be set if the EDX:EAX operand was shifted left above.
        !           200:        bsrl    %ecx, %edx              / hi bit position 0-21 to ECX
        !           201:        jz      ?lotest                 / hi mantissa is 0
        !           202:        cmpl    %ecx, $21
        !           203:        jz      ?rshift                 / right shift one bit if hi bit was 21
        !           204: 
        !           205: ?normalize:
        !           206:        negl    %ecx
        !           207:        addl    %ecx, $20               / shift count to CL
        !           208: 
        !           209:        / Shift left by CL bits to normalize, adjusting exponent accordingly.
        !           210: ?lshift:
        !           211:        shld    %edx, %eax, %cl
        !           212:        shll    %eax, %cl               / shift EDX:EAX left CL bits
        !           213:        subl    %ebp, %ecx              / adjust result exponent
        !           214:        jg      ?pack                   / pack and go home
        !           215:        subl    %eax, %eax              / exponent underflow, return 0.0
        !           216:        subl    %edx, %edx
        !           217:        jmp     ?done
        !           218:        
        !           219:        / Hi mantissa EDX is 0, check lo mantissa EAX.
        !           220: ?lotest:
        !           221:        bsrl    %ecx, %eax              / hi bit position 0-31 to ECX
        !           222:        jz      ?done                   / mantissa difference 0, return 0.0
        !           223:        cmpl    %ecx, $20
        !           224:        jle     ?losmall                / hi bit of lo mantissa is 0 to 20
        !           225:                                        / else hi bit is bit 21 to 31
        !           226:        negl    %ecx                    / -21 to -31
        !           227:        addl    %ecx, $52               / 31 to 21 shift count in CL
        !           228:        jmp     ?lshift                 / shift EDX:EAX as above
        !           229: 
        !           230:        / EDX is 0, hi bit of EAX is bit 0 to 20; move EAX up to EDX.
        !           231: ?losmall:
        !           232:        xchgl   %edx, %eax              / shift left 32 bits; EAX becomes 0
        !           233:        subl    %ebp, $32               / adjust exponent
        !           234:        jmp     ?normalize              / and shift as above
        !           235: 
        !           236:        / Compare mantissas, exchange args to force |EDX:EAX| >= |ESI:EDI|.
        !           237: ?argtest:
        !           238:        cmpl    %edx, %esi              / compare hi mantissas
        !           239:        ja      ?dosub
        !           240:        jb      ?flip
        !           241:        cmpl    %eax, %edi              / compare lo mantissas
        !           242:        jae     ?dosub
        !           243: 
        !           244:        / |ESI:EDI| > |EDX:EAX|, flip args and change result sign.
        !           245: ?flip:
        !           246:        xchgl   %esi, %edx
        !           247:        xchgl   %edi, %eax              / force |EDX:EAX| >= |ESI:EDI| (for sure)
        !           248:        xorl    %ebx, $SGNMASK          / and change the result sign
        !           249:        jmp     ?dosub
        !           250: 
        !           251: / end of libc/crt/i386/daddsub.s

unix.superglobalmegacorp.com

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.