Commit Graph

116 Commits

Author SHA1 Message Date
Stefan Krah f21587e3a8 mpd_qpowmod(): calculate result with zero-exponent for compatibility with
decimal.py. The hack to remove the ideal exponent is no longer required.
2012-08-23 15:05:29 +02:00
Stefan Krah 17183a2972 Use zero bottom margin. 2012-08-22 19:28:12 +02:00
Stefan Krah ad5b43995e In the 32-bit build, dec_hash() raised InvalidOperation if the operand
had a coefficient with MAX_PREC=425000000 digits and a negative exponent.
Increasing the context limits above the official values fixes the issue
and is safe (in this case!).
2012-08-22 19:11:50 +02:00
Stefan Krah 2fd502f6a1 1) Use _mpd_basedivmod() regardless of the length of the dividend. This is
required for a corner case in dec_hash() in the following commit and also
   usually faster. dec_hash() needs some extra precision above MPD_MAX_PREC,
   and _mpd_base_ndivmod() is not audited for that.

2) Use _mpd_basemul() if the length of the smaller operand is less than
   or equal to 256. While this is technically an optimization, it is
   required for *testing* corner cases in dec_hash() in reasonable time.
2012-08-22 18:54:37 +02:00
Stefan Krah b9e36b97f2 Add unused parameter to a couple of METH_NOARGS functions. The previous
form is used in many places in the source tree, but was found to be incorrect
in a recent tracker discussion.
2012-07-28 13:53:47 +02:00
Stefan Krah 26a1c7a905 Fix Visual Studio warning. 2012-07-20 12:34:18 +02:00
Stefan Krah e574402bd6 Issue #7652: Clean up _mpd_qinvroot() and mark it LIBMPDEC_ONLY. Use the
algorithm from decimal.py for mpd_qsqrt().
2012-07-12 21:17:59 +02:00
Stefan Krah d57caf36bd Remove ISSUES.txt. 2012-07-01 12:24:20 +02:00
Stefan Krah 5431e30853 After 79d2eb29c755 it is no longer necessary to zero the output array:
None of the _mpd_shortadd() or _mpd_shortmul() functions read uninitialized
values. Previously zeroing was required since _mpd_real_size() was called
on the output array.
2012-06-30 21:57:49 +02:00
Stefan Krah c35a8e5c98 Proactive reliability fix for broken FPUs: The base conversion functions
use log10() to calculate the size of the output array. The current code
has been tested on x86/amd64 (and to a lesser extent on qemu-mips qemu-sparc)
and produces sufficiently large values for all inputs tested so far (coefficient
sizes of 10**18 - 1 are hard to test exhaustively).

The new code does not rely on the correctness of log10() and resizes
the output arrays if the allocated space is insufficient.
2012-06-30 18:05:33 +02:00
Stefan Krah 1edab78859 Update test script to Visual Studio 2010. 2012-06-25 14:41:37 +02:00
Stefan Krah 39e810eb6c Make the benchmark more fair for _decimal/decimal.py by setting context.prec
only once (float obviously doesn't set any context at all).
2012-06-24 14:10:49 +02:00
Stefan Krah 78f075636c Speed up _decimal by another 10-15% by caching the thread local context
that was last accessed. In the pi benchmark (64-bit platform, prec=9),
_decimal is now only 1.5x slower than float.
2012-06-24 12:20:03 +02:00
Stefan Krah 3077ab8237 Whitespace. 2012-06-23 00:31:04 +02:00
Stefan Krah 50b0a365ba Fix comment. 2012-06-20 23:38:51 +02:00
Stefan Krah 22385011ed Many cleanups of redundant code in mpd_qrem_near():
1) _mpd_qdivmod() uses the context precision only in two places, and
     the new code should be exactly equivalent to the previous code.

  2) Remove misleading comment.

  3) The quotient *is* an integer with exponent 0, so calling mpd_qtrunc()
     is pointless.

  4) Replace two instances of identical code by a single one.

  5) Use _mpd_cmp_abs() instead of mpd_cmp_total_mag(): the operands
     are not special.

  6) Don't clear MPD_Rounded in the status (with the current code it should
     not be set within the function).
2012-06-20 23:34:58 +02:00
Stefan Krah 9c1feb88f3 Add comments to the power functions, in particular to _mpd_qpow_real(). 2012-06-18 19:57:23 +02:00
Stefan Krah c62bd13cb2 1) State the relative errors of the power functions for integer exponents.
2) _mpd_qpow_mpd(): Abort the loop for all specials, not only infinity.

