Lines Matching +full:half +full:- +full:precision
3 M68000 Hi-Performance Microprocessor Division
5 Production Release P1.00 -- October 10, 1994
276 set LV, -LOCAL_SIZE # stack offset
285 set EXC_AREGS, -68 # offset of all address regs
286 set EXC_DREGS, -100 # offset of all data regs
287 set EXC_FPREGS, -36 # offset of all fp regs
372 set FTEMP_EX, 0 # extended precision
379 set LOCAL_EX, 0 # extended precision
386 set DST_EX, 0 # extended precision
391 set SRC_EX, 0 # extended precision
402 set EXT_BIAS, 0x3fff # extended precision bias
403 set SGL_BIAS, 0x007f # single precision bias
404 set DBL_BIAS, 0x03ff # double precision bias
499 set x_mode, 0x0 # extended precision
500 set s_mode, 0x4 # single precision
501 set d_mode, 0x8 # double precision
503 set rn_mode, 0x0 # round-to-nearest
504 set rz_mode, 0x1 # round-to-zero
505 set rm_mode, 0x2 # round-tp-minus-infinity
506 set rp_mode, 0x3 # round-to-plus-infinity
528 set mda7_flg, 0x08 # flag bit: -(a7) <ea>
539 # TRANSCENDENTAL "LAST-OP" FLAGS #
563 link %a6,&-LOCAL_SIZE
565 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
612 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
620 link %a6,&-LOCAL_SIZE
622 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
670 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
678 link %a6,&-LOCAL_SIZE
680 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
728 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
740 link %a6,&-LOCAL_SIZE
742 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
789 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
797 link %a6,&-LOCAL_SIZE
799 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
847 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
855 link %a6,&-LOCAL_SIZE
857 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
905 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
917 link %a6,&-LOCAL_SIZE
919 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
966 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
974 link %a6,&-LOCAL_SIZE
976 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1024 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1032 link %a6,&-LOCAL_SIZE
1034 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1082 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1094 link %a6,&-LOCAL_SIZE
1096 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1143 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1151 link %a6,&-LOCAL_SIZE
1153 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1201 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1209 link %a6,&-LOCAL_SIZE
1211 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1259 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1271 link %a6,&-LOCAL_SIZE
1273 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1320 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1328 link %a6,&-LOCAL_SIZE
1330 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1378 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1386 link %a6,&-LOCAL_SIZE
1388 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1436 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1448 link %a6,&-LOCAL_SIZE
1450 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1497 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1505 link %a6,&-LOCAL_SIZE
1507 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1555 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1563 link %a6,&-LOCAL_SIZE
1565 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1613 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1625 link %a6,&-LOCAL_SIZE
1627 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1674 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1682 link %a6,&-LOCAL_SIZE
1684 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1732 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1740 link %a6,&-LOCAL_SIZE
1742 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1790 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1802 link %a6,&-LOCAL_SIZE
1804 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1851 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1859 link %a6,&-LOCAL_SIZE
1861 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1909 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1917 link %a6,&-LOCAL_SIZE
1919 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
1967 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
1979 link %a6,&-LOCAL_SIZE
1981 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2028 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2036 link %a6,&-LOCAL_SIZE
2038 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2086 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2094 link %a6,&-LOCAL_SIZE
2096 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2144 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2156 link %a6,&-LOCAL_SIZE
2158 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2205 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2213 link %a6,&-LOCAL_SIZE
2215 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2263 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2271 link %a6,&-LOCAL_SIZE
2273 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2321 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2333 link %a6,&-LOCAL_SIZE
2335 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2382 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2390 link %a6,&-LOCAL_SIZE
2392 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2440 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2448 link %a6,&-LOCAL_SIZE
2450 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2498 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2510 link %a6,&-LOCAL_SIZE
2512 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2559 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2567 link %a6,&-LOCAL_SIZE
2569 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2617 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2625 link %a6,&-LOCAL_SIZE
2627 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2675 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2687 link %a6,&-LOCAL_SIZE
2689 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2736 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2744 link %a6,&-LOCAL_SIZE
2746 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2794 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2802 link %a6,&-LOCAL_SIZE
2804 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2852 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2864 link %a6,&-LOCAL_SIZE
2866 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2913 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2921 link %a6,&-LOCAL_SIZE
2923 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
2971 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
2979 link %a6,&-LOCAL_SIZE
2981 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3029 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3041 link %a6,&-LOCAL_SIZE
3043 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3090 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3098 link %a6,&-LOCAL_SIZE
3100 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3148 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3156 link %a6,&-LOCAL_SIZE
3158 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3206 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3218 link %a6,&-LOCAL_SIZE
3220 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3267 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3275 link %a6,&-LOCAL_SIZE
3277 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3325 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3333 link %a6,&-LOCAL_SIZE
3335 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3383 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3395 link %a6,&-LOCAL_SIZE
3397 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3444 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3452 link %a6,&-LOCAL_SIZE
3454 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3502 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3510 link %a6,&-LOCAL_SIZE
3512 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3560 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3572 link %a6,&-LOCAL_SIZE
3574 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3621 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3629 link %a6,&-LOCAL_SIZE
3631 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3679 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3687 link %a6,&-LOCAL_SIZE
3689 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3737 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3749 link %a6,&-LOCAL_SIZE
3751 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3798 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3806 link %a6,&-LOCAL_SIZE
3808 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3856 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3864 link %a6,&-LOCAL_SIZE
3866 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3914 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3926 link %a6,&-LOCAL_SIZE
3928 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
3975 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
3983 link %a6,&-LOCAL_SIZE
3985 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4033 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4041 link %a6,&-LOCAL_SIZE
4043 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4091 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4103 link %a6,&-LOCAL_SIZE
4105 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4152 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4154 fmovm.x &0x03,-(%sp) # store off fp0/fp1
4162 link %a6,&-LOCAL_SIZE
4164 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4212 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4214 fmovm.x &0x03,-(%sp) # store off fp0/fp1
4222 link %a6,&-LOCAL_SIZE
4224 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4272 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4274 fmovm.x &0x03,-(%sp) # store off fp0/fp1
4286 link %a6,&-LOCAL_SIZE
4288 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4344 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4352 link %a6,&-LOCAL_SIZE
4354 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4410 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4418 link %a6,&-LOCAL_SIZE
4420 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4478 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4490 link %a6,&-LOCAL_SIZE
4492 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4548 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4556 link %a6,&-LOCAL_SIZE
4558 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4614 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4622 link %a6,&-LOCAL_SIZE
4624 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4682 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4694 link %a6,&-LOCAL_SIZE
4696 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4752 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4760 link %a6,&-LOCAL_SIZE
4762 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4818 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4826 link %a6,&-LOCAL_SIZE
4828 movm.l &0x0303,EXC_DREGS(%a6) # save d0-d1/a0-a1
4886 movm.l EXC_DREGS(%a6),&0x0303 # restore d0-d1/a0-a1
4902 # a0 = pointer to extended precision input #
4903 # d0 = round precision,mode #
4915 # rounded to double precision. The result is provably monotonic #
4916 # in double precision. #
4923 # 2. If |X| >= 15Pi or |X| < 2**(-40), go to 7. #
4931 # 5. (k is odd) Set j := (k-1)/2, sgn := (-1)**j. #
4937 # 6. (k is even) Set j := k/2, sgn := (-1)**j. Return sgn*sin(r) #
4944 # 8. (|X|<2**(-40)) If SIN is invoked, return X; #
4951 # 1. If |X| >= 15Pi or |X| < 2**(-40), go to 6. #
4958 # 4. (k is odd) Set j1 := (k-1)/2, j2 := j1 (EOR) (k mod 2), ie. #
4960 # sgn1 := (-1)**j1, sgn2 := (-1)**j2. #
4965 # 5. (k is even) Set j1 := k/2, sgn1 := (-1)**j1. #
4972 # 7. (|X|<2**(-40)) SIN(X) = X and COS(X) = 1. Exit. #
5026 #--SAVE FPCR, FP1. CHECK IF |X| IS TOO SMALL OR LARGE
5036 cmpi.l %d1,&0x3FD78000 # is |X| >= 2**(-40)?
5045 #--THIS IS THE USUAL CASE, |X| <= 15 PI.
5046 #--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
5051 lea PITBL+0x200(%pc),%a1 # TABLE OF N*PI/2, N = -32,...,32
5061 fsub.x (%a1)+,%fp0 # X-Y1
5062 fsub.s (%a1),%fp0 # fp0 = R = (X-Y1)-Y2
5065 #--continuation from REDUCEX
5067 #--GET N+ADJN AND SEE IF SIN(R) OR COS(R) IS NEEDED
5074 #--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
5075 #--THEN WE RETURN SGN*SIN(R). SGN*SIN(R) IS COMPUTED BY
5076 #--R' + R'*S*(A1 + S(A2 + S(A3 + S(A4 + ... + SA7)))), WHERE
5077 #--R' = SGN*R, S=R*R. THIS CAN BE REWRITTEN AS
5078 #--R' + R'*S*( [A1+T(A3+T(A5+TA7))] + [S(A2+T(A4+TA6))])
5079 #--WHERE T=S*S.
5080 #--NOTE THAT A3 THROUGH A7 ARE STORED IN DOUBLE PRECISION
5081 #--WHILE A1 AND A2 ARE IN DOUBLE-EXTENDED FORMAT.
5083 fmovm.x &0x0c,-(%sp) # save fp2/fp3
5119 fmul.x %fp1,%fp0 # SIN(R')-R'
5124 fadd.x X(%a6),%fp0 # last inst - possible exception set
5127 #--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
5128 #--THEN WE RETURN SGN*COS(R). SGN*COS(R) IS COMPUTED BY
5129 #--SGN + S'*(B1 + S(B2 + S(B3 + S(B4 + ... + SB8)))), WHERE
5130 #--S=R*R AND S'=SGN*S. THIS CAN BE REWRITTEN AS
5131 #--SGN + S'*([B1+T(B3+T(B5+TB7))] + [S(B2+T(B4+T(B6+TB8)))])
5132 #--WHERE T=S*S.
5133 #--NOTE THAT B4 THROUGH B8 ARE STORED IN DOUBLE PRECISION
5134 #--WHILE B2 AND B3 ARE IN DOUBLE-EXTENDED FORMAT, B1 IS -1/2
5135 #--AND IS THEREFORE STORED AS SINGLE PRECISION.
5137 fmovm.x &0x0c,-(%sp) # save fp2/fp3
5186 fadd.s POSNEG1(%a6),%fp0 # last inst - possible exception set
5192 #--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION.
5193 #--IF |X| < 2**(-40), RETURN X OR 1.
5203 # here, the operation may underflow iff the precision is sgl or dbl.
5210 fmov.x X(%a6),%fp0 # last inst - possible exception set
5216 fadd.s &0x80800000,%fp0 # last inst - possible exception set
5221 #--SIN(X) = X FOR DENORMALIZED X
5227 #--COS(X) = 1 FOR DENORMALIZED X
5236 #--SET ADJN TO 4
5246 cmp.l %d1,&0x3FD78000 # |X| >= 2**(-40)?
5256 #--THIS IS THE USUAL CASE, |X| <= 15 PI.
5257 #--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
5263 lea PITBL+0x200(%pc),%a1 # TABLE OF N*PI/2, N = -32,...,32
5271 fsub.x (%a1)+,%fp0 # X-Y1
5272 fsub.s (%a1),%fp0 # FP0 IS R = (X-Y1)-Y2
5275 #--continuation point from REDUCEX
5283 #--REGISTERS SAVED SO FAR: D0, A0, FP2.
5284 fmovm.x &0x04,-(%sp) # save fp2
5293 mov.l %d2,-(%sp)
5357 #--REGISTERS SAVED SO FAR: FP2.
5358 fmovm.x &0x04,-(%sp) # save fp2
5452 #--SIN AND COS OF X FOR DENORMALIZED X
5454 mov.l %d0,-(%sp) # save d0
5462 #--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW.
5463 #--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING
5464 #--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE.
5466 fmovm.x &0x3c,-(%sp) # save {fp2-fp5}
5467 mov.l %d2,-(%sp) # save d2
5470 #--If compact form of abs(arg) in d0=$7ffeffff, argument is so large that
5471 #--there is a danger of unwanted overflow in first LOOP iteration. In this
5472 #--case, reduce argument by one remainder step to make subsequent reduction
5473 #--safe.
5482 # create low half of 2**16383*PI/2 at FP_SCR1
5499 #--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4.
5500 #--integer quotient will be stored in N
5501 #--Intermeditate remainder is 66-bit long; (R,r) in (FP0,FP1)
5503 fmov.x %fp0,INARG(%a6) # +-2**K * F, 1 <= F < 2
5511 sub.l &27,%d1 # d0 = L := K-27
5519 #--FIND THE REMAINDER OF (R,r) W.R.T. 2**L * (PI/2). L IS SO CHOSEN
5520 #--THAT INT( X * (2/PI) / 2**(L) ) < 2**29.
5522 #--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63),
5523 #--2**L * (PIby2_1), 2**L * (PIby2_2)
5526 sub.l %d1,%d2 # BIASED EXP OF 2**(-L)*(2/PI)
5530 mov.w %d2,FP_SCR0_EX(%a6) # FP_SCR0 = 2**(-L)*(2/PI)
5533 fmul.x FP_SCR0(%a6),%fp2 # fp2 = X * 2**(-L)*(2/PI)
5535 #--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN
5536 #--FLOATING POINT FORMAT, THE TWO FMOVE'S FMOVE.L FP <--> N
5537 #--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT
5538 #--(SIGN(INARG)*2**63 + FP2) - SIGN(INARG)*2**63 WILL GIVE
5539 #--US THE DESIRED VALUE IN FLOATING POINT.
5549 #--CREATING 2**(L)*Piby2_1 and 2**(L)*Piby2_2
5564 #--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and
5565 #--P2 = 2**(L) * Piby2_2
5572 #--we want P+p = W+w but |p| <= half ulp of P
5573 #--Then, we need to compute A := R-P and a := r-p
5575 fsub.x %fp3,%fp4 # fp4 = W-P
5577 fsub.x %fp3,%fp0 # fp0 = A := R - P
5578 fadd.x %fp5,%fp4 # fp4 = p = (W-P)+w
5581 fsub.x %fp4,%fp1 # fp1 = a := r - p
5583 #--Now we need to normalize (A,a) to "new (R,r)" where R+r = A+a but
5584 #--|r| <= half ulp of R.
5586 #--No need to calculate r if this is the last loop
5590 #--Need to calculate r
5591 fsub.x %fp0,%fp3 # fp3 = A-R
5592 fadd.x %fp3,%fp1 # fp1 = r := (A-R)+a
5598 fmovm.x (%sp)+,&0x3c # restore {fp2-fp5}
5611 # a0 = pointer to extended precision input #
5612 # d0 = round precision,mode #
5620 # rounded to double precision. The result is provably monotonic #
5621 # in double precision. #
5625 # 1. If |X| >= 15Pi or |X| < 2**(-40), go to 6. #
5638 # 4. (k is odd) Tan(X) = -cot(r). Since tan(r) is approximated by #
5642 # -Cot(r) = -V/U. Exit. #
5646 # 7. (|X|<2**(-40)) Tan(X) = X. Exit. #
5681 #--N*PI/2, -32 <= N <= 32, IN A LEADING TERM IN EXT. AND TRAILING
5682 #--TERM IN SGL. NOTE THAT PI IS 64-BIT LONG, THUS N*PI/2 IS AT
5683 #--MOST 69 BITS LONG.
5766 cmp.l %d1,&0x3FD78000 # |X| >= 2**(-40)?
5775 #--THIS IS THE USUAL CASE, |X| <= 15 PI.
5776 #--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
5780 lea.l PITBL+0x200(%pc),%a1 # TABLE OF N*PI/2, N = -32,...,32
5787 fsub.x (%a1)+,%fp0 # X-Y1
5789 fsub.s (%a1),%fp0 # FP0 IS R = (X-Y1)-Y2
5795 fmovm.x &0x0c,-(%sp) # save fp2,fp3
5833 fdiv.x %fp1,%fp0 # last inst - possible exception set
5868 fmov.x %fp1,-(%sp)
5872 fdiv.x (%sp)+,%fp0 # last inst - possible exception set
5876 #--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION.
5877 #--IF |X| < 2**(-40), RETURN X OR 1.
5882 fmov.x %fp0,-(%sp)
5885 fmov.x (%sp)+,%fp0 # last inst - posibble exception set
5889 #--TAN(X) = X FOR DENORMALIZED X
5893 #--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW.
5894 #--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING
5895 #--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE.
5897 fmovm.x &0x3c,-(%sp) # save {fp2-fp5}
5898 mov.l %d2,-(%sp) # save d2
5901 #--If compact form of abs(arg) in d0=$7ffeffff, argument is so large that
5902 #--there is a danger of unwanted overflow in first LOOP iteration. In this
5903 #--case, reduce argument by one remainder step to make subsequent reduction
5904 #--safe.
5913 # create low half of 2**16383*PI/2 at FP_SCR1
5930 #--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4.
5931 #--integer quotient will be stored in N
5932 #--Intermeditate remainder is 66-bit long; (R,r) in (FP0,FP1)
5934 fmov.x %fp0,INARG(%a6) # +-2**K * F, 1 <= F < 2
5942 sub.l &27,%d1 # d0 = L := K-27
5950 #--FIND THE REMAINDER OF (R,r) W.R.T. 2**L * (PI/2). L IS SO CHOSEN
5951 #--THAT INT( X * (2/PI) / 2**(L) ) < 2**29.
5953 #--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63),
5954 #--2**L * (PIby2_1), 2**L * (PIby2_2)
5957 sub.l %d1,%d2 # BIASED EXP OF 2**(-L)*(2/PI)
5961 mov.w %d2,FP_SCR0_EX(%a6) # FP_SCR0 = 2**(-L)*(2/PI)
5964 fmul.x FP_SCR0(%a6),%fp2 # fp2 = X * 2**(-L)*(2/PI)
5966 #--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN
5967 #--FLOATING POINT FORMAT, THE TWO FMOVE'S FMOVE.L FP <--> N
5968 #--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT
5969 #--(SIGN(INARG)*2**63 + FP2) - SIGN(INARG)*2**63 WILL GIVE
5970 #--US THE DESIRED VALUE IN FLOATING POINT.
5980 #--CREATING 2**(L)*Piby2_1 and 2**(L)*Piby2_2
5995 #--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and
5996 #--P2 = 2**(L) * Piby2_2
6003 #--we want P+p = W+w but |p| <= half ulp of P
6004 #--Then, we need to compute A := R-P and a := r-p
6006 fsub.x %fp3,%fp4 # fp4 = W-P
6008 fsub.x %fp3,%fp0 # fp0 = A := R - P
6009 fadd.x %fp5,%fp4 # fp4 = p = (W-P)+w
6012 fsub.x %fp4,%fp1 # fp1 = a := r - p
6014 #--Now we need to normalize (A,a) to "new (R,r)" where R+r = A+a but
6015 #--|r| <= half ulp of R.
6017 #--No need to calculate r if this is the last loop
6021 #--Need to calculate r
6022 fsub.x %fp0,%fp3 # fp3 = A-R
6023 fadd.x %fp3,%fp1 # fp1 = r := (A-R)+a
6029 fmovm.x (%sp)+,&0x3c # restore {fp2-fp5}
6041 # a0 = pointer to extended precision input #
6042 # d0 = round precision,mode #
6050 # rounded to double precision. The result is provably monotonic #
6051 # in double precision. #
6057 # Note that k = -4, -3,..., or 3. #
6059 # significant bits of X with a bit-1 attached at the 6-th #
6060 # bit position. Define u to be u = (X-F) / (1 + X*F). #
6071 # Step 7. Define X' = -1/X. Approximate arctan(X') by an odd #
6240 #--ENTRY POINT FOR ATAN(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S
6258 #--THE MOST LIKELY CASE, |X| IN [1/16, 16). WE USE TABLE TECHNIQUE
6259 #--THE IDEA IS ATAN(X) = ATAN(F) + ATAN( [X-F] / [1+XF] ).
6260 #--SO IF F IS CHOSEN TO BE CLOSE TO X AND ATAN(F) IS STORED IN
6261 #--A TABLE, ALL WE NEED IS TO APPROXIMATE ATAN(U) WHERE
6262 #--U = (X-F)/(1+XF) IS SMALL (REMEMBER F IS CLOSE TO X). IT IS
6263 #--TRUE THAT A DIVIDE IS NOW NEEDED, BUT THE APPROXIMATION FOR
6264 #--ATAN(U) IS A VERY SHORT POLYNOMIAL AND THE INDEXING TO
6265 #--FETCH F AND SAVING OF REGISTERS CAN BE ALL HIDED UNDER THE
6266 #--DIVIDE. IN THE END THIS METHOD IS MUCH FASTER THAN A TRADITIONAL
6267 #--ONE. NOTE ALSO THAT THE TRADITIONAL SCHEME THAT APPROXIMATE
6268 #--ATAN(X) DIRECTLY WILL NEED TO USE A RATIONAL APPROXIMATION
6269 #--(DIVISION NEEDED) ANYWAY BECAUSE A POLYNOMIAL APPROXIMATION
6270 #--WILL INVOLVE A VERY LONG POLYNOMIAL.
6272 #--NOW WE SEE X AS +-2^K * 1.BBBBBBB....B <- 1. + 63 BITS
6273 #--WE CHOSE F TO BE +-2^K * 1.BBBB1
6274 #--THAT IS IT MATCHES THE EXPONENT AND FIRST 5 BITS OF X, THE
6275 #--SIXTH BITS IS SET TO BE 1. SINCE K = -4, -3, ..., 3, THERE
6276 #--ARE ONLY 8 TIMES 16 = 2^7 = 128 |F|'S. SINCE ATAN(-|F|) IS
6277 #-- -ATAN(|F|), WE NEED TO STORE ONLY ATAN(|F|).
6282 or.l &0x04000000,XFRAC(%a6) # SET 6-TH BIT TO 1
6287 fsub.x X(%a6),%fp0 # FP0 IS X-F
6289 fdiv.x %fp1,%fp0 # FP0 IS U = (X-F)/(1+X*F)
6291 #--WHILE THE DIVISION IS TAKING ITS TIME, WE FETCH ATAN(|F|)
6292 #--CREATE ATAN(F) AND STORE IT IN ATANF, AND
6293 #--SAVE REGISTERS FP2.
6295 mov.l %d2,-(%sp) # SAVE d2 TEMPORARILY
6313 #--THAT'S ALL I HAVE TO DO FOR NOW,
6314 #--BUT ALAS, THE DIVIDE IS STILL CRANKING!
6316 #--U IN FP0, WE ARE NOW READY TO COMPUTE ATAN(U) AS
6317 #--U + A1*U*V*(A2 + V*(A3 + V)), V = U*U
6318 #--THE POLYNOMIAL MAY LOOK STRANGE, BUT IS NEVERTHELESS CORRECT.
6319 #--THE NATURAL FORM IS U + U*V*(A1 + V*(A2 + V*A3))
6320 #--WHAT WE HAVE HERE IS MERELY A1 = A3, A2 = A1/A3, A3 = A2/A3.
6321 #--THE REASON FOR THIS REARRANGEMENT IS TO MAKE THE INDEPENDENT
6322 #--PARTS A1*U*V AND (A2 + ... STUFF) MORE LOAD-BALANCED
6324 fmovm.x &0x04,-(%sp) # save fp2
6344 #--|X| IS IN d0 IN COMPACT FORM. FP1, d0 SAVED.
6345 #--FP0 IS X AND |X| <= 1/16 OR |X| >= 16.
6350 #--|X| <= 1/16
6351 #--IF |X| < 2^(-40), RETURN X AS ANSWER. OTHERWISE, APPROXIMATE
6352 #--ATAN(X) BY X + X*Y*(B1+Y*(B2+Y*(B3+Y*(B4+Y*(B5+Y*B6)))))
6353 #--WHICH IS X + X*Y*( [B1+Z*(B3+Z*B5)] + [Y*(B2+Z*(B4+Z*B6)] )
6354 #--WHERE Y = X*X, AND Z = Y*Y.
6359 #--COMPUTE POLYNOMIAL
6360 fmovm.x &0x0c,-(%sp) # save fp2/fp3
6396 #--|X| < 2^(-40), ATAN(X) = X
6400 fmov.x X(%a6),%fp0 # last inst - possible exception set
6405 #--IF |X| > 2^(100), RETURN SIGN(X)*(PI/2 - TINY). OTHERWISE,
6406 #--RETURN SIGN(X)*PI/2 + ATAN(-1/X).
6410 #--APPROXIMATE ATAN(-1/X) BY
6411 #--X'+X'*Y*(C1+Y*(C2+Y*(C3+Y*(C4+Y*C5)))), X' = -1/X, Y = X'*X'
6412 #--THIS CAN BE RE-WRITTEN AS
6413 #--X'+X'*Y*( [C1+Z*(C3+Z*C5)] + [Y*(C2+Z*C4)] ), Z = Y*Y.
6415 fmovm.x &0x0c,-(%sp) # save fp2/fp3
6417 fmov.s &0xBF800000,%fp1 # LOAD -1
6418 fdiv.x %fp0,%fp1 # FP1 IS -1/X
6420 #--DIVIDE IS STILL CRANKING
6465 #--RETURN SIGN(X)*(PIBY2 - TINY) = SIGN(X)*PIBY2 - SIGN(X)*TINY
6482 #--ENTRY POINT FOR ATAN(X) FOR DENORMALIZED ARGUMENT
6491 # a0 = pointer to extended precision input #
6492 # d0 = round precision,mode #
6500 # rounded to double precision. The result is provably monotonic #
6501 # in double precision. #
6509 # z := sqrt( [1-X][1+X] ) #
6540 #--THIS IS THE USUAL CASE, |X| < 1
6541 #--ASIN(X) = ATAN( X / SQRT( (1-X)(1+X) ) )
6545 fsub.x %fp0,%fp1 # 1-X
6546 fmovm.x &0x4,-(%sp) # {fp2}
6549 fmul.x %fp2,%fp1 # (1+X)(1-X)
6551 fsqrt.x %fp1 # SQRT([1-X][1+X])
6552 fdiv.x %fp1,%fp0 # X/SQRT([1-X][1+X])
6553 fmovm.x &0x01,-(%sp) # save X/SQRT(...)
6564 #--|X| = 1, ASIN(X) = +- PI/2.
6569 or.l &0x3F800000,%d1 # +-1 IN SGL FORMAT
6570 mov.l %d1,-(%sp) # push SIGN(X) IN SGL-FMT
6575 #--|X| < 2^(-40), ATAN(X) = X
6579 fmov.x (%a0),%fp0 # last inst - possible exception
6583 #--ASIN(X) = X FOR DENORMALIZED X
6592 # a0 = pointer to extended precision input #
6593 # d0 = round precision,mode #
6601 # rounded to double precision. The result is provably monotonic #
6602 # in double precision. #
6610 # z := (1-X) / (1+X) #
6633 #--THIS IS THE USUAL CASE, |X| < 1
6634 #--ACOS(X) = 2 * ATAN( SQRT( (1-X)/(1+X) ) )
6639 fneg.x %fp0 # -X
6640 fadd.s &0x3F800000,%fp0 # 1-X
6641 fdiv.x %fp1,%fp0 # (1-X)/(1+X)
6642 fsqrt.x %fp0 # SQRT((1-X)/(1+X))
6643 mov.l %d0,-(%sp) # save original users fpcr
6645 fmovm.x &0x01,-(%sp) # save SQRT(...) to stack
6647 bsr satan # ATAN(SQRT([1-X]/[1+X]))
6659 #--|X| = 1, ACOS(X) = 0 OR PI
6663 #--X = -1
6675 #--ACOS(X) = PI/2 FOR DENORMALIZED X
6688 # a0 = pointer to extended precision input #
6689 # d0 = round precision,mode #
6692 # fp0 = exp(X) or exp(X)-1 #
6697 # rounded to double precision. The result is provably monotonic #
6698 # in double precision. #
6703 # ------ #
6706 # Step 2. Return ans := ans + sign(X)*2^(-126). Exit. #
6707 # Notes: This will always generate one exception -- inexact. #
6711 # ----- #
6714 # 1.1 If |X| >= 2^(-65), go to Step 1.3. #
6718 # Notes: The usual case should take the branches 1.1 -> 1.3 -> 2.#
6719 # To avoid the use of floating-point comparisons, a #
6721 # 32-bit integer, the upper (more significant) 16 bits #
6733 # Step 2. Calculate N = round-to-nearest-int( X * 64/log2 ). #
6734 # 2.1 Set AdjFlag := 0 (indicates the branch 1.3 -> 2 #
6736 # 2.2 N := round-to-nearest-integer( X * 64/log2 ). #
6739 # 2.4 Calculate M = (N - J)/64; so N = 64M + J. #
6745 # N := round-to-nearest-integer(Z) #
6747 # constant := single-precision( 64/log 2 ). #
6749 # Using a single-precision constant avoids memory #
6750 # access. Another effect of using a single-precision #
6753 # Z = X*(64/log2)*(1+eps), |eps| <= 2^(-24). #
6757 # Step 3. Calculate X - N*log2/64. #
6759 # where L1 := single-precision(-log2/64). #
6761 # L2 := extended-precision(-log2/64 - L1).#
6763 # approximate the value -log2/64 to 88 bits of accuracy. #
6772 # N = rnd-to-int( X*64/log2 (1+eps) ), |eps|<=2^(-24) #
6774 # X*64/log2 - N = f - eps*X 64/log2 #
6775 # X - N*log2/64 = f*log2/64 - eps*X #
6780 # |X - N*log2/64| <= (0.5 + 16446/2^(18))*log2/64 #
6784 # Step 4. Approximate exp(R)-1 by a polynomial #
6788 # and A5 are single precision; A2 and A3 are double #
6789 # precision. #
6791 # |p - (exp(R)-1)| < 2^(-68.8) for all |R| <= 0.0062. #
6803 # 2^(J/64) to roughly 85 bits; T is in extended precision #
6804 # and t is in single precision. Note also that T is #
6806 # zero. The reason for such a special form is that T-1, #
6807 # T-2, and T-8 will all be exact --- a property that will #
6835 # Notes: For non-zero X, the inexact exception will always be #
6846 # (mimic 2.2 - 2.6) #
6847 # 8.2 N := round-to-integer( X * 64/log2 ) #
6849 # 8.4 K := (N-J)/64, M1 := truncate(K/2), M = K-M1, #
6855 # Notes: Refer to notes for 2.2 - 2.6. #
6865 # extended-precision numbers whose square over/underflow #
6870 # -------- #
6876 # precision prescribed by the user FPCR. #
6879 # ------- #
6886 # Notes: The usual case should take the branches 1.1 -> 1.3 -> 2.#
6888 # because EXPM1 is intended to evaluate exp(X)-1 #
6892 # Step 2. Calculate N = round-to-nearest-int( X * 64/log2 ). #
6893 # 2.1 N := round-to-nearest-integer( X * 64/log2 ). #
6896 # 2.3 Calculate M = (N - J)/64; so N = 64M + J. #
6900 # OnebySc := -2^(-M). #
6903 # Step 3. Calculate X - N*log2/64. #
6905 # where L1 := single-precision(-log2/64). #
6907 # L2 := extended-precision(-log2/64 - L1).#
6912 # Step 4. Approximate exp(R)-1 by a polynomial #
6916 # and A6 are single precision; A2, A3 and A4 are double #
6917 # precision. #
6919 # |p - (exp(R)-1)| < |R| * 2^(-72.7) #
6931 # 2^(J/64) to roughly 85 bits; T is in extended precision #
6932 # and t is in single precision. Note also that T is #
6934 # zero. The reason for such a special form is that T-1, #
6935 # T-2, and T-8 will all be exact --- a property that will #
6937 # in p is no bigger than 2^(-67.7) compared to the final #
6940 # Step 6. Reconstruction of exp(X)-1 #
6941 # exp(X)-1 = 2^M * ( 2^(J/64) + p - 2^(-M) ). #
6944 # 6.3 If M >= -3, go to 6.5. #
6952 # Step 7. exp(X)-1 for |X| < 1/4. #
6953 # 7.1 If |X| >= 2^(-65), go to Step 9. #
6956 # Step 8. Calculate exp(X)-1, |X| < 2^(-65). #
6957 # 8.1 If |X| < 2^(-16312), goto 8.3 #
6958 # 8.2 Restore FPCR; return ans := X - 2^(-16382). #
6961 # 8.4 Restore FPCR; ans := ans - 2^(-16382). #
6963 # Notes: The idea is to return "X - tiny" under the user #
6964 # precision and rounding modes. To avoid unnecessary #
6966 # the best we can. For |X| >= 2^(-16312), the #
6970 # Step 9. Calculate exp(X)-1, |X| < 1/4, by a polynomial #
6974 # to B12 are single precision; B3 to B8 are double #
6975 # precision; and B2 is double extended. #
6977 # |p - (exp(X)-1)| < |X| 2^(-70.6) #
6989 # Step 10. Calculate exp(X)-1 for |X| >= 70 log 2. #
6990 # 10.1 If X >= 70log2 , exp(X) - 1 = exp(X) for all #
6992 # 10.2 If X <= -70log2, exp(X) - 1 = -1 for all practical #
6994 # ans := -1 #
6996 # Return ans := ans + 2^(-126). Exit. #
6997 # Notes: 10.2 will always create an inexact and return -1 + tiny #
6998 # in the user rounding precision and mode. #
7102 #--entry point for EXP(X), here X is finite, non-zero, and not NaN's
7104 #--Step 1.
7107 cmp.l %d1,&0x3FBE0000 # 2^(-65)
7112 #--The case |X| >= 2^(-65)
7119 #--Step 2.
7120 #--This is the normal branch: 2^(-65) <= |X| < 16380 log2.
7125 fmovm.x &0xc,-(%sp) # save fp2 {%fp2/%fp3}
7129 fmov.l %d1,%fp0 # convert to floating-format
7141 #--Step 3.
7142 #--fp1,fp2 saved on the stack. fp0 is N, fp1 is X,
7143 #--a0 points to 2^(J/64), D0 is biased expo. of 2^(M)
7145 fmul.s &0xBC317218,%fp0 # N * L1, L1 = lead(-log2/64)
7146 fmul.x L2(%pc),%fp2 # N * L2, L1+L2 = -log2/64
7150 #--Step 4.
7151 #--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL
7152 #-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*A5))))
7153 #--TO FULLY UTILIZE THE PIPELINE, WE COMPUTE S = R*R
7154 #--[R+R*S*(A2+S*A4)] + [S*(A1+S*(A3+S*A5))]
7182 fadd.x %fp2,%fp0 # fp0 is EXP(R) - 1
7184 #--Step 5
7185 #--final reconstruction process
7186 #--EXP(X) = 2^M * ( 2^(J/64) + 2^(J/64)*(EXP(R)-1) )
7188 fmul.x %fp1,%fp0 # 2^(J/64)*(Exp(R)-1)
7195 #--Step 6
7207 #--Step 7
7214 #--Step 8
7217 #--Steps 8.2 -- 8.6
7222 fmovm.x &0xc,-(%sp) # save fp2 {%fp2/%fp3}
7226 fmov.l %d1,%fp0 # convert to floating-format
7245 #--Step 9
7252 #--entry point for EXP(X), X is denormalized
7253 mov.l (%a0),-(%sp)
7255 ori.l &0x00800000,(%sp) # sign(X)*2^(-126)
7265 #--entry point for EXPM1(X), here X is finite, non-zero, non-NaN
7267 #--Step 1.
7268 #--Step 1.1
7276 #--Step 1.3
7277 #--The case |X| >= 1/4
7284 #--Step 2.
7285 #--This is the case: 1/4 <= |X| <= 70 log2.
7290 fmovm.x &0xc,-(%sp) # save fp2 {%fp2/%fp3}
7293 fmov.l %d1,%fp0 # convert to floating-format
7303 #--Step 3.
7304 #--fp1,fp2 saved on the stack. fp0 is N, fp1 is X,
7305 #--a0 points to 2^(J/64), D0 and a1 both contain M
7307 fmul.s &0xBC317218,%fp0 # N * L1, L1 = lead(-log2/64)
7308 fmul.x L2(%pc),%fp2 # N * L2, L1+L2 = -log2/64
7313 #--Step 4.
7314 #--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL
7315 #-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*(A5 + R*A6)))))
7316 #--TO FULLY UTILIZE THE PIPELINE, WE COMPUTE S = R*R
7317 #--[R*S*(A2+S*(A4+S*A6))] + [R+S*(A1+S*(A3+S*A5))]
7336 neg.w %d1 # D0 is -M
7338 add.w &0x3FFF,%d1 # biased expo. of 2^(-M)
7343 or.w &0x8000,%d1 # signed/expo. of -2^(-M)
7344 mov.w %d1,ONEBYSC(%a6) # OnebySc is -2^(-M)
7352 fadd.x %fp2,%fp0 # fp0 IS EXP(R)-1
7356 #--Step 5
7357 #--Compute 2^(J/64)*p
7359 fmul.x (%a1),%fp0 # 2^(J/64)*(Exp(R)-1)
7361 #--Step 6
7362 #--Step 6.1
7366 #--Step 6.2 M >= 64
7373 #--Step 6.3 M <= 63
7374 cmp.l %d1,&-3
7377 #--Step 6.4 M <= -4
7383 #--Step 6.5 -3 <= M <= 63
7390 #--Step 6.6
7396 #--Step 7 |X| < 1/4.
7397 cmp.l %d1,&0x3FBE0000 # 2^(-65)
7401 #--Step 8 |X| < 2^(-65)
7402 cmp.l %d1,&0x00330000 # 2^(-16312)
7404 #--Step 8.2
7405 mov.l &0x80010000,SC(%a6) # SC is -2^(-16382)
7415 #--Step 8.3
7428 #--Step 9 exp(X)-1 by a simple polynomial
7431 fmovm.x &0xc,-(%sp) # save fp2 {%fp2/%fp3}
7479 #--Step 10 |X| > 70 log2
7483 #--Step 10.2
7484 fmov.s &0xBF800000,%fp0 # fp0 is -1
7486 fadd.s &0x00800000,%fp0 # -1 + 2^(-126)
7491 #--entry point for EXPM1(X), here X is denormalized
7492 #--Step 0.
7498 # returned as an extended precision number in fp0. #
7502 # mantissa is converted to an extended precision number w/ #
7504 # the result is [1.0 - 2.0). #
7508 # a0 = pointer to extended precision input #
7531 neg.w %d0 # new exp = -(shft amt)
7541 bclr &0xe,%d0 # make it the new exp +-3fff
7556 # For denormalized numbers, shift the mantissa until the j-bit = 1,
7569 # a0 = pointer to extended precision input #
7570 # d0 = round precision,mode #
7578 # rounded to double precision. The result is provably monotonic #
7579 # in double precision. #
7599 # Y' := Y - 16381 log2 #
7623 #--THIS IS THE USUAL CASE, |X| < 16380 LOG2
7624 #--COSH(X) = (1/2) * ( EXP(X) + 1/EXP(X) )
7628 mov.l %d0,-(%sp)
7630 fmovm.x &0x01,-(%sp) # save |X| to stack
7650 fsub.d T1(%pc),%fp0 # (|X|-16381LOG2_LEAD)
7651 fsub.d T2(%pc),%fp0 # |X| - 16381 LOG2, ACCURATE
7653 mov.l %d0,-(%sp)
7655 fmovm.x &0x01,-(%sp) # save fp0 to stack
7670 #--COSH(X) = 1 FOR DENORMALIZED X
7683 # a0 = pointer to extended precision input #
7684 # d0 = round precision,mode #
7692 # rounded to double precision. The result is provably monotonic #
7693 # in double precision. #
7714 # Y' := Y - 16381 log2 #
7736 #--THIS IS THE USUAL CASE, |X| < 16380 LOG2
7737 #--Y = |X|, Z = EXPM1(Y), SINH(X) = SIGN(X)*(1/2)*( Z + Z/(1+Z) )
7741 movm.l &0x8040,-(%sp) # {a1/d0}
7742 fmovm.x &0x01,-(%sp) # save Y on stack
7752 fmov.x %fp0,-(%sp)
7758 mov.l %d1,-(%sp)
7762 fmul.s (%sp)+,%fp0 # last fp inst - possible exceptions set
7769 fsub.d T1(%pc),%fp0 # (|X|-16381LOG2_LEAD)
7770 mov.l &0,-(%sp)
7771 mov.l &0x80000000,-(%sp)
7775 mov.l %d1,-(%sp) # EXTENDED FMT
7776 fsub.d T2(%pc),%fp0 # |X| - 16381 LOG2, ACCURATE
7778 mov.l %d0,-(%sp)
7780 fmovm.x &0x01,-(%sp) # save fp0 on stack
7792 #--SINH(X) = X FOR DENORMALIZED X
7801 # a0 = pointer to extended precision input #
7802 # d0 = round precision,mode #
7810 # rounded to double precision. The result is provably monotonic #
7811 # in double precision. #
7816 # 1. If |X| >= (5/2) log2 or |X| <= 2**(-40), go to 3. #
7818 # 2. (2**(-40) < |X| < (5/2) log2) Calculate tanh(X) by #
7823 # 3. (|X| <= 2**(-40) or |X| >= (5/2) log2). If |X| < 1, #
7830 # tanh(X) = sgn - [ sgn*2/(1+z) ]. #
7833 # 6. (|X| >= 50 log2) Tanh(X) = +-1 (round to nearest). Thus, we #
7835 # sgn := sign(X), Tiny := 2**(-126), #
7836 # tanh(X) := sgn - sgn*Tiny. #
7839 # 7. (|X| < 2**(-40)). Tanh(X) = X. Exit. #
7859 cmp.l %d1, &0x3fd78000 # is |X| < 2^(-40)?
7864 #--THIS IS THE USUAL CASE
7865 #--Y = 2|X|, Z = EXPM1(Y), TANH(X) = SIGN(X) * Z / (Z+2).
7875 mov.l %d0,-(%sp)
7877 fmovm.x &0x1,-(%sp) # save Y on stack
7900 #-- (5/2) LOG2 < |X| < 50 LOG2,
7901 #--TANH(X) = 1 - (2/[EXP(2X)+1]). LET Y = 2|X|, SGN = SIGN(X),
7902 #--TANH(X) = SGN - SGN*2/[EXP(Y)+1].
7913 mov.l %d0,-(%sp)
7915 fmovm.x &0x01,-(%sp) # save Y on stack
7923 eor.l &0xC0000000,%d1 # -SIGN(X)*2
7924 fmov.s %d1,%fp1 # -SIGN(X)*2 IN SGL FMT
7925 fdiv.x %fp0,%fp1 # -SIGN(X)2 / [EXP(Y)+1 ]
7939 fmov.x X(%a6),%fp0 # last inst - possible exception set
7942 #---RETURN SGN(X) - SGN(X)EPS
7949 eor.l &0x80800000,%d1 # -SIGN(X)*EPS
7956 #--TANH(X) = X FOR DENORMALIZED X
7967 # a0 = pointer to extended precision input #
7968 # d0 = round precision,mode #
7976 # rounded to double precision. The result is provably monotonic #
7977 # in double precision. #
7981 # Step 1. If |X-1| < 1/16, approximate log(X) by an odd #
7982 # polynomial in u, where u = 2(X-1)/(X+1). Otherwise, #
7986 # seven significant bits of Y plus 2**(-7), i.e. #
7988 # of Y. Note that |Y-F| <= 2**(-7). #
7990 # Step 3. Define u = (Y-F)/F. Approximate log(1+u) by a #
8006 # approximates log(1+u), u = (Y-F)/F. #
8011 # 1/F are also tabulated so that the division in (Y-F)/F #
8015 # the value Y-F has to be calculated carefully when #
8212 #--ENTRY POINT FOR LOG(X) FOR X FINITE, NON-ZERO, NOT NAN'S
8218 #--FPCR SAVED AND CLEARED, INPUT IS 2^(ADJK)*FP0, FP0 CONTAINS
8219 #--A FINITE, NON-ZERO, NORMALIZED NUMBER.
8237 #--THIS SHOULD BE THE USUAL CASE, X NOT VERY CLOSE TO 1
8239 #--X = 2^(K) * Y, 1 <= Y < 2. THUS, Y = 1.XXXXXXXX....XX IN BINARY.
8240 #--WE DEFINE F = 1.XXXXXX1, I.E. FIRST 7 BITS OF Y AND ATTACH A 1.
8241 #--THE IDEA IS THAT LOG(X) = K*LOG2 + LOG(Y)
8242 #-- = K*LOG2 + LOG(F) + LOG(1 + (Y-F)/F).
8243 #--NOTE THAT U = (Y-F)/F IS VERY SMALL AND THUS APPROXIMATING
8244 #--LOG(1+U) CAN BE VERY EFFICIENT.
8245 #--ALSO NOTE THAT THE VALUE 1/F IS STORED IN A TABLE SO THAT NO
8246 #--DIVISION IS NEEDED TO CALCULATE (Y-F)/F.
8248 #--GET K, Y, F, AND ADDRESS OF 1/F.
8254 fmov.l %d1,%fp1 # CONVERT K TO FLOATING-POINT FORMAT
8256 #--WHILE THE CONVERSION IS GOING ON, WE GET F AND ADDRESS OF 1/F
8257 mov.l &0x3FFF0000,X(%a6) # X IS NOW Y, I.E. 2^(-K)*X
8271 fsub.x F(%a6),%fp0 # Y-F
8272 fmovm.x &0xc,-(%sp) # SAVE FP2-3 WHILE FP0 IS NOT READY
8273 #--SUMMARY: FP0 IS Y-F, A0 IS ADDRESS OF 1/F, FP1 IS K
8274 #--REGISTERS SAVED: FPCR, FP1, FP2
8277 #--AN RE-ENTRY POINT FOR LOGNP1
8278 fmul.x (%a0),%fp0 # FP0 IS U = (Y-F)/F
8284 #--LOG(1+U) IS APPROXIMATED BY
8285 #--U + V*(A1+U*(A2+U*(A3+U*(A4+U*(A5+U*A6))))) WHICH IS
8286 #--[U + V*(A1+V*(A3+V*A5))] + [U*V*(A2+V*(A4+V*A6))]
8311 fmovm.x (%sp)+,&0x30 # RESTORE FP2-3
8327 #--REGISTERS SAVED: FPCR, FP1. FP0 CONTAINS THE INPUT.
8329 fsub.s one(%pc),%fp1 # FP1 IS X-1
8331 fadd.x %fp1,%fp1 # FP1 IS 2(X-1)
8332 #--LOG(X) = LOG(1+U/2)-LOG(1-U/2) WHICH IS AN ODD POLYNOMIAL
8333 #--IN U, U = 2(X-1)/(X+1) = FP1/FP0
8336 #--THIS IS AN RE-ENTRY POINT FOR LOGNP1
8338 fmovm.x &0xc,-(%sp) # SAVE FP2-3
8339 #--REGISTERS SAVED ARE NOW FPCR,FP1,FP2,FP3
8340 #--LET V=U*U, W=V*V, CALCULATE
8341 #--U + U*V*(B1 + V*(B2 + V*(B3 + V*(B4 + V*B5)))) BY
8342 #--U + U*V*( [B1 + W*(B3 + W*B5)] + [V*(B2 + W*B4)] )
8366 fmovm.x (%sp)+,&0x30 # FP2-3 RESTORED
8374 #--REGISTERS SAVED FPCR. LOG(-VE) IS INVALID
8380 #--ENTRY POINT FOR LOG(X) FOR DENORMALIZED INPUT
8382 mov.l &-100,ADJK(%a6) # INPUT = 2^(ADJK) * FP0
8384 #----normalize the input value by left shifting k bits (k to be determined
8385 #----below), adjusting exponent and storing -k to ADJK
8386 #----the value TWOTO100 is no longer needed.
8387 #----Note that this code assumes the denormalized input is NON-ZERO.
8389 movm.l &0x3f00,-(%sp) # save some registers {d2-d7}
8413 movm.l (%sp)+,&0xfc # restore registers {d2-d7}
8435 movm.l (%sp)+,&0xfc # restore registers {d2-d7}
8440 #--ENTRY POINT FOR LOG(1+X) FOR X FINITE, NON-ZERO, NOT NAN'S
8460 ble.w LP1NEG0 # LOG OF ZERO OR -VE
8465 #--IF 1+Z > 3/2 OR 1+Z < 1/2, THEN X, WHICH IS ROUNDING 1+Z,
8466 #--CONTAINS AT LEAST 63 BITS OF INFORMATION OF Z. IN THAT CASE,
8467 #--SIMPLY INVOKE LOG(X) FOR LOG(1+Z).
8470 #--NEXT SEE IF EXP(-1/16) < X < EXP(1/16)
8477 #--EXP(-1/16) < X < EXP(1/16). LOG(1+Z) = LOG(1+U/2) - LOG(1-U/2)
8478 #--WHERE U = 2Z/(2+Z) = 2Z/(1+X).
8481 #--U = FP1/FP0
8485 #--HERE WE USE THE USUAL TABLE DRIVEN APPROACH. CARE HAS TO BE
8486 #--TAKEN BECAUSE 1+Z CAN HAVE 67 BITS OF INFORMATION AND WE MUST
8487 #--PRESERVE ALL THE INFORMATION. BECAUSE 1+Z IS IN [1/2,3/2],
8488 #--THERE ARE ONLY TWO CASES.
8489 #--CASE 1: 1+Z < 1, THEN K = -1 AND Y-F = (2-F) + 2Z
8490 #--CASE 2: 1+Z > 1, THEN K = 0 AND Y-F = (1-F) + Z
8491 #--ON RETURNING TO LP1CONT1, WE MUST HAVE K IN FP1, ADDRESS OF
8492 #--(1/F) IN A0, Y-F IN FP0, AND FP2 SAVED.
8504 fsub.x F(%a6),%fp0 # 2-F
8511 fmovm.x &0xc,-(%sp) # SAVE FP2 {%fp2/%fp3}
8512 fadd.x %fp1,%fp0 # FP0 IS Y-F = (2-F)+2Z
8515 fmov.s negone(%pc),%fp1 # FP1 IS K = -1
8522 fsub.x F(%a6),%fp0 # 1-F
8528 fadd.x %fp1,%fp0 # FP0 IS Y-F
8529 fmovm.x &0xc,-(%sp) # FP2 SAVED {%fp2/%fp3}
8536 #--FPCR SAVED. D0 IS X IN COMPACT FORM.
8552 #--ENTRY POINT FOR LOG(1+Z) FOR DENORMALIZED INPUT
8562 # a0 = pointer to extended precision input #
8563 # d0 = round precision,mode #
8571 # rounded to double precision. The result is provably monotonic #
8572 # in double precision. #
8582 # z := 2y/(1-y) #
8589 # divide-by-zero by #
8607 #--THIS IS THE USUAL CASE, |X| < 1
8608 #--Y = |X|, Z = 2Y/(1-Y), ATANH(X) = SIGN(X) * (1/2) * LOG1P(Z).
8612 fneg.x %fp1 # -Y
8614 fadd.s &0x3F800000,%fp1 # 1-Y
8615 fdiv.x %fp1,%fp0 # 2Y/(1-Y)
8618 or.l &0x3F000000,%d1 # SIGN(X)*HALF
8619 mov.l %d1,-(%sp)
8621 mov.l %d0,-(%sp) # save rnd prec,mode
8623 fmovm.x &0x01,-(%sp) # save Z on stack
8641 #--ATANH(X) = X FOR DENORMALIZED X
8646 # slog10(): computes the base-10 logarithm of a normalized input #
8647 # slog10d(): computes the base-10 logarithm of a denormalized input #
8648 # slog2(): computes the base-2 logarithm of a normalized input #
8649 # slog2d(): computes the base-2 logarithm of a denormalized input #
8652 # a0 = pointer to extended precision input #
8653 # d0 = round precision,mode #
8661 # rounded to double precision. The result is provably monotonic #
8662 # in double precision. #
8670 # Notes: Default means round-to-nearest mode, no floating-point #
8671 # traps, and precision control = double extended. #
8684 # Notes: Default means round-to-nearest mode, no floating-point #
8685 # traps, and precision control = double extended. #
8697 # Notes: Default means round-to-nearest mode, no floating-point #
8698 # traps, and precision control = double extended. #
8711 # Notes: Default means round-to-nearest mode, no floating-point #
8712 # traps, and precision control = double extended. #
8720 # 2.3 Return ans := convert-to-double-extended(k). #
8737 #--entry point for Log10(X), X is normalized
8745 mov.l %d0,-(%sp)
8753 #--entry point for Log10(X), X is denormalized
8757 mov.l %d0,-(%sp)
8765 #--entry point for Log2(X), X is normalized
8777 #--X = 2^k.
8787 mov.l %d0,-(%sp)
8798 #--entry point for Log2(X), X is denormalized
8802 mov.l %d0,-(%sp)
8816 # a0 = pointer to extended precision input #
8817 # d0 = round precision,mode #
8825 # rounded to double precision. The result is provably monotonic #
8826 # in double precision. #
8833 # 2. If |X| < 2**(-70), go to ExpSm. #
8846 # 2. If |X| < 2**(-70), go to ExpSm. #
8849 # N := round-to-int(y). Decompose N as #
8853 # r := ((X - N*L1)-N*L2) * L10 #
8984 #--ENTRY POINT FOR 2**(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S
8993 cmp.l %d1,&0x3FB98000 # |X| >= 2**(-70)?
9003 #--USUAL CASE, 2^(-70) <= |X| <= 16480
9007 fmov.l %fp1,INT(%a6) # N = ROUND-TO-INT(64 X)
9008 mov.l %d2,-(%sp)
9010 fmov.l INT(%a6),%fp1 # N --> FLOATING FMT
9022 #--SUMMARY: a1 IS ADDRESS FOR THE LEADING PORTION OF 2^(J/64),
9023 #--D0 IS M WHERE N = 64(M+M') + J. NOTE THAT |M| <= 16140 BY DESIGN.
9024 #--ADJFACT = 2^(M').
9025 #--REGISTERS SAVED SO FAR ARE (IN ORDER) FPCR, D0, FP1, a1, AND FP2.
9027 fmovm.x &0x0c,-(%sp) # save fp2/fp3
9035 fsub.x %fp1,%fp0 # X - (1/64)*INT(64 X)
9047 #--FPCR, D0 SAVED
9051 #--|X| IS SMALL, RETURN 1 + X
9058 #--|X| IS LARGE, GENERATE OVERFLOW IF X > 0; ELSE GENERATE UNDERFLOW
9059 #--REGISTERS SAVE SO FAR ARE FPCR AND D0
9071 #--ENTRY POINT FOR 2**(X) FOR DENORMALIZED ARGUMENT
9073 fmov.l %d0,%fpcr # set user's rounding mode/precision
9081 #--ENTRY POINT FOR 10**(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S
9090 cmp.l %d1,&0x3FB98000 # |X| >= 2**(-70)?
9100 #--USUAL CASE, 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10
9105 mov.l %d2,-(%sp)
9107 fmov.l INT(%a6),%fp1 # N --> FLOATING FMT
9119 #--SUMMARY: a1 IS ADDRESS FOR THE LEADING PORTION OF 2^(J/64),
9120 #--D0 IS M WHERE N = 64(M+M') + J. NOTE THAT |M| <= 16140 BY DESIGN.
9121 #--ADJFACT = 2^(M').
9122 #--REGISTERS SAVED SO FAR ARE (IN ORDER) FPCR, D0, FP1, a1, AND FP2.
9123 fmovm.x &0x0c,-(%sp) # save fp2/fp3
9134 fsub.x %fp1,%fp0 # X - N L_LEAD
9137 fsub.x %fp2,%fp0 # X - N L_TRAIL
9148 #--FPCR, FP2, FP3 ARE SAVED IN ORDER AS SHOWN.
9149 #--ADJFACT CONTAINS 2**(M'), FACT1 + FACT2 = 2**(M) * 2**(J/64).
9150 #--FP0 IS R. THE FOLLOWING CODE COMPUTES
9151 #-- 2**(M'+M) * 2**(J/64) * EXP(R)
9173 fadd.x %fp2,%fp0 # FP0 IS EXP(R) - 1
9177 #--FINAL RECONSTRUCTION PROCESS
9178 #--EXP(X) = 2^M*2^(J/64) + 2^M*2^(J/64)*(EXP(R)-1) - (1 OR 0)
9195 #--ENTRY POINT FOR 10**(X) FOR DENORMALIZED ARGUMENT
9197 fmov.l %d0,%fpcr # set user's rounding mode/precision
9210 # a0 = pointer to double-extended source operand X #
9211 # a1 = pointer to double-extended destination operand Y #
9222 mov.l %d0,-(%sp) # store off ctrl bits for now
9251 mov.l %d0,-(%sp) # save src for now
9264 cmpi.w %d0,&-0x3fff # is the shft amt really low?
9275 subi.l &-0x3fff,%d0 # how many should we shift?
9280 clr.l -(%sp) # insert zero low mantissa
9281 mov.l %d1,-(%sp) # insert new high mantissa
9282 clr.l -(%sp) # make zero exponent
9287 mov.l %d1,-(%sp) # insert new low mantissa
9288 clr.l -(%sp) # insert zero high mantissa
9289 clr.l -(%sp) # make zero exponent
9301 clr.l -(%sp) # insert new exponent
9302 mov.l &0x80000000,-(%sp) # insert new high mantissa
9303 mov.l %d0,-(%sp) # insert new lo mantissa
9345 # a0 = pointer to extended precision input X #
9346 # a1 = pointer to extended precision input Y #
9347 # d0 = round precision,mode #
9362 # Step 2. Set L := expo(X)-expo(Y), k := 0, Q := 0. #
9366 # R := 2^(-L)X, j := L. #
9371 # 3.2 If R > Y, then { R := R - Y, Q := Q + 1} #
9373 # 3.4 k := k + 1, j := j - 1, Q := 2Q, R := 2R. Go to #
9376 # Step 4. At this point, R = X - QY = MOD(X,Y). Set #
9386 # then { Q := Q + 1, signX := -signX }. #
9390 # Step 7. If Last_Subtract = true, R := R - Y. #
9394 # Step 9. At this point, R = 2^(-j)*X - Q Y = Y. Thus, #
9421 mov.l %d0,-(%sp) # save ctrl bits
9428 mov.l %d0,-(%sp) # save ctrl bits
9433 movm.l &0x3f00,-(%sp) # save data registers
9526 mov.l %d0,-(%sp) # save biased exp(X)
9527 sub.l %d3,%d0 # L := expo(X)-expo(Y)
9529 clr.l %d6 # D6 := carry <- 0
9544 #..At this point R = 2^(-L)X; Q = 0; k = 0; and k+j = L
9564 #..and Y < (D1,D2) < 2Y. Either way, perform R - Y
9565 sub.l %d5,%d2 # lo(R) - lo(Y)
9566 subx.l %d4,%d1 # hi(R) - hi(Y)
9571 #..At this point, Carry=0, R < Y. R = 2^(k-L)X - QY; k+j = L; j >= 0.
9580 subq.l &1,%d0 # j := j - 1
9581 #..At this point, R=(Carry,D1,D2) = 2^(k-L)X - QY, j+k=L, j >= 0, R < 2Y.
9586 #..k = L, j = 0, Carry = 0, R = (D1,D2) = X - QY, R < Y.
9696 movm.l (%sp)+,&0xfc # {%d2-%d7}
9717 #..R = 2^(-j)X - Q Y = Y, thus R = 0 and quotient = 2^j (Q+1)
9740 #..Q is odd, Q := Q + 1, signX := -signX
9757 # a0 = pointer to extended precision operand #
9764 # Simply test the exponent, j-bit, and mantissa values to #
9831 # a0 = pointer to extended precision source operand. #
9858 # dz is disabled. return a -INF.
9859 fmov.s &0xff800000,%fp0 # return -INF
9867 fmov.s &0xbf800000,%fp1 # load -1
9868 fdiv.s &0x00000000,%fp1 # -1 / 0
9931 fmovm.x &0x04,-(%sp) # save fp2
9958 # a0 = pointer to extended precision source operand #
10015 # a0 = pointer to extended precision source operand #
10038 # result would be inexact for the given precision. make a copy of the
10044 movm.l &0xc080,-(%sp) # save d0-d1/a0
10046 movm.l (%sp)+,&0x0103 # restore d0-d1/a0
10048 cmpi.b %d1,&0x40 # is precision sgl?
10143 # t_minx2(): Handle inexact 060FPLSP exception for "-" results. #
10204 # a0 = pointer to extended precision input operand #
10248 # dst_qnan --- force result when destination is a NaN
10263 # src_qnan --- force result when source is a NaN
10288 fmov.l %fpcr,-(%sp) # save fpcr
10297 fmov.l %fpcr,-(%sp) # save fpcr
10312 fmov.l %fpcr,-(%sp) # save fpcr
10321 fmov.l %fpcr,-(%sp) # save fpcr
10336 fmov.l %fpcr,-(%sp) # save fpcr
10345 fmov.l %fpcr,-(%sp) # save fpcr
10360 fmov.l %fpcr,-(%sp) # save fpcr
10369 fmov.l %fpcr,-(%sp) # save fpcr
10479 fmov.s &0x80000000,%fp0 # load -0
10514 fmov.s &0xff800000,%fp0 # load -INF
10582 fmov.s &0xbf800000,%fp0 # load -1
10612 fmov.x mpiby2(%pc),%fp0 # load -pi/2
10710 mov.l %d0,-(%sp)
10817 # norm(): normalize the mantissa of an extended precision input. the #
10827 # a0 = pointer fp extended precision operand to normalize #
10837 mov.l %d2, -(%sp) # create some temp regs
10838 mov.l %d3, -(%sp)
10879 # unnorm_fix(): - changes an UNNORM to one of NORM, DENORM, or ZERO #
10880 # - returns corresponding optype tag #
10886 # norm() - normalize the mantissa #
10889 # a0 = pointer to unnormalized extended precision number #
10892 # d0 = optype tag - is corrected to one of NORM, DENORM, or ZERO #