![[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.