3) _mpd_qpow_mpd(): Make the function more general and distinguish between
   zero clamping and folding down the exponent. The latter case is currently
   handled by setting context->clamp to 0 before calling the function.

4) _mpd_qpow_int(): Add one to the work precision in case of a negative
   exponent. This is to get the same relative error (0.1 * 10**-prec)
   for both positive and negative exponents. The previous relative
   error for negative exponents was (0.2 * 10**-prec).

   Both errors are _before_ the final rounding to the context precision.
2012-06-16 19:45:35 +02:00
Stefan Krah b7832939c7 1) Fix signature of _mpd_qpow_uint(): contrary to the comment base is constant.
2) Abort the loop for all specials, not only infinity.

3) Make the function more general and distinguish between zero clamping
   and folding down the exponent. The latter case is currently handled
   by setting context->clamp to 0 before calling the function.
2012-06-12 21:06:06 +02:00
Stefan Krah 88e19779ad 1) Replace long-winded abort() construct by assert().
2) Remove micro optimization (inline checking for NaN before calling
   mpd_qcheck_nans()) that probably has no benefit in this case.
2012-06-11 08:57:17 +02:00
Stefan Krah 9253862f45 1) State restrictions for the transform length.
2) Switch argument order to match the function signature of mpd_calloc()
   (cosmetic change, since the order is irrelevant).
2012-06-10 16:50:55 +02:00
Stefan Krah afc0c77b42 Add one extra comparison to the _mpd_shortmul() case to avoid repetitive code. 2012-06-09 15:28:36 +02:00
Stefan Krah 5248a2d3c1 Enumerate all cases in the overflow detection strategy in mpd_qlog10(). 2012-06-09 00:01:28 +02:00
Stefan Krah 1cf6dfc8b2 1) List relative error for _mpd_qln10().
2) Add rigorous error analysis to _mpd_qlog10 (ACL2 proofs exist).

3) Use the relative error as a basis for the interval generation in the
   correction loop (same as in _mpd_qln()).
2012-06-08 18:41:33 +02:00
Stefan Krah 7bda265662 1) The overflow detection in mpd_qln() has a surprising number of case splits.
List all of them in the comment.

2) Use the recently stated relative error of _mpd_qln() to generate the
   interval for the exact value of ln(x). See also the comment in mpd_qexp().
2012-06-07 17:48:47 +02:00
Stefan Krah a3394bce33 1) Add error analysis comments to mpd_qln10() and _mpd_qln().
2) Simplify the precision adjustment code for values in [0.900, 1.15].
2012-06-06 15:57:18 +02:00
Stefan Krah 67ee1d05dd word.digits are always initialized before use in the Taylor series loop,
but this is more readable.
2012-06-01 10:58:16 +02:00
Stefan Krah 0271766c88 Use workctx instead of ctx for cosmetic reasons. Also zero-pad the result
in the simple path (not correctly rounded but faster).
2012-05-31 20:49:24 +02:00
Stefan Krah 4d3e0a695a Improve Underflow handling in the correct-rounding loop. The case for
Underflow to zero hasn't changed: _mpd_qexp() internally uses MIN_EMIN,
so the result would also underflow to zero for all emin > MIN_EMIN.

In case digits are left, the informal argument is as follows: Underflow can
occur only once in the last multiplication of the power stage (in the Horner
stage Underflow provably cannot occur, and if Underflow occurred twice in
the power stage, the result would underflow to zero on the second occasion).

Since there is no double rounding during Underflow, the effective work
precision is now 1 <= result->digits < prec. It can be shown by a somewhat
tedious argument that abs(result - e**x) < ulp(result, result->digits).

