Browse Source

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.
master
Richard Kreckel 19 years ago
parent
commit
e528307482
  1. 68
      ChangeLog
  2. 29
      NEWS
  3. 4
      doc/cln.tex
  4. 4
      examples/pi.cc
  5. 13
      include/cln/complex_class.h
  6. 2
      include/cln/dfloat.h
  7. 2
      include/cln/ffloat.h
  8. 52
      include/cln/float.h
  9. 3
      include/cln/floatformat.h
  10. 44
      include/cln/integer.h
  11. 12
      include/cln/integer_class.h
  12. 26
      include/cln/lfloat.h
  13. 43
      include/cln/number.h
  14. 56
      include/cln/rational.h
  15. 12
      include/cln/rational_class.h
  16. 76
      include/cln/real.h
  17. 12
      include/cln/real_class.h
  18. 2
      include/cln/sfloat.h
  19. 12
      include/cln/types.h
  20. 25
      src/base/cl_low.h
  21. 2
      src/base/cl_macros.h
  22. 35
      src/base/low/cl_low_isqrt.cc
  23. 18
      src/complex/algebraic/cl_LF_hypot.cc
  24. 18
      src/complex/elem/division/cl_C_LF_recip.cc
  25. 2
      src/float/cl_F.h
  26. 2
      src/float/dfloat/misc/cl_DF_exponent.cc
  27. 2
      src/float/ffloat/misc/cl_FF_exponent.cc
  28. 4
      src/float/input/cl_F_read.cc
  29. 4
      src/float/lfloat/algebraic/cl_LF_sqrt.cc
  30. 10
      src/float/lfloat/cl_LF.h
  31. 22
      src/float/lfloat/cl_LF_impl.h
  32. 14
      src/float/lfloat/elem/cl_LF_1plus.cc
  33. 4
      src/float/lfloat/elem/cl_LF_I_div.cc
  34. 4
      src/float/lfloat/elem/cl_LF_I_mul.cc
  35. 8
      src/float/lfloat/elem/cl_LF_compare.cc
  36. 6
      src/float/lfloat/elem/cl_LF_div.cc
  37. 6
      src/float/lfloat/elem/cl_LF_from_I.cc
  38. 12
      src/float/lfloat/elem/cl_LF_fround.cc
  39. 12
      src/float/lfloat/elem/cl_LF_ftrunc.cc
  40. 12
      src/float/lfloat/elem/cl_LF_futrunc.cc
  41. 6
      src/float/lfloat/elem/cl_LF_mul.cc
  42. 8
      src/float/lfloat/elem/cl_LF_scale.cc
  43. 4
      src/float/lfloat/elem/cl_LF_scale_I.cc
  44. 6
      src/float/lfloat/elem/cl_LF_square.cc
  45. 2
      src/float/lfloat/elem/cl_LF_to_I.cc
  46. 4
      src/float/lfloat/misc/cl_LF_decode.cc
  47. 6
      src/float/lfloat/misc/cl_LF_exponent.cc
  48. 2
      src/float/lfloat/misc/cl_LF_idecode.cc
  49. 6
      src/float/lfloat/misc/cl_LF_shortenrel.cc
  50. 6
      src/float/lfloat/misc/cl_LF_shortenwith.cc
  51. 4
      src/float/misc/cl_F_decode.cc
  52. 2
      src/float/misc/cl_F_exponent.cc
  53. 6
      src/float/misc/cl_F_shortenrel.cc
  54. 19
      src/float/misc/cl_float_format.cc
  55. 65
      src/float/output/cl_F_dprint.cc
  56. 2
      src/float/sfloat/misc/cl_SF_exponent.cc
  57. 10
      src/float/transcendental/cl_F_atanhx.cc
  58. 8
      src/float/transcendental/cl_F_atanx.cc
  59. 2
      src/float/transcendental/cl_F_cosh.cc
  60. 16
      src/float/transcendental/cl_F_expx.cc
  61. 8
      src/float/transcendental/cl_F_lnx.cc
  62. 14
      src/float/transcendental/cl_F_sinhx.cc
  63. 14
      src/float/transcendental/cl_F_sinx.cc
  64. 4
      src/float/transcendental/cl_LF_pi.cc
  65. 33
      src/integer/cl_I.h
  66. 196
      src/integer/conv/cl_I_from_Q2.cc

68
ChangeLog

@ -1,3 +1,71 @@
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.

29
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
==========================

4
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).

4
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;

13
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)

2
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)

2
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)

52
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)

3
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

44
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);

12
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

26
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)

43
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)

56
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; }

12
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

76
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<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; }
@ -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<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; }
@ -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

12
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)

2
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)

12
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

25
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)

2
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))

35
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

18
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'

18
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'

2
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:

2
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; }

2
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; }

4
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;
}
}

4
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:

10
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; }

22
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

14
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

4
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(); }

4
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();

8
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
if (x_uexp > 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| -> x<y
}

6
src/float/lfloat/elem/cl_LF_div.cc

@ -41,9 +41,9 @@ 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ä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

6
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

12
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, <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
@ -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

12
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, <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
@ -48,12 +48,12 @@ const cl_LF ftruncate (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

12
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, <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)
@ -64,12 +64,12 @@ const cl_LF futruncate (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_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

6
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:

8
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(); }

4
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

6
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.

2
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;

4
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
);
}

6
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

2
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;

6
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?

6
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.

4
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
);
}

2
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);

6
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?

19
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

65
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 {

2
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; }

10
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.

8
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.

2
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))

16
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:

8
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)

14
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;
}}

14
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;
}}

4
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;

33
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<y ? -1 : 0), x-y );
}
#endif
// Umwandlungsroutinen Digit sequence <--> Longword:

196
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
|||||||
100:0
Loading…
Cancel
Save