+2006-12-11 Richard B. Kreckel <kreckel@ginac.de>
+
+ 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 <kreckel@ginac.de>
* src/base/digitseq/cl_DS.h: #undef DS, needed for i386-Solaris.
+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
==========================
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
@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).
int
main (int argc, char * argv[])
{
- int digits = 100;
+ long digits = 100;
if (argc > 1) {
if (argc == 2 && !strcmp(argv[1],"--help")) {
usage(cout);
return 0;
}
if (argc == 2 && isdigit(argv[1][0])) {
- digits = atoi(argv[1]);
+ digits = atol(argv[1]);
} else {
usage(cerr);
return 1;
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
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.
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)
// 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)
// 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)
// 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.
{ 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)
{ 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)
{ 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)
{ 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)
{ 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)
{ 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)
{ 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)
{ 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)
// 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)
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
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:
{ 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)
{ 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);
{ 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)
{ 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);
{ 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)
{ 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);
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.
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
{ 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)
{ 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);
{ 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)
{ 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);
// 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)
}
#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 {
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
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.
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)
{ 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)
{ 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);
{ 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)
{ 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);
{ 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)
{ 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);
{ 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)
{ 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.
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; }
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; }
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.
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
{ 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)
{ 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)
{ 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)
{ 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)
{ 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)
{ 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)
{ 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)
{ 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)
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<cl_F&>(x = x + y); }
inline cl_F& operator+= (cl_R& x, const double y) { return static_cast<cl_F&>(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; }
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<cl_F&>(x = x - y); }
inline cl_F& operator-= (cl_R& x, const double y) { return static_cast<cl_F&>(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; }
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; }
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
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
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.
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)
// 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)
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
return (uint32)high << 16;
}
-#ifdef HAVE_FAST_LONGLONG
+#ifdef HAVE_LONGLONG
// High-Word einer 64-Bit-Zahl bestimmen
// high32(wert)
return (uint64)high << 32;
}
-#endif /* HAVE_FAST_LONGLONG */
+#endif /* HAVE_LONGLONG */
// Multipliziert zwei 16-Bit-Zahlen miteinander und liefert eine 32-Bit-Zahl:
// < 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)
// 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))
}}
}
+#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
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'
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'
// 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:
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; }
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; }
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)
);
(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;
}
}
// 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:
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:
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
};
// 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; }
// 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;
// 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))); \
// 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 */
// 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
// 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:
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 ?
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
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
}
// 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(); }
}
// 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();
{ 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
if (x_uexp > y_uexp) return signean_plus; // x>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| -> x<y
}
var uintC len1 = TheLfloat(x1)->len;
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:
}
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
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);
{ 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
// 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, <intDsize
+ var uintC count = floor((uintE)exp,intDsize); // zu kopierende Digits, < mantlen
+ var uintC bitcount = ((uintE)exp) % intDsize; // zu kopierende Bits danach, >=0, <intDsize
var uintD mask = minus_bit(intDsize-bitcount-1); // Maske mit bitcount+1 Bits
var const uintD* mantptr = mantMSDptr mspop count;
if ((mspref(mantptr,0) & -mask) ==0) goto ab; // Bit 16n-e-1 =0 -> abrunden
}
#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
// 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
{ 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, <intDsize
+ { var uintC count = floor((uintE)exp,intDsize); // zu kopierende Digits, < mantlen
+ var uintC bitcount = ((uintE)exp) % intDsize; // zu kopierende Bits danach, >=0, <intDsize
var uintD* ptr =
copy_loop_msp(mantMSDptr,MSDptr,count); // count ganze Digits kopieren
msprefnext(ptr) = mspref(mantMSDptr,count) & minus_bitm(intDsize-bitcount); // dann bitcount Bits kopieren
}
#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
// 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, <intDsize
+ var uintC count = floor((uintE)exp,intDsize); // zu kopierende Digits, < mantlen
+ var uintC bitcount = ((uintE)exp) % intDsize; // zu kopierende Bits danach, >=0, <intDsize
var uintD mask = minus_bitm(intDsize-bitcount); // Maske mit bitcount Bits
var uintD* mantptr = mantMSDptr mspop count;
if ( ((mspref(mantptr,0) & ~mask) ==0)
}
#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_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
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:
} }
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:
// 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?
)
{ 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(); }
// 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
{
// 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)
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.
// 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;
{
// 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
);
}
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
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;
// 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?
// 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.
{
// 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
);
}
namespace cln {
-sintL float_exponent (const cl_F& x)
+sintE float_exponent (const cl_F& x)
{
floatcase(x
, return float_exponent(x);
// 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?
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
// 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);
// 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
// |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.
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.
// 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),
}
}
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);
}
{
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
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 {
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; }
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
// 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.
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.
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
// 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.
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.
// 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))
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:
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:
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);
// 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)
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)
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);
while (e > e_limit) {
z = z + x2 * square(z);
x2 = scale_float(x2,2); // x^2 := x^2*4
- e_limit++;
+ e--;
}
return z;
}}
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);
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;
}}
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);
// 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);
while (e > e_limit) {
z = z - x2 * square(z);
x2 = scale_float(x2,2); // x^2 := x^2*4
- e_limit++;
+ e--;
}
return z;
}}
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);
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;
}}
// (/ (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);
// 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;
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.
#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.
#endif
}
+#ifdef intQsize
+
+inline const cl_I minus (uintQ x, uintQ y)
+{
+ return Q2_to_I( (x<y ? -1 : 0), x-y );
+}
+
+#endif
// Umwandlungsroutinen Digit sequence <--> Longword:
--- /dev/null
+// 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
git://www.ginac.de / cln.git / commitdiff