Therefore the correct rounding loop now uses ulp(result, result->digits)
to generate the bounds for e**x in case of Underflow.
2012-05-31 20:01:05 +02:00
Stefan Krah 9a5beece1b Improve comments. 2012-05-31 16:21:34 +02:00
Stefan Krah 5ddbcfc53e Pad the result with zeros just before the final rounding. 2012-05-31 16:00:21 +02:00
Stefan Krah 30c35e8154 Do not clobber existing flags. 2012-05-31 15:09:27 +02:00
Stefan Krah e34a209584 Fix Visual Studio warning. 2012-05-16 20:20:03 +02:00
Stefan Krah 696d10f1bb Changes in _mpd_qexp():
-----------------------

  1) Reduce the number of iterations in the Horner scheme for operands with
     a negative adjusted exponent. Previously the number was overestimated
     quite generously.

  2) The function _mpd_get_exp_iterations() now has an ACL2 proof and
     is rewritten accordingly.

  3) The proof relies on abs(op) > 9 * 10**(-prec-1), so operands without
     that property are now handled by the new function _mpd_qexp_check_one().

  4) The error analysis for the evaluation of the truncated Taylor series
     in Hull&Abrham's paper relies on the fact that the reduced operand
     'r' has fewer than context.prec digits.

     Since the operands may have more than context.prec digits, a new ACL2
     proof covers the case that r.digits > context.prec. To facilitate the
     proof, the Horner step now uses fma instead of rounding twice in
     multiply/add.


Changes in mpd_qexp():
----------------------

  1) Fix a bound in the correct rounding loop that was too optimistic. In
     practice results were always correctly rounded, because it is unlikely
     that the error in _mpd_qexp() ever reaches the theoretical maximum.
2012-05-16 20:10:21 +02:00
Stefan Krah 9d3a5aeabe Defensive programming: mpd_isspecial(r) already implies mpd_isspecial(q), but
this is more readable.
2012-04-20 21:00:31 +02:00
Stefan Krah 3c23a87e58 The divmod function for large numbers now has an ACL2 proof. Related changes:
1) Rename _mpd_qbarrett_divmod into _mpd_base_ndivmod: The function is
     only marginally related to either Barrett's algorithm or to the version
     in Hasselstrom's paper.

  2) In places where the proof assumes exact operations, use new versions of
     add/sub/multiply that set NaN/Invalid_operation if this condition is
     not met. According to the proof this cannot happen, so this should be
     regarded as an extra safety net.

  3) Raise Division_impossible for operands with a number of digits greater
     than MPD_MAX_PREC. This facilitates the audit of the function and can
     practically only occur in the 32-bit version under conditions where
     a MemoryError is already imminent.

  4) Use _mpd_qmul() in places where the result can exceed MPD_MAX_PREC in
     a well defined manner.

  5) Test for mpd_isspecial(qq) in a place where the addition of one
     can theoretically trigger a Malloc_error.

  6) Remove redundant code in _mpd_qdivmod().

  7) Add many comments.
2012-04-20 19:59:20 +02:00
Stefan Krah c51b7fd65b 1) Simplify comment -- one has to read the complete proof (available in ACL2)
in order to understand the algorithm anyway.

2) v->exp == -v->digits may be assumed.

3) Fix comment (v always shares data with a).
2012-04-18 19:27:32 +02:00
Stefan Krah 5d0d2e2b04 Explain the strategy to avoid huge alignment shifts in _mpd_qadd() in detail. 2012-04-18 18:59:56 +02:00
Stefan Krah ed4b21ff4f Cosmetic change: initialize digits to 1 (redundant). 2012-04-18 18:45:22 +02:00
Stefan Krah bc771e9b19 Remove redundant finalization of the result. 2012-04-18 18:25:37 +02:00
Stefan Krah aecaf0b663 Fix comments and whitespace. 2012-04-18 18:08:20 +02:00
Stefan Krah 6369f77d20 Support mythical ones' complement machines. 2012-04-18 17:57:56 +02:00
Stefan Krah 140893cbaa The previous code is correct, but hard to verify: The libmpdec documentation
rightfully states that an mpd_t with a coefficient flagged as MPD_CONST_DATA
must not be in the position of the result operand. In this particular case
several assumptions guarantee that a resize will never occur in all possible
code paths, which was the reason for using MPD_CONST_DATA and saving an
instruction by omitting the initialization of tmp.alloc.

For readability, tmp is now flagged as MPD_STATIC_DATA and tmp.alloc
is initialized.
2012-04-18 17:48:34 +02:00
Stefan Krah ec766a6179 1) Remove claim of an input invariant that is only true for static mpd_t.
Resizing is used _inside_ libmpdec functions, and it is permitted to
   change x->alloc several times while setting x->len at the end of the
   function. Therefore, for dynamic mpd_t x->alloc can _temporarily_ drop
   below x->len. Of course the final result always has x->len <= x->alloc.

   For static mpd_t this cannot happen, since resizing to a smaller
   coefficient is a no-op.

