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

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

VxWorks

开发平台：

C/C++

ssin.s：源码内容
  1. /* ssin.s - Motorola 68040 FP sine routines (EXC) */
  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  added modification history; 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
  25. ssinsa 3.2 12/18/90
  27.  WIND RIVER MODIFICATION HISTORY
  28.  01a,31jul91,kdl from Motorola FPSP v2.0.
  30. The entry point sSIN computes the sine of an input argument
  31. sCOS computes the cosine, and sSINCOS computes both. The
  32. corresponding entry points with a "d" computes the same
  33. corresponding function values for denormalized inputs.
  35. Input: Double-extended number X in location pointed to
  36. by address register a0.
  38. Output: The funtion value sin(X) or cos(X) returned in Fp0 if SIN or
  39. COS is requested. Otherwise, for SINCOS, sin(X) is returned
  40. in Fp0, and cos(X) is returned in Fp1.
  42. Modifies: Fp0 for SIN or COS|  both Fp0 and Fp1 for SINCOS.
  44. Accuracy and Monotonicity: The returned result is within 1 ulp in
  45. 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
  46. result is subsequently rounded to double precision. The
  47. result is provably monotonic in double precision.
  49. Speed: The programs sSIN and sCOS take approximately 150 cycles for
  50. input argument X such that |X| < 15Pi, which is the the usual
  51. situation. The speed for sSINCOS is approximately 190 cycles.
  53. Algorithm:
  55. SIN and COS:
  56. 1. If SIN is invoked, set AdjN := 0|  otherwise, set AdjN := 1.
  58. 2. If |X| >= 15Pi or |X| < 2**(-40), go to 7.
  60. 3. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
  61. k = N mod 4, so in particular, k = 0,1,2,or 3. Overwirte
  62. k by k := k + AdjN.
  64. 4. If k is even, go to 6.
  66. 5. (k is odd) Set j := (k-1)/2, sgn := (-1)**j. Return sgn*cos(r)
  67. where cos(r) is approximated by an even polynomial in r,
  68. 1 + r*r*(B1+s*(B2+ |... + s*B8)), s = r*r.
  69. Exit.
  71. 6. (k is even) Set j := k/2, sgn := (-1)**j. Return sgn*sin(r)
  72. where sin(r) is approximated by an odd polynomial in r
  73. r + r*s*(A1+s*(A2+ |... + s*A7)), s = r*r.
  74. Exit.
  76. 7. If |X| > 1, go to 9.
  78. 8. (|X|<2**(-40)) If SIN is invoked, return X|  otherwise return 1.
  80. 9. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 3.
  82. SINCOS:
  83. 1. If |X| >= 15Pi or |X| < 2**(-40), go to 6.
  85. 2. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
  86. k = N mod 4, so in particular, k = 0,1,2,or 3.
  88. 3. If k is even, go to 5.
  90. 4. (k is odd) Set j1 := (k-1)/2, j2 := j1 (EOR) (k mod 2), i.e.
  91. j1 exclusive or with the ls.b. of k.
  92. sgn1 := (-1)**j1, sgn2 := (-1)**j2.
  93. SIN(X) = sgn1 * cos(r) and COS(X) = sgn2*sin(r) where
  94. sin(r) and cos(r) are computed as odd and even polynomials
  95. in r, respectively. Exit
  97. 5. (k is even) Set j1 := k/2, sgn1 := (-1)**j1.
  98. SIN(X) = sgn1 * sin(r) and COS(X) = sgn1*cos(r) where
  99. sin(r) and cos(r) are computed as odd and even polynomials
  100. in r, respectively. Exit
  102. 6. If |X| > 1, go to 8.
  104. 7. (|X|<2**(-40)) SIN(X) = X and COS(X) = 1. Exit.
  106. 8. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 2.
  109. Copyright (C) Motorola, Inc. 1990
  110. All Rights Reserved
  112. THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF MOTOROLA
  113. The copyright notice above does not evidence any
  114. actual or intended publication of such source code.
  116. SSIN idnt 2,1 Motorola 040 Floating Point Software Package
  118. section 8
  120. NOMANUAL
  121. */
  123. #include "fpsp040E.h"
  125. BOUNDS1: .long 0x3FD78000,0x4004BC7E
  126. TWOBYPI: .long 0x3FE45F30,0x6DC9C883
  128. SINA7: .long 0xBD6AAA77,0xCCC994F5
  129. SINA6: .long 0x3DE61209,0x7AAE8DA1
  131. SINA5: .long 0xBE5AE645,0x2A118AE4
  132. SINA4: .long 0x3EC71DE3,0xA5341531
  134. SINA3: .long 0xBF2A01A0,0x1A018B59,0x00000000,0x00000000
  136. SINA2: .long 0x3FF80000,0x88888888,0x888859AF,0x00000000
  138. SINA1: .long 0xBFFC0000,0xAAAAAAAA,0xAAAAAA99,0x00000000
  140. COSB8: .long 0x3D2AC4D0,0xD6011EE3
  141. COSB7: .long 0xBDA9396F,0x9F45AC19
  143. COSB6: .long 0x3E21EED9,0x0612C972
  144. COSB5: .long 0xBE927E4F,0xB79D9FCF
  146. COSB4: .long 0x3EFA01A0,0x1A01D423,0x00000000,0x00000000
  148. COSB3: .long 0xBFF50000,0xB60B60B6,0x0B61D438,0x00000000
  150. COSB2: .long 0x3FFA0000,0xAAAAAAAA,0xAAAAAB5E
  151. COSB1: .long 0xBF000000
  153. INVTWOPI: .long 0x3FFC0000,0xA2F9836E,0x4E44152A
  155. TWOPI1: .long 0x40010000,0xC90FDAA2,0x00000000,0x00000000
  156. TWOPI2: .long 0x3FDF0000,0x85A308D4,0x00000000,0x00000000
  158. | xref __x_PITBL
  160. #define INARG FP_SCR4
  162. #define X FP_SCR5
  163. #define XDCARE X+2
  164. #define XFRAC X+4
  166. #define RPRIME FP_SCR1
  167. #define SPRIME FP_SCR2
  169. #define POSNEG1 L_SCR1
  170. #define TWOTO63 L_SCR1
  172. #define ENDFLAG L_SCR2
  173. #define N L_SCR2
  175. #define ADJN L_SCR3
  177. | xref __x_t_frcinx
  178. | xref __x_t_extdnrm
  179. | xref __x_sto_cos
  182. .text
  184. .globl __x_ssind
  185. __x_ssind:
  186. |--SIN(X) = X FOR DENORMALIZED X
  187. jra  __x_t_extdnrm
  189. .globl __x_scosd
  190. __x_scosd:
  191. |--COS(X) = 1 FOR DENORMALIZED X
  193. /* fmoves &0x3F800000,fp0 */  .long 0xf23c4400,0x3f800000
  194. |
  195. | 9D25B Fix: Sometimes the previous fmoves sets fpsr bits
  196. |
  197. fmovel #0,fpsr
  198. |
  199. jra  __x_t_frcinx
  201. .globl __x_ssin
  202. __x_ssin:
  203. |--SET ADJN TO 0
  204. movel #0,a6@(ADJN)
  205. jra  SINBGN
  207. .globl __x_scos
  208. __x_scos:
  209. |--SET ADJN TO 1
  210. movel #1,a6@(ADJN)
  212. SINBGN:
  213. |--SAVE fpcr, FP1. CHECK IF |X| IS TOO SMALL OR LARGE
  215. fmovex a0@,fp0 |...lOAD INPUT
  217. movel A0@,d0
  218. movew A0@(4),d0
  219. fmovex fp0,a6@(X)
  220. andil #0x7FFFFFFF,d0 |...COMPACTIFY X
  222. cmpil #0x3FD78000,d0 |...|X| >= 2**(-40)?
  223. jge  SOK1
  224. jra  SINSM
  226. SOK1:
  227. cmpil #0x4004BC7E,d0 |...|X| < 15 PI?
  228. jlt  SINMAIN
  229. jra  REDUCEX
  231. SINMAIN:
  232. |--THIS IS THE USUAL CASE, |X| <= 15 PI.
  233. |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
  234. fmovex fp0,fp1
  235. fmuld TWOBYPI,fp1 |...X*2/PI
  237. |--HIDE THE NEXT THREE INSTRUCTIONS
  238. lea __x_PITBL+0x200,a1 |...TABLE OF N*PI/2, N = -32,...,32
  241. |--FP1 IS NOW READY
  242. fmovel fp1,a6@(N) |...CONVERT TO INTEGER
  244. movel a6@(N),d0
  245. asll #4,d0
  246. addal d0,a1 |...A1 IS THE ADDRESS OF N*PIBY2
  247. | |...wHICH IS IN TWO PIECES Y1 # Y2
  249. fsubx A1@+,fp0 |...X-Y1
  250. |--HIDE THE NEXT ONE
  251. fsubs A1@,fp0 |...FP0 IS R = (X-Y1)-Y2
  253. SINCONT:
  254. |--continuation from REDUCEX
  256. |--GET N+ADJN AND SEE IF SIN(R) OR COS(R) IS NEEDED
  257. movel a6@(N),d0
  258. addl a6@(ADJN),d0 |...SEE IF d0 IS ODD OR EVEN
  259. rorl #1,d0 |...D0 WAS ODD IFF d0 IS NEGATIVE
  260. cmpil #0,d0
  261. jlt  COSPOLY
  263. SINPOLY:
  264. |--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
  265. |--THEN WE RETURN SGN*SIN(R). SGN*SIN(R) IS COMPUTED BY
  266. /* |--R' + R'*S*(A1 + S(A2 + S(A3 + S(A4 + |... + SA7)))), WHERE */
  267. /* |--R' = SGN*R, S=R*R. THIS CAN BE REWRITTEN AS */
  268. /* |--R' + R'*S*( [A1+T(A3+T(A5+TA7))] + [S(A2+T(A4+TA6))]) */
  269. |--WHERE T=S*S.
  270. |--NOTE THAT A3 THROUGH A7 ARE STORED IN DOUBLE PRECISION
  271. |--WHILE A1 AND A2 ARE IN DOUBLE-EXTENDED FORMAT.
  272. fmovex fp0,a6@(X) |...X IS R
  273. fmulx fp0,fp0 |...FP0 IS S
  274. |---HIDE THE NEXT TWO WHILE WAITING FOR FP0
  275. fmoved SINA7,fp3
  276. fmoved SINA6,fp2
  277. |--FP0 IS NOW READY
  278. fmovex fp0,fp1
  279. fmulx fp1,fp1 |...FP1 IS T
  280. |--HIDE THE NEXT TWO WHILE WAITING FOR FP1
  282. rorl #1,d0
  283. andil #0x80000000,d0
  284. | |...lEAST SIG. BIT OF D0 IN SIGN POSITION
  285. eorl d0,a6@(X) /* |...X IS NOW R'= SGN*R */
  287. fmulx fp1,fp3 |...TA7
  288. fmulx fp1,fp2 |...TA6
  290. faddd SINA5,fp3 |...A5+TA7
  291. faddd SINA4,fp2 |...A4+TA6
  293. fmulx fp1,fp3 |...T(A5+TA7)
  294. fmulx fp1,fp2 |...T(A4+TA6)
  296. faddd SINA3,fp3 |...A3+T(A5+TA7)
  297. faddx SINA2,fp2 |...A2+T(A4+TA6)
  299. fmulx fp3,fp1 |...T(A3+T(A5+TA7))
  301. fmulx fp0,fp2 |...S(A2+T(A4+TA6))
  302. faddx SINA1,fp1 |...A1+T(A3+T(A5+TA7))
  303. fmulx a6@(X),fp0 /* |...R'*S */
  305. faddx fp2,fp1 |...[A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))]
  306. |--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
  307. |--FP2 RELEASED, RESTORE NOW AND TAKE FULL ADVANTAGE OF HIDING
  310. fmulx fp1,fp0 /* |...SIN(R')-R' */
  311. |--FP1 RELEASED.
  313. fmovel d1,fpcr | restore users exceptions
  314. faddx a6@(X),fp0 | last inst - possible exception set
  315. jra  __x_t_frcinx
  318. COSPOLY:
  319. |--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
  320. |--THEN WE RETURN SGN*COS(R). SGN*COS(R) IS COMPUTED BY
  321. /* |--SGN + S'*(B1 + S(B2 + S(B3 + S(B4 + |... + SB8)))), WHERE */
  322. /* |--S=R*R AND S'=SGN*S. THIS CAN BE REWRITTEN AS */
  323. /* |--SGN + S'*([B1+T(B3+T(B5+TB7))] + [S(B2+T(B4+T(B6+TB8)))]) */
  324. |--WHERE T=S*S.
  325. |--NOTE THAT B4 THROUGH B8 ARE STORED IN DOUBLE PRECISION
  326. |--WHILE B2 AND B3 ARE IN DOUBLE-EXTENDED FORMAT, B1 IS -1/2
  327. |--AND IS THEREFORE STORED AS SINGLE PRECISION.
  329. fmulx fp0,fp0 |...FP0 IS S
  330. |---HIDE THE NEXT TWO WHILE WAITING FOR FP0
  331. fmoved COSB8,fp2
  332. fmoved COSB7,fp3
  333. |--FP0 IS NOW READY
  334. fmovex fp0,fp1
  335. fmulx fp1,fp1 |...FP1 IS T
  336. |--HIDE THE NEXT TWO WHILE WAITING FOR FP1
  337. fmovex fp0,a6@(X) |...X IS S
  338. rorl #1,d0
  339. andil #0x80000000,d0
  340. | |...lEAST SIG. BIT OF D0 IN SIGN POSITION
  342. fmulx fp1,fp2 |...TB8
  343. |--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
  344. eorl d0,a6@(X) /* |...X IS NOW S'= SGN*S */
  345. andil #0x80000000,d0
  347. fmulx fp1,fp3 |...TB7
  348. |--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
  349. oril #0x3F800000,d0 |...D0 IS SGN IN SINGLE
  350. movel d0,a6@(POSNEG1)
  352. faddd COSB6,fp2 |...B6+TB8
  353. faddd COSB5,fp3 |...B5+TB7
  355. fmulx fp1,fp2 |...T(B6+TB8)
  356. fmulx fp1,fp3 |...T(B5+TB7)
  358. faddd COSB4,fp2 |...B4+T(B6+TB8)
  359. faddx COSB3,fp3 |...B3+T(B5+TB7)
  361. fmulx fp1,fp2 |...T(B4+T(B6+TB8))
  362. fmulx fp3,fp1 |...T(B3+T(B5+TB7))
  364. faddx COSB2,fp2 |...B2+T(B4+T(B6+TB8))
  365. fadds COSB1,fp1 |...B1+T(B3+T(B5+TB7))
  367. fmulx fp2,fp0 |...S(B2+T(B4+T(B6+TB8)))
  368. |--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
  369. |--FP2 RELEASED.
  372. faddx fp1,fp0
  373. |--FP1 RELEASED
  375. fmulx a6@(X),fp0
  377. fmovel d1,fpcr | restore users exceptions
  378. fadds a6@(POSNEG1),fp0 | last inst - possible exception set
  379. jra  __x_t_frcinx
  382. SINBORS:
  383. |--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION.
  384. |--IF |X| < 2**(-40), RETURN X OR 1.
  385. cmpil #0x3FFF8000,d0
  386. jgt  REDUCEX
  389. SINSM:
  390. movel a6@(ADJN),d0
  391. cmpil #0,d0
  392. jgt  COSTINY
  394. SINTINY:
  395. movew #0x0000,a6@(XDCARE) |...JUST IN CASE
  396. fmovel d1,fpcr | restore users exceptions
  397. fmovex a6@(X),fp0 | last inst - possible exception set
  398. jra  __x_t_frcinx
  401. COSTINY:
  402. .long 0xf23c4400,0x3f800000 /*  fmoves  &0x3F800000,fp0 */
  404. fmovel d1,fpcr | restore users exceptions
  405. .long 0xf23c4428,0x00800000 /*  fsubs  &0x00800000,fp0 */
  406. jra  __x_t_frcinx
  409. REDUCEX:
  410. |--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW.
  411. |--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING
  412. |--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE.
  414. fmovemx fp2-fp5,A7@- |...save fp2 through fp5
  415. movel d2,A7@-
  416. /* fmoves &0x00000000,fp1 */  .long 0xf23c4480,0x00000000
  418. |--If compact form of abs(arg) in d0=0x7ffeffff, argument is so large that
  419. |--there is a danger of unwanted overflow in first LOOP iteration.  In this
  420. |--case, reduce argument by one remainder step to make subsequent reduction
  421. |--safe.
  422. cmpil #0x7ffeffff,d0 | is argument dangerously large?
  423. jne  LOOP
  424. movel #0x7ffe0000,a6@(FP_SCR2) | yes
  425. | | create 2**16383*PI/2
  426. movel #0xc90fdaa2,a6@(FP_SCR2+4)
  427. clrl a6@(FP_SCR2+8)
  428. ftstx fp0 | test sign of argument
  429. movel #0x7fdc0000,a6@(FP_SCR3) | create low half of 2**16383*
  430. | | PI/2 at FP_SCR3
  431. movel #0x85a308d3,a6@(FP_SCR3+4)
  432. clrl   a6@(FP_SCR3+8)
  433. fblt red_neg
  434. orw #0x8000,a6@(FP_SCR2) | positive arg
  435. orw #0x8000,a6@(FP_SCR3)
  436. red_neg:
  437. faddx  a6@(FP_SCR2),fp0 | high part of reduction is exact
  438. fmovex  fp0,fp1 | save high result in fp1
  439. faddx  a6@(FP_SCR3),fp0 | low part of reduction
  440. fsubx   fp0,fp1 | determine low component of result
  441. faddx  a6@(FP_SCR3),fp1 | fp0/fp1 are reduced argument.
  443. |--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4.
  444. |--integer quotient will be stored in N
  445. |--Intermeditate remainder is 66-bit long|  (R,r) in (FP0,FP1)
  447. LOOP:
  448. fmovex fp0,a6@(INARG) |...+-2**K * F, 1 <= F < 2
  449. movew a6@(INARG),d0
  450.         movel          d0,a1 |...save a copy of d0
  451. andil #0x00007FFF,d0
  452. subil #0x00003FFF,d0 |...D0 IS K
  453. cmpil #28,d0
  454. jle  LASTLOOP
  455. CONTLOOP:
  456. subil #27,d0  |...D0 IS L := K-27
  457. movel #0,a6@(ENDFLAG)
  458. jra  WORK
  459. LASTLOOP:
  460. clrl d0 |...D0 IS L := 0
  461. movel #1,a6@(ENDFLAG)
  463. WORK:
  464. |--FIND THE REMAINDER OF (R,r) W.R.T. 2**L * (PI/2). L IS SO CHOSEN
  465. |--THAT INT( X * (2/PI) / 2**(L) ) < 2**29.
  467. |--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63),
  468. |--2**L * (PIby2_1), 2**L * (PIby2_2)
  470. movel #0x00003FFE,d2 |...BIASED EXPO OF 2/PI
  471. subl d0,d2 |...BIASED EXPO OF 2**(-L)*(2/PI)
  473. movel #0xA2F9836E,a6@(FP_SCR1+4)
  474. movel #0x4E44152A,a6@(FP_SCR1+8)
  475. movew d2,a6@(FP_SCR1) |...FP_SCR1 is 2**(-L)*(2/PI)
  477. fmovex fp0,fp2
  478. fmulx a6@(FP_SCR1),fp2
  479. |--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN
  480. /* |--FLOATING POINT FORMAT, THE TWO FMOVE'S FMOVE.l FP <--> N */
  481. |--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT
  482. |--(SIGN(INARG)*2**63 + FP2) - SIGN(INARG)*2**63 WILL GIVE
  483. |--US THE DESIRED VALUE IN FLOATING POINT.
  485. |--HIDE SIX CYCLES OF INSTRUCTION
  486.         movel a1,d2
  487.         swap d2
  488. andil #0x80000000,d2
  489. oril #0x5F000000,d2 |...D2 IS SIGN(INARG)*2**63 IN SGL
  490. movel d2,a6@(TWOTO63)
  492. movel d0,d2
  493. addil #0x00003FFF,d2 |...BIASED EXPO OF 2**L * (PI/2)
  495. |--FP2 IS READY
  496. fadds a6@(TWOTO63),fp2 |...THE FRACTIONAL PART OF fp1 IS ROUNDED
  498. |--HIDE 4 CYCLES OF INSTRUCTION|  creating 2**(L)*Piby2_1  and  2**(L)*Piby2_2
  499.         movew d2,a6@(FP_SCR2)
  500. clrw           a6@(FP_SCR2+2)
  501. movel #0xC90FDAA2,a6@(FP_SCR2+4)
  502. clrl a6@(FP_SCR2+8) |...FP_SCR2 is  2**(L) * Piby2_1
  504. |--FP2 IS READY
  505. fsubs a6@(TWOTO63),fp2 |...FP2 is N
  507. addil #0x00003FDD,d0
  508.         movew d0,a6@(FP_SCR3)
  509. clrw           a6@(FP_SCR3+2)
  510. movel #0x85A308D3,a6@(FP_SCR3+4)
  511. clrl a6@(FP_SCR3+8) |...FP_SCR3 is 2**(L) * Piby2_2
  513. movel a6@(ENDFLAG),d0
  515. |--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and
  516. |--P2 = 2**(L) * Piby2_2
  517. fmovex fp2,fp4
  518. fmulx a6@(FP_SCR2),fp4 |...w = N*P1
  519. fmovex fp2,fp5
  520. fmulx a6@(FP_SCR3),fp5 |...w = N*P2
  521. fmovex fp4,fp3
  522. |--we want P+p = W+w  but  |p| <= half ulp of P
  523. |--Then, we need to compute  A := R-P   and  a := r-p
  524. faddx fp5,fp3 |...FP3 is P
  525. fsubx fp3,fp4 |...w-P
  527. fsubx fp3,fp0 |...FP0 is A := R - P
  528.         faddx fp5,fp4 |...FP4 is p = (W-P)+w
  530. fmovex fp0,fp3 |...FP3 A
  531. fsubx fp4,fp1 |...FP1 is a := r - p
  533. |--Now we need to normalize (A,a) to  "new (R,r)" where R+r = A+a but
  534. |--|r| <= half ulp of R.
  535. faddx fp1,fp0 |...FP0 is R := A+a
  536. |--No need to calculate r if this is the last loop
  537. cmpil #0,d0
  538. jgt  RESTORE
  540. |--Need to calculate r
  541. fsubx fp0,fp3 |...A-R
  542. faddx fp3,fp1 |...FP1 is r := (A-R)+a
  543. jra  LOOP
  545. RESTORE:
  546.         fmovel fp2,a6@(N)
  547. movel A7@+,d2
  548. fmovemx A7@+,fp2-fp5
  551. movel a6@(ADJN),d0
  552. cmpil #4,d0
  554. jlt  SINCONT
  555. jra  SCCONT
  557. .globl __x_ssincosd
  558. __x_ssincosd:
  559. |--SIN AND COS OF X FOR DENORMALIZED X
  561. .long 0xf23c4480,0x3f800000 /*  fmoves  &0x3F800000,fp1 */
  562. bsrl __x_sto_cos | store cosine result
  563. jra  __x_t_extdnrm
  565. .globl __x_ssincos
  566. __x_ssincos:
  567. |--SET ADJN TO 4
  568. movel #4,a6@(ADJN)
  570. fmovex a0@,fp0 |...lOAD INPUT
  572. movel A0@,d0
  573. movew A0@(4),d0
  574. fmovex fp0,a6@(X)
  575. andil #0x7FFFFFFF,d0 |...COMPACTIFY X
  577. cmpil #0x3FD78000,d0 |...|X| >= 2**(-40)?
  578. jge  SCOK1
  579. jra  SCSM
  581. SCOK1:
  582. cmpil #0x4004BC7E,d0 |...|X| < 15 PI?
  583. jlt  SCMAIN
  584. jra  REDUCEX
  587. SCMAIN:
  588. |--THIS IS THE USUAL CASE, |X| <= 15 PI.
  589. |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
  590. fmovex fp0,fp1
  591. fmuld TWOBYPI,fp1 |...X*2/PI
  593. |--HIDE THE NEXT THREE INSTRUCTIONS
  594. lea __x_PITBL+0x200,a1 |...TABLE OF N*PI/2, N = -32,...,32
  597. |--FP1 IS NOW READY
  598. fmovel fp1,a6@(N) |...CONVERT TO INTEGER
  600. movel a6@(N),d0
  601. asll #4,d0
  602. addal d0,a1 |...ADDRESS OF N*PIBY2, IN Y1, Y2
  604. fsubx A1@+,fp0 |...X-Y1
  605.         fsubs A1@,fp0 |...FP0 IS R = (X-Y1)-Y2
  607. SCCONT:
  608. |--continuation point from REDUCEX
  610. |--HIDE THE NEXT TWO
  611. movel a6@(N),d0
  612. rorl #1,d0
  614. cmpil #0,d0 |...D0 < 0 IFF N IS ODD
  615. jge  NEVEN
  617. NODD:
  618. |--REGISTERS SAVED SO FAR: D0, A0, FP2.
  620. fmovex fp0,a6@(RPRIME)
  621. fmulx fp0,fp0   |...FP0 IS S = R*R
  622. fmoved SINA7,fp1 |...A7
  623. fmoved COSB8,fp2 |...B8
  624. fmulx fp0,fp1   |...SA7
  625. movel d2,A7@-
  626. movel d0,d2
  627. fmulx fp0,fp2   |...SB8
  628. rorl #1,d2
  629. andil #0x80000000,d2
  631. faddd SINA6,fp1 |...A6+SA7
  632. eorl d0,d2
  633. andil #0x80000000,d2
  634. faddd COSB7,fp2 |...B7+SB8
  636. fmulx fp0,fp1   |...S(A6+SA7)
  637. eorl d2,a6@(RPRIME)
  638. movel A7@+,d2
  639. fmulx fp0,fp2   |...S(B7+SB8)
  640. rorl #1,d0
  641. andil #0x80000000,d0
  643. faddd SINA5,fp1 |...A5+S(A6+SA7)
  644. movel #0x3F800000,a6@(POSNEG1)
  645. eorl d0,a6@(POSNEG1)
  646. faddd COSB6,fp2 |...B6+S(B7+SB8)
  648. fmulx fp0,fp1   |...S(A5+S(A6+SA7))
  649. fmulx fp0,fp2   |...S(B6+S(B7+SB8))
  650. fmovex fp0,a6@(SPRIME)
  652. faddd SINA4,fp1 |...A4+S(A5+S(A6+SA7))
  653. eorl d0,a6@(SPRIME)
  654. faddd COSB5,fp2 |...B5+S(B6+S(B7+SB8))
  656. fmulx fp0,fp1   |...S(A4+...)
  657. fmulx fp0,fp2   |...S(B5+...)
  659. faddd SINA3,fp1 |...A3+S(A4+...)
  660. faddd COSB4,fp2 |...B4+S(B5+...)
  662. fmulx fp0,fp1   |...S(A3+...)
  663. fmulx fp0,fp2   |...S(B4+...)
  665. faddx SINA2,fp1 |...A2+S(A3+...)
  666. faddx COSB3,fp2 |...B3+S(B4+...)
  668. fmulx fp0,fp1   |...S(A2+...)
  669. fmulx fp0,fp2   |...S(B3+...)
  671. faddx SINA1,fp1 |...A1+S(A2+...)
  672. faddx COSB2,fp2 |...B2+S(B3+...)
  674. fmulx fp0,fp1   |...S(A1+...)
  675. fmulx fp2,fp0   |...S(B2+...)
  679. fmulx a6@(RPRIME),fp1 /* |...R'S(A1+...) */
  680. fadds COSB1,fp0 |...B1+S(B2...)
  681. fmulx a6@(SPRIME),fp0 /* |...S'(B1+S(B2+...)) */
  683. movel d1,a7@- | restore users mode # precision
  684. andil #0xff,d1 | mask off all exceptions
  685. fmovel d1,fpcr
  686. faddx a6@(RPRIME),fp1 |...COS(X)
  687. bsrl __x_sto_cos | store cosine result
  688. fmovel a7@+,fpcr | restore users exceptions
  689. fadds a6@(POSNEG1),fp0 |...SIN(X)
  691. jra  __x_t_frcinx
  694. NEVEN:
  695. |--REGISTERS SAVED SO FAR: FP2.
  697. fmovex fp0,a6@(RPRIME)
  698. fmulx fp0,fp0   |...FP0 IS S = R*R
  699. fmoved COSB8,fp1 |...B8
  700. fmoved SINA7,fp2 |...A7
  701. fmulx fp0,fp1   |...SB8
  702. fmovex fp0,a6@(SPRIME)
  703. fmulx fp0,fp2   |...SA7
  704. rorl #1,d0
  705. andil #0x80000000,d0
  706. faddd COSB7,fp1 |...B7+SB8
  707. faddd SINA6,fp2 |...A6+SA7
  708. eorl d0,a6@(RPRIME)
  709. eorl d0,a6@(SPRIME)
  710. fmulx fp0,fp1   |...S(B7+SB8)
  711. oril #0x3F800000,d0
  712. movel d0,a6@(POSNEG1)
  713. fmulx fp0,fp2   |...S(A6+SA7)
  715. faddd COSB6,fp1 |...B6+S(B7+SB8)
  716. faddd SINA5,fp2 |...A5+S(A6+SA7)
  718. fmulx fp0,fp1   |...S(B6+S(B7+SB8))
  719. fmulx fp0,fp2   |...S(A5+S(A6+SA7))
  721. faddd COSB5,fp1 |...B5+S(B6+S(B7+SB8))
  722. faddd SINA4,fp2 |...A4+S(A5+S(A6+SA7))
  724. fmulx fp0,fp1   |...S(B5+...)
  725. fmulx fp0,fp2   |...S(A4+...)
  727. faddd COSB4,fp1 |...B4+S(B5+...)
  728. faddd SINA3,fp2 |...A3+S(A4+...)
  730. fmulx fp0,fp1   |...S(B4+...)
  731. fmulx fp0,fp2   |...S(A3+...)
  733. faddx COSB3,fp1 |...B3+S(B4+...)
  734. faddx SINA2,fp2 |...A2+S(A3+...)
  736. fmulx fp0,fp1   |...S(B3+...)
  737. fmulx fp0,fp2   |...S(A2+...)
  739. faddx COSB2,fp1 |...B2+S(B3+...)
  740. faddx SINA1,fp2 |...A1+S(A2+...)
  742. fmulx fp0,fp1   |...S(B2+...)
  743. fmulx fp2,fp0   |...s(a1+...)
  747. fadds COSB1,fp1 |...B1+S(B2...)
  748. fmulx a6@(RPRIME),fp0 /* |...R'S(A1+...) */
  749. fmulx a6@(SPRIME),fp1 /* |...S'(B1+S(B2+...)) */
  751. movel d1,a7@- | save users mode # precision
  752. andil #0xff,d1 | mask off all exceptions
  753. fmovel d1,fpcr
  754. fadds a6@(POSNEG1),fp1 |...COS(X)
  755. bsrl __x_sto_cos | store cosine result
  756. fmovel a7@+,fpcr | restore users exceptions
  757. faddx a6@(RPRIME),fp0 |...SIN(X)
  759. jra  __x_t_frcinx
  761. SCBORS:
  762. cmpil #0x3FFF8000,d0
  763. jgt  REDUCEX
  766. SCSM:
  767. movew #0x0000,a6@(XDCARE)
  768. .long 0xf23c4480,0x3f800000 /*  fmoves  &0x3F800000,fp1 */
  770. movel d1,a7@- | save users mode # precision
  771. andil #0xff,d1 | mask off all exceptions
  772. fmovel d1,fpcr
  773. .long 0xf23c44a8,0x00800000 /*  fsubs  &0x00800000,fp1 */
  774. bsrl __x_sto_cos | store cosine result
  775. fmovel a7@+,fpcr | restore users exceptions
  776. fmovex a6@(X),fp0
  777. jra  __x_t_frcinx
  779. | end