#10042: Doctest failure in sage/rings/polynomial/polynomial_element.pyx
-------------------------------+--------------------------------------------
Reporter: mpatel | Owner: mvngu
Type: defect | Status: needs_work
Priority: blocker | Milestone: sage-4.6
Component: doctest | Keywords:
Author: | Upstream: N/A
Reviewer: | Merged:
Work_issues: report upstream |
-------------------------------+--------------------------------------------
Comment(by leif):
Replying to [comment:44 zimmerma]:
> [...] it seems the computation is done by the lapack {{{dgeev}}}
routine. From the source file
> {{{dgeev.f}}} in the Sage distribution, it does not call {{{POW}}}
however it calls {{{SQRT}}}. But the square root is a standard IEEE
function. What could happen however is that on some machines lapack could
use double extended precision (i.e., 64-bit instead of 53-bit
significand).
Yes, the difference is caused by LAPACK, not !NumPy (see below).
I'm not sure but I doubt PPCs have (Intels!) extended double precision,
i.e. {{{-mfpmath=sse}}} vs. {{{-mfpmath=387}}} doesn't hold there. (But
IBM uses other floating-point formats beyond 64-bit double precision.)
> One way to debug this is to add print statements in the {{{dgeev.f}}}
file, recompile lapack and see the difference between two different
computers. How can one (easily) recompile and build a spkg (I guess
{{{sage -br}}} will not do that)?
For hacking the Fortran code (avoiding to create new LAPACK spkgs), it is
sufficient to extract the corresponding files (e.g. {{{dgeev.f}}}) from
the spkg and put modified versions into the [static] LAPACK library, which
is {{{$SAGE_ROOT/local/lib/liblapack.a}}}, and to afterwards rebuild /
reinstall the !NumPy spkg, which links this library into its dynamic ones.
Assuming one has copied {{{lapack-20071123.p1/src/SRC/dgeev.f}}} to
{{{$SAGE_ROOT}}}, one can do e.g. the following to include a modified
version:
{{{
#!sh
~/sage-x.y.z$ gfortran -fPIC -c dgeev.f # or
sage_fortran -c dgeev.f
~/sage-x.y.z$ ar rs local/lib/liblapack.a dgeev.o # "s" saves
running ranlib separately
~/sage-x.y.z$ ./sage -f spkg/standard/numpy-1.3.0.p4.spkg # force
reinstallation
}}}
and then run Sage.
(Dave: I've ''not'' verified if the "s" option is POSIX... ;-) The
{{{ar}}} command above is just a short-cut for
{{{
#!sh
~/sage-x.y.z$ ar r local/lib/liblapack.a dgeev.o # replace the existing
dgeev module
~/sage-x.y.z$ ranlib local/lib/liblapack.a # rebuild the symbol
table
}}}
)
On a 32-bit Pentium4 (Prescott) and a 64-bit Core2 I get, respectively:
{{{
sage: R.<x>=RealField(200)[]
sage: f=x^2 - R(pi)
sage: import numpy
sage: na=numpy.array([float(c) for c in reversed(f.list())],
dtype="double")
sage: n_roots=numpy.roots(na)
Workspace query only:
A( 1, 1)= -0.0000000000000000000000000D+00
A( 1, 2)= 0.3141592653589793115998000D+01
A( 2, 1)= 0.1000000000000000000000000D+01
A( 2, 2)= 0.0000000000000000000000000D+00
Actual computation:
A( 1, 1)= -0.0000000000000000000000000D+00
A( 1, 2)= 0.3141592653589793115998000D+01
A( 2, 1)= 0.1000000000000000000000000D+01
A( 2, 2)= 0.0000000000000000000000000D+00
Eigenvalues only...
WR( 1)= -0.1772453850905515881919400D+01
WR( 2)= 0.1772453850905516103964000D+01
sage:
}}}
{{{
sage: R.<x> = RealField(200)[]
sage: f = x^2 - R(pi)
sage: import numpy
sage: na=numpy.array([float(c) for c in reversed(f.list())],
dtype='double')
sage: n_roots=numpy.roots(na)
Workspace query only:
A( 1, 1)= -0.0000000000000000000000000D+00
A( 1, 2)= 0.3141592653589793115998000D+01
A( 2, 1)= 0.1000000000000000000000000D+01
A( 2, 2)= 0.0000000000000000000000000D+00
Actual computation:
A( 1, 1)= -0.0000000000000000000000000D+00
A( 1, 2)= 0.3141592653589793115998000D+01
A( 2, 1)= 0.1000000000000000000000000D+01
A( 2, 2)= 0.0000000000000000000000000D+00
Eigenvalues only...
WR( 1)= -0.1772453850905515881919400D+01
WR( 2)= 0.1772453850905515881919400D+01
sage:
}}}
(The machine parameters queried by LAPACK's {{{DLAMCH()}}} are identical;
the chosen precision / format for the {{{PRINT}}} statements,
{{{D35.25}}}, should be sufficient.)
What about (temporarily?) tagging the doctest result(s) with {{{#
32-bit}}} / {{{# 64-bit}}}? (It seems the results are identical on
'''all''' 32-bit platforms...?)
I've not yet tracked this down further (perhaps someone more knowledgeable
of LAPACK / Fortran can do this), nor have I tried the more recent LAPACK
version.
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/10042#comment:56>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