2) Remove micro optimization in mpd_switch_to_dyn(): Previously only the
   valid initialized part of the existing coefficient up to x->len was
   copied to the new dynamic memory area. Now copying does the same as
   realloc() and the entire old memory area is copied.

   The rationale for this change is that it is no longer needed to memorize
   the explanation given in 1).
2012-04-10 23:11:54 +02:00
Stefan Krah 7b544ca08d Fix stale comment. 2012-04-10 23:08:29 +02:00
Stefan Krah cc74b6ab14 Issue #14478: Cache the hash of a Decimal in the C version. 2012-04-10 16:27:58 +02:00
Stefan Krah e37f8b29fc Issue #14520: Add __sizeof__() method to the Decimal object. 2012-04-09 21:27:20 +02:00
Stefan Krah f69aef747a Resize the coefficient to MPD_MINALLOC also if the requested size is below
MPD_MINALLOC. Previously the resize was skipped as a micro optimization.
2012-04-09 20:47:57 +02:00
Stefan Krah dd159ce606 Speed up _decimal by 30-40% for numerical workloads by improving the cache
locality for regularly sized coefficients.
2012-04-09 20:24:57 +02:00
Stefan Krah 94ef3e4c1d Use the MPD() accessor macro. 2012-04-09 19:20:46 +02:00
Stefan Krah dc36efa368 1) Fix comment.
2) Assert that the source operand is not special. Prevent resulting assert
   failure (harmless) by initializing flags before calling mpd_qshiftr_inplace.

3) Save a couple of instructions (mpd_zerocoeff already sets digits and len).
   Reorder initialization to match the order in the mpd_t struct.
2012-04-07 15:57:59 +02:00
Stefan Krah 4771cca817 Whitespace. 2012-04-05 16:15:25 +02:00
Stefan Krah 871b96bd5a Reduce array size. 2012-04-05 16:07:22 +02:00
Stefan Krah a6169484c2 Formatting. 2012-04-05 15:48:59 +02:00
Stefan Krah ff3eca0cc3 Allow printing a leading '-' and the maximum number of exponent digits
rather than raising RuntimeError (allocated space is sufficient for the
additional character).
2012-04-05 15:46:19 +02:00
Stefan Krah 0774e9b9f5 Raise InvalidOperation if exponents of zeros are clamped during exact
conversion in the Decimal constructor. Exact here refers to the
representation and not to the value (clamping does not change the value).
2012-04-05 15:21:58 +02:00
Stefan Krah 91c0274bc4 Improve comments. 2012-04-02 20:51:08 +02:00
Stefan Krah 5100171d81 Clear the context flags if a context is initialized from the DefaultContext. 2012-04-02 15:02:21 +02:00
Stefan Krah 41e031004b Fix Overflow exception in the bignum factorial benchmark that is due to
the recent change of the default value for context.Emax.
2012-04-01 23:25:34 +02:00
Stefan Krah 0e41981cd5 Use abort() rather than exit() to appease tools like rpmlint. abort() is used
in libmpdec to prevent undefined behavior if an invalid context is used. This
cannot occur for the _decimal module since user input for the context is
validated.
2012-03-30 14:12:20 +02:00
Stefan Krah fe17b2bc77 Raise MemoryError instead of InvalidOperation/MallocError for compatibility
with decimal.py. The standard specifies InsufficientStorage (MallocError) as
a sub-condition of InvalidOperation. This allows a calculation to continue
with NaN results when allocation fails.
2012-03-25 18:59:21 +02:00
Stefan Krah c64150bcac Fix formatting after removing tabs. 2012-03-23 16:34:41 +01:00
Stefan Krah cd9e1d0205 Whitespace. 2012-03-23 16:22:05 +01:00
Stefan Krah b6405efd1b Use the same exception hierarchy as decimal.py. FloatOperation now also
inherits from TypeError. Cleanup in module initialization to make repeated
import failures robust.
2012-03-23 14:46:48 +01:00
Stefan Krah 7cc5521d40 Whitespace. 2012-03-21 20:21:20 +01:00
Stefan Krah 1919b7e72b Issue #7652: Integrate the decimal floating point libmpdec library to speed
up the decimal module. Performance gains of the new C implementation are
between 12x and 80x, depending on the application.
2012-03-21 18:25:23 +01:00