search

sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.46k stars 485 forks source link

series yields wrong result #12589

Open dkrenn opened 12 years ago

dkrenn commented 12 years ago

The following was posted on the public bug reports from the notebook interface by Clemens Heuberger on 1/4/2011 and on sage-support. In the latter other code examples were posted.

f.series(q,2) for the f defined below yields Order(q^2) which is incorrect, as f.subs(q=0) equals 1 (which is correct).

sage: var('q')
sage: f=(q^13362120470/((q - 1)*(q^5 - 1)*(q^21 - 1)*(q^85 - 1)*(q^341 -
1)*(q^1365 - 1)*(q^5461 - 1)*(q^21845 - 1)*(q^87381 - 1)*(q^349525 -
1)*(q^1398101 - 1)*(q^5592405 - 1)*(q^22369621 - 1)*(q^89478485 -
1)*(q^357913941 - 1)*(q^1431655765 - 1)*(q^5726623061 - 1)) +
1)/(q^7635497409/((q - 1)*(q^5 - 1)*(q^21 - 1)*(q^85 - 1)*(q^341 -
1)*(q^1365 - 1)*(q^5461 - 1)*(q^21845 - 1)*(q^87381 - 1)*(q^349525 -
1)*(q^1398101 - 1)*(q^5592405 - 1)*(q^22369621 - 1)*(q^89478485 -
1)*(q^357913941 - 1)*(q^1431655765 - 1)*(q^5726623061 - 1)) + 1)
sage: f.subs(q=0)
1
sage: f.series(q,2)

CC: @cheuberg @rwst

Component: symbolics

Keywords: series, order, symbolics

Stopgaps: todo

Issue created by migration from https://trac.sagemath.org/ticket/12589

dkrenn commented 12 years ago

Description changed:

