diff options
Diffstat (limited to 'newlib/libm/machine/spu/headers/divd2.h')
| -rw-r--r-- | newlib/libm/machine/spu/headers/divd2.h | 279 |
1 files changed, 142 insertions, 137 deletions
diff --git a/newlib/libm/machine/spu/headers/divd2.h b/newlib/libm/machine/spu/headers/divd2.h index 7bcf366eb..005194b86 100644 --- a/newlib/libm/machine/spu/headers/divd2.h +++ b/newlib/libm/machine/spu/headers/divd2.h @@ -51,65 +51,115 @@ * * DESCRIPTION * _divd2 divides the vector dividend a by the vector divisor b and - * returns the resulting vector quotient. Maximum error 0.5 ULPS for - * normalized results, 1ulp for denorm results, over entire double - * range including denorms, compared to true result in round-to-nearest - * rounding mode. Handles Inf or NaN operands and results correctly. + * returns the resulting vector quotient. Maximum error about 0.5 ulp + * over entire double range including denorms, compared to true result + * in round-to-nearest rounding mode. Handles Inf or NaN operands and + * results correctly. */ -static __inline vector double _divd2(vector double a, vector double b) +static __inline vector double _divd2(vector double a_in, vector double b_in) { + /* Variables */ + vec_int4 exp, exp_bias; + vec_uint4 no_underflow, overflow; + vec_float4 mant_bf, inv_bf; + vec_ullong2 exp_a, exp_b; + vec_ullong2 a_nan, a_zero, a_inf, a_denorm; + vec_ullong2 b_nan, b_zero, b_inf, b_denorm; + vec_ullong2 nan; + vec_double2 a, b; + vec_double2 mant_a, mant_b, inv_b, q0, q1, q2, mult; + /* Constants */ + vec_float4 onef = spu_splats(1.0f); + vec_ullong2 exp_mask = spu_splats(0x7FF0000000000000ULL); + vec_double2 one = spu_splats(1.0); - /* Variables - */ - vec_float4 inv_bf, mant_bf; - vec_double2 mant_a, mant_b, inv_b, q0, q1, q2, mult; - vec_int4 exp, tmp; - vec_uint4 exp_a, exp_b, exp_q1, overflow, nounderflow, normal, utmp, - sign_a, sign_b, a_frac, b_frac, a_frac_0, b_frac_0, a_exp_0, b_exp_0, - a_exp_ones, b_exp_ones, a_nan, b_nan, a_inf, b_inf, a_zero, b_zero, - res_nan, sign_res; +#ifdef __SPU_EDP__ + vec_double2 denorm_scale = (vec_double2)spu_splats(0x4330000000000000ULL); - /* Constants + /* Identify all possible special values that must be accomodated including: + * +-0, +-infinity, +-denorm, and NaNs. */ - vec_float4 onef = spu_splats(1.0f); - vec_double2 one = spu_splats(1.0); - vec_uint4 exp_mask = (vec_uint4) { 0x7FF00000, 0, 0x7FF00000, 0 }; - vec_uint4 sign_mask = (vec_uint4) { 0x80000000, 0, 0x80000000, 0}; - vec_uint4 sign_exp_mask = (vec_uint4) { 0xFFF00000, 0, 0xFFF00000,0}; - vec_uint4 frac_mask =(vec_uint4) { 0x000FFFFF, 0xFFFFFFFF, 0x000FFFFF, 0xFFFFFFFF }; - vec_uchar16 swap32 = (vec_uchar16) ((vec_uint4) { 0x04050607, 0x00010203, 0x0C0D0E0F, 0x08090A0B} ); - vec_uint4 zero = (vec_uint4) { 0, 0, 0, 0 }; - vec_int4 e1022 = (vec_int4) { 0x000003FE, 0, 0x000003FE, 0 }; - vec_int4 emax = (vec_int4) { 0x000007FE, 0, 0x000007FE, 0 }; - vec_int4 e1 = (vec_int4) { 0x00000001, 0, 0x00000001, 0 }; - - vec_uint4 nan = (vec_uint4) { 0x7FF80000, 0, 0x7FF80000, 0}; - - /* Extract exponents and underflow denorm arguments to signed zero. + a_nan = spu_testsv(a_in, (SPU_SV_NAN)); + a_zero = spu_testsv(a_in, (SPU_SV_NEG_ZERO | SPU_SV_POS_ZERO)); + a_inf = spu_testsv(a_in, (SPU_SV_NEG_INFINITY | SPU_SV_POS_INFINITY)); + a_denorm = spu_testsv(a_in, (SPU_SV_NEG_DENORM | SPU_SV_POS_DENORM)); + + b_nan = spu_testsv(b_in, (SPU_SV_NAN)); + b_zero = spu_testsv(b_in, (SPU_SV_NEG_ZERO | SPU_SV_POS_ZERO)); + b_inf = spu_testsv(b_in, (SPU_SV_NEG_INFINITY | SPU_SV_POS_INFINITY)); + b_denorm = spu_testsv(b_in, (SPU_SV_NEG_DENORM | SPU_SV_POS_DENORM)); + + /* Scale denorm inputs to into normalized numbers by conditionally scaling the + * input parameters. + */ + a = spu_sel(a_in, spu_mul(a_in, denorm_scale), a_denorm); + b = spu_sel(b_in, spu_mul(b_in, denorm_scale), b_denorm); + +#else /* !__SPU_EDP__ */ + vec_uint4 a_exp, b_exp; + vec_ullong2 a_mant_0, b_mant_0; + vec_ullong2 a_exp_1s, b_exp_1s; + vec_ullong2 sign_exp_mask; + + vec_uint4 exp_mask_u32 = spu_splats((unsigned int)0x7FF00000); + vec_uchar16 splat_hi = (vec_uchar16){0,1,2,3, 0,1,2,3, 8, 9,10,11, 8,9,10,11}; + vec_uchar16 swap_32 = (vec_uchar16){4,5,6,7, 0,1,2,3, 12,13,14,15, 8,9,10,11}; + vec_ullong2 sign_mask = spu_splats(0x8000000000000000ULL); + vec_double2 exp_53 = (vec_double2)spu_splats(0x0350000000000000ULL); + + sign_exp_mask = spu_or(sign_mask, exp_mask); + + /* Extract the floating point components from each of the operands including + * exponent and mantissa. */ - exp_a = spu_and((vec_uint4)a, exp_mask); - exp_b = spu_and((vec_uint4)b, exp_mask); + a_exp = (vec_uint4)spu_and((vec_uint4)a_in, exp_mask_u32); + a_exp = spu_shuffle(a_exp, a_exp, splat_hi); + b_exp = (vec_uint4)spu_and((vec_uint4)b_in, exp_mask_u32); + b_exp = spu_shuffle(b_exp, b_exp, splat_hi); + + a_mant_0 = (vec_ullong2)spu_cmpeq((vec_uint4)spu_andc((vec_ullong2)a_in, sign_exp_mask), 0); + a_mant_0 = spu_and(a_mant_0, spu_shuffle(a_mant_0, a_mant_0, swap_32)); + + b_mant_0 = (vec_ullong2)spu_cmpeq((vec_uint4)spu_andc((vec_ullong2)b_in, sign_exp_mask), 0); + b_mant_0 = spu_and(b_mant_0, spu_shuffle(b_mant_0, b_mant_0, swap_32)); - sign_a = spu_and((vec_uint4)a, sign_mask); - sign_b = spu_and((vec_uint4)b, sign_mask); + a_exp_1s = (vec_ullong2)spu_cmpeq(a_exp, exp_mask_u32); + b_exp_1s = (vec_ullong2)spu_cmpeq(b_exp, exp_mask_u32); - a_exp_0 = spu_cmpeq (exp_a, 0); - utmp = spu_shuffle (a_exp_0, a_exp_0, swap32); - a_exp_0 = spu_and (a_exp_0, utmp); - b_exp_0 = spu_cmpeq (exp_b, 0); - utmp = spu_shuffle (b_exp_0, b_exp_0, swap32); - b_exp_0 = spu_and (b_exp_0, utmp); + /* Identify all possible special values that must be accomodated including: + * +-denorm, +-0, +-infinity, and NaNs. + */ + a_denorm = (vec_ullong2)spu_cmpeq(a_exp, 0); /* really is a_exp_0 */ + a_nan = spu_andc(a_exp_1s, a_mant_0); + a_zero = spu_and (a_denorm, a_mant_0); + a_inf = spu_and (a_exp_1s, a_mant_0); + + b_denorm = (vec_ullong2)spu_cmpeq(b_exp, 0); /* really is b_exp_0 */ + b_nan = spu_andc(b_exp_1s, b_mant_0); + b_zero = spu_and (b_denorm, b_mant_0); + b_inf = spu_and (b_exp_1s, b_mant_0); + + /* Scale denorm inputs to into normalized numbers by conditionally scaling the + * input parameters. + */ + a = spu_sub(spu_or(a_in, exp_53), spu_sel(exp_53, a_in, sign_mask)); + a = spu_sel(a_in, a, a_denorm); - a = spu_sel(a, (vec_double2)sign_a, (vec_ullong2)a_exp_0); - b = spu_sel(b, (vec_double2)sign_b, (vec_ullong2)b_exp_0); + b = spu_sub(spu_or(b_in, exp_53), spu_sel(exp_53, b_in, sign_mask)); + b = spu_sel(b_in, b, b_denorm); - /* Force the divisor and dividend into the range [1.0,2.0). - (Unless they're zero.) - */ - mant_a = spu_sel(a, one, (vec_ullong2)sign_exp_mask); - mant_b = spu_sel(b, one, (vec_ullong2)sign_exp_mask); +#endif /* __SPU_EDP__ */ + + /* Extract the divisor and dividend exponent and force parameters into the signed + * range [1.0,2.0) or [-1.0,2.0). + */ + exp_a = spu_and((vec_ullong2)a, exp_mask); + exp_b = spu_and((vec_ullong2)b, exp_mask); + mant_a = spu_sel(a, one, (vec_ullong2)exp_mask); + mant_b = spu_sel(b, one, (vec_ullong2)exp_mask); + /* Approximate the single reciprocal of b by using * the single precision reciprocal estimate followed by one * single precision iteration of Newton-Raphson. @@ -118,112 +168,67 @@ static __inline vector double _divd2(vector double a, vector double b) inv_bf = spu_re(mant_bf); inv_bf = spu_madd(spu_nmsub(mant_bf, inv_bf, onef), inv_bf, inv_bf); - /* Perform 2 more Newton-Raphson iterations in double precision. + /* Perform 2 more Newton-Raphson iterations in double precision. The + * result (q1) is in the range (0.5, 2.0). */ inv_b = spu_extend(inv_bf); inv_b = spu_madd(spu_nmsub(mant_b, inv_b, one), inv_b, inv_b); q0 = spu_mul(mant_a, inv_b); q1 = spu_madd(spu_nmsub(mant_b, q0, mant_a), inv_b, q0); - /* Compute the quotient's expected exponent. If the exponent - * is out of range, then force the resulting exponent to 0. - * (1023 with the bias). We correct for the out of range - * values by computing a multiplier (mult) that will force the - * result to the correct out of range value and set the - * correct exception flag (UNF, OVF, or neither). - */ - exp_q1 = spu_and((vec_uint4)q1, exp_mask); - exp = spu_sub((vec_int4)exp_a, (vec_int4)exp_b); - exp = spu_rlmaska(exp, -20); // shift right to allow enough bits for working - tmp = spu_rlmaska((vec_int4)exp_q1, -20); - exp = spu_add(exp, tmp); // biased exponent of result (right justified) - - /* The default multiplier is 1.0. If an underflow is detected (the computed - * exponent is less than or equal to a biased 0), force the multiplier to 0.0. - * If exp<=0 set mult = 2**(unbiased exp + 1022) and unbiased exp = -1022 - * = biased 1, the smallest normalized exponent. If exp<-51 set - * mult = 2**(-1074) to ensure underflowing result. Otherwise mult=1. - */ - normal = spu_cmpgt(exp, 0); - nounderflow = spu_cmpgt(exp, -52); - tmp = spu_add(exp, e1022); - mult = (vec_double2)spu_sl(tmp, 20); - mult = spu_sel(mult, one, (vec_ullong2)normal); - mult = spu_sel((vec_double2)e1, mult, (vec_ullong2)nounderflow); - exp = spu_sel(e1, exp, normal); // unbiased -1022 is biased 1 - - /* Force the multiplier to positive infinity (exp_mask) and the biased - * exponent to 1022, if the computed biased exponent is > emax. + + /* Determine the exponent correction factor that must be applied + * to q1 by taking into account the exponent of the normalized inputs + * and the scale factors that were applied to normalize them. */ - overflow = spu_cmpgt(exp, (vec_int4)emax); - exp = spu_sel(exp, (vec_int4)e1022, overflow); - mult = spu_sel(mult, (vec_double2)exp_mask, (vec_ullong2)overflow); - - /* Determine if a, b are Inf, NaN, or zero. - * Since these are rare, it would improve speed if these could be detected - * quickly and a branch used to avoid slowing down the main path. However - * most of the work seems to be in the detection. + exp = spu_rlmaska(spu_sub((vec_int4)exp_a, (vec_int4)exp_b), -20); + exp = spu_add(exp, (vec_int4)spu_add(spu_and((vec_int4)a_denorm, -0x34), spu_and((vec_int4)b_denorm, 0x34))); + + /* Bias the quotient exponent depending on the sign of the exponent correction + * factor so that a single multiplier will ensure the entire double precision + * domain (including denorms) can be achieved. + * + * exp bias q1 adjust exp + * ===== ======== ========== + * positive 2^+65 -65 + * negative 2^-64 +64 */ - a_exp_ones = spu_cmpeq (exp_a, exp_mask); - utmp = spu_shuffle (a_exp_ones, a_exp_ones, swap32); - a_exp_ones = spu_and (a_exp_ones, utmp); - - a_frac = spu_and ((vec_uint4)a, frac_mask); - a_frac_0 = spu_cmpeq (a_frac, 0); - utmp = spu_shuffle (a_frac_0, a_frac_0, swap32); - a_frac_0 = spu_and (a_frac_0, utmp); - - a_zero = spu_and (a_exp_0, a_frac_0); - a_inf = spu_and (a_exp_ones, a_frac_0); - a_nan = spu_andc (a_exp_ones, a_frac_0); - - b_exp_ones = spu_cmpeq (exp_b, exp_mask); - utmp = spu_shuffle (b_exp_ones, b_exp_ones, swap32); - b_exp_ones = spu_and (b_exp_ones, utmp); + exp_bias = spu_xor(spu_rlmaska(exp, -31), 64); - b_frac = spu_and ((vec_uint4)b, frac_mask); - b_frac_0 = spu_cmpeq (b_frac, 0); - utmp = spu_shuffle (b_frac_0, b_frac_0, swap32); - b_frac_0 = spu_and (b_frac_0, utmp); - b_zero = spu_and (b_exp_0, b_frac_0); - b_inf = spu_and (b_exp_ones, b_frac_0); - b_nan = spu_andc (b_exp_ones, b_frac_0); + exp = spu_sub(exp, exp_bias); - /* Handle exception cases */ + q1 = spu_sel(q1, (vec_double2)spu_add((vec_int4)q1, spu_sl(exp_bias, 20)), exp_mask); - /* Result is 0 for 0/x, x!=0, or x/Inf, x!=Inf. - * Set mult=0 for 0/0 or Inf/Inf now, since it will be replaced - * with NaN later. + /* Compute a multiplier (mult) to applied to the quotient (q1) to produce the + * expected result. */ - utmp = spu_or (a_zero, b_inf); - mult = spu_sel(mult, (vec_double2)zero, (vec_ullong2)utmp); - - /* Result is Inf for x/0, x!=0. Set mult=Inf for 0/0 now, since it - * will be replaced with NaN later. + exp = spu_add(exp, 0x3FF); + no_underflow = spu_cmpgt(exp, 0); + overflow = spu_cmpgt(exp, 0x7FF); + exp = spu_and(spu_sl(exp, 20), (vec_int4)no_underflow); + exp = spu_and(exp, (vec_int4)exp_mask); + mult = spu_sel((vec_double2)exp, (vec_double2)exp_mask, (vec_ullong2)overflow); + + /* Handle special value conditions. These include: + * + * 1) IF either operand is a NaN OR both operands are 0 or INFINITY THEN a NaN + * results. + * 2) ELSE IF the dividend is an INFINITY OR the divisor is 0 THEN a INFINITY results. + * 3) ELSE IF the dividend is 0 OR the divisor is INFINITY THEN a 0 results. */ - mult = spu_sel(mult, (vec_double2)exp_mask, (vec_ullong2)b_zero); + mult = spu_andc(mult, (vec_double2)spu_or(a_zero, b_inf)); + mult = spu_sel(mult, (vec_double2)exp_mask, spu_or(a_inf, b_zero)); - /* Result is NaN if either operand is, or Inf/Inf, or 0/0. - */ - res_nan = spu_or (a_nan, b_nan); - utmp = spu_and (a_inf, b_inf); - res_nan = spu_or (res_nan, utmp); - utmp = spu_and (a_zero, b_zero); - res_nan = spu_or (res_nan, utmp); - mult = spu_sel(mult, (vec_double2)nan, (vec_ullong2)res_nan); - - /* Insert sign of result into mult. - */ - sign_res = spu_xor (sign_a, sign_b); - mult = spu_or (mult, (vec_double2)sign_res); + nan = spu_or(a_nan, b_nan); + nan = spu_or(nan, spu_and(a_zero, b_zero)); + nan = spu_or(nan, spu_and(a_inf, b_inf)); - /* Insert the sign and exponent into the result and perform the - * final multiplication. - */ - exp = spu_sl(exp, 20); - q2 = spu_sel(q1, (vec_double2)exp, (vec_ullong2)exp_mask); - q2 = spu_mul(q2, mult); + mult = spu_or(mult, (vec_double2)nan); + + /* Scale the final quotient */ + + q2 = spu_mul(q1, mult); return (q2); } |