ssin.S 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745
  1. |
  2. | ssin.sa 3.3 7/29/91
  3. |
  4. | The entry point sSIN computes the sine of an input argument
  5. | sCOS computes the cosine, and sSINCOS computes both. The
  6. | corresponding entry points with a "d" computes the same
  7. | corresponding function values for denormalized inputs.
  8. |
  9. | Input: Double-extended number X in location pointed to
  10. | by address register a0.
  11. |
  12. | Output: The function value sin(X) or cos(X) returned in Fp0 if SIN or
  13. | COS is requested. Otherwise, for SINCOS, sin(X) is returned
  14. | in Fp0, and cos(X) is returned in Fp1.
  15. |
  16. | Modifies: Fp0 for SIN or COS; both Fp0 and Fp1 for SINCOS.
  17. |
  18. | Accuracy and Monotonicity: The returned result is within 1 ulp in
  19. | 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
  20. | result is subsequently rounded to double precision. The
  21. | result is provably monotonic in double precision.
  22. |
  23. | Speed: The programs sSIN and sCOS take approximately 150 cycles for
  24. | input argument X such that |X| < 15Pi, which is the usual
  25. | situation. The speed for sSINCOS is approximately 190 cycles.
  26. |
  27. | Algorithm:
  28. |
  29. | SIN and COS:
  30. | 1. If SIN is invoked, set AdjN := 0; otherwise, set AdjN := 1.
  31. |
  32. | 2. If |X| >= 15Pi or |X| < 2**(-40), go to 7.
  33. |
  34. | 3. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
  35. | k = N mod 4, so in particular, k = 0,1,2,or 3. Overwrite
  36. | k by k := k + AdjN.
  37. |
  38. | 4. If k is even, go to 6.
  39. |
  40. | 5. (k is odd) Set j := (k-1)/2, sgn := (-1)**j. Return sgn*cos(r)
  41. | where cos(r) is approximated by an even polynomial in r,
  42. | 1 + r*r*(B1+s*(B2+ ... + s*B8)), s = r*r.
  43. | Exit.
  44. |
  45. | 6. (k is even) Set j := k/2, sgn := (-1)**j. Return sgn*sin(r)
  46. | where sin(r) is approximated by an odd polynomial in r
  47. | r + r*s*(A1+s*(A2+ ... + s*A7)), s = r*r.
  48. | Exit.
  49. |
  50. | 7. If |X| > 1, go to 9.
  51. |
  52. | 8. (|X|<2**(-40)) If SIN is invoked, return X; otherwise return 1.
  53. |
  54. | 9. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 3.
  55. |
  56. | SINCOS:
  57. | 1. If |X| >= 15Pi or |X| < 2**(-40), go to 6.
  58. |
  59. | 2. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
  60. | k = N mod 4, so in particular, k = 0,1,2,or 3.
  61. |
  62. | 3. If k is even, go to 5.
  63. |
  64. | 4. (k is odd) Set j1 := (k-1)/2, j2 := j1 (EOR) (k mod 2), i.e.
  65. | j1 exclusive or with the l.s.b. of k.
  66. | sgn1 := (-1)**j1, sgn2 := (-1)**j2.
  67. | SIN(X) = sgn1 * cos(r) and COS(X) = sgn2*sin(r) where
  68. | sin(r) and cos(r) are computed as odd and even polynomials
  69. | in r, respectively. Exit
  70. |
  71. | 5. (k is even) Set j1 := k/2, sgn1 := (-1)**j1.
  72. | SIN(X) = sgn1 * sin(r) and COS(X) = sgn1*cos(r) where
  73. | sin(r) and cos(r) are computed as odd and even polynomials
  74. | in r, respectively. Exit
  75. |
  76. | 6. If |X| > 1, go to 8.
  77. |
  78. | 7. (|X|<2**(-40)) SIN(X) = X and COS(X) = 1. Exit.
  79. |
  80. | 8. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 2.
  81. |
  82. | Copyright (C) Motorola, Inc. 1990
  83. | All Rights Reserved
  84. |
  85. | For details on the license for this file, please see the
  86. | file, README, in this same directory.
  87. |SSIN idnt 2,1 | Motorola 040 Floating Point Software Package
  88. |section 8
  89. #include "fpsp.h"
  90. BOUNDS1: .long 0x3FD78000,0x4004BC7E
  91. TWOBYPI: .long 0x3FE45F30,0x6DC9C883
  92. SINA7: .long 0xBD6AAA77,0xCCC994F5
  93. SINA6: .long 0x3DE61209,0x7AAE8DA1
  94. SINA5: .long 0xBE5AE645,0x2A118AE4
  95. SINA4: .long 0x3EC71DE3,0xA5341531
  96. SINA3: .long 0xBF2A01A0,0x1A018B59,0x00000000,0x00000000
  97. SINA2: .long 0x3FF80000,0x88888888,0x888859AF,0x00000000
  98. SINA1: .long 0xBFFC0000,0xAAAAAAAA,0xAAAAAA99,0x00000000
  99. COSB8: .long 0x3D2AC4D0,0xD6011EE3
  100. COSB7: .long 0xBDA9396F,0x9F45AC19
  101. COSB6: .long 0x3E21EED9,0x0612C972
  102. COSB5: .long 0xBE927E4F,0xB79D9FCF
  103. COSB4: .long 0x3EFA01A0,0x1A01D423,0x00000000,0x00000000
  104. COSB3: .long 0xBFF50000,0xB60B60B6,0x0B61D438,0x00000000
  105. COSB2: .long 0x3FFA0000,0xAAAAAAAA,0xAAAAAB5E
  106. COSB1: .long 0xBF000000
  107. INVTWOPI: .long 0x3FFC0000,0xA2F9836E,0x4E44152A
  108. TWOPI1: .long 0x40010000,0xC90FDAA2,0x00000000,0x00000000
  109. TWOPI2: .long 0x3FDF0000,0x85A308D4,0x00000000,0x00000000
  110. |xref PITBL
  111. .set INARG,FP_SCR4
  112. .set X,FP_SCR5
  113. .set XDCARE,X+2
  114. .set XFRAC,X+4
  115. .set RPRIME,FP_SCR1
  116. .set SPRIME,FP_SCR2
  117. .set POSNEG1,L_SCR1
  118. .set TWOTO63,L_SCR1
  119. .set ENDFLAG,L_SCR2
  120. .set N,L_SCR2
  121. .set ADJN,L_SCR3
  122. | xref t_frcinx
  123. |xref t_extdnrm
  124. |xref sto_cos
  125. .global ssind
  126. ssind:
  127. |--SIN(X) = X FOR DENORMALIZED X
  128. bra t_extdnrm
  129. .global scosd
  130. scosd:
  131. |--COS(X) = 1 FOR DENORMALIZED X
  132. fmoves #0x3F800000,%fp0
  133. |
  134. | 9D25B Fix: Sometimes the previous fmove.s sets fpsr bits
  135. |
  136. fmovel #0,%fpsr
  137. |
  138. bra t_frcinx
  139. .global ssin
  140. ssin:
  141. |--SET ADJN TO 0
  142. movel #0,ADJN(%a6)
  143. bras SINBGN
  144. .global scos
  145. scos:
  146. |--SET ADJN TO 1
  147. movel #1,ADJN(%a6)
  148. SINBGN:
  149. |--SAVE FPCR, FP1. CHECK IF |X| IS TOO SMALL OR LARGE
  150. fmovex (%a0),%fp0 | ...LOAD INPUT
  151. movel (%a0),%d0
  152. movew 4(%a0),%d0
  153. fmovex %fp0,X(%a6)
  154. andil #0x7FFFFFFF,%d0 | ...COMPACTIFY X
  155. cmpil #0x3FD78000,%d0 | ...|X| >= 2**(-40)?
  156. bges SOK1
  157. bra SINSM
  158. SOK1:
  159. cmpil #0x4004BC7E,%d0 | ...|X| < 15 PI?
  160. blts SINMAIN
  161. bra REDUCEX
  162. SINMAIN:
  163. |--THIS IS THE USUAL CASE, |X| <= 15 PI.
  164. |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
  165. fmovex %fp0,%fp1
  166. fmuld TWOBYPI,%fp1 | ...X*2/PI
  167. |--HIDE THE NEXT THREE INSTRUCTIONS
  168. lea PITBL+0x200,%a1 | ...TABLE OF N*PI/2, N = -32,...,32
  169. |--FP1 IS NOW READY
  170. fmovel %fp1,N(%a6) | ...CONVERT TO INTEGER
  171. movel N(%a6),%d0
  172. asll #4,%d0
  173. addal %d0,%a1 | ...A1 IS THE ADDRESS OF N*PIBY2
  174. | ...WHICH IS IN TWO PIECES Y1 & Y2
  175. fsubx (%a1)+,%fp0 | ...X-Y1
  176. |--HIDE THE NEXT ONE
  177. fsubs (%a1),%fp0 | ...FP0 IS R = (X-Y1)-Y2
  178. SINCONT:
  179. |--continuation from REDUCEX
  180. |--GET N+ADJN AND SEE IF SIN(R) OR COS(R) IS NEEDED
  181. movel N(%a6),%d0
  182. addl ADJN(%a6),%d0 | ...SEE IF D0 IS ODD OR EVEN
  183. rorl #1,%d0 | ...D0 WAS ODD IFF D0 IS NEGATIVE
  184. cmpil #0,%d0
  185. blt COSPOLY
  186. SINPOLY:
  187. |--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
  188. |--THEN WE RETURN SGN*SIN(R). SGN*SIN(R) IS COMPUTED BY
  189. |--R' + R'*S*(A1 + S(A2 + S(A3 + S(A4 + ... + SA7)))), WHERE
  190. |--R' = SGN*R, S=R*R. THIS CAN BE REWRITTEN AS
  191. |--R' + R'*S*( [A1+T(A3+T(A5+TA7))] + [S(A2+T(A4+TA6))])
  192. |--WHERE T=S*S.
  193. |--NOTE THAT A3 THROUGH A7 ARE STORED IN DOUBLE PRECISION
  194. |--WHILE A1 AND A2 ARE IN DOUBLE-EXTENDED FORMAT.
  195. fmovex %fp0,X(%a6) | ...X IS R
  196. fmulx %fp0,%fp0 | ...FP0 IS S
  197. |---HIDE THE NEXT TWO WHILE WAITING FOR FP0
  198. fmoved SINA7,%fp3
  199. fmoved SINA6,%fp2
  200. |--FP0 IS NOW READY
  201. fmovex %fp0,%fp1
  202. fmulx %fp1,%fp1 | ...FP1 IS T
  203. |--HIDE THE NEXT TWO WHILE WAITING FOR FP1
  204. rorl #1,%d0
  205. andil #0x80000000,%d0
  206. | ...LEAST SIG. BIT OF D0 IN SIGN POSITION
  207. eorl %d0,X(%a6) | ...X IS NOW R'= SGN*R
  208. fmulx %fp1,%fp3 | ...TA7
  209. fmulx %fp1,%fp2 | ...TA6
  210. faddd SINA5,%fp3 | ...A5+TA7
  211. faddd SINA4,%fp2 | ...A4+TA6
  212. fmulx %fp1,%fp3 | ...T(A5+TA7)
  213. fmulx %fp1,%fp2 | ...T(A4+TA6)
  214. faddd SINA3,%fp3 | ...A3+T(A5+TA7)
  215. faddx SINA2,%fp2 | ...A2+T(A4+TA6)
  216. fmulx %fp3,%fp1 | ...T(A3+T(A5+TA7))
  217. fmulx %fp0,%fp2 | ...S(A2+T(A4+TA6))
  218. faddx SINA1,%fp1 | ...A1+T(A3+T(A5+TA7))
  219. fmulx X(%a6),%fp0 | ...R'*S
  220. faddx %fp2,%fp1 | ...[A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))]
  221. |--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
  222. |--FP2 RELEASED, RESTORE NOW AND TAKE FULL ADVANTAGE OF HIDING
  223. fmulx %fp1,%fp0 | ...SIN(R')-R'
  224. |--FP1 RELEASED.
  225. fmovel %d1,%FPCR |restore users exceptions
  226. faddx X(%a6),%fp0 |last inst - possible exception set
  227. bra t_frcinx
  228. COSPOLY:
  229. |--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
  230. |--THEN WE RETURN SGN*COS(R). SGN*COS(R) IS COMPUTED BY
  231. |--SGN + S'*(B1 + S(B2 + S(B3 + S(B4 + ... + SB8)))), WHERE
  232. |--S=R*R AND S'=SGN*S. THIS CAN BE REWRITTEN AS
  233. |--SGN + S'*([B1+T(B3+T(B5+TB7))] + [S(B2+T(B4+T(B6+TB8)))])
  234. |--WHERE T=S*S.
  235. |--NOTE THAT B4 THROUGH B8 ARE STORED IN DOUBLE PRECISION
  236. |--WHILE B2 AND B3 ARE IN DOUBLE-EXTENDED FORMAT, B1 IS -1/2
  237. |--AND IS THEREFORE STORED AS SINGLE PRECISION.
  238. fmulx %fp0,%fp0 | ...FP0 IS S
  239. |---HIDE THE NEXT TWO WHILE WAITING FOR FP0
  240. fmoved COSB8,%fp2
  241. fmoved COSB7,%fp3
  242. |--FP0 IS NOW READY
  243. fmovex %fp0,%fp1
  244. fmulx %fp1,%fp1 | ...FP1 IS T
  245. |--HIDE THE NEXT TWO WHILE WAITING FOR FP1
  246. fmovex %fp0,X(%a6) | ...X IS S
  247. rorl #1,%d0
  248. andil #0x80000000,%d0
  249. | ...LEAST SIG. BIT OF D0 IN SIGN POSITION
  250. fmulx %fp1,%fp2 | ...TB8
  251. |--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
  252. eorl %d0,X(%a6) | ...X IS NOW S'= SGN*S
  253. andil #0x80000000,%d0
  254. fmulx %fp1,%fp3 | ...TB7
  255. |--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
  256. oril #0x3F800000,%d0 | ...D0 IS SGN IN SINGLE
  257. movel %d0,POSNEG1(%a6)
  258. faddd COSB6,%fp2 | ...B6+TB8
  259. faddd COSB5,%fp3 | ...B5+TB7
  260. fmulx %fp1,%fp2 | ...T(B6+TB8)
  261. fmulx %fp1,%fp3 | ...T(B5+TB7)
  262. faddd COSB4,%fp2 | ...B4+T(B6+TB8)
  263. faddx COSB3,%fp3 | ...B3+T(B5+TB7)
  264. fmulx %fp1,%fp2 | ...T(B4+T(B6+TB8))
  265. fmulx %fp3,%fp1 | ...T(B3+T(B5+TB7))
  266. faddx COSB2,%fp2 | ...B2+T(B4+T(B6+TB8))
  267. fadds COSB1,%fp1 | ...B1+T(B3+T(B5+TB7))
  268. fmulx %fp2,%fp0 | ...S(B2+T(B4+T(B6+TB8)))
  269. |--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
  270. |--FP2 RELEASED.
  271. faddx %fp1,%fp0
  272. |--FP1 RELEASED
  273. fmulx X(%a6),%fp0
  274. fmovel %d1,%FPCR |restore users exceptions
  275. fadds POSNEG1(%a6),%fp0 |last inst - possible exception set
  276. bra t_frcinx
  277. SINBORS:
  278. |--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION.
  279. |--IF |X| < 2**(-40), RETURN X OR 1.
  280. cmpil #0x3FFF8000,%d0
  281. bgts REDUCEX
  282. SINSM:
  283. movel ADJN(%a6),%d0
  284. cmpil #0,%d0
  285. bgts COSTINY
  286. SINTINY:
  287. movew #0x0000,XDCARE(%a6) | ...JUST IN CASE
  288. fmovel %d1,%FPCR |restore users exceptions
  289. fmovex X(%a6),%fp0 |last inst - possible exception set
  290. bra t_frcinx
  291. COSTINY:
  292. fmoves #0x3F800000,%fp0
  293. fmovel %d1,%FPCR |restore users exceptions
  294. fsubs #0x00800000,%fp0 |last inst - possible exception set
  295. bra t_frcinx
  296. REDUCEX:
  297. |--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW.
  298. |--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING
  299. |--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE.
  300. fmovemx %fp2-%fp5,-(%a7) | ...save FP2 through FP5
  301. movel %d2,-(%a7)
  302. fmoves #0x00000000,%fp1
  303. |--If compact form of abs(arg) in d0=$7ffeffff, argument is so large that
  304. |--there is a danger of unwanted overflow in first LOOP iteration. In this
  305. |--case, reduce argument by one remainder step to make subsequent reduction
  306. |--safe.
  307. cmpil #0x7ffeffff,%d0 |is argument dangerously large?
  308. bnes LOOP
  309. movel #0x7ffe0000,FP_SCR2(%a6) |yes
  310. | ;create 2**16383*PI/2
  311. movel #0xc90fdaa2,FP_SCR2+4(%a6)
  312. clrl FP_SCR2+8(%a6)
  313. ftstx %fp0 |test sign of argument
  314. movel #0x7fdc0000,FP_SCR3(%a6) |create low half of 2**16383*
  315. | ;PI/2 at FP_SCR3
  316. movel #0x85a308d3,FP_SCR3+4(%a6)
  317. clrl FP_SCR3+8(%a6)
  318. fblt red_neg
  319. orw #0x8000,FP_SCR2(%a6) |positive arg
  320. orw #0x8000,FP_SCR3(%a6)
  321. red_neg:
  322. faddx FP_SCR2(%a6),%fp0 |high part of reduction is exact
  323. fmovex %fp0,%fp1 |save high result in fp1
  324. faddx FP_SCR3(%a6),%fp0 |low part of reduction
  325. fsubx %fp0,%fp1 |determine low component of result
  326. faddx FP_SCR3(%a6),%fp1 |fp0/fp1 are reduced argument.
  327. |--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4.
  328. |--integer quotient will be stored in N
  329. |--Intermediate remainder is 66-bit long; (R,r) in (FP0,FP1)
  330. LOOP:
  331. fmovex %fp0,INARG(%a6) | ...+-2**K * F, 1 <= F < 2
  332. movew INARG(%a6),%d0
  333. movel %d0,%a1 | ...save a copy of D0
  334. andil #0x00007FFF,%d0
  335. subil #0x00003FFF,%d0 | ...D0 IS K
  336. cmpil #28,%d0
  337. bles LASTLOOP
  338. CONTLOOP:
  339. subil #27,%d0 | ...D0 IS L := K-27
  340. movel #0,ENDFLAG(%a6)
  341. bras WORK
  342. LASTLOOP:
  343. clrl %d0 | ...D0 IS L := 0
  344. movel #1,ENDFLAG(%a6)
  345. WORK:
  346. |--FIND THE REMAINDER OF (R,r) W.R.T. 2**L * (PI/2). L IS SO CHOSEN
  347. |--THAT INT( X * (2/PI) / 2**(L) ) < 2**29.
  348. |--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63),
  349. |--2**L * (PIby2_1), 2**L * (PIby2_2)
  350. movel #0x00003FFE,%d2 | ...BIASED EXPO OF 2/PI
  351. subl %d0,%d2 | ...BIASED EXPO OF 2**(-L)*(2/PI)
  352. movel #0xA2F9836E,FP_SCR1+4(%a6)
  353. movel #0x4E44152A,FP_SCR1+8(%a6)
  354. movew %d2,FP_SCR1(%a6) | ...FP_SCR1 is 2**(-L)*(2/PI)
  355. fmovex %fp0,%fp2
  356. fmulx FP_SCR1(%a6),%fp2
  357. |--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN
  358. |--FLOATING POINT FORMAT, THE TWO FMOVE'S FMOVE.L FP <--> N
  359. |--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT
  360. |--(SIGN(INARG)*2**63 + FP2) - SIGN(INARG)*2**63 WILL GIVE
  361. |--US THE DESIRED VALUE IN FLOATING POINT.
  362. |--HIDE SIX CYCLES OF INSTRUCTION
  363. movel %a1,%d2
  364. swap %d2
  365. andil #0x80000000,%d2
  366. oril #0x5F000000,%d2 | ...D2 IS SIGN(INARG)*2**63 IN SGL
  367. movel %d2,TWOTO63(%a6)
  368. movel %d0,%d2
  369. addil #0x00003FFF,%d2 | ...BIASED EXPO OF 2**L * (PI/2)
  370. |--FP2 IS READY
  371. fadds TWOTO63(%a6),%fp2 | ...THE FRACTIONAL PART OF FP1 IS ROUNDED
  372. |--HIDE 4 CYCLES OF INSTRUCTION; creating 2**(L)*Piby2_1 and 2**(L)*Piby2_2
  373. movew %d2,FP_SCR2(%a6)
  374. clrw FP_SCR2+2(%a6)
  375. movel #0xC90FDAA2,FP_SCR2+4(%a6)
  376. clrl FP_SCR2+8(%a6) | ...FP_SCR2 is 2**(L) * Piby2_1
  377. |--FP2 IS READY
  378. fsubs TWOTO63(%a6),%fp2 | ...FP2 is N
  379. addil #0x00003FDD,%d0
  380. movew %d0,FP_SCR3(%a6)
  381. clrw FP_SCR3+2(%a6)
  382. movel #0x85A308D3,FP_SCR3+4(%a6)
  383. clrl FP_SCR3+8(%a6) | ...FP_SCR3 is 2**(L) * Piby2_2
  384. movel ENDFLAG(%a6),%d0
  385. |--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and
  386. |--P2 = 2**(L) * Piby2_2
  387. fmovex %fp2,%fp4
  388. fmulx FP_SCR2(%a6),%fp4 | ...W = N*P1
  389. fmovex %fp2,%fp5
  390. fmulx FP_SCR3(%a6),%fp5 | ...w = N*P2
  391. fmovex %fp4,%fp3
  392. |--we want P+p = W+w but |p| <= half ulp of P
  393. |--Then, we need to compute A := R-P and a := r-p
  394. faddx %fp5,%fp3 | ...FP3 is P
  395. fsubx %fp3,%fp4 | ...W-P
  396. fsubx %fp3,%fp0 | ...FP0 is A := R - P
  397. faddx %fp5,%fp4 | ...FP4 is p = (W-P)+w
  398. fmovex %fp0,%fp3 | ...FP3 A
  399. fsubx %fp4,%fp1 | ...FP1 is a := r - p
  400. |--Now we need to normalize (A,a) to "new (R,r)" where R+r = A+a but
  401. |--|r| <= half ulp of R.
  402. faddx %fp1,%fp0 | ...FP0 is R := A+a
  403. |--No need to calculate r if this is the last loop
  404. cmpil #0,%d0
  405. bgt RESTORE
  406. |--Need to calculate r
  407. fsubx %fp0,%fp3 | ...A-R
  408. faddx %fp3,%fp1 | ...FP1 is r := (A-R)+a
  409. bra LOOP
  410. RESTORE:
  411. fmovel %fp2,N(%a6)
  412. movel (%a7)+,%d2
  413. fmovemx (%a7)+,%fp2-%fp5
  414. movel ADJN(%a6),%d0
  415. cmpil #4,%d0
  416. blt SINCONT
  417. bras SCCONT
  418. .global ssincosd
  419. ssincosd:
  420. |--SIN AND COS OF X FOR DENORMALIZED X
  421. fmoves #0x3F800000,%fp1
  422. bsr sto_cos |store cosine result
  423. bra t_extdnrm
  424. .global ssincos
  425. ssincos:
  426. |--SET ADJN TO 4
  427. movel #4,ADJN(%a6)
  428. fmovex (%a0),%fp0 | ...LOAD INPUT
  429. movel (%a0),%d0
  430. movew 4(%a0),%d0
  431. fmovex %fp0,X(%a6)
  432. andil #0x7FFFFFFF,%d0 | ...COMPACTIFY X
  433. cmpil #0x3FD78000,%d0 | ...|X| >= 2**(-40)?
  434. bges SCOK1
  435. bra SCSM
  436. SCOK1:
  437. cmpil #0x4004BC7E,%d0 | ...|X| < 15 PI?
  438. blts SCMAIN
  439. bra REDUCEX
  440. SCMAIN:
  441. |--THIS IS THE USUAL CASE, |X| <= 15 PI.
  442. |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
  443. fmovex %fp0,%fp1
  444. fmuld TWOBYPI,%fp1 | ...X*2/PI
  445. |--HIDE THE NEXT THREE INSTRUCTIONS
  446. lea PITBL+0x200,%a1 | ...TABLE OF N*PI/2, N = -32,...,32
  447. |--FP1 IS NOW READY
  448. fmovel %fp1,N(%a6) | ...CONVERT TO INTEGER
  449. movel N(%a6),%d0
  450. asll #4,%d0
  451. addal %d0,%a1 | ...ADDRESS OF N*PIBY2, IN Y1, Y2
  452. fsubx (%a1)+,%fp0 | ...X-Y1
  453. fsubs (%a1),%fp0 | ...FP0 IS R = (X-Y1)-Y2
  454. SCCONT:
  455. |--continuation point from REDUCEX
  456. |--HIDE THE NEXT TWO
  457. movel N(%a6),%d0
  458. rorl #1,%d0
  459. cmpil #0,%d0 | ...D0 < 0 IFF N IS ODD
  460. bge NEVEN
  461. NODD:
  462. |--REGISTERS SAVED SO FAR: D0, A0, FP2.
  463. fmovex %fp0,RPRIME(%a6)
  464. fmulx %fp0,%fp0 | ...FP0 IS S = R*R
  465. fmoved SINA7,%fp1 | ...A7
  466. fmoved COSB8,%fp2 | ...B8
  467. fmulx %fp0,%fp1 | ...SA7
  468. movel %d2,-(%a7)
  469. movel %d0,%d2
  470. fmulx %fp0,%fp2 | ...SB8
  471. rorl #1,%d2
  472. andil #0x80000000,%d2
  473. faddd SINA6,%fp1 | ...A6+SA7
  474. eorl %d0,%d2
  475. andil #0x80000000,%d2
  476. faddd COSB7,%fp2 | ...B7+SB8
  477. fmulx %fp0,%fp1 | ...S(A6+SA7)
  478. eorl %d2,RPRIME(%a6)
  479. movel (%a7)+,%d2
  480. fmulx %fp0,%fp2 | ...S(B7+SB8)
  481. rorl #1,%d0
  482. andil #0x80000000,%d0
  483. faddd SINA5,%fp1 | ...A5+S(A6+SA7)
  484. movel #0x3F800000,POSNEG1(%a6)
  485. eorl %d0,POSNEG1(%a6)
  486. faddd COSB6,%fp2 | ...B6+S(B7+SB8)
  487. fmulx %fp0,%fp1 | ...S(A5+S(A6+SA7))
  488. fmulx %fp0,%fp2 | ...S(B6+S(B7+SB8))
  489. fmovex %fp0,SPRIME(%a6)
  490. faddd SINA4,%fp1 | ...A4+S(A5+S(A6+SA7))
  491. eorl %d0,SPRIME(%a6)
  492. faddd COSB5,%fp2 | ...B5+S(B6+S(B7+SB8))
  493. fmulx %fp0,%fp1 | ...S(A4+...)
  494. fmulx %fp0,%fp2 | ...S(B5+...)
  495. faddd SINA3,%fp1 | ...A3+S(A4+...)
  496. faddd COSB4,%fp2 | ...B4+S(B5+...)
  497. fmulx %fp0,%fp1 | ...S(A3+...)
  498. fmulx %fp0,%fp2 | ...S(B4+...)
  499. faddx SINA2,%fp1 | ...A2+S(A3+...)
  500. faddx COSB3,%fp2 | ...B3+S(B4+...)
  501. fmulx %fp0,%fp1 | ...S(A2+...)
  502. fmulx %fp0,%fp2 | ...S(B3+...)
  503. faddx SINA1,%fp1 | ...A1+S(A2+...)
  504. faddx COSB2,%fp2 | ...B2+S(B3+...)
  505. fmulx %fp0,%fp1 | ...S(A1+...)
  506. fmulx %fp2,%fp0 | ...S(B2+...)
  507. fmulx RPRIME(%a6),%fp1 | ...R'S(A1+...)
  508. fadds COSB1,%fp0 | ...B1+S(B2...)
  509. fmulx SPRIME(%a6),%fp0 | ...S'(B1+S(B2+...))
  510. movel %d1,-(%sp) |restore users mode & precision
  511. andil #0xff,%d1 |mask off all exceptions
  512. fmovel %d1,%FPCR
  513. faddx RPRIME(%a6),%fp1 | ...COS(X)
  514. bsr sto_cos |store cosine result
  515. fmovel (%sp)+,%FPCR |restore users exceptions
  516. fadds POSNEG1(%a6),%fp0 | ...SIN(X)
  517. bra t_frcinx
  518. NEVEN:
  519. |--REGISTERS SAVED SO FAR: FP2.
  520. fmovex %fp0,RPRIME(%a6)
  521. fmulx %fp0,%fp0 | ...FP0 IS S = R*R
  522. fmoved COSB8,%fp1 | ...B8
  523. fmoved SINA7,%fp2 | ...A7
  524. fmulx %fp0,%fp1 | ...SB8
  525. fmovex %fp0,SPRIME(%a6)
  526. fmulx %fp0,%fp2 | ...SA7
  527. rorl #1,%d0
  528. andil #0x80000000,%d0
  529. faddd COSB7,%fp1 | ...B7+SB8
  530. faddd SINA6,%fp2 | ...A6+SA7
  531. eorl %d0,RPRIME(%a6)
  532. eorl %d0,SPRIME(%a6)
  533. fmulx %fp0,%fp1 | ...S(B7+SB8)
  534. oril #0x3F800000,%d0
  535. movel %d0,POSNEG1(%a6)
  536. fmulx %fp0,%fp2 | ...S(A6+SA7)
  537. faddd COSB6,%fp1 | ...B6+S(B7+SB8)
  538. faddd SINA5,%fp2 | ...A5+S(A6+SA7)
  539. fmulx %fp0,%fp1 | ...S(B6+S(B7+SB8))
  540. fmulx %fp0,%fp2 | ...S(A5+S(A6+SA7))
  541. faddd COSB5,%fp1 | ...B5+S(B6+S(B7+SB8))
  542. faddd SINA4,%fp2 | ...A4+S(A5+S(A6+SA7))
  543. fmulx %fp0,%fp1 | ...S(B5+...)
  544. fmulx %fp0,%fp2 | ...S(A4+...)
  545. faddd COSB4,%fp1 | ...B4+S(B5+...)
  546. faddd SINA3,%fp2 | ...A3+S(A4+...)
  547. fmulx %fp0,%fp1 | ...S(B4+...)
  548. fmulx %fp0,%fp2 | ...S(A3+...)
  549. faddx COSB3,%fp1 | ...B3+S(B4+...)
  550. faddx SINA2,%fp2 | ...A2+S(A3+...)
  551. fmulx %fp0,%fp1 | ...S(B3+...)
  552. fmulx %fp0,%fp2 | ...S(A2+...)
  553. faddx COSB2,%fp1 | ...B2+S(B3+...)
  554. faddx SINA1,%fp2 | ...A1+S(A2+...)
  555. fmulx %fp0,%fp1 | ...S(B2+...)
  556. fmulx %fp2,%fp0 | ...s(a1+...)
  557. fadds COSB1,%fp1 | ...B1+S(B2...)
  558. fmulx RPRIME(%a6),%fp0 | ...R'S(A1+...)
  559. fmulx SPRIME(%a6),%fp1 | ...S'(B1+S(B2+...))
  560. movel %d1,-(%sp) |save users mode & precision
  561. andil #0xff,%d1 |mask off all exceptions
  562. fmovel %d1,%FPCR
  563. fadds POSNEG1(%a6),%fp1 | ...COS(X)
  564. bsr sto_cos |store cosine result
  565. fmovel (%sp)+,%FPCR |restore users exceptions
  566. faddx RPRIME(%a6),%fp0 | ...SIN(X)
  567. bra t_frcinx
  568. SCBORS:
  569. cmpil #0x3FFF8000,%d0
  570. bgt REDUCEX
  571. SCSM:
  572. movew #0x0000,XDCARE(%a6)
  573. fmoves #0x3F800000,%fp1
  574. movel %d1,-(%sp) |save users mode & precision
  575. andil #0xff,%d1 |mask off all exceptions
  576. fmovel %d1,%FPCR
  577. fsubs #0x00800000,%fp1
  578. bsr sto_cos |store cosine result
  579. fmovel (%sp)+,%FPCR |restore users exceptions
  580. fmovex X(%a6),%fp0
  581. bra t_frcinx
  582. |end