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

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

VxWorks

开发平台：

C/C++

uss_dpfncs.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/DPFNCS: IEEE Double Precision Functions
  25. |DPFNCS idnt    1,0             ; IEEE Double Precision Functions
  26. |                               ; DPFNCS.A68
  27. |
  28. |* * * * * * * * * *
  29. |
  30. |       68000 DPAC  --  Double 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. |  Double Precision Floating Point Function library.
  55. |
  56. |  These functions operate on double precision floating point values,
  57. |  yielding double 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        -708.40 < a < 709.78    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 forward branches to SHORT
  94. |
  95.         .globl  DPEXP
  96.         .globl  DPLN
  97.         .globl  DPLOG
  98.         .globl  DPXTOI
  99.         .globl  DPSQRT
  100.         .globl  DPATN
  101.         .globl  DPCOS
  102.         .globl  DPSIN
  103.         .globl  DPTAN
  104. |
  105. |
  106. |       xref    DPMUL           ;Floating point operations
  107. |       xref    DPDIV
  108. |       xref    DPRDIV
  109. |       xref    DPADD
  110. |       xref    DAINT           ;Basic floating point functions
  111. |       xref    DFLOAT
  112. |       xref    DINT
  113. |
  114. |       xref    DNANRS
  115. |       xref    DINFRS
  116. |       xref    DUNFRS
  117. |       xref    DZERRS
  118. |
  119. |
  120. |
  121. |
  122. |       .set    DBIAS,1023
  123. |
  124. |
  125. |       xref    FPERR
  126. |
  127. |
  128. |       .set    ERUNF,1         |FPERR UNDERFLOW ERROR CODE
  129. |       .set    EROVF,2         |FPERR OVERFLOW ERROR CODE
  130. |       .set    ERNAN,3         |FPERR INVALID OPERATION ERROR CODE
  131. /*
  132. |
  133. |
  134. |
  135. | ###   SUBTTL  DPEXP: Exponentiation of e Function
  136. |       page
  137. |
  138. |
  139. |  EXP:   Calculate e to the power specified by the stack argument.
  140. |         The result is returned on the stack.
  141. |         placed at DI.
  142. |
  143. |         FPERR is set upon overflow, underflow, or NaN argument.
  144. |
  145. |  Algorithm:
  146. |  The argument is multiplied by 1/LN 2.  The integer portion of the
  147. |  product becomes the two*s exponent.  A Taylor series approximation
  148. |  about zero is used to compute 2**FRAC(ARG * 1/LN 2).
  149. |
  150. */
  151.         .text
  152. |
  153. DEXPCN:
  154.         |dsw    0               | EXP CONSTANT LIST (2**X EVALUATION)
  155. |
  156. | ???     (LN 2)^16/16!   =   1.3570 24794 87551 47193 D-16
  157. | ???   DC.L    $3CA38E89,$AE79F8B4
  158. |         (LN 2)^15/15!   =   3.1324 36707 08842 86216 D-15
  159.         .long   0x3CEC36E8,0x43B04022
  160. |         (LN 2)^14/14!   =   6.7787 26354 82254 56335 D-14
  161.         .long   0x3D331496,0x4D5878A9
  162. |         (LN 2)^13/13!   =   1.3691 48885 39041 28881 D-12
  163.         .long   0x3D781619,0x3166D0F9
  164. |         (LN 2)^12/12!   =   2.5678 43599 34882 05142 D-11
  165.         .long   0x3DBC3BD6,0x50FC2986
  166. |         (LN 2)^11/11!   =   4.4455 38271 87081 14976 D-10
  167.         .long   0x3DFE8CAC,0x7351BB25
  168. |         (LN 2)^10/10!   =   7.0549 11620 80112 33299 D-09
  169.         .long   0x3E3E4CF5,0x158B8ECA
  170. |         (LN 2)^ 9/ 9!   =   1.0178 08600 92396 99728 D-07
  171.         .long   0x3E7B5253,0xD395E7C4
  172. |         (LN 2)^ 8/ 8!   =   1.3215 48679 01443 09488 D-06
  173.         .long   0x3EB62C02,0x23A5C824
  174. |         (LN 2)^ 7/ 7!   =   1.5252 73380 40598 40280 D-05
  175.         .long   0x3EEFFCBF,0xC588B0C7
  176. |         (LN 2)^ 6/ 6!   =   1.5403 53039 33816 09954 D-04
  177.         .long   0x3F243091,0x2F86C787
  178. |         (LN 2)^ 5/ 5!   =   1.3333 55814 64284 43423 D-03
  179.         .long   0x3F55D87F,0xE78A6731
  180. |         (LN 2)^ 4/ 4!   =   9.6181 29107 62847 71620 D-03
  181.         .long   0x3F83B2AB,0x6FBA4E77
  182. |         (LN 2)^ 3/ 3!   =   5.5504 10866 48215 79953 D-02
  183.         .long   0x3FAC6B08,0xD704A0C0
  184. |         (LN 2)^ 2/ 2!   =   2.4022 65069 59100 71233 D-01
  185.         .long   0x3FCEBFBD,0xFF82C58F
  186. |         (LN 2)^1/ 1!   =   6.9314 71805 59945 30942 D-01
  187. DLN2:
  188.         .long   0x3FE62E42,0xFEFA39EF
  189. |                             1.0000 00000 00000 00000 D+00
  190.         .long   0x3FF00000,0x00000000
  191. |
  192. |       .set    DNEXPC,16               | Number of constants in DEXP poly
  193. |
  194. |         1 / LN 2  =  1.4426 95040 88896 34074 D+00
  195. DINLN2:
  196.         .long   0x3FF71547,0x652B82FE
  197. |
  198. |
  199. |
  200. | ###   PUBLIC  DPEXP
  201. |
  202. DPEXP:
  203.         |dsw    0
  204.         bsr     DPFADJ          | Set-up for DP functions
  205. |
  206.         jvs     DFNANR          | NaN parm -> NaN result
  207.         jcc     DEXP01          | J/ arg not INF
  208. |
  209.         jmi     DFUNFR          | -INF parm -> underflow result
  210.         jra     DFPINR          | +INF parm -> positive INF result
  211. |
  212. DEXP01:
  213.         movel   DINLN2+4,sp@-   | 1/LN 2 on stack
  214.         movel   DINLN2+0,sp@-
  215.         jsr     DPMUL           | Convert to 2^ function
  216. |
  217.         movew   sp@,d1
  218.         andiw   #0x7FF0,d1
  219.         cmpiw   #16*DBIAS+144,d1
  220.         jhi     DEXP10          | J/ overflow or underflow  >= +/- 2^10
  221.         jne     DEXP20          | J/ within range            < +/- 2^09
  222.         tstb    sp@
  223.         jpl     DEXP20          | J/ result >= 1.0
  224. |
  225.         movel   sp@,d2          | Examine top 12 bits of mantissa
  226.         rorl    #8,d2
  227.         andiw   #0x0FFF,d2
  228.         cmpiw   #0x0FF8,d2
  229.         jcs     DEXP20          | J/ result >= 2^-1022
  230. |
  231. DEXP10:
  232.         tstb    sp@
  233.         jpl     DFPINR          | J/ DP function +INF result
  234.         jra     DFUNFR          | J/ DP function underflow result
  235. |
  236. DEXP20:
  237.         movel   a7@(4),sp@-
  238.         movel   a7@(4),sp@-     | Double the argument on the stack
  239. |
  240.         jsr     DINT
  241. |
  242.         movew   d1,a7@(8)       | Save result of DINT
  243. |
  244.         negl    d1
  245.         negxl   d0              | Negate integer
  246.         jsr     DFLOAT
  247.         jsr     DPADD           | Special DFRAC
  248. |
  249.         pea     DEXPCN
  250.         movel   sp@+,a6 | Point a6 to constant list
  251.         moveq   #DNEXPC ,d0      | d0 holds the number of constants
  252.         moveq   #DNEXPC, d0      | d0 holds the number of constants
  253.         moveq   #16,d0      | XXX DEBUG DEBUG
  254.         bsr     DXSER           | Polynomial series
  255. |
  256.         movew   a7@(8),d0       | Fetch power of two scaling
  257.         lslw    #4,d0
  258.         swap    d0
  259.         clrw    d0
  260.         addl    sp@+,d0 | Combine with DXSER result
  261. |
  262.         movel   sp@,a7@(4)      | Downshift result
  263.         movel   d0,sp@
  264.         jmp     a4@             | Return
  266. /*
  267. |
  268. | ###   SUBTTL  DPLOG, DPLN: Logarithm Functions
  269. |       page
  270. |
  271. |  LOGARITHM (BASE E) FUNCTION.
  272. |
  273. |  The natural logarithm of the double precision floating point
  274. |  value on the stack is computed using a split domain polynomial
  275. |  approximation.  The common log is computed by scaling the
  276. |  result of the natural log computation.
  277. |
  278. |  TWOS = FLOAT(ARG.EXPONENT) * LOG(2)
  279. |
  280. |  The argument, (after the two's power scaling, 2.0 > arg' >= 1.0),
  281. |  is sorted into one of five classes.  Within that class, a center
  282. |  point with a known log value is used in combination with a 9th
  283. |  degree Chebyshey polynomial approximation, to calculate a
  284. |  logarithm.
  285. |
  286. |  The center point of each class is adjusted so that its natural
  287. |  logarithm expressed in double precision form is accurate to
  288. |  more than 80 bits.
  289. |
  290. |  log(arg) = two's log + center point log + polynomial approximation
  291. |
  292. |  Polynomial approximation:
  293. |
  294. |  log((1-t)/(1+t)) = t*(c1*t^8 + c2*t^6 + c3*t^4 + c4*t^2 + c5)
  295. |         (accurate with 1.0E-17 for t within 0.1 of 1.0)
  296. |
  297. |  t = - (center - arg')/(center + arg')  [to calc log of arg'/center]
  298. |
  299. */
  301. DLNCNS:
  302.         |dsw    0               | LN constants list
  303. |
  304. |         c1  =   0.22330 022
  305.         .long   0x3FCC951A,0x03000000
  306. |         c2  =   0.28571 20487
  307.         .long   0x3FD2491B,0x30450000
  308. |         c3  =   0.40000 00018 947
  309.         .long   0x3FD99999,0x9BA26900
  310. |         c4  =   0.66666 66666 66105
  311.         .long   0x3FE55555,0x55554192
  312. |         c5  =   2.00000 00000 00000 00000
  313.         .long   0x40000000,0x00000000
  314. |
  315. |       .set    DNLNCN,5                | Number of constants in DLN poly
  316. |
  317. |
  318. DLNCEN:
  319.         |dsw    0               | Center Points
  320. |         Center  #1 = 1.0000 00000 00000 00000
  321.         .long   0x3FF00000,0x00000000
  322. |         Center  #2 = 1.1892 07115 00272 05647  (close to 2^(1/4))
  323.         .long   0x3FF306FE,0x0A31B713
  324. |         Center  #3 = 1.4142 13562 37309 38548  (close to 2^(1/2))
  325.         .long   0x3FF6A09E,0x667F3BC7
  326. |         Center  #4 = 1.6817 92830 50742 54625  (close to 2^(3/4))
  327.         .long   0x3FFAE89F,0x995AD39D
  328. |         Center  #5 = 2.0000 00000 00000 00000
  329. | ---     shifted to 1.0 center..
  330. |
  331. DLNLOG:
  332.         |dsw    0               | DLNCEN values adj for exact ln values
  333. |         LN(Center #2) = 0.17328 67951 39985 90522
  334.         .long   0x3FC62E42,0xFEFA39E0
  335. |         LN(Center #3) = 0.34657 35902 79971 81045
  336.         .long   0x3FD62E42,0xFEFA39E0
  337. |         LN(Center #4) = 0.51986 03854 19956 82749
  338.         .long   0x3FE0A2B2,0x3F3BAB60
  339. |
  340. |
  341. |         1 / LN 10  =  0.43429 44819 03251 82765
  342. DILN10:
  343.         .long   0x3FDBCB7B,0x1526E50E
  344. |
  345. |
  346. |
  347. | ###   PUBLIC  DPLN
  348. | ###   PUBLIC  DPLOG
  349. |
  350. DPLOG:
  351.         |dsw    0
  352.         moveq   #-1,d0          | Signal common log scaling
  353.         jra     DLN000
  354. |
  355. DPLN:
  356.         |dsw    0
  357.         clrb    d0              | Signal no post calculation scaling
  358. |
  359. DLN000:
  360.         bsr     DPFADJ          | Adjust parameter on stack
  361. |
  362.         jvs     DFNANR          | J/ NaN arg -> NaN result
  363.         jmi     DFNANR          | J/ Neg arg -> NaN result
  364.         jeq     DFMINR          | J/ 0.0 arg -> -INF result
  365.         jcs     DFPINR          | J/ +INF arg -> +INF result
  366. |
  367.         moveb   d0,a7@(10)      | Save natural/common flag
  368. |
  369.         movew   sp@,d1          | Prepare to calc parm sign and exp
  370.         lsrw    #4,d1           | Position exponent (right justified)
  371.         subiw   #DBIAS,d1
  372.         movew   d1,a7@(8)       | /* Save two's exponent value */
  373. |
  374.         lslw    #4,d1           | Scale parameter
  375.         subw    d1,sp@
  376. |
  377.         clrw    d0              | Set class number to 0
  378.         movel   sp@,d1
  379.         lsrl    #8,d1           | d1.W = 0xEMMM
  380. |
  381.         cmpiw   #0xF000+371,d1
  382.         jcs     DLN050          | J/ < 1 +  371/4096  (1.090576)
  383.         addqw   #8,d0
  384.         cmpiw   #0xF000+1216,d1
  385.         jcs     DLN050          | J/ < 1 + 1216/4096  (1.296875)
  386.         addqw   #8,d0
  387.         cmpiw   #0xF000+2221,d1
  388.         jcs     DLN050          | J/ < 1 + 2221/4096  (1.542236)
  389.         addqw   #8,d0
  390.         cmpiw   #0xF000+3416,d1
  391.         jcs     DLN050          | J/ < 1 + 3416/4096  (1.833984)
  392.         addqw   #1,a7@(8)       | /* Class 4: Bump two's exp, then class 0 */
  393.         clrw    d0
  394.         subib   #0x10,a7@(1)    | Reduce scaled exponent
  395. |
  396. DLN050:
  397.         moveb   d0,a7@(11)      | Save center number
  398. |
  399.         pea     DLNCEN
  400.         movel   sp@+,a2
  401.         movel   a2@(4,d0:W),sp@- |Place center on stack
  402.         movel   a2@(0,d0:W),sp@-
  403. |
  404.         movel   a7@(12),sp@-    | Copy scaled argument
  405.         movel   a7@(12),sp@-
  406. |
  407.         movel   a2@(4,d0:W),sp@- |Place center on stack
  408.         movel   a2@(0,d0:W),sp@-
  409. |
  410.         jsr     DPADD           | /* Calc  (CENTER + PARM') */
  411. |
  412.         movel   sp@,d0          | Exch top stack item w/ 3rd
  413.         movel   a7@(16),d1
  414.         movel   d1,sp@
  415.         movel   d0,a7@(16)
  416.         movel   a7@(4),d0
  417.         movel   a7@(20),d1
  418.         movel   d1,a7@(4)
  419.         movel   d0,a7@(20)
  420. |
  421.         bset    #7,a7@(8)       | Negate CENTER
  422. |
  423.         jsr     DPADD           | /* Calc  - (CENTER - PARM') */
  424.         jsr     DPRDIV          | /* T = - (CENTER-PARM')/(CENTER+PARM) */
  425. |
  426.         movel   a7@(4),sp@-     | Duplicate t
  427.         movel   a7@(4),sp@-
  428. |
  429.         pea     DLNCNS
  430.         movel   sp@+,a6 | Poly approx
  431.         moveq   #DNLNCN,d0
  432.         bsr     DX2SER
  433. |
  434.         jsr     DPMUL           | Times t
  435.         clrw    d0
  436.         moveb   a7@(11),d0
  437.         jeq     DLN060          | J/ center = 0 -> no additive log val
  438. |
  439.         pea     DLNLOG-8
  440.         movel   sp@+,a6
  441.         movel   a6@(4,d0:W),sp@-
  442.         movel   a6@(0,d0:W),sp@-
  443.         jsr     DPADD           | Add in center log value
  444. |
  445. DLN060:
  446.         movew   a7@(8),d1       | /* Get two's exponent value */
  447.         clrl    d0
  448.         extl    d1
  449.         jpl     DLN061
  450.         moveq   #-1,d0
  451. DLN061:
  452.         |dsw    0
  453. |
  454.         jsr     DFLOAT
  455. |
  456.         movel   DLN2+4,sp@-     | Log of 2 on stack
  457.         movel   DLN2+0,sp@-
  458.         jsr     DPMUL
  459.         jsr     DPADD
  460. |
  461.         tstb    a7@(10)
  462.         jeq     DLN070          | J/ natural log
  463. |
  464.         movel   DILN10+4,sp@-   | Scaling value
  465.         movel   DILN10+0,sp@-
  466.         jsr     DPMUL
  467. |
  468. DLN070:
  469.         movel   a7@(4),a7@(8)
  470.         movel   sp@+,sp@
  471.         jmp     a4@
  473. /*
  474. | ###   SUBTTL  DPXTOI: Floating Point Number to Integer Power Function
  475. |       page
  476. |
  477. |  X to I power Function.
  478. |
  479. |  The double precision floating point value on the stack is
  480. |  raised to the power specified in D0.W then returned on stack.
  481. |
  482. |  A shift and multiply technique is used (possibly with a trailing
  483. |  recipication.
  484. |
  485. |
  486. | ###   PUBLIC  DPXTOI
  487. |
  488. */
  490. DPXTOI:
  491.         |dsw    0
  492.         bsr     DPFADJ
  493. |
  494.         jvs     DFNANR          | J/ NaN arg -> NaN result
  495.         jcc     DPXT10          | J/ arg is not INF
  496. |
  497.         tstw    d0
  498.         jeq     DFNANR          | J/ arg is +/- INF, I is 0 -> NaN
  499.         jmi     DFUNFR          | J/ x is +/- INF, I < 0 -> underflow
  500. |
  501.         rorw    #1,d0
  502.         andw    sp@,d0
  503.         jmi     DFMINR          | J/ arg is -INF, I is +odd -> -INF
  504.         jra     DFPINR          | else -> +INF
  505. |
  506. |
  507. DPXT10:
  508.         jne     DPXT15          | J/ parm is a number <> 0.0
  509. |
  510.         tstw    d0
  511.         jmi     DFPINR          | J/ 0.0 to -int -> +INF
  512.         jeq     DFNANR          | J/ 0.0 to 0 -> NaN
  513.         jra     DFZERR          | J/ 0.0 to +int -> 0.0
  514. |
  515. DPXT15:
  516.         tstw    d0
  517.         jeq     DFONER          | J/ num to 0 -> 1.0
  518. |
  519.         movew   d0,a7@(10)      | Save int as its own sign flag
  520.         jpl     DPXT20
  521.         negw    d0
  522. DPXT20:
  523.         |dsw    0
  524. |
  525.         movel   a7@(4),sp@-     | Result init to parm value
  526.         movel   a7@(4),sp@-
  527. |
  528.         moveq   #16,d1          | Find MS bit of power
  529. DPXT21:
  530.         lslw    #1,d0
  531.         jcs     DPXT22          | J/ MS bit moved into carry
  532.         dbra    d1,DPXT21       | Dec d1 and jump (will always jump)
  533. |
  534. DPXT22:
  535.         movew   d0,a7@(16)      | Power pattern on stack
  536.         moveb   d1,a7@(19)      | Bit slots left count on stack
  537. |
  538. |
  539. DPXT30:
  540.         subqb   #1,a7@(19)      | Decrement bit slots left
  541.         jeq     DPXT35          | J/ evaluation complete
  542. |
  543.         movel   a7@(4),sp@-     | Square result value
  544.         movel   a7@(4),sp@-
  545.         jsr     DPMUL           | Square it
  546. |
  547.         lslw    a7@(16)         | Shift power pattern
  548.         jcc     DPXT30          | J/ this product bit not set
  549. |
  550.         movel   a7@(12),sp@-    | Copy parm value
  551.         movel   a7@(12),sp@-
  552.         jsr     DPMUL
  553.         jra     DPXT30
  554. |
  555. DPXT35:
  556.         tstb    a7@(18)         | Check for recipocation
  557.         jpl     DPXT36          | J/ no recipocation
  558. |
  559.         clrl    sp@-            | Place 1.0 on stack
  560.         movel   #0x3FF00000,sp@-
  561.         jsr     DPRDIV
  562. |
  563. DPXT36:
  564.         movel   sp@+,a7@(8)     | Shift result value
  565.         movel   sp@+,a7@(8)
  566.         addql   #4,sp           | Delete excess stack area
  567.         jmp     a4@             | Return.
  568. /*
  569. |
  570. | ###   SUBTTL  DPSQRT: Square Root Function
  571. |       page
  572. |
  573. |  Square Root Function.
  574. |
  575. |  Take square root of the double precision floating point value
  576. |  on the top of the stack.
  577. |
  578. |  Use the Newton iteration technique to compute the square root.
  579. |
  580. |      X(n+1) = (X(n) + Z/X(n)) / 2
  581. |
  582. |  The two*s exponent is scaled to restrict the solution domain to 1.0
  583. |  through 4.0.  A linear approximation to the square root is used to
  584. |  produce a first guess with greater than 4 bits of accuracy.  Three
  585. |  successive iterations are performed in registers to obtain accuracy
  586. |  of about 30 bits.  The final iteration is performed in the floating
  587. |  point domain.
  588. |
  589. | ###   PUBLIC  DPSQRT
  590. |
  591. */
  593. DPSQRT:
  594.         |dsw    0
  595.         bsr     DPFADJ
  596. |
  597.         jvs     DFNANR          | J/ NaN arg -> NaN result
  598.         jmi     DFNANR          | J/ neg arg -> NaN result
  599.         jcs     DFPINR          | J/ +INF arg -> +INF result
  600.         jeq     DFZERR          | J/ 0.0 arg -> 0.0 result
  601. |
  602.         movew   sp@,d1          | Get S/E/M word
  603.         subiw   #16*DBIAS,d1    | Extract argument's two's exp
  604.         andib   #0xE0,d1                | Make it a factor of two
  605.         subw    d1,sp@          | /* Scale arg. range to 4.0 > arg' >= 1.0 */
  606.         asrw    #1,d1           | Square root of scaled two power
  607.         movew   d1,a7@(8)       | /* Save two's exp of result on stack */
  608. |
  609.         movel   sp@,d1          | Create fixed point integer for approx
  610.         movew   a7@(4),d2
  611.         lsll    #8,d1           | /* Produce arg' * 2^30 in d1 */
  612.         lsll    #3,d1
  613.         lsrl    #5,d2
  614.         orw     d2,d1
  615.         bset    #31,d1          | Set implicit bit
  616.         jeq     DPSQ10          | /* J/ arg' >= 2.0 */
  617. |
  618.         lsrl    #1,d1           | Adjust d1
  619. |
  620. DPSQ10:
  621.         movew   #42720-65536,d2 | d2 = 0.325926 * 2^17
  622.         swap    d1              | /* d1.W = arg' * 2^14 */
  623.         mulu    d1,d2           | /* d2 = arg' * 0.325926 * 2^31 */
  624.         swap    d2
  625.         addiw   #23616,d2       | + 0.7207 * 2^15 - to 4+ bits
  626.         subxw   d3,d3
  627.         orw     d3,d2           | Top out approximation at 1.99997
  628. |
  629.         swap    d1
  630.         lsrl    #1,d1           | /* Arg' * 2^29 in d1 (prevent overflow) */
  631. |
  632.         movel   d1,d3           | Copy into d3
  633.         divu    d2,d3           | /* Arg'/X0 * 2^14 in d3 */
  634.         lsrw    #1,d2
  635.         addw    d3,d2           | X1 in d2 - to 8 bits
  636. |
  637.         movel   d1,d3           | Second in-register iteration
  638.         divu    d2,d3
  639.         lsrw    #1,d2
  640.         addw    d3,d2           | X2 in d2 - to 16 bits
  641. |
  642.         movel   d1,d3
  643.         divu    d2,d3
  644.         movew   d3,d4
  645.         clrw    d3
  646.         swap    d4
  647.         divu    d2,d3
  648.         movew   d3,d4           | 32 bit division result
  649.         swap    d2
  650.         clrw    d2
  651.         lsrl    #1,d2
  652.         addl    d4,d2           | X3 in d2 - to 29 bits
  653.         subxl   d4,d4
  654.         orl     d4,d2           | Top out at 1.9999999995
  655. |
  656.         movel   a7@(4),sp@-     | /* Down shift arg' */
  657.         movel   a7@(4),sp@-
  658. |
  659.         lsll    #1,d2           | Create DP of X3 (good to 29 bits) ...
  660.         clrl    d3              | ... in d2:d3
  661.         movew   d2,d3
  662.         andiw   #0xFFF,d3
  663.         eorw    d3,d2
  664.         oriw    #DBIAS,d2       | Scale to floating point
  665.         moveq   #12,d0          | Shift count in d0
  666.         rorl    d0,d2           | Position bits
  667.         rorl    d0,d3
  668. |
  669.         movel   d2,a7@(8)       | /* On stack: X4, arg', X4, <flags> */
  670.         movel   d3,a7@(12)
  671.         movel   d3,sp@-
  672.         movel   d2,sp@-
  673.         jsr     DPDIV           | Last iter - to X4 - in DP domain
  674.         jsr     DPADD
  675.         subiw   #0x0010,sp@     | "Divide" by 2
  676. |
  677.         movew   a7@(8),d7
  678.         addw    d7,sp@          | Scale result
  679.         movel   a7@(4),a7@(8)   | Down shift result
  680.         movel   sp@+,sp@
  681.         jmp     a4@
  683. /*
  684. | ###   SUBTTL  DPATN: Arctangent Function
  685. |       page
  686. |
  687. |  ARCTANGENT Function.
  688. |
  689. |  The arctangent of the double precision floating point value at SI is
  690. |  computed by using a split domain with a polynomial approximation.
  691. |  A principal range radian value is returned.
  692. |
  693. |  The domain is split at 0.125 (eight) intervals.  A polynomial is
  694. |  used to approximate the arctangent for magnitudes less than 1/16.
  695. |
  696. |  Using the trigonometric identity:
  697. |         ARCTAN(y) + ARCTAN(z) = ARCTAN((y+z)/(1+yz))
  698. |         If z = (x-y)/(1+xy) then ARCTAN((y+z)/(1+yz)) = ARCTAN(x).
  699. |
  700. |         ARCTAN(-v) = -ARCTAN(v)        * make argument positive
  701. |         ARCTAN(1/v) = PI/2 - ARCTAN(v) * reduce argument to <= 1.0
  702. |
  703. |
  704. |  ARCTANGENT approximation polynomical coefficients
  705. |
  706. |         C4  =    1.1022 81616 12614 90000E-01
  707. */
  709. DATNCN:
  710.         .long   0x3FBC37E9,0xAD397134
  711. |         C3  =   -1.4285 41305 08745 00000E-01
  712.         .long   0xBFC2490B,0x4D511901
  713. |         C2  =    1.9999 99958 01446 40000E-01
  714.         .long   0x3FC99999,0x90956BB6
  715. |         C1  =   -3.3333 33333 31284 50000E-01
  716.         .long   0xBFD55555,0x5554C529
  717. |         C0  =    1.0000 00000 00000 00000E+00
  718.         .long   0x3FF00000,0x00000000
  719. |
  720. |       .set    NDATNC,5
  721. |
  722. |
  723. |
  724. |  Table of ARCTANGENT values at 0.125 intervals
  725. |
  726. |         ATAN(1/8)  =  1.2435 49945 46761 43503E-01
  727. DATNTB:
  728.         .long   0x3FBFD5BA,0x9AAC2F6E
  729. |         ATAN(2/8)  =  2.4497 86631 26864 15417E-01
  730.         .long   0x3FCF5B75,0xF92C80DE
  731. |         ATAN(3/8)  =  3.5877 06702 70572 22040E-01
  732.         .long   0x3FD6F619,0x41E4DEF1
  733. |         ATAN(4/8)  =  4.6364 76090 00806 11621E-01
  734.         .long   0x3FDDAC67,0x0561BB4F
  735. |         ATAN(5/8)  =  5.5859 93153 43562 43597E-01
  736.         .long   0x3FE1E00B,0xABDEFEB4
  737. |         ATAN(6/8)  =  6.4350 11087 93284 38680E-01
  738.         .long   0x3FE4978F,0xA3269EE1
  739. |         ATAN(7/8)  =  7.1882 99996 21624 50542E-01
  740.         .long   0x3FE700A7,0xC5784634
  741. |         ATAN(8/8)  =  7.8539 81633 97448 30962E-01   ( = PI/4)
  742.         .long   0x3FE921FB,0x54442D18
  743. |
  744. |
  745. |  DP PI/2 (duplicate of DPIO2 because of assembler bug)
  746. |
  747. DXPIO2:
  748.         .long   0x3FF921FB,0x54442D18
  749. |
  750. | ###   PUBLIC  DPATN
  751. |
  752. DPATN:
  753.         |dsw    0
  754.         bsr     DPFADJ
  755. |
  756.         jvs     DFNANR          | J/ NaN arg -> NaN result
  757.         jcc     DPAT10          | J/ not INF
  758. |
  759.         movel   DXPIO2+0,d1     | Get top long word of PI/2
  760.         roxll   #1,d1
  761.         roxlw   sp@             | INF sign bit into X
  762.         roxrl   #1,d1           | PI/2 given sign of INF
  763.         addql   #4,sp           | Delete four bytes from the stack
  764.         movel   DXPIO2+4,a7@(4)
  765.         movel   d1,sp@
  766.         jmp     a4@
  767. |
  768. DPAT10:
  769.         jeq     DFZERR          | J/ 0.0 arg -> 0.0 result
  770. |
  771.         bclr    #7,sp@          | Insure argument positive
  772.         sne     d0              | Create flag byte (0xFF iff negative)
  773.         andib   #0x80,d0                | Keep sign bit only
  774.         moveb   d0,a7@(10)      | Save flag byte
  775. |
  776.         movew   sp@,d1
  777.         cmpiw   #16*DBIAS,d1
  778.         jle     DPAT20          | J/ arg < 1 + 1/16
  779. |
  780.         addqb   #1,a7@(10)
  781. |
  782.         clrl    sp@-            | Place 1.0 onstack
  783.         movel   #0x3FF00000,sp@-
  784.         jsr     DPRDIV          | Invert the number
  785. |
  786. DPAT20:
  787.         movew   sp@,d1
  788.         cmpiw   #16*DBIAS-64,d1
  789.         jlt     DPAT30          | J/ arg < 1/16
  790. |
  791.         moveb   d1,d2           | Number of sixteenths in d2
  792.         andib   #0x0F,d2
  793.         orib    #0x10,d2                | Implicit bit
  794.         lsrw    #4,d1
  795.         moveq   #-1,d3
  796.         subb    d1,d3
  797.         lsrb    d3,d2
  798.         addqb   #1,d2
  799.         lsrb    #1,d2           | Rounded eighths in d2.B
  800. |
  801.         clrl    d0
  802.         clrl    d1
  803.         moveb   d2,d1           | 64 bit integer in d0:d1
  804. |
  805.         lslb    #3,d2
  806.         addb    d2,a7@(10)      | Save 8 * eigths on stack
  807. |
  808.         jsr     DFLOAT
  809.         subiw   #16*3,sp@       | Produce y, floating point eighths
  810. |
  811.         moveq   #4-1,d0
  812. DPAT25:
  813.         movel   a7@(12),sp@-
  814.         dbra    d0,DPAT25       | On stack:  y, arg', y, arg', <temps>
  815. |
  816.         jsr     DPMUL           | /* arg'*y */
  817.         clrl    sp@-            | Place 1.0 on stack
  818.         movel   #0x3FF00000,sp@-
  819.         jsr     DPADD           | /* 1.0 + arg'*y */
  820. |
  821.         movel   a7@(16),d0      | Exchange stack item
  822.         movel   a7@(20),d1
  823.         movel   sp@,a7@(16)
  824.         movel   a7@(4),a7@(20)
  825.         movel   d0,sp@
  826.         movel   d1,a7@(4)       | On stack: arg', y, (1+arg'*y), <tmps>
  827.         bset    #7,a7@(8)       | Negate y
  828.         jsr     DPADD           | /* (arg'-y) */
  829.         jsr     DPRDIV          | /* (arg'-y) / (1 + arg'*y) */
  830. |
  831. DPAT30:
  832.         |dsw    0
  833.         movel   a7@(4),sp@-     | Duplicate z
  834.         movel   a7@(4),sp@-
  835. |
  836.         pea     DATNCN
  837.         movel   sp@+,a6 | Polynomial approximation to small ATN
  838.         moveq   #NDATNC,d0
  839.         bsr     DX2SER
  840.         jsr     DPMUL           | Complete approximation
  841. |
  842.         moveb   a7@(10),d0
  843.         andiw   #0x0078,d0      | Trim to table index
  844.         jeq     DPAT40          | J/ y = 0.0, ARCTAN(y) = 0.0
  845. |
  846.         pea     DATNTB-8
  847.         movel   sp@+,a6
  848.         movel   a6@(4,d0:W),sp@-
  849.         movel   a6@(0,d0:W),sp@-
  850.         jsr     DPADD           | Add in ARCTAN(y)
  851. |
  852. DPAT40:
  853.         btst    #0,a7@(10)      | Check for inversion
  854.         jeq     DPAT50          | J/ no inversion
  855. |
  856.         bset    #7,sp@          | Negate ARCTAN
  857.         movel   DXPIO2+4,sp@-
  858.         movel   DXPIO2+0,sp@-
  859.         jsr     DPADD           | Inversion via subtraction
  860. |
  861. DPAT50:
  862.         tstb    a7@(10)         | Check sign of result
  863.         jpl     DPAT60          | J/ positive
  864. |
  865.         bset    #7,sp@          | Negate result
  866. |
  867. DPAT60:
  868.         movel   a7@(4),a7@(8)   | Downshift result
  869.         movel   sp@+,sp@
  870.         jmp     a4@
  872. /*
  873. |
  874. | ###   SUBTTL  DPCOS, DPSIN, DPTAN: Trigonometric Functions
  875. |       page
  876. |
  877. |  TRIG ROUTINES.
  878. |
  879. |  The support routine DTRGSV converts the radian mode argument to
  880. |  an quadrant value between -0.5 and 0.5 (quadrants).  The sign of the
  881. |  is saved (the argument is forced non-negative).  Computations are
  882. |  performed for values quadrant values of -0.5 to 0.5* other values
  883. |  are transformed as per the following transformation formulae
  884. |  and table:
  885. |
  886. |     Let:     z = shifted quadrant value
  887. |     Recall:  sin(-x) = -sin(x)    tan(-x) = -tan(x)
  888. |              cos(-x) =  cos(x)
  889. |
  890. |     l========l=================l=================l=================l
  891. |     l  QUAD  l       SIN       l       COS       l       TAN       l
  892. |     l========l=================l=================l=================l
  893. |     l    0   l      sin(z)     l      cos(z)     l      tan(z)     l
  894. |     l    1   l      cos(z)     l     -sin(z)     l   -1/tan(z)     l
  895. |     l    2   l     -sin(z)     l     -cos(z)     l      tan(z)     l
  896. |     l    3   l     -cos(z)     l      sin(z)     l   -1/tan(z)     l
  897. |     l========l=================l=================l=================l
  898. |
  899. |  The sign of the argument is, by the equations above, important only
  900. |  when computing SIN and TAN.
  901. |
  902. |  Chebyshev Polynomials are used to approximate the trigonometric
  903. |  value in the range -0.5 to 0.5 quandrants.  (Note: constants are
  904. |  adjusted to reflect calculations done in quadrant mode instead
  905. |  of radian mode).
  906. |
  907. |
  908. |
  909. |         C16  =   6.5659 63114 97947 23622E-11
  910. */
  912. DCOSCN:
  913.         .long   0x3DD20C62,0xC2F2D7F5
  914. |         C14  =  -6.3866 03083 79185 22411E-09
  915.         .long   0xBE3B6E24,0xF44B128F
  916. |         C12  =   4.7108 74778 81817 15037E-07
  917.         .long   0x3E9F9D38,0xA3763CC3
  918. |         C10  =  -2.5202 04237 30606 05481E-05
  919.         .long   0xBEFA6D1F,0x2A204A8C
  920. |         C08  =   9.1926 02748 39426 58024E-04
  921.         .long   0x3F4E1F50,0x6891BABB
  922. |         C06  =  -2.0863 48076 33529 96087E-02
  923.         .long   0xBF955D3C,0x7E3CBFFA
  924. |         C04  =   2.5366 95079 01048 01364E-01
  925.         .long   0x3FD03C1F,0x081B5AC4
  926. |         C02  =  -1.2337 00550 13616 98274E+00
  927.         .long   0xBFF3BD3C,0xC9BE45DE
  928. |         C00  =   1.0000 00000 00000 00000E+00
  929.         .long   0x3FF00000,0x00000000
  930. |
  931. |       .set    NDCOSC,9
  932. |
  933. |
  934. |
  935. |         C15  =  -6.6880 35109 81146 72325E-10
  936. DSINCN:
  937.         .long   0xBE06FADB,0x9F155744
  938. |         C13  =   5.6921 72921 96792 68118E-08
  939.         .long   0x3E6E8F43,0x4D018D63
  940. |         C11  =  -3.5988 43235 21208 53405E-06
  941.         .long   0xBECE3074,0xFDE8871F
  942. |         C09  =   1.6044 11847 87359 82187E-04
  943.         .long   0x3F250783,0x487EE782
  944. |         C07 =   -4.6817 54135 31868 81007E-03
  945.         .long   0xBF732D2C,0xCE62BD86
  946. |         C05  =   7.9692 62624 61670 45121E-02
  947.         .long   0x3FB466BC,0x6775AAE2
  948. |         C03  =  -6.4596 40975 06246 25366E-01
  949.         .long   0xBFE4ABBC,0xE625BE53
  950. |         C01  =   1.5707 96326 79489 66192E+00
  951. DPIO2:
  952.         .long   0x3FF921FB,0x54442D18
  953. |
  954. |       .set    NDSINC,8
  955. |
  956. |
  957. |
  958. |         Q3  =    1.7732 32244 08118 84863E+01
  959. DTNQCN:
  960.         .long   0x4031BB79,0x7BC569F8
  961. |         Q2  =   -8.0485 82645 07638 77427E+02
  962.         .long   0xC08926DD,0xB9C83D02
  963. |         Q1  =    6.4501 55566 23337 83845E+03
  964.         .long   0x40B93227,0xD3304CB9
  965. |         Q0  =   -5.6630 29630 56808 22115E+03
  966.         .long   0xC0B61F07,0x95DE70E0
  967. |
  968. |       .set    NDTNQC,4
  969. |
  970. |
  971. |         P3  =    1.0000 00000 00000 00000E+00
  972. DTNPCN:
  973.         .long   0x3FF00000,0x00000000
  974. |         P2  =   -1.5195 78275 66504 79235E+02
  975.         .long   0xC062FEA6,0x85FF2B0D
  976. |         P1  =    2.8156 53021 77302 44048E+03
  977.         .long   0x40A5FF4E,0x58DEAD6E
  978. |         P0  =   -8.8954 66142 22700 41047E+03
  979.         .long   0xC0C15FBB,0xAA8C6A22
  980. |
  981. |       .set    NDTNPC,4
  982. |
  983. |
  984. |         INV2PI  =  1 / 2*PI  =   1.5915 49430 91895 33577E-01
  985. DIN2PI:
  986.         .long   0x3FC45F30,0x6DC9C883
  987. |
  988. |
  989. | ###   PUBLIC  DPCOS
  990. |
  991. DPCOS:
  992.         bsr     DTRGSV          | Prepare argument for operation
  993. |
  994.         moveb   d0,d1           | Copy into d1
  995.         andib   #3,d0           | Strip sign bit
  996.         addqb   #1,d0
  997.         rorb    #2,d0           | Move sign bit to B7
  998.         moveb   d0,a7@(8)       | Save it
  999.         rorb    #1,d1
  1000.         jcs     DSINOP          | Compute sine
  1001. |
  1002. |
  1003. DCOSOP:
  1004.         pea     DCOSCN
  1005.         movel   sp@+,a6
  1006.         moveq   #NDCOSC,d0
  1007.         bsr     DX2SER
  1008. |
  1009. DTRGFN:
  1010.         moveb   a7@(8),d0
  1011.         andib   #0x80,d0                | Isolate sign bit
  1012.         eorb    d0,sp@
  1013. |
  1014.         movel   a7@(4),a7@(8)
  1015.         movel   sp@+,sp@
  1016.         jmp     a4@
  1017. |
  1018. |
  1019. |
  1020. | ###   PUBLIC  DPSIN
  1021. |
  1022. DPSIN:
  1023.         bsr     DTRGSV          | Prepare argument
  1024. |
  1025.         addib   #0x7E,a7@(8)    | Compute sign
  1026.         rorb    #1,d0
  1027.         jcs     DCOSOP
  1028. |
  1029. |
  1030. DSINOP:
  1031.         movel   a7@(4),sp@-     | Duplicate argument
  1032.         movel   a7@(4),sp@-
  1033. |
  1034.         pea     DSINCN
  1035.         movel   sp@+,a6
  1036.         moveq   #NDSINC,d0
  1037.         bsr     DX2SER
  1038.         jsr     DPMUL
  1039. |
  1040.         jra     DTRGFN
  1041. |
  1042. |
  1043. |
  1044. | ###   PUBLIC  DPTAN
  1045. |
  1046. DPTAN:
  1047.         bsr     DTRGSV          | Prepare argument
  1048. |
  1049.         rorb    #1,d0
  1050.         andib   #0x80,d0
  1051.         eorb    d0,a7@(8)
  1052. |
  1053.         moveq   #4-1,d1         | Double duplication of argument
  1054. DPTN01:
  1055.         movel   a7@(4),sp@-
  1056.         dbra    d1,DPTN01
  1057. |
  1058.         pea     DTNPCN
  1059.         movel   sp@+,a6
  1060.         moveq   #NDTNPC,d0
  1061.         bsr     DX2SER
  1062.         jsr     DPMUL
  1063. |
  1064.         movel   a7@(8),d0
  1065.         movel   a7@(12),d1
  1066.         movel   sp@,a7@(8)
  1067.         movel   a7@(4),a7@(12)
  1068.         movel   d0,sp@
  1069.         movel   d1,a7@(4)
  1070. |
  1071.         pea     DTNQCN
  1072.         movel   sp@+,a6
  1073.         moveq   #NDTNQC,d0
  1074.         bsr     DX2SER
  1075. |
  1076.         btst    #0,a7@(16)
  1077.         jne     DPTN20          | J/ reverse division
  1078. |
  1079.         jsr     DPDIV
  1080.         jra     DTRGFN
  1081. |
  1082. DPTN20:
  1083.         jsr     DPRDIV
  1084.         jra     DTRGFN
  1085. /*
  1086. |
  1087. |       page
  1088. |
  1089. |  Trigonometric Service Routine
  1090. |
  1091. |  The double precision floating point value to be processed is on the
  1092. |  stack.  Excess 2 PI's are removed from the argument.  The range of
  1093. |  the argument is then scaled to within PI/4 of 0.  The shift size
  1094. |  (in number of quadrants) is returned in D0.  The original sign bit
  1095. |  is returned in the D0.B sign bit.
  1096. |
  1097. |  If the argument is NaN, +INF, -INF, or too large (>= 2^16), this
  1098. |  routine causes a return to the caller with a NaN result.
  1099. |
  1100. */
  1102. DTRGSV:
  1103.         moveal  sp@+,a6 | DTRGSV return addr
  1104. |
  1105.         bsr     DPFADJ
  1106.         jcs     DFNANR          | Return NAN if +/- INF
  1107.         jvs     DFNANR          | Return NaN if NaN
  1108. |
  1109.         clrw    a7@(8)          | Create arg type return byte
  1110. |
  1111.         roxlw   sp@             | Sign bit into X
  1112.         roxrw   a7@(8)          | Sign bit into temp byte
  1113.         roxrw   sp@             | Restore DP value (with sign = +)
  1114. |
  1115.         cmpiw   #16*DBIAS+256,sp@
  1116.         jge     DFNANR          | J/ return NaN if ABS(arg) >= 2^16
  1117. |
  1118.         movel   DIN2PI+4,sp@-   | Scale arg by 1/(2*PI)
  1119.         movel   DIN2PI+0,sp@-
  1120.         jsr     DPMUL
  1121. |
  1122.         cmpiw   #16*DBIAS,sp@
  1123.         jlt     DTS10           | /* J/ no excess two PI's */
  1124. |
  1125.         movel   a7@(4),sp@-     | Double the scaled argument
  1126.         movel   a7@(4),sp@-
  1127.         jsr     DINT
  1128.         negl    d1
  1129.         negxl   d0
  1130.         jsr     DFLOAT
  1131.         jsr     DPADD
  1132. |
  1133. DTS10:
  1134.         tstw    sp@
  1135.         jeq     DTS20
  1136. |
  1137.         addiw   #16*3,sp@       | Multiply value by 8
  1138.         cmpiw   #16*DBIAS,sp@
  1139.         jlt     DTS20           | J/ < 1
  1140. |
  1141.         movel   a7@(4),sp@-     | Double parameter
  1142.         movel   a7@(4),sp@-
  1143.         jsr     DINT
  1144.         addqb   #1,d1
  1145.         lsrb    #1,d1
  1146.         addb    d1,a7@(8)       | Mix quadrant number with sign bit
  1147. |
  1148.         lslb    #1,d1
  1149.         negl    d1
  1150.         negxl   d0
  1151.         jsr     DFLOAT
  1152.         jsr     DPADD
  1153. |
  1154. DTS20:
  1155.         tstw    sp@
  1156.         jeq     DTS30           | J/ arg = 0.0
  1157. |
  1158.         subiw   #16*1,sp@       | Range result to -0.5 to 0.5
  1159. |
  1160. DTS30:
  1161.         andib   #0x83,a7@(8)    | Limit quadrant shift to 0 thru 3
  1162.         moveb   a7@(8),d0       | Return cntl byte in d0 and at a7@(8)
  1163.         jmp     a6@
  1164. /*
  1165. |
  1166. | ###   SUBTTL  Polynomial Evaluation Routine
  1167. |       page
  1168. |
  1169. |  Polynomial Evaluation Routines.
  1170. |
  1171. |  A list of constants is pointed to by A6.  X is on the stack.  D0
  1172. |  contains the number of constants (the polynomial degree plus one).
  1173. |
  1174. |  A6 enters pointing to C[1].  Upon return, the value on stack has
  1175. |  been replaced with:
  1176. |  DXSER   = C[1]*X^(N-1)  + C[2]*X^(N-2)  + ... + C[N-1]*X   + C[N]
  1177. |
  1178. |  DX2SER  = C[1]*X^(2N-2) + C[2]*X^(2N-4) + ... + C[N-1]*X^2 + C[N]
  1179. |
  1180. |
  1181. */
  1183. DX2SER:
  1184.         moveal  a4,a5           | Save a4 in a5
  1185.         bsr     DPFADJ
  1186. |
  1187.         moveb   d0,a7@(8)       | Save count
  1188. |
  1189.         movel   a7@(4),sp@-     | Square the argument
  1190.         movel   a7@(4),sp@-
  1191.         jsr     DPMUL
  1192.         jra     DXSR01          | J/ join DXSER routine
  1193. |
  1194. DXSER:
  1195.         moveal  a4,a5
  1196.         bsr     DPFADJ
  1197. |
  1198.         moveb   d0,a7@(8)
  1199. |
  1200. DXSR01:
  1201.         movel   a6@(4),sp@-     | Place C[1] on stack as accum
  1202.         movel   a6@,sp@-
  1203.         addql   #8,a6
  1204. |
  1205.         subqb   #1,a7@(16)
  1206. |
  1207. DXSR02:
  1208.         movel   a7@(12),sp@-    | Copy X
  1209.         movel   a7@(12),sp@-
  1210.         jsr     DPMUL
  1211.         movel   a6@(4),sp@-     | Add in next coefficient
  1212.         movel   a6@,sp@-
  1213.         addql   #8,a6           | Advance to next coefficient
  1214.         jsr     DPADD
  1215. |
  1216.         subqb   #1,a7@(16)
  1217.         jne     DXSR02          | J/ more coefficients
  1218. |
  1219.         movel   a7@(4),a7@(16)  | Adjust return
  1220.         movel   sp@+,a7@(8)
  1221.         addql   #8,sp
  1222. |
  1223.         moveal  a4,a3
  1224.         moveal  a5,a4           | Restore return address
  1225.         jmp     a3@
  1226. /*
  1227. |
  1228. | ###   SUBTTL  Function Completion Segments
  1229. |       page
  1230. |
  1231. |  DPFADJ - Adjust the stack for double precision functions.
  1232. |          Move the return address to the function caller into A4.
  1233. |          Shift the double precision floating point argument down
  1234. |          four bytes, leaving a temporary area of four bytes at
  1235. |          8(A7) thru 11(A7).  Set the condition code bits N, Z,
  1236. |          C (INF), and V (NaN) to quickly type the argument.
  1237. |
  1238. */
  1240. DPFADJ:
  1241.         moveal  sp@+,a3 | Function return addr
  1242.         moveal  sp@,a4          | Caller return addr
  1243.         movel   a7@(4),sp@      | Down shift parameter
  1244.         movel   a7@(8),a7@(4)
  1245. |
  1246.         movew   sp@,d2          | Get SEM word
  1247.         rolw    #1,d2
  1248.         cmpiw   #32*0x7FF,d2
  1249.         jcs     DPFA10          | J/ not INF or NaN
  1250. |
  1251.         andiw   #0x1E,d2
  1252.         jeq     DPFA05
  1253. |
  1254.         orib    #0x02,ccr       | Set -V- bit -> NaN argument
  1255.         jmp     a3@
  1256. |
  1257. DPFA05:
  1258.         tstb    sp@             | /* Set N bit as req'd, clear Z bit. */
  1259.         orib    #0x01,ccr       | Set -C- bit -> INF argument
  1260.         jmp     a3@
  1261. |
  1262. DPFA10:
  1263.         tstw    sp@             | /* Set N, Z bits as req'd, clear V, C */
  1264.         jmp     a3@
  1265. |
  1266. |
  1267. |
  1268. DFPINR:
  1269.         |dsw    0               | Result is positive overflow
  1270.         moveq   #0,d0           | Sign is positive
  1271.         jra     DFIN01
  1272. |
  1273. DFMINR:
  1274.         |dsw    0               | Result in negative overflow
  1275.         moveq   #-1,d0          | Sign is negative
  1276. |
  1277. DFIN01:
  1278.         moveaw  d0,a2           | Sign into a2
  1279.         addql   #8,sp           | Delete twelve bytes from the stack
  1280.         addql   #4,sp
  1281.         moveal  a4,a0           | Return addr in a0
  1282.         jmp     DINFRS          | (Use DPOPNS exception return)
  1283. |
  1284. |
  1285. DFNANR:
  1286.         |dsw    0               | Result is NaN
  1287.         addql   #8,sp           | Delete twelve bytes from the stack
  1288.         addql   #4,sp
  1289.         moveal  a4,a0
  1290.         jmp     DNANRS
  1291. |
  1292. |
  1293. DFZERR:
  1294.         |dsw    0               | Result is zero
  1295.         addql   #8,sp
  1296.         addql   #4,sp
  1297.         moveal  a4,a0
  1298.         jmp     DZERRS
  1299. |
  1300. |
  1301. DFUNFR:
  1302.         |dsw    0               | Underflow result
  1303.         addql   #8,sp           | Delete twelve bytes from the stack
  1304.         addql   #4,sp
  1305.         moveal  a4,a0           | Return address into a0
  1306.         jmp     DUNFRS
  1307. |
  1308. |
  1309. DFONER:
  1310.         |dsw    0               | Result is 1.0
  1311.         addql   #4,sp
  1312.         movel   #0x3FF00000,sp@
  1313.         clrl    a7@(4)
  1314.         clrb    FPERR
  1315.         jmp     a4@
  1316. |
  1317. |       end