--- 
+++ 
@@ -1,4 +1,4 @@
-The following was posted by Clemens Heuberger on [the public bug reports from the notebook interface](https://spreadsheets.google.com/pub?key=pCwvGVwSMxTzT6E2xNdo5fA) and on [sage-support](http://groups.google.com/group/sage-support/browse_thread/thread/2949e7c3d0e84ef/852fabf0229526ed). In the latter other code examples were posted.
+The following was posted on [the public bug reports from the notebook interface](https://spreadsheets.google.com/pub?key=pCwvGVwSMxTzT6E2xNdo5fA) by Clemens Heuberger on 1/4/2011 and on [sage-support](http://groups.google.com/group/sage-support/browse_thread/thread/2949e7c3d0e84ef/852fabf0229526ed). In the latter other code examples were posted.

 `f.series(q,2)` for the `f` defined below yields `Order(q^2)` which is incorrect, as `f.subs(q=0)` equals `1` (which is correct).
dkrenn commented 9 years ago
comment:6

Still there in 6.6.

ea1d0bf8-c27a-4548-8cb7-de0b1d02441a commented 9 years ago

Stopgaps: todo

rwst commented 9 years ago
comment:9

This looks like an int overflow problem. A more minimal case is:

sage: f=(q^13362120470/((q - 1)) +1)/(q^7635497409/((q - 1)) + 1)
sage: f.series(q,2)
Order(q^2)

sage: f=(q^1336212047/((q - 1)) +1)/(q^763549741/((q - 1)) + 1)
sage: f.series(q,2)
1 + Order(q^2)

EDIT: Note that log_2(7635497409) = 32.8...

rwst commented 8 years ago
comment:10

Note that in Pynac/GiNaC the degree/ldegree member functions of basic returns an int and since that is virtual =0, all other objects do as well. So changing this would be major enhancement with performance repercussions.

dkrenn commented 8 years ago
comment:11

Replying to @rwst:

Note that in Pynac/GiNaC the degree/ldegree member functions of basic returns an int and since that is virtual =0, all other objects do as well. So changing this would be major enhancement with performance repercussions.

Being stucked with this bug is not an option either... What else can we do?

rwst commented 8 years ago
comment:12

Replying to @dkrenn:

Being stucked with this bug is not an option either... What else can we do?

Mathematically you are dealing with so-called "lacunar" or "super-sparse" series. Changing from 32bit to 64bit does not resolve it in principle if it were possible with existing Sage functionality, try this:

sage: R.<a>=QQ[[]]
sage: 2^33
8589934592
sage: a^8589934592-1

versus

sage: R.<x>=PolynomialRing(QQ,sparse=True)
sage: 2^65
36893488147419103232
sage: x^36893488147419103232-1
x^36893488147419103232 - 1

versus

sage: R.<x>=PowerSeriesRing(QQ,sparse=True)
sage: x^36893488147419103232-1
-1 + x^36893488147419103232
sage: x^36893488147419103232-1+O(x^5555555555555555555555)
OverflowError: Python int too large to convert to C long

sage: R.<q>=PowerSeriesRing(QQ,sparse=True)
sage: f=(q^13362120470/((q - 1)) +1)/(q^7635497409/((q - 1)) + 1)
MemoryError

so apparently Sage's sparse power series ring can represent monomials with bigint degree but not a bigint series order term. Representing or computing a full series lastly can involve memory problems for whatever reason.

However, if you request f.series(q,2) you do not want a full series so maybe lazy series could be your solution.

rwst commented 8 years ago
comment:13

I wrote earlier:

Note that in Pynac/GiNaC the degree/ldegree member functions of basic returns an int and since that is virtual =0, all other objects do as well. So changing this would be major enhancement with performance repercussions.

Maybe it is possible to check the expression for bigint exponents and, if so, call a special algorithm like lazy series. There are likely to be special cases of series expansion algorithms in the future, anyway, e.g. a call to fast univariate expansion via flint.

nbruin commented 8 years ago
comment:14

Replying to @rwst:

Note that in Pynac/GiNaC the degree/ldegree member functions of basic returns an int and since that is virtual =0, all other objects do as well. So changing this would be major enhancement with performance repercussions.

Presently, we are silently getting wrong results, though. A good step forward would be if we'd get an "overflow error" rather than a series expansion that seems to have 0-terms in order 0, 1.

rwst commented 8 years ago

Description changed:

--- 
+++ 
@@ -17,3 +17,4 @@
 sage: f.subs(q=0)
 1

+The pynac ticket implementing an overflow check is https://github.com/pynac/pynac/issues/121. When this is done this ticket should doctest it.

rwst commented 8 years ago
comment:16

Replying to @dkrenn:

Being stucked with this bug is not an option either... What else can we do?

We could adapt a polynomial/series package that supports unlimited exponents: http://bluescarni.github.io/piranha/sphinx/polynomials.html

rwst commented 7 years ago
comment:17

The behaviour of this test case changed somewhat. With Sage-8.1.beta6 it takes minutes to output (-1)*q^(-1) + 1 + Order(q^2). The time is spent in FLINT manipulating huge dense polynomials. The change is because Pynac no longer uses differentiation for these kind of series. It looks like the way to go here is therefore to 1. recognize rational expressions; 2. normalize them (this needs a new algorithm for sparse polynomials); 3. add a new algorithm to compute the sparse series. As to 2. this is also needed for correct handling of the test case of #23925.

rwst commented 7 years ago
comment:18

Replying to @dkrenn:

Replying to @rwst:

Note that in Pynac/GiNaC the degree/ldegree member functions of basic returns an int and since that is virtual =0, all other objects do as well. So changing this would be major enhancement with performance repercussions.

Being stucked with this bug is not an option either... What else can we do?

Actually, Pynac master now supports taking any non-symbolic real degree. I'm not sure atm how this can be used to make progress with this ticket.

rwst commented 7 years ago
comment:19

Replying to @rwst:

The behaviour of this test case changed somewhat. With Sage-8.1.beta6 it takes minutes to output (-1)*q^(-1) + 1 + Order(q^2). The time is spent in FLINT manipulating huge dense polynomials. The change is because Pynac no longer uses differentiation for these kind of series. It looks like the way to go here is therefore to 1. recognize rational expressions; 2. normalize them (this needs a new algorithm for sparse polynomials); 3. add a new algorithm to compute the sparse series. As to 2. this is also needed for correct handling of the test case of #23925.

However 2. can only be resolved with a dedicated sparse polynomial package. Actually if 1. #23925 is done; and 2. Pynac normalizes rational functions before developing series; then this ticket can be considered done because we will get a SIGABRT exception because of memory error from Singular+FLINT. This would avoid setting a hard limit on exponents, shifting the responsibility to the polynomial package.

rwst commented 7 years ago

Description changed:

--- 
+++ 
@@ -12,9 +12,7 @@
 1)*(q^1365 - 1)*(q^5461 - 1)*(q^21845 - 1)*(q^87381 - 1)*(q^349525 -
 1)*(q^1398101 - 1)*(q^5592405 - 1)*(q^22369621 - 1)*(q^89478485 -
 1)*(q^357913941 - 1)*(q^1431655765 - 1)*(q^5726623061 - 1)) + 1)
-sage: f.series(q,2)
-Order(q^2)
 sage: f.subs(q=0)
 1
+sage: f.series(q,2)

-The pynac ticket implementing an overflow check is https://github.com/pynac/pynac/issues/121. When this is done this ticket should doctest it.

rwst commented 7 years ago
comment:21

There is need for a hard limit.

FLINT allows only long exponents in functions accessing their fmpz_t polynomials, so there is a size restriction when using FLINT, i.e. in Pynac series and polynomial manipulation via Singular (which uses FLINT as default).

Note that expansion of polynomials only uses Singular above a certain size, so e.g. expansion of small products with less than 400 terms overall works fine:

sage: ((1+x^(2^100))*(1-x^(2^100))).expand()
-x^2535301200456458802993406410752 + 1
sage: 2^101
2535301200456458802993406410752
rwst commented 7 years ago
comment:22

With pynac-0.7.12 we have now:

sage: f.series(q,2)
Exception (FLINT memory_manager). Unable to allocate memory.
...
RuntimeError: Aborted

sage: (1/(x^3689348814741910323-1)).series(x)
...
RuntimeError: exponent too big

Note that f.series(q,2) crashes inside the step where rational functions are normalized. This step is however not necessary in principle, and a smarter series functionality is needed to resolve this ticket. At least there is an error message now.