diff --git a/ChangeLog b/ChangeLog index f00389e..e1c81e8 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,18 @@ +2007-05-31 Richard B. Kreckel + + * include/cln/integer.h (cl_I_to_E, cl_I_to_UE): New functions. + * src/float/transcendental/cl_LF_exp_aux.cc: Make lq argument an uintE. + * src/float/transcendental/cl_LF_cossin_aux.cc: Likewise. + * src/float/transcendental/cl_LF_coshsinh_aux.cc: Likewise. + * src/float/transcendental/cl_F_tran.h: Change declaration of lq. + * src/float/transcendental/cl_F_lnx.cc: Fix some exponent types. + * src/float/transcendental/cl_F_expx.cc: Likewise. + * src/float/transcendental/cl_F_sinh.cc: Likewise. + * src/float/transcendental/cl_F_atanx.cc: Likewise. + * src/float/transcendental/cl_F_coshsinh.cc: Likewise. + * src/float/transcendental/cl_LF_cossin.cc: Likewise. + * src/float/transcendental/cl_LF_coshsinh.cc: Likewise. + 2007-04-09 Richard B. Kreckel More memory efficient constants: diff --git a/include/cln/integer.h b/include/cln/integer.h index 14012c5..28d5164 100644 --- a/include/cln/integer.h +++ b/include/cln/integer.h @@ -32,17 +32,22 @@ CL_DEFINE_AS_CONVERSION(cl_I) inline unsigned int cl_I_to_uint (const cl_I& x) { return cl_I_to_UL(x); } #endif +// Convert an integer to a 64-bit 'quad' type. +#ifdef intQsize + extern uint64 cl_I_to_UQ (const cl_I& obj); + extern sint64 cl_I_to_Q (const cl_I& obj); +#endif + // Convert an integer to a C `long' or `unsigned long'. #if (long_bitsize==32) inline long cl_I_to_long (const cl_I& x) { return cl_I_to_L(x); } inline unsigned long cl_I_to_ulong (const cl_I& x) { return cl_I_to_UL(x); } #elif (long_bitsize==64) - extern uint64 cl_I_to_UQ (const cl_I& obj); - extern sint64 cl_I_to_Q (const cl_I& obj); inline long cl_I_to_long (const cl_I& x) { return cl_I_to_Q(x); } inline unsigned long cl_I_to_ulong (const cl_I& x) { return cl_I_to_UQ(x); } #endif +// Convert an integer to a counter type. #if (intCsize==long_bitsize) inline uintC cl_I_to_UC (const cl_I& x) { return cl_I_to_ulong(x); } inline sintC cl_I_to_C (const cl_I& x) { return cl_I_to_long(x); } @@ -51,6 +56,15 @@ CL_DEFINE_AS_CONVERSION(cl_I) inline sintC cl_I_to_C (const cl_I& x) { return cl_I_to_int(x); } #endif +// Convert an integer to an exponent type. +#if (intEsize==intLsize) + inline uintE cl_I_to_UE (const cl_I& x) { return cl_I_to_UL(x); } + inline sintE cl_I_to_E (const cl_I& x) { return cl_I_to_L(x); } +#elif (intEsize==intQsize) + inline uintE cl_I_to_UE (const cl_I& x) { return cl_I_to_UQ(x); } + inline sintE cl_I_to_E (const cl_I& x) { return cl_I_to_Q(x); } +#endif + // Logische Operationen auf Integers: diff --git a/src/float/transcendental/cl_F_atanx.cc b/src/float/transcendental/cl_F_atanx.cc index b670853..38fd109 100644 --- a/src/float/transcendental/cl_F_atanx.cc +++ b/src/float/transcendental/cl_F_atanx.cc @@ -195,15 +195,15 @@ static const cl_LF atanx_ratseries (const cl_LF& t) var cl_idecoded_float y_ = integer_decode_float(y); // y = (-1)^sign * 2^exponent * mantissa var uintC lm = integer_length(y_.mantissa); - var uintL me = cl_I_to_UL(- y_.exponent); + var uintE me = cl_I_to_UE(- y_.exponent); var cl_I p; - var uintC lq; + var uintE lq; var cl_boolean last_step = cl_false; if (lm >= me) { // |y| >= 1/2 ? p = y_.sign; // 1 or -1 lq = 1; } else { - var uintL n = me - lm; // |y| < 2^-n with n maximal + var uintE n = me - lm; // |y| < 2^-n with n maximal // Set p to the first n bits of |y|: if (lm > n) { p = y_.mantissa >> (lm - n); @@ -221,7 +221,7 @@ static const cl_LF atanx_ratseries (const cl_LF& t) if (2*n >= lm) last_step = cl_true; } - z = z + scale_float(cl_I_to_LF(p,len),-(sintC)lq); + z = z + scale_float(cl_I_to_LF(p,len),-(sintE)lq); if (last_step) break; var cl_LF_cos_sin_t cis_z = cl_cossin_aux(-p,lq,len); diff --git a/src/float/transcendental/cl_F_coshsinh.cc b/src/float/transcendental/cl_F_coshsinh.cc index 7888285..55b9087 100644 --- a/src/float/transcendental/cl_F_coshsinh.cc +++ b/src/float/transcendental/cl_F_coshsinh.cc @@ -33,7 +33,7 @@ const cosh_sinh_t cosh_sinh (const cl_F& x) // (scale-float (+ y (/ y)) -1) und (scale-float (- y (/ y)) -1) bilden. // Genauigkeit wieder verringern. - var sintL e = float_exponent(x); + var sintE e = float_exponent(x); if (e < 0) { // Exponent e abtesten // e<0 if (zerop(x) || (e <= (1-(sintC)float_digits(x))>>1)) @@ -54,7 +54,7 @@ const cosh_sinh_t cosh_sinh (const cl_F& x) #endif if (TheLfloat(x)->len >= 585) { // verwende exp(x), schneller als cl_coshsinh_ratseries - var cl_LF xx = extend(x,TheLfloat(x)->len+ceiling((uintL)(-e),intDsize)); + var cl_LF xx = extend(x,TheLfloat(x)->len+ceiling((uintE)(-e),intDsize)); var cl_F y = exp(xx); var cl_F y_inv = recip(y); return cosh_sinh_t( diff --git a/src/float/transcendental/cl_F_expx.cc b/src/float/transcendental/cl_F_expx.cc index c7956a3..3f3dfec 100644 --- a/src/float/transcendental/cl_F_expx.cc +++ b/src/float/transcendental/cl_F_expx.cc @@ -139,7 +139,7 @@ const cl_LF expx_ratseries (const cl_LF& x) var uintC len = TheLfloat(x)->len; var cl_idecoded_float x_ = integer_decode_float(x); // x = (-1)^sign * 2^exponent * mantissa - var uintL lq = cl_I_to_UL(- x_.exponent); + var uintE lq = cl_I_to_UE(- x_.exponent); var const cl_I& p = x_.mantissa; // Compute exp(p/2^lq) by splitting into pieces. // Each piece gives rise to a factor exp(pk/2^lqk). @@ -153,12 +153,12 @@ const cl_LF expx_ratseries (const cl_LF& x) // Values 2.5 <= c <= 3.2 seem best. Let's choose c = 2.875. var cl_boolean first_factor = cl_true; var cl_LF product; - var uintL b1; - var uintL b2; + var uintE b1; + var uintE b2; for (b1 = 0, b2 = 1; b1 < lq; b1 = b2, b2 = ceiling(b2*23,8)) { // Piece containing bits b1+1..b2 after "decimal point" // in the binary representation of (p/2^lq). - var uintL lqk = (lq >= b2 ? b2 : lq); + var uintE lqk = (lq >= b2 ? b2 : lq); var cl_I pk = ldb(p,cl_byte(lqk-b1,lq-lqk)); // Compute exp(pk/2^lqk). if (!zerop(pk)) { diff --git a/src/float/transcendental/cl_F_lnx.cc b/src/float/transcendental/cl_F_lnx.cc index e6b2ea7..b00a5d1 100644 --- a/src/float/transcendental/cl_F_lnx.cc +++ b/src/float/transcendental/cl_F_lnx.cc @@ -201,15 +201,15 @@ const cl_LF lnx_ratseries (const cl_LF& x) if (zerop(x1_.mantissa)) break; var uintC lm = integer_length(x1_.mantissa); - var uintL me = cl_I_to_UL(- x1_.exponent); + var uintE me = cl_I_to_UE(- x1_.exponent); var cl_I p; - var uintC lq; + var uintE lq; var cl_boolean last_step = cl_false; if (lm >= me) { // |x'| >= 1/2 ? p = x1_.sign; // 1 or -1 lq = 1; } else { - var uintL n = me - lm; // |x'| < 2^-n with n maximal + var uintE n = me - lm; // |x'| < 2^-n with n maximal // Set p to the first n bits of |x'|: if (lm > n) { p = x1_.mantissa >> (lm - n); @@ -227,7 +227,7 @@ const cl_LF lnx_ratseries (const cl_LF& x) if (2*n >= lm) last_step = cl_true; } - y = y + scale_float(cl_I_to_LF(p,len),-(sintC)lq); + y = y + scale_float(cl_I_to_LF(p,len),-(sintE)lq); if (last_step) break; x = x * cl_exp_aux(-p,lq,len); diff --git a/src/float/transcendental/cl_F_sinh.cc b/src/float/transcendental/cl_F_sinh.cc index b1f0bfc..21d37db 100644 --- a/src/float/transcendental/cl_F_sinh.cc +++ b/src/float/transcendental/cl_F_sinh.cc @@ -41,7 +41,7 @@ const cl_F sinh (const cl_F& x) // verwende exp(x), schneller als cl_coshsinh_ratseries // (aber nur bei 0 > e > -d/2, denn wir müssen, um // Auslöschung zu verhindern, |e| Bits dazunehmen) - var cl_LF xx = extend(x,TheLfloat(x)->len+ceiling((uintL)(-float_exponent(x)),intDsize)); + var cl_LF xx = extend(x,TheLfloat(x)->len+ceiling((uintE)(-float_exponent(x)),intDsize)); var cl_F y = exp(xx); var cl_F z = scale_float(y - recip(y), -1); // (/ (- y (/ y)) 2) return cl_float(z,x); diff --git a/src/float/transcendental/cl_F_tran.h b/src/float/transcendental/cl_F_tran.h index 739226b..29a7a87 100644 --- a/src/float/transcendental/cl_F_tran.h +++ b/src/float/transcendental/cl_F_tran.h @@ -19,7 +19,7 @@ extern const cl_LF pi (uintC len); // computes it even further // cl_exp_aux(p,lq,len) liefert die Zahl exp(p/2^lq) mit len Digits. // 0 < |p| < 2^lq. // Es sollte |p|^2 < 2^lq sein, sonst ist das nicht effizient. -extern const cl_LF cl_exp_aux (const cl_I& p, uintC lq, uintC len); +extern const cl_LF cl_exp_aux (const cl_I& p, uintE lq, uintC len); // cl_cossin_aux(p,lq,len) liefert cos(p/2^lq) und sin(p/2^lq) mit len Digits. // 0 < |p| < 2^lq. @@ -31,7 +31,7 @@ struct cl_LF_cos_sin_t { cl_LF_cos_sin_t (const cl_LF& u, const cl_LF& v) : cos (u), sin (v) {} cl_LF_cos_sin_t () {} }; -extern const cl_LF_cos_sin_t cl_cossin_aux (const cl_I& p, uintC lq, uintC len); +extern const cl_LF_cos_sin_t cl_cossin_aux (const cl_I& p, uintE lq, uintC len); // cl_coshsinh_aux(p,lq,len) liefert cosh(p/2^lq) und sinh(p/2^lq) mit len Digits. // 0 < |p| < 2^lq. @@ -43,7 +43,7 @@ struct cl_LF_cosh_sinh_t { cl_LF_cosh_sinh_t (const cl_LF& u, const cl_LF& v) : cosh (u), sinh (v) {} cl_LF_cosh_sinh_t () {} }; -extern const cl_LF_cosh_sinh_t cl_coshsinh_aux (const cl_I& p, uintC lq, uintC len); +extern const cl_LF_cosh_sinh_t cl_coshsinh_aux (const cl_I& p, uintE lq, uintC len); // atanhx(x) liefert zu einem Float x (betragsmäßig <1/2) atanh(x) als Float. extern const cl_F atanhx (const cl_F& x); diff --git a/src/float/transcendental/cl_LF_coshsinh.cc b/src/float/transcendental/cl_LF_coshsinh.cc index da15148..7a926b7 100644 --- a/src/float/transcendental/cl_LF_coshsinh.cc +++ b/src/float/transcendental/cl_LF_coshsinh.cc @@ -26,17 +26,17 @@ const cl_LF_cosh_sinh_t cl_coshsinh_ratseries (const cl_LF& x) var uintC len = TheLfloat(x)->len; var cl_idecoded_float x_ = integer_decode_float(x); // x = (-1)^sign * 2^exponent * mantissa - var uintL lq = cl_I_to_UL(- x_.exponent); + var uintE lq = cl_I_to_UE(- x_.exponent); var const cl_I& p = x_.mantissa; // Compute sinh(p/2^lq) and cosh(p/2^lq) by splitting into pieces. var cl_boolean first_factor = cl_true; var cl_LF_cosh_sinh_t product; - var uintL b1; - var uintL b2; + var uintE b1; + var uintE b2; for (b1 = 0, b2 = 1; b1 < lq; b1 = b2, b2 = 2*b2) { // Piece containing bits b1+1..b2 after "decimal point" // in the binary representation of (p/2^lq). - var uintL lqk = (lq >= b2 ? b2 : lq); + var uintE lqk = (lq >= b2 ? b2 : lq); var cl_I pk = ldb(p,cl_byte(lqk-b1,lq-lqk)); // Compute sinh(pk/2^lqk) and cosh(pk/2^lqk). if (!zerop(pk)) { diff --git a/src/float/transcendental/cl_LF_coshsinh_aux.cc b/src/float/transcendental/cl_LF_coshsinh_aux.cc index ee68f82..5aebbf1 100644 --- a/src/float/transcendental/cl_LF_coshsinh_aux.cc +++ b/src/float/transcendental/cl_LF_coshsinh_aux.cc @@ -26,10 +26,10 @@ namespace cln { // by a power series evaluation brings 20% speedup, even more for small lengths. #define TRIVIAL_SPEEDUP -const cl_LF_cosh_sinh_t cl_coshsinh_aux (const cl_I& p, uintC lq, uintC len) +const cl_LF_cosh_sinh_t cl_coshsinh_aux (const cl_I& p, uintE lq, uintC len) { { Mutable(cl_I,p); - var uintC lp = integer_length(p); // now |p| < 2^lp. + var uintE lp = integer_length(p); // now |p| < 2^lp. if (!(lp <= lq)) cl_abort(); lp = lq - lp; // now |p/2^lq| < 2^-lp. // Minimize lq (saves computation time). diff --git a/src/float/transcendental/cl_LF_cossin.cc b/src/float/transcendental/cl_LF_cossin.cc index 624e3c6..15d5a81 100644 --- a/src/float/transcendental/cl_LF_cossin.cc +++ b/src/float/transcendental/cl_LF_cossin.cc @@ -26,17 +26,17 @@ const cl_LF_cos_sin_t cl_cossin_ratseries (const cl_LF& x) var uintC len = TheLfloat(x)->len; var cl_idecoded_float x_ = integer_decode_float(x); // x = (-1)^sign * 2^exponent * mantissa - var uintL lq = cl_I_to_UL(- x_.exponent); + var uintE lq = cl_I_to_UE(- x_.exponent); var const cl_I& p = x_.mantissa; // Compute sin(p/2^lq) and cos(p/2^lq) by splitting into pieces. var cl_boolean first_factor = cl_true; var cl_LF_cos_sin_t product; - var uintL b1; - var uintL b2; + var uintE b1; + var uintE b2; for (b1 = 0, b2 = 1; b1 < lq; b1 = b2, b2 = 2*b2) { // Piece containing bits b1+1..b2 after "decimal point" // in the binary representation of (p/2^lq). - var uintL lqk = (lq >= b2 ? b2 : lq); + var uintE lqk = (lq >= b2 ? b2 : lq); var cl_I pk = ldb(p,cl_byte(lqk-b1,lq-lqk)); // Compute sin(pk/2^lqk) and cos(pk/2^lqk). if (!zerop(pk)) { diff --git a/src/float/transcendental/cl_LF_cossin_aux.cc b/src/float/transcendental/cl_LF_cossin_aux.cc index 32d0b77..943158b 100644 --- a/src/float/transcendental/cl_LF_cossin_aux.cc +++ b/src/float/transcendental/cl_LF_cossin_aux.cc @@ -26,10 +26,10 @@ namespace cln { // by a power series evaluation brings 20% speedup, even more for small lengths. #define TRIVIAL_SPEEDUP -const cl_LF_cos_sin_t cl_cossin_aux (const cl_I& p, uintC lq, uintC len) +const cl_LF_cos_sin_t cl_cossin_aux (const cl_I& p, uintE lq, uintC len) { { Mutable(cl_I,p); - var uintC lp = integer_length(p); // now |p| < 2^lp. + var uintE lp = integer_length(p); // now |p| < 2^lp. if (!(lp <= lq)) cl_abort(); lp = lq - lp; // now |p/2^lq| < 2^-lp. // Minimize lq (saves computation time). diff --git a/src/float/transcendental/cl_LF_exp_aux.cc b/src/float/transcendental/cl_LF_exp_aux.cc index 1dda6d2..3401768 100644 --- a/src/float/transcendental/cl_LF_exp_aux.cc +++ b/src/float/transcendental/cl_LF_exp_aux.cc @@ -22,10 +22,10 @@ namespace cln { -const cl_LF cl_exp_aux (const cl_I& p, uintC lq, uintC len) +const cl_LF cl_exp_aux (const cl_I& p, uintE lq, uintC len) { { Mutable(cl_I,p); - var uintC lp = integer_length(p); // now |p| < 2^lp. + var uintE lp = integer_length(p); // now |p| < 2^lp. if (!(lp <= lq)) cl_abort(); lp = lq - lp; // now |p/2^lq| < 2^-lp. // Minimize lq (saves computation time).