![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to daddsub.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 daddsub.s CVS log | Up to [MW Coherent from dump] / coherent / b / lib / libc / crt / i386 |
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
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.