[PATCH 2/2] chinrem_gcd: return correct integer content when GCD is a number.
This patch fixes a silly typo and makes chinrem_gcd work correctly with
polynomials which are "almost" (up to an integer coefficient) relatively
prime.
Thanks to Ernst Moritz Hahn for a bug report and a test case.
[PATCH 1/2] extract_integer_content: check for rational numbers properly.
Commit 3d09388a (titled as `[bugfix] chinrem_gcd: handle polynomials over
rationals properly.') broke extract_integer_content: now it always returns 1.
The check for rational `integer_contnent' introduced by that commit is wrong
(since integers is a subset of rationals). Rewrite the check proprerly.
[bugfix] chinrem_gcd: handle polynomials over rationals properly.
* extract_integer_content():
integer_content() can also return a rational number, e.g if the expression
is a polynomial over rationals (in this case expr/content is polynomial
over integers with integer content being 1). Therefore check if
integer_content() is really an integer (and if it's not just return 1,
GCD for polynomials over a field is defined up to arbitrary element of
the field).
This fixes possible segfault when computing GCD of polynomials over
rationals (this is not theoretical, see the added test case).
Thanks to Ernst Moritz Hahn for reporting this bug.
mul: algebraic_subs_mul(), has(): don't write beyond the end of array
algebraic_match_mul_with_mul() iterates over operands of mul, that is
for (size_t i=0; i<e.nops(); ++i)
However, the size of arrays (`vectors' in STL speak) passed to this
function is seq.size(), which is nops() - 1 for any mul object. Thus
algebraic_match_mul_with_mul() accesses beyond the arrays limit. Usually
it's not a problem, since any reasonable implementation of std::vector<bool>
packs booleans into ints (or longs). However, some STL implementations
(in particular, the one shipped with msvc) are more picky, and access
beyond the vector<bool> limits results in a segfault. Therefore let's
play safe and allocate proper number of elements (that is, nops()) for
those arrays (subsed and currsubsed).
Avoid infinite loop when unarchiving realsymbol and possymbol.
symbol::read_archive(): explicitly set status_flags::evaluated (and
status_flags::expanded) on object being unarchived. These flags get
reset by basic::operator=(const basic&) for realsymbol and possymbol,
and nothing sets (except symbol ctor), so automatic evaluation never
terminates (or rather, terminates due to a stack overflow). Therefore
it's necessary need to set status_flags::evaluated explicitly.
Thanks to Markus Fröb for a bugreport and a test case.
Richard Kreckel [Wed, 22 Sep 2010 22:40:38 +0000 (00:40 +0200)]
Make sure add::eval() collects all numeric terms.
Apparently, add::eval() assumed that none of the elements of its epvector
has a numeric rest. However, nothing guarantees that -- in particular
evalchildren() doesn't (and actually cannot) do so. Since there are many
places where a new add is constructed directly from an epvector, enforcing
this doesn't make sense either. One example where it did fail was found by
Burgin Erocal: real_part(1+2*(sqrt(2)+1)*(sqrt(2)-1)) returned 1+2*1, not 3.
Richard Kreckel [Wed, 15 Sep 2010 07:11:57 +0000 (09:11 +0200)]
Be more careful about final top-level substitution.
Substituting x==log(x) in exp(x) erroneously returned log(x) because of a
final subst(x==log(x)) after having eval'ed exp(log(x)) -> x. This final
substitution is wrong in the general case. On the other hand, the intent
is to syntactically substitute functions of a given kind etc. This patch
suppresses the final top-level substitution unless the intermediate result
is a container.
Thanks to Burcin Erocal for reporting this bug (originally described by
Kees van Schaik on sage-support@googlegroops.com).
Richard Kreckel [Mon, 13 Sep 2010 20:39:05 +0000 (22:39 +0200)]
Fix autoconf warning about AC_INIT use.
Apparently, autoconf doesn't like the angle brackets in the bug-report
argument any more. It comlains "not a literal: <ginac-list@ginac.de>".
So let's remove it. And while at it, also provide tarname and url
arguments.
While working on fsolve bug I've noticed the following in valgrind log:
==17455== 136 (56 direct, 80 indirect) bytes in 1 blocks are definitely lost in loss record 16 of 19
==17455== at 0x4C249C7: operator new(unsigned long) (vg_replace_malloc.c:220)
==17455== by 0x516CA70: GiNaC::power::eval(int) const (power.cpp:540)
==17455== by 0x4FC1E39: GiNaC::ex::construct_from_basic(GiNaC::basic const&) (ex.cpp:310)
==17455== by 0x406FBF: main (ex.h:255)
Heap allocated objects definitely need the status_flags::dyncallocated flag.
fsolve: avoid useless numerical evaluation of the function
Don't compute f(x) if new x is outside of the interval. We don't need that
value anyway, and the function might be difficult to compute numerically or
even ill defined outside the interval.
As a result fsolve is able to find root(s) of some weird functions.
For example
In general, ex_to is unsafe, and should be used only after proper checks.
evalf() may return non-numeric expression for various reasons (bad
convergence, floating point over- or underflow, out of memory, etc).
So let's add missing checks.
Richard Kreckel [Sat, 22 May 2010 20:34:42 +0000 (22:34 +0200)]
Be more careful with conjugate(f(x)) -> f(conjugate(x)).
That identity is correct for holomorphic functions, but we have to be
careful with some functions not being holomorphic everywhere. On branch
cuts, it is wrong. For a discussion, see:
<http://www.ginac.de/pipermail/ginac-list/2010-April/001601.html>.
This patch changes the default behavior of user-defined functions. From
now on, a user-defined conjugate_func must be registered, in order to
evaluate conjugate(f(x)) -> f(conjugate(x)). This patch also adds such
functions for most initially known functions.
pgcd(), chinrem_gcd(): use appropriate definition of the degree.
Effect: fixes (rare) GCD miscalculation.
pgcd() treats polynomials Z_p[x_0, ..., x_n] as R[x_0, ..., x_{n - 1}], where
R = Z_p[x_n]. Therefore one should use correct definition of the degree
(i.e. compute it w.r.t. x_0, ..., x_{n-1}, but NOT w.r.t. x_n!).
One should use appropriate definition of degree (that is, w.r.t. x_0, ..., x_n)
in chinrem_gcd() too.
pgcd(): avoid infinite loop if the first GCD candidate coincides with GCD
136 if (H_lcoeff.is_equal(lc_gcd)) {
137 if ((Hprev-H).expand().smod(pn).is_zero()) // (*)
138 continue;
The check (*) can result in infinite loop if the (very) first GCD candidate
gives the correct GCD: any subsequent iterations give the same result, so
Hprev and H will be equal (hence the infinite loop). That check is not very
useful (AFAIR I was trying to avoid extra division checks), so let's remove it.
Richard Kreckel [Thu, 13 May 2010 16:30:40 +0000 (18:30 +0200)]
Explicit function disambiguation.
This is necessary for compilers where the hyperbolic and the tgamma function
templates are pulled in from <cmath> and compete with GiNaC's function
templates. GCC 4.4 and 4.5 need this, when compiling with -std=cxx0x.
Richard Kreckel [Thu, 25 Mar 2010 22:08:54 +0000 (23:08 +0100)]
Really fixed the atan2(-Pi, 0) bug.
The problem was that atan2_eval assumed that if y is real and not
positive, it must be negative. But this neglects the fact that in
symbolic compution it may just not be possible to determine the
sign. Now, the expression is returned as it is. Ugly, but correct.
Jens Vollinga [Thu, 25 Mar 2010 09:36:41 +0000 (10:36 +0100)]
Fixed a bug in atan2. It gave a division by zero error for calls like
atan2(-Pi,0), because arguments like -Pi were not recognized (via
info_flags) as negative but as real nevertheless.
pgcd(): detect relatively prime polynomials properly...
... so pgcd() doesn't loop forever any more. Division test does not handle
relatively prime polynomials (because C = 1 divides any polynomial). So we
should stop interpolation (and return gcd of contents) if we got relatively
prime images (we should do that for efficiency reasons too).
Richard Kreckel [Tue, 29 Sep 2009 20:25:16 +0000 (22:25 +0200)]
Remove debian/ directory.
The lifecycle of the Debian packaging files is better maintained in
the Debian pool's diff file. Removing these files makes the life of
the GiNaC releaser easier.
Richard Kreckel [Tue, 29 Sep 2009 20:08:11 +0000 (22:08 +0200)]
Properly document how to deal with libtool versioning.
The old description was referring to ginac_lt_age which should always
stay 0, because we promised to mark incompatible binary interfaces by
bumping ginac_minor_version (or even ginac_major_version). Let's get
rid of ginac_lt_age.
modular_matrix: don't use STL iterators in {mul, sub}_col and friends.
Not using STL iterators makes the code simpler and cleaner. Also we don't
have to argue whether the standard ([lib.random.access.iterators]) imposes
any pre/post-condictions on operator+=.
shaker_sort, permutation_sign: fix invalid use of STL iterators.
According to the standard incrementing end(), and decrementing begin()
are ill-defined operations (see [lib.forward.iterators],
[lib.bidirectional.iterators]).
According to [lib.bidirectional.iterators] it's not OK to decrement
an iterator pointing to the beginning of the sequence. Fortunately random
access iterators provided by (current versions of) gcc/libstdc++ don't have
this silly limitation, so the code which works with pointers works with
iterators without any changes at all. However,
- there's no guarantee that the current behavior won't change in the future
- some non-GNU compilers are not that smart (i.e. the program segfaults
upon when decrementing beginning-of-the-sequence iterator).
Changed the parser such that it understands all defined functions
including the user defined ones.
To this end a method has been added to class function to allow the
modified get_default_reader() function to build up a complete prototype
table. The autogen tool is no longer required.
Added include for cstdio header.
gcc 4.0 does not include cstdio as part of other headers like string
anymore. As a result, gcc 4.0 complained about EOF being undefined.
Due to the strange (although permitted by the standard) behaviour of C++
RTTI on woe32 calchash() returns different hash values for equal objects.
As a result automatic evaluation gets spectacularly broken:
examining clifford objects.....({1+t,2+x,3+y,4+z}) - ({1+t,2+x,3+y,4+z}) erroneously returned -{1+t,2+x,3+y,4+z}+{1+t,2+x,3+y,4+z} instead of 0
({1+t,2+x,3+y,4+z}) - ({1+t,2+x,3+y,4+z}) erroneously returned -{1+t,2+x,3+y,4+z}+{1+t,2+x,3+y,4+z} instead of 0
.({1+t,2+x,3+y,4+z}) - ({1+t,2+x,3+y,4+z}) erroneously returned -{1+t,2+x,3+y,4+z}+{1+t,2+x,3+y,4+z} instead of 0
({1+t,2+x,3+y,4+z}) - ({1+t,2+x,3+y,4+z}) erroneously returned {1+t,2+x,3+y,4+z}-{1+t,2+x,3+y,4+z} instead of 0
[skipped]
.......FAIL: exam_clifford.exe
This patch works around `features' of woe32 RTTI, so calchash() works
properly.
Univariate GCD timing: use sr_gcd when appropriate.
The benchmark consists of two parts:
1) timing of different GCD calculation methods (i.e. subresultant PRS,
heuristic, chinese remaindering).
2) timing of different implementations of the same method. The purpose
is to find out how (in)efficient GiNaC::ex is as a representation
of (univariate) polynomials (as a side note, the result is a bit
depressing -- using coefficient vector instead of GiNaC:ex makes
GCD calculation 50x -- 1000x faster).
Now GiNaC uses modular (chinese remaindering) GCD by default, so part 2)
got broken, i.e. instead of (intended) timings
a) (heuristic, GiNaC::ex) versus (heuristic, coefficient vector)
b) (PRS, GiNaC::ex) versus (PRS, coefficient vector)
one gets
a') (heuristic, GiNaC::ex) versus (heuristic, coefficient vector)
b') (chinese remaindering, GiNaC::ex) versus (PRS, coefficient vector)
Set the gcd_options::use_sr_gcd to restore the meaning of the benchmark.
Note: this patch does not affect the library proper.
Jens Vollinga [Fri, 6 Feb 2009 15:07:10 +0000 (16:07 +0100)]
Prettified source code.
- Added copyright and GPL licencing to new files.
- Increased year to 2009.
- Changed guarding macros in header to uniform pattern without leading or
trailing __ (double underscores).
- Put includes of system wide header consistently below own includes (help
a tiny bit to clarify dependencies).
Implement modular multivariate GCD (based on chinese remaindering algorithm).
Both algorithm and implementation are a bit inefficient:
* algorithm does not make use of sparseness of inputs (and most of
multivariate polynomials are quite sparse)
* GiNaC's expressions are essentially immutable. Thus some simple
operations (i.e. multiplying a polynomial by an element of the base ring)
are prohibitively expensive.
* All numbers (i.e. GiNaC::numeric objects) are heap allocated.
Still it's much faster (~5x on bivariate polynomials, ~100x on 3-variate
ones) than (subresultant) PRS algorithm, so gcd() uses modular algorithm
by default.
Jens Vollinga [Thu, 18 Dec 2008 10:10:50 +0000 (11:10 +0100)]
Fixed bug in multivariate factorization. Fractional numerical content was not
extracted and led to a not so finite running behavior. Added new checks for that
bug.
Added static keywords to hide debugging output operators.
The libtool naming scheme cannot guarantee that on all systems, the numbering
is consecutive. It only guarantees that it is increasing. This doesn't matter,
though: there is not incurred cost for numbers that are omitted, except for
shrinking the available space of leftover numbers. Not something we need to
worry about yet. ;-) On the other hand, tricks which we were using to make
version numbers consecutive on GNU/Linux break versioning on other OSes.
Jens Vollinga [Thu, 13 Nov 2008 10:13:05 +0000 (11:13 +0100)]
Fixed lots of bugs in factor_multivariate().
factor_univariate() takes now the content of the polynomial into accout when
choosing a prime. It also returns the chosen prime. This supports the
factor_multivariate() code.
Polynomials are expanded before calling factor_sqrfree(). This avoids problems
with some input polynomials.
Jens Vollinga [Mon, 10 Nov 2008 12:38:11 +0000 (13:38 +0100)]
Added modular square free factorization.
Completed distinct degree factorization.
Univariate polynomial factorization uses now upoly.
Merged class Partition and function split into class factor_partition.
Jens Vollinga [Wed, 5 Nov 2008 10:22:19 +0000 (11:22 +0100)]
Changed code from using cl_UP_MI to using umodpoly. The cl_UP_MI interface was
inconvenient to use and caused several very difficult bugs (some were still
unresolved before changing to umodpoly!).
Fixed severe bug in multivariate factorization. The condition that all modular
factors should be relatively prime in the base ring was violated. This caused
the factorization sometimes to go into an infinite loop.
symbol: remove return_type_tinfo() and return_type() (shrink symbol by 8 bytes)
Bloating symbol to the state sizeof(symbol) == 64 is not appreciated at all.
If someone needs/wants non-commutative symbols or other fancy stuff, please
derive from class symbol and do whatever you want.
return_type_tinfo() is used in order to distingish between non-commutative
objects of different type. However, often it's necessary to distingish
between non-commutative objects of the same type, for example, between
gamma matrices with different representation label. return_type_tinfo()
does not provide any clean way to do that. Hence, introduce return_type_t
type which holds representation label along with type information, and
use it for return_type_tinfo().