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

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

VxWorks

开发平台：

C/C++

l_ssin.s：源码内容
  1. /* l_ssin.s - Motorola 68040 FP sine/cosine 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. ssinsa 3.2 12/18/90
  26. The entry point sSIN computes the sine of an input argument
  27. sCOS computes the cosine, and sSINCOS computes both. The
  28. corresponding entry points with a "d" computes the same
  29. corresponding function values for denormalized inputs.
  31. Input: Double-extended number X in location pointed to
  32. by address register a0.
  34. Output: The funtion value sin(X) or cos(X) returned in Fp0 if SIN or
  35. COS is requested. Otherwise, for SINCOS, sin(X) is returned
  36. in Fp0, and cos(X) is returned in Fp1.
  38. Modifies: Fp0 for SIN or COS|  both Fp0 and Fp1 for SINCOS.
  40. Accuracy and Monotonicity: The returned result is within 1 ulp in
  41. 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
  42. result is subsequently rounded to double precision. The
  43. result is provably monotonic in double precision.
  45. Speed: The programs sSIN and sCOS take approximately 150 cycles for
  46. input argument X such that |X| < 15Pi, which is the the usual
  47. situation. The speed for sSINCOS is approximately 190 cycles.
  49. Algorithm:
  51. SIN and COS:
  52. 1. If SIN is invoked, set AdjN := 0|  otherwise, set AdjN := 1.
  54. 2. If |X| >= 15Pi or |X| < 2**(-40), go to 7.
  56. 3. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
  57. k = N mod 4, so in particular, k = 0,1,2,or 3. Overwirte
  58. k by k := k + AdjN.
  60. 4. If k is even, go to 6.
  62. 5. (k is odd) Set j := (k-1)/2, sgn := (-1)**j. Return sgn*cos(r)
  63. where cos(r) is approximated by an even polynomial in r,
  64. 1 + r*r*(B1+s*(B2+ |... + s*B8)), s = r*r.
  65. Exit.
  67. 6. (k is even) Set j := k/2, sgn := (-1)**j. Return sgn*sin(r)
  68. where sin(r) is approximated by an odd polynomial in r
  69. r + r*s*(A1+s*(A2+ |... + s*A7)), s = r*r.
  70. Exit.
  72. 7. If |X| > 1, go to 9.
  74. 8. (|X|<2**(-40)) If SIN is invoked, return X|  otherwise return 1.
  76. 9. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 3.
  78. SINCOS:
  79. 1. If |X| >= 15Pi or |X| < 2**(-40), go to 6.
  81. 2. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
  82. k = N mod 4, so in particular, k = 0,1,2,or 3.
  84. 3. If k is even, go to 5.
  86. 4. (k is odd) Set j1 := (k-1)/2, j2 := j1 (EOR) (k mod 2), i.e.
  87. j1 exclusive or with the ls.b. of k.
  88. sgn1 := (-1)**j1, sgn2 := (-1)**j2.
  89. SIN(X) = sgn1 * cos(r) and COS(X) = sgn2*sin(r) where
  90. sin(r) and cos(r) are computed as odd and even polynomials
  91. in r, respectively. Exit
  93. 5. (k is even) Set j1 := k/2, sgn1 := (-1)**j1.
  94. SIN(X) = sgn1 * sin(r) and COS(X) = sgn1*cos(r) where
  95. sin(r) and cos(r) are computed as odd and even polynomials
  96. in r, respectively. Exit
  98. 6. If |X| > 1, go to 8.
  100. 7. (|X|<2**(-40)) SIN(X) = X and COS(X) = 1. Exit.
  102. 8. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 2.
  105. Copyright (C) Motorola, Inc. 1990
  106. All Rights Reserved
  108. THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF MOTOROLA
  109. The copyright notice above does not evidence any
  110. actual or intended publication of such source code.
  112. SSIN idnt 2,1 Motorola 040 Floating Point Software Package
  114. section 8
  116. NOMANUAL
  117. */
  119. #include "fpsp040L.h"
  121. BOUNDS1: .long 0x3FD78000,0x4004BC7E
  122. TWOBYPI: .long 0x3FE45F30,0x6DC9C883
  124. SINA7: .long 0xBD6AAA77,0xCCC994F5
  125. SINA6: .long 0x3DE61209,0x7AAE8DA1
  127. SINA5: .long 0xBE5AE645,0x2A118AE4
  128. SINA4: .long 0x3EC71DE3,0xA5341531
  130. SINA3: .long 0xBF2A01A0,0x1A018B59,0x00000000,0x00000000
  132. SINA2: .long 0x3FF80000,0x88888888,0x888859AF,0x00000000
  134. SINA1: .long 0xBFFC0000,0xAAAAAAAA,0xAAAAAA99,0x00000000
  136. COSB8: .long 0x3D2AC4D0,0xD6011EE3
  137. COSB7: .long 0xBDA9396F,0x9F45AC19
  139. COSB6: .long 0x3E21EED9,0x0612C972
  140. COSB5: .long 0xBE927E4F,0xB79D9FCF
  142. COSB4: .long 0x3EFA01A0,0x1A01D423,0x00000000,0x00000000
  144. COSB3: .long 0xBFF50000,0xB60B60B6,0x0B61D438,0x00000000
  146. COSB2: .long 0x3FFA0000,0xAAAAAAAA,0xAAAAAB5E
  147. COSB1: .long 0xBF000000
  149. INVTWOPI: .long 0x3FFC0000,0xA2F9836E,0x4E44152A
  151. TWOPI1: .long 0x40010000,0xC90FDAA2,0x00000000,0x00000000
  152. TWOPI2: .long 0x3FDF0000,0x85A308D4,0x00000000,0x00000000
  154. | xref __l_PITBL
  156. #define INARG FP_SCR4
  158. #define X FP_SCR5
  159. #define XDCARE X+2
  160. #define XFRAC X+4
  162. #define RPRIME FP_SCR1
  163. #define SPRIME FP_SCR2
  165. #define POSNEG1 L_SCR1
  166. #define TWOTO63 L_SCR1
  168. #define ENDFLAG L_SCR2
  169. #define N L_SCR2
  171. #define ADJN L_SCR3
  173. | xref __l_t_frcinx
  174. | xref __l_t_extdnrm
  175. | xref __l_sto_cos
  178. .text
  180. .globl __l_ssind
  181. __l_ssind:
  182. |--SIN(X) = X FOR DENORMALIZED X
  183. jra  __l_t_extdnrm
  185. .globl __l_scosd
  186. __l_scosd:
  187. |--COS(X) = 1 FOR DENORMALIZED X
  189. .long 0xf23c4400,0x3f800000 /* fmoves &0x3F800000,fp0 */
  190. |
  191. | 9D25B Fix: Sometimes the previous fmoves sets fpsr bits
  192. |
  193. fmovel #0,fpsr
  194. |
  195. jra  __l_t_frcinx
  197. .globl __l_ssin
  198. __l_ssin:
  199. |--SET ADJN TO 0
  200. movel #0,a6@(ADJN)
  201. jra  SINBGN
  203. .globl __l_scos
  204. __l_scos:
  205. |--SET ADJN TO 1
  206. movel #1,a6@(ADJN)
  208. SINBGN:
  209. |--SAVE fpcr, FP1. CHECK IF |X| IS TOO SMALL OR LARGE
  211. fmovex a0@,fp0 |...lOAD INPUT
  213. movel A0@,d0
  214. movew A0@(4),d0
  215. fmovex fp0,a6@(X)
  216. andil #0x7FFFFFFF,d0 |...COMPACTIFY X
  218. cmpil #0x3FD78000,d0 |...|X| >= 2**(-40)?
  219. jge  SOK1
  220. jra  SINSM
  222. SOK1:
  223. cmpil #0x4004BC7E,d0 |...|X| < 15 PI?
  224. jlt  SINMAIN
  225. jra  REDUCEX
  227. SINMAIN:
  228. |--THIS IS THE USUAL CASE, |X| <= 15 PI.
  229. |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
  230. fmovex fp0,fp1
  231. fmuld TWOBYPI,fp1 |...X*2/PI
  233. |--HIDE THE NEXT THREE INSTRUCTIONS
  234. lea __l_PITBL+0x200,a1 |...TABLE OF N*PI/2, N = -32,...,32
  237. |--FP1 IS NOW READY
  238. fmovel fp1,a6@(N) |...CONVERT TO INTEGER
  240. movel a6@(N),d0
  241. asll #4,d0
  242. addal d0,a1 |...A1 IS THE ADDRESS OF N*PIBY2
  243. | |...wHICH IS IN TWO PIECES Y1 # Y2
  245. fsubx A1@+,fp0 |...X-Y1
  246. |--HIDE THE NEXT ONE
  247. fsubs A1@,fp0 |...FP0 IS R = (X-Y1)-Y2
  249. SINCONT:
  250. |--continuation from REDUCEX
  252. |--GET N+ADJN AND SEE IF SIN(R) OR COS(R) IS NEEDED
  253. movel a6@(N),d0
  254. addl a6@(ADJN),d0 |...SEE IF d0 IS ODD OR EVEN
  255. rorl #1,d0 |...D0 WAS ODD IFF d0 IS NEGATIVE
  256. cmpil #0,d0
  257. jlt  COSPOLY
  259. SINPOLY:
  260. |--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
  261. |--THEN WE RETURN SGN*SIN(R). SGN*SIN(R) IS COMPUTED BY
  262. /* |--R' + R'*S*(A1 + S(A2 + S(A3 + S(A4 + |... + SA7)))), WHERE */
  263. /* |--R' = SGN*R, S=R*R. THIS CAN BE REWRITTEN AS */
  264. /* |--R' + R'*S*( [A1+T(A3+T(A5+TA7))] + [S(A2+T(A4+TA6))]) */
  265. |--WHERE T=S*S.
  266. |--NOTE THAT A3 THROUGH A7 ARE STORED IN DOUBLE PRECISION
  267. |--WHILE A1 AND A2 ARE IN DOUBLE-EXTENDED FORMAT.
  268. fmovex fp0,a6@(X) |...X IS R
  269. fmulx fp0,fp0 |...FP0 IS S
  270. |---HIDE THE NEXT TWO WHILE WAITING FOR FP0
  271. fmoved SINA7,fp3
  272. fmoved SINA6,fp2
  273. |--FP0 IS NOW READY
  274. fmovex fp0,fp1
  275. fmulx fp1,fp1 |...FP1 IS T
  276. |--HIDE THE NEXT TWO WHILE WAITING FOR FP1
  278. rorl #1,d0
  279. andil #0x80000000,d0
  280. | |...lEAST SIG. BIT OF D0 IN SIGN POSITION
  281. eorl d0,a6@(X) /* |...X IS NOW R'= SGN*R */
  283. fmulx fp1,fp3 |...TA7
  284. fmulx fp1,fp2 |...TA6
  286. faddd SINA5,fp3 |...A5+TA7
  287. faddd SINA4,fp2 |...A4+TA6
  289. fmulx fp1,fp3 |...T(A5+TA7)
  290. fmulx fp1,fp2 |...T(A4+TA6)
  292. faddd SINA3,fp3 |...A3+T(A5+TA7)
  293. faddx SINA2,fp2 |...A2+T(A4+TA6)
  295. fmulx fp3,fp1 |...T(A3+T(A5+TA7))
  297. fmulx fp0,fp2 |...S(A2+T(A4+TA6))
  298. faddx SINA1,fp1 |...A1+T(A3+T(A5+TA7))
  299. fmulx a6@(X),fp0 /* |...R'*S */
  301. faddx fp2,fp1 |...[A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))]
  302. |--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
  303. |--FP2 RELEASED, RESTORE NOW AND TAKE FULL ADVANTAGE OF HIDING
  306. fmulx fp1,fp0 /* |...SIN(R')-R' */
  307. |--FP1 RELEASED.
  309. fmovel d1,fpcr | restore users exceptions
  310. faddx a6@(X),fp0 | last inst - possible exception set
  311. jra  __l_t_frcinx
  314. COSPOLY:
  315. |--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
  316. |--THEN WE RETURN SGN*COS(R). SGN*COS(R) IS COMPUTED BY
  317. /* |--SGN + S'*(B1 + S(B2 + S(B3 + S(B4 + |... + SB8)))), WHERE */
  318. /* |--S=R*R AND S'=SGN*S. THIS CAN BE REWRITTEN AS */
  319. /* |--SGN + S'*([B1+T(B3+T(B5+TB7))] + [S(B2+T(B4+T(B6+TB8)))]) */
  320. |--WHERE T=S*S.
  321. |--NOTE THAT B4 THROUGH B8 ARE STORED IN DOUBLE PRECISION
  322. |--WHILE B2 AND B3 ARE IN DOUBLE-EXTENDED FORMAT, B1 IS -1/2
  323. |--AND IS THEREFORE STORED AS SINGLE PRECISION.
  325. fmulx fp0,fp0 |...FP0 IS S
  326. |---HIDE THE NEXT TWO WHILE WAITING FOR FP0
  327. fmoved COSB8,fp2
  328. fmoved COSB7,fp3
  329. |--FP0 IS NOW READY
  330. fmovex fp0,fp1
  331. fmulx fp1,fp1 |...FP1 IS T
  332. |--HIDE THE NEXT TWO WHILE WAITING FOR FP1
  333. fmovex fp0,a6@(X) |...X IS S
  334. rorl #1,d0
  335. andil #0x80000000,d0
  336. | |...lEAST SIG. BIT OF D0 IN SIGN POSITION
  338. fmulx fp1,fp2 |...TB8
  339. |--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
  340. eorl d0,a6@(X) /* |...X IS NOW S'= SGN*S */
  341. andil #0x80000000,d0
  343. fmulx fp1,fp3 |...TB7
  344. |--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
  345. oril #0x3F800000,d0 |...D0 IS SGN IN SINGLE
  346. movel d0,a6@(POSNEG1)
  348. faddd COSB6,fp2 |...B6+TB8
  349. faddd COSB5,fp3 |...B5+TB7
  351. fmulx fp1,fp2 |...T(B6+TB8)
  352. fmulx fp1,fp3 |...T(B5+TB7)
  354. faddd COSB4,fp2 |...B4+T(B6+TB8)
  355. faddx COSB3,fp3 |...B3+T(B5+TB7)
  357. fmulx fp1,fp2 |...T(B4+T(B6+TB8))
  358. fmulx fp3,fp1 |...T(B3+T(B5+TB7))
  360. faddx COSB2,fp2 |...B2+T(B4+T(B6+TB8))
  361. fadds COSB1,fp1 |...B1+T(B3+T(B5+TB7))
  363. fmulx fp2,fp0 |...S(B2+T(B4+T(B6+TB8)))
  364. |--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
  365. |--FP2 RELEASED.
  368. faddx fp1,fp0
  369. |--FP1 RELEASED
  371. fmulx a6@(X),fp0
  373. fmovel d1,fpcr | restore users exceptions
  374. fadds a6@(POSNEG1),fp0 | last inst - possible exception set
  375. jra  __l_t_frcinx
  378. SINBORS:
  379. |--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION.
  380. |--IF |X| < 2**(-40), RETURN X OR 1.
  381. cmpil #0x3FFF8000,d0
  382. jgt  REDUCEX
  385. SINSM:
  386. movel a6@(ADJN),d0
  387. cmpil #0,d0
  388. jgt  COSTINY
  390. SINTINY:
  391. movew #0x0000,a6@(XDCARE) |...JUST IN CASE
  392. fmovel d1,fpcr | restore users exceptions
  393. fmovex a6@(X),fp0 | last inst - possible exception set
  394. jra  __l_t_frcinx
  397. COSTINY:
  398. .long 0xf23c4400,0x3f800000  /* fmoves &0x3F800000,fp0 */
  401. fmovel d1,fpcr | restore users exceptions
  402. .long 0xf23c4428,0x00800000  /* fsubs &0x00800000,fp0 */
  403. jra  __l_t_frcinx
  406. REDUCEX:
  407. |--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW.
  408. |--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING
  409. |--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE.
  411. fmovemx fp2-fp5,A7@- |...save fp2 through fp5
  412. movel d2,A7@-
  413. .long 0xf23c4480,0x00000000 /* fmoves &0x00000000,fp1 */
  415. |--If compact form of abs(arg) in d0=0x7ffeffff, argument is so large that
  416. |--there is a danger of unwanted overflow in first LOOP iteration.  In this
  417. |--case, reduce argument by one remainder step to make subsequent reduction
  418. |--safe.
  419. cmpil #0x7ffeffff,d0 | is argument dangerously large?
  420. jne  LOOP
  421. movel #0x7ffe0000,a6@(FP_SCR2) | yes
  422. | | create 2**16383*PI/2
  423. movel #0xc90fdaa2,a6@(FP_SCR2+4)
  424. clrl a6@(FP_SCR2+8)
  425. ftstx fp0 | test sign of argument
  426. movel #0x7fdc0000,a6@(FP_SCR3) | create low half of 2**16383*
  427. | | PI/2 at FP_SCR3
  428. movel #0x85a308d3,a6@(FP_SCR3+4)
  429. clrl   a6@(FP_SCR3+8)
  430. fblt red_neg
  431. orw #0x8000,a6@(FP_SCR2) | positive arg
  432. orw #0x8000,a6@(FP_SCR3)
  433. red_neg:
  434. faddx  a6@(FP_SCR2),fp0 | high part of reduction is exact
  435. fmovex  fp0,fp1 | save high result in fp1
  436. faddx  a6@(FP_SCR3),fp0 | low part of reduction
  437. fsubx   fp0,fp1 | determine low component of result
  438. faddx  a6@(FP_SCR3),fp1 | fp0/fp1 are reduced argument.
  440. |--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4.
  441. |--integer quotient will be stored in N
  442. |--Intermeditate remainder is 66-bit long|  (R,r) in (FP0,FP1)
  444. LOOP:
  445. fmovex fp0,a6@(INARG) |...+-2**K * F, 1 <= F < 2
  446. movew a6@(INARG),d0
  447.         movel          d0,a1 |...save a copy of d0
  448. andil #0x00007FFF,d0
  449. subil #0x00003FFF,d0 |...D0 IS K
  450. cmpil #28,d0
  451. jle  LASTLOOP
  452. CONTLOOP:
  453. subil #27,d0  |...D0 IS L := K-27
  454. movel #0,a6@(ENDFLAG)
  455. jra  WORK
  456. LASTLOOP:
  457. clrl d0 |...D0 IS L := 0
  458. movel #1,a6@(ENDFLAG)
  460. WORK:
  461. |--FIND THE REMAINDER OF (R,r) W.R.T. 2**L * (PI/2). L IS SO CHOSEN
  462. |--THAT INT( X * (2/PI) / 2**(L) ) < 2**29.
  464. |--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63),
  465. |--2**L * (PIby2_1), 2**L * (PIby2_2)
  467. movel #0x00003FFE,d2 |...BIASED EXPO OF 2/PI
  468. subl d0,d2 |...BIASED EXPO OF 2**(-L)*(2/PI)
  470. movel #0xA2F9836E,a6@(FP_SCR1+4)
  471. movel #0x4E44152A,a6@(FP_SCR1+8)
  472. movew d2,a6@(FP_SCR1) |...FP_SCR1 is 2**(-L)*(2/PI)
  474. fmovex fp0,fp2
  475. fmulx a6@(FP_SCR1),fp2
  476. |--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN
  477. /* |--FLOATING POINT FORMAT, THE TWO FMOVE'S FMOVE.l FP <--> N */
  478. |--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT
  479. |--(SIGN(INARG)*2**63 + FP2) - SIGN(INARG)*2**63 WILL GIVE
  480. |--US THE DESIRED VALUE IN FLOATING POINT.
  482. |--HIDE SIX CYCLES OF INSTRUCTION
  483.         movel a1,d2
  484.         swap d2
  485. andil #0x80000000,d2
  486. oril #0x5F000000,d2 |...D2 IS SIGN(INARG)*2**63 IN SGL
  487. movel d2,a6@(TWOTO63)
  489. movel d0,d2
  490. addil #0x00003FFF,d2 |...BIASED EXPO OF 2**L * (PI/2)
  492. |--FP2 IS READY
  493. fadds a6@(TWOTO63),fp2 |...THE FRACTIONAL PART OF fp1 IS ROUNDED
  495. |--HIDE 4 CYCLES OF INSTRUCTION|  creating 2**(L)*Piby2_1  and  2**(L)*Piby2_2
  496.         movew d2,a6@(FP_SCR2)
  497. clrw           a6@(FP_SCR2+2)
  498. movel #0xC90FDAA2,a6@(FP_SCR2+4)
  499. clrl a6@(FP_SCR2+8) |...FP_SCR2 is  2**(L) * Piby2_1
  501. |--FP2 IS READY
  502. fsubs a6@(TWOTO63),fp2 |...FP2 is N
  504. addil #0x00003FDD,d0
  505.         movew d0,a6@(FP_SCR3)
  506. clrw           a6@(FP_SCR3+2)
  507. movel #0x85A308D3,a6@(FP_SCR3+4)
  508. clrl a6@(FP_SCR3+8) |...FP_SCR3 is 2**(L) * Piby2_2
  510. movel a6@(ENDFLAG),d0
  512. |--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and
  513. |--P2 = 2**(L) * Piby2_2
  514. fmovex fp2,fp4
  515. fmulx a6@(FP_SCR2),fp4 |...w = N*P1
  516. fmovex fp2,fp5
  517. fmulx a6@(FP_SCR3),fp5 |...w = N*P2
  518. fmovex fp4,fp3
  519. |--we want P+p = W+w  but  |p| <= half ulp of P
  520. |--Then, we need to compute  A := R-P   and  a := r-p
  521. faddx fp5,fp3 |...FP3 is P
  522. fsubx fp3,fp4 |...w-P
  524. fsubx fp3,fp0 |...FP0 is A := R - P
  525.         faddx fp5,fp4 |...FP4 is p = (W-P)+w
  527. fmovex fp0,fp3 |...FP3 A
  528. fsubx fp4,fp1 |...FP1 is a := r - p
  530. |--Now we need to normalize (A,a) to  "new (R,r)" where R+r = A+a but
  531. |--|r| <= half ulp of R.
  532. faddx fp1,fp0 |...FP0 is R := A+a
  533. |--No need to calculate r if this is the last loop
  534. cmpil #0,d0
  535. jgt  RESTORE
  537. |--Need to calculate r
  538. fsubx fp0,fp3 |...A-R
  539. faddx fp3,fp1 |...FP1 is r := (A-R)+a
  540. jra  LOOP
  542. RESTORE:
  543.         fmovel fp2,a6@(N)
  544. movel A7@+,d2
  545. fmovemx A7@+,fp2-fp5
  548. movel a6@(ADJN),d0
  549. cmpil #4,d0
  551. jlt  SINCONT
  552. jra  SCCONT
  554. .globl __l_ssincosd
  555. __l_ssincosd:
  556. |--SIN AND COS OF X FOR DENORMALIZED X
  558. .long 0xf23c4480,0x3f800000 /* fmoves &0x3F800000,fp1 */
  559. bsrl __l_sto_cos | store cosine result
  560. jra  __l_t_extdnrm
  562. .globl __l_ssincos
  563. __l_ssincos:
  564. |--SET ADJN TO 4
  565. movel #4,a6@(ADJN)
  567. fmovex a0@,fp0 |...lOAD INPUT
  569. movel A0@,d0
  570. movew A0@(4),d0
  571. fmovex fp0,a6@(X)
  572. andil #0x7FFFFFFF,d0 |...COMPACTIFY X
  574. cmpil #0x3FD78000,d0 |...|X| >= 2**(-40)?
  575. jge  SCOK1
  576. jra  SCSM
  578. SCOK1:
  579. cmpil #0x4004BC7E,d0 |...|X| < 15 PI?
  580. jlt  SCMAIN
  581. jra  REDUCEX
  584. SCMAIN:
  585. |--THIS IS THE USUAL CASE, |X| <= 15 PI.
  586. |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
  587. fmovex fp0,fp1
  588. fmuld TWOBYPI,fp1 |...X*2/PI
  590. |--HIDE THE NEXT THREE INSTRUCTIONS
  591. lea __l_PITBL+0x200,a1 |...TABLE OF N*PI/2, N = -32,...,32
  594. |--FP1 IS NOW READY
  595. fmovel fp1,a6@(N) |...CONVERT TO INTEGER
  597. movel a6@(N),d0
  598. asll #4,d0
  599. addal d0,a1 |...ADDRESS OF N*PIBY2, IN Y1, Y2
  601. fsubx A1@+,fp0 |...X-Y1
  602.         fsubs A1@,fp0 |...FP0 IS R = (X-Y1)-Y2
  604. SCCONT:
  605. |--continuation point from REDUCEX
  607. |--HIDE THE NEXT TWO
  608. movel a6@(N),d0
  609. rorl #1,d0
  611. cmpil #0,d0 |...D0 < 0 IFF N IS ODD
  612. jge  NEVEN
  614. NODD:
  615. |--REGISTERS SAVED SO FAR: D0, A0, FP2.
  617. fmovex fp0,a6@(RPRIME)
  618. fmulx fp0,fp0   |...FP0 IS S = R*R
  619. fmoved SINA7,fp1 |...A7
  620. fmoved COSB8,fp2 |...B8
  621. fmulx fp0,fp1   |...SA7
  622. movel d2,A7@-
  623. movel d0,d2
  624. fmulx fp0,fp2   |...SB8
  625. rorl #1,d2
  626. andil #0x80000000,d2
  628. faddd SINA6,fp1 |...A6+SA7
  629. eorl d0,d2
  630. andil #0x80000000,d2
  631. faddd COSB7,fp2 |...B7+SB8
  633. fmulx fp0,fp1   |...S(A6+SA7)
  634. eorl d2,a6@(RPRIME)
  635. movel A7@+,d2
  636. fmulx fp0,fp2   |...S(B7+SB8)
  637. rorl #1,d0
  638. andil #0x80000000,d0
  640. faddd SINA5,fp1 |...A5+S(A6+SA7)
  641. movel #0x3F800000,a6@(POSNEG1)
  642. eorl d0,a6@(POSNEG1)
  643. faddd COSB6,fp2 |...B6+S(B7+SB8)
  645. fmulx fp0,fp1   |...S(A5+S(A6+SA7))
  646. fmulx fp0,fp2   |...S(B6+S(B7+SB8))
  647. fmovex fp0,a6@(SPRIME)
  649. faddd SINA4,fp1 |...A4+S(A5+S(A6+SA7))
  650. eorl d0,a6@(SPRIME)
  651. faddd COSB5,fp2 |...B5+S(B6+S(B7+SB8))
  653. fmulx fp0,fp1   |...S(A4+...)
  654. fmulx fp0,fp2   |...S(B5+...)
  656. faddd SINA3,fp1 |...A3+S(A4+...)
  657. faddd COSB4,fp2 |...B4+S(B5+...)
  659. fmulx fp0,fp1   |...S(A3+...)
  660. fmulx fp0,fp2   |...S(B4+...)
  662. faddx SINA2,fp1 |...A2+S(A3+...)
  663. faddx COSB3,fp2 |...B3+S(B4+...)
  665. fmulx fp0,fp1   |...S(A2+...)
  666. fmulx fp0,fp2   |...S(B3+...)
  668. faddx SINA1,fp1 |...A1+S(A2+...)
  669. faddx COSB2,fp2 |...B2+S(B3+...)
  671. fmulx fp0,fp1   |...S(A1+...)
  672. fmulx fp2,fp0   |...S(B2+...)
  676. fmulx a6@(RPRIME),fp1 /* |...R'S(A1+...) */
  677. fadds COSB1,fp0 |...B1+S(B2...)
  678. fmulx a6@(SPRIME),fp0 /* |...S'(B1+S(B2+...)) */
  680. movel d1,a7@- | restore users mode # precision
  681. andil #0xff,d1 | mask off all exceptions
  682. fmovel d1,fpcr
  683. faddx a6@(RPRIME),fp1 |...COS(X)
  684. bsrl __l_sto_cos | store cosine result
  685. fmovel a7@+,fpcr | restore users exceptions
  686. fadds a6@(POSNEG1),fp0 |...SIN(X)
  688. jra  __l_t_frcinx
  691. NEVEN:
  692. |--REGISTERS SAVED SO FAR: FP2.
  694. fmovex fp0,a6@(RPRIME)
  695. fmulx fp0,fp0   |...FP0 IS S = R*R
  696. fmoved COSB8,fp1 |...B8
  697. fmoved SINA7,fp2 |...A7
  698. fmulx fp0,fp1   |...SB8
  699. fmovex fp0,a6@(SPRIME)
  700. fmulx fp0,fp2   |...SA7
  701. rorl #1,d0
  702. andil #0x80000000,d0
  703. faddd COSB7,fp1 |...B7+SB8
  704. faddd SINA6,fp2 |...A6+SA7
  705. eorl d0,a6@(RPRIME)
  706. eorl d0,a6@(SPRIME)
  707. fmulx fp0,fp1   |...S(B7+SB8)
  708. oril #0x3F800000,d0
  709. movel d0,a6@(POSNEG1)
  710. fmulx fp0,fp2   |...S(A6+SA7)
  712. faddd COSB6,fp1 |...B6+S(B7+SB8)
  713. faddd SINA5,fp2 |...A5+S(A6+SA7)
  715. fmulx fp0,fp1   |...S(B6+S(B7+SB8))
  716. fmulx fp0,fp2   |...S(A5+S(A6+SA7))
  718. faddd COSB5,fp1 |...B5+S(B6+S(B7+SB8))
  719. faddd SINA4,fp2 |...A4+S(A5+S(A6+SA7))
  721. fmulx fp0,fp1   |...S(B5+...)
  722. fmulx fp0,fp2   |...S(A4+...)
  724. faddd COSB4,fp1 |...B4+S(B5+...)
  725. faddd SINA3,fp2 |...A3+S(A4+...)
  727. fmulx fp0,fp1   |...S(B4+...)
  728. fmulx fp0,fp2   |...S(A3+...)
  730. faddx COSB3,fp1 |...B3+S(B4+...)
  731. faddx SINA2,fp2 |...A2+S(A3+...)
  733. fmulx fp0,fp1   |...S(B3+...)
  734. fmulx fp0,fp2   |...S(A2+...)
  736. faddx COSB2,fp1 |...B2+S(B3+...)
  737. faddx SINA1,fp2 |...A1+S(A2+...)
  739. fmulx fp0,fp1   |...S(B2+...)
  740. fmulx fp2,fp0   |...s(a1+...)
  744. fadds COSB1,fp1 |...B1+S(B2...)
  745. fmulx a6@(RPRIME),fp0 /*  |...R'S(A1+...) */
  746. fmulx a6@(SPRIME),fp1 /*  |...S'(B1+S(B2+...)) */
  748. movel d1,a7@- | save users mode # precision
  749. andil #0xff,d1 | mask off all exceptions
  750. fmovel d1,fpcr
  751. fadds a6@(POSNEG1),fp1 |...COS(X)
  752. bsrl __l_sto_cos | store cosine result
  753. fmovel a7@+,fpcr | restore users exceptions
  754. faddx a6@(RPRIME),fp0 |...SIN(X)
  756. jra  __l_t_frcinx
  758. SCBORS:
  759. cmpil #0x3FFF8000,d0
  760. jgt  REDUCEX
  763. SCSM:
  764. movew #0x0000,a6@(XDCARE)
  765. .long 0xf23c4480,0x3f800000 /* fmoves &0x3F800000,fp1 */
  767. movel d1,a7@- | save users mode # precision
  768. andil #0xff,d1 | mask off all exceptions
  769. fmovel d1,fpcr
  770. .long 0xf23c44a8,0x00800000 /* fsubs &0x00800000,fp1 */
  771. bsrl __l_sto_cos | store cosine result
  772. fmovel a7@+,fpcr | restore users exceptions
  773. fmovex a6@(X),fp0
  774. jra  __l_t_frcinx
  776. | end