开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 嵌入式/单片机编程 > VxWorks > 查看源码
emuOperALib.s
资源名称：vxwork_src.rar　[点击查看]
上传用户：nvosite88
上传日期：2007-01-17
资源大小：4983k
文件大小：34k
源码类别：

VxWorks

开发平台：

C/C++

emuOperALib.s：源码内容
  1. /* Copyright 1991-2001 Wind River Systems, Inc. */
  3. /*
  4. modification history
  5. --------------------
  6. 01c,14jul99,tdl  added FUNC / FUNC_LABEL
  7. 01b,24oct96,yp   stuck in a # in USSC mailing addr so cpp will work
  8. 01a,24jan95,hdn  original US Software version.
  9. */
  11. #* * * * * * * * * *
  12. #       Filename:   EMUOPER.ASM
  13. #
  14. #       Copyright (C) 1990,1991 By
  15. #       # United States Software Corporation
  16. #       # 14215 N.W. Science Park Drive
  17. #       # Portland, Oregon 97229
  18. #
  19. #       This software is furnished under a license and may be used
  20. #       and copied only in accordance with the terms of such license
  21. #       and with the inclusion of the above copyright notice.
  22. #       This software or any other copies thereof may not be provided
  23. #       or otherwise made available to any other person.  No title to
  24. #       and ownership of the software is hereby transferred.
  25. #
  26. #       The information in this software is subject to change without
  27. #       notice and should not be construed as a commitment by United
  28. #       States Software Corporation.
  29. #
  30. #       Version:        V3.09  01/19/95 GOFAST  WRS GNU version of GF-PROT
  31. #       Released:       1 March 1991
  32. #
  33. #* * * * * * * * * *
  35. #define _ASMLANGUAGE
  36. #include "vxWorks.h"
  37. #include "asm.h"
  39. .data
  40. .globl  FUNC(copyright_wind_river)
  41. .long   FUNC(copyright_wind_river)
  43. .include "emuIncALib.s"
  44. .text
  47. #        extrn   roundi:near,rounde:near,chkarg:near
  48. #        extrn   retnan:near,tnormal:near,renorm:near,getqn:near
  49. #        extrn   mul10:near,shiftr:near,exproc:near,exprop:near
  50. #        extrn   exret:near
  52. #DGROUP  GROUP   hwseg
  53. #CGROUP  GROUP   emuseg
  55. #hwseg   segment RW para public 'DATA'
  56. #        assume  ds:DGROUP, es:DGROUP, ss:DGROUP
  57. #        extrn   h_ctrl:word,h_stat:word,s_iptr:dword,s_memp:dword
  58. #hwseg   ends
  60. #emuseg  segment public USE32 'CODE'
  61. #        assume  cs:CGROUP
  63. __lib:
  65. #        subttl  Convert To Internal Floating-Point Format
  66. #        page
  68. # convert 16-bit or 32-bit integer to extended floating-point
  69. # input:        *eax = signed integer
  70. # output:       EP in ecx:edx:eax
  72. .globl witoep,wiep
  73. witoep:
  74. movswl (%eax),%edx
  75. jmp wiep2
  77. .globl sitoep
  78. sitoep:
  79. movl (%eax),%edx #pick up number
  81. wiep2: #SAVEII
  82. movl %esi,s_iptr
  83. movl %eax,s_memp
  85. wiep: movl $31,%ecx #initial sign -> cx
  86. xorl %eax,%eax
  88. orl %edx,%edx #weed out zero
  89. jz liep20
  91. jns siep5 #if negative, set sign, negate
  92. orl $0x80000000,%ecx
  93. negl %edx
  95. siep5: bsrl %edx,%esi
  96. subl %esi,%ecx
  97. shll %cl,%edx
  98. negw %cx
  99. addw $tBIAS+31,%cx
  100. ret
  102. # convert 64-bit integer to extended floating-point
  103. # input:        *eax = signed integer
  104. # output:       EP in ecx:edx:eax
  106. .globl litoep
  107. litoep:
  108. #SAVEII
  109. movl %esi,s_iptr
  110. movl %eax,s_memp
  112. movl 4(%eax),%edx #pick up number
  113. movl (%eax),%eax
  115. movl %eax,%ecx #jump if all zero
  116. orl %edx,%ecx
  117. jz liep20
  119. movl $tBIAS+63,%ecx #initial exponent -> cx
  121. orl %edx,%edx #if negative, set sign, negate
  122. jns liep5
  123. orl $0x80000000,%ecx
  124. notl %edx
  125. negl %eax
  126. sbbl $-1,%edx
  127. js liep7
  129. liep5: decl %ecx #shift left until high bit = 1
  130. addl %eax,%eax
  131. adcl %edx,%edx
  132. jns liep5
  134. liep7: ret
  136. liep20: movl $0x00010000,%ecx #zero
  137. ret
  139. # convert 10-byte packed decimal integer to extended floating-point
  140. # input:        *eax = signed integer
  141. # output:       EP in ecx:edx:eax
  143. .globl pdtoep
  144. pdtoep:
  145. movl $9,%ecx #set up loop counter and pointer
  146. lea 9(%eax),%edi
  148. xorl %edx,%edx #initialize mantissa
  149. xorl %eax,%eax
  151. pdep3: call mul10 #multiply mantissa by 10
  153. decl %edi #pick up next 2-digit
  154. movb (%edi),%bl
  156. pushl %ebx #add low digit to mantissa
  157. shrb $4,%bl
  158. andl $0x0f,%ebx
  159. addl %ebx,%eax
  160. adcl $0,%edx
  162. call mul10 #again multiply mantissa by 10
  164. popl %ebx #add high digit to mantissa
  165. andl $0x0f,%ebx
  166. addl %ebx,%eax
  167. adcl $0,%edx
  169. loop pdep3 #loop for all digits
  171. movl %eax,%ecx #jump if all zero
  172. orl %edx,%ecx
  173. jz pdep20
  175. movl $tBIAS+63,%ecx #initial exponent -> cx
  177. pdep5: decl %ecx #shift left until high bit = 1
  178. addl %eax,%eax
  179. adcl %edx,%edx
  180. jns pdep5
  182. pdep7: movb 9(%edi),%bl #set the sign bit
  183. andb $0x80,%bl
  184. shll $24,%ebx
  185. orl %ebx,%ecx
  186. ret
  188. pdep20: movl $0x00010000,%ecx #zero
  189. jmp pdep7
  191. # convert single-precision to extended floating-point
  192. # input:        *eax = single-precision
  193. # output:       EP in ecx:edx:eax
  195. .globl fptoep
  196. fptoep:
  197. #SAVEII
  198. movl %esi,s_iptr
  199. movl %eax,s_memp
  201. movl (%eax),%edx #pick up number
  203. movl %edx,%ecx #sign, exponent to ecx
  204. sarl $23,%ecx
  205. andl $0x800000ff,%ecx
  207. shll $9,%edx #mantissa -> edx:eax
  208. xorl %eax,%eax
  210. cmpb $0,%cl #go handle zero and denormal
  211. jz fpep20
  213. cmpb $-1,%cl #go handle INF and NaN
  214. jz dpep30
  215. rcrl $1,%edx
  217. fpep7: addw $tBIAS-sBIAS,%cx #adjust bias
  218. ret
  220. fpep20: shrl $1,%edx #this is either zero or denormal
  221. jz dpep29 #   zero if mantissa = 0
  222. orb $denorm,h_stat
  223. call exprop
  224. incl %ecx
  225. fpep21: decw %cx
  226. addl %edx,%edx
  227. jns fpep21
  228. jmp fpep7
  230. # convert double-precision to extended floating-point
  231. # input:        eax -> double-precision
  232. # output:       EP in ecx:edx:eax
  234. .globl dptoep
  235. dptoep:
  236. #SAVEII
  237. movl %esi,s_iptr
  238. movl %eax,s_memp
  240. movl 4(%eax),%edx #pick up number
  241. movl (%eax),%eax
  243. movl %edx,%ecx #sign, exponent to ecx
  244. sarl $20,%ecx
  245. andl $0x800007ff,%ecx
  247. shldl $12,%eax,%edx #mantissa -> edx:eax
  248. shll $11,%eax
  249. orw %cx,%cx
  250. jz dpep20
  251. cmpw $0x7ff,%cx
  252. jz dpep30
  253. rcrl $1,%edx
  255. dpep7: addw $tBIAS-lBIAS,%cx #adjust exponent bias
  256. dpep9: ret
  258. dpep30: orl $0x00027fff,%ecx #this is INF or NaN
  259. movl %edx,%ebx #   INF if mantissa 800..00
  260. shrdl $1,%ecx,%edx
  261. orl %eax,%ebx
  262. jz dpep9
  263. nanret: btl $30,%edx
  264. jc dpep32
  265. orw $invop,h_stat
  266. call exproc
  267. btsl $30,%edx
  268. dpep32: xorl %ebx,%ebx
  269. ret
  271. dpep20: shrl $1,%edx #this is either zero or denormal
  272. movl %edx,%ebx #   zero if mantissa = 0
  273. orl %eax,%ebx
  274. jz dpep29
  275. orb $denorm,h_stat
  276. call exprop
  277. incl %ecx
  278. dpep21: decw %cx
  279. addl %eax,%eax
  280. adcl %edx,%edx
  281. jns dpep21
  282. jmp dpep7
  284. dpep29: btsl $16,%ecx #zero
  285. ret
  287. # convert extended floating-point number to internal formal
  288. # input:        EP in *ss:eax
  289. # output:       EP in ecx:edx:eax
  291. .globl eptoep
  292. eptoep:
  293. movl 4(%eax),%edx #pick up the number
  294. movzwl 8(%eax),%ecx
  295. movl (%eax),%eax
  296. addw %cx,%cx
  297. rcrl $1,%ecx
  299. jz epe11 #tags: 10 = INF, denormal, illegal
  300. cmpw $0x7fff,%cx #      01 = zero
  301. jz epe6
  302. orl %edx,%edx
  303. jns epe6
  304. movl %edx,%ebx
  305. orl %eax,%ebx
  306. jz epe6
  307. ret
  309. epe6: btsl $17,%ecx #illegal etc
  310. ret
  312. epe11: movl %edx,%ebx #this is zero or denormal
  313. orl %eax,%ebx
  314. jnz epe6
  315. btsl $16,%ecx
  316. ret
  318. #        subttl  Convert From Internal Floating-Point Format
  319. #        page
  321. # convert extended floating-point to 16 or 32 bit integer
  322. # input:        EP in ecx:edx:eax
  323. # output:       eax = signed integer
  325. .globl eptowi
  326. eptowi:
  327. movl $15,%ebx #integer size
  328. jmp epsi1
  330. .globl eptosi
  331. eptosi:
  332. movl $31,%ebx #integer size
  334. epsi1: testl $0x00020000,%ecx #weed out specials
  335. jnz epsi20
  337. epsi3: subw $tBIAS,%cx #subtract the bias
  339. cmpw %bx,%cx #if exponent too large, return INF
  340. jge epsi32
  342. subw $63,%cx #here shift right
  343. negw %cx
  344. call shiftr
  346. call roundi #IEEE rounding
  347. adcl $0,%eax
  349. orl %ecx,%ecx #if negative take 2s complement
  350. jns epsi11
  351. negl %eax
  353. epsi11: ret #return
  355. epsi20: rorl $8,%ecx #check for stack underflow
  356. cmpb $3,%ch
  357. roll $8,%ecx
  358. jz epsi80
  360. orw %cx,%cx #specials, but denormals go back
  361. jz epsi3
  362. jmp epsi30
  364. epsi32: jg epsi30 #overflow
  365. addl %edx,%edx
  366. orl %eax,%edx
  367. jz epsi34
  369. epsi30: orb $invop,h_stat #either +INF or -INF
  370. call exproc
  371. epsi34: movl $1,%eax
  372. movb %bl,%cl
  373. shll %cl,%eax
  374. ret
  376. epsi80: orb $stacku,h_stat #stack underflow
  377. call exproc
  378. jmp epsi34
  380. # convert extended floating-point to 64-bit integer
  381. # input:        EP in ecx:edx:eax
  382. # output:       edx:eax = signed integer
  384. .globl eptoli
  385. eptoli:
  386. testl $0x00020000,%ecx #weed out specials
  387. jnz epli20
  389. epli3: subw $tBIAS,%cx #subtract the bias
  391. cmpw $63,%cx #if exponent too large, return INF
  392. je epli32
  393. jg epli30
  395. subw $63,%cx #here shift right
  396. negw %cx
  397. call shiftr
  399. call roundi #do IEEE rounding
  400. adcl $0,%eax
  401. adcl $0,%edx
  403. orl %ecx,%ecx #if negative take 2s complement
  404. jns epli11
  405. notl %edx
  406. negl %eax
  407. sbbl $-1,%edx
  409. epli11: ret #return
  411. epli20: shldl $16,%ecx,%ebx #check for stack underflow
  412. cmpb $3,%bl
  413. jz epli80
  415. orw %cx,%cx #specials, but denormals go back
  416. jz epli3
  417. jmp epli30
  419. epli32: addl %edx,%edx #treat 0x80000000 as special case
  420. orl %eax,%edx
  421. jz epli34
  423. epli30: orb $invop,h_stat #either +INF or -INF
  424. call exproc
  425. epli34: movl $0x80000000,%edx
  426. xorl %eax,%eax
  427. jmp epli11
  429. epli80: orb $stacku,h_stat #stack underflow
  430. call exproc
  431. jmp epli34
  433. # convert EP to 10-byte packed-decimal integer
  434. # input:        EP in ecx:edx:eax
  435. #               edi pointer to result area
  436. # output:       signed integer in result area
  438. .globl eptopd
  440. eppd20: shldl $16,%ecx,%ebx #check for stack underflow
  441. cmpb $3,%bl
  442. jz eppd80
  444. orw %cx,%cx #specials, but denormals go back
  445. jz eppd3
  447. eppd30: orb $invop,h_stat #either +INF or -INF
  448. call exproc
  449. eppd34: movl $0xffff8000,%eax
  450. movl %eax,6(%edi)
  451. movl $6,%ecx
  452. repz
  453. stosb
  454. ret
  456. eppd80: orb $stacku,h_stat #stack underflow
  457. call exproc
  458. jmp eppd34
  460. eppd33: subl $9,%edi #number too big
  461. jmp eppd30
  463. eptopd:
  464. movl %ecx,6(%edi) #store the sign
  466. btl $17,%ecx #weed out specials
  467. jc eppd20
  469. eppd3: subw $tBIAS+63,%cx #subtract the bias
  470. jge eppd30
  471. negw %cx
  472. call shiftr
  474. call roundi #do IEEE rounding
  475. adcl $0,%eax
  476. adcl $0,%edx
  478. movl %edx,%esi #get the digits by dividing with 10
  479. movl $9,%ecx
  480. pushl %ebp
  481. movl $10,%ebp
  483. eppd11: xchgl %esi,%eax
  484. xorl %edx,%edx
  485. divl %ebp
  486. xchgl %esi,%eax
  487. divl %ebp
  488. movb %dl,%bl
  489. xchgl %esi,%eax
  490. xorl %edx,%edx
  491. divl %ebp
  492. xchgl %esi,%eax
  493. divl %ebp
  494. shlb $4,%dl
  495. orb %bl,%dl
  496. movb %dl,(%edi)
  497. incl %edi
  498. loop eppd11
  500. popl %ebp
  501. orl %eax,%eax
  502. jnz eppd33
  503. ret
  505. # convert extended floating-point to single-precision floating-point
  506. # input:        EP in ecx:edx:eax
  507. # output:       eax = floating-point
  509. .globl eptofp
  511. epfp2: shldl $16,%ecx,%ebx #check for stack underflow
  512. cmpb $3,%bl
  513. jz epfp38
  515. btsl $16,%ecx #check zero, denormal
  516. jc epfp9
  517. orw %cx,%cx
  518. jz epfp3
  519. orl %edx,%edx
  520. jns epfp36
  521. jmp epfp30
  523. eptofp:
  524. testl $0x00030000,%ecx #jump if special value
  525. jnz epfp2
  527. epfp3: addw $sBIAS-tBIAS,%cx #adjust the bias
  529. cmpw $0,%cx #weed out underflow and overflow
  530. jle epfp20
  531. cmpw $0x0ff,%cx
  532. jnc epfp32
  534. epfp6: addl %edx,%edx #shift left to delete hidden bit
  536. epfp7: movl %edx,%ebx #set rounding bits
  537. shll $23,%ebx
  538. addl $-1,%eax
  539. adcl %ebx,%ebx
  540. jz epfp71
  541. movb $0x80,%bl
  542. epfp71: rcrb $1,%bl
  544. epfp8: shldl $9,%ecx,%eax #combine sign, exponent, mantissa
  545. movb %cl,%al
  546. shldl $23,%edx,%eax
  548. orb %cl,%cl #set underflow if necessary
  549. jnz epfp11
  550. testb $underf,h_ctrl
  551. jz epfp13
  552. orb %bl,%bl
  553. jz epfp11
  554. epfp13: orb $underf,h_stat
  555. call exproc
  557. epfp11: call roundi #IEEE rounding
  558. adcl $0,%eax
  559. ret
  561. epfp9: addl %ecx,%ecx #return signed zero
  562. rcrl $1,%eax
  563. ret
  565. epfp20: negw %cx #underflow, we return zero or denormal
  566. jz epfp7
  567. cmpw $25,%cx
  568. jc epfp23
  569. movw $25,%cx
  570. epfp23: xorl %ebx,%ebx
  571. shrdl %cl,%edx,%ebx
  572. shrl %cl,%edx
  573. xorw %cx,%cx
  574. orl %ebx,%ebx
  575. jz epfp7
  576. orb $1,%dl
  577. jmp epfp7
  579. epfp30: btcl $31,%edx #this is NaN or INF
  580. orl %edx,%eax
  581. jnz epfp34
  582. movw $2*sBIAS+1,%cx
  583. jmp epfp7
  585. epfp32: orb $overf,h_stat #this is overflow
  586. call exproc
  587. movw $2*sBIAS,%cx
  588. movl $-1,%edx
  589. jmp epfp7
  591. epfp38: call getqn
  592. orb $stacku,h_stat
  593. jmp epfp39
  594. epfp36: call getqn #this is NaN
  595. jmp epfp37
  596. epfp34: btl $30,%edx
  597. jc epfp35
  598. epfp37: orb $invop,h_stat
  599. epfp39: call exproc
  600. btsl $30,%edx
  601. epfp35: movw $2*sBIAS+1,%cx
  602. addl %edx,%edx
  603. xorl %ebx,%ebx
  604. jmp epfp8
  606. # convert extended floating-point to double-precision floating-point
  607. # input:        EP in ecx:edx:eax
  608. # output:       edx:eax = double-precision
  610. .globl eptodp
  612. epdp32: orb $overf,h_stat #this is overflow
  613. call exproc
  614. movw $0x7fe,%cx
  615. movl $-1,%edx
  616. movl %edx,%eax
  617. jmp epdp6
  619. epdp2: shldl $16,%ecx,%ebx #check for stack underflow
  620. cmpb $3,%bl
  621. jz epdp38
  623. btsl $16,%ecx #check zero, denormal
  624. jc epdp3
  625. orw %cx,%cx
  626. jz epdp4
  627. orl %edx,%edx
  628. jns epdp36
  629. jmp epdp30
  631. epdp3: movl %ecx,%edx #return zero
  632. andl $0x80000000,%edx
  633. ret
  635. eptodp:
  636. testl $0x00030000,%ecx #jump if special value
  637. jnz epdp2
  639. epdp4: addw $lBIAS-tBIAS,%cx #adjust the bias
  640. jle epdp20 #   weed out underflow and overflow
  641. cmpw $0x07ff,%cx
  642. jnc epdp32
  644. epdp6: addl %eax,%eax #shift left mantissa to delete hidden bit
  645. adcl %edx,%edx
  647. xorl %ebx,%ebx #discarded bits into ebx
  648. shrdl $12,%eax,%ebx
  650. shrdl $12,%edx,%eax #shift right everything 12 bits
  651. shrdl $11,%ecx,%edx
  652. btl $31,%ecx
  653. rcrl $1,%edx
  655. epdp8: orl %ebx,%ebx #exact = no rounding
  656. jz epdp19
  658. testb $0x0c,h_ctrl+1 #quick round
  659. jnz epdp14
  660. btl $0,%eax
  661. adcl $-1,%ebx
  662. jns epdp17
  663. epdp18: #PREEX   precis+rup
  664. orw $precis+rup,h_stat
  665. testb $precis,h_ctrl
  666. jnz LL1
  667. orw $0x8080,h_stat
  668. movl $exret,-4(%ebp)
  669. LL1:
  670. addl $1,%eax
  671. adcl $0,%edx
  672. ret
  673. epdp17: #PREEXS  precis
  674. orb $precis,h_stat
  675. testb $precis,h_ctrl
  676. jnz LL2
  677. orw $0x8080,h_stat
  678. movl $exret,-4(%ebp)
  679. LL2:
  680. epdp19: andb $0x0fd,h_stat+1
  681. ret
  683. epdp14: roll $8,%ebx #IEEE rounding
  684. cmpl $0x0ff,%ebx
  685. jna epdp11
  686. orb $0x40,%bl
  687. epdp11: call roundi
  688. adcl $0,%eax
  689. adcl $0,%edx
  690. ret
  692. epdp20: negw %cx #this is zero or denormal
  693. jz epdp7
  694. cmpw $53,%cx
  695. jc epdp21
  696. movw $53,%cx
  697. epdp21: shrl $1,%edx
  698. rcrl $1,%eax
  699. jnc epdp22
  700. orl $1,%eax
  701. epdp22: decw %cx
  702. jnz epdp21
  704. epdp7: xorl %ebx,%ebx #discarded bits into ebx
  705. shrdl $12,%eax,%ebx
  707. shrdl $12,%edx,%eax #shift right everything 12 bits
  708. shrdl $11,%ecx,%edx
  709. btl $31,%ecx
  710. rcrl $1,%edx
  712. testb $underf,h_ctrl #set underflow if necessary
  713. jz epdp13
  714. orl %ebx,%ebx
  715. jz epdp19
  716. epdp13: orb $underf,h_stat
  717. call exproc
  718. jmp epdp8
  720. epdp30: movl %edx,%ebx #this is INF or NaN
  721. addl %ebx,%ebx
  722. orl %eax,%ebx
  723. jnz epdp34
  724. movw $0x7ff,%cx
  725. jmp epdp6
  727. epdp38: call getqn
  728. orb $stacku,h_stat
  729. jmp epdp39
  730. epdp36: call getqn #this is NaN
  731. jmp epdp37
  732. epdp34: btl $30,%edx
  733. jc epdp35
  734. epdp37: orb $invop,h_stat
  735. epdp39: call exproc
  736. btsl $30,%edx
  737. epdp35: movw $2*lBIAS+1,%cx
  738. andw $0xf800,%ax
  739. jmp epdp6
  741. # round extended floating-point to whole floating-point
  742. # input:        EP in ecx:edx:eax
  743. # output:       EP in ecx:edx:eax no decimals
  745. .globl eptown
  747. epwn2: call chkarg #check for stack underflow, NaN
  748. jc epwn82
  750. btl $16,%ecx #check zero, denormal
  751. jc epwn82
  752. orw %cx,%cx
  753. jz epwn85
  754. cmpw $1,%cx
  755. jz epwn43
  756. orl %edx,%edx
  757. jns stdnan
  759. epwn82: ret #INF
  761. epwn85: orw $denorm+precis,h_stat
  762. call exproc
  763. epwn43: orb $1,%bl
  765. epwn40: movb $0,%bl #here answer is 0 or 1 or -1
  766. jnz epwn41
  767. shldl $8,%edx,%ebx
  768. addl %edx,%edx
  769. epwn41: orl %eax,%edx
  770. jz epwn42
  771. orb $0x40,%bl
  772. epwn42: xorl %edx,%edx
  773. xorl %eax,%eax
  774. andl $0x80000000,%ecx
  775. movw $0x3fff,%cx
  776. call roundi
  777. jc epwn16
  778. addl $0xc001,%ecx
  779. ret
  781. eptown:
  782. testl $0x00030000,%ecx #jump if special value
  783. jnz epwn2
  785. epwn3: movl %ecx,%esi #number of bits to clear
  786. movw $tBIAS+63,%si
  787. subw %cx,%si
  788. jle epwn8
  790. cmpw $64,%si #over not needed
  791. jnc epwn40
  793. xchgl %esi,%ecx
  794. xorl %ebx,%ebx
  796. cmpb $32,%cl #clear decimal bits
  797. jc epwn21
  798. xchgl %ebx,%eax
  799. xchgl %eax,%edx
  800. jz epwn22
  801. cmpl %ebx,%edx
  802. rcrl $1,%ebx
  803. epwn21: shrdl %cl,%eax,%ebx
  804. shrdl %cl,%edx,%eax
  805. shrl %cl,%edx
  806. epwn22: roll $8,%ebx
  807. cmpl $0x100,%ebx
  808. jc epwn23
  809. orb $0x20,%bl
  810. epwn23: call roundi
  811. adcl $0,%eax
  812. adcl $0,%edx
  813. cmpb $32,%cl
  814. jc epwn24
  815. xchgl %eax,%edx
  816. shrl $1,%eax
  817. adcl $0,%esi
  818. epwn24: shldl %cl,%eax,%edx
  819. adcl $0,%esi
  820. shll %cl,%eax
  822. xchgl %esi,%ecx
  824. epwn16: btsl $31,%edx #set high mantissa bit
  825. epwn8: ret
  827. # extract exponent and mantissa from extended floating-point
  828. # input:        EP in ecx:edx:eax
  829. # output:       mantissa as EP in ecx:edx:eax
  830. #               [edi] exponent as EP
  832. .globl xtract
  834. xtr2: btl $16,%ecx #check zero, denormal
  835. jc xtr80
  836. orw %cx,%cx
  837. jz xtr23
  838. orl %edx,%edx
  839. jns xtr92
  840. jmp xtr90
  842. xtr23: call tnormal #normalize
  843. jmp xtr3
  845. xtract:
  846. movl %ecx,%ebx #save original sign
  847. andl $0x80000000,%ebx
  849. testl $0x00030000,%ecx #jump if special value
  850. jnz xtr2
  852. xtr3: pushl %edx #save the argument
  853. pushl %eax
  855. movswl %cx,%edx #convert real exponent to EP floating-point
  856. subl $tBIAS,%edx
  857. call wiep
  858. orw %cx,%cx
  859. jnz xtr7
  860. orl %ebx,%ecx
  861. xtr7: #PUTEP   
  862. movl %eax,(%edi)
  863. movl %edx,4(%edi)
  864. movl %ecx,8(%edi)
  866. popl %eax #restore original argument
  867. popl %edx
  869. movl $0x3fff,%ecx #make exponent 1
  870. orl %ebx,%ecx
  871. ret
  873. xtr80: orb $zerod,h_stat #argument was 0, store -INF
  874. call exproc
  875. movl %eax,(%edi)
  876. movl $0x80000000,4(%edi)
  877. movl $0x80007fff,8(%edi)
  878. ret
  880. xtr92: call getqn #invalid
  881. jmp xtr97
  883. xtr90: movl %edx,%ebx #either NaN or INF
  884. addl %ebx,%ebx
  885. orl %eax,%ebx
  886. jnz xtr95
  888. movl %eax,(%edi) #argument was INF, store INF
  889. movl $0x7fff,8(%edi)
  890. movl $0x80000000,%edx
  891. movl %edx,4(%edi)
  892. ret
  894. xtr95: btl $30,%edx #NaN, make into quiet
  895. jc xtr96
  896. xtr97: orb $invop,h_stat
  897. call exproc
  898. btsl $30,%edx
  899. xtr96: #PUTEP   
  900. movl %eax,(%edi)
  901. movl %edx,4(%edi)
  902. movl %ecx,8(%edi)
  903. ret
  905. # add integer to exponent or EP
  906. # input:        integer EP in ecx:edx:eax
  907. #               [edi] EP
  908. # output:       EP scaled in ecx:edx:eax
  910. .globl scale
  912. scaw: call retnan #remove NaNs and invalids
  913. movb $0,%bl
  915. movl %ecx,%esi #set denormal exception if either
  916. addl %esi,%esi #   argument denormal
  917. cmpl $0x00040000,%esi
  918. jz sca24
  919. movl 8(%edi),%esi
  920. addl %esi,%esi
  921. cmpl $0x00040000,%esi
  922. jnz sca26
  923. movb $1,8(%edi)
  925. sca24: orb $denorm,h_stat
  926. call exproc
  928. sca26: testb $0x01,10(%edi) #cases x = 0
  929. jz sca28
  930. cmpl $0x00027fff,%ecx
  931. jz stdnan
  933. scazer: xorl %edx,%edx #return zero 
  934. xorl %eax,%eax
  935. movl $0x00010000,%ecx
  936. jmp sca17
  938. sca28: cmpw $0x7fff,8(%edi) #cases x = INF
  939. jnz sca30
  940. cmpl $0x80027fff,%ecx
  941. jz stdnan
  943. #GETEP                   #return x
  944. movl (%edi),%eax
  945. movl 4(%edi),%edx
  946. movl 8(%edi),%ecx
  947. ret
  949. sca30: cmpl $0x80027fff,%ecx #cases n = INF
  950. jz scazer
  951. cmpl $0x00027fff,%ecx
  952. jz sca17
  953. jmp sca4
  955. sca40: movl $-2,%ebx #this is special cases for integer
  956. jmp sca9
  957. sca44: xorl %ebx,%ebx
  958. jmp sca9
  960. scale:
  961. andb $0x0fd,h_stat+1
  963. shldl $16,%ecx,%ebx #weed out special arguments
  964. orb 10(%edi),%bl
  965. jnz scaw
  967. sca4: movl %ecx,%esi #save the sign
  969. subw $tBIAS+15,%cx #calculate shift count
  970. ja sca40
  971. negw %cx
  972. cmpw $16,%cx
  973. jnc sca44
  974. addb $16,%cl
  976. shrl %cl,%edx #integer -> bx
  977. movl %edx,%ebx
  979. sca9: #GETEP                   #load the EP
  980. movl (%edi),%eax
  981. movl 4(%edi),%edx
  982. movl 8(%edi),%ecx
  984. andl $0xffff,%ebx #we need 32-bit integers
  985. andl $0xffff,%ecx
  987. sarl $31,%esi #make twos complement of integer
  988. xorl %esi,%ebx
  989. subl %esi,%ebx
  991. addl %ebx,%ecx #add integer to exponent
  992. jle sca70
  993. cmpl $0x7fff,%ecx
  994. jge sca74
  996. orl %edx,%edx #also check for denormal
  997. jns sca78
  999. movb $0,%bl #no extra bits here
  1001. sca17: roll $16,%ecx #put sign in place
  1002. movb 11(%edi),%ch
  1003. roll $16,%ecx
  1004. ret
  1006. sca70: cmpl $0xffff8000,%ecx #renormalization needed
  1007. jge sca71
  1008. movw $0x8000,%cx
  1009. sca71: clc
  1010. jmp sca78
  1011. sca74: stc
  1012. sca78: movb $0,%bl
  1013. call renorm
  1014. jmp sca17
  1016. itsstu: orb $stacku,h_stat
  1017. jmp nan3
  1018. itsnan: btl $30,%edx
  1019. jc nan9
  1020. stdnan: orb $invop,h_stat
  1021. nan3: call getqn
  1022. call exproc
  1023. btsl $30,%edx
  1024. nan9: xorl %ebx,%ebx
  1025. ret
  1027. #        subttl  CCALL Arithmetic Operations
  1028. #        page
  1030. # examine an EP operand
  1031. # input:        ecx:edx:eax
  1032. # output:       condition bits in al
  1034. .globl epexam
  1036. epexam:
  1037. xorl %ebx,%ebx #first set sign bit
  1038. shldl $2,%ecx,%ebx
  1040. testl $0x00030000,%ecx #check for specials
  1041. jnz exa10
  1043. exa8: orb $4,%bl #regular bit
  1045. exa9: movb %bl,%al #return bits
  1046. ret
  1048. exa10: btl $16,%ecx #check for zero, NaN, unsupported
  1049. jc exa12
  1050. orw %cx,%cx
  1051. jz exa14
  1052. orl %edx,%edx
  1053. jns exa9
  1055. exa20: orb $0x01,%bl #INF or NaN
  1056. addl %edx,%edx
  1057. orl %eax,%edx
  1058. jnz exa9
  1059. jmp exa8
  1061. exa14: orb $0x4,%bl #denormal case
  1062. exa12: orb $0x40,%bl #zero case
  1063. jmp exa9
  1065. # compare two EP operands
  1066. # input:        cx:edx:eax
  1067. #               in epcomu quiet NaN will not cause exception
  1068. # output:       ST - input -> condition bits in al
  1070. .globl epcomp, epcomu, xcmp
  1072. epcomu:
  1073. movb $1,%bl #unordered entry
  1074. jmp cmp1
  1076. epcomp:
  1077. movb $0,%bl #Intel entry
  1079. cmp1: movl 8(%edi),%esi #sign, exponent of B
  1081. shldl $16,%ecx,%ebx #weed out special arguments
  1082. orb 10(%edi),%bl
  1083. jnz cmpw
  1085. cmp7: movl %ecx,%ebx #set bit in bh if both negative
  1086. andl %esi,%ebx
  1087. shrl $23,%ebx
  1089. subl (%edi),%eax #compare hook line and sinker
  1090. sbbl 4(%edi),%edx
  1091. sbbl %esi,%ecx
  1092. jz cmp10
  1093. setnl %al
  1095. cmp8: xorb %bh,%al #if both negative, complement answer
  1097. cmp9: andb $0x0b8,h_stat+1
  1098. orb %al,h_stat+1
  1099. ret
  1101. cmp10: orl %eax,%edx #possibly equal
  1102. movb $1,%al
  1103. jnz cmp8
  1104. cmpeq: movb $0x40,%al
  1105. jmp cmp9
  1107. cmpw: andb $0x0b8,h_stat+1
  1108. orb $0x45,h_stat+1
  1109. call retnan
  1111. xcmp: shldl $16,%ecx,%ebx #two zeroes are equal regardless of sign
  1112. andb 10(%edi),%bl
  1113. testb $1,%bl
  1114. jnz cmpeq
  1116. movl %ecx,%ebx #set denormal exception if either
  1117. addl %ebx,%ebx #   argument denormal
  1118. cmpl $0x00040000,%ebx
  1119. jz cmpw2
  1120. movl %esi,%ebx
  1121. addl %ebx,%ebx
  1122. cmpl $0x00040000,%ebx
  1123. jnz cmpw4
  1124. cmpw2: orb $denorm,h_stat
  1125. testb $denorm,h_ctrl
  1126. jnz cmpw4
  1127. orw $0x8080,h_stat
  1128. movl $exret,-4(%ebp)
  1130. cmpw4: orw %cx,%cx #adjust exponent of pseudo-denormal
  1131. jnz cmpw5
  1132. orl %edx,%edx
  1133. jns cmpw5
  1134. incl %ecx
  1135. cmpw5: orw %si,%si
  1136. jnz cmpw6
  1137. testb $0x80,7(%edi)
  1138. jz cmpw6
  1139. incl %esi
  1141. cmpw6: andl $0x8000ffff,%ecx
  1142. andl $0x8000ffff,%esi
  1143. jmp cmp7
  1145. # add/subtract from registers to stack
  1146. # input:        ecx:edx:eax = EP number
  1147. # result:       registers = ST + number
  1149. .globl epadd,epadd1,epsub,epsubr
  1151. epsub:
  1152. shldl $16,%ecx,%ebx #weed out special arguments
  1153. orb 10(%edi),%bl
  1154. jnz sub2
  1155. movl 8(%edi),%esi
  1156. xorl $0x80000000,%esi
  1157. jmp epadd2
  1159. sub2: call retnan #handle special arguments
  1160. xorb $0x80,11(%edi)
  1161. jmp addw
  1163. epsubr:
  1164. shldl $16,%ecx,%ebx #weed out special arguments
  1165. orb 10(%edi),%bl
  1166. jnz subr2
  1167. xorl $0x80000000,%ecx
  1168. jmp epadd1
  1170. subr2: call retnan #handle special arguments
  1171. xorl $0x80000000,%ecx
  1173. addw: testb $0x02,%bl #zero is handled normally
  1174. jz epadd1
  1176. call retnan #first remove NaNs and invalids
  1178. cmpw $0x7fff,%cx #weed out A = INF
  1180. jnz addw3
  1181. cmpw %cx,8(%edi)
  1182. jnz addrta
  1183. movl %ecx,%esi
  1184. xorl 8(%edi),%esi
  1185. jns addrta
  1186. jmp stdnan
  1188. addrtb: #GETEP                   #return B or A as such
  1189. movl (%edi),%eax
  1190. movl 4(%edi),%edx
  1191. movl 8(%edi),%ecx
  1192. addrta: movb $0,%bl
  1193. ret
  1195. addw3: cmpw $0x7fff,8(%edi) #weed out B = INF
  1196. jz addrtb
  1198. call tnormal #normalize arguments
  1199. xchgl %eax,(%edi)
  1200. xchgl %edx,4(%edi)
  1201. xchgl %ecx,8(%edi)
  1202. call tnormal
  1203. jmp epadd1
  1205. epadd:
  1206. shldl $16,%ecx,%ebx #weed out special arguments
  1207. orb 10(%edi),%bl
  1208. jnz addw
  1210. epadd1: movl 8(%edi),%esi #the other exponent
  1212. epadd2: cmpw %si,%cx #if exponent of A > exponent of B,
  1213. jl add7 #   we swap arguments
  1214. movl $0,%ebx
  1215. je add34
  1216. xchgl %eax,(%edi)
  1217. xchgl %edx,4(%edi)
  1218. xchgl %ecx,%esi
  1219. add7: subw %si,%cx
  1220. negw %cx
  1222. call shiftr #shift mantissa right by cl
  1224. movw %si,%cx #load the exponent
  1226. add34: andl $0x8000ffff,%ecx #clear special bits
  1228. xorl %ecx,%esi #go subtract if different signs
  1229. js add20
  1231. addl (%edi),%eax #add mantissas
  1232. adcl 4(%edi),%edx
  1233. jnc add15
  1234. rcrl $1,%edx
  1235. rcrl $1,%eax
  1236. rcrb $1,%bl
  1237. incl %ecx
  1238. cmpw $0x7fff,%cx
  1239. je add19
  1241. add15: orw %cx,%cx #renormalize if needed
  1242. jle add28
  1243. ret
  1245. add20: sbbl (%edi),%eax #subtract mantissas
  1246. sbbl 4(%edi),%edx
  1247. jnc add24
  1248. xorl $0x80000000,%ecx
  1249. notl %edx
  1250. notl %eax
  1251. negb %bl
  1252. sbbl $-1,%eax
  1253. sbbl $-1,%edx
  1255. add24: jns add25 #weed out not-normal, possible all-zero
  1256. orw %cx,%cx
  1257. jle add19
  1258. ret
  1260. add25: je add27 #here needs normalizing, or is zero             
  1261. add26: decw %cx
  1262. addb %bl,%bl
  1263. adcl %eax,%eax
  1264. adcl %edx,%edx
  1265. jns add26
  1266. orw %cx,%cx
  1267. jle add19
  1268. ret
  1270. add28: orl %edx,%edx #this may be 0+0
  1271. jnz add19
  1272. movb %bl,%dl
  1273. orl %eax,%edx
  1274. jnz add22
  1275. btsl $16,%ecx
  1276. ret
  1278. add27: orb %bl,%dl #here may be all zero
  1279. orl %eax,%edx
  1280. jz addzer
  1281. add22: xorl %edx,%edx
  1283. add19: call renorm #normalize
  1284. ret
  1286. addzer: movl $0x00010000,%ecx #result went all the way to zero
  1287. movb h_ctrl+1,%bh
  1288. andb $0x0c,%bh
  1289. cmpb $0x04,%bh
  1290. jnz add18
  1291. btsl $31,%ecx
  1292. add18: ret
  1294. # multiply registers to stack
  1295. # input:        cx:edx:eax = EP number
  1296. # result:       registers = ST * number
  1298. .globl epmul,epmul1
  1300. mulzer: shrl $16,%ecx
  1301. movb $1,%cl
  1302. shll $16,%ecx
  1303. xorl %edx,%edx
  1304. jmp mulin2
  1306. mulinf: orl $0x00027fff,%ecx #return INF
  1307. movl $0x80000000,%edx
  1308. mulin2: xorl %eax,%eax
  1309. xorb %bl,%bl
  1310. jmp mul21
  1312. mulw: call retnan #weed out NaN arguments
  1314. movl %ecx,%esi #weed out special cases
  1315. addl 8(%edi),%esi
  1316. btrl $31,%esi
  1317. cmpl $0x37fff,%esi
  1318. jz stdnan
  1319. testb $1,%bl
  1320. jnz mulzer
  1321. cmpw $0x7fff,%si
  1322. jnc mulinf
  1324. call tnormal #normalize arguments if needed
  1325. xchgl %eax,(%edi)
  1326. xchgl %edx,4(%edi)
  1327. xchgl %ecx,8(%edi)
  1328. call tnormal
  1329. jmp epmul1
  1331. epmul:
  1332. shldl $16,%ecx,%ebx #weed out special arguments
  1333. orb 10(%edi),%bl
  1334. jnz mulw
  1336. epmul1: pushl %ebp #save exponent, sign, FP stack pointer
  1337. pushl %ecx
  1339. movl %eax,%esi #save A0123
  1340. movl %edx,%ecx
  1342. mull (%edi) # A23 * B23
  1343. xchgl %eax,%esi
  1344. movl %edx,%ebx
  1346. mull 4(%edi) # A23 * B01
  1347. addl %eax,%ebx
  1348. adcl $0,%edx
  1349. movl %edx,%ebp
  1351. movl (%edi),%eax # A01 * B23
  1352. mull %ecx
  1353. addl %eax,%ebx
  1354. adcl %edx,%ebp
  1355. movl $0,%eax
  1356. adcl $0,%eax
  1357. xchgl %ecx,%eax
  1359. mull 4(%edi) # A01 * B01
  1360. addl %ebp,%eax
  1361. adcl %ecx,%edx
  1363. popl %ecx #exponent
  1365. js mul42 #possibly normalize
  1366. addl %ebx,%ebx
  1367. adcl %eax,%eax
  1368. adcl %edx,%edx
  1369. decl %ecx
  1371. mul42: addl $-1,%esi #set rounding bits
  1372. adcl %ebx,%ebx
  1373. jz mul43
  1374. movb $0x80,%bl
  1375. mul43: rcrb $1,%bl
  1377. popl %ebp #restore base pointer
  1379. subw $tBIAS-2,%cx #calculate new exponent
  1380. addw 8(%edi),%cx
  1381. jo mul44
  1382. decw %cx
  1383. jle mul45
  1385. mul21: roll $8,%ecx #put in sign
  1386. xorb 11(%edi),%cl
  1387. rorl $8,%ecx
  1388. ret #return
  1390. mul44: stc #here re-normalize
  1391. jmp mul46
  1392. mul45: clc
  1393. mul46: call renorm
  1394. jmp mul21
  1396. # divide registers to stack
  1397. # input:        cx:edx:eax = EP number
  1398. # result:       registers = ST / number
  1400. .globl epdiv,epdivr,epdiv1
  1402. divw61: orb $zerod,h_stat
  1403. call exproc
  1404. jmp mulinf
  1406. divw: call retnan #weed out NaN arguments
  1408. cmpw $0x7fff,%cx #weed out A = INF
  1409. jnz divw3
  1410. cmpw %cx,8(%edi)
  1411. jz stdnan
  1412. jmp mulinf
  1414. divw3: cmpw $0x7fff,8(%edi) #weed out B = INF
  1415. jz mulzer
  1417. btl $16,%ecx #if A = 0 and B = 0 return NaN
  1418. jnc divw7
  1419. testb $1,10(%edi)
  1420. jnz stdnan
  1421. jmp mulzer
  1423. divw7: testb $1,10(%edi) #if B = 0 "zero divide"
  1424. jnz divw61
  1426. call tnormal #normalize arguments
  1427. xchgl %eax,(%edi)
  1428. xchgl %edx,4(%edi)
  1429. xchgl %ecx,8(%edi)
  1430. call tnormal
  1431. xchgl %eax,(%edi)
  1432. xchgl %edx,4(%edi)
  1433. xchgl %ecx,8(%edi)
  1434. jmp epdiv1
  1436. epdivr:
  1437. xchgl (%edi),%eax
  1438. xchgl 4(%edi),%edx
  1439. xchgl 8(%edi),%ecx
  1441. epdiv:
  1442. shldl $16,%ecx,%ebx #weed out special arguments
  1443. orb 10(%edi),%bl
  1444. jnz divw
  1446. epdiv1: pushl %ebp #save exponent, sign
  1447. pushl %ecx
  1449. movl $1,%ebp #first bit = 0
  1451. cmpl 4(%edi),%edx #jump if x < y
  1452. jc div50
  1454. subl 0(%edi),%eax #get first bit
  1455. sbbl 4(%edi),%edx
  1456. jnc div49
  1458. movl %eax,%edx #prevent overflow
  1459. xorl %eax,%eax
  1460. movl %eax,%ecx
  1461. jmp div31
  1463. div49: decl %ebp #first bit = 1
  1465. div50: divl 4(%edi) # q = A0123 / B01 -> ebx
  1466. movl %eax,%ecx #     r -> esi
  1467. movl %edx,%esi
  1469. mull 0(%edi) # x * B23 -> edx:eax
  1471. xchgl %edx,%esi #subtract that from qA23
  1472. negl %eax
  1473. sbbl %esi,%edx
  1474. jnc div41
  1475. div31: decl %ecx
  1476. addl 0(%edi),%eax
  1477. adcl 4(%edi),%edx
  1478. jnc div31
  1480. div41: cmpl 4(%edi),%edx #prevent overflow
  1481. jnz div51
  1482. movl %eax,%esi
  1483. xorl %eax,%eax
  1484. movl %eax,%ebx
  1485. subl 0(%edi),%esi
  1486. jmp div32
  1488. div51: divl 4(%edi) # x = A0123 / B01 -> bx
  1489. movl %eax,%ebx
  1490. movl %edx,%esi
  1492. mull 0(%edi) #then q * B23 -> edx:eax
  1494. negl %eax #subtract that from qA23
  1495. sbbl %edx,%esi
  1496. jnc div42
  1497. div32: decl %ebx
  1498. addl 0(%edi),%eax
  1499. adcl 4(%edi),%esi
  1500. jnc div32
  1502. div42: movl %esi,%edx #set rounding bit, first bit flag
  1503. orl %eax,%edx
  1504. xchgl %ebx,%edx
  1505. jz div54
  1506. movb $0x80,%bl
  1507. addl %eax,%eax
  1508. adcl %esi,%esi
  1509. jc div43
  1510. subl 0(%edi),%eax
  1511. sbbl 4(%edi),%esi
  1512. cmc
  1513. div43: rcrb $1,%bl
  1515. div54: movl %edx,%eax #load result to edx:eax
  1516. movl %ecx,%edx
  1518. popl %ecx #restore registers
  1520. cmpl $1,%ebp #handle first bit
  1521. jnc div57
  1522. rcrl $1,%edx
  1523. rcrl $1,%eax
  1524. rcrb $1,%bl
  1525. incw %cx
  1527. div57: popl %ebp
  1529. subw 8(%edi),%cx #calculate new exponent
  1530. addw $tBIAS-0,%cx
  1531. jo div23
  1532. decw %cx
  1533. jle div24
  1535. div21: roll $8,%ecx #put in sign
  1536. xorb 11(%edi),%cl
  1537. rorl $8,%ecx
  1538. ret
  1540. div23: stc #renormalize if needed
  1541. jmp div26
  1542. div24: clc
  1543. div26: call renorm
  1544. jmp div21
  1546. # partial remainder, Intel version
  1547. # input:        cx:edx:eax = EP number
  1548. # result:       registers = ST % number
  1550. .globl eprem
  1552. res88: cmpw $-1,%cx #here exponent of A < exponent of B
  1553. jnz rem83 #   return A unless A > B/2
  1554. pushl %eax #   in which case return A-B
  1555. pushl %edx
  1556. shll $1,(%edi)
  1557. rcll $1,4(%edi)
  1558. btrl $31,%edx
  1559. jmp res17
  1561. rem87: btl $16,%ebx #branch if IEEE remainder
  1562. jc res88
  1563. rem83: jmp rem16
  1565. remw: andb $0x0fb,h_stat+1 #clear C2
  1566. call retnan #remove NaNs and invalids
  1568. andb $0x0b8,h_stat+1 #clear condition codes
  1570. btl $16,%ecx #if B = 0 or A = INF return NaN
  1571. jc stdnan
  1572. cmpw $0x7fff,%si
  1573. jz stdnan
  1575. call tnormal #normalize arguments
  1576. xchgl %eax,(%edi)
  1577. xchgl %edx,4(%edi)
  1578. xchgl %ecx,%esi
  1579. cmpw $0x7fff,%si
  1580. jz rem17
  1581. btl $16,%ecx
  1582. jc rem19
  1583. call tnormal
  1584. jmp rem3
  1586. eprem:
  1587. movl 8(%edi),%esi
  1589. shldl $16,%ecx,%ebx #weed out special arguments
  1590. orb 10(%edi),%bl
  1591. jnz remw
  1593. andb $0x0b8,h_stat+1 #clear condition codes
  1595. xchgl (%edi),%eax #swap registers
  1596. xchgl 4(%edi),%edx
  1597. xchgl %esi,%ecx
  1599. rem3: movb $0,%bl #initialize quotient
  1601. subw %si,%cx #get the shift count
  1602. js rem87
  1604. cmpw $64,%cx #maximum shift count 64
  1605. jl rem4
  1606. addw %cx,%si
  1607. orb $32,%cl
  1608. andw $63,%cx
  1609. subw %cx,%si
  1610. orb $0x04,h_stat+1
  1611. rem4: movw %si,8(%edi)
  1613. pushl %ebp #save registers
  1614. movl (%edi),%ebp
  1615. movl 4(%edi),%esi
  1617. rem5: subl %ebp,%eax #divide using shift-subtract
  1618. sbbl %esi,%edx
  1619. jnc rem9
  1620. addl %ebp,%eax
  1621. adcl %esi,%edx
  1622. rem9: cmc
  1623. rem11: adcb %bl,%bl
  1624. decb %cl
  1625. js rem15
  1626. addl %eax,%eax
  1627. adcl %edx,%edx
  1628. jnc rem5
  1629. subl %ebp,%eax
  1630. sbbl %esi,%edx
  1631. jmp rem11
  1633. rem15: popl %ebp #restore registers
  1635. testb $0x04,h_stat+1 #if C2 is set (incomplete)
  1636. jnz rem16 #   leave Q zero
  1638. testl $0x10000,%ebx #branch if IEEE remainder
  1639. jnz res18
  1641. rem41: andl $7,%ebx #massage quotient for the required form
  1642. rorw $1,%bx #   bits are discontinuous and reversed
  1643. shrb $4,%bh
  1644. rorw $1,%bx
  1645. shrb $2,%bh
  1646. shlb $7,%bl
  1647. shll $1,%ebx
  1648. orb %bh,h_stat+1
  1650. rem16: movw 8(%edi),%cx #normalize
  1651. rem17: xorb %bl,%bl
  1652. call renorm
  1653. rem19: ret
  1655. res18: pushl %eax #save remainder
  1656. pushl %edx
  1658. res17: addl %eax,%eax #if next quotent bit would be 0, we are
  1659. adcl %edx,%edx #   done
  1660. setc %bh
  1661. subl (%edi),%eax
  1662. sbbl 4(%edi),%edx
  1663. sbbb $0,%bh
  1664. jc res42
  1666. movl %ebx,%esi #next bit is one, if all the rest are 0,
  1667. andl $1,%esi #   and Q is even, we are likewise done
  1668. orl %edx,%esi
  1669. orl %eax,%esi
  1670. jz res42
  1672. incl %ebx #here we round up the quotient
  1673. popl %edx
  1674. popl %eax
  1675. subl (%edi),%eax
  1676. sbbl 4(%edi),%edx
  1677. notl %edx
  1678. negl %eax
  1679. sbbl $-1,%edx
  1681. btcl $31,%ecx
  1682. jmp rem41
  1684. res42: popl %edx
  1685. popl %eax
  1686. jmp rem41
  1688. # square-root
  1689. # input:        ecx:edx:eax = EP number
  1690. # result:       ST = sqrt(number)
  1692. .globl epsqrt
  1694. sqr89: ret
  1695. sqrw: xorl %ebx,%ebx
  1696. btl $16,%ecx #zero returns as such
  1697. jc sqr89
  1699. orw %cx,%cx #check for unsupported
  1700. jz sqrw5
  1701. orl %edx,%edx
  1702. jns stdnan
  1704. cmpw $0x7fff,%cx #check for NaN
  1705. jnz sqrw5
  1706. movl %edx,%ebx
  1707. addl %ebx,%ebx
  1708. orl %eax,%ebx
  1709. jnz nanret
  1711. sqrw5: orl %ecx,%ecx #negative is illegal
  1712. js stdnan
  1714. cmpw $0x7fff,%cx #+INF returns as such
  1715. jz sqr89
  1717. call tnormal #normalize
  1718. movl %ecx,8(%edi)
  1719. jmp sqr2
  1721. epsqrt:
  1722. testl $0x80030000,%ecx #jump if special value
  1723. jnz sqrw
  1725. sqr2:
  1726. pushl %ebp #save registers
  1727. pushl %edi
  1729. xorl %ebx,%ebx #save x, possibly shifted right by 1
  1730. movl %eax,%esi
  1731. andb $1,%cl
  1732. shrdl %cl,%esi,%ebx
  1733. shrdl %cl,%edx,%esi
  1734. pushl %esi
  1735. pushl %ebx
  1737. incl %ecx #load high 30 or 31 bits -> esi:edi
  1738. shrdl %cl,%edx,%eax
  1739. shrl %cl,%edx
  1740. movl %edx,%esi
  1741. movl %eax,%edi
  1743. testb $1,%cl #get first estimate with a*x+b
  1744. movl $0x4d12cfe4,%ebx
  1745. jz sqr3
  1746. movl $0x6cff9300,%ebx
  1747. sqr3: shldl $17,%esi,%eax
  1748. mulw %bx
  1749. shrl $16,%ebx
  1750. addw %dx,%bx
  1752. movl %esi,%eax #16-bit iteration y = (x/y + y)/2
  1753. shldl $16,%esi,%edx
  1754. divw %bx
  1755. shrw $1,%bx
  1756. addw %ax,%bx
  1757. sbbw $0,%bx
  1758. shll $16,%ebx
  1760. movl %esi,%edx #32-bit iteration y = (x/y + y)/2
  1761. movl %edi,%eax
  1762. divl %ebx
  1763. shrl $1,%ebx
  1764. addl %eax,%ebx
  1765. sbbl $0,%ebx
  1767. movl %esi,%edx #64-bit iteration y = (x/y + y)/2
  1768. movl %edi,%eax
  1769. divl %ebx
  1770. movl %eax,%ecx
  1771. xorl %eax,%eax
  1772. divl %ebx
  1773. xorl %edx,%edx
  1774. shrl $1,%ebx
  1775. rcrl $1,%edx
  1776. addl %edx,%eax
  1777. adcl %ebx,%ecx
  1779. movl %eax,%ebx #mantissa of y -> ecx:ebx
  1781. mull %eax #calculate y^2 -> (none):si:edi:ebp
  1782. movl %eax,%ebp
  1783. movl %edx,%edi
  1784. movl %ecx,%eax
  1785. mull %eax
  1786. movl %eax,%esi
  1787. movl %ebx,%eax
  1788. mull %ecx
  1789. addl %eax,%edi
  1790. adcl %edx,%esi
  1791. addl %eax,%edi
  1792. adcl %edx,%esi
  1794. movl %ebx,%eax #mantissa of y -> ecx:eax
  1796. popl %ebx #now y^2 - x into si:edi:ebp
  1797. subl %ebx,%edi
  1798. popl %ebx
  1799. sbbl %ebx,%esi
  1800. js sqr37
  1801. jz sqr35
  1803. sqr31: subl $1,%eax #try out next smallest y 
  1804. sbbl $0,%ecx
  1806. stc #new residue = old - 2y + 1
  1807. sbbl %eax,%ebp
  1808. sbbl %ecx,%edi
  1809. sbbl $0,%esi
  1810. subl %eax,%ebp
  1811. sbbl %ecx,%edi
  1812. sbbl $0,%esi
  1814. jns sqr31 #more if we are still above real sqrt(x)
  1815. cmc
  1817. sqr32: sbbb %bl,%bl #1/2 bit here, never exactly half
  1818. orb $0x40,%bl
  1820. sqr33: movl %ecx,%edx
  1822. popl %edi #restore registers
  1823. popl %ebp
  1825. movl 8(%edi),%ecx #calculate new exponent
  1826. shrw $1,%cx
  1827. adcw $0x1fff,%cx
  1828. ret
  1830. sqr37: stc #try next largest x
  1831. adcl %eax,%ebp
  1832. adcl %ecx,%edi
  1833. adcl $0,%esi
  1834. addl %eax,%ebp
  1835. adcl %ecx,%edi
  1836. adcl $0,%esi
  1838. jns sqr32 #quit if differential > 0
  1839. addl $1,%eax
  1840. adcl $0,%ecx
  1841. jmp sqr37
  1843. sqr35: movl %edi,%ebx #see if we have an exact answer
  1844. orl %ebp,%ebx
  1845. jz sqr33
  1846. jmp sqr31
  1848. #__lib   endp
  1850. #emuseg  ends
  1852. #        end