![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ddiv.s CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [MW Coherent from dump] / coherent / b / lib / libc / crt / i386|
Return to ddiv.s CVS log | Up to [MW Coherent from dump] / coherent / b / lib / libc / crt / i386 |
1.1 ! root 1: ////////// ! 2: / libc/crt/i386/ddiv.s ! 3: / i386 C runtime library. ! 4: / IEEE software floating point support. ! 5: ////////// ! 6: ! 7: ////////// ! 8: / double _ddiv(double d) ! 9: / Return d / %edx:eax in %edx:%eax. ! 10: / ! 11: / double _drdiv(double d) ! 12: / Return %edx:eax / d in %edx:%eax. ! 13: / ! 14: / The hard part of floating point division is computing the result mantissa, ! 15: / i.e. the quotient of the numerator mantissa by the denominator mantissa; ! 16: / all are 53-bit quantities. ! 17: / We shift the denominator mantissa to make it as large as possible, ! 18: / then use i386 "divl" (64-bit by 32-bit unsigned integer divide) ! 19: / twice to find the two quotient dwords efficiently. ! 20: / ! 21: / This does not handle denormals, though it could without much trouble. ! 22: ////////// ! 23: ! 24: d .equ 8 ! 25: BIAS .equ 1023 ! 26: EXPMASK .equ 0x7FF00000 ! 27: MANMASK .equ 0x000FFFFF ! 28: SGNMASK .equ 0x80000000 ! 29: HIDDEN .equ 0x00100000 ! 30: ! 31: .globl _ddiv ! 32: .globl _drdiv ! 33: ! 34: _ddiv: ! 35: xchgl %edx, 8(%esp) ! 36: xchgl %eax, 4(%esp) / exchange arg order ! 37: / jmp _drdiv / and fall through to divide ! 38: ! 39: / Numerator is in EDX:EAX, call it A = hiA:loA. ! 40: / Denominator is on stack, call it B = hiB:loB. ! 41: / Compute the quotient A/B, call it Q = hiQ:loQ. ! 42: _drdiv: ! 43: push %ebp ! 44: movl %ebp, %esp ! 45: push %esi ! 46: push %edi ! 47: push %ebx ! 48: push %ecx ! 49: ! 50: movl %esi, d+4(%ebp) ! 51: movl %edi, d(%ebp) / B to ESI:EDI ! 52: / now done with EBP as index register ! 53: ! 54: / Compute result sign. ! 55: movl %ebx, %edx ! 56: xorl %ebx, %esi ! 57: andl %ebx, $SGNMASK ! 58: push %ebx / save result sign bit ! 59: ! 60: / Check for special cases +-0.0, +-infinity, NaN on each side. ! 61: movl %ecx, %esi ! 62: andl %ecx, $EXPMASK ! 63: movl %ebx, %edx ! 64: andl %ebx, $EXPMASK ! 65: jz ?lhszero / A is 0.0; ignore denormal ! 66: cmpl %ebx, $EXPMASK ! 67: jz ?lhsmax / A is +-infinity or NaN ! 68: orl %ecx, %ecx ! 69: jz ?inf / normal/0.0, return +-infinity; ignore denormal ! 70: cmpl %ecx, $EXPMASK ! 71: jz ?rhsmax / normal/+-infinity or normal/NaN ! 72: ! 73: / Compute probable result exponent in EBX. ! 74: / It might get incremented below (so 0 might become 1). ! 75: shrl %ebx, $20 / A biased exp in EBX ! 76: shrl %ecx, $20 / B biased exp in ECX ! 77: subl %ecx, $BIAS-1 / unbiased, fudged to get right result ! 78: subl %ebx, %ecx / form biased result exponent ! 79: jl ?zero / underflow, return 0.0 ! 80: cmpl %ebx, $EXPMASK>>20 ! 81: jge ?inf / overflow, return +-infinity ! 82: ! 83: / Extract the mantissas. ! 84: / Shift the denominator mantissa to range 2^63 <= B < 2^64. ! 85: andl %edx, $MANMASK ! 86: orl %edx, $HIDDEN / extract A mantissa, restore hidden bit ! 87: andl %esi, $MANMASK ! 88: orl %esi, $HIDDEN / extract B mantissa, restore hidden bit ! 89: shld %esi, %edi, $11 ! 90: shll %edi, $11 / shift left, B bit 31 is now 1 ! 91: jnz ?hard / loB is nonzero ! 92: ! 93: / Division is easy when the lo divisor dword is zero: ! 94: / perform a 96-bit by 32-bit divide of hiA:loA:0 by hiB to get ! 95: / a 64-bit quotient and a 32-bit remainder, ! 96: / then zero-extend the remainder to 64 bits. ! 97: / hiQ = q1 = A/hiB and r1 = A%hiB, then loQ = q2 = r1:0/hiB. ! 98: / This is a special case for efficiency; division by any double with ! 99: / up to 32 mantissa bits, notably any int or unsigned, is quite fast. ! 100: divl %esi / q1 = A/hiB to EAX, r1 = A%hiB to EDX ! 101: movl %ebp, %eax / save q1 ! 102: subl %eax, %eax / r1:0 in EDX:EAX ! 103: divl %esi / q2 = r1:0/hiB to EAX, r to EDX ! 104: xchgl %ebp, %edx / result quotient in EDX:EAX ! 105: subl %ecx, %ecx / r:0 in EBP:ECX ! 106: ! 107: / The quotient is in EDX:EAX, the remainder is in EBP:ECX. ! 108: / Round up the quotient when remainder >= B / 2, i.e. 2*r >= B. ! 109: ?remtest: ! 110: shld %ebp, %ecx, $1 / 2*hiR ! 111: jc ?roundup / too big ! 112: cmpl %ebp, %esi ! 113: jb ?position / 2*r < B ! 114: ja ?roundup / 2*r > B, round up ! 115: shll %ecx, $1 / 2*loR ! 116: cmpl %ecx, %edi / hi(2*R) == hiB, compare lo dword ! 117: jb ?position / 2*r < B ! 118: ! 119: / Round up the quotient in EDX:EAX. ! 120: ?roundup: ! 121: addl %eax, $1 ! 122: adcl %edx, $0 ! 123: ! 124: / The quotient in EDX:EAX might be correctly positioned as it stands, ! 125: / or it might require a 1-bit right shift. ! 126: ?position: ! 127: testl %edx, $HIDDEN<<1 / check if shift required ! 128: jz ?pack / no shift required ! 129: ?rshift: ! 130: incl %ebx / bump exponent ! 131: shrd %eax, %edx, $1 ! 132: pushfl / save CF for rounding ! 133: shrl %edx, $1 / shift EDX:EAX right by 1 bit ! 134: popfl / restore CF ! 135: adcl %eax, $0 ! 136: adcl %edx, $0 / round if appropriate ! 137: testl %edx, $HIDDEN<<1 / watch out for carry past hidden bit ! 138: jnz ?rshift ! 139: ! 140: / Pack result mantissa in EDX:EAX with exponent from EBX and sign from stack. ! 141: ?pack: ! 142: orl %ebx, %ebx ! 143: jle ?zero / exponent underflow, return 0.0 ! 144: cmp %ebx, $EXPMASK>>20 ! 145: jge ?inf / exponent overflow, return infinity ! 146: shll %ebx, $20 / position exponent ! 147: andl %edx, $MANMASK / mask off hidden bit ! 148: orl %edx, %ebx / pack mantissa and exponent ! 149: pop %ecx ! 150: orl %edx, %ecx / pack with sign ! 151: ! 152: ?done: ! 153: pop %ecx ! 154: pop %ebx ! 155: pop %edi ! 156: pop %esi ! 157: pop %ebp ! 158: ret ! 159: ! 160: / Numerator is 0.0 (or denormal, ignored here). ! 161: ?lhszero: ! 162: jecxz ?NaN / 0/0, return NaN; ignore denormal ! 163: cmpl %ecx, $EXPMASK ! 164: jnz ?zero / 0/normal, return 0.0 ! 165: ?rhsmax: ! 166: / Numerator is normal or 0, denominator is +-infinity or NaN. ! 167: andl %esi, $MANMASK ! 168: jnz ?NaN / A/NaN, return NaN ! 169: orl %edi, %edi ! 170: jnz ?NaN / A/NaN, return NaN ! 171: / jmp ?zero / 0/+-infinity or normal/+-infinity, return 0.0 ! 172: / fall through... ! 173: / Return +0.0. ! 174: ?zero: ! 175: pop %edx / pop sign bit and ignore ! 176: subl %edx, %edx / return 0.0 ! 177: ?zeroeax: ! 178: subl %eax, %eax ! 179: jmp ?done ! 180: ! 181: / Numerator is +-infinity or NaN. ! 182: ?lhsmax: ! 183: andl %edx, $MANMASK ! 184: jnz ?NaN / NaN/B, return NaN ! 185: orl %eax, %eax ! 186: jnz ?NaN / NaN/B, return NaN ! 187: cmpl %ecx, $EXPMASK ! 188: jz ?NaN / +-infinity/NaN or +-infinity/+-infinity ! 189: / jmp ?inf / +-infinity/normal or +-infinity/0, return infinity ! 190: / fall through... ! 191: ! 192: / Return +-infinity. ! 193: ?inf: ! 194: pop %edx / pop result sign bit ! 195: orl %edx, $EXPMASK / max exp, zero mantissa for infinity ! 196: jmp ?zeroeax ! 197: ! 198: / Return NaN. ! 199: ?NaN: ! 200: pop %edx / pop sign bit and ignore ! 201: movl %edx, $EXPMASK|MANMASK / max exp, nonzero mantissa for NaN ! 202: jmp ?zeroeax ! 203: ! 204: / The hard case performs a 128-bit by 64-bit division of hiA:loA:0:0 ! 205: / by hiB:loB to get a 64-bit quotient hiQ:loQ and a 64-bit remainder hiR:loR. ! 206: / Execute the following code twice to compute the two quotient dwords, ! 207: / saving hiQ in EBP and using it as a termination flag. ! 208: / Each iteration does a 96-bit by 64-bit divide to get a 32-bit quotient ! 209: / and a 64-bit remainder. ! 210: / Use q = A / hiB as a guess at the quotient, ! 211: / then decrement the remainder r = A % hiB by loB*q. ! 212: / We may need to decrement q once or twice to get the right quotient. ! 213: ! 214: ?hard: ! 215: push %ebx / save result exponent ! 216: subl %ebp, %ebp / clear flag for first pass ! 217: ?divide: ! 218: divl %esi / q = A/hiB to EAX, r = A%hiB to EDX ! 219: movl %ecx, %eax / save q in ECX ! 220: movl %ebx, %edx / and save r in EBX ! 221: mull %edi / loB*q to EDX:EAX ! 222: xchgl %edx, %ebx / r to EDX, hi(loB*q) to EBX ! 223: negl %eax / 0-lo(loB*q) to EAX ! 224: sbbl %edx, %ebx / EDX:EAX gets r-loB*q ! 225: jnc ?gotcha / q is the right quotient ! 226: ! 227: / q is too big a guess, adjust to q' = q-1 and adjust the remainder. ! 228: / This gets executed at most twice (because the high bit of B is 1). ! 229: ?adjust: ! 230: decl %ecx / decrement the quotient ! 231: addl %eax, %edi ! 232: adcl %edx, %esi / add B back to remainder ! 233: jnc ?adjust / repeat if still negative ! 234: ! 235: / The correct dword q is in ECX, the remainder r is in EDX:EAX. ! 236: ?gotcha: ! 237: orl %ebp, %ebp ! 238: jz ?loQ / repeat to find loQ ! 239: xchgl %ebp, %edx ! 240: xchgl %ecx, %eax / q to EDX:EAX, r to EBP:ECX ! 241: pop %ebx / restore result exponent ! 242: jmp ?remtest ! 243: ! 244: / Compute the lo dword of the quotient. ! 245: / r is in EDX:EAX, hiQ is in ECX. ! 246: ?loQ: ! 247: movl %ebp, %ecx / save hiQ (nonzero) in EBP ! 248: cmpl %edx, %esi / be wary of overflow on divide ! 249: jb ?divide / no overflow, proceed as above ! 250: ! 251: / Since r is the remainder of the first division by B, ! 252: / it must be strictly less than B. But it is possible that hiR = hiB, ! 253: / in which case the divide using "divl" would overflow. ! 254: / For this case, use q2 = 0xFFFFFFFF = 2^32 - 1 as the quotient guess. ! 255: / Compute the adjusted remainder r2 = r - q2*B as r - 2^32 * B + B. ! 256: movl %ecx, $0xFFFFFFFF / q2 ! 257: subl %eax, %edi / loR - loB, must be negative ! 258: / sbbl %edx, %esi / hiR = hiB, so hi dword EDX becomes 0 ! 259: movl %edx, %eax / EDX:EAX *= 2^32; EDX is now negative ! 260: movl %eax, %edi / 0 + loB ! 261: addl %edx, %esi / add back B to EDX:EAX ! 262: jc ?gotcha / remainder became positive, q2 is ok ! 263: jmp ?adjust / remainder stayed negative, adjust q2 ! 264: ! 265: / end of libc/crt/i386/ddiv.s
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.