|
|
1.1 ! root 1: ////////// ! 2: / libc/crt/i386/dmul.s ! 3: / i386 C runtime library. ! 4: / IEEE software floating point support. ! 5: ////////// ! 6: ! 7: ////////// ! 8: / double _dmul(double d) ! 9: / Return d * %edx:eax in %edx:%eax. ! 10: / ! 11: / In 32-bit dwords: ! 12: / A = hiA:loA ! 13: / B = hiB:loB ! 14: / A*B = hi(hiA*hiB):lo(hiA*hiB) ! 15: / + hi(hiA*loB):lo(hiA*loB) ! 16: / + hi(loA*hiB):lo(loA*hiB) ! 17: / + lo(loA*loB):lo(loA*loB) ! 18: / Multiplying two n bit quantities produces a 2*n-1 or 2*n bit result. ! 19: / Thus, multiplying two 53 bit quantities produces a 105 or 106 bit result. ! 20: / We want a normalized 53-bit result, so the lo-order dword above is irrelevant. ! 21: ////////// ! 22: ! 23: d .equ 8 ! 24: BIAS .equ 1023 ! 25: EXPMASK .equ 0x7FF00000 ! 26: MANMASK .equ 0x000FFFFF ! 27: SGNMASK .equ 0x80000000 ! 28: HIDDEN .equ 0x00100000 ! 29: ! 30: .globl _dmul ! 31: ! 32: _dmul: ! 33: push %ebp ! 34: movl %ebp, %esp ! 35: push %esi ! 36: push %edi ! 37: push %ebx ! 38: push %ecx ! 39: ! 40: movl %esi, d+4(%ebp) ! 41: movl %edi, d(%ebp) / d to ESI:EDI, call it A ! 42: / now done with EBP as index register ! 43: ! 44: / Check for special cases +-0.0 and +-infinity on each side. ! 45: movl %ebx, %esi ! 46: andl %ebx, $EXPMASK ! 47: jz ?retlhs / A is 0.0, return it; ignore denormal ! 48: cmpl %ebx, $EXPMASK ! 49: jz ?retlhs / A is +-infinity or NaN, return it ! 50: movl %ecx, %edx ! 51: andl %ecx, $EXPMASK ! 52: jz ?done / B is 0.0, return it; ignore denormal ! 53: cmpl %ecx, $EXPMASK ! 54: jz ?done / B is +-infinity or NaN, return it ! 55: ! 56: / Compute result sign. ! 57: movl %ebp, %edx ! 58: xorl %ebp, %esi ! 59: andl %ebp, $SGNMASK ! 60: push %ebp / save result sign bit ! 61: ! 62: / Compute result exponent. ! 63: shrl %ebx, $20 / A biased exp in EBX ! 64: shrl %ecx, $20 / B biased exp in ECX ! 65: subl %ecx, $BIAS / unbiased ! 66: addl %ecx, %ebx / form biased result exponent ! 67: jl ?zero / underflow, return 0.0 ! 68: cmpl %ecx, $EXPMASK>>20 ! 69: jge ?inf / overflow, return +-infinity ! 70: movl %ebp, %ecx / save result exponent in EBP ! 71: ! 72: / Perform the 64*64 bit multiply. ! 73: andl %esi, $MANMASK ! 74: orl %esi, $HIDDEN / extract A mantissa, restore hidden bit ! 75: andl %edx, $MANMASK ! 76: orl %edx, $HIDDEN / extract B mantissa, restore hidden bit ! 77: movl %ecx, %edx ! 78: movl %ebx, %eax / B mantissa to ECX:EBX ! 79: ! 80: mul %edi / loA*loB ! 81: push %edx / save hi(loA*loB), ignore lo dword ! 82: movl %eax, %ebx / loB ! 83: mul %esi / hiA*loB ! 84: xchgl %edi, %eax / lo(hiA*loB) to EDI, loA to EAX ! 85: movl %ebx, %edx / hi(hiA*loB) to EBX ! 86: mul %ecx / loA*hiB ! 87: xchgl %esi, %eax / lo(loA*hiB) to ESI, hiA to EAX ! 88: xchgl %ecx, %edx / hi(loA*hiB) to ECX, hiB to EDX ! 89: mul %edx / hiA*hiB, leave in EDX:EAX ! 90: ! 91: / Add partial products to obtain 96 bit product in EDX:EAX:ESI. ! 92: / This is really a 105-32=73 or 106-32=74 bit value, higher bits are 0. ! 93: addl %esi, %edi ! 94: adcl %eax, %ebx ! 95: adcl %edx, $0 ! 96: pop %edi / hi(loA*loB) ! 97: addl %esi, %edi ! 98: adcl %eax, %ecx ! 99: adcl %edx, $0 ! 100: ! 101: / Normalize the result mantissa. ! 102: / The product presently contains 73 or 74 bits. ! 103: / The normalized result contains 53+32 = 85 bits. ! 104: movb %cl, $12 / shift count to CL ! 105: testl %edx, $0x200 / examine product bit 74 ! 106: jz ?L0 / 73 bit product ! 107: decb %cl / 74 bit product, decrement shift count ! 108: incl %ebp / and bump the exponent ! 109: ! 110: ?L0: ! 111: / Shift EDX:EAX:ESI left by CL bits. ! 112: shld %edx, %eax, %cl ! 113: shld %eax, %esi, %cl ! 114: shll %esi, %cl ! 115: ! 116: / Round up result if appropriate. ! 117: shll %esi, $1 / high bit of lo 64 bits to CF ! 118: ?addcarry: ! 119: adcl %eax, $0 ! 120: adcl %edx, $0 ! 121: testl %edx, $HIDDEN<<1 ! 122: jnz ?carry ! 123: ! 124: / Pack result mantissa from EDX:EAX, exponent from EBP, sign from stack. ! 125: ?pack: ! 126: andl %edx, $MANMASK / mask off hidden bit ! 127: orl %ebp, %ebp ! 128: jle ?zero / exponent underflow, return 0.0 ! 129: cmp %ebp, $EXPMASK>>20 ! 130: jge ?inf / exponent overflow, return infinity ! 131: shll %ebp, $20 / position exponent ! 132: orl %edx, %ebp / pack mantissa and exponent ! 133: pop %ecx ! 134: orl %edx, %ecx / pack with sign ! 135: ! 136: / Return the value from EDX:EAX. ! 137: ?done: ! 138: pop %ecx ! 139: pop %ebx ! 140: pop %edi ! 141: pop %esi ! 142: pop %ebp ! 143: ret ! 144: ! 145: / Rounding generated carry past hidden bit, shift one bit right. ! 146: ?carry: ! 147: incl %ebp / bump the exponent ! 148: shrd %eax, %edx, $1 ! 149: pushfl ! 150: shrl %edx, $1 ! 151: popfl ! 152: jmp ?addcarry ! 153: ! 154: / Return the value from ESI:EDI. ! 155: ?retlhs: ! 156: xchgl %esi, %edx ! 157: xchgl %edi, %eax ! 158: jmp ?done ! 159: ! 160: ?inf: ! 161: pop %edx / pop result sign bit ! 162: orl %edx, $EXPMASK / max exp, zero mantissa for infinity ! 163: jmp ?L2 ! 164: ! 165: ?zero: ! 166: pop %edx / pop sign bit and ignore ! 167: subl %edx, %edx / return 0.0 ! 168: ?L2: ! 169: subl %eax, %eax ! 170: jmp ?done ! 171: ! 172: / end of libc/crt/i386/dmul.s
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.