From: Richard Kreckel Date: Mon, 11 Dec 2006 17:31:06 +0000 (+0000) Subject: Extend the exponent range from 32 bits to 64 bits on selected platforms. X-Git-Tag: cln_1-2-0~59 X-Git-Url: https://www.ginac.de/CLN/cln.git//cln.git?a=commitdiff_plain;h=e52830748240df6b0ab51d8a78243c88132f6c75;p=cln.git Extend the exponent range from 32 bits to 64 bits on selected platforms. * include/cln/number.h: Add signatures for operations with long long. * include/cln/complex_class.h: Likewise. * include/cln/real_class.h: Likewise. * include/cln/real.h: Likewise. * include/cln/rational_class.h: Likewise. * include/cln/rational.h: Likewise. * include/cln/integer_class.h: Likewise. * include/cln/integer.h: Likewise. * include/cln/float.h: Likewise. * include/cln/lfloat.h: Likewise. * include/cln/types.h (sintE and uintE): New types for exponents. * include/cln/*float.h: Use the new types for exponents. * include/cln/floatformat.h (float_format_t): Make underlying type compatible with sintE. * doc/cln.tex: Document changed float_exponent return value. * src/float/cl_F.h: Likewise. * src/float/ffloat/misc/cl_FF_exponent.cc: Likewise. * src/float/input/cl_F_read.cc: Likewise. * src/float/lfloat/cl_LF.h: Likewise. * src/float/lfloat/cl_LF_impl.h: Likewise. * src/float/lfloat/algebraic/cl_LF_sqrt.cc: Likewise. * src/float/lfloat/elem/cl_LF_1plus.cc: Likewise. * src/float/lfloat/elem/cl_LF_I_div.cc: Likewise. * src/float/lfloat/elem/cl_LF_I_mul.cc: Likewise. * src/float/lfloat/elem/cl_LF_compare.cc: Likewise. * src/float/lfloat/elem/cl_LF_div.cc: Likewise. * src/float/lfloat/elem/cl_LF_from_I.cc: Likewise. * src/float/lfloat/elem/cl_LF_fround.cc: Likewise. * src/float/lfloat/elem/cl_LF_ftrunc.cc: Likewise. * src/float/lfloat/elem/cl_LF_futrunc.cc: Likewise. * src/float/lfloat/elem/cl_LF_mul.cc: Likewise. * src/float/lfloat/elem/cl_LF_scale.cc: Likewise. * src/float/lfloat/elem/cl_LF_scale_I.cc: Likewise. * src/float/lfloat/elem/cl_LF_square.cc: Likewise. * src/float/lfloat/elem/cl_LF_to_I.cc: Likewise. * src/float/lfloat/misc/cl_LF_decode.cc: Likewise. * src/float/lfloat/misc/cl_LF_exponent.cc: Likewise. * src/float/lfloat/misc/cl_LF_idecode.cc: Likewise. * src/float/lfloat/misc/cl_LF_shortenrel.cc: Likewise. * src/float/lfloat/misc/cl_LF_shortenwith.cc: Likewise. * src/float/misc/cl_F_decode.cc: Likewise. * src/float/misc/cl_F_exponent.cc: Likewise. * src/float/misc/cl_F_shortenrel.cc: Likewise. * src/float/misc/cl_float_format.cc: Likewise. * src/float/output/cl_F_dprint.cc: Likewise. * src/float/sfloat/misc/cl_SF_exponent.cc: Likewise. * src/float/transcendental/cl_F_atanhx.cc: Likewise. * src/float/transcendental/cl_F_atanx.cc: Likewise. * src/float/transcendental/cl_F_cosh.cc: Likewise. * src/float/transcendental/cl_F_expx.cc: Likewise. * src/float/transcendental/cl_F_lnx.cc: Likewise. * src/float/transcendental/cl_F_sinhx.cc: Likewise. * src/float/transcendental/cl_F_sinx.cc: Likewise. * src/float/transcendental/cl_LF_pi.cc: Likewise. * src/integer/cl_I.h: Likewise. * src/complex/algebraic/cl_LF_hypot.cc: Likewise. * src/complex/elem/division/cl_C_LF_recip.cc: Likewise. * src/float/dfloat/misc/cl_DF_exponent.cc: Likewise. * src/integer/conv/cl_I_from_Q2.cc: Added. * src/base/cl_low.h (isqrtC): New function, for 64 bit falls back to... * src/base/low/cl_low_isqrt.cc (isqrt): ...this new implementation. * src/base/cl_macros.h (bitc): Make sure 64 bit is used if required by exponent operations. * examples/pi.cc: Support more than 646456614 decimal digits. --- diff --git a/ChangeLog b/ChangeLog index e536e84..5fe28d7 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,71 @@ +2006-12-11 Richard B. Kreckel + + Extend the exponent range from 32 bits to 64 bits on selected platforms. + * include/cln/number.h: Add signatures for operations with long long. + * include/cln/complex_class.h: Likewise. + * include/cln/real_class.h: Likewise. + * include/cln/real.h: Likewise. + * include/cln/rational_class.h: Likewise. + * include/cln/rational.h: Likewise. + * include/cln/integer_class.h: Likewise. + * include/cln/integer.h: Likewise. + * include/cln/float.h: Likewise. + * include/cln/lfloat.h: Likewise. + * include/cln/types.h (sintE and uintE): New types for exponents. + * include/cln/*float.h: Use the new types for exponents. + * include/cln/floatformat.h (float_format_t): Make underlying type + compatible with sintE. + * doc/cln.tex: Document changed float_exponent return value. + * src/float/cl_F.h: Likewise. + * src/float/ffloat/misc/cl_FF_exponent.cc: Likewise. + * src/float/input/cl_F_read.cc: Likewise. + * src/float/lfloat/cl_LF.h: Likewise. + * src/float/lfloat/cl_LF_impl.h: Likewise. + * src/float/lfloat/algebraic/cl_LF_sqrt.cc: Likewise. + * src/float/lfloat/elem/cl_LF_1plus.cc: Likewise. + * src/float/lfloat/elem/cl_LF_I_div.cc: Likewise. + * src/float/lfloat/elem/cl_LF_I_mul.cc: Likewise. + * src/float/lfloat/elem/cl_LF_compare.cc: Likewise. + * src/float/lfloat/elem/cl_LF_div.cc: Likewise. + * src/float/lfloat/elem/cl_LF_from_I.cc: Likewise. + * src/float/lfloat/elem/cl_LF_fround.cc: Likewise. + * src/float/lfloat/elem/cl_LF_ftrunc.cc: Likewise. + * src/float/lfloat/elem/cl_LF_futrunc.cc: Likewise. + * src/float/lfloat/elem/cl_LF_mul.cc: Likewise. + * src/float/lfloat/elem/cl_LF_scale.cc: Likewise. + * src/float/lfloat/elem/cl_LF_scale_I.cc: Likewise. + * src/float/lfloat/elem/cl_LF_square.cc: Likewise. + * src/float/lfloat/elem/cl_LF_to_I.cc: Likewise. + * src/float/lfloat/misc/cl_LF_decode.cc: Likewise. + * src/float/lfloat/misc/cl_LF_exponent.cc: Likewise. + * src/float/lfloat/misc/cl_LF_idecode.cc: Likewise. + * src/float/lfloat/misc/cl_LF_shortenrel.cc: Likewise. + * src/float/lfloat/misc/cl_LF_shortenwith.cc: Likewise. + * src/float/misc/cl_F_decode.cc: Likewise. + * src/float/misc/cl_F_exponent.cc: Likewise. + * src/float/misc/cl_F_shortenrel.cc: Likewise. + * src/float/misc/cl_float_format.cc: Likewise. + * src/float/output/cl_F_dprint.cc: Likewise. + * src/float/sfloat/misc/cl_SF_exponent.cc: Likewise. + * src/float/transcendental/cl_F_atanhx.cc: Likewise. + * src/float/transcendental/cl_F_atanx.cc: Likewise. + * src/float/transcendental/cl_F_cosh.cc: Likewise. + * src/float/transcendental/cl_F_expx.cc: Likewise. + * src/float/transcendental/cl_F_lnx.cc: Likewise. + * src/float/transcendental/cl_F_sinhx.cc: Likewise. + * src/float/transcendental/cl_F_sinx.cc: Likewise. + * src/float/transcendental/cl_LF_pi.cc: Likewise. + * src/integer/cl_I.h: Likewise. + * src/complex/algebraic/cl_LF_hypot.cc: Likewise. + * src/complex/elem/division/cl_C_LF_recip.cc: Likewise. + * src/float/dfloat/misc/cl_DF_exponent.cc: Likewise. + * src/integer/conv/cl_I_from_Q2.cc: Added. + * src/base/cl_low.h (isqrtC): New function, for 64 bit falls back to... + * src/base/low/cl_low_isqrt.cc (isqrt): ...this new implementation. + * src/base/cl_macros.h (bitc): Make sure 64 bit is used if required by + exponent operations. + * examples/pi.cc: Support more than 646456614 decimal digits. + 2006-11-02 Richard B. Kreckel * src/base/digitseq/cl_DS.h: #undef DS, needed for i386-Solaris. diff --git a/NEWS b/NEWS index eded1a9..9a9f7ca 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,32 @@ +2006-mm-dd, version 1.2.0 +========================= + +Implementation changes +---------------------- + +* Added support for huge numbers... + + +2006-08-08, version 1.1.13 +========================== + +* Compilation fixes for 64-bit brokenness introduced in last release. + + +2006-08-06, version 1.1.12 +========================== + +Implementation changes +---------------------- + +* Fix rare assertion when printing quite large floats. + +Other changes +------------- + +* Compilation fixes for several platforms: *BSD, Intel Mac, and MinGW. + + 2005-11-23, version 1.1.11 ========================== diff --git a/doc/cln.tex b/doc/cln.tex index a747bef..82c8e99 100644 --- a/doc/cln.tex +++ b/doc/cln.tex @@ -2025,7 +2025,7 @@ The following functions provide an abstract interface to the underlying representation of floating-point numbers. @table @code -@item sintL float_exponent (const @var{type}& x) +@item sintE float_exponent (const @var{type}& x) @cindex @code{float_exponent ()} Returns the exponent @code{e} of @code{x}. For @code{x = 0.0}, this is 0. For @code{x} non-zero, this is the unique @@ -2118,7 +2118,7 @@ The type @code{float_format_t} describes a floating-point format. @cindex @code{float_format_t} @table @code -@item float_format_t float_format (uintL n) +@item float_format_t float_format (uintE n) @cindex @code{float_format ()} Returns the smallest float format which guarantees at least @code{n} decimal digits in the mantissa (after the decimal point). diff --git a/examples/pi.cc b/examples/pi.cc index a2768b5..815aea4 100644 --- a/examples/pi.cc +++ b/examples/pi.cc @@ -23,7 +23,7 @@ usage (ostream &os) int main (int argc, char * argv[]) { - int digits = 100; + long digits = 100; if (argc > 1) { if (argc == 2 && !strcmp(argv[1],"--help")) { usage(cout); @@ -46,7 +46,7 @@ main (int argc, char * argv[]) return 0; } if (argc == 2 && isdigit(argv[1][0])) { - digits = atoi(argv[1]); + digits = atol(argv[1]); } else { usage(cerr); return 1; diff --git a/include/cln/complex_class.h b/include/cln/complex_class.h index b5dec56..4e7a774 100644 --- a/include/cln/complex_class.h +++ b/include/cln/complex_class.h @@ -21,6 +21,10 @@ public: cl_N (const unsigned int); // argument must be < 2^29 cl_N (const long); cl_N (const unsigned long); +#ifdef HAVE_LONGLONG + cl_N (const long long); + cl_N (const unsigned long long); +#endif cl_N (const float); cl_N (const double); cl_N& operator= (const int); // |argument| must be < 2^29 @@ -29,6 +33,10 @@ public: cl_N& operator= (const unsigned long); cl_N& operator= (const float); cl_N& operator= (const double); +#ifdef HAVE_LONGLONG + cl_N& operator= (const long long); + cl_N& operator= (const unsigned long long); +#endif // Other constructors. cl_N (const char *); // Private constructor. @@ -56,7 +64,10 @@ CL_DEFINE_INT_CONSTRUCTORS(cl_N) CL_DEFINE_INT_ASSIGNMENT_OPERATORS(cl_N) CL_DEFINE_LONG_CONSTRUCTORS(cl_N) CL_DEFINE_LONG_ASSIGNMENT_OPERATORS(cl_N) -// Constructors and assignment operators from C numeric types. +#ifdef HAVE_LONGLONG +CL_DEFINE_LONGLONG_CONSTRUCTORS(cl_N) +CL_DEFINE_LONGLONG_ASSIGNMENT_OPERATORS(cl_N) +#endif CL_DEFINE_FLOAT_CONSTRUCTOR(cl_N) CL_DEFINE_DOUBLE_CONSTRUCTOR(cl_N) diff --git a/include/cln/dfloat.h b/include/cln/dfloat.h index 2b9d0c7..30534b5 100644 --- a/include/cln/dfloat.h +++ b/include/cln/dfloat.h @@ -243,7 +243,7 @@ extern const decoded_dfloat decode_float (const cl_DF& x); // den Exponenten von (decode-float x). // x = 0.0 liefert 0. // x = (-1)^s * 2^e * m liefert e. -extern sintL float_exponent (const cl_DF& x); +extern sintE float_exponent (const cl_DF& x); // float_radix(x) liefert (float-radix x), wo x ein Float ist. inline sintL float_radix (const cl_DF& x) diff --git a/include/cln/ffloat.h b/include/cln/ffloat.h index d690327..a9ade1c 100644 --- a/include/cln/ffloat.h +++ b/include/cln/ffloat.h @@ -243,7 +243,7 @@ extern const decoded_ffloat decode_float (const cl_FF& x); // den Exponenten von (decode-float x). // x = 0.0 liefert 0. // x = (-1)^s * 2^e * m liefert e. -extern sintL float_exponent (const cl_FF& x); +extern sintE float_exponent (const cl_FF& x); // float_radix(x) liefert (float-radix x), wo x ein Float ist. inline sintL float_radix (const cl_FF& x) diff --git a/include/cln/float.h b/include/cln/float.h index 4f9945c..abee49b 100644 --- a/include/cln/float.h +++ b/include/cln/float.h @@ -59,7 +59,7 @@ extern float_format_t default_float_format; // Returns the smallest float format which guarantees at least n decimal digits // in the mantissa (after the decimal point). -extern float_format_t float_format (uintL n); +extern float_format_t float_format (uintE n); // cl_float(x,y) wandelt ein Float x in das Float-Format des Floats y um // und rundet dabei nötigenfalls. @@ -172,6 +172,12 @@ inline const cl_F operator+ (const long x, const cl_F& y) { return cl_I(x) + y; } inline const cl_F operator+ (const unsigned long x, const cl_F& y) { return cl_I(x) + y; } +#ifdef HAVE_LONGLONG +inline const cl_F operator+ (const long long x, const cl_F& y) + { return cl_I(x) + y; } +inline const cl_F operator+ (const unsigned long long x, const cl_F& y) + { return cl_I(x) + y; } +#endif inline const cl_F operator+ (const float x, const cl_F& y) { return cl_F(x) + y; } inline const cl_F operator+ (const double x, const cl_F& y) @@ -184,6 +190,12 @@ inline const cl_F operator+ (const cl_F& x, const long y) { return x + cl_I(y); } inline const cl_F operator+ (const cl_F& x, const unsigned long y) { return x + cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_F operator+ (const cl_F& x, const long long y) + { return x + cl_I(y); } +inline const cl_F operator+ (const cl_F& x, const unsigned long long y) + { return x + cl_I(y); } +#endif inline const cl_F operator+ (const cl_F& x, const float y) { return x + cl_F(y); } inline const cl_F operator+ (const cl_F& x, const double y) @@ -209,6 +221,12 @@ inline const cl_F operator- (const long x, const cl_F& y) { return cl_I(x) - y; } inline const cl_F operator- (const unsigned long x, const cl_F& y) { return cl_I(x) - y; } +#ifdef HAVE_LONGLONG +inline const cl_F operator- (const long long x, const cl_F& y) + { return cl_I(x) - y; } +inline const cl_F operator- (const unsigned long long x, const cl_F& y) + { return cl_I(x) - y; } +#endif inline const cl_F operator- (const float x, const cl_F& y) { return cl_F(x) - y; } inline const cl_F operator- (const double x, const cl_F& y) @@ -221,6 +239,12 @@ inline const cl_F operator- (const cl_F& x, const long y) { return x - cl_I(y); } inline const cl_F operator- (const cl_F& x, const unsigned long y) { return x - cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_F operator- (const cl_F& x, const long long y) + { return x - cl_I(y); } +inline const cl_F operator- (const cl_F& x, const unsigned long long y) + { return x - cl_I(y); } +#endif inline const cl_F operator- (const cl_F& x, const float y) { return x - cl_F(y); } inline const cl_F operator- (const cl_F& x, const double y) @@ -258,6 +282,12 @@ inline const cl_R operator* (const long x, const cl_F& y) { return cl_I(x) * y; } inline const cl_R operator* (const unsigned long x, const cl_F& y) { return cl_I(x) * y; } +#ifdef HAVE_LONGLONG +inline const cl_R operator* (const long long x, const cl_F& y) + { return cl_I(x) * y; } +inline const cl_R operator* (const unsigned long long x, const cl_F& y) + { return cl_I(x) * y; } +#endif inline const cl_F operator* (const float x, const cl_F& y) { return cl_F(x) * y; } inline const cl_F operator* (const double x, const cl_F& y) @@ -270,6 +300,12 @@ inline const cl_R operator* (const cl_F& x, const long y) { return x * cl_I(y); } inline const cl_R operator* (const cl_F& x, const unsigned long y) { return x * cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_R operator* (const cl_F& x, const long long y) + { return x * cl_I(y); } +inline const cl_R operator* (const cl_F& x, const unsigned long long y) + { return x * cl_I(y); } +#endif inline const cl_F operator* (const cl_F& x, const float y) { return x * cl_F(y); } inline const cl_F operator* (const cl_F& x, const double y) @@ -294,6 +330,12 @@ inline const cl_F operator/ (const cl_F& x, const long y) { return x / cl_I(y); } inline const cl_F operator/ (const cl_F& x, const unsigned long y) { return x / cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_F operator/ (const cl_F& x, const long long y) + { return x / cl_I(y); } +inline const cl_F operator/ (const cl_F& x, const unsigned long long y) + { return x / cl_I(y); } +#endif inline const cl_F operator/ (const cl_F& x, const float y) { return x / cl_F(y); } inline const cl_F operator/ (const cl_F& x, const double y) @@ -306,6 +348,12 @@ inline const cl_R operator/ (const long x, const cl_F& y) { return cl_I(x) / y; } inline const cl_R operator/ (const unsigned long x, const cl_F& y) { return cl_I(x) / y; } +#ifdef HAVE_LONGLONG +inline const cl_R operator/ (const long long x, const cl_F& y) + { return cl_I(x) / y; } +inline const cl_R operator/ (const unsigned long long x, const cl_F& y) + { return cl_I(x) / y; } +#endif inline const cl_F operator/ (const float x, const cl_F& y) { return cl_F(x) / y; } inline const cl_F operator/ (const double x, const cl_F& y) @@ -451,7 +499,7 @@ extern const decoded_float decode_float (const cl_F& x); // den Exponenten von (decode-float x). // x = 0.0 liefert 0. // x = (-1)^s * 2^e * m liefert e. -extern sintL float_exponent (const cl_F& x); +extern sintE float_exponent (const cl_F& x); // float_radix(x) liefert (float-radix x), wo x ein Float ist. inline sintL float_radix (const cl_F& x) diff --git a/include/cln/floatformat.h b/include/cln/floatformat.h index f41200b..7a84c2e 100644 --- a/include/cln/floatformat.h +++ b/include/cln/floatformat.h @@ -12,7 +12,8 @@ enum float_format_t { float_format_sfloat = 17, float_format_ffloat = 24, float_format_dfloat = 53, - float_format_lfloat_min = ((53+intDsize-1)/intDsize)*intDsize // = round_up(53,intDsize) + float_format_lfloat_min = ((53+intDsize-1)/intDsize)*intDsize, // = round_up(53,intDsize) + float_format_lfloat_max = ~((sintE)(1) << intEsize-1) // force correct underlying type of float_format_t }; } // namespace cln diff --git a/include/cln/integer.h b/include/cln/integer.h index 72b848f..14012c5 100644 --- a/include/cln/integer.h +++ b/include/cln/integer.h @@ -43,6 +43,14 @@ CL_DEFINE_AS_CONVERSION(cl_I) inline unsigned long cl_I_to_ulong (const cl_I& x) { return cl_I_to_UQ(x); } #endif +#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); } +#elif (intCsize==int_bitsize) + inline uintC cl_I_to_UC (const cl_I& x) { return cl_I_to_uint(x); } + inline sintC cl_I_to_C (const cl_I& x) { return cl_I_to_int(x); } +#endif + // Logische Operationen auf Integers: @@ -191,6 +199,12 @@ inline const cl_I operator+ (const long x, const cl_I& y) { return cl_I(x) + y; } inline const cl_I operator+ (const unsigned long x, const cl_I& y) { return cl_I(x) + y; } +#ifdef HAVE_LONGLONG +inline const cl_I operator+ (const long long x, const cl_I& y) + { return cl_I(x) + y; } +inline const cl_I operator+ (const unsigned long long x, const cl_I& y) + { return cl_I(x) + y; } +#endif inline const cl_I operator+ (const cl_I& x, const int y) { return x + cl_I(y); } inline const cl_I operator+ (const cl_I& x, const unsigned int y) @@ -199,6 +213,12 @@ inline const cl_I operator+ (const cl_I& x, const long y) { return x + cl_I(y); } inline const cl_I operator+ (const cl_I& x, const unsigned long y) { return x + cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_I operator+ (const cl_I& x, const long long y) + { return x + cl_I(y); } +inline const cl_I operator+ (const cl_I& x, const unsigned long long y) + { return x + cl_I(y); } +#endif // (- x), wenn x ein Integer ist. Ergebnis Integer. extern const cl_I operator- (const cl_I& x); @@ -214,6 +234,12 @@ inline const cl_I operator- (const long x, const cl_I& y) { return cl_I(x) - y; } inline const cl_I operator- (const unsigned long x, const cl_I& y) { return cl_I(x) - y; } +#ifdef HAVE_LONGLONG +inline const cl_I operator- (const long long x, const cl_I& y) + { return cl_I(x) - y; } +inline const cl_I operator- (const unsigned long long x, const cl_I& y) + { return cl_I(x) - y; } +#endif inline const cl_I operator- (const cl_I& x, const int y) { return x - cl_I(y); } inline const cl_I operator- (const cl_I& x, const unsigned int y) @@ -222,6 +248,12 @@ inline const cl_I operator- (const cl_I& x, const long y) { return x - cl_I(y); } inline const cl_I operator- (const cl_I& x, const unsigned long y) { return x - cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_I operator- (const cl_I& x, const long long y) + { return x - cl_I(y); } +inline const cl_I operator- (const cl_I& x, const unsigned long long y) + { return x - cl_I(y); } +#endif // (abs x), wenn x ein Integer ist. Ergebnis Integer. extern const cl_I abs (const cl_I& x); @@ -310,6 +342,12 @@ inline const cl_I operator* (const long x, const cl_I& y) { return cl_I(x) * y; } inline const cl_I operator* (const unsigned long x, const cl_I& y) { return cl_I(x) * y; } +#ifdef HAVE_LONGLONG +inline const cl_I operator* (const long long x, const cl_I& y) + { return cl_I(x) * y; } +inline const cl_I operator* (const unsigned long long x, const cl_I& y) + { return cl_I(x) * y; } +#endif inline const cl_I operator* (const cl_I& x, const int y) { return x * cl_I(y); } inline const cl_I operator* (const cl_I& x, const unsigned int y) @@ -318,6 +356,12 @@ inline const cl_I operator* (const cl_I& x, const long y) { return x * cl_I(y); } inline const cl_I operator* (const cl_I& x, const unsigned long y) { return x * cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_I operator* (const cl_I& x, const long long y) + { return x * cl_I(y); } +inline const cl_I operator* (const cl_I& x, const unsigned long long y) + { return x * cl_I(y); } +#endif // (EXPT x 2), wo x Integer ist. extern const cl_I square (const cl_I& x); diff --git a/include/cln/integer_class.h b/include/cln/integer_class.h index db2068e..23f5632 100644 --- a/include/cln/integer_class.h +++ b/include/cln/integer_class.h @@ -21,10 +21,18 @@ public: cl_I (const unsigned int); // argument must be < 2^29 cl_I (const long); cl_I (const unsigned long); +#ifdef HAVE_LONGLONG + cl_I (const long long); + cl_I (const unsigned long long); +#endif cl_I& operator= (const int); // |argument| must be < 2^29 cl_I& operator= (const unsigned int); // argument must be < 2^29 cl_I& operator= (const long); cl_I& operator= (const unsigned long); +#ifdef HAVE_LONGLONG + cl_I& operator= (const long long); + cl_I& operator= (const unsigned long long); +#endif // Other constructors. cl_I (const char *); // Private constructor. @@ -51,6 +59,10 @@ CL_DEFINE_INT_CONSTRUCTORS(cl_I) CL_DEFINE_INT_ASSIGNMENT_OPERATORS(cl_I) CL_DEFINE_LONG_CONSTRUCTORS(cl_I) CL_DEFINE_LONG_ASSIGNMENT_OPERATORS(cl_I) +#ifdef HAVE_LONGLONG +CL_DEFINE_LONGLONG_CONSTRUCTORS(cl_I) +CL_DEFINE_LONGLONG_ASSIGNMENT_OPERATORS(cl_I) +#endif } // namespace cln diff --git a/include/cln/lfloat.h b/include/cln/lfloat.h index 24041d8..58d1e89 100644 --- a/include/cln/lfloat.h +++ b/include/cln/lfloat.h @@ -83,6 +83,12 @@ inline const cl_R operator* (const long x, const cl_LF& y) { return cl_I(x) * y; } inline const cl_R operator* (const unsigned long x, const cl_LF& y) { return cl_I(x) * y; } +#ifdef HAVE_LONGLONG +inline const cl_R operator* (const long long x, const cl_LF& y) + { return cl_I(x) * y; } +inline const cl_R operator* (const unsigned long long x, const cl_LF& y) + { return cl_I(x) * y; } +#endif inline const cl_R operator* (const cl_LF& x, const int y) { return x * cl_I(y); } inline const cl_R operator* (const cl_LF& x, const unsigned int y) @@ -91,6 +97,12 @@ inline const cl_R operator* (const cl_LF& x, const long y) { return x * cl_I(y); } inline const cl_R operator* (const cl_LF& x, const unsigned long y) { return x * cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_R operator* (const cl_LF& x, const long long y) + { return x * cl_I(y); } +inline const cl_R operator* (const cl_LF& x, const unsigned long long y) + { return x * cl_I(y); } +#endif // Spezialfall x = y. // Liefert zu einem Long-Float x : (* x x), ein LF. extern const cl_LF square (const cl_LF& x); @@ -127,6 +139,12 @@ inline const cl_LF operator/ (const cl_LF& x, const long y) { return x / cl_I(y); } inline const cl_LF operator/ (const cl_LF& x, const unsigned long y) { return x / cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_LF operator/ (const cl_LF& x, const long long y) + { return x / cl_I(y); } +inline const cl_LF operator/ (const cl_LF& x, const unsigned long long y) + { return x / cl_I(y); } +#endif inline const cl_R operator/ (const int x, const cl_LF& y) { return cl_I(x) / y; } inline const cl_R operator/ (const unsigned int x, const cl_LF& y) @@ -135,6 +153,12 @@ inline const cl_R operator/ (const long x, const cl_LF& y) { return cl_I(x) / y; } inline const cl_R operator/ (const unsigned long x, const cl_LF& y) { return cl_I(x) / y; } +#ifdef HAVE_LONGLONG +inline const cl_R operator/ (const long long x, const cl_LF& y) + { return cl_I(x) / y; } +inline const cl_R operator/ (const unsigned long long x, const cl_LF& y) + { return cl_I(x) / y; } +#endif // Liefert zu einem Long-Float x>=0 : (sqrt x), ein LF. extern const cl_LF sqrt (const cl_LF& x); @@ -333,7 +357,7 @@ extern const decoded_lfloat decode_float (const cl_LF& x); // den Exponenten von (decode-float x). // x = 0.0 liefert 0. // x = (-1)^s * 2^e * m liefert e. -extern sintL float_exponent (const cl_LF& x); +extern sintE float_exponent (const cl_LF& x); // float_radix(x) liefert (float-radix x), wo x ein Float ist. inline sintL float_radix (const cl_LF& x) diff --git a/include/cln/number.h b/include/cln/number.h index aff318b..da0b2fd 100644 --- a/include/cln/number.h +++ b/include/cln/number.h @@ -109,6 +109,37 @@ inline _class_& _class_::operator= (const unsigned long wert) \ } #endif +#ifdef HAVE_LONGLONG +#if (long_long_bitsize==64) +// `long' == `sintQ', `unsigned long' == `uintQ'. +#define CL_DEFINE_LONGLONG_CONSTRUCTORS(_class_) \ +inline _class_::_class_ (const long long wert) \ +{ \ + extern cl_private_thing cl_I_constructor_from_Q (sint64 wert); \ + pointer = cl_I_constructor_from_Q(wert); \ +} \ +inline _class_::_class_ (const unsigned long long wert) \ +{ \ + extern cl_private_thing cl_I_constructor_from_UQ (uint64 wert); \ + pointer = cl_I_constructor_from_UQ(wert); \ +} +#define CL_DEFINE_LONGLONG_ASSIGNMENT_OPERATORS(_class_) \ +inline _class_& _class_::operator= (const long long wert) \ +{ \ + extern cl_private_thing cl_I_constructor_from_Q (sint64 wert); \ + cl_dec_refcount(*this); \ + pointer = cl_I_constructor_from_Q(wert); \ + return *this; \ +} \ +inline _class_& _class_::operator= (const unsigned long long wert) \ +{ \ + extern cl_private_thing cl_I_constructor_from_UQ (uint64 wert); \ + cl_dec_refcount(*this); \ + pointer = cl_I_constructor_from_UQ(wert); \ + return *this; \ +} +#endif +#endif namespace cln { @@ -164,6 +195,10 @@ public: cl_number (const unsigned int); // argument must be < 2^29 cl_number (const long); cl_number (const unsigned long); +#ifdef HAVE_LONGLONG + cl_number (const long long); + cl_number (const unsigned long long); +#endif cl_number (const float); cl_number (const double); cl_number& operator= (const int); // |argument| must be < 2^29 @@ -172,6 +207,10 @@ public: cl_number& operator= (const unsigned long); cl_number& operator= (const float); cl_number& operator= (const double); +#ifdef HAVE_LONGLONG + cl_number& operator= (const long long); + cl_number& operator= (const unsigned long long); +#endif // Other constructors. // cl_number (const char *); // Private pointer manipulations. @@ -192,6 +231,10 @@ CL_DEFINE_INT_CONSTRUCTORS(cl_number) CL_DEFINE_INT_ASSIGNMENT_OPERATORS(cl_number) CL_DEFINE_LONG_CONSTRUCTORS(cl_number) CL_DEFINE_LONG_ASSIGNMENT_OPERATORS(cl_number) +#ifdef HAVE_LONGLONG +CL_DEFINE_LONGLONG_CONSTRUCTORS(cl_number) +CL_DEFINE_LONGLONG_ASSIGNMENT_OPERATORS(cl_number) +#endif CL_DEFINE_FLOAT_CONSTRUCTOR(cl_number) CL_DEFINE_DOUBLE_CONSTRUCTOR(cl_number) diff --git a/include/cln/rational.h b/include/cln/rational.h index fdc5beb..6b50855 100644 --- a/include/cln/rational.h +++ b/include/cln/rational.h @@ -33,6 +33,12 @@ inline const cl_RA operator+ (const long x, const cl_RA& y) { return cl_I(x) + y; } inline const cl_RA operator+ (const unsigned long x, const cl_RA& y) { return cl_I(x) + y; } +#ifdef HAVE_LONGLONG +inline const cl_RA operator+ (const long long x, const cl_RA& y) + { return cl_I(x) + y; } +inline const cl_RA operator+ (const unsigned long long x, const cl_RA& y) + { return cl_I(x) + y; } +#endif inline const cl_RA operator+ (const cl_RA& x, const int y) { return x + cl_I(y); } inline const cl_RA operator+ (const cl_RA& x, const unsigned int y) @@ -41,6 +47,12 @@ inline const cl_RA operator+ (const cl_RA& x, const long y) { return x + cl_I(y); } inline const cl_RA operator+ (const cl_RA& x, const unsigned long y) { return x + cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_RA operator+ (const cl_RA& x, const long long y) + { return x + cl_I(y); } +inline const cl_RA operator+ (const cl_RA& x, const unsigned long long y) + { return x + cl_I(y); } +#endif // (- r s), wo r und s rationale Zahlen sind. extern const cl_RA operator- (const cl_RA& r, const cl_RA& s); @@ -53,6 +65,12 @@ inline const cl_RA operator- (const long x, const cl_RA& y) { return cl_I(x) - y; } inline const cl_RA operator- (const unsigned long x, const cl_RA& y) { return cl_I(x) - y; } +#ifdef HAVE_LONGLONG +inline const cl_RA operator- (const long long x, const cl_RA& y) + { return cl_I(x) - y; } +inline const cl_RA operator- (const unsigned long long x, const cl_RA& y) + { return cl_I(x) - y; } +#endif inline const cl_RA operator- (const cl_RA& x, const int y) { return x - cl_I(y); } inline const cl_RA operator- (const cl_RA& x, const unsigned int y) @@ -61,6 +79,12 @@ inline const cl_RA operator- (const cl_RA& x, const long y) { return x - cl_I(y); } inline const cl_RA operator- (const cl_RA& x, const unsigned long y) { return x - cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_RA operator- (const cl_RA& x, const long long y) + { return x - cl_I(y); } +inline const cl_RA operator- (const cl_RA& x, const unsigned long long y) + { return x - cl_I(y); } +#endif // (1+ r), wo r eine rationale Zahl ist. extern const cl_RA plus1 (const cl_RA& r); @@ -116,6 +140,12 @@ inline const cl_RA operator* (const long x, const cl_RA& y) { return cl_I(x) * y; } inline const cl_RA operator* (const unsigned long x, const cl_RA& y) { return cl_I(x) * y; } +#ifdef HAVE_LONGLONG +inline const cl_RA operator* (const long long x, const cl_RA& y) + { return cl_I(x) * y; } +inline const cl_RA operator* (const unsigned long long x, const cl_RA& y) + { return cl_I(x) * y; } +#endif inline const cl_RA operator* (const cl_RA& x, const int y) { return x * cl_I(y); } inline const cl_RA operator* (const cl_RA& x, const unsigned int y) @@ -124,6 +154,12 @@ inline const cl_RA operator* (const cl_RA& x, const long y) { return x * cl_I(y); } inline const cl_RA operator* (const cl_RA& x, const unsigned long y) { return x * cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_RA operator* (const cl_RA& x, const long long y) + { return x * cl_I(y); } +inline const cl_RA operator* (const cl_RA& x, const unsigned long long y) + { return x * cl_I(y); } +#endif // Quadrat (* r r), wo r eine rationale Zahl ist. extern const cl_RA square (const cl_RA& r); @@ -139,6 +175,12 @@ inline const cl_RA operator/ (const long x, const cl_RA& y) { return cl_I(x) / y; } inline const cl_RA operator/ (const unsigned long x, const cl_RA& y) { return cl_I(x) / y; } +#ifdef HAVE_LONGLONG +inline const cl_RA operator/ (const long long x, const cl_RA& y) + { return cl_I(x) / y; } +inline const cl_RA operator/ (const unsigned long long x, const cl_RA& y) + { return cl_I(x) / y; } +#endif inline const cl_RA operator/ (const cl_RA& x, const int y) { return x / cl_I(y); } inline const cl_RA operator/ (const cl_RA& x, const unsigned int y) @@ -147,6 +189,12 @@ inline const cl_RA operator/ (const cl_RA& x, const long y) { return x / cl_I(y); } inline const cl_RA operator/ (const cl_RA& x, const unsigned long y) { return x / cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_RA operator/ (const cl_RA& x, const long long y) + { return x / cl_I(y); } +inline const cl_RA operator/ (const cl_RA& x, const unsigned long long y) + { return x / cl_I(y); } +#endif // Return type for rounding operators. // x / y --> (q,r) with x = y*q+r. @@ -273,6 +321,10 @@ inline cl_RA& operator+= (cl_RA& x, const int y) { return x = x + y; } inline cl_RA& operator+= (cl_RA& x, const unsigned int y) { return x = x + y; } inline cl_RA& operator+= (cl_RA& x, const long y) { return x = x + y; } inline cl_RA& operator+= (cl_RA& x, const unsigned long y) { return x = x + y; } +#ifdef HAVE_LONGLONG +inline cl_RA& operator+= (cl_RA& x, const long long y) { return x = x + y; } +inline cl_RA& operator+= (cl_RA& x, const unsigned long long y) { return x = x + y; } +#endif inline cl_RA& operator++ /* prefix */ (cl_RA& x) { return x = plus1(x); } inline void operator++ /* postfix */ (cl_RA& x, int dummy) { (void)dummy; x = plus1(x); } inline cl_RA& operator-= (cl_RA& x, const cl_RA& y) { return x = x - y; } @@ -280,6 +332,10 @@ inline cl_RA& operator-= (cl_RA& x, const int y) { return x = x - y; } inline cl_RA& operator-= (cl_RA& x, const unsigned int y) { return x = x - y; } inline cl_RA& operator-= (cl_RA& x, const long y) { return x = x - y; } inline cl_RA& operator-= (cl_RA& x, const unsigned long y) { return x = x - y; } +#ifdef HAVE_LONGLONG +inline cl_RA& operator-= (cl_RA& x, const long long y) { return x = x - y; } +inline cl_RA& operator-= (cl_RA& x, const unsigned long long y) { return x = x - y; } +#endif inline cl_RA& operator-- /* prefix */ (cl_RA& x) { return x = minus1(x); } inline void operator-- /* postfix */ (cl_RA& x, int dummy) { (void)dummy; x = minus1(x); } inline cl_RA& operator*= (cl_RA& x, const cl_RA& y) { return x = x * y; } diff --git a/include/cln/rational_class.h b/include/cln/rational_class.h index 3546598..633f7cf 100644 --- a/include/cln/rational_class.h +++ b/include/cln/rational_class.h @@ -22,10 +22,18 @@ public: cl_RA (const unsigned int); // argument must be < 2^29 cl_RA (const long); cl_RA (const unsigned long); +#ifdef HAVE_LONGLONG + cl_RA (const long long); + cl_RA (const unsigned long long); +#endif cl_RA& operator= (const int); // |argument| must be < 2^29 cl_RA& operator= (const unsigned int); // argument must be < 2^29 cl_RA& operator= (const long); cl_RA& operator= (const unsigned long); +#ifdef HAVE_LONGLONG + cl_RA& operator= (const long long); + cl_RA& operator= (const unsigned long long); +#endif // Other constructors. cl_RA (const char *); // Private constructor. @@ -53,6 +61,10 @@ CL_DEFINE_INT_CONSTRUCTORS(cl_RA) CL_DEFINE_INT_ASSIGNMENT_OPERATORS(cl_RA) CL_DEFINE_LONG_CONSTRUCTORS(cl_RA) CL_DEFINE_LONG_ASSIGNMENT_OPERATORS(cl_RA) +#ifdef HAVE_LONGLONG +CL_DEFINE_LONGLONG_CONSTRUCTORS(cl_RA) +CL_DEFINE_LONGLONG_ASSIGNMENT_OPERATORS(cl_RA) +#endif } // namespace cln diff --git a/include/cln/real.h b/include/cln/real.h index 5c2a2bd..d89cfef 100644 --- a/include/cln/real.h +++ b/include/cln/real.h @@ -84,6 +84,12 @@ inline const cl_R operator+ (const long x, const cl_R& y) { return cl_I(x) + y; } inline const cl_R operator+ (const unsigned long x, const cl_R& y) { return cl_I(x) + y; } +#ifdef HAVE_LONGLONG +inline const cl_R operator+ (const long long x, const cl_R& y) + { return cl_I(x) + y; } +inline const cl_R operator+ (const unsigned long long x, const cl_R& y) + { return cl_I(x) + y; } +#endif inline const cl_F operator+ (const float x, const cl_R& y) { return The(cl_F)(cl_R(x) + y); } inline const cl_F operator+ (const double x, const cl_R& y) @@ -96,6 +102,12 @@ inline const cl_R operator+ (const cl_R& x, const long y) { return x + cl_I(y); } inline const cl_R operator+ (const cl_R& x, const unsigned long y) { return x + cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_R operator+ (const cl_R& x, const long long y) + { return x + cl_I(y); } +inline const cl_R operator+ (const cl_R& x, const unsigned long long y) + { return x + cl_I(y); } +#endif inline const cl_F operator+ (const cl_R& x, const float y) { return The(cl_F)(x + cl_R(y)); } inline const cl_F operator+ (const cl_R& x, const double y) @@ -117,6 +129,12 @@ inline const cl_R operator- (const long x, const cl_R& y) { return cl_I(x) - y; } inline const cl_R operator- (const unsigned long x, const cl_R& y) { return cl_I(x) - y; } +#ifdef HAVE_LONGLONG +inline const cl_R operator- (const long long x, const cl_R& y) + { return cl_I(x) - y; } +inline const cl_R operator- (const unsigned long long x, const cl_R& y) + { return cl_I(x) - y; } +#endif inline const cl_F operator- (const float x, const cl_R& y) { return The(cl_F)(cl_R(x) - y); } inline const cl_F operator- (const double x, const cl_R& y) @@ -129,6 +147,12 @@ inline const cl_R operator- (const cl_R& x, const long y) { return x - cl_I(y); } inline const cl_R operator- (const cl_R& x, const unsigned long y) { return x - cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_R operator- (const cl_R& x, const long long y) + { return x - cl_I(y); } +inline const cl_R operator- (const cl_R& x, const unsigned long long y) + { return x - cl_I(y); } +#endif inline const cl_F operator- (const cl_R& x, const float y) { return The(cl_F)(x - cl_R(y)); } inline const cl_F operator- (const cl_R& x, const double y) @@ -145,6 +169,12 @@ inline const cl_R operator* (const long x, const cl_R& y) { return cl_I(x) * y; } inline const cl_R operator* (const unsigned long x, const cl_R& y) { return cl_I(x) * y; } +#ifdef HAVE_LONGLONG +inline const cl_R operator* (const long long x, const cl_R& y) + { return cl_I(x) * y; } +inline const cl_R operator* (const unsigned long long x, const cl_R& y) + { return cl_I(x) * y; } +#endif inline const cl_R operator* (const float x, const cl_R& y) { return cl_R(x) * y; } inline const cl_R operator* (const double x, const cl_R& y) @@ -157,6 +187,12 @@ inline const cl_R operator* (const cl_R& x, const long y) { return x * cl_I(y); } inline const cl_R operator* (const cl_R& x, const unsigned long y) { return x * cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_R operator* (const cl_R& x, const long long y) + { return x * cl_I(y); } +inline const cl_R operator* (const cl_R& x, const unsigned long long y) + { return x * cl_I(y); } +#endif inline const cl_R operator* (const cl_R& x, const float y) { return x * cl_R(y); } inline const cl_R operator* (const cl_R& x, const double y) @@ -179,6 +215,12 @@ inline const cl_R operator/ (const long x, const cl_R& y) { return cl_I(x) / y; } inline const cl_R operator/ (const unsigned long x, const cl_R& y) { return cl_I(x) / y; } +#ifdef HAVE_LONGLONG +inline const cl_R operator/ (const long long x, const cl_R& y) + { return cl_I(x) / y; } +inline const cl_R operator/ (const unsigned long long x, const cl_R& y) + { return cl_I(x) / y; } +#endif inline const cl_F operator/ (const float x, const cl_R& y) { return The(cl_F)(cl_R(x) / y); } inline const cl_F operator/ (const double x, const cl_R& y) @@ -191,6 +233,12 @@ inline const cl_R operator/ (const cl_R& x, const long y) { return x / cl_I(y); } inline const cl_R operator/ (const cl_R& x, const unsigned long y) { return x / cl_I(y); } +#ifdef HAVE_LONGLONG +inline const cl_R operator/ (const cl_R& x, const long long y) + { return x / cl_I(y); } +inline const cl_R operator/ (const cl_R& x, const unsigned long long y) + { return x / cl_I(y); } +#endif inline const cl_R operator/ (const cl_R& x, const float y) { return x / cl_R(y); } inline const cl_R operator/ (const cl_R& x, const double y) @@ -487,12 +535,20 @@ inline cl_R& operator+= (cl_R& x, const int y) { return x = x + y; } inline cl_R& operator+= (cl_R& x, const unsigned int y) { return x = x + y; } inline cl_R& operator+= (cl_R& x, const long y) { return x = x + y; } inline cl_R& operator+= (cl_R& x, const unsigned long y) { return x = x + y; } +#ifdef HAVE_LONGLONG +inline cl_R& operator+= (cl_R& x, const long long y) { return x = x + y; } +inline cl_R& operator+= (cl_R& x, const unsigned long long y) { return x = x + y; } +#endif inline cl_F& operator+= (cl_R& x, const float y) { return static_cast(x = x + y); } inline cl_F& operator+= (cl_R& x, const double y) { return static_cast(x = x + y); } inline cl_F& operator+= (cl_F& x, const int y) { return x = x + y; } inline cl_F& operator+= (cl_F& x, const unsigned int y) { return x = x + y; } inline cl_F& operator+= (cl_F& x, const long y) { return x = x + y; } inline cl_F& operator+= (cl_F& x, const unsigned long y) { return x = x + y; } +#ifdef HAVE_LONGLONG +inline cl_F& operator+= (cl_F& x, const long long y) { return x = x + y; } +inline cl_F& operator+= (cl_F& x, const unsigned long long y) { return x = x + y; } +#endif inline cl_R& operator++ /* prefix */ (cl_R& x) { return x = plus1(x); } inline void operator++ /* postfix */ (cl_R& x, int dummy) { (void)dummy; x = plus1(x); } inline cl_R& operator-= (cl_R& x, const cl_R& y) { return x = x - y; } @@ -503,12 +559,20 @@ inline cl_R& operator-= (cl_R& x, const int y) { return x = x - y; } inline cl_R& operator-= (cl_R& x, const unsigned int y) { return x = x - y; } inline cl_R& operator-= (cl_R& x, const long y) { return x = x - y; } inline cl_R& operator-= (cl_R& x, const unsigned long y) { return x = x - y; } +#ifdef HAVE_LONGLONG +inline cl_R& operator-= (cl_R& x, const long long y) { return x = x - y; } +inline cl_R& operator-= (cl_R& x, const unsigned long long y) { return x = x - y; } +#endif inline cl_F& operator-= (cl_R& x, const float y) { return static_cast(x = x - y); } inline cl_F& operator-= (cl_R& x, const double y) { return static_cast(x = x - y); } inline cl_F& operator-= (cl_F& x, const int y) { return x = x - y; } inline cl_F& operator-= (cl_F& x, const unsigned int y) { return x = x - y; } inline cl_F& operator-= (cl_F& x, const long y) { return x = x - y; } inline cl_F& operator-= (cl_F& x, const unsigned long y) { return x = x - y; } +#ifdef HAVE_LONGLONG +inline cl_F& operator-= (cl_F& x, const long long y) { return x = x - y; } +inline cl_F& operator-= (cl_F& x, const unsigned long long y) { return x = x - y; } +#endif inline cl_R& operator-- /* prefix */ (cl_R& x) { return x = minus1(x); } inline void operator-- /* postfix */ (cl_R& x, int dummy) { (void)dummy; x = minus1(x); } inline cl_R& operator*= (cl_R& x, const cl_R& y) { return x = x * y; } @@ -516,6 +580,10 @@ inline cl_R& operator*= (cl_R& x, const int y) { return x = x * y; } inline cl_R& operator*= (cl_R& x, const unsigned int y) { return x = x * y; } inline cl_R& operator*= (cl_R& x, const long y) { return x = x * y; } inline cl_R& operator*= (cl_R& x, const unsigned long y) { return x = x * y; } +#ifdef HAVE_LONGLONG +inline cl_R& operator*= (cl_R& x, const long long y) { return x = x * y; } +inline cl_R& operator*= (cl_R& x, const unsigned long long y) { return x = x * y; } +#endif inline cl_R& operator*= (cl_R& x, const float y) { return x = x * y; } inline cl_R& operator*= (cl_R& x, const double y) { return x = x * y; } inline cl_R& operator/= (cl_R& x, const cl_R& y) { return x = x / y; } @@ -526,12 +594,20 @@ inline cl_R& operator/= (cl_R& x, const int y) { return x = x / y; } inline cl_R& operator/= (cl_R& x, const unsigned int y) { return x = x / y; } inline cl_R& operator/= (cl_R& x, const long y) { return x = x / y; } inline cl_R& operator/= (cl_R& x, const unsigned long y) { return x = x / y; } +#ifdef HAVE_LONGLONG +inline cl_R& operator/= (cl_R& x, const long long y) { return x = x / y; } +inline cl_R& operator/= (cl_R& x, const unsigned long long y) { return x = x / y; } +#endif inline cl_R& operator/= (cl_R& x, const float y) { return x = x / y; } inline cl_R& operator/= (cl_R& x, const double y) { return x = x / y; } inline cl_F& operator/= (cl_F& x, const int y) { return x = x / y; } inline cl_F& operator/= (cl_F& x, const unsigned int y) { return x = x / y; } inline cl_F& operator/= (cl_F& x, const long y) { return x = x / y; } inline cl_F& operator/= (cl_F& x, const unsigned long y) { return x = x / y; } +#ifdef HAVE_LONGLONG +inline cl_F& operator/= (cl_F& x, const long long y) { return x = x / y; } +inline cl_F& operator/= (cl_F& x, const unsigned long long y) { return x = x / y; } +#endif #endif diff --git a/include/cln/real_class.h b/include/cln/real_class.h index bae1e51..44f25e3 100644 --- a/include/cln/real_class.h +++ b/include/cln/real_class.h @@ -22,6 +22,10 @@ public: cl_R (const unsigned int); // argument must be < 2^29 cl_R (const long); cl_R (const unsigned long); +#ifdef HAVE_LONGLONG + cl_R (const long long); + cl_R (const unsigned long long); +#endif cl_R (const float); cl_R (const double); cl_R& operator= (const int); // |argument| must be < 2^29 @@ -30,6 +34,10 @@ public: cl_R& operator= (const unsigned long); cl_R& operator= (const float); cl_R& operator= (const double); +#ifdef HAVE_LONGLONG + cl_R& operator= (const long long); + cl_R& operator= (const unsigned long long); +#endif // Other constructors. cl_R (const char *); // Private constructor. @@ -56,6 +64,10 @@ CL_DEFINE_INT_CONSTRUCTORS(cl_R) CL_DEFINE_INT_ASSIGNMENT_OPERATORS(cl_R) CL_DEFINE_LONG_CONSTRUCTORS(cl_R) CL_DEFINE_LONG_ASSIGNMENT_OPERATORS(cl_R) +#ifdef HAVE_LONGLONG +CL_DEFINE_LONGLONG_CONSTRUCTORS(cl_R) +CL_DEFINE_LONGLONG_ASSIGNMENT_OPERATORS(cl_R) +#endif CL_DEFINE_FLOAT_CONSTRUCTOR(cl_R) CL_DEFINE_DOUBLE_CONSTRUCTOR(cl_R) diff --git a/include/cln/sfloat.h b/include/cln/sfloat.h index a95d65f..4ca1bf6 100644 --- a/include/cln/sfloat.h +++ b/include/cln/sfloat.h @@ -223,7 +223,7 @@ extern const decoded_sfloat decode_float (const cl_SF& x); // den Exponenten von (decode-float x). // x = 0.0 liefert 0. // x = (-1)^s * 2^e * m liefert e. -extern sintL float_exponent (const cl_SF& x); +extern sintE float_exponent (const cl_SF& x); // float_radix(x) liefert (float-radix x), wo x ein Float ist. inline sintL float_radix (const cl_SF& x) diff --git a/include/cln/types.h b/include/cln/types.h index d23fcc7..067411e 100644 --- a/include/cln/types.h +++ b/include/cln/types.h @@ -96,6 +96,18 @@ typedef unsigned int uintC; #endif +// Integer type used for lfloat exponents. +// Constraint: sizeof(uintE) >= sizeof(uintC) +#if (defined(HAVE_LONGLONG) && (defined(__alpha__) || defined(__ia64__) || defined(__powerpc64__) || defined(__x86_64__) || defined(__i386__))) + #define intEsize 64 + typedef sint64 sintE; + typedef uint64 uintE; + #else + #define intEsize 32 + typedef sint32 sintE; + typedef uint32 uintE; + #endif + // Integer type as large as a pointer. // Assumption: sizeof(long) == sizeof(void*) #define intPsize long_bitsize diff --git a/src/base/cl_low.h b/src/base/cl_low.h index 98b371e..46485b9 100644 --- a/src/base/cl_low.h +++ b/src/base/cl_low.h @@ -77,7 +77,7 @@ inline uint32 highlow32_0 (uint16 high) return (uint32)high << 16; } -#ifdef HAVE_FAST_LONGLONG +#ifdef HAVE_LONGLONG // High-Word einer 64-Bit-Zahl bestimmen // high32(wert) @@ -107,7 +107,7 @@ inline uint64 highlow64_0 (uint32 high) return (uint64)high << 32; } -#endif /* HAVE_FAST_LONGLONG */ +#endif /* HAVE_LONGLONG */ // Multipliziert zwei 16-Bit-Zahlen miteinander und liefert eine 32-Bit-Zahl: @@ -1211,6 +1211,27 @@ inline uint32 mulu32_unchecked (uint32 arg1, uint32 arg2) // < uintL ergebnis : Wurzel, >=0, <2^16 extern uintL isqrt (uintL x); +#ifdef HAVE_LONGLONG +// Extracts integer root of a 64-bit number and returns a 32-bit number. +// isqrt(x) +// > uintQ x : radicand, >=0, <2^64 +// < uintL result : square root, >=0, <2^32 + extern uintL isqrt (uintQ x); +#endif + +// Sorry for this. We need an isqrt function taking uintC arguments but we +// cannot use overloading since this would lead to ambiguities with any of the +// two signatures above. + inline uintL isqrtC (uintC x) + { +#if (intCsize==32) + return isqrt((uintL)x); +#else + return isqrt((uintQ)x); +#endif + } + + // Zieht die Ganzzahl-Wurzel aus einer 64-Bit-Zahl und // liefert eine 32-Bit-Wurzel. // isqrt(x1,x0) diff --git a/src/base/cl_macros.h b/src/base/cl_macros.h index 8c160e9..ec49995 100644 --- a/src/base/cl_macros.h +++ b/src/base/cl_macros.h @@ -152,7 +152,7 @@ namespace cln { // Return 2^n, n a constant expression. // Same as bit(n), but undefined if n<0 or n>={long_}long_bitsize. - #ifdef HAVE_FAST_LONGLONG + #if defined(HAVE_FAST_LONGLONG) || defined(intQsize) #define bitc(n) (1ULL << (((n) >= 0 && (n) < long_long_bitsize) ? (n) : 0)) #else #define bitc(n) (1UL << (((n) >= 0 && (n) < long_bitsize) ? (n) : 0)) diff --git a/src/base/low/cl_low_isqrt.cc b/src/base/low/cl_low_isqrt.cc index b4cf1cc..2588c0a 100644 --- a/src/base/low/cl_low_isqrt.cc +++ b/src/base/low/cl_low_isqrt.cc @@ -61,4 +61,39 @@ uintL isqrt (uintL x) }} } +#ifdef HAVE_LONGLONG +uintL isqrt (uintQ x) +{ + // As isqrt (uintL) above, but with 64-bit numbers. + if (x==0) return 0; // x=0 -> y=0 + { var uintC k2; integerlength64(x,k2=); // 2^(k2-1) <= x < 2^k2 + {var uintC k1 = floor(k2-1,2); // k1 = k-1, k as above + if (k1 < 32-1) + // k < 32 + { var uintL y = (x >> (k1+2)) | bit(k1); // always 2^(k-1) <= y < 2^k + loop + { var uintL z; + divu_6432_3232(high32(x),low32(x),y, z=,); // z := x/y (works, since x/y < 2^(2k)/2^(k-1) = 2^(k+1) <= 2^32) + if (z >= y) break; + y = floor(z+y,2); // geht, da z+y < 2*y < 2^(k+1) <= 2^32 + } + return y; + } + else + // k = 32, careful! + { var uintL x1 = high32(x); + var uintL y = (x >> (32+1)) | bit(32-1); // stets 2^(k-1) <= y < 2^k + loop + { var uintL z; + if (x1 >= y) break; // division x/y would overflow -> z > y + divu_6432_3232(high32(x),low32(x),y, z=,); // divide x/y + if (z >= y) break; + y = floor(z+y,2); + } + return y; + } + }} +} +#endif + } // namespace cln diff --git a/src/complex/algebraic/cl_LF_hypot.cc b/src/complex/algebraic/cl_LF_hypot.cc index 70f51a0..cec82df 100644 --- a/src/complex/algebraic/cl_LF_hypot.cc +++ b/src/complex/algebraic/cl_LF_hypot.cc @@ -46,29 +46,29 @@ const cl_LF cl_hypot (const cl_LF& a, const cl_LF& b) a = shorten(a,b_len); } } - var sintL a_exp; - var sintL b_exp; + var sintE a_exp; + var sintE b_exp; { // Exponenten von a holen: - var uintL uexp = TheLfloat(a)->expo; + var uintE uexp = TheLfloat(a)->expo; if (uexp == 0) // a=0.0 -> liefere (abs b) : return (minusp(b) ? -b : b); - a_exp = (sintL)(uexp - LF_exp_mid); + a_exp = (sintE)(uexp - LF_exp_mid); } { // Exponenten von b holen: - var uintL uexp = TheLfloat(b)->expo; + var uintE uexp = TheLfloat(b)->expo; if (uexp == 0) // b=0.0 -> liefere (abs a) : return (minusp(a) ? -a : a); - b_exp = (sintL)(uexp - LF_exp_mid); + b_exp = (sintE)(uexp - LF_exp_mid); } // Nun a_exp = float_exponent(a), b_exp = float_exponent(b). - var sintL e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten + var sintE e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten // a und b durch 2^e dividieren: - var cl_LF na = ((b_exp > a_exp) && ((uintL)(b_exp-a_exp) > (uintL)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(a)->len) : scale_float(a,-e)); - var cl_LF nb = ((a_exp > b_exp) && ((uintL)(a_exp-b_exp) > (uintL)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(b)->len) : scale_float(b,-e)); + var cl_LF na = ((b_exp > a_exp) && ((uintE)(b_exp-a_exp) > (uintE)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(a)->len) : scale_float(a,-e)); + var cl_LF nb = ((a_exp > b_exp) && ((uintE)(a_exp-b_exp) > (uintE)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(b)->len) : scale_float(b,-e)); // c' := a'*a'+b'*b' berechnen: var cl_LF nc = square(na) + square(nb); return scale_float(sqrt(nc),e); // c' := sqrt(c'), 2^e*c' diff --git a/src/complex/elem/division/cl_C_LF_recip.cc b/src/complex/elem/division/cl_C_LF_recip.cc index c62ef08..0bf2d7f 100644 --- a/src/complex/elem/division/cl_C_LF_recip.cc +++ b/src/complex/elem/division/cl_C_LF_recip.cc @@ -43,29 +43,29 @@ const cl_C_LF cl_C_recip (const cl_LF& a, const cl_LF& b) a = shorten(a,b_len); } } - var sintL a_exp; - var sintL b_exp; + var sintE a_exp; + var sintE b_exp; { // Exponenten von a holen: - var uintL uexp = TheLfloat(a)->expo; + var uintE uexp = TheLfloat(a)->expo; if (uexp == 0) // a=0.0 -> liefere (complex a (- (/ b))) : return cl_C_LF(a,-recip(b)); - a_exp = (sintL)(uexp - LF_exp_mid); + a_exp = (sintE)(uexp - LF_exp_mid); } { // Exponenten von b holen: - var uintL uexp = TheLfloat(b)->expo; + var uintE uexp = TheLfloat(b)->expo; if (uexp == 0) // b=0.0 -> liefere (complex (/ a) b) : return cl_C_LF(recip(a),b); - b_exp = (sintL)(uexp - LF_exp_mid); + b_exp = (sintE)(uexp - LF_exp_mid); } // Nun a_exp = float_exponent(a), b_exp = float_exponent(b). - var sintL e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten + var sintE e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten // a und b durch 2^e dividieren: - var cl_LF na = ((b_exp > a_exp) && ((uintL)(b_exp-a_exp) > (uintL)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(a)->len) : scale_float(a,-e)); - var cl_LF nb = ((a_exp > b_exp) && ((uintL)(a_exp-b_exp) > (uintL)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(b)->len) : scale_float(b,-e)); + var cl_LF na = ((b_exp > a_exp) && ((uintE)(b_exp-a_exp) > (uintE)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(a)->len) : scale_float(a,-e)); + var cl_LF nb = ((a_exp > b_exp) && ((uintE)(a_exp-b_exp) > (uintE)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(b)->len) : scale_float(b,-e)); // c' := a'*a'+b'*b' berechnen: var cl_LF nc = square(na) + square(nb); // 2^(-e)*a'/c' + i * -2^(-e)*b'/c' diff --git a/src/float/cl_F.h b/src/float/cl_F.h index 9dd7b20..c140642 100644 --- a/src/float/cl_F.h +++ b/src/float/cl_F.h @@ -24,7 +24,7 @@ nonreturning_function(extern, cl_error_floating_point_underflow, (void)); // SF 1 8 16 // FF 1 8 23 // DF 1 11 52 -// LF 1 32 intDsize*n >= 53 +// LF 1 32 or 64 intDsize*n >= 53 // Konversionen ohne Rundung: diff --git a/src/float/dfloat/misc/cl_DF_exponent.cc b/src/float/dfloat/misc/cl_DF_exponent.cc index db94aaa..c4058e3 100644 --- a/src/float/dfloat/misc/cl_DF_exponent.cc +++ b/src/float/dfloat/misc/cl_DF_exponent.cc @@ -14,7 +14,7 @@ namespace cln { MAYBE_INLINE -sintL float_exponent (const cl_DF& x) +sintE float_exponent (const cl_DF& x) { var uintL uexp = DF_uexp(TheDfloat(x)->dfloat_value_semhi); if (uexp==0) { return 0; } diff --git a/src/float/ffloat/misc/cl_FF_exponent.cc b/src/float/ffloat/misc/cl_FF_exponent.cc index 72f3cee..f9f6b67 100644 --- a/src/float/ffloat/misc/cl_FF_exponent.cc +++ b/src/float/ffloat/misc/cl_FF_exponent.cc @@ -14,7 +14,7 @@ namespace cln { MAYBE_INLINE -sintL float_exponent (const cl_FF& x) +sintE float_exponent (const cl_FF& x) { var uintL uexp = FF_uexp(cl_ffloat_value(x)); if (uexp==0) { return 0; } diff --git a/src/float/input/cl_F_read.cc b/src/float/input/cl_F_read.cc index 3ed96bf..627652d 100644 --- a/src/float/input/cl_F_read.cc +++ b/src/float/input/cl_F_read.cc @@ -131,7 +131,7 @@ const cl_F read_float (const cl_read_flags& flags, const char * string, const ch ptr_after_prec = skip_digits(ptr,string_limit,10); if (ptr_after_prec == ptr) goto not_float_syntax; var cl_I prec1 = digits_to_I(ptr,ptr_after_prec-ptr,10); - var uintC prec2 = cl_I_to_UL(prec1); + var uintC prec2 = cl_I_to_UC(prec1); prec = (float_base==10 ? float_format(prec2) : (float_format_t)((uintC)((1+prec2)*::log((double)float_base)*1.442695041)+1) ); @@ -155,7 +155,7 @@ const cl_F read_float (const cl_read_flags& flags, const char * string, const ch (float_base==10 ? float_format(prec2) : (float_format_t)((uintC)((1+prec2)*::log((double)float_base)*1.442695041)+1) ); - if ((uintC)precx > (uintC)prec) + if ((uintE)precx > (uintE)prec) prec = precx; } } diff --git a/src/float/lfloat/algebraic/cl_LF_sqrt.cc b/src/float/lfloat/algebraic/cl_LF_sqrt.cc index 87d974c..8d53c6a 100644 --- a/src/float/lfloat/algebraic/cl_LF_sqrt.cc +++ b/src/float/lfloat/algebraic/cl_LF_sqrt.cc @@ -35,7 +35,7 @@ const cl_LF sqrt (const cl_LF& x) // sonst aufrunden. // Bei rounding overflow Mantisse um 1 Bit nach rechts schieben // und Exponent incrementieren. - var uintL uexp = TheLfloat(x)->expo; + var uintE uexp = TheLfloat(x)->expo; if (uexp==0) { return x; } // x=0.0 -> 0.0 als Ergebnis var uintC len = TheLfloat(x)->len; // Radikanden bilden: @@ -59,7 +59,7 @@ const cl_LF sqrt (const cl_LF& x) clear_loop_msp(ptr,len+1); // n+1 Nulldigits anh�gen } // Compute ((uexp - LF_exp_mid + 1) >> 1) + LF_exp_mid without risking - // uintL overflow. + // uintE overflow. uexp = ((uexp - ((LF_exp_mid - 1) & 1)) >> 1) - ((LF_exp_mid - 1) >> 1) + LF_exp_mid; // Ergebnis allozieren: diff --git a/src/float/lfloat/cl_LF.h b/src/float/lfloat/cl_LF.h index 8da8d40..df93d75 100644 --- a/src/float/lfloat/cl_LF.h +++ b/src/float/lfloat/cl_LF.h @@ -12,7 +12,7 @@ namespace cln { struct cl_heap_lfloat : cl_heap { unsigned int len; // length of mantissa (in digits) int sign; // sign (0 or -1) - uint32 expo; // exponent + uintE expo; // exponent uintD data[1]; // mantissa }; @@ -20,11 +20,15 @@ struct cl_heap_lfloat : cl_heap { // so that a LF has not fewer mantissa bits than a DF. #define LF_minlen ceiling(53,intDsize) // Exponent. -// Define as 'unsigned int', not 'unsigned long', so that -// LF_exp_high+1 wraps around to 0 just like the 'expo' field does. +#if (intEsize==64) + #define LF_exp_low 1 + #define LF_exp_mid 0x8000000000000000ULL + #define LF_exp_high 0xFFFFFFFFFFFFFFFFULL +#else #define LF_exp_low 1 #define LF_exp_mid 0x80000000U #define LF_exp_high 0xFFFFFFFFU +#endif inline cl_heap_lfloat* TheLfloat (cl_heap_lfloat* p) { return p; } diff --git a/src/float/lfloat/cl_LF_impl.h b/src/float/lfloat/cl_LF_impl.h index 7194fc9..617a5ba 100644 --- a/src/float/lfloat/cl_LF_impl.h +++ b/src/float/lfloat/cl_LF_impl.h @@ -16,10 +16,10 @@ extern cl_class cl_class_lfloat; // Builds a long-float, without filling the mantissa. // allocate_lfloat(len,expo,sign) // > uintC len: length of mantissa (in digits) -// > uint32 expo: exponent +// > uintE expo: exponent // > cl_signean sign: sign (0 = +, -1 = -) // The long-float is only complete when the mantissa has been filled in! -inline cl_heap_lfloat* allocate_lfloat (uintC len, uint32 expo, cl_signean sign) +inline cl_heap_lfloat* allocate_lfloat (uintC len, uintE expo, cl_signean sign) { cl_heap_lfloat* p = (cl_heap_lfloat*) malloc_hook(offsetofa(cl_heap_lfloat,data)+sizeof(uintD)*len); p->refcount = 1; @@ -60,17 +60,17 @@ inline cl_LF::cl_LF (cl_heap_lfloat* ptr) : cl_F ((cl_private_thing) ptr) {} // zerlegt ein Long-Float obj. // Ist obj=0.0, wird zero_statement ausgeführt. // Sonst: cl_signean sign = Vorzeichen (0 = +, -1 = -), -// sintL exp = Exponent (vorzeichenbehaftet), +// sintE exp = Exponent (vorzeichenbehaftet), // UDS mantMSDptr/mantlen/mantLSDptr = Mantisse // (>= 2^(intDsize*mantlen-1), < 2^(intDsize*mantlen)), // mit mantlen>=LF_minlen. #define LF_decode(obj, zero_statement, sign_zuweisung,exp_zuweisung,mantMSDptr_zuweisung,mantlen_zuweisung,mantLSDptr_zuweisung) \ { var Lfloat _x = TheLfloat(obj); \ - var uintL uexp = _x->expo; \ + var uintE uexp = _x->expo; \ if (uexp==0) \ { mantlen_zuweisung _x->len; zero_statement } /* e=0 -> Zahl 0.0 */\ else \ - { exp_zuweisung (sintL)(uexp - LF_exp_mid); /* Exponent */ \ + { exp_zuweisung (sintE)(uexp - LF_exp_mid); /* Exponent */ \ sign_zuweisung _x->sign; /* Vorzeichen */\ unused (mantMSDptr_zuweisung arrayMSDptr(_x->data, (uintP)(mantlen_zuweisung _x->len))); /* Mantissen-UDS */\ unused (mantLSDptr_zuweisung arrayLSDptr(_x->data, (uintP)(mantlen_zuweisung _x->len))); \ @@ -112,12 +112,12 @@ inline const cl_LF encode_LF1 (uintC len) // Einpacken eines Long-Float: // encode_LFu(sign,uexp,mantMSDptr,mantlen) liefert ein Long-Float // > cl_signean sign: Vorzeichen -// > uintL exp: Exponent + LF_exp_mid +// > uintE exp: Exponent + LF_exp_mid // > uintD* mantMSDptr: Pointer auf eine NUDS mit gesetztem höchstem Bit // > uintC mantlen: Anzahl der Digits, >= LF_minlen // < cl_LF erg: neues Long-Float mit der UDS mantMSDptr/mantlen/.. als Mantisse // Der Exponent wird nicht auf Überlauf/Unterlauf getestet. -inline const cl_LF encode_LFu (cl_signean sign, uintL uexp, const uintD* mantMSDptr, uintC mantlen) +inline const cl_LF encode_LFu (cl_signean sign, uintE uexp, const uintD* mantMSDptr, uintC mantlen) { var Lfloat erg = allocate_lfloat(mantlen,uexp,sign); /* Exponent */ copy_loop_msp(mantMSDptr,arrayMSDptr(TheLfloat(erg)->data,mantlen),mantlen); /* Mantisse übertragen */ @@ -127,20 +127,20 @@ inline const cl_LF encode_LFu (cl_signean sign, uintL uexp, const uintD* mantMSD // Einpacken eines Long-Float: // encode_LF(sign,exp,mantMSDptr,mantlen) liefert ein Long-Float // > cl_signean sign: Vorzeichen -// > sintL exp: Exponent +// > sintE exp: Exponent // > uintD* mantMSDptr: Pointer auf eine NUDS mit gesetztem höchstem Bit // > uintC mantlen: Anzahl der Digits, >= LF_minlen // < cl_LF erg: neues Long-Float mit der UDS mantMSDptr/mantlen/.. als Mantisse // Der Exponent wird nicht auf Überlauf/Unterlauf getestet. -inline const cl_LF encode_LF (cl_signean sign, sintL exp, const uintD* mantMSDptr, uintC mantlen) +inline const cl_LF encode_LF (cl_signean sign, sintE exp, const uintD* mantMSDptr, uintC mantlen) { - return encode_LFu(sign,LF_exp_mid+(uintL)exp,mantMSDptr,mantlen); + return encode_LFu(sign,LF_exp_mid+(uintE)exp,mantMSDptr,mantlen); } // Einpacken eines Long-Float: // encode_LF_array(sign,exp,mantarr,mantlen) liefert ein Long-Float // > cl_signean sign: Vorzeichen -// > sintL exp: Exponent +// > sintE exp: Exponent // > uintD mantarr[]: NUDS mit gesetztem höchstem Bit // > uintC mantlen: Anzahl der Digits, >= LF_minlen // < cl_LF erg: neues Long-Float mit der UDS mantarr[] als Mantisse diff --git a/src/float/lfloat/elem/cl_LF_1plus.cc b/src/float/lfloat/elem/cl_LF_1plus.cc index 00fcfe6..2114521 100644 --- a/src/float/lfloat/elem/cl_LF_1plus.cc +++ b/src/float/lfloat/elem/cl_LF_1plus.cc @@ -35,15 +35,15 @@ const cl_LF LF_LF_plus_LF (const cl_LF& arg1, const cl_LF& arg2) // Normalisiere, fertig. var cl_LF x1 = arg1; var cl_LF x2 = arg2; - var uintL uexp1 = TheLfloat(arg1)->expo; - var uintL uexp2 = TheLfloat(arg2)->expo; + var uintE uexp1 = TheLfloat(arg1)->expo; + var uintE uexp2 = TheLfloat(arg2)->expo; if (uexp1 < uexp2) // x1 und x2 vertauschen - { x1 = arg2; x2 = arg1; swap(uintL, uexp1,uexp2); } + { x1 = arg2; x2 = arg1; swap(uintE, uexp1,uexp2); } // uexp1 >= uexp2 if (uexp2==0) { return x1; } // x2=0.0 -> x1 als Ergebnis var uintC len = TheLfloat(x1)->len; // Länge n von x1 und x2 - var uintL expdiff = uexp1-uexp2; // e1-e2 + var uintE expdiff = uexp1-uexp2; // e1-e2 if ((expdiff == 0) && (TheLfloat(x1)->sign != TheLfloat(x2)->sign)) // verschiedene Vorzeichen, aber gleicher Exponent { // Vorzeichen des Ergebnisses festlegen: @@ -54,7 +54,7 @@ const cl_LF LF_LF_plus_LF (const cl_LF& arg1, const cl_LF& arg2) if (erg<0) // |x1| < |x2| // x1 und x2 vertauschen, expdiff bleibt =0 { x1.pointer = arg2.pointer; x2.pointer = arg1.pointer; - swap(uintL, uexp1,uexp2); + swap(uintE, uexp1,uexp2); } } if (expdiff >= intDsize * len + 2) // e1-e2 >= 16n+2 ? @@ -172,7 +172,7 @@ const cl_LF LF_LF_plus_LF (const cl_LF& arg1, const cl_LF& arg2) rounding_bits = 0; // und keine weiteren Rundungsbits // Exponenten um intDsize*k erniedrigen: k = intDsize*k; - {var uintL uexp = TheLfloat(y)->expo; + {var uintE uexp = TheLfloat(y)->expo; #if !(LF_exp_low==1) if (uexp < k+LF_exp_low) #else @@ -205,7 +205,7 @@ const cl_LF LF_LF_plus_LF (const cl_LF& arg1, const cl_LF& arg2) rounding_bits = 0; // = rounding_bits << s; } // Exponenten um s erniedrigen: - {var uintL uexp = TheLfloat(y)->expo; + {var uintE uexp = TheLfloat(y)->expo; #if !(LF_exp_low==1) if (uexp < s+LF_exp_low) #else diff --git a/src/float/lfloat/elem/cl_LF_I_div.cc b/src/float/lfloat/elem/cl_LF_I_div.cc index 188f089..139f3e1 100644 --- a/src/float/lfloat/elem/cl_LF_I_div.cc +++ b/src/float/lfloat/elem/cl_LF_I_div.cc @@ -84,8 +84,8 @@ const cl_LF cl_LF_I_div (const cl_LF& x, const cl_I& y) } // Quotient MSDptr/len/.. ist nun normalisiert: höchstes Bit =1. // exponent := exponent(x) - intDsize*y_len + shiftcount - var uintL uexp = TheLfloat(x)->expo; - var uintL dexp = intDsize*y_len - shiftcount; // >= 0 ! + var uintE uexp = TheLfloat(x)->expo; + var uintE dexp = intDsize*y_len - shiftcount; // >= 0 ! if ((uexp < dexp) || ((uexp = uexp - dexp) < LF_exp_low)) { if (underflow_allowed()) { cl_error_floating_point_underflow(); } diff --git a/src/float/lfloat/elem/cl_LF_I_mul.cc b/src/float/lfloat/elem/cl_LF_I_mul.cc index 927c268..d83424a 100644 --- a/src/float/lfloat/elem/cl_LF_I_mul.cc +++ b/src/float/lfloat/elem/cl_LF_I_mul.cc @@ -60,8 +60,8 @@ const cl_R cl_LF_I_mul (const cl_LF& x, const cl_I& y) } // Produkt ist nun normalisiert: höchstes Bit =1. // exponent := exponent(x) + intDsize*y_len - shiftcount - var uintL uexp = TheLfloat(x)->expo; - var uintL iexp = intDsize*y_len - shiftcount; // >= 0 ! + var uintE uexp = TheLfloat(x)->expo; + var uintE iexp = intDsize*y_len - shiftcount; // >= 0 ! uexp = uexp + iexp; if ((uexp < iexp) || (uexp > LF_exp_high)) cl_error_floating_point_overflow(); diff --git a/src/float/lfloat/elem/cl_LF_compare.cc b/src/float/lfloat/elem/cl_LF_compare.cc index 797f403..9d64f59 100644 --- a/src/float/lfloat/elem/cl_LF_compare.cc +++ b/src/float/lfloat/elem/cl_LF_compare.cc @@ -32,8 +32,8 @@ cl_signean compare (const cl_LF& x, const cl_LF& y) { if (!minusp(x)) // y>=0, x>=0 { // Vergleiche Exponenten und Mantissen: - { var uintL x_uexp = TheLfloat(x)->expo; - var uintL y_uexp = TheLfloat(y)->expo; + { var uintE x_uexp = TheLfloat(x)->expo; + var uintE y_uexp = TheLfloat(y)->expo; if (x_uexp < y_uexp) return signean_minus; // x y_uexp) return signean_plus; // x>y } @@ -72,8 +72,8 @@ cl_signean compare (const cl_LF& x, const cl_LF& y) else // y<0, x<0 { // Vergleiche Exponenten und Mantissen: - { var uintL x_uexp = TheLfloat(x)->expo; - var uintL y_uexp = TheLfloat(y)->expo; + { var uintE x_uexp = TheLfloat(x)->expo; + var uintE y_uexp = TheLfloat(y)->expo; if (x_uexp < y_uexp) return signean_plus; // |x|<|y| -> x>y if (x_uexp > y_uexp) return signean_minus; // |x|>|y| -> xlen; var uintC len2 = TheLfloat(x2)->len; var uintC len = (len1 < len2 ? len1 : len2); // min. Länge n von x1 und x2 - var uintL uexp2 = TheLfloat(x2)->expo; + var uintE uexp2 = TheLfloat(x2)->expo; if (uexp2==0) { cl_error_division_by_0(); } // x2=0.0 -> Error - var uintL uexp1 = TheLfloat(x1)->expo; + var uintE uexp1 = TheLfloat(x1)->expo; if (uexp1==0) // x1=0.0 -> Ergebnis 0.0 { if (len < len1) return shorten(x1,len); else return x1; } // Exponenten subtrahieren: @@ -55,7 +55,7 @@ const cl_LF operator/ (const cl_LF& x1, const cl_LF& x2) } else { uexp1 = uexp1 - uexp2; // Carry - if (uexp1 < (uintL)(LF_exp_low-1-LF_exp_mid)) + if (uexp1 < (uintE)(LF_exp_low-1-LF_exp_mid)) { if (underflow_allowed()) { cl_error_floating_point_underflow(); } else diff --git a/src/float/lfloat/elem/cl_LF_from_I.cc b/src/float/lfloat/elem/cl_LF_from_I.cc index bb9c483..c18e329 100644 --- a/src/float/lfloat/elem/cl_LF_from_I.cc +++ b/src/float/lfloat/elem/cl_LF_from_I.cc @@ -39,11 +39,11 @@ const cl_LF cl_I_to_LF (const cl_I& x, uintC len) var uintC exp = integer_length(abs_x); // (integer-length x) < intDsize*2^intCsize // Teste, ob exp <= LF_exp_high-LF_exp_mid : if ( (log2_intDsize+intCsize < 32) - && ((uintL)(intDsize*bitc(intCsize)-1) <= (uintL)(LF_exp_high-LF_exp_mid)) + && ((uintE)(intDsize*bitc(intCsize)-1) <= (uintE)(LF_exp_high-LF_exp_mid)) ) {} // garantiert exp <= intDsize*2^intCsize-1 <= LF_exp_high-LF_exp_mid else - { if (!(exp <= (uintL)(LF_exp_high-LF_exp_mid))) { cl_error_floating_point_overflow(); } } + { if (!(exp <= (uintE)(LF_exp_high-LF_exp_mid))) { cl_error_floating_point_overflow(); } } // Long-Float bauen: var Lfloat y = allocate_lfloat(len,exp+LF_exp_mid,sign); var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len); @@ -91,7 +91,7 @@ const cl_LF cl_I_to_LF (const cl_I& x, uintC len) { mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0 // Exponenten incrementieren: if ( (log2_intDsize+intCsize < 32) - && ((uintL)(intDsize*bitc(intCsize)-1) < (uintL)(LF_exp_high-LF_exp_mid)) + && ((uintE)(intDsize*bitc(intCsize)-1) < (uintE)(LF_exp_high-LF_exp_mid)) ) // garantiert exp < intDsize*2^intCsize-1 <= LF_exp_high-LF_exp_mid { (TheLfloat(y)->expo)++; } // jetzt exp <= LF_exp_high-LF_exp_mid diff --git a/src/float/lfloat/elem/cl_LF_fround.cc b/src/float/lfloat/elem/cl_LF_fround.cc index 8decf62..979a2e7 100644 --- a/src/float/lfloat/elem/cl_LF_fround.cc +++ b/src/float/lfloat/elem/cl_LF_fround.cc @@ -24,18 +24,18 @@ const cl_LF fround (const cl_LF& x) // e>=16n -> Ergebnis x #if 0 var cl_signean sign; - var sintL exp; + var sintE exp; var const uintD* mantMSDptr; var uintC mantlen; LF_decode(x, { return x; }, sign=,exp=,mantMSDptr=,mantlen=,); if (exp<0) { return encode_LF0(mantlen); } // e<0 -> Ergebnis 0.0 - if ((uintL)exp >= intDsize*mantlen) // e>=16n -> x als Ergebnis + if ((uintE)exp >= intDsize*mantlen) // e>=16n -> x als Ergebnis { return x; } else // 0 <= e < 16n { // alle hinteren 16n-e Bits wegrunden: - var uintC count = floor((uintL)exp,intDsize); // zu kopierende Digits, < mantlen - var uintC bitcount = ((uintL)exp) % intDsize; // zu kopierende Bits danach, >=0, =0, abrunden @@ -75,12 +75,12 @@ const cl_LF fround (const cl_LF& x) } #else var uintC len = TheLfloat(x)->len; - var uintL uexp = TheLfloat(x)->expo; + var uintE uexp = TheLfloat(x)->expo; if (uexp < LF_exp_mid) { if (uexp == 0) { return x; } // x=0.0 -> Ergebnis 0.0 return encode_LF0(len); // e<0 -> Ergebnis 0.0 } - var uintL exp = uexp - LF_exp_mid; + var uintE exp = uexp - LF_exp_mid; if (exp >= intDsize*len) // e>=16n -> x als Ergebnis { return x; } // 0 <= e < 16n diff --git a/src/float/lfloat/elem/cl_LF_ftrunc.cc b/src/float/lfloat/elem/cl_LF_ftrunc.cc index 764d732..f0f204e 100644 --- a/src/float/lfloat/elem/cl_LF_ftrunc.cc +++ b/src/float/lfloat/elem/cl_LF_ftrunc.cc @@ -24,12 +24,12 @@ const cl_LF ftruncate (const cl_LF& x) // e>=16n -> Ergebnis x #if 0 var cl_signean sign; - var sintL exp; + var sintE exp; var const uintD* mantMSDptr; var uintC mantlen; LF_decode(x, { return x; }, sign=,exp=,mantMSDptr=,mantlen=,); if (exp<=0) { return encode_LF0(mantlen); } // e<=0 -> Ergebnis 0.0 - if ((uintL)exp >= intDsize*mantlen) // e>=16n -> x als Ergebnis + if ((uintE)exp >= intDsize*mantlen) // e>=16n -> x als Ergebnis { return x; } else // 0 < e < 16n @@ -37,8 +37,8 @@ const cl_LF ftruncate (const cl_LF& x) { CL_ALLOCA_STACK; var uintD* MSDptr; num_stack_alloc(mantlen, MSDptr=,); - { var uintC count = floor((uintL)exp,intDsize); // zu kopierende Digits, < mantlen - var uintC bitcount = ((uintL)exp) % intDsize; // zu kopierende Bits danach, >=0, =0, len; - var uintL uexp = TheLfloat(x)->expo; + var uintE uexp = TheLfloat(x)->expo; if (uexp <= LF_exp_mid) { if (uexp == 0) { return x; } // x=0.0 -> Ergebnis 0.0 return encode_LF0(len); // e<=0 -> Ergebnis 0.0 } - var uintL exp = uexp - LF_exp_mid; + var uintE exp = uexp - LF_exp_mid; if (exp >= intDsize*len) // e>=16n -> x als Ergebnis { return x; } // 0 < e < 16n diff --git a/src/float/lfloat/elem/cl_LF_futrunc.cc b/src/float/lfloat/elem/cl_LF_futrunc.cc index fe31e90..0afe85c 100644 --- a/src/float/lfloat/elem/cl_LF_futrunc.cc +++ b/src/float/lfloat/elem/cl_LF_futrunc.cc @@ -29,18 +29,18 @@ const cl_LF futruncate (const cl_LF& x) // e>=16n -> Ergebnis x. #if 0 var cl_signean sign; - var sintL exp; + var sintE exp; var const uintD* mantMSDptr; var uintC mantlen; LF_decode(x, { return x; }, sign=,exp=,mantMSDptr=,mantlen=,); if (exp<=0) { return encode_LF1s(sign,mantlen); } // e<=0 -> Ergebnis +-1.0 - if ((uintL)exp >= intDsize*mantlen) // e>=16n -> x als Ergebnis + if ((uintE)exp >= intDsize*mantlen) // e>=16n -> x als Ergebnis { return x; } else // 0 < e < 16n { // Testen, ob alle hinteren 16n-e Bits =0 sind: - var uintC count = floor((uintL)exp,intDsize); // zu kopierende Digits, < mantlen - var uintC bitcount = ((uintL)exp) % intDsize; // zu kopierende Bits danach, >=0, =0, len; - var uintL uexp = TheLfloat(x)->expo; + var uintE uexp = TheLfloat(x)->expo; if (uexp <= LF_exp_mid) { if (uexp == 0) { return x; } // x=0.0 -> Ergebnis 0.0 return encode_LF1s(TheLfloat(x)->sign,len); // e<=0 -> Ergebnis +-1.0 } - var uintL exp = uexp - LF_exp_mid; + var uintE exp = uexp - LF_exp_mid; if (exp >= intDsize*len) // e>=16n -> x als Ergebnis { return x; } // 0 < e < 16n diff --git a/src/float/lfloat/elem/cl_LF_mul.cc b/src/float/lfloat/elem/cl_LF_mul.cc index 8790695..5429a6c 100644 --- a/src/float/lfloat/elem/cl_LF_mul.cc +++ b/src/float/lfloat/elem/cl_LF_mul.cc @@ -30,10 +30,10 @@ const cl_LF operator* (const cl_LF& x1, const cl_LF& x2) var uintC len1 = TheLfloat(x1)->len; var uintC len2 = TheLfloat(x2)->len; var uintC len = (len1 < len2 ? len1 : len2); // min. L�ge n von x1 und x2 - var uintL uexp1 = TheLfloat(x1)->expo; + var uintE uexp1 = TheLfloat(x1)->expo; if (uexp1==0) // x1=0.0 -> Ergebnis 0.0 { if (len < len1) return shorten(x1,len); else return x1; } - var uintL uexp2 = TheLfloat(x2)->expo; + var uintE uexp2 = TheLfloat(x2)->expo; if (uexp2==0) // x2=0.0 -> Ergebnis 0.0 { if (len < len2) return shorten(x2,len); else return x2; } // Exponenten addieren: @@ -49,7 +49,7 @@ const cl_LF operator* (const cl_LF& x1, const cl_LF& x2) } } else // Carry - { if (uexp1 > (uintL)(LF_exp_mid+LF_exp_high+1)) { cl_error_floating_point_overflow(); } } + { if (uexp1 > (uintE)(LF_exp_mid+LF_exp_high+1)) { cl_error_floating_point_overflow(); } } uexp1 = uexp1 - LF_exp_mid; // Nun ist LF_exp_low <= uexp1 <= LF_exp_high+1. // neues Long-Float allozieren: diff --git a/src/float/lfloat/elem/cl_LF_scale.cc b/src/float/lfloat/elem/cl_LF_scale.cc index 3ec4909..7a2738f 100644 --- a/src/float/lfloat/elem/cl_LF_scale.cc +++ b/src/float/lfloat/elem/cl_LF_scale.cc @@ -23,9 +23,9 @@ const cl_LF scale_float (const cl_LF& x, sintC delta) // delta muß ein Integer betragsmäßig <= LF_exp_high-LF_exp_low sein. // Neues LF mit um delta vergrößertem Exponenten bilden. if (delta == 0) { return x; } // delta=0 -> x als Ergebnis - var uintL uexp = TheLfloat(x)->expo; + var uintE uexp = TheLfloat(x)->expo; if (uexp==0) { return x; } - var uintC udelta = delta; + var uintE udelta = delta; if (delta >= 0) { // udelta = delta >=0 if ( ((uexp = uexp+udelta) < udelta) // Exponent-Überlauf? @@ -33,8 +33,8 @@ const cl_LF scale_float (const cl_LF& x, sintC delta) ) { cl_error_floating_point_overflow(); } } else { - // delta <0, udelta = 2^intCsize+delta - if ( ((uintL)(-(uexp = uexp+udelta)) <= (uintC)(-udelta)) // oder Exponent-Unterlauf? + // delta <0, udelta = 2^intEsize+delta + if ( ((uintE)(-(uexp = uexp+udelta)) <= (uintE)(-udelta)) // oder Exponent-Unterlauf? || (uexp < LF_exp_low) // oder Exponent zu klein? ) { cl_error_floating_point_underflow(); } diff --git a/src/float/lfloat/elem/cl_LF_scale_I.cc b/src/float/lfloat/elem/cl_LF_scale_I.cc index daf92fd..6691fd7 100644 --- a/src/float/lfloat/elem/cl_LF_scale_I.cc +++ b/src/float/lfloat/elem/cl_LF_scale_I.cc @@ -24,9 +24,9 @@ const cl_LF scale_float (const cl_LF& x, const cl_I& delta) // delta muß ein Integer betragsmäßig <= LF_exp_high-LF_exp_low sein. // Neues LF mit um delta vergrößertem Exponenten bilden. if (eq(delta,0)) { return x; } // delta=0 -> x als Ergebnis - var uintL uexp = TheLfloat(x)->expo; + var uintE uexp = TheLfloat(x)->expo; if (uexp==0) { return x; } - var uintV udelta; + var uintE udelta; // |delta| muß <= LF_exp_high-LF_exp_low < 2^32 sein. Wie bei I_to_UL: if (fixnump(delta)) { // Fixnum diff --git a/src/float/lfloat/elem/cl_LF_square.cc b/src/float/lfloat/elem/cl_LF_square.cc index 53fb6f8..72d5827 100644 --- a/src/float/lfloat/elem/cl_LF_square.cc +++ b/src/float/lfloat/elem/cl_LF_square.cc @@ -20,12 +20,12 @@ const cl_LF square (const cl_LF& x) { // Methode: wie operator*(x,x). var uintC len = TheLfloat(x)->len; - var uintL uexp = TheLfloat(x)->expo; + var uintE uexp = TheLfloat(x)->expo; if (uexp==0) // x=0.0 -> Ergebnis 0.0 { return x; } // Exponenten addieren: // (uexp-LF_exp_mid) + (uexp-LF_exp_mid) = (2*uexp-LF_exp_mid)-LF_exp_mid - if ((sintL)uexp >= 0) + if ((sintE)uexp >= 0) // kein Carry { uexp = 2*uexp; if (uexp < LF_exp_mid+LF_exp_low) @@ -37,7 +37,7 @@ const cl_LF square (const cl_LF& x) else // Carry { uexp = 2*uexp; - if (uexp > (uintL)(LF_exp_mid+LF_exp_high+1)) { cl_error_floating_point_overflow(); } + if (uexp > (uintE)(LF_exp_mid+LF_exp_high+1)) { cl_error_floating_point_overflow(); } } uexp = uexp - LF_exp_mid; // Nun ist LF_exp_low <= uexp <= LF_exp_high+1. diff --git a/src/float/lfloat/elem/cl_LF_to_I.cc b/src/float/lfloat/elem/cl_LF_to_I.cc index d91421e..fb5655c 100644 --- a/src/float/lfloat/elem/cl_LF_to_I.cc +++ b/src/float/lfloat/elem/cl_LF_to_I.cc @@ -24,7 +24,7 @@ const cl_I cl_LF_to_I (const cl_LF& x) // Methode: // Falls x=0.0, Ergebnis 0. // Sonst (ASH Vorzeichen*Mantisse (e-16n)). - var uintL uexp = TheLfloat(x)->expo; + var uintE uexp = TheLfloat(x)->expo; if (uexp==0) { return 0; } // x=0.0 -> Ergebnis 0 // Mantisse zu einem Integer machen: CL_ALLOCA_STACK; diff --git a/src/float/lfloat/misc/cl_LF_decode.cc b/src/float/lfloat/misc/cl_LF_decode.cc index f35f058..a65e142 100644 --- a/src/float/lfloat/misc/cl_LF_decode.cc +++ b/src/float/lfloat/misc/cl_LF_decode.cc @@ -19,14 +19,14 @@ const decoded_lfloat decode_float (const cl_LF& x) { // x entpacken: var cl_signean sign; - var sintL exp; + var sintE exp; var uintC mantlen; var const uintD* mantMSDptr; LF_decode(x, { return decoded_lfloat(x, 0, encode_LF1(mantlen)); }, sign=,exp=,mantMSDptr=,mantlen=,); return decoded_lfloat( encode_LFu(0,0+LF_exp_mid,mantMSDptr,mantlen), // (-1)^0 * 2^0 * m erzeugen - L_to_I(exp), // e als Fixnum + E_to_I(exp), // e als Fixnum encode_LF1s(sign,mantlen) // (-1)^s erzeugen ); } diff --git a/src/float/lfloat/misc/cl_LF_exponent.cc b/src/float/lfloat/misc/cl_LF_exponent.cc index 7fd12b1..d6f7346 100644 --- a/src/float/lfloat/misc/cl_LF_exponent.cc +++ b/src/float/lfloat/misc/cl_LF_exponent.cc @@ -14,11 +14,11 @@ namespace cln { MAYBE_INLINE -sintL float_exponent (const cl_LF& x) +sintE float_exponent (const cl_LF& x) { - var uintL uexp = TheLfloat(x)->expo; + var uintE uexp = TheLfloat(x)->expo; if (uexp==0) { return 0; } - return (sintL)(uexp - LF_exp_mid); + return (sintE)(uexp - LF_exp_mid); } } // namespace cln diff --git a/src/float/lfloat/misc/cl_LF_idecode.cc b/src/float/lfloat/misc/cl_LF_idecode.cc index b21490c..033e4bd 100644 --- a/src/float/lfloat/misc/cl_LF_idecode.cc +++ b/src/float/lfloat/misc/cl_LF_idecode.cc @@ -19,7 +19,7 @@ MAYBE_INLINE const cl_idecoded_float integer_decode_float (const cl_LF& x) { // x entpacken: - var uintL uexp = TheLfloat(x)->expo; + var uintE uexp = TheLfloat(x)->expo; if (uexp == 0) { return cl_idecoded_float(0, 0, 1); } var cl_signean sign = TheLfloat(x)->sign; diff --git a/src/float/lfloat/misc/cl_LF_shortenrel.cc b/src/float/lfloat/misc/cl_LF_shortenrel.cc index 3207bda..77d5f75 100644 --- a/src/float/lfloat/misc/cl_LF_shortenrel.cc +++ b/src/float/lfloat/misc/cl_LF_shortenrel.cc @@ -28,15 +28,15 @@ const cl_LF cl_LF_shortenrelative (const cl_LF& x, const cl_LF& y) // dx := float_digits(x), dy := float_digits(y). // 1 ulp(x) = 2^(ex-dx), 1 ulp(y) = 2^(ey-dy). // Falls ex-dx < ey-dy, x von Precision dx auf dy-ey+ex verkürzen. - var sintL ey = float_exponent(y); + var sintE ey = float_exponent(y); var sintC dy = float_precision(y); if (dy==0) // zerop(y) ? cl_abort(); - var sintL ex = float_exponent(x); + var sintE ex = float_exponent(x); var sintC dx = float_precision(x); if (dx==0) // zerop(x) ? return x; - var sintL d = ex - ey; + var sintE d = ex - ey; if (ex>=0 && ey<0 && d<0) // d overflow? return x; if (ex<0 && ey>=0 && d>=0) // d underflow? diff --git a/src/float/lfloat/misc/cl_LF_shortenwith.cc b/src/float/lfloat/misc/cl_LF_shortenwith.cc index f0ed744..b12156a 100644 --- a/src/float/lfloat/misc/cl_LF_shortenwith.cc +++ b/src/float/lfloat/misc/cl_LF_shortenwith.cc @@ -27,12 +27,12 @@ const cl_LF cl_LF_shortenwith (const cl_LF& x, const cl_LF& y) // ex := float_exponent(x), dx := float_digits(x), 1 ulp(x) = 2^(ex-dx). // ey := float_exponent(y). // Falls ex-dx < ey, x von Precision dx auf ex-ey verkürzen. - var sintL ey = float_exponent(y); - var sintL ex = float_exponent(x); + var sintE ey = float_exponent(y); + var sintE ex = float_exponent(x); var uintC dx = float_precision(x); if (dx==0) // zerop(x) ? return x; - var sintL ulpx = ex - dx; + var sintE ulpx = ex - dx; if ((ex<0 && ulpx>=0) // underflow? || (ulpx < ey) ) { // Now ex-dx < ey, hence ex-ey < dx. diff --git a/src/float/misc/cl_F_decode.cc b/src/float/misc/cl_F_decode.cc index 0c16e7c..42f9771 100644 --- a/src/float/misc/cl_F_decode.cc +++ b/src/float/misc/cl_F_decode.cc @@ -84,14 +84,14 @@ inline const decoded_float decode_float (const cl_LF& x) { // x entpacken: var cl_signean sign; - var sintL exp; + var sintE exp; var uintC mantlen; var const uintD* mantMSDptr; LF_decode(x, { return decoded_float(x, 0, encode_LF1(mantlen)); }, sign=,exp=,mantMSDptr=,mantlen=,); return decoded_float( encode_LFu(0,0+LF_exp_mid,mantMSDptr,mantlen), // (-1)^0 * 2^0 * m erzeugen - L_to_I(exp), // e als Fixnum + E_to_I(exp), // e als Fixnum encode_LF1s(sign,mantlen) // (-1)^s erzeugen ); } diff --git a/src/float/misc/cl_F_exponent.cc b/src/float/misc/cl_F_exponent.cc index 10c89ba..8f94561 100644 --- a/src/float/misc/cl_F_exponent.cc +++ b/src/float/misc/cl_F_exponent.cc @@ -20,7 +20,7 @@ namespace cln { -sintL float_exponent (const cl_F& x) +sintE float_exponent (const cl_F& x) { floatcase(x , return float_exponent(x); diff --git a/src/float/misc/cl_F_shortenrel.cc b/src/float/misc/cl_F_shortenrel.cc index d4bea6d..c0b292f 100644 --- a/src/float/misc/cl_F_shortenrel.cc +++ b/src/float/misc/cl_F_shortenrel.cc @@ -22,15 +22,15 @@ const cl_F cl_F_shortenrelative (const cl_F& x, const cl_F& y) // dx := float_digits(x), dy := float_digits(y). // 1 ulp(x) = 2^(ex-dx), 1 ulp(y) = 2^(ey-dy). // Falls ex-dx < ey-dy, x von Precision dx auf dy-ey+ex verkürzen. - var sintL ey = float_exponent(y); + var sintE ey = float_exponent(y); var sintC dy = float_precision(y); if (dy==0) // zerop(y) ? cl_abort(); - var sintL ex = float_exponent(x); + var sintE ex = float_exponent(x); var sintC dx = float_precision(x); if (dx==0) // zerop(x) ? return x; - var sintL d = ex - ey; + var sintE d = ex - ey; if (ex>=0 && ey<0 && d<0) // d overflow? return x; if (ex<0 && ey>=0 && d>=0) // d underflow? diff --git a/src/float/misc/cl_float_format.cc b/src/float/misc/cl_float_format.cc index 5bfbf60..8c6e4df 100644 --- a/src/float/misc/cl_float_format.cc +++ b/src/float/misc/cl_float_format.cc @@ -11,21 +11,26 @@ namespace cln { -float_format_t float_format (uintL n) +float_format_t float_format (uintE n) { // Methode: // Mindestens 1+n Dezimalstellen (inklusive Vorkommastelle) // bedeutet mindestens ceiling((1+n)*ln(10)/ln(2)) Binärstellen. -// ln(10)/ln(2) = 3.321928095 = (binär) 11.01010010011010011110000100101111... -// = (binär) 100 - 0.10101101100101100001111011010001 +// ln(10)/ln(2) = 3.321928095 = (binär) 11.0101001001101001111000010010111100110100... +// = (binär) 100 - 0.1010110110010110000111101101000111001011... // Durch diese Berechnungsmethode wird das Ergebnis sicher >= (1+n)*ln(10)/ln(2) -// sein, evtl. um ein paar Bit zu groß, aber nicht zu klein. +// sein, evtl. um ein paar Bit zu groß aber nicht zu klein. n = 1+n; return (float_format_t) ((n << 2) - - (n >> 1) - (n >> 3) - (n >> 5) - (n >> 6) - (n >> 8) - - (n >> 9) - (n >> 12) - (n >> 14) - (n >> 15) - ); + - (n >> 1) - (n >> 3) - (n >> 5) - (n >> 6) + - (n >> 8) - (n >> 9) - (n >> 12) - (n >> 14) + - (n >> 15) - (n >> 20) - (n >> 21) - (n >> 22) + - (n >> 23) - (n >> 25) - (n >> 26) - (n >> 28) +#if (intEsize>32) + - (n >> 32) - (n >> 33) - (n >> 34) - (n >> 35) +#endif + ); } } // namespace cln diff --git a/src/float/output/cl_F_dprint.cc b/src/float/output/cl_F_dprint.cc index 79d182c..9f871a0 100644 --- a/src/float/output/cl_F_dprint.cc +++ b/src/float/output/cl_F_dprint.cc @@ -61,13 +61,14 @@ namespace cln { // typedef struct cl_decimal_decoded_float { char * a; - uintL k; + uintC k; cl_I e; cl_I s; // Constructor. - cl_decimal_decoded_float (char * ap, uintL kp, const cl_I& ep, const cl_I& sp) : a(ap), k(kp), e(ep), s(sp) {} + cl_decimal_decoded_float (char * ap, uintC kp, const cl_I& ep, const cl_I& sp) : a(ap), k(kp), e(ep), s(sp) {} }; + static const cl_decimal_decoded_float decode_float_decimal (const cl_F& x) { var cl_idecoded_float x_idecoded = integer_decode_float(x); @@ -120,21 +121,41 @@ static const cl_decimal_decoded_float decode_float_decimal (const cl_F& x) // 1838395/6107016 1936274/6432163 13456039/44699994 // 15392313/51132157 44240665/146964308 59632978/198096465 // 103873643/345060773 475127550/1578339557 579001193/1923400330 + // 24793177656/82361153417 149338067129/496090320832 + // 174131244785/578451474249 845863046269/2809896217828 + // 1865857337323/6198243909905 6443435058238/21404627947543 // ) // e>=0 : wähle lg(2) < a/b < lg(2) + 1/e, // dann ist d <= floor(e*a/b) <= d+1 . // e<0 : wähle lg(2) - 1/abs(e) < a/b < lg(2), // dann ist d <= floor(e*a/b) <= d+1 . - // Es ist bekannt, daß abs(e) <= 2^31 + 2^20 . + // Es ist bekannt, dass abs(e) <= 2^31 + 2^32*64, falls intEsize == 32, + // bzw. dass abs(e) <= 2^63 + 2^64*64, falls intEsize == 64. + // (Hierbei steht 64 für die maximale intDsize und es wurde benutzt, + // dass intEsize >= intCsize.) // Unser d sei := floor(e*a/b)-1. (d /= 0, da abs(e) >= 7.) d = minus1(minusp(e) ? (e >= -970 ? floor1(e*3,10) // Näherungsbruch 3/10 - : floor1(e*21306,70777) // Näherungsbruch 21306/70777 +#if (intEsize==32) + : floor1(e*97879,325147) // Näherungsbruch 97879/325147 +#else + : (e >= -1800000000LL + ? floor1(e*8651,28738) // Näherungsbruch 8651/28738 + : floor1(e*24793177656LL,82361153417LL) // Näherungsbruch 24793177656/82361153417 + ) +#endif ) : (e <= 22000 ? floor1(e*28,93) // Näherungsbruch 28/93 - : floor1(e*12655,42039) // Näherungsbruch 12655/42039 +#if (intEsize==32) + : floor1(e*1838395,6107016) // Näherungsbruch 1838395/6107016 +#else + : (e <= 3300000000LL + ? floor1(e*12655,42039) // Näherungsbruch 12655/42039 + : floor1(e*149338067129LL,496090320832LL) // Näherungsbruch 149338067129/496090320832 + ) +#endif ) ); // Das wahre d wird durch diese Schätzung entweder getroffen @@ -248,10 +269,18 @@ static const cl_decimal_decoded_float decode_float_decimal (const cl_F& x) // |e| ist recht klein -> man kann 2^e und 10^d exakt ausrechnen if (!minusp(e)) { // e >= 0. Schätze d = floor(e*lg(2)) wie oben. - // Es ist e<=2*l<2^21. + // Es ist e<=2*l<2^39, falls intCsize == 32, + // bzw. e<=2*l<2^71, falls intCsize == 64. d = (e <= 22000 ? floor1(e*28,93) // Näherungsbruch 28/93 - : floor1(e*4004,13301) // Näherungsbruch 4004/13301 +#if (intCsize==32) + : floor1(e*1838395,6107016) // Näherungsbruch 1838395/6107016 +#else + : (e <= 3300000000LL + ? floor1(e*12655,42039) // Näherungsbruch 12655/42039 + : floor1(e*149338067129LL,496090320832LL) // Näherungsbruch 149338067129/496090320832 + ) +#endif ); // Das wahre d wird durch diese Schätzung entweder getroffen // oder um 1 überschätzt, aber das können wir leicht feststellen. @@ -267,10 +296,18 @@ static const cl_decimal_decoded_float decode_float_decimal (const cl_F& x) a2 = floor1(minus1(ash(oben,e)),zehn_d); } else { // e < 0. Schätze d = floor(e*lg(2)) wie oben. - // Es ist |e|<=2*l<2^21. + // Es ist |e|<=2*l<2^39, falls intCsize == 32, + // bzw. |e|<=2*l<2^71, falls intCsize == 64. d = (e >= -970 ? floor1(e*3,10) // Näherungsbruch 3/10 - : floor1(e*643,2136) // Näherungsbruch 643/2136 +#if (intCsize==32) + : floor1(e*97879,325147) // Näherungsbruch 97879/325147 +#else + : (e >= -1800000000LL + ? floor1(e*8651,28738) // Näherungsbruch 8651/28738 + : floor1(e*24793177656LL,82361153417LL) // Näherungsbruch 24793177656/82361153417 + ) +#endif ); // Das wahre d wird durch diese Schätzung entweder getroffen // oder um 1 überschätzt, aber das können wir leicht feststellen. @@ -317,8 +354,8 @@ static const cl_decimal_decoded_float decode_float_decimal (const cl_F& x) // Nun a in einen Dezimalstring umwandeln // und dann Nullen am Schluß streichen: var char* as = cl_decimal_string(a); // Ziffernfolge zu a>0 - var uintL las = ::strlen(as); // Länge der Ziffernfolge - var uintL k = las; // Länge ohne die gestrichenen Nullen am Schluß + var uintC las = ::strlen(as); // Länge der Ziffernfolge + var uintC k = las; // Länge ohne die gestrichenen Nullen am Schluß var cl_I ee = k+d; // a * 10^d = a * 10^(-k+ee) while (as[k-1] == '0') // eine 0 am Schluß? { // ja -> a := a / 10 (wird aber nicht mehr gebraucht), @@ -355,7 +392,7 @@ static const cl_decimal_decoded_float decode_float_decimal (const cl_F& x) } } var char* as = cl_decimal_string(a); // Ziffernfolge zu a>0 - var uintL k = ::strlen(as); + var uintC k = ::strlen(as); ASSERT(as[k-1] != '0'); return cl_decimal_decoded_float(as,k,k+d,sign); } @@ -365,7 +402,7 @@ void print_float (std::ostream& stream, const cl_print_float_flags& flags, const { var cl_decimal_decoded_float z_decoded = decode_float_decimal(z); var char * & mantstring = z_decoded.a; - var uintL& mantlen = z_decoded.k; + var uintC& mantlen = z_decoded.k; var cl_I& expo = z_decoded.e; var cl_I& sign = z_decoded.s; // arg in Dezimaldarstellung: +/- 0.mant * 10^expo, wobei @@ -405,7 +442,7 @@ void print_float (std::ostream& stream, const cl_print_float_flags& flags, const fprintchar(stream,mantstring[i]); } fprintchar(stream,'.'); - { for (uintL i = scale; i < mantlen; i++) + { for (uintC i = scale; i < mantlen; i++) fprintchar(stream,mantstring[i]); } } else { diff --git a/src/float/sfloat/misc/cl_SF_exponent.cc b/src/float/sfloat/misc/cl_SF_exponent.cc index d9c8f3b..54adfe5 100644 --- a/src/float/sfloat/misc/cl_SF_exponent.cc +++ b/src/float/sfloat/misc/cl_SF_exponent.cc @@ -14,7 +14,7 @@ namespace cln { MAYBE_INLINE -sintL float_exponent (const cl_SF& x) +sintE float_exponent (const cl_SF& x) { var uintL uexp = SF_uexp(x); if (uexp==0) { return 0; } diff --git a/src/float/transcendental/cl_F_atanhx.cc b/src/float/transcendental/cl_F_atanhx.cc index 2d15ede..3faf4b9 100644 --- a/src/float/transcendental/cl_F_atanhx.cc +++ b/src/float/transcendental/cl_F_atanhx.cc @@ -53,12 +53,12 @@ const cl_LF atanhx (const cl_LF& x) return x; var uintC actuallen = TheLfloat(x)->len; var uintC d = float_digits(x); - var sintL e = float_exponent(x); + var sintE e = float_exponent(x); if (e <= (sintC)(-d)>>1) // e <= -d/2 <==> e <= -ceiling(d/2) return x; // ja -> x als Ergebnis if (actuallen >= 34) { DeclareType(cl_LF,x); - 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)); return cl_float(scale_float(ln((1+xx)/(1-xx)),-1),x); } var uintL k = 0; // Rekursionszähler k:=0 @@ -67,7 +67,7 @@ const cl_LF atanhx (const cl_LF& x) // schlecht). Ein guter Wert ist: // für naive1: limit_scope = 0.625 = 5/8, // für naive2: limit_scope = 0.4 = 13/32. - var uintL sqrt_d = floor(isqrt(d)*13,32); // limit_slope*floor(sqrt(d)) + var uintL sqrt_d = floor(isqrtC(d)*13,32); // limit_slope*floor(sqrt(d)) var cl_LF xx = x; if (e >= (sintL)(-sqrt_d)) { // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern. @@ -130,14 +130,14 @@ const cl_F atanhx (const cl_F& x) if (zerop(x)) return x; var uintC d = float_digits(x); - var sintL e = float_exponent(x); + var sintE e = float_exponent(x); if (e <= (sintC)(-d)>>1) // e <= -d/2 <==> e <= -ceiling(d/2) return x; // ja -> x als Ergebnis var uintL k = 0; // Rekursionszähler k:=0 // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe // angewandt werden. limit_slope = 1.0 ist schlecht (ca. 15% zu // schlecht). Ein guter Wert ist limit_scope = 0.625 = 5/8. - var uintL sqrt_d = floor(isqrt(d)*5,8); // limit_slope*floor(sqrt(d)) + var uintL sqrt_d = floor(isqrtC(d)*5,8); // limit_slope*floor(sqrt(d)) var cl_F xx = x; if (e >= (sintL)(-sqrt_d)) { // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern. diff --git a/src/float/transcendental/cl_F_atanx.cc b/src/float/transcendental/cl_F_atanx.cc index ad852a1..b670853 100644 --- a/src/float/transcendental/cl_F_atanx.cc +++ b/src/float/transcendental/cl_F_atanx.cc @@ -51,7 +51,7 @@ static const cl_LF atanx_naive (const cl_LF& x) return x; var uintC actuallen = TheLfloat(x)->len; var uintC d = float_digits(x); - var sintL e = float_exponent(x); + var sintE e = float_exponent(x); if (e <= (sintC)(-d)>>1) // e <= -d/2 <==> e <= -ceiling(d/2) return x; // ja -> x als Ergebnis var uintL k = 0; // Rekursionszähler k:=0 @@ -61,7 +61,7 @@ static const cl_LF atanx_naive (const cl_LF& x) // Für naive1: limit_scope = 0.5. // Für naive2: limit_scope = 0.375 (ca. 0.5 für kleine len, 0.35 für // große len). - var uintL sqrt_d = floor(isqrt(d)*3,8); // limit_slope*floor(sqrt(d)) + var uintL sqrt_d = floor(isqrtC(d)*3,8); // limit_slope*floor(sqrt(d)) var cl_LF xx = x; if (e >= (sintL)(-sqrt_d)) { // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern. @@ -120,11 +120,11 @@ static const cl_F atanx_naive (const cl_F& x) if (zerop(x)) return x; var uintC d = float_digits(x); - var sintL e = float_exponent(x); + var sintE e = float_exponent(x); if (e <= (sintC)(-d)>>1) // e <= -d/2 <==> e <= -ceiling(d/2) return x; // ja -> x als Ergebnis var uintL k = 0; // Rekursionszähler k:=0 - var uintL sqrt_d = floor(isqrt(d),2); // limit_slope*floor(sqrt(d)) + var uintL sqrt_d = floor(isqrtC(d),2); // limit_slope*floor(sqrt(d)) // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe // angewandt werden. limit_slope = 1.0 ist schlecht (ca. 20% zu // schlecht). Ein guter Wert ist limit_scope = 0.5. diff --git a/src/float/transcendental/cl_F_cosh.cc b/src/float/transcendental/cl_F_cosh.cc index da29b15..a64b733 100644 --- a/src/float/transcendental/cl_F_cosh.cc +++ b/src/float/transcendental/cl_F_cosh.cc @@ -29,7 +29,7 @@ const cl_F cosh (const cl_F& x) // cosh(x) = 1+x*y*(sinh(y)/y)^2 errechnen. // falls e>=0: y:=exp(x) errechnen, (scale-float (+ y (/ y)) -1) bilden. - var sintL e = float_exponent(x); + var sintE e = float_exponent(x); if (e < 0) { // Exponent e abtesten // e<0 if (zerop(x)) diff --git a/src/float/transcendental/cl_F_expx.cc b/src/float/transcendental/cl_F_expx.cc index 1695723..c7956a3 100644 --- a/src/float/transcendental/cl_F_expx.cc +++ b/src/float/transcendental/cl_F_expx.cc @@ -48,19 +48,19 @@ const cl_LF expx_naive (const cl_LF& x) return cl_float(1,x); var uintL actuallen = TheLfloat(x)->len; var uintC d = float_digits(x); - var sintL e = float_exponent(x); + var sintE e = float_exponent(x); if (e < -(sintC)d) // e < -d ? return cl_float(1,x); // ja -> 1.0 als Ergebnis { Mutable(cl_LF,x); - var uintL k = 0; // Rekursionszähler k:=0 + var uintE k = 0; // Rekursionszähler k:=0 // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe // angewandt werden. limit_slope = 1.0 ist nicht schlecht, // auch im Bereich d = ca. 800. - var sintL e_limit = -1-isqrt(d); // -1-floor(sqrt(d)) + var sintL e_limit = -1-isqrtC(d); // -1-floor(sqrt(d)) if (e > e_limit) { // e > -1-floor(sqrt(d)) -> muß |x| verkleinern. k = e - e_limit; - x = scale_float(x,-(sintL)k); // x := x/2^k + x = scale_float(x,-(sintE)k); // x := x/2^k // Neuer Exponent = e-k = e_limit. } // Potenzreihe anwenden: @@ -94,22 +94,22 @@ const cl_F expx_naive (const cl_F& x) if (zerop(x)) return cl_float(1,x); var uintC d = float_digits(x); - var sintL e = float_exponent(x); + var sintE e = float_exponent(x); if (e < -(sintC)d) // e < -d ? return cl_float(1,x); // ja -> 1.0 als Ergebnis { Mutable(cl_F,x); - var uintL k = 0; // Rekursionszähler k:=0 + var uintE k = 0; // Rekursionszähler k:=0 // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe // angewandt werden. limit_slope = 1.0 ist nicht schlecht. Für // d > 1600 scheint der Bereich 2.0 <= limit_slope <= 2.6 am besten // zu sein (mit bis zu 15% Beschleunigung gegenüber limit_slope = 1.0), // aber in diesem Bereich rechnen wir gar nicht. // Wir wählen limit_slope = 1.5. - var sintL e_limit = -1-floor(isqrt(d)*3,2); // -1-floor(sqrt(d)) + var sintL e_limit = -1-floor(isqrtC(d)*3,2); // -1-floor(sqrt(d)) if (e > e_limit) { // e > -1-floor(sqrt(d)) -> muß |x| verkleinern. k = e - e_limit; - x = scale_float(x,-(sintL)k); // x := x/2^k + x = scale_float(x,-(sintE)k); // x := x/2^k // Neuer Exponent = e-k = e_limit. } // Potenzreihe anwenden: diff --git a/src/float/transcendental/cl_F_lnx.cc b/src/float/transcendental/cl_F_lnx.cc index 73af930..e6b2ea7 100644 --- a/src/float/transcendental/cl_F_lnx.cc +++ b/src/float/transcendental/cl_F_lnx.cc @@ -49,7 +49,7 @@ const cl_LF lnx_naive (const cl_LF& x) return y; var uintL actuallen = TheLfloat(x)->len; var uintC d = float_digits(x); - var sintL e = float_exponent(y); + var sintE e = float_exponent(y); if (e <= -(sintC)d) // e <= -d ? return y; // ja -> y als Ergebnis { Mutable(cl_LF,x); @@ -60,7 +60,7 @@ const cl_LF lnx_naive (const cl_LF& x) // für ln(1+y), naive2: limit_slope = 11/16 = 0.7, // für atanh(z), naive1: limit_slope = 0.6, // für atanh(z), naive1: limit_slope = 0.5. - var sintL e_limit = -1-floor(isqrt(d),2); // -1-floor(sqrt(d)) + var sintL e_limit = -1-floor(isqrtC(d),2); // -1-floor(sqrt(d)) while (e > e_limit) { // e > -1-floor(sqrt(d)) -> muß |y| verkleinern. x = sqrt(x); // x := (sqrt x) @@ -147,13 +147,13 @@ const cl_F lnx_naive (const cl_F& x) if (zerop(y)) // y=0.0 -> y als Ergebnis return y; var uintC d = float_digits(x); - var sintL e = float_exponent(y); + var sintE e = float_exponent(y); if (e <= -(sintC)d) // e <= -d ? return y; // ja -> y als Ergebnis { Mutable(cl_F,x); var uintL k = 0; // Rekursionszähler k:=0 // Bei e <= -1-floor(sqrt(d)) kann die Potenzreihe angewandt werden. - var sintL e_limit = -1-isqrt(d); // -1-floor(sqrt(d)) + var sintL e_limit = -1-isqrtC(d); // -1-floor(sqrt(d)) while (e > e_limit) { // e > -1-floor(sqrt(d)) -> muß |y| verkleinern. x = sqrt(x); // x := (sqrt x) diff --git a/src/float/transcendental/cl_F_sinhx.cc b/src/float/transcendental/cl_F_sinhx.cc index e5fd97c..d2b61e6 100644 --- a/src/float/transcendental/cl_F_sinhx.cc +++ b/src/float/transcendental/cl_F_sinhx.cc @@ -52,13 +52,13 @@ const cl_F sinhxbyx_naive (const cl_F& x) if (zerop(x)) return cl_float(1,x); var uintC d = float_digits(x); - var sintL e = float_exponent(x); + var sintE e = float_exponent(x); if (e <= (1-(sintC)d)>>1) // e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ? return cl_float(1,x); // ja -> 1.0 als Ergebnis { Mutable(cl_F,x); // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe // angewandt werden. Wähle limit_slope = 13/32 = 0.4. - var sintL e_limit = -1-floor(isqrt(d)*13,32); // -1-floor(sqrt(d)) + var sintL e_limit = -1-floor(isqrtC(d)*13,32); // -1-floor(sqrt(d)) if (e > e_limit) { // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern. x = scale_float(x,e_limit-e); @@ -82,7 +82,7 @@ const cl_F sinhxbyx_naive (const cl_F& x) while (e > e_limit) { z = z + x2 * square(z); x2 = scale_float(x2,2); // x^2 := x^2*4 - e_limit++; + e--; } return z; }} @@ -116,15 +116,15 @@ const cl_LF sinhx_naive (const cl_LF& x) return x; var uintL actuallen = TheLfloat(x)->len; var uintC d = float_digits(x); - var sintL e = float_exponent(x); + var sintE e = float_exponent(x); if (e <= (1-(sintC)d)>>1) // e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ? return square(x); // ja -> x^2 als Ergebnis { Mutable(cl_LF,x); - var sintL ee = e; + var sintE ee = e; // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe // angewandt werden. Ein guter Wert für naive1 ist limit_slope = 0.6, // für naive3 aber limit_slope = 0.5. - var sintL e_limit = -1-floor(isqrt(d),2); // -1-floor(sqrt(d)) + var sintL e_limit = -1-floor(isqrtC(d),2); // -1-floor(sqrt(d)) if (e > e_limit) { // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern. x = scale_float(x,e_limit-e); @@ -182,7 +182,7 @@ const cl_LF sinhx_naive (const cl_LF& x) var cl_LF z = square(powser_value); // sinh^2 als Ergebnis while (e > e_limit) { z = square(cl_float(1,x) + scale_float(z,1)) - cl_float(1,x); // z := (1+2*z)^2-1 - e_limit++; + e--; } return z; }} diff --git a/src/float/transcendental/cl_F_sinx.cc b/src/float/transcendental/cl_F_sinx.cc index 5dc12b1..430cfed 100644 --- a/src/float/transcendental/cl_F_sinx.cc +++ b/src/float/transcendental/cl_F_sinx.cc @@ -52,7 +52,7 @@ const cl_F sinxbyx_naive (const cl_F& x) if (zerop(x)) return cl_float(1,x); var uintC d = float_digits(x); - var sintL e = float_exponent(x); + var sintE e = float_exponent(x); if (e <= (-(sintC)d)>>1) // e <= (-d)/2 <==> e <= -ceiling(d/2) ? return cl_float(1,x); // ja -> 1.0 als Ergebnis { Mutable(cl_F,x); @@ -67,7 +67,7 @@ const cl_F sinxbyx_naive (const cl_F& x) // 100 0.40-0.45 // 200 0.40-0.45 // Wähle limit_slope = 13/32 = 0.4. - var sintL e_limit = -1-floor(isqrt(d)*13,32); // -1-floor(sqrt(d)) + var sintL e_limit = -1-floor(isqrtC(d)*13,32); // -1-floor(sqrt(d)) if (e > e_limit) { // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern. x = scale_float(x,e_limit-e); @@ -91,7 +91,7 @@ const cl_F sinxbyx_naive (const cl_F& x) while (e > e_limit) { z = z - x2 * square(z); x2 = scale_float(x2,2); // x^2 := x^2*4 - e_limit++; + e--; } return z; }} @@ -125,16 +125,16 @@ const cl_LF sinx_naive (const cl_LF& x) return x; var uintL actuallen = TheLfloat(x)->len; var uintC d = float_digits(x); - var sintL e = float_exponent(x); + var sintE e = float_exponent(x); if (e <= (-(sintC)d)>>1) // e <= (-d)/2 <==> e <= -ceiling(d/2) ? return square(x); // ja -> x^2 als Ergebnis { Mutable(cl_LF,x); - var sintL ee = e; + var sintE ee = e; // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe // angewandt werden. limit_slope = 1.0 ist schlecht (ca. 10% zu // schlecht). Ein guter Wert für naive1 ist limit_slope = 0.6, // für naive3 aber limit_slope = 0.5. - var sintL e_limit = -1-floor(isqrt(d),2); // -1-floor(sqrt(d)) + var sintL e_limit = -1-floor(isqrtC(d),2); // -1-floor(sqrt(d)) if (e > e_limit) { // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern. x = scale_float(x,e_limit-e); @@ -192,7 +192,7 @@ const cl_LF sinx_naive (const cl_LF& x) var cl_LF z = square(powser_value); // sin^2 als Ergebnis while (e > e_limit) { z = cl_float(1,x) - square(cl_float(1,x) - scale_float(z,1)); // z := 1-(1-2*z)^2 - e_limit++; + e--; } return z; }} diff --git a/src/float/transcendental/cl_LF_pi.cc b/src/float/transcendental/cl_LF_pi.cc index 311ce3b..ba8d2cd 100644 --- a/src/float/transcendental/cl_LF_pi.cc +++ b/src/float/transcendental/cl_LF_pi.cc @@ -62,7 +62,7 @@ const cl_LF compute_pi_brent_salamin (uintC len) // (/ (expt a 2) t) // ) var uintC actuallen = len + 1; // 1 Schutz-Digit - var uintC uexp_limit = LF_exp_mid - intDsize*len; + var uintE uexp_limit = LF_exp_mid - intDsize*len; // Ein Long-Float ist genau dann betragsmäßig <2^-n, wenn // sein Exponent < LF_exp_mid-n = uexp_limit ist. var cl_LF a = cl_I_to_LF(1,actuallen); @@ -112,7 +112,7 @@ const cl_LF compute_pi_brent_salamin_quartic (uintC len) // Hence, // pi = AGM(1,1/sqrt(2))^2 * 1/(1/2 - sum(k even, 2^k*[....])). var uintC actuallen = len + 1; // 1 Schutz-Digit - var uintC uexp_limit = LF_exp_mid - intDsize*len; + var uintE uexp_limit = LF_exp_mid - intDsize*len; var cl_LF one = cl_I_to_LF(1,actuallen); var cl_LF a = one; var cl_LF wa = one; diff --git a/src/integer/cl_I.h b/src/integer/cl_I.h index 52c137d..282aece 100644 --- a/src/integer/cl_I.h +++ b/src/integer/cl_I.h @@ -234,6 +234,11 @@ inline sint64 FN_to_Q (const cl_I& x) return cl_I(cl_I_constructor_from_UQ(wert)); } + extern cl_private_thing cl_I_constructor_from_Q2 (sint64 wert_hi, uint64 wert_lo ); + inline const cl_I Q2_to_I( sint64 wert_hi, uint64 wert_lo) + { + return cl_I(cl_I_constructor_from_Q2(wert_hi, wert_lo)); + } #endif // Wandelt Doppel-Longword in Integer um. @@ -290,6 +295,26 @@ inline sint64 FN_to_Q (const cl_I& x) #define UV_to_I(wert) UQ_to_I(wert) #endif +// Wandelt sintE in Integer um. +// E_to_I(wert) +// > wert: Wert des Integers, ein sintE. +// < ergebnis: Integer mit diesem Wert. +#if (intEsize<=32) + #define E_to_I(wert) L_to_I(wert) +#else + #define E_to_I(wert) Q_to_I(wert) +#endif + +// Wandelt uintE in Integer >=0 um. +// UE_to_I(wert) +// > wert: Wert des Integers, ein uintE. +// < ergebnis: Integer mit diesem Wert. +#if (intEsize<=32) + #define UE_to_I(wert) UL_to_I(wert) +#else + #define UE_to_I(wert) UQ_to_I(wert) +#endif + // Wandelt uintD in Integer >=0 um. // UD_to_I(wert) // > wert: Wert des Integers, ein uintD. @@ -313,6 +338,14 @@ inline const cl_I minus (uintL x, uintL y) #endif } +#ifdef intQsize + +inline const cl_I minus (uintQ x, uintQ y) +{ + return Q2_to_I( (x Longword: diff --git a/src/integer/conv/cl_I_from_Q2.cc b/src/integer/conv/cl_I_from_Q2.cc new file mode 100644 index 0000000..e75830e --- /dev/null +++ b/src/integer/conv/cl_I_from_Q2.cc @@ -0,0 +1,196 @@ +// Q2_to_I() helper. + +// General includes. +#include "cl_sysdep.h" + +// Specification. +#include "cl_I.h" + + +// Implementation. + +#include "cln/number.h" + +#ifdef intQsize + +#include "cl_DS.h" + +namespace cln { + +cl_private_thing cl_I_constructor_from_Q2 (sint64 wert_hi, uint64 wert_lo) +{ + if (wert_hi == 0) { + if ((wert_lo & minus_bit(cl_value_len-1)) == 0) + return (cl_private_thing)(cl_combine(cl_FN_tag,wert_lo)); + } + elif (wert_hi == ~(sint64)0) { + if ((~wert_lo & minus_bit(cl_value_len-1)) == 0) + return (cl_private_thing)(cl_combine(cl_FN_tag,(sint64)wert_lo)); + } + // Create bignum with length n, where: + // bn_minlength <= n <= ceiling(128/intDsize) + #define FILL_1_DIGIT(l,i,from) \ + arrayLSref(ptr->data,l,i) = (uintD)from; + #define FILL_2_DIGIT(l,i,from) \ + arrayLSref(ptr->data,l,i) = (uintD)from; \ + arrayLSref(ptr->data,l,i+1) = (uintD)(from>>intDsize); + #define FILL_3_DIGIT(l,i,from) \ + arrayLSref(ptr->data,l,i) = (uintD)from; from>>=intDsize; \ + arrayLSref(ptr->data,l,i+1) = (uintD)from; \ + arrayLSref(ptr->data,l,i+2) = (uintD)(from>>intDsize); + #define FILL_4_DIGIT(l,i,from) \ + arrayLSref(ptr->data,l,i) = (uintD)from; from>>=intDsize; \ + arrayLSref(ptr->data,l,i+1) = (uintD)from; from>>=intDsize; \ + arrayLSref(ptr->data,l,i+2) = (uintD)from; \ + arrayLSref(ptr->data,l,i+3) = (uintD)(from>>intDsize); + #define FILL_8_DIGIT(l,i,from) \ + arrayLSref(ptr->data,l,i) = (uintD)from; from >>=intDsize; \ + arrayLSref(ptr->data,l,i+1) = (uintD)from; from >>=intDsize; \ + arrayLSref(ptr->data,l,i+2) = (uintD)from; from >>=intDsize; \ + arrayLSref(ptr->data,l,i+3) = (uintD)from; from >>=intDsize; \ + arrayLSref(ptr->data,l,i+4) = (uintD)from; from >>=intDsize; \ + arrayLSref(ptr->data,l,i+5) = (uintD)from; from >>=intDsize; \ + arrayLSref(ptr->data,l,i+6) = (uintD)from; \ + arrayLSref(ptr->data,l,i+7) = (uintD)from>>intDsize; + #if (intDsize==64) + #define FILL_1 FILL_1_DIGIT(1,0,wert_lo); + #define FILL_2 FILL_1_DIGIT(2,1,wert_hi); FILL_1_DIGIT(2,0,wert_lo); + #endif + #if (32/intDsize==1) + #define FILL_1 FILL_1_DIGIT(1,0,wert_lo); + #define FILL_2 FILL_2_DIGIT(2,0,wert_lo); + #define FILL_3 FILL_1_DIGIT(3,2,wert_hi); FILL_2_DIGIT(3,0,wert_lo); + #define FILL_4 FILL_2_DIGIT(4,2,wert_hi); FILL_2_DIGIT(4,0,wert_lo); + #endif + #if (32/intDsize==2) + #define FILL_1 FILL_1_DIGIT(1,0,wert_lo); + #define FILL_2 FILL_2_DIGIT(2,0,wert_lo); + #define FILL_3 FILL_3_DIGIT(3,0,wert_lo); + #define FILL_4 FILL_4_DIGIT(4,0,wert_lo); + #define FILL_5 FILL_1_DIGIT(5,4,wert_hi); FILL_4_DIGIT(5,0,wert_lo); + #define FILL_6 FILL_2_DIGIT(6,4,wert_hi); FILL_4_DIGIT(6,0,wert_lo); + #define FILL_7 FILL_3_DIGIT(7,4,wert_hi); FILL_4_DIGIT(7,0,wert_lo); + #define FILL_8 FILL_4_DIGIT(8,4,wert_hi); FILL_4_DIGIT(8,0,wert_lo); + #endif + #if (32/intDsize==4) + #define FILL_1 FILL_1_DIGIT(1,0,wert_lo); + #define FILL_2 FILL_2_DIGIT(2,0,wert_lo); + #define FILL_3 FILL_3_DIGIT(3,0,wert_lo); + #define FILL_4 FILL_4_DIGIT(4,0,wert_lo); + #define FILL_5 FILL_5_DIGIT(5,0,wert_lo); + #define FILL_6 FILL_6_DIGIT(6,0,wert_lo); + #define FILL_7 FILL_7_DIGIT(7,0,wert_lo); + #define FILL_8 FILL_8_DIGIT(8,0,wert_lo); + #define FILL_9 FILL_1_DIGIT(9,8,wert_hi); FILL_8_DIGIT(9,0,wert_lo); + #define FILL_10 FILL_2_DIGIT(10,8,wert_hi); FILL_8_DIGIT(10,0,wert_lo); + #define FILL_11 FILL_3_DIGIT(11,8,wert_hi); FILL_8_DIGIT(11,0,wert_lo); + #define FILL_12 FILL_4_DIGIT(12,8,wert_hi); FILL_8_DIGIT(12,0,wert_lo); + #define FILL_13 FILL_5_DIGIT(13,8,wert_hi); FILL_8_DIGIT(13,0,wert_lo); + #define FILL_14 FILL_6_DIGIT(14,8,wert_hi); FILL_8_DIGIT(14,0,wert_lo); + #define FILL_15 FILL_7_DIGIT(15,8,wert_hi); FILL_8_DIGIT(15,0,wert_lo); + #define FILL_16 FILL_8_DIGIT(16,8,wert_hi); FILL_8_DIGIT(16,0,wert_lo); + #endif + if (wert_hi >= 0) { + #define IF_LENGTH(i) \ + if ((bn_minlength <= i) && (i*intDsize <= 128)) \ + if (!((i+1)*intDsize <= 128) \ + || (i*intDsize-1 < 64 \ + ? ((wert_hi == 0) && (wert_lo < (uint64)bitc(i*intDsize-1))) \ + : ((uint64)wert_hi < (uint64)bitc(i*intDsize-1-64)) \ + ) ) + #define ALLOC(i) \ + var cl_heap_bignum* ptr = allocate_bignum(i); + #define OK \ + return (cl_private_thing)(ptr); + IF_LENGTH(1) + bignum1: { ALLOC(1); FILL_1; OK; } + IF_LENGTH(2) + bignum2: { ALLOC(2); FILL_2; OK; } + #if (intDsize <= 32) + IF_LENGTH(3) + bignum3: { ALLOC(3); FILL_3; OK; } + IF_LENGTH(4) + bignum4: { ALLOC(4); FILL_4; OK; } + #if (intDsize <= 16) + IF_LENGTH(5) + bignum5: { ALLOC(5); FILL_5; OK; } + IF_LENGTH(6) + bignum6: { ALLOC(6); FILL_6; OK; } + IF_LENGTH(7) + bignum7: { ALLOC(7); FILL_7; OK; } + IF_LENGTH(8) + bignum8: { ALLOC(8); FILL_8; OK; } + #if (intDsize <= 8) + IF_LENGTH(9) + bignum9: { ALLOC(9); FILL_9; OK; } + IF_LENGTH(10) + bignum10: { ALLOC(10); FILL_10; OK; } + IF_LENGTH(11) + bignum11: { ALLOC(11); FILL_11; OK; } + IF_LENGTH(12) + bignum12: { ALLOC(12); FILL_12; OK; } + IF_LENGTH(13) + bignum13: { ALLOC(13); FILL_13; OK; } + IF_LENGTH(14) + bignum14: { ALLOC(14); FILL_14; OK; } + IF_LENGTH(15) + bignum15: { ALLOC(15); FILL_15; OK; } + IF_LENGTH(16) + bignum16: { ALLOC(16); FILL_16; OK; } + #endif + #endif + #endif + #undef IF_LENGTH + } else { + #define IF_LENGTH(i) \ + if ((bn_minlength <= i) && (i*intDsize <= 128)) \ + if (!((i+1)*intDsize <= 128) \ + || (i*intDsize-1 < 64 \ + ? ((wert_hi == ~(sint64)0) && (wert_lo >= (uint64)(-bitc(i*intDsize-1)))) \ + : ((uint64)wert_hi >= (uint64)(-bitc(i*intDsize-1-64))) \ + ) ) + IF_LENGTH(1) + goto bignum1; + IF_LENGTH(2) + goto bignum2; + #if (intDsize <= 32) + IF_LENGTH(3) + goto bignum3; + IF_LENGTH(4) + goto bignum4; + #if (intDsize <= 16) + IF_LENGTH(5) + goto bignum5; + IF_LENGTH(6) + goto bignum6; + IF_LENGTH(7) + goto bignum7; + IF_LENGTH(8) + goto bignum8; + #if (intDsize <= 8) + IF_LENGTH(9) + goto bignum9; + IF_LENGTH(10) + goto bignum10; + IF_LENGTH(11) + goto bignum11; + IF_LENGTH(12) + goto bignum12; + IF_LENGTH(13) + goto bignum13; + IF_LENGTH(14) + goto bignum14; + IF_LENGTH(15) + goto bignum15; + IF_LENGTH(16) + goto bignum16; + #endif + #endif + #endif + #undef IF_LENGTH + } +} + +} // namespace cln + +#endif