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

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

VxWorks

开发平台：

C/C++

uss_fpfncs.s：源码内容
  1. /* Copyright 1991-1992 Wind River Systems, Inc. */
  2. .data
  3. .globl _copyright_wind_river
  4. .long _copyright_wind_river
  6. /*
  7. modification history
  8. --------------------
  9. 01f,23aug92,jcf  changed bxxx to jxx.
  10. 01e,05jun92,kdl  changed external branches to jumps (and bsr's to jsr's).
  11. 01d,26may92,rrr  the tree shuffle
  12. 01c,30mar92,kdl  added include of "uss_fp.h"; commented-out ".set" directives
  13.  (SPR #1398).
  14. 01b,04oct91,rrr  passed through the ansification filter
  15.   -changed ASMLANGUAGE to _ASMLANGUAGE
  16.   -changed copyright notice
  17. 01a,28jan91,kdl  original US Software version.
  18. */
  21. /*
  22. DESCRIPTION
  24. |       ttl     FPAC 68K/FPFNCS: IEEE Single Precision Functions
  25. |FPFNCS idnt    1,0             ; IEEE Single Precision Functions
  26. |                               ; FPFNCS.A68
  27. |
  28. |* * * * * * * * * *
  29. |
  30. |       68000 FPAC  --  Single Precision Functions
  31. |
  32. |       Copyright (C) 1984 - 1990 By
  33. |       United States Software Corporation
  34. |       14215 N.W. Science Park Drive
  35. |       Portland, Oregon 97229
  36. |
  37. |       This software is furnished under a license and may be used
  38. |       and copied only in accordance with the terms of such license
  39. |       and with the inclusion of the above copyright notice.
  40. |       This software or any other copies thereof may not be provided
  41. |       or otherwise made available to any other person.  No title to
  42. |       and ownership of the software is hereby transferred.
  43. |
  44. |       The information in this software is subject to change without
  45. |       notice and should not be construed as a commitment by United
  46. |       States Software Corporation.
  47. |
  48. |       Version:        See VSNLOG.TXT
  49. |       Released:       01 Jan 1984
  50. |
  51. |* * * * * * * * * *
  52. |
  53. |
  54. |  Single Precision Floating Point Function library.
  55. |
  56. |  These functions operate on single precision floating point values,
  57. |  yielding single precision floating point results.  All functions
  58. |  process the value on top of the stack and return a result on the
  59. |  stack.  The XTOI function expects its integer argument in D0.
  60. |
  61. |
  62. |  NAME       ARGUMENT                DESCRIPTION
  63. |  ========== =====================   =================================
  64. |  EXP        -87.43 < a < 88.72      e to the power specified.
  65. |                                     (if out of range, 0.0 or +INF)
  66. |
  67. |  LN         a > 0.0                 natural log of value specified.
  68. |                                     (invalid opn if neg or zero)
  69. |
  70. |  XTOI       all value pairs         number raised to integer power.
  71. |
  72. |  ATN        any numeric value       arc whose tanget is given value
  73. |                                     (no domain restriction, result in
  74. |                                     the range is -PI/2 to +PI/2)
  75. |
  76. |  SQRT       a >= 0.0                square root of the given value
  77. |                                     (invalid operation if negative)
  78. |
  79. |  COS                                All trigonometric functions
  80. |  SIN                                expect a radian mode argument
  81. |  TAN                                whose magnitude is < 65536 (2^16)
  82. |
  84. NOMANUAL
  85. */
  87. #define _ASMLANGUAGE
  89. #include "uss_fp.h"
  91. |       page
  92. |
  93. |       opt     BRS             ; Default to forward branches SHORT
  94. |
  95.         .globl  FPEXP
  96.         .globl  FPLN
  97.         .globl  FPLOG
  98.         .globl  FPXTOI
  99.         .globl  FPATN
  100.         .globl  FPSQRT
  101.         .globl  FPCOS
  102.         .globl  FPSIN
  103.         .globl  FPTAN
  104. |
  105. |       xref    FPMUL   ;Floating point operations
  106. |       xref    FPDIV
  107. |       xref    FPRDIV
  108. |       xref    FPADD
  109. |       xref    FLOAT
  110. |       xref    INT
  111. |
  112. |       xref    FNANRS
  113. |       xref    FINFRS
  114. |       xref    FUNFRS
  115. |       xref    FZERRS
  116. |
  117. |       xref    FPERR
  118. |
  119. |
  120. |       .set    FBIAS,127
  121. |
  122. |       .set    ERUNF,1 |FPERR UNDERFLOW ERROR CODE
  123. |       .set    EROVF,2 |FPERR OVERFLOW ERROR CODE
  124. |       .set    ERNAN,3 |FPERR INVALID OPERATION ERROR CODE
  125. /*
  126. |
  127. |
  128. |
  129. | ###   SUBTTL  FPEXP: Exponentiation of e Function
  130. |       page
  131. |
  132. |
  133. |  EXP:   Calculate e to the power specified by the stack argument.
  134. |         The result is returned on the stack.
  135. |
  136. |         FPERR is set upon overflow, underflow, or NaN argument.
  137. |
  138. |  Algorithm:
  139. |  The argument is multiplied by 1/LN 2.  The integer portion of the
  140. |  product becomes the two*s exponent.  A Taylor series approximation
  141. |  about zero is used to compute 2**FRAC(ARG * 1/LN 2).
  142. |
  143. */
  145.         .text
  146. |
  147. FEXPCN:
  148.         |dsw    0       | EXP CONSTANT LIST (2**X EVALUATION)
  149. |
  150. |         (LN 2)^ 9/ 9!   =   1.0178 08601E-07
  151.         .long   0x33DA929E
  152. |         (LN 2)^ 8/ 8!   =   1.3215 48679E-06
  153.         .long   0x35B16011
  154. |         (LN 2)^ 7/ 7!   =   1.5252 73380E-05
  155.         .long   0x377FE5FE
  156. |         (LN 2)^ 6/ 6!   =   1.5403 53039E-04
  157.         .long   0x39218489
  158. |         (LN 2)^ 5/ 5!   =   1.3333 55815E-03
  159.         .long   0x3AAEC3FF
  160. |         (LN 2)^ 4/ 4!   =   9.6181 29108E-03
  161.         .long   0x3C1D955B
  162. |         (LN 2)^ 3/ 3!   =   5.5504 10866E-02
  163.         .long   0x3D635847
  164. |         (LN 2)^ 2/ 2!   =   2.4022 65070E-01
  165.         .long   0x3E75FDF0
  166. |         (LN 2)^ 1/ 1!   =   6.9314 71806E-01
  167. FLN2:
  168.         .long   0x3F317218
  169. |                             1.0000 00000E+00
  170.         .long   0x3F800000
  171. |
  172. |       .set    FNEXPC,10               | Number of constants in FEXP poly
  173. |
  174. |         1 / LN 2  =  1.4426 95040 88896 34074 D+00
  175. FINLN2:
  176.         .long   0x3FB8AA3B
  177. |
  178. |
  179. |
  180. | ###   PUBLIC  FPEXP
  181. |
  182. FPEXP:
  183.         |dsw    0
  184.         bsr     FPFADJ          | Set-up for FP functions
  185. |
  186.         jvs     FFNANR          | NaN parm -> NaN result
  187.         jcc     FEXP01          | J/ arg not INF
  188. |
  189.         jmi     FFUNFR          | -INF parm -> underflow result
  190.         jra     FFPINR          | +INF parm -> positive INF result
  191. |
  192. FEXP01:
  193.         movel   FINLN2,sp@-     | 1/LN 2 on stack
  194.         jsr     FPMUL           | Convert to 2^ function
  195. |
  196.         movew   sp@,d1
  197.         rorw    #7,d1           | Exponent in d1.B
  198.         cmpib   #FBIAS+6,d1
  199.         jhi     FEXP10          | J/ overflow or underflow  >= +/- 2^7
  200.         jne     FEXP20          | J/ within range            < +/- 2^6
  201.         tstb    sp@
  202.         jpl     FEXP20          | J/ result >= 1.0
  203. |
  204.         cmpiw   #0xF800+0x0100+FBIAS+6,d1       | Check ms 7 man bits (+s+exp)
  205.         jcs     FEXP20                  | J/ result >= 2^-126
  206. |
  207. FEXP10:
  208.         tstb    sp@
  209.         jpl     FFPINR          | J/ FP function +INF result
  210.         jra     FFUNFR          | J/ FP function underflow result
  211. |
  212. FEXP20:
  213.         movel   sp@,sp@-        | Duplicate the argument on the stack
  214. |
  215.         jsr     INT
  216. |
  217.         movew   d0,a7@(4)       | Save result of INT
  218. |
  219.         negl    d0              | Negate integer
  220.         jsr     FLOAT
  221.         jsr     FPADD           | Special FFRAC
  222. |
  223.         pea     FEXPCN
  224.         movel   sp@+,a6 | Point a6 to constant list
  225.         moveq   #FNEXPC,d0      | d0 holds the number of constants
  226.         bsr     FXSER           | Polynomial series
  227. |
  228.         movew   a7@(4),d0       | Fetch power of two scaling
  229.         lslw    #7,d0
  230.         addw    d0,sp@          | Combine with FXSER result
  231. |
  232.         movel   sp@+,sp@        | Downshift result
  233.         jmp     a4@             | Return
  235. /*
  236. |
  237. | ###   SUBTTL  FPLOG, FPLN: Logarithm Functions
  238. |       page
  239. |
  240. |  LOGARITHM (BASE E) FUNCTION.
  241. |
  242. |  The natural logarithm of the single precision floating point
  243. |  value on the stack is computed using a split domain polynomial
  244. |  approximation.  The common log is computed by scaling the
  245. |  result of the natural log computation.
  246. |
  247. |  TWOS = FLOAT(ARG.EXPONENT) * LOG(2)
  248. |
  249. |  The argument, (after the two's power scaling, 2.0 > arg' >= 1.0),
  250. |  is sorted into one of five classes.  Within that class, a center
  251. |  point with a known log value is used in combination with a 7th
  252. |  degree Chebyshey polynomial approximation, to calculate a
  253. |  logarithm.
  254. |
  255. |  The center point of each class is adjusted so that its natural
  256. |  logarithm expressed in single precision form is accurate to
  257. |  more than 40 bits.
  258. |
  259. |  log(arg) = two's log + center point log + polynomial approximation
  260. |
  261. |  Polynomial approximation:
  262. |
  263. |  log((1-t)/(1+t)) = t*(c1*t^8 + c2*t^6 + c3*t^4 + c4*t^2 + c5)
  264. |         (accurate with 1.0E-17 for t within 0.1 of 1.0)
  265. |
  266. |  t = - (center - arg')/(center + arg')  [to calc log of arg'/center]
  267. |
  268. */
  270. FLNCNS:
  271.         |dsw    0               | LN constants list
  272. |
  273. |         c1  =   0.28571 20487
  274.         .long   0x3E9248DA
  275. |         c2  =   0.40000 00019
  276.         .long   0x3ECCCCCD
  277. |         c3  =   0.66666 66667
  278.         .long   0x3F2AAAAB
  279. |         c4  =   2.00000 00000
  280.         .long   0x40000000
  281. |
  282. |       .set    FNLNCN,4                | Number of constants in FLN poly
  283. |
  284. |
  285. FLNCEN:
  286.         |dsw    0               | Center Points
  287. |         Center  #1 = 1.0000 00000 00000 00000
  288.         .long   0x3F800000
  289. |         Center  #2 = 1.1892 07098     (close to 2^(1/4))
  290.         .long   0x3F9837F0
  291. |         Center  #3 = 1.4142 13521     (close to 2^(1/2))
  292.         .long   0x3FB504F3
  293. |         Center  #4 = 1.6817 92733     (close to 2^(3/4))
  294.         .long   0x3FD744FC
  295. |         Center  #5 = 2.0000 00000
  296. | ---     shifted to 1.0 center..
  297. |
  298. FLNLOG:
  299.         |dsw    0               | FLNCEN values adj for exact ln values
  300. |         LN(Center #2) = 0.17328 67807
  301.         .long   0x3E317217
  302. |         LN(Center #3) = 0.34657 35614
  303.         .long   0x3EB17217
  304. |         LN(Center #4) = 0.51986 03272
  305.         .long   0x3F051591
  306. |
  307. |
  308. |         1 / LN 10  =  0.43429 44819
  309. FILN10:
  310.         .long   0x3EDE5BD9
  311. |
  312. |
  313. |
  314. | ###   PUBLIC  FPLN
  315. | ###   PUBLIC  FPLOG
  316. |
  317. FPLOG:
  318.         |dsw    0
  319.         moveq   #-1,d0          | Signal common log scaling
  320.         jra     FLN000
  321. |
  322. FPLN:
  323.         |dsw    0
  324.         clrb    d0              | Signal no post calculation scaling
  325. |
  326. FLN000:
  327.         bsr     FPFADJ          | Adjust parameter on stack
  328. |
  329.         jvs     FFNANR          | J/ NaN arg -> NaN result
  330.         jmi     FFNANR          | J/ Neg arg -> NaN result
  331.         jeq     FFMINR          | J/ 0.0 arg -> -INF result
  332.         jcs     FFPINR          | J/ +INF arg -> +INF result
  333. |
  334.         moveb   d0,a7@(6)       | Save natural/common flag
  335. |
  336.         movew   sp@,d1          | Prepare to calc parm sign and exp
  337.         lsrw    #7,d1           | Position exponent (right justified)
  338.         subiw   #FBIAS,d1
  339.         movew   d1,a7@(4)       | /* Save two's exponent value */
  340. |
  341.         lslw    #7,d1           | Scale parameter
  342.         subw    d1,sp@
  343. |
  344.         clrw    d0              | Set class number to 0
  345.         movel   sp@,d1
  346.         lsrl    #8,d1           | d1.W = 0xEMMM
  347. |
  348.         cmpiw   #0x8000+2966,d1
  349.         jcs     FLN050          | J/ < 1 +  2966/32768  (1.090515)
  350.         addqw   #4,d0
  351.         cmpiw   #0x8000+9727,d1
  352.         jcs     FLN050          | J/ < 1 +  9727/32768  (1.296844)
  353.         addqw   #4,d0
  354.         cmpiw   #0x8000+17767,d1
  355.         jcs     FLN050          | J/ < 1 + 17767/32768  (1.542206)
  356.         addqw   #4,d0
  357.         cmpiw   #0x8000+27329,d1
  358.         jcs     FLN050          | J/ < 1 + 27329/32768  (1.834015)
  359.         addqw   #1,a7@(4)       | /* Class 4: Bump two's exp, then class 0 */
  360.         clrw    d0
  361.         subib   #0x80,a7@(1)    | Reduce scaled exponent
  362. |
  363. FLN050:
  364.         moveb   d0,a7@(7)       | Save center number
  365. |
  366.         pea     FLNCEN
  367.         movel   sp@+,a2
  368.         movel   a2@(0,d0:W),sp@- |Place center on stack
  369. |
  370.         movel   a7@(4),sp@-     | Copy scaled argument
  371. |
  372.         movel   a2@(0,d0:W),sp@- |Place center on stack
  373. |
  374.         jsr     FPADD           | /* Calc  (CENTER + PARM') */
  375. |
  376.         movel   sp@,d0          | Exch top stack item w/ 3rd
  377.         movel   a7@(8),d1
  378.         movel   d1,sp@
  379.         movel   d0,a7@(8)
  380. |
  381.         bset    #7,a7@(4)       | Negate CENTER
  382. |
  383.         jsr     FPADD           | /* Calc  - (CENTER - PARM') */
  384.         jsr     FPRDIV          | /* T = - (CENTER-PARM')/(CENTER+PARM) */
  385. |
  386.         movel   sp@,sp@-        | Duplicate t
  387. |
  388.         pea     FLNCNS
  389.         movel   sp@+,a6 | Poly approx
  390.         moveq   #FNLNCN,d0
  391.         bsr     FX2SER
  392. |
  393.         jsr     FPMUL           | Times t
  394.         clrw    d0
  395.         moveb   a7@(7),d0
  396.         jeq     FLN060          | J/ center = 0 -> no additive log val
  397. |
  398.         pea     FLNLOG-4
  399.         movel   sp@+,a6
  400.         movel   a6@(0,d0:W),sp@-
  401.         jsr     FPADD           | Add in center log value
  402. |
  403. FLN060:
  404.         movew   a7@(4),d0       | /* Get two's exponent value */
  405.         extl    d0
  406. | ##    DC.W    $48C0           ; (Assembled  EXT.L D0 instruction)
  407.         jsr     FLOAT
  408. |
  409.         movel   FLN2,sp@-       | Log of 2 on stack
  410.         jsr     FPMUL
  411.         jsr     FPADD
  412. |
  413.         tstb    a7@(6)
  414.         jeq     FLN070          | J/ natural log
  415. |
  416.         movel   FILN10,sp@-     | Scaling value
  417.         jsr     FPMUL
  418. |
  419. FLN070:
  420.         movel   sp@+,sp@
  421.         jmp     a4@
  422. /*
  423. | ###   SUBTTL  FPXTOI: Floating Point Number to Integer Power Function
  424. |       page
  425. |
  426. |  X to I power Function.
  427. |
  428. |  The single precision floating point value on the stack is
  429. |  raised to the power specified in D0.W then returned on stack.
  430. |
  431. |  A shift and multiply technique is used (possibly with a trailing
  432. |  recipication.
  433. |
  434. |
  435. | ###   PUBLIC  FPXTOI
  436. |
  437. */
  439. FPXTOI:
  440.         |dsw    0
  441.         bsr     FPFADJ
  442. |
  443.         jvs     FFNANR          | J/ NaN arg -> NaN result
  444.         jcc     FPXT10          | J/ arg is not INF
  445. |
  446.         tstw    d0
  447.         jeq     FFNANR          | J/ arg is +/- INF, I is 0 -> NaN
  448.         jmi     FFUNFR          | J/ x is +/- INF, I < 0 -> underflow
  449. |
  450.         rorb    #1,d0
  451.         andb    sp@,d0
  452.         jmi     FFMINR          | J/ arg is -INF, I is +odd -> -INF
  453.         jra     FFPINR          | else -> +INF
  454. |
  455. |
  456. FPXT10:
  457.         jne     FPXT15          | J/ parm is a number <> 0.0
  458. |
  459.         tstw    d0
  460.         jmi     FFPINR          | J/ 0.0 to -int -> +INF
  461.         jeq     FFNANR          | J/ 0.0 to 0 -> NaN
  462.         jra     FFZERR          | J/ 0.0 to +int -> 0.0
  463. |
  464. FPXT15:
  465.         tstw    d0
  466.         jeq     FFONER          | J/ num to 0 -> 1.0
  467. |
  468.         movew   d0,a7@(6)       | Save int as its own sign flag
  469.         jpl     FPXT20
  470.         negw    d0
  471. FPXT20:
  472.         |dsw    0
  473. |
  474.         movel   sp@,sp@-        | Result init to parm value
  475. |
  476.         moveq   #16,d1          | Find MS bit of power
  477. FPXT21:
  478.         lslw    #1,d0
  479.         jcs     FPXT22          | J/ MS bit moved into carry
  480.         dbra    d1,FPXT21       | Dec d1 and jump (will always jump)
  481. |
  482. FPXT22:
  483.         movew   d0,a7@(8)       | Power pattern on stack
  484.         moveb   d1,a7@(11)      | Bit slots left count on stack
  485. |
  486. |
  487. FPXT30:
  488.         subqb   #1,a7@(11)      | Decrement bit slots left
  489.         jeq     FPXT35          | J/ evaluation complete
  490. |
  491.         movel   sp@,sp@-        | Square result value
  492.         jsr     FPMUL           | Square it
  493. |
  494.         lslw    a7@(8)          | Shift power pattern
  495.         jcc     FPXT30          | J/ this product bit not set
  496. |
  497.         movel   a7@(4),sp@-     | Copy parm value
  498.         jsr     FPMUL
  499.         jra     FPXT30
  500. |
  501. FPXT35:
  502.         tstb    a7@(10)         | Check for recipocation
  503.         jpl     FPXT36          | J/ no recipocation
  504. |
  505.         movel   #0x3F800000,sp@- |Place 1.0 on stack
  506.         jsr     FPRDIV
  507. |
  508. FPXT36:
  509.         movel   sp@+,a7@(4)     | Shift result value
  510.         addql   #4,sp           | Delete excess stack area
  511.         jmp     a4@             | Return.
  512. /*
  513. |
  514. | ###   SUBTTL  FPSQRT: Square Root Function
  515. |       page
  516. |
  517. |  Square Root Function.
  518. |
  519. |  Take square root of the single precision floating point value
  520. |  on the top of the stack.
  521. |
  522. |  Use the Newton iteration technique to compute the square root.
  523. |
  524. |      X(n+1) = (X(n) + Z/X(n)) / 2
  525. |
  526. |  The two*s exponent is scaled to restrict the solution domain to 1.0
  527. |  through 4.0.  A linear approximation to the square root is used to
  528. |  produce a first guess with greater than 4 bits of accuracy.  Three
  529. |  successive iterations are performed in registers to obtain accuracy
  530. |  of about 30 bits.  The final iteration is performed in the floating
  531. |  point domain.
  532. |
  533. | ###   PUBLIC  FPSQRT
  534. |
  535. */
  537. FPSQRT:
  538.         |dsw    0
  539.         bsr     FPFADJ
  540. |
  541.         jvs     FFNANR          | J/ NaN arg -> NaN result
  542.         jmi     FFNANR          | J/ neg arg -> NaN result
  543.         jcs     FFPINR          | J/ +INF arg -> +INF result
  544.         jeq     FFZERR          | J/ 0.0 arg -> 0.0 result
  545. |
  546.         movew   sp@,d1          | Get S/E/M word
  547.         subiw   #128*FBIAS,d1   | Extract argument's two's exp
  548.         clrb    d1              | Make it a factor of two
  549.         subw    d1,sp@          | /* Scale arg. range to 4.0 > arg' >= 1.0 */
  550.         asrw    #8,d1           | Square root of scaled two power
  551.         moveb   d1,a7@(4)       | /* Save two's exp of result on stack */
  552. |
  553.         movel   sp@,d1          | Create fixed point integer for approx
  554.         lsll    #8,d1           | /* Produce arg' * 2^30 in d1 */
  555.         bset    #31,d1          | Set implicit bit
  556.         jeq     FPSQ10          | /* J/ arg' >= 2.0 */
  557. |
  558.         lsrl    #1,d1           | Adjust d1
  559. |
  560. FPSQ10:
  561.         movew   #42720-65536,d2 | d2 = 0.325926 * 2^17
  562.         swap    d1              | /* d1.W = arg' * 2^14 */
  563.         mulu    d1,d2           | /* d2 = arg' * 0.325926 * 2^31 */
  564.         swap    d2
  565.         addiw   #23616,d2       | + 0.7207 * 2^15 - to 4+ bits
  566.         subxw   d3,d3
  567.         orw     d3,d2           | Top out approximation at 1.99997
  568. |
  569.         swap    d1
  570.         lsrl    #1,d1           | /* Arg' * 2^29 in d1 (prevent overflow) */
  571. |
  572.         movel   d1,d3           | Copy into d3
  573.         divu    d2,d3           | /* Arg'/X0 * 2^14 in d3 */
  574.         lsrw    #1,d2
  575.         addw    d3,d2           | X1 in d2 - to 8 bits
  576. |
  577.         movel   d1,d3           | Second in-register iteration
  578.         divu    d2,d3
  579.         lsrw    #1,d2
  580.         addw    d3,d2           | X2 in d2 - to 16 bits
  581. |
  582.         movel   d1,d3
  583.         divu    d2,d3
  584.         movew   d3,d4
  585.         clrw    d3
  586.         swap    d4
  587.         divu    d2,d3
  588.         movew   d3,d4           | 32 bit division result
  589.         swap    d2
  590.         clrw    d2
  591.         lsrl    #1,d2
  592.         addl    d4,d2           | X3 in d2 - to 24+ bits
  593. |
  594.         moveq   #0x7F,d1
  595.         addl    d1,d2           | Round value in d2 to 24 bits
  596.         roll    #1,d2           | Strip implicit bit
  597. |       .set    XBIAS,127               | *** for assembler bug ***
  598.         moveb   #XBIAS,d2       | Place bias value
  599.         addb    a7@(4),d2       | Scale result
  600.         rorl    #8,d2           | Position FP value
  601.         lsrl    #1,d2           | Force sign bit to 0
  602.         addql   #4,sp           | Delete four bytes from the stack
  603.         movel   d2,sp@
  604.         jmp     a4@
  605. /*
  606. |
  607. | ###   SUBTTL  FPATN: Arctangent Function
  608. |       page
  609. |
  610. |  ARCTANGENT Function.
  611. |
  612. |  The arctangent of the single precision floating point value on
  613. |  the stack computed by using a split domain with a polynomial
  614. |  approximation.  A principal range radian value is returned.
  615. |
  616. |  The domain is split at 0.25 (four) intervals.  A polynomial is
  617. |  used to approximate the arctangent for magnitudes less than 1/8.
  618. |
  619. |  Using the trigonometric identity:
  620. |         ARCTAN(y) + ARCTAN(z) = ARCTAN((y+z)/(1+yz))
  621. |         If z = (x-y)/(1+xy) then ARCTAN((y+z)/(1+yz)) = ARCTAN(x).
  622. |
  623. |         ARCTAN(-v) = -ARCTAN(v)        * make argument positive
  624. |         ARCTAN(1/v) = PI/2 - ARCTAN(v) * reduce argument to <= 1.0
  625. |
  626. |
  627. |  ARCTANGENT approximation polynomical coefficients
  628. |
  629. |         C3  =   -1.4285 71429E-01     = -1/7
  630. */
  632. FATNCN:
  633.         .long   0xBE124925
  634. |         C2  =    2.0000 00000E-01     =  1/5
  635.         .long   0x3E4CCCCD
  636. |         C1  =   -3.3333 33333E-01     = -1/3
  637.         .long   0xBEAAAAAB
  638. |         C0  =    1.0000 00000E+00     =  1/1
  639.         .long   0x3F800000
  640. |
  641. |       .set    NFATNC,4
  642. |
  643. |
  644. |
  645. |  Table of ARCTANGENT values at 0.125 intervals
  646. |
  647. |         ATAN(1/4)  =  2.4497 86631E-01
  648. FATNTB:
  649.         .long   0x3E7ADBB0
  650. |         ATAN(2/4)  =  4.6364 76090E-01
  651.         .long   0x3EED6338
  652. |         ATAN(3/4)  =  6.4350 11088E-01
  653.         .long   0x3F24BC7D
  654. |         ATAN(1/1)  =  7.8539 81634E-01   ( = PI/4)
  655.         .long   0x3F490FDB
  656. |
  657. |
  658. |  FP PI/2 (duplicate of FPIO2 because of assembler bug)
  659. |
  660. FXPIO2:
  661.         .long   0x3FC90FDB
  662. |
  663. | ###   PUBLIC  FPATN
  664. |
  665. FPATN:
  666.         |dsw    0
  667.         bsr     FPFADJ
  668. |
  669.         jvs     FFNANR          | J/ NaN arg -> NaN result
  670.         jcc     FPAT10          | J/ not INF
  671. |
  672.         movel   FXPIO2,d1       | Get PI/2
  673.         roxll   #1,d1
  674.         roxlw   sp@             | INF sign bit into X
  675.         roxrl   #1,d1           | PI/2 given sign of INF
  676.         addql   #4,sp           | Delete four bytes from the stack
  677.         movel   d1,sp@
  678.         jmp     a4@
  679. |
  680. FPAT10:
  681.         jeq     FFZERR          | J/ 0.0 arg -> 0.0 result
  682. |
  683.         bclr    #7,sp@          | Insure argument positive
  684.         sne     d0              | Create flag byte (0xFF iff negative)
  685.         andib   #0x80,d0                | Keep sign bit only
  686.         moveb   d0,a7@(6)       | Save flag byte
  687. |
  688.         movew   sp@,d1
  689.         cmpiw   #128*FBIAS+0x10,d1      | (top word of 1.125)
  690.         jlt     FPAT20                  | J/ arg < 1 + 1/8
  691. |
  692.         addqb   #1,a7@(6)
  693. |
  694.         movel   #0x3F800000,sp@- |Place 1.0 onstack
  695.         jsr     FPRDIV          | Invert the number
  696. |
  697. FPAT20:
  698.         movew   sp@,d1
  699.         cmpiw   #128*FBIAS-256,d1
  700.         jlt     FPAT30          | J/ arg < 1/4
  701. |
  702.         moveb   d1,d2           | Number of eights in d2
  703.         orib    #0x80,d2                | Implicit bit
  704.         lsrw    #7,d1
  705.         moveq   #FBIAS+4-256,d3
  706.         subb    d1,d3
  707.         lsrb    d3,d2
  708.         addqb   #1,d2
  709.         lsrb    #1,d2           | Rounded quarters in d2.B
  710. |
  711.         clrl    d0
  712.         moveb   d2,d0           | 32 bit integer in d0
  713. |
  714.         lslb    #2,d2
  715.         addb    d2,a7@(6)       | Save 4 * quarters on stack
  716. |
  717.         jsr     FLOAT
  718.         subiw   #128*2,sp@      | Produce y, floating point quarters
  719. |
  720.         movel   a7@(4),sp@-
  721.         movel   a7@(4),sp@-     | On stack:  y, arg', y, arg', <temps>
  722. |
  723.         jsr     FPMUL           | /* arg'*y */
  724.         movel   #0x3F800000,sp@-
  725.         jsr     FPADD           | /* 1.0 + arg'*y */
  726. |
  727.         movel   a7@(8),d0       | Exchange stack items
  728.         movel   sp@,a7@(8)
  729.         movel   d0,sp@          | On stack: arg', y, (1+arg'*y), <tmps>
  730.         bset    #7,a7@(4)       | Negate y
  731.         jsr     FPADD           | /* (arg'-y) */
  732.         jsr     FPRDIV          | (arg'-y) / (1 + arg'*y)
  733. |
  734. FPAT30:
  735.         |dsw    0
  736.         movel   sp@,sp@-        | Duplicate z
  737. |
  738.         pea     FATNCN
  739.         movel   sp@+,a6 | Polynomial approximation to small ATN
  740.         moveq   #NFATNC,d0
  741.         bsr     FX2SER
  742.         jsr     FPMUL           | Complete approximation
  743. |
  744.         moveb   a7@(6),d0
  745.         andiw   #0x001C,d0      | Trim to table index
  746.         jeq     FPAT40          | J/ y = 0.0, ARCTAN(y) = 0.0
  747. |
  748.         pea     FATNTB-4
  749.         movel   sp@+,a6
  750.         movel   a6@(0,d0:W),sp@-
  751.         jsr     FPADD           | Add in ARCTAN(y)
  752. |
  753. FPAT40:
  754.         btst    #0,a7@(6)       | Check for inversion
  755.         jeq     FPAT50          | J/ no inversion
  756. |
  757.         bset    #7,sp@          | Negate ARCTAN
  758.         movel   FXPIO2,sp@-
  759.         jsr     FPADD           | Inversion via subtraction
  760. |
  761. FPAT50:
  762.         tstb    a7@(6)          | Check sign of result
  763.         jpl     FPAT60          | J/ positive
  764. |
  765.         bset    #7,sp@          | Negate result
  766. |
  767. FPAT60:
  768.         movel   sp@+,sp@        | Downshift result
  769.         jmp     a4@
  770. /*
  771. |
  772. | ###   SUBTTL  FPCOS, FPSIN, FPTAN: Trigonometric Functions
  773. |       page
  774. |
  775. |  TRIG ROUTINES.
  776. |
  777. |  The support routine FTRGSV converts the radian mode argument to
  778. |  an quadrant value between -0.5 and 0.5 (quadrants).  The sign of the
  779. |  is saved (the argument is forced non-negative).  Computations are
  780. |  performed for values quadrant values of -0.5 to 0.5* other values
  781. |  are transformed as per the following transformation formulae
  782. |  and table:
  783. |
  784. |     Let:     z = shifted quadrant value
  785. |     Recall:  sin(-x) = -sin(x)    tan(-x) = -tan(x)
  786. |              cos(-x) =  cos(x)
  787. |
  788. |     l========l=================l=================l=================l
  789. |     l  QUAD  l       SIN       l       COS       l       TAN       l
  790. |     l========l=================l=================l=================l
  791. |     l    0   l      sin(z)     l      cos(z)     l      tan(z)     l
  792. |     l    1   l      cos(z)     l     -sin(z)     l   -1/tan(z)     l
  793. |     l    2   l     -sin(z)     l     -cos(z)     l      tan(z)     l
  794. |     l    3   l     -cos(z)     l      sin(z)     l   -1/tan(z)     l
  795. |     l========l=================l=================l=================l
  796. |
  797. |  The sign of the argument is, by the equations above, important only
  798. |  when computing SIN and TAN.
  799. |
  800. |  Chebyshev Polynomials are used to approximate the trigonometric
  801. |  value in the range -0.5 to 0.5 quandrants.  (Note: constants are
  802. |  adjusted to reflect calculations done in quadrant mode instead
  803. |  of radian mode).
  804. |
  805. |
  806. |
  807. |         C08  =   9.1926 02748E-04
  808. */
  810. FCOSCN:
  811.         .long   0x3A70FA83
  812. |         C06  =  -2.0863 48076E-02
  813.         .long   0xBCAAE9E4
  814. |         C04  =   2.5366 95079E-01
  815.         .long   0x3E81E0F8
  816. |         C02  =  -1.2337 00550E+00
  817.         .long   0xBF9DE9E6
  818. |         C00  =   1.0000 00000E+00
  819.         .long   0x3F800000
  820. |
  821. |       .set    NFCOSC,5
  822. |
  823. |
  824. |
  825. |         C09  =   1.6044 11848E-04
  826. FSINCN:
  827.         .long   0x39283C1A
  828. |         C07 =   -4.6817 54135E-03
  829.         .long   0xBB996966
  830. |         C05  =   7.9692 62625E-02
  831.         .long   0x3DA335E3
  832. |         C03  =  -6.4596 40975E-01
  833.         .long   0xBF255DE7
  834. |         C01  =   1.5707 96327E+00
  835. FPIO2:
  836.         .long   0x3FC90FDB
  837. |
  838. |       .set    NFSINC,5
  839. |
  840. |
  841. |
  842. |         Q3  =    1.7732 32244E+01
  843. FTNQCN:
  844.         .long   0x418DDBCC
  845. |         Q2  =   -8.0485 82645E+02
  846.         .long   0xC44936EE
  847. |         Q1  =    6.4501 55566E+03
  848.         .long   0x45C9913F
  849. |         Q0  =   -5.6630 29631E+03
  850.         .long   0xC5B0F83D
  851. |
  852. |       .set    NFTNQC,4
  853. |
  854. |
  855. |         P3  =    1.0000 00000E+00
  856. FTNPCN:
  857.         .long   0x3F800000
  858. |         P2  =   -1.5195 78276E+02
  859.         .long   0xC317F534
  860. |         P1  =    2.8156 53022E+03
  861.         .long   0x452FFA73
  862. |         P0  =   -8.8954 66142E+03
  863.         .long   0xC60AFDDD
  864. |
  865. |       .set    NFTNPC,4
  866. |
  867. |
  868. |         INV2PI  =  1 / 2*PI  =   1.5915 49431E-01
  869. FIN2PI:
  870.         .long   0x3E22F983
  871. |
  872. |
  873. | ###   PUBLIC  FPCOS
  874. |
  875. FPCOS:
  876.         bsr     FTRGSV          | Prepare argument for operation
  877. |
  878.         moveb   d0,d1           | Copy into d1
  879.         andib   #3,d0           | Strip sign bit
  880.         addqb   #1,d0
  881.         rorb    #2,d0           | Move sign bit to B7
  882.         moveb   d0,a7@(4)       | Save it
  883.         rorb    #1,d1
  884.         jcs     FSINOP          | Compute sine
  885. |
  886. |
  887. FCOSOP:
  888.         pea     FCOSCN
  889.         movel   sp@+,a6
  890.         moveq   #NFCOSC,d0
  891.         bsr     FX2SER
  892. |
  893. FTRGFN:
  894.         moveb   a7@(4),d0
  895.         andib   #0x80,d0                | Isolate sign bit
  896.         eorb    d0,sp@
  897. |
  898.         movel   sp@+,sp@
  899.         jmp     a4@
  900. |
  901. |
  902. |
  903. | ###   PUBLIC  FPSIN
  904. |
  905. FPSIN:
  906.         bsr     FTRGSV          | Prepare argument
  907. |
  908.         addib   #0x7E,a7@(4)    | Compute sign
  909.         rorb    #1,d0
  910.         jcs     FCOSOP
  911. |
  912. |
  913. FSINOP:
  914.         movel   sp@,sp@-        | Duplicate argument
  915. |
  916.         pea     FSINCN
  917.         movel   sp@+,a6
  918.         moveq   #NFSINC,d0
  919.         bsr     FX2SER
  920.         jsr     FPMUL
  921. |
  922.         jra     FTRGFN
  923. |
  924. |
  925. |
  926. | ###   PUBLIC  FPTAN
  927. |
  928. FPTAN:
  929.         bsr     FTRGSV          | Prepare argument
  930. |
  931.         rorb    #1,d0
  932.         andib   #0x80,d0
  933.         eorb    d0,a7@(4)
  934. |
  935.         movel   sp@,sp@-        | Double duplication of argument
  936.         movel   sp@,sp@-
  937. |
  938.         pea     FTNPCN
  939.         movel   sp@+,a6
  940.         moveq   #NFTNPC,d0
  941.         bsr     FX2SER
  942.         jsr     FPMUL
  943. |
  944.         movel   a7@(4),d0
  945.         movel   sp@,a7@(4)
  946.         movel   d0,sp@
  947. |
  948.         pea     FTNQCN
  949.         movel   sp@+,a6
  950.         moveq   #NFTNQC,d0
  951.         bsr     FX2SER
  952. |
  953.         btst    #0,a7@(8)
  954.         jne     FPTN20          | J/ reverse division
  955. |
  956.         jsr     FPDIV
  957.         jra     FTRGFN
  958. |
  959. FPTN20:
  960.         jsr     FPRDIV
  961.         jra     FTRGFN
  962. /*
  963. |
  964. |       page
  965. |
  966. |  Trigonometric Service Routine
  967. |
  968. |  The single precision floating point value to be processed is on the
  969. |  stack.  Excess 2 PI's are removed from the argument.  The range of
  970. |  the argument is then scaled to within PI/4 of 0.  The shift size
  971. |  (in number of quadrants) is returned in D0.  The original sign bit
  972. |  is returned in the D0.B sign bit.
  973. |
  974. |  If the argument is NaN, +INF, -INF, or too large (>= 2^16), this
  975. |  routine causes a return to the caller with a NaN result.
  976. |
  977. */
  979. FTRGSV:
  980.         moveal  sp@+,a6 | FTRGSV return addr
  981. |
  982.         bsr     FPFADJ
  983.         jcs     FFNANR          | Return NAN if +/- INF
  984.         jvs     FFNANR          | Return NaN if NaN
  985. |
  986.         clrw    a7@(4)          | Create arg type return byte
  987. |
  988.         roxlw   sp@             | Sign bit into X
  989.         roxrw   a7@(4)          | Sign bit into temp byte
  990.         roxrw   sp@             | Restore FP value (with sign = +)
  991. |
  992.         cmpiw   #128*FBIAS+2048,sp@     | (2048 = 128*16)
  993.         jge     FFNANR                  | J/ return NaN if ABS(arg) >= 2^16
  994. |
  995.         movel   FIN2PI,sp@-     | Scale arg by 1/(2*PI)
  996.         jsr     FPMUL
  997. |
  998.         cmpiw   #128*FBIAS,sp@
  999.         jlt     FTS10           | /* J/ no excess two PI's */
  1000. |
  1001.         movel   sp@,sp@-        | Double the scaled argument
  1002.         jsr     INT
  1003.         negl    d0
  1004.         jsr     FLOAT
  1005.         jsr     FPADD
  1006. |
  1007. FTS10:
  1008.         tstw    sp@
  1009.         jeq     FTS20
  1010. |
  1011.         addiw   #128*3,sp@      | Multiply value by 8
  1012.         cmpiw   #128*FBIAS,sp@
  1013.         jlt     FTS20           | J/ < 1
  1014. |
  1015.         movel   sp@,sp@-        | Double parameter
  1016.         jsr     INT
  1017.         addqb   #1,d0
  1018.         lsrb    #1,d0
  1019.         addb    d0,a7@(4)       | Mix quadrant number with sign bit
  1020. |
  1021.         lslb    #1,d0
  1022.         negl    d0
  1023.         jsr     FLOAT
  1024.         jsr     FPADD
  1025. |
  1026. FTS20:
  1027.         tstw    sp@
  1028.         jeq     FTS30           | J/ arg = 0.0
  1029. |
  1030.         subiw   #128*1,sp@      | Range result to -0.5 to 0.5
  1031. |
  1032. FTS30:
  1033.         andib   #0x83,a7@(4)    | Limit quadrant shift to 0 thru 3
  1034.         moveb   a7@(4),d0       | Return cntl byte in d0 and at a7@(4)
  1035.         jmp     a6@
  1036. /*
  1037. |
  1038. | ###   SUBTTL  Polynomial Evaluation Routine
  1039. |       page
  1040. |
  1041. |  Polynomial Evaluation Routines.
  1042. |
  1043. |  A list of constants is pointed to by A6.  X is on the stack.  D0
  1044. |  contains the number of constants (the polynomial degree plus one).
  1045. |
  1046. |  A6 enters pointing to C[1].  Upon return, the value on stack has
  1047. |  been replaced with:
  1048. |  FXSER   = C[1]*X^(N-1)  + C[2]*X^(N-2)  + ... + C[N-1]*X   + C[N]
  1049. |
  1050. |  FX2SER  = C[1]*X^(2N-2) + C[2]*X^(2N-4) + ... + C[N-1]*X^2 + C[N]
  1051. |
  1052. |
  1053. */
  1055. FX2SER:
  1056.         moveal  a4,a5           | Save a4 in a5
  1057.         bsr     FPFADJ
  1058. |
  1059.         moveb   d0,a7@(4)       | Save count
  1060. |
  1061.         movel   sp@,sp@-        | Square the argument
  1062.         jsr     FPMUL
  1063.         jra     FXSR01          | J/ join FXSER routine
  1064. |
  1065. FXSER:
  1066.         moveal  a4,a5
  1067.         bsr     FPFADJ
  1068. |
  1069.         moveb   d0,a7@(4)
  1070. |
  1071. FXSR01:
  1072.         movel   a6@+,sp@-       | Place C[1] on stack as accum
  1073. |
  1074.         subqb   #1,a7@(8)
  1075. |
  1076. FXSR02:
  1077.         movel   a7@(4),sp@-     | Copy X
  1078.         jsr     FPMUL
  1079.         movel   a6@+,sp@-       | Add in next coefficient
  1080.         jsr     FPADD
  1081. |
  1082.         subqb   #1,a7@(8)
  1083.         jne     FXSR02          | J/ more coefficients
  1084. |
  1085.         movel   sp@,a7@(8)      | Adjust return
  1086.         addql   #8,sp
  1087. |
  1088.         moveal  a4,a3
  1089.         moveal  a5,a4           | Restore return address
  1090.         jmp     a3@
  1091. /*
  1092. |
  1093. | ###   SUBTTL  Function Completion Segments
  1094. |       page
  1095. |
  1096. |  FPFADJ - Adjust the stack for single precision functions.
  1097. |          Move the return address to the function caller into A4.
  1098. |          Shift the single precision floating point argument down
  1099. |          four bytes, leaving a temporary area of four bytes at
  1100. |          4(A7) thru 7(A7).  Set the condition code bits N, Z,
  1101. |          C (INF), and V (NaN) to quickly type the argument.
  1102. |
  1103. */
  1105. FPFADJ:
  1106.         moveal  sp@+,a3 | Function return addr
  1107.         moveal  sp@,a4          | Caller return addr
  1108.         movel   a7@(4),sp@      | Down shift parameter
  1109. |
  1110.         movew   sp@,d2          | Get SEM word
  1111.         rolw    #1,d2
  1112.         cmpiw   #256*0xFF,d2
  1113.         jcs     FPFA10          | J/ not INF or NaN
  1114. |
  1115.         andiw   #0xFE,d2
  1116.         jeq     FPFA05
  1117. |
  1118.         orib    #0x02,ccr       | Set -V- bit -> NaN argument
  1119. | ##    DC.W    $003C,$0002
  1120.         jmp     a3@
  1121. |
  1122. FPFA05:
  1123.         tstb    sp@             | /* Set N bit as req'd, clear Z bit. */
  1124.         orib    #0x01,ccr       | Set -C- bit -> INF argument
  1125. | ##    DC.W    $003C,$0001
  1126.         jmp     a3@
  1127. |
  1128. FPFA10:
  1129.         tstw    sp@             | /* Set N, Z bits as req'd, clear V, C */
  1130.         jmp     a3@
  1131. |
  1132. |
  1133. |
  1134. FFPINR:
  1135.         |dsw    0               | Result is positive overflow
  1136.         moveq   #0,d0           | Sign is positive
  1137.         jra     FFIN01
  1138. |
  1139. FFMINR:
  1140.         |dsw    0               | Result in negative overflow
  1141.         moveq   #-1,d0          | Sign is negative
  1142. |
  1143. FFIN01:
  1144.         moveaw  d0,a2           | Sign into a2
  1145.         addql   #8,sp           | Delete eight bytes from the stack
  1146.         moveal  a4,a0           | Return addr in a0
  1147.         jmp     FINFRS          | (Use FPOPNS exception return)
  1148. |
  1149. |
  1150. FFNANR:
  1151.         |dsw    0               | Result is NaN
  1152.         addql   #8,sp           | Delete eight bytes from the stack
  1153.         moveal  a4,a0
  1154.         jmp     FNANRS
  1155. |
  1156. |
  1157. FFZERR:
  1158.         |dsw    0               | Result is zero
  1159.         addql   #8,sp
  1160.         moveal  a4,a0
  1161.         jmp     FZERRS
  1162. |
  1163. |
  1164. FFUNFR:
  1165.         |dsw    0               | Underflow result
  1166.         addql   #8,sp           | Delete eight bytes from the stack
  1167.         moveal  a4,a0           | Return address into a0
  1168.         jmp     FUNFRS
  1169. |
  1170. |
  1171. FFONER:
  1172.         |dsw    0               | Result is 1.0
  1173.         addql   #4,sp
  1174.         movel   #0x3F800000,sp@
  1175.         clrb    FPERR
  1176.         jmp     a4@
  1177. |
  1178. |
  1179. |       end