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

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

VxWorks

开发平台：

C/C++

l_stan.s：源码内容
  1. /* l_stan.s - Motorola 68040 FP tangent routines (LIB) */
  3. /* Copyright 1991-1993 Wind River Systems, Inc. */
  4. .data
  5. .globl _copyright_wind_river
  6. .long _copyright_wind_river
  8. /*
  9. modification history
  10. --------------------
  11. 01f,21jul93,kdl  added .text (SPR #2372).
  12. 01e,23aug92,jcf  changed bxxx to jxx.
  13. 01d,26may92,rrr  the tree shuffle
  14. 01c,10jan92,kdl  general cleanup.
  15. 01b,17dec91,kdl  put in changes from Motorola v3.3 (from FPSP 2.1):
  16.  reduce argument by one step before general reduction
  17.  loop.
  18. 01a,15aug91,kdl  original version, from Motorola FPSP v2.0.
  19. */
  21. /*
  22. DESCRIPTION
  24. stansa 3.2 12/18/90
  26. The entry point __l_stan computes the tangent of
  27. an input argument;
  28. __l_stand does the same except for denormalized input.
  30. Input: Double-extended number X in location pointed to
  31. by address register a0.
  33. Output: The value tan(X) returned in floating-point register Fp0.
  35. Accuracy and Monotonicity: The returned result is within 3 ulp in
  36. 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
  37. result is subsequently rounded to double precision. The
  38. result is provably monotonic in double precision.
  40. Speed: The program sTAN takes approximately 170 cycles for
  41. input argument X such that |X| < 15Pi, which is the the usual
  42. situation.
  44. Algorithm:
  46. 1. If |X| >= 15Pi or |X| < 2**(-40), go to 6.
  48. 2. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
  49. k = N mod 2, so in particular, k = 0 or 1.
  51. 3. If k is odd, go to 5.
  53. 4. (k is even) Tan(X) = tan(r) and tan(r) is approximated by a
  54. rational function U/V where
  55. U = r + r*s*(P1 + s*(P2 + s*P3)), and
  56. V = 1 + s*(Q1 + s*(Q2 + s*(Q3 + s*Q4))),  s = r*r.
  57. Exit.
  59. 5. (k is odd) Tan(X) = -cot(r). Since tan(r) is approximated by a
  60. rational function U/V where
  61. U = r + r*s*(P1 + s*(P2 + s*P3)), and
  62. V = 1 + s*(Q1 + s*(Q2 + s*(Q3 + s*Q4))), s = r*r,
  63. -Cot(r) = -V/U. Exit.
  65. 6. If |X| > 1, go to 8.
  67. 7. (|X|<2**(-40)) Tan(X) = X. Exit.
  69. 8. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 2.
  72. Copyright (C) Motorola, Inc. 1990
  73. All Rights Reserved
  75. THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF MOTOROLA
  76. The copyright notice above does not evidence any
  77. actual or intended publication of such source code.
  79. STAN idnt 2,1 Motorola 040 Floating Point Software Package
  81. section 8
  83. NOMANUAL
  84. */
  86. #include "fpsp040L.h"
  88. BOUNDS1: .long 0x3FD78000,0x4004BC7E
  89. TWOBYPI: .long 0x3FE45F30,0x6DC9C883
  91. TANQ4: .long 0x3EA0B759,0xF50F8688
  92. TANP3: .long 0xBEF2BAA5,0xA8924F04
  94. TANQ3: .long 0xBF346F59,0xB39BA65F,0x00000000,0x00000000
  96. TANP2: .long 0x3FF60000,0xE073D3FC,0x199C4A00,0x00000000
  98. TANQ2: .long 0x3FF90000,0xD23CD684,0x15D95FA1,0x00000000
  100. TANP1: .long 0xBFFC0000,0x8895A6C5,0xFB423BCA,0x00000000
  102. TANQ1: .long 0xBFFD0000,0xEEF57E0D,0xA84BC8CE,0x00000000
  104. INVTWOPI: .long 0x3FFC0000,0xA2F9836E,0x4E44152A,0x00000000
  106. TWOPI1: .long 0x40010000,0xC90FDAA2,0x00000000,0x00000000
  107. TWOPI2: .long 0x3FDF0000,0x85A308D4,0x00000000,0x00000000
  109. |--N*PI/2, -32 <= N <= 32, IN A LEADING TERM IN EXT. AND TRAILING
  110. |--TERM IN SGL. NOTE THAT PI IS 64-BIT LONG, THUS N*PI/2 IS AT
  111. |--MOST 69 BITS LONG.
  112. .globl __l_PITBL
  113. __l_PITBL:
  114.   .long  0xC0040000,0xC90FDAA2,0x2168C235,0x21800000
  115.   .long  0xC0040000,0xC2C75BCD,0x105D7C23,0xA0D00000
  116.   .long  0xC0040000,0xBC7EDCF7,0xFF523611,0xA1E80000
  117.   .long  0xC0040000,0xB6365E22,0xEE46F000,0x21480000
  118.   .long  0xC0040000,0xAFEDDF4D,0xDD3BA9EE,0xA1200000
  119.   .long  0xC0040000,0xA9A56078,0xCC3063DD,0x21FC0000
  120.   .long  0xC0040000,0xA35CE1A3,0xBB251DCB,0x21100000
  121.   .long  0xC0040000,0x9D1462CE,0xAA19D7B9,0xA1580000
  122.   .long  0xC0040000,0x96CBE3F9,0x990E91A8,0x21E00000
  123.   .long  0xC0040000,0x90836524,0x88034B96,0x20B00000
  124.   .long  0xC0040000,0x8A3AE64F,0x76F80584,0xA1880000
  125.   .long  0xC0040000,0x83F2677A,0x65ECBF73,0x21C40000
  126.   .long  0xC0030000,0xFB53D14A,0xA9C2F2C2,0x20000000
  127.   .long  0xC0030000,0xEEC2D3A0,0x87AC669F,0x21380000
  128.   .long  0xC0030000,0xE231D5F6,0x6595DA7B,0xA1300000
  129.   .long  0xC0030000,0xD5A0D84C,0x437F4E58,0x9FC00000
  130.   .long  0xC0030000,0xC90FDAA2,0x2168C235,0x21000000
  131.   .long  0xC0030000,0xBC7EDCF7,0xFF523611,0xA1680000
  132.   .long  0xC0030000,0xAFEDDF4D,0xDD3BA9EE,0xA0A00000
  133.   .long  0xC0030000,0xA35CE1A3,0xBB251DCB,0x20900000
  134.   .long  0xC0030000,0x96CBE3F9,0x990E91A8,0x21600000
  135.   .long  0xC0030000,0x8A3AE64F,0x76F80584,0xA1080000
  136.   .long  0xC0020000,0xFB53D14A,0xA9C2F2C2,0x1F800000
  137.   .long  0xC0020000,0xE231D5F6,0x6595DA7B,0xA0B00000
  138.   .long  0xC0020000,0xC90FDAA2,0x2168C235,0x20800000
  139.   .long  0xC0020000,0xAFEDDF4D,0xDD3BA9EE,0xA0200000
  140.   .long  0xC0020000,0x96CBE3F9,0x990E91A8,0x20E00000
  141.   .long  0xC0010000,0xFB53D14A,0xA9C2F2C2,0x1F000000
  142.   .long  0xC0010000,0xC90FDAA2,0x2168C235,0x20000000
  143.   .long  0xC0010000,0x96CBE3F9,0x990E91A8,0x20600000
  144.   .long  0xC0000000,0xC90FDAA2,0x2168C235,0x1F800000
  145.   .long  0xBFFF0000,0xC90FDAA2,0x2168C235,0x1F000000
  146.   .long  0x00000000,0x00000000,0x00000000,0x00000000
  147.   .long  0x3FFF0000,0xC90FDAA2,0x2168C235,0x9F000000
  148.   .long  0x40000000,0xC90FDAA2,0x2168C235,0x9F800000
  149.   .long  0x40010000,0x96CBE3F9,0x990E91A8,0xA0600000
  150.   .long  0x40010000,0xC90FDAA2,0x2168C235,0xA0000000
  151.   .long  0x40010000,0xFB53D14A,0xA9C2F2C2,0x9F000000
  152.   .long  0x40020000,0x96CBE3F9,0x990E91A8,0xA0E00000
  153.   .long  0x40020000,0xAFEDDF4D,0xDD3BA9EE,0x20200000
  154.   .long  0x40020000,0xC90FDAA2,0x2168C235,0xA0800000
  155.   .long  0x40020000,0xE231D5F6,0x6595DA7B,0x20B00000
  156.   .long  0x40020000,0xFB53D14A,0xA9C2F2C2,0x9F800000
  157.   .long  0x40030000,0x8A3AE64F,0x76F80584,0x21080000
  158.   .long  0x40030000,0x96CBE3F9,0x990E91A8,0xA1600000
  159.   .long  0x40030000,0xA35CE1A3,0xBB251DCB,0xA0900000
  160.   .long  0x40030000,0xAFEDDF4D,0xDD3BA9EE,0x20A00000
  161.   .long  0x40030000,0xBC7EDCF7,0xFF523611,0x21680000
  162.   .long  0x40030000,0xC90FDAA2,0x2168C235,0xA1000000
  163.   .long  0x40030000,0xD5A0D84C,0x437F4E58,0x1FC00000
  164.   .long  0x40030000,0xE231D5F6,0x6595DA7B,0x21300000
  165.   .long  0x40030000,0xEEC2D3A0,0x87AC669F,0xA1380000
  166.   .long  0x40030000,0xFB53D14A,0xA9C2F2C2,0xA0000000
  167.   .long  0x40040000,0x83F2677A,0x65ECBF73,0xA1C40000
  168.   .long  0x40040000,0x8A3AE64F,0x76F80584,0x21880000
  169.   .long  0x40040000,0x90836524,0x88034B96,0xA0B00000
  170.   .long  0x40040000,0x96CBE3F9,0x990E91A8,0xA1E00000
  171.   .long  0x40040000,0x9D1462CE,0xAA19D7B9,0x21580000
  172.   .long  0x40040000,0xA35CE1A3,0xBB251DCB,0xA1100000
  173.   .long  0x40040000,0xA9A56078,0xCC3063DD,0xA1FC0000
  174.   .long  0x40040000,0xAFEDDF4D,0xDD3BA9EE,0x21200000
  175.   .long  0x40040000,0xB6365E22,0xEE46F000,0xA1480000
  176.   .long  0x40040000,0xBC7EDCF7,0xFF523611,0x21E80000
  177.   .long  0x40040000,0xC2C75BCD,0x105D7C23,0x20D00000
  178.   .long  0x40040000,0xC90FDAA2,0x2168C235,0xA1800000
  180. #define INARG FP_SCR4
  182. #define TWOTO63      L_SCR1
  183. #define ENDFLAG L_SCR2
  184. #define N            L_SCR3
  186. | xref __l_t_frcinx
  187. | xref __l_t_extdnrm
  190. .text
  192. .globl __l_stand
  193. __l_stand:
  194. |--TAN(X) = X FOR DENORMALIZED X
  196. jra  __l_t_extdnrm
  198. .globl __l_stan
  199. __l_stan:
  200. fmovex a0@,fp0 |...lOAD INPUT
  202. movel A0@,d0
  203. movew A0@(4),d0
  204. andil #0x7FFFFFFF,d0
  206. cmpil #0x3FD78000,d0 |...|X| >= 2**(-40)?
  207. jge  TANOK1
  208. jra  TANSM
  209. TANOK1:
  210. cmpil #0x4004BC7E,d0 |...|X| < 15 PI?
  211. jlt  TANMAIN
  212. jra  REDUCEX
  215. TANMAIN:
  216. |--THIS IS THE USUAL CASE, |X| <= 15 PI.
  217. |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
  218. fmovex fp0,fp1
  219. fmuld TWOBYPI,fp1 |...X*2/PI
  221. |--HIDE THE NEXT TWO INSTRUCTIONS
  222. lea __l_PITBL+0x200,a1 |...TABLE OF N*PI/2, N = -32,...,32
  224. |--FP1 IS NOW READY
  225. fmovel fp1,d0 |...CONVERT TO INTEGER
  227. asll #4,d0
  228. addal d0,a1 |...ADDRESS N*PIBY2 IN Y1, Y2
  230. fsubx a1@+,fp0 |...X-Y1
  231. |--HIDE THE NEXT ONE
  233. fsubs a1@,fp0 |...FP0 IS R = (X-Y1)-Y2
  235. rorl #5,d0
  236. andil #0x80000000,d0 |...D0 WAS ODD IFF d0 < 0
  238. TANCONT:
  240. cmpil #0,d0
  241. jlt  NODD
  243. fmovex fp0,fp1
  244. fmulx fp1,fp1   |...S = R*R
  246. fmoved TANQ4,fp3
  247. fmoved TANP3,fp2
  249. fmulx fp1,fp3   |...SQ4
  250. fmulx fp1,fp2   |...SP3
  252. faddd TANQ3,fp3 |...Q3+SQ4
  253. faddx TANP2,fp2 |...P2+SP3
  255. fmulx fp1,fp3   |...S(Q3+SQ4)
  256. fmulx fp1,fp2   |...S(P2+SP3)
  258. faddx TANQ2,fp3 |...Q2+S(Q3+SQ4)
  259. faddx TANP1,fp2 |...P1+S(P2+SP3)
  261. fmulx fp1,fp3   |...S(Q2+S(Q3+SQ4))
  262. fmulx fp1,fp2   |...S(P1+S(P2+SP3))
  264. faddx TANQ1,fp3 |...Q1+S(Q2+S(Q3+SQ4))
  265. fmulx fp0,fp2   |...RS(P1+S(P2+SP3))
  267. fmulx fp3,fp1   |...S(Q1+S(Q2+S(Q3+SQ4)))
  270. faddx fp2,fp0   |...R+RS(P1+S(P2+SP3))
  273. .long 0xf23c44a2,0x3f800000 /* fadds &0x3F800000,fp1 */
  275. fmovel d1,fpcr | restore users exceptions
  276. fdivx fp1,fp0 | last inst - possible exception set
  278. jra  __l_t_frcinx
  280. NODD:
  281. fmovex fp0,fp1
  282. fmulx fp0,fp0   |...S = R*R
  284. fmoved TANQ4,fp3
  285. fmoved TANP3,fp2
  287. fmulx fp0,fp3   |...SQ4
  288. fmulx fp0,fp2   |...SP3
  290. faddd TANQ3,fp3 |...Q3+SQ4
  291. faddx TANP2,fp2 |...P2+SP3
  293. fmulx fp0,fp3   |...S(Q3+SQ4)
  294. fmulx fp0,fp2   |...S(P2+SP3)
  296. faddx TANQ2,fp3 |...Q2+S(Q3+SQ4)
  297. faddx TANP1,fp2 |...P1+S(P2+SP3)
  299. fmulx fp0,fp3   |...S(Q2+S(Q3+SQ4))
  300. fmulx fp0,fp2   |...S(P1+S(P2+SP3))
  302. faddx TANQ1,fp3 |...Q1+S(Q2+S(Q3+SQ4))
  303. fmulx fp1,fp2   |...RS(P1+S(P2+SP3))
  305. fmulx fp3,fp0   |...S(Q1+S(Q2+S(Q3+SQ4)))
  308. faddx fp2,fp1   |...R+RS(P1+S(P2+SP3))
  309. .long 0xf23c4422,0x3f800000 /* fadds  &0x3F800000,fp0 */
  312. fmovex fp1,a7@-
  313. eoril #0x80000000,a7@
  315. fmovel d1,fpcr   | restore users exceptions
  316. fdivx a7@+,fp0 | last inst - possible exception set
  318. jra  __l_t_frcinx
  320. TANBORS:
  321. |--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION.
  322. |--IF |X| < 2**(-40), RETURN X OR 1.
  323. cmpil #0x3FFF8000,d0
  324. jgt  REDUCEX
  326. TANSM:
  328. fmovex fp0,a7@-
  329. fmovel d1,fpcr  | restore users exceptions
  330. fmovex a7@+,fp0 | last inst - posibble exception set
  332. jra  __l_t_frcinx
  335. REDUCEX:
  336. |--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW.
  337. |--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING
  338. |--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE.
  340. fmovemx fp2-fp5,A7@- |...save fp2 through fp5
  341. movel d2,A7@-
  342. .long 0xf23c4480,0x00000000 /* fmoves  &0x00000000,fp1 */
  344. |--If compact form of abs(arg) in d0=0x7ffeffff, argument is so large that
  345. |--there is a danger of unwanted overflow in first LOOP iteration.  In this
  346. |--case, reduce argument by one remainder step to make subsequent reduction
  347. |--safe.
  348. cmpil #0x7ffeffff,d0 | is argument dangerously large?
  349. jne  LOOP
  350. movel #0x7ffe0000,a6@(FP_SCR2) | yes
  351. | | create 2**16383*PI/2
  352. movel #0xc90fdaa2,a6@(FP_SCR2+4)
  353. clrl a6@(FP_SCR2+8)
  354. ftstx fp0 | test sign of argument
  355. movel #0x7fdc0000,a6@(FP_SCR3) | create low half of 2**16383*
  356. | | PI/2 at FP_SCR3
  357. movel #0x85a308d3,a6@(FP_SCR3+4)
  358. clrl   a6@(FP_SCR3+8)
  359. fblt red_neg
  360. orw #0x8000,a6@(FP_SCR2) | positive arg
  361. orw #0x8000,a6@(FP_SCR3)
  362. red_neg:
  363. faddx  a6@(FP_SCR2),fp0 | high part of reduction is exact
  364. fmovex  fp0,fp1 | save high result in fp1
  365. faddx  a6@(FP_SCR3),fp0 | low part of reduction
  366. fsubx   fp0,fp1 | determine low component of result
  367. faddx  a6@(FP_SCR3),fp1 | fp0/fp1 are reduced argument.
  369. |--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4.
  370. |--integer quotient will be stored in N
  371. |--Intermeditate remainder is 66-bit long|  (R,r) in (FP0,FP1)
  373. LOOP:
  374. fmovex fp0,a6@(INARG) |...+-2**K * F, 1 <= F < 2
  375. movew a6@(INARG),d0
  376.         movel          d0,a1 |...save a copy of d0
  377. andil #0x00007FFF,d0
  378. subil #0x00003FFF,d0 |...D0 IS K
  379. cmpil #28,d0
  380. jle  LASTLOOP
  381. CONTLOOP:
  382. subil #27,d0   |...D0 IS L := K-27
  383. movel #0,a6@(ENDFLAG)
  384. jra  WORK
  385. LASTLOOP:
  386. clrl d0 |...D0 IS L := 0
  387. movel #1,a6@(ENDFLAG)
  389. WORK:
  390. |--FIND THE REMAINDER OF (R,r) W.R.T. 2**L * (PI/2). L IS SO CHOSEN
  391. |--THAT INT( X * (2/PI) / 2**(L) ) < 2**29.
  393. |--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63),
  394. |--2**L * (PIby2_1), 2**L * (PIby2_2)
  396. movel #0x00003FFE,d2 |...BIASED EXPO OF 2/PI
  397. subl d0,d2 |...BIASED EXPO OF 2**(-L)*(2/PI)
  399. movel #0xA2F9836E,a6@(FP_SCR1+4)
  400. movel #0x4E44152A,a6@(FP_SCR1+8)
  401. movew d2,a6@(FP_SCR1) |...FP_SCR1 is 2**(-L)*(2/PI)
  403. fmovex fp0,fp2
  404. fmulx a6@(FP_SCR1),fp2
  405. |--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN
  406. /* |--FLOATING POINT FORMAT, THE TWO FMOVE'S FMOVE.l FP <--> N */
  407. |--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT
  408. |--(SIGN(INARG)*2**63 + FP2) - SIGN(INARG)*2**63 WILL GIVE
  409. |--US THE DESIRED VALUE IN FLOATING POINT.
  411. |--HIDE SIX CYCLES OF INSTRUCTION
  412.         movel a1,d2
  413.         swap d2
  414. andil #0x80000000,d2
  415. oril #0x5F000000,d2 |...D2 IS SIGN(INARG)*2**63 IN SGL
  416. movel d2,a6@(TWOTO63)
  418. movel d0,d2
  419. addil #0x00003FFF,d2 |...BIASED EXPO OF 2**L * (PI/2)
  421. |--FP2 IS READY
  422. fadds a6@(TWOTO63),fp2 |...THE FRACTIONAL PART OF fp1 IS ROUNDED
  424. |--HIDE 4 CYCLES OF INSTRUCTION|  creating 2**(L)*Piby2_1  and  2**(L)*Piby2_2
  425.         movew d2,a6@(FP_SCR2)
  426. clrw           a6@(FP_SCR2+2)
  427. movel #0xC90FDAA2,a6@(FP_SCR2+4)
  428. clrl a6@(FP_SCR2+8) |...FP_SCR2 is  2**(L) * Piby2_1
  430. |--FP2 IS READY
  431. fsubs a6@(TWOTO63),fp2 |...FP2 is N
  433. addil #0x00003FDD,d0
  434.         movew d0,a6@(FP_SCR3)
  435. clrw           a6@(FP_SCR3+2)
  436. movel #0x85A308D3,a6@(FP_SCR3+4)
  437. clrl a6@(FP_SCR3+8) |...FP_SCR3 is 2**(L) * Piby2_2
  439. movel a6@(ENDFLAG),d0
  441. |--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and
  442. |--P2 = 2**(L) * Piby2_2
  443. fmovex fp2,fp4
  444. fmulx a6@(FP_SCR2),fp4 |...w = N*P1
  445. fmovex fp2,fp5
  446. fmulx a6@(FP_SCR3),fp5 |...w = N*P2
  447. fmovex fp4,fp3
  448. |--we want P+p = W+w  but  |p| <= half ulp of P
  449. |--Then, we need to compute  A := R-P   and  a := r-p
  450. faddx fp5,fp3 |...FP3 is P
  451. fsubx fp3,fp4 |...w-P
  453. fsubx fp3,fp0 |...FP0 is A := R - P
  454.         faddx fp5,fp4 |...FP4 is p = (W-P)+w
  456. fmovex fp0,fp3 |...FP3 A
  457. fsubx fp4,fp1 |...FP1 is a := r - p
  459. |--Now we need to normalize (A,a) to  "new (R,r)" where R+r = A+a but
  460. |--|r| <= half ulp of R.
  461. faddx fp1,fp0 |...FP0 is R := A+a
  462. |--No need to calculate r if this is the last loop
  463. cmpil #0,d0
  464. jgt  RESTORE
  466. |--Need to calculate r
  467. fsubx fp0,fp3 |...A-R
  468. faddx fp3,fp1 |...FP1 is r := (A-R)+a
  469. jra  LOOP
  471. RESTORE:
  472.         fmovel fp2,a6@(N)
  473. movel A7@+,d2
  474. fmovemx A7@+,fp2-fp5
  477. movel a6@(N),d0
  478.         rorl #1,d0
  481. jra  TANCONT
  483. | end