Lines Matching +full:ulp +full:- +full:0
9 | Input: Double-extended number X in location pointed to
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
30 | 1. If SIN is invoked, set AdjN := 0; otherwise, set AdjN := 1.
32 | 2. If |X| >= 15Pi or |X| < 2**(-40), go to 7.
35 | k = N mod 4, so in particular, k = 0,1,2,or 3. Overwrite
40 | 5. (k is odd) Set j := (k-1)/2, sgn := (-1)**j. Return sgn*cos(r)
45 | 6. (k is even) Set j := k/2, sgn := (-1)**j. Return sgn*sin(r)
52 | 8. (|X|<2**(-40)) If SIN is invoked, return X; otherwise return 1.
57 | 1. If |X| >= 15Pi or |X| < 2**(-40), go to 6.
60 | k = N mod 4, so in particular, k = 0,1,2,or 3.
64 | 4. (k is odd) Set j1 := (k-1)/2, j2 := j1 (EOR) (k mod 2), i.e.
66 | sgn1 := (-1)**j1, sgn2 := (-1)**j2.
71 | 5. (k is even) Set j1 := k/2, sgn1 := (-1)**j1.
78 | 7. (|X|<2**(-40)) SIN(X) = X and COS(X) = 1. Exit.
95 BOUNDS1: .long 0x3FD78000,0x4004BC7E
96 TWOBYPI: .long 0x3FE45F30,0x6DC9C883
98 SINA7: .long 0xBD6AAA77,0xCCC994F5
99 SINA6: .long 0x3DE61209,0x7AAE8DA1
101 SINA5: .long 0xBE5AE645,0x2A118AE4
102 SINA4: .long 0x3EC71DE3,0xA5341531
104 SINA3: .long 0xBF2A01A0,0x1A018B59,0x00000000,0x00000000
106 SINA2: .long 0x3FF80000,0x88888888,0x888859AF,0x00000000
108 SINA1: .long 0xBFFC0000,0xAAAAAAAA,0xAAAAAA99,0x00000000
110 COSB8: .long 0x3D2AC4D0,0xD6011EE3
111 COSB7: .long 0xBDA9396F,0x9F45AC19
113 COSB6: .long 0x3E21EED9,0x0612C972
114 COSB5: .long 0xBE927E4F,0xB79D9FCF
116 COSB4: .long 0x3EFA01A0,0x1A01D423,0x00000000,0x00000000
118 COSB3: .long 0xBFF50000,0xB60B60B6,0x0B61D438,0x00000000
120 COSB2: .long 0x3FFA0000,0xAAAAAAAA,0xAAAAAB5E
121 COSB1: .long 0xBF000000
123 INVTWOPI: .long 0x3FFC0000,0xA2F9836E,0x4E44152A
125 TWOPI1: .long 0x40010000,0xC90FDAA2,0x00000000,0x00000000
126 TWOPI2: .long 0x3FDF0000,0x85A308D4,0x00000000,0x00000000
153 |--SIN(X) = X FOR DENORMALIZED X
158 |--COS(X) = 1 FOR DENORMALIZED X
160 fmoves #0x3F800000,%fp0
164 fmovel #0,%fpsr
170 |--SET ADJN TO 0
171 movel #0,ADJN(%a6)
176 |--SET ADJN TO 1
180 |--SAVE FPCR, FP1. CHECK IF |X| IS TOO SMALL OR LARGE
187 andil #0x7FFFFFFF,%d0 | ...COMPACTIFY X
189 cmpil #0x3FD78000,%d0 | ...|X| >= 2**(-40)?
194 cmpil #0x4004BC7E,%d0 | ...|X| < 15 PI?
199 |--THIS IS THE USUAL CASE, |X| <= 15 PI.
200 |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
204 |--HIDE THE NEXT THREE INSTRUCTIONS
205 lea PITBL+0x200,%a1 | ...TABLE OF N*PI/2, N = -32,...,32
208 |--FP1 IS NOW READY
216 fsubx (%a1)+,%fp0 | ...X-Y1
217 |--HIDE THE NEXT ONE
218 fsubs (%a1),%fp0 | ...FP0 IS R = (X-Y1)-Y2
221 |--continuation from REDUCEX
223 |--GET N+ADJN AND SEE IF SIN(R) OR COS(R) IS NEEDED
227 cmpil #0,%d0
231 |--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
232 |--THEN WE RETURN SGN*SIN(R). SGN*SIN(R) IS COMPUTED BY
233 |--R' + R'*S*(A1 + S(A2 + S(A3 + S(A4 + ... + SA7)))), WHERE
234 |--R' = SGN*R, S=R*R. THIS CAN BE REWRITTEN AS
235 |--R' + R'*S*( [A1+T(A3+T(A5+TA7))] + [S(A2+T(A4+TA6))])
236 |--WHERE T=S*S.
237 |--NOTE THAT A3 THROUGH A7 ARE STORED IN DOUBLE PRECISION
238 |--WHILE A1 AND A2 ARE IN DOUBLE-EXTENDED FORMAT.
241 |---HIDE THE NEXT TWO WHILE WAITING FOR FP0
244 |--FP0 IS NOW READY
247 |--HIDE THE NEXT TWO WHILE WAITING FOR FP1
250 andil #0x80000000,%d0
273 |--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
274 |--FP2 RELEASED, RESTORE NOW AND TAKE FULL ADVANTAGE OF HIDING
277 fmulx %fp1,%fp0 | ...SIN(R')-R'
278 |--FP1 RELEASED.
281 faddx X(%a6),%fp0 |last inst - possible exception set
286 |--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
287 |--THEN WE RETURN SGN*COS(R). SGN*COS(R) IS COMPUTED BY
288 |--SGN + S'*(B1 + S(B2 + S(B3 + S(B4 + ... + SB8)))), WHERE
289 |--S=R*R AND S'=SGN*S. THIS CAN BE REWRITTEN AS
290 |--SGN + S'*([B1+T(B3+T(B5+TB7))] + [S(B2+T(B4+T(B6+TB8)))])
291 |--WHERE T=S*S.
292 |--NOTE THAT B4 THROUGH B8 ARE STORED IN DOUBLE PRECISION
293 |--WHILE B2 AND B3 ARE IN DOUBLE-EXTENDED FORMAT, B1 IS -1/2
294 |--AND IS THEREFORE STORED AS SINGLE PRECISION.
297 |---HIDE THE NEXT TWO WHILE WAITING FOR FP0
300 |--FP0 IS NOW READY
303 |--HIDE THE NEXT TWO WHILE WAITING FOR FP1
306 andil #0x80000000,%d0
310 |--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
312 andil #0x80000000,%d0
315 |--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
316 oril #0x3F800000,%d0 | ...D0 IS SGN IN SINGLE
335 |--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
336 |--FP2 RELEASED.
340 |--FP1 RELEASED
345 fadds POSNEG1(%a6),%fp0 |last inst - possible exception set
350 |--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION.
351 |--IF |X| < 2**(-40), RETURN X OR 1.
352 cmpil #0x3FFF8000,%d0
358 cmpil #0,%d0
362 movew #0x0000,XDCARE(%a6) | ...JUST IN CASE
364 fmovex X(%a6),%fp0 |last inst - possible exception set
369 fmoves #0x3F800000,%fp0
372 fsubs #0x00800000,%fp0 |last inst - possible exception set
377 |--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW.
378 |--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING
379 |--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE.
381 fmovemx %fp2-%fp5,-(%a7) | ...save FP2 through FP5
382 movel %d2,-(%a7)
383 fmoves #0x00000000,%fp1
384 |--If compact form of abs(arg) in d0=$7ffeffff, argument is so large that
385 |--there is a danger of unwanted overflow in first LOOP iteration. In this
386 |--case, reduce argument by one remainder step to make subsequent reduction
387 |--safe.
388 cmpil #0x7ffeffff,%d0 |is argument dangerously large?
390 movel #0x7ffe0000,FP_SCR2(%a6) |yes
392 movel #0xc90fdaa2,FP_SCR2+4(%a6)
395 movel #0x7fdc0000,FP_SCR3(%a6) |create low half of 2**16383*
397 movel #0x85a308d3,FP_SCR3+4(%a6)
400 orw #0x8000,FP_SCR2(%a6) |positive arg
401 orw #0x8000,FP_SCR3(%a6)
409 |--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4.
410 |--integer quotient will be stored in N
411 |--Intermediate remainder is 66-bit long; (R,r) in (FP0,FP1)
414 fmovex %fp0,INARG(%a6) | ...+-2**K * F, 1 <= F < 2
417 andil #0x00007FFF,%d0
418 subil #0x00003FFF,%d0 | ...D0 IS K
422 subil #27,%d0 | ...D0 IS L := K-27
423 movel #0,ENDFLAG(%a6)
426 clrl %d0 | ...D0 IS L := 0
430 |--FIND THE REMAINDER OF (R,r) W.R.T. 2**L * (PI/2). L IS SO CHOSEN
431 |--THAT INT( X * (2/PI) / 2**(L) ) < 2**29.
433 |--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63),
434 |--2**L * (PIby2_1), 2**L * (PIby2_2)
436 movel #0x00003FFE,%d2 | ...BIASED EXPO OF 2/PI
437 subl %d0,%d2 | ...BIASED EXPO OF 2**(-L)*(2/PI)
439 movel #0xA2F9836E,FP_SCR1+4(%a6)
440 movel #0x4E44152A,FP_SCR1+8(%a6)
441 movew %d2,FP_SCR1(%a6) | ...FP_SCR1 is 2**(-L)*(2/PI)
445 |--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN
446 |--FLOATING POINT FORMAT, THE TWO FMOVE'S FMOVE.L FP <--> N
447 |--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT
448 |--(SIGN(INARG)*2**63 + FP2) - SIGN(INARG)*2**63 WILL GIVE
449 |--US THE DESIRED VALUE IN FLOATING POINT.
451 |--HIDE SIX CYCLES OF INSTRUCTION
454 andil #0x80000000,%d2
455 oril #0x5F000000,%d2 | ...D2 IS SIGN(INARG)*2**63 IN SGL
459 addil #0x00003FFF,%d2 | ...BIASED EXPO OF 2**L * (PI/2)
461 |--FP2 IS READY
464 |--HIDE 4 CYCLES OF INSTRUCTION; creating 2**(L)*Piby2_1 and 2**(L)*Piby2_2
467 movel #0xC90FDAA2,FP_SCR2+4(%a6)
470 |--FP2 IS READY
473 addil #0x00003FDD,%d0
476 movel #0x85A308D3,FP_SCR3+4(%a6)
481 |--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and
482 |--P2 = 2**(L) * Piby2_2
488 |--we want P+p = W+w but |p| <= half ulp of P
489 |--Then, we need to compute A := R-P and a := r-p
491 fsubx %fp3,%fp4 | ...W-P
493 fsubx %fp3,%fp0 | ...FP0 is A := R - P
494 faddx %fp5,%fp4 | ...FP4 is p = (W-P)+w
497 fsubx %fp4,%fp1 | ...FP1 is a := r - p
499 |--Now we need to normalize (A,a) to "new (R,r)" where R+r = A+a but
500 |--|r| <= half ulp of R.
502 |--No need to calculate r if this is the last loop
503 cmpil #0,%d0
506 |--Need to calculate r
507 fsubx %fp0,%fp3 | ...A-R
508 faddx %fp3,%fp1 | ...FP1 is r := (A-R)+a
514 fmovemx (%a7)+,%fp2-%fp5
525 |--SIN AND COS OF X FOR DENORMALIZED X
527 fmoves #0x3F800000,%fp1
533 |--SET ADJN TO 4
541 andil #0x7FFFFFFF,%d0 | ...COMPACTIFY X
543 cmpil #0x3FD78000,%d0 | ...|X| >= 2**(-40)?
548 cmpil #0x4004BC7E,%d0 | ...|X| < 15 PI?
554 |--THIS IS THE USUAL CASE, |X| <= 15 PI.
555 |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
559 |--HIDE THE NEXT THREE INSTRUCTIONS
560 lea PITBL+0x200,%a1 | ...TABLE OF N*PI/2, N = -32,...,32
563 |--FP1 IS NOW READY
570 fsubx (%a1)+,%fp0 | ...X-Y1
571 fsubs (%a1),%fp0 | ...FP0 IS R = (X-Y1)-Y2
574 |--continuation point from REDUCEX
576 |--HIDE THE NEXT TWO
580 cmpil #0,%d0 | ...D0 < 0 IFF N IS ODD
584 |--REGISTERS SAVED SO FAR: D0, A0, FP2.
591 movel %d2,-(%a7)
595 andil #0x80000000,%d2
599 andil #0x80000000,%d2
607 andil #0x80000000,%d0
610 movel #0x3F800000,POSNEG1(%a6)
649 movel %d1,-(%sp) |restore users mode & precision
650 andil #0xff,%d1 |mask off all exceptions
661 |--REGISTERS SAVED SO FAR: FP2.
671 andil #0x80000000,%d0
677 oril #0x3F800000,%d0
717 movel %d1,-(%sp) |save users mode & precision
718 andil #0xff,%d1 |mask off all exceptions
728 cmpil #0x3FFF8000,%d0
733 movew #0x0000,XDCARE(%a6)
734 fmoves #0x3F800000,%fp1
736 movel %d1,-(%sp) |save users mode & precision
737 andil #0xff,%d1 |mask off all exceptions
739 fsubs #0x00800000,%fp1