|
|
1.1 root 1: /*
2: * Copyright (c) 1999 Apple Computer, Inc. All rights reserved.
3: *
4: * @APPLE_LICENSE_HEADER_START@
5: *
6: * Portions Copyright (c) 1999 Apple Computer, Inc. All Rights
7: * Reserved. This file contains Original Code and/or Modifications of
8: * Original Code as defined in and that are subject to the Apple Public
9: * Source License Version 1.1 (the "License"). You may not use this file
10: * except in compliance with the License. Please obtain a copy of the
11: * License at http://www.apple.com/publicsource and read it before using
12: * this file.
13: *
14: * The Original Code and all software distributed under the License are
15: * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
16: * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
17: * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
18: * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
19: * License for the specific language governing rights and limitations
20: * under the License.
21: *
22: * @APPLE_LICENSE_HEADER_END@
23: */
24:
25: .file "lipsq.s"
26:
27: .ident "@(#)kern-fp:lipsq.s 1.1"
28:
29: // ***************************************************************************
30: //
31: // l i p s q . m o d
32: // =================
33: //
34: // ===============================================================
35: // intel corporation proprietary information
36: // this software is supplied under the terms of a license
37: // agreement or non-disclosure agreement with intel corporation
38: // and may not be copied nor disclosed except in accordance with
39: // the terms of that agreement.
40: // ===============================================================
41: //
42: // functions:
43: // implements the loading of the constants:
44: // one, log base 2 of ten, log base 2 of e, pi,
45: // log base 10 of 2, log base e of 2, and zero.
46: // implements the 80387 square root instruction.
47: // implements the 80387 integer part instruction.
48: //
49: // public procedures:
50: // load_con sqrt intpt
51: //
52: // *****************************************************************************
53: //
54: //...september 16, 1983...
55: //
56: // .file *a_mli*
57: //
58: //$nolist
59: #include "fp_e80387.h"
60: //$list
61: .text //a_med segment er public
62: //
63: // extrn put_si_result,sticky_right_shift,round
64: // extrn addition_normalize,left_shift,getx
65: // extrn set_up_indefinite,clear_6w,set_6w
66: // extrn test_4w,left_shift_result_cl
67: // extrn left_shift_frac1_cl,left_shift_frac2_cl
68: // extrn add_to_frac_2,gradual_underflow
69: // extrn subtraction_normalize,put_indefinite
70: // extrn i_masked_,set_i_masked_,set_d_masked_
71: // extrn affine_infinity_,get_precision
72: //
73: .globl load_con
74: .globl intpt
75: .globl sqrt
76: //
77: // temp real floating point numbers for push constant instructions
78: //
79: treal_table:
80: .value 0x00000,0x00000,0x00000,0x00000,0x08000,0x03fff //treal_one
81: .value 0x04000,0x08afe,0x0cd1b,0x0784b,0x0d49a,0x04000 //treal_l2t
82: .value 0x0c000,0x0f0bb,0x05c17,0x03b29,0x0b8aa,0x03fff //treal_l2e
83: .value 0x0c000,0x0c234,0x02168,0x0daa2,0x0c90f,0x04000 //treal_pi
84: .value 0x0a000,0x0f798,0x0fbcf,0x09a84,0x09a20,0x03ffd //treal_lg2
85: .value 0x0e000,0x079ab,0x0d1cf,0x017f7,0x0b172,0x03ffe //treal_ln2
86: .value 0x00000,0x00000,0x00000,0x00000,0x00000,0x00000 //treal_0
87: //$eject
88: // ***************************************************************
89: // load_con
90: // ******
91: // function:
92: // implements the load constants instructions.
93: //
94: // inputs:
95: // offset of fpn pointer in bx register.
96: //
97: // outputs:
98: // constant value on top of stack
99: //
100: // data accessed:
101: // - offset_result
102: //
103: // data changed:
104: // - mem_operand_pointer
105: //
106: // procedures called:
107: // getx
108: //
109: // *******************************************************************
110: // all load constant instructions use the same entry point.
111: //
112: ALIGN
113: load_con: //proc near
114: //
115: // the contents of bx is used as an index to treal_table.
116: // getx unpacks the extended floating point number and put_result
117: // pushes the value to the top of the 80387 stack.
118: //
119: shrb $1,%bl //take care of extra shi
120: mov $0x06,%eax //treal_table ptr = cs:
121: mulb %bl // index * 6 - (6 * 2)
122: add $(treal_table-(load_1_op*12)),%eax // * (load_1_op)
123: mov %eax, %esi
124: movl $0, result_dword_frac(%ebp)
125: movw %cs:(%esi), %ax
126: movw %ax, result_dword_frac+frac80(%ebp)
127: mov %cs:2(%esi), %eax
128: mov %eax, result_dword_frac+frac64(%ebp)
129: mov %cs:6(%esi), %eax
130: mov %eax, result_dword_frac+frac32(%ebp)
131: movzwl %cs:10(%esi), %eax
132: mov %eax, dword_result_expon
133: xor %eax, %eax
134: mov %eax, result_sign(%ebp)
135: movb prec64, %dl
136: mov $offset_result, %edi
137: call round
138: mov $offset_result, %edi
139: call addition_normalize
140: andb $~a_mask, %gs:sr_flags
141: xor %eax, %eax
142: mov before_error_signals(%ebp), %eax
143: lahf
144: mov $offset_result, %edi
145: jmp put_si_result
146: //load_con endp
147: //$eject
148: // ***********************************************************************
149: // intpt:
150: // *****
151: // function:
152: // implements 80387 integer part instruction
153: //
154: // inputs:
155: //
156: // outputs:
157: //
158: // data accessed:
159: // - offset_result offset_operand1
160: // - tag1 expon1
161: // - word_frac1
162: //
163: // data changed:
164: // - tag1 expon1
165: // - result
166: //
167: // procedures called:
168: // get_operand gradual_underflow
169: // round subtraction_normalize
170: // put_si_result test_4w
171: // i_masked? set_i_masked?
172: // set_d_masked?
173: //
174: // **************************************************************************
175: ALIGN
176: intpt: //proc
177: // movb true,%ah / load constant
178: jz separate_cases // branch if no stack error
179: call set_stk_u_error
180: testb invalid_mask, %gs:sr_masks
181: jz intpt_done
182: jmp put_op1_result //math stack error, so return indef
183: ALIGN
184: separate_cases:
185: andb $~a_mask, %gs:sr_flags //initialize a-bit to zero
186: movb tag1(%ebp), %al // load op1 tag
187: cmpb valid,%al
188: je valid_case
189: cmpb denormd,%al
190: je denormalized_operand
191: cmpb unsupp,%al
192: jne check_nan
193: movl $offset_operand1, %edi
194: call set_up_indefinite
195: set_i_err:
196: orb invalid_mask, %gs:sr_errors
197: testb invalid_mask, %gs:sr_masks
198: jz intpt_done
199: jmp put_op1_result // masked_i error so return indef
200: ALIGN
201: intpt_done:
202: ret
203: ALIGN
204: check_nan:
205: cmpb inv, %al
206: je kind_of_nan_
207: jmp put_op1_result //infinity or zero, same answer
208: ALIGN
209: kind_of_nan_:
210: testb $0x40, msb_frac1
211: jz make_qnan
212: jmp put_op1_result // op1 is a qnan so pass it through
213: ALIGN
214: make_qnan:
215: orb $0x40, msb_frac1 //op1 is an snan so make it a qnan
216: jmp set_i_err //and signal i_error
217: ALIGN
218: denormalized_operand:
219: orb denorm_mask, %gs:sr_errors
220: testb denorm_mask, %gs:sr_masks
221: jz intpt_done
222: movw $0x0001,expon1(%ebp) //if masked d-error, make valid
223: movb valid,tag1(%ebp)
224: valid_case:
225: movl $0x403e,%eax
226: cmpl %eax,dword_expon1 //if expon >=63, then number
227: jge give_op1 // is already an integer
228: mov $offset_operand1,%edi //gradual uflow until expon=63
229: push %edi
230: call gradual_underflow
231: pop %edi //round to precision 64
232: movb prec64,%dl
233: movb false,%al
234: call round
235: cmpb true, rnd1_inexact
236: jne detect_zero
237: orb $inexact_mask, %gs:sr_errors
238: cmpb true,added_one
239: jne detect_zero
240: orb $a_mask, %gs:sr_flags
241: detect_zero:
242: xor %eax,%eax
243: movl $dword_frac1+4,%edi //if fraction = 0, result = 0
244: call test_4w
245: jz zero_result
246: mov $offset_operand1,%edi //normalize
247: call subtraction_normalize
248: jmp put_op1_result
249: ALIGN
250: zero_result:
251: movl %eax,dword_expon1 //set result to true zero
252: movb special,tag1(%ebp)
253: give_op1:
254: jmp put_op1_result
255: //intpt endp
256: //$eject
257: // ***************************************************************************
258: // sqrtx:
259: //
260: // function:
261: // fractional square root routine.
262: //
263: // inputs:
264: // assumes valid, non-zero,positive, normalized
265: // fraction is in frac1.
266: //
267: // outputs
268: // leaves fractional square root in result.
269: //
270: // data accessed:
271: // - offset_operand1 lsb_frac1
272: // - offset_operand2 word_frac2
273: // - offset_result lsb_frac2
274: // - result_word_frac lsb_result
275: // - msb_frac1
276: //
277: // data changed:
278: // - word_frac2 lsb_frac2
279: // - result_word_frac lsb_result
280: //
281: // procedures called:
282: // left_shift_result_cl left_shift_frac1_cl
283: // left_shift_frac2_cl clear_6w
284: // set_6w add_to_frac_2
285: //
286: // ***********************************************************************
287: ALIGN
288: sqrtx: //proc
289: movl $result_dword_frac,%edi //during this computation,
290: call clear_6w //the lsb of the result will
291: mov $dword_frac2,%edi // hold g and s, and the msb
292: call set_6w // will hold carry-out bits
293: and $0x000000ff,%ecx //clear high word
294: stc //cf holds the quotient bit
295: push $65 // iterate 65 times
296: pushf // stack the quotient bit
297: jmp enter_sqrt_loop
298: ALIGN
299: sqrt_loop:
300: push %ecx
301: pushf // stack the q_bit
302: movb $1,%cl
303: call left_shift_result_cl //shift result left one bit
304: popf // inject the new q_bit into
305: pushf // the least significant byte
306: adcb %cl,%cl
307: orb %cl,1+lsb_result
308: movb $2,%cl // into the lsb
309: call left_shift_frac2_cl //shift frac2 left 2 bits
310: enter_sqrt_loop:
311: movb msb_frac1,%al //g and s bits of frac2 <--
312: andb $0x0c0,%al // top 2 bits of frac1
313: movb %al,lsb_frac2
314: movb $2,%cl // shift frac1 left 2
315: call left_shift_frac1_cl
316: movb $0x0c0,%al
317: movb 1+lsb_result,%ah
318: popf //test q_bit
319: jc q_bit_set
320: shl $16,%eax
321: call add_to_frac2 //frac2.gs <-- frac2.gs +
322: jmp set_q_bit //result.11
323: q_bit_set:
324: notb %ah //frac2.gs <-- frac2.gs +
325: shl $16,%eax //take care of extra low word
326: add %eax,dword_frac2(%ebp) //not(result).11
327: mov result_dword_frac+frac64(%ebp),%eax
328: not %eax
329: adc %eax,dword_frac2+frac64(%ebp)
330: mov result_dword_frac+frac32(%ebp),%eax
331: not %eax
332: adc %eax,dword_frac2+frac32(%ebp)
333: set_q_bit:
334: pop %ecx // reload loop count
335: loop sqrt_loop // loop until done
336: rcrb $1,%cl // set g bit of result to q_bit
337: orb %cl,lsb_result
338: movw result_word_frac,%ax //frac2 <- frac2 +
339: incb %ah // result + 1
340: shl $16,%eax
341: call add_to_frac2
342: mov $dword_frac2+frac64,%edi //if frac2 = 0 then
343: xor %eax,%eax //s_bit of result = 0,
344: call test_4w //otherwise 1.
345: orb 1+lsb_frac2,%al
346: jz left_adjust_result
347: orb $0x40,lsb_result
348: left_adjust_result:
349: movb $8,%cl //shift result left 8 bits
350: jmp left_shift_result_cl // to eliminate carry
351: add_to_frac2:
352: mov $result_dword_frac,%esi // si points to the addend
353: mov $dword_frac2,%edi // di points to the result
354: jmp add_to_frac_2 // add result frac to frac2
355: //sqrtx endp
356: //$eject
357: // ***********************************************************************
358: // sqrt:
359: //
360: // function:
361: // implements the 80387 sqrt instruction
362: //
363: // inputs:
364: // assume op1 is set up
365: //
366: // outputs:
367: // result
368: //
369: // data accessed:
370: // - result_record_offset result_expon
371: // - offset_operand1 tag1
372: // - sign1 expon1
373: // - msb_frac1 offset_operand1
374: // - offset_result result_sign
375: // - result_tag
376: //
377: // data changed:
378: // - expon1 result_sign
379: // - result_tag result_expon
380: //
381: // procedures called:
382: // set_up_indefinite sticky_right_shift
383: // sqrtx round
384: // addition_normalize affine_infinity?
385: // set_i_masked? get_precision
386: // set_d_masked? put_si_result
387: //
388: // ***********************************************************************
389: ALIGN
390: sqrt: //proc
391: jz sqrt_cont // if stack error, sqrt done
392: call set_stk_u_error
393: testb invalid_mask, %gs:sr_masks
394: jnz put_op1 //masked stack error, so return indef
395: sqrt_done:
396: ret
397: ALIGN
398: sqrt_cont:
399: andb $~a_mask, %gs:sr_flags
400: movb tag1(%ebp),%al // load tag for op1
401: cmpb valid,%al
402: jne op1_denorm_
403: cmpb positive, sign1(%ebp)
404: jne i_error
405: jmp sqrt_valid_case
406: op1_denorm_:
407: cmpb denormd,%al //if op1 denormalized, then give
408: jne op1_zero_
409: cmpb positive, sign1(%ebp)
410: jne i_error
411: jmp d_error
412: op1_zero_:
413: cmpb special,%al //if op1 = 0, then give 0 as the
414: je put_op1 //result
415: cmpb unsupp,%al
416: je i_error
417: cmpb inv,%al //if op1 inv, then give i_error
418: jne inf_op1
419: jmp kind_of_nan_
420: inf_op1:
421: cmpb positive, sign1(%ebp) //infinity case
422: je put_op1 // -infinity is invalid
423: i_error:
424: call set_i_masked_ //denormalized, unnormalized,
425: jz sqrt_done // negative, or proj +infinity
426: mov $offset_result,%edi
427: call set_up_indefinite //if masked i_error, then give
428: jmp sqrt_give_result //indefinite
429: put_op1: //if op1 = zero, nan, or +inf,
430: jmp give_op1 // then give op1 as the result
431: d_error:
432: call set_d_masked_ //op1 is denormalized
433: jz sqrt_done
434: movl $offset_operand1, %edi
435: call norm_denorm //if d_error masked, make valid
436: sqrt_valid_case:
437: subl $exponent_bias,dword_expon1
438: testl $0x0001,dword_expon1
439: jz even_expon
440: decl dword_expon1 //if expon1 odd, then expon1 <--
441: jmp halve_exponent // expon1 - 1
442: even_expon:
443: movb $1,%cl //if expon1 even, then shift
444: xorb %al,%al // frac1 right one bit
445: mov $offset_operand1,%edi
446: call sticky_right_shift
447: halve_exponent:
448: movl dword_expon1,%eax
449: sarl $1,%eax
450: addl $exponent_bias,%eax
451: movl %eax,dword_result_expon
452: call sqrtx // calculate fraction
453: sqrt_round_result:
454: xorl %eax,%eax // not second rounding
455: movl %eax,result_sign(%ebp) // set sign and tag
456: movl $offset_result,%edi // round result
457: call get_precision
458: call round
459: mov $offset_result,%edi
460: call addition_normalize //(possible renormalize)
461: cmpb true, rnd1_inexact
462: jne sqrt_give_result
463: orb $inexact_mask, %gs:sr_errors
464: cmpb true, added_one
465: jne sqrt_give_result
466: orb $a_mask, %gs:sr_flags
467: sqrt_give_result:
468: mov $offset_result,%edi
469: jmp put_si_result
470: //sqrt endp
471: //
472: //a_med ends
473: //
474: // end
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.