multivariate polynomial multiplication, division, GCD, factorization

43 views
Skip to first unread message

parisse

unread,
Jul 9, 2010, 8:32:27 AM7/9/10
to sage-devel
I have updated the benchmarks page of giac, with comparison with magma
and trip there:
http://www-fourier.ujf-grenoble.fr/~parisse/giac/benchmarks/benchmarks.html
I'm curious to know if this will raise interest in integrating giac in
sage one day, now that giac has about the same timings than magma for
GCD...
Message has been deleted

parisse

unread,
Jul 10, 2010, 1:29:57 AM7/10/10
to sage-devel
> It does raise interest. But then there are questions of encouraging
> people to get involved in maintaining the interface between giac and
> Sage. Someone needs to step up and maintain that interface for at
> least 3 years.
>

I don't think it would be a problem, because I could probably maintain
it myself once it is done, I believe that the hard work is to
implement the initial interface.

William Stein

unread,
Jul 15, 2010, 6:28:17 PM7/15/10
to sage-...@googlegroups.com
On Fri, Jul 9, 2010 at 2:32 PM, parisse <bernard...@ujf-grenoble.fr> wrote:
> I have updated the benchmarks page of giac, with comparison with magma
> and trip there:
> http://www-fourier.ujf-grenoble.fr/~parisse/giac/benchmarks/benchmarks.html
> I'm curious to know if this will raise interest in integrating giac in
> sage one day,

Yes, it raises interest.

> now that giac has about the same timings than magma for
> GCD...

We have all run into some problems trying to build the source code you
posted here:

http://xcas.svn.sourceforge.net/viewvc/xcas/

I wonder if you might use the account I gave you on
redhawk.math.washington.edu (username: parisse) to
update giac so that it will build fine on there. It would also be
nice if you could login to redhawk.math.washington.edu (same
credentials as I gave you for sage.math -- I made an account for you),
go to /scratch/ make yourself a directory, and build giac there as
well, then checkin any changes needed to do the build to your svn
repo. Let me know what you need me to apt-get install.

Thanks!

-- William


--
William Stein
Professor of Mathematics
University of Washington
http://wstein.org

William Stein

unread,
Jul 15, 2010, 6:35:50 PM7/15/10
to sage-...@googlegroups.com
On Fri, Jul 9, 2010 at 2:32 PM, parisse <bernard...@ujf-grenoble.fr> wrote:

Hi,

Is your factoring code (and corresponding benchmarks) multithreaded?

William Stein

unread,
Jul 15, 2010, 7:28:58 PM7/15/10
to sage-...@googlegroups.com
On Fri, Jul 16, 2010 at 12:28 AM, William Stein <wst...@gmail.com> wrote:
> On Fri, Jul 9, 2010 at 2:32 PM, parisse <bernard...@ujf-grenoble.fr> wrote:
>> I have updated the benchmarks page of giac, with comparison with magma
>> and trip there:
>> http://www-fourier.ujf-grenoble.fr/~parisse/giac/benchmarks/benchmarks.html
>> I'm curious to know if this will raise interest in integrating giac in
>> sage one day,
>
> Yes, it raises interest.
>
>>  now that giac has about the same timings than magma for
>> GCD...
>
> We have all run into some problems trying to build the source code you
> posted here:
>
>     http://xcas.svn.sourceforge.net/viewvc/xcas/
...
>
> Thanks!

I just wanted to add that we had no trouble at all building a giac
package for Sage using your release 0.9.0 version of Giac. So
whatever you were doing for that build was good.

In case you're curious, the giac-0.9.0 I made is here:

http://sage.math.washington.edu/home/wstein/tmp/giac-0.9.0.spkg

It's just a bzip2'd tarball:

tar jxvf giac-0.9.0.spkg

to see inside.

parisse

unread,
Jul 16, 2010, 2:12:05 AM7/16/10
to sage-devel


> Hi,
>
> Is your factoring code (and corresponding benchmarks) multithreaded?
>

The factoring code is not multithreaded, but it calls multivariate
polynomial multiplication several times during the Hensel lift phase,
and that code is multithreaded. The timings reported are when you add
time for all threads, real time (as reported by the bash time command)
is lower on multi-CPU. I'm still optimizing the multiplication code
and also discovered a small bug which slowed down the factorization,
the timings I get now are 11s for the last example instead of 27s.
The svn repository at sourceforge is only a backup (Makefile,
configure not inside or not up to date), as you have done, the right
source is the tarball from my website. I'm on vacation with low
Internet bandwith, I will give you more feedback on your giac spkg
later, what does it do currently? Build giac inside sage but no
interface from sage prompt I guess?

William Stein

unread,
Jul 16, 2010, 2:51:39 AM7/16/10
to sage-...@googlegroups.com
On Fri, Jul 16, 2010 at 8:12 AM, parisse
<bernard...@ujf-grenoble.fr> wrote:
>
>
>> Hi,
>>
>> Is your factoring code (and corresponding benchmarks) multithreaded?
>>
>
> The factoring code is not multithreaded, but it calls multivariate
> polynomial multiplication several times during the Hensel lift phase,
> and that code is multithreaded. The timings reported are when you add
> time for all threads, real time (as reported by the bash time command)
> is lower on multi-CPU. I'm still optimizing the multiplication code
> and also discovered a small bug which slowed down the factorization,
> the timings I get now are 11s for the last example instead of 27s.
> The svn repository at sourceforge is only a backup (Makefile,
> configure not inside or not up to date), as you have done, the right
> source is the tarball from my website.

I'm *REALLY*, really, really confused about which tarball is the
source code for the
latest version of giac. Just to clarify, please send me a direct link
in email here, please.

Thanks!

> I'm on vacation with low
> Internet bandwith, I will give you more feedback on your giac spkg
> later, what does it do currently? Build giac inside sage but no
> interface from sage prompt I guess?
>

> --
> To post to this group, send an email to sage-...@googlegroups.com
> To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/sage-devel
> URL: http://www.sagemath.org

parisse

unread,
Jul 16, 2010, 8:50:55 AM7/16/10
to sage-devel
> I'm *REALLY*, really, really confused about which tarball is the
> source code for the
> latest version of giac.  Just to clarify, please send me a direct link
> in email here, please.
>

Should be there
http://www-fourier.ujf-grenoble.fr/~parisse/giac/giac_unstable.tgz
or there
ftp://ftp-fourier.ujf-grenoble.fr/xcas/giac_unstable.tgz
with the latest changes.
I realize I did not update the http link since 3 months, should be up
to date now.
On x64, you should run
export CXXFLAGS='-O2 -DSMARTPTR64 -D_I386_ '
before calling configure to get the best performances

William Stein

unread,
Jul 17, 2010, 11:13:16 AM7/17/10
to sage-...@googlegroups.com
On Fri, Jul 16, 2010 at 2:50 PM, parisse
<bernard...@ujf-grenoble.fr> wrote:
>> I'm *REALLY*, really, really confused about which tarball is the
>> source code for the
>> latest version of giac.  Just to clarify, please send me a direct link
>> in email here, please.
>>
>
> Should be there
> http://www-fourier.ujf-grenoble.fr/~parisse/giac/giac_unstable.tgz

After installing enough latex packages so that the docs seemed to
build, my build still fails here:


...
make[3]: Entering directory
`/scratch/wstein/build/sage-4.4.4/spkg/build/giac-2010-07-17/src/doc'
make install-exec-hook
make[4]: Entering directory
`/scratch/wstein/build/sage-4.4.4/spkg/build/giac-2010-07-17/src/doc'
/bin/bash /scratch/wstein/build/sage-4.4.4/spkg/build/giac-2010-07-17/src/config/install-sh
-d /scratch/wstein/build/sage-4.4.4/local/shar
e/giac/doc/dxcas
/bin/bash /scratch/wstein/build/sage-4.4.4/spkg/build/giac-2010-07-17/src/config/install-sh
-d /scratch/wstein/build/sage-4.4.4/local/shar
e/giac/doc//pari
/usr/bin/install -c -m 644 ../doc/dxcas/*
/scratch/wstein/build/sage-4.4.4/local/share/giac/doc/dxcas
/usr/bin/install -c -m 644 ../doc/pari/*
/scratch/wstein/build/sage-4.4.4/local/share/giac/doc/pari
ln -s ../aide_cas /scratch/wstein/build/sage-4.4.4/local/share/giac/doc/aide_cas
ln: creating symbolic link
`/scratch/wstein/build/sage-4.4.4/local/share/giac/doc/aide_cas': File
exists
make[4]: *** [install-exec-hook] Error 1
make[4]: Leaving directory
`/scratch/wstein/build/sage-4.4.4/spkg/build/giac-2010-07-
----


I manually deleted the relevant file

/scratch/wstein/build/sage-4.4.4/local/share/giac/doc/aide_cas

and was able to get the build to finish. So very close to "works first try".

How do I surpress printing when using giac? E.g., how do I keep giac
from printing massive amounts when I do:

8>> f1:=(x+y+z+1)^20+1; p1:=normal(f1*(f1+1));
9>> factor(p1);

By the way, I tried the above and it took only 0.63 seconds, which is
very impressive (!). It's much faster than Sage-4.4.4 and Magma.

sage: R.<x,y,z> = QQ[]
sage: f1 = (x+y+z+1)^20+1; p1=f1*(f1+1)
sage: time a = p1.factor()
CPU times: user 19.65 s, sys: 0.00 s, total: 19.65 s
Wall time: 19.68 s


sage: k=magma(p1)
sage: magma.eval('time a := Factorization(%s);'%k.name())
'Time: 9.480'

William Stein

unread,
Jul 17, 2010, 11:18:45 AM7/17/10
to sage-...@googlegroups.com
Hi,

I put my experimental giac spkg here:

http://sage.math.washington.edu/home/wstein/patches/giac-20100717.spkg

To build it you have to "sage -f -m giac-20100717". Then "sage -sh",
then delete the file it complains about (repeat, lather), type "make
install", etc.

Note that there are a *ton* of system-wide dependencies needed to get
these to build.
I installed them all on redhawk.math.washington.edu, which
sagemath-users have access
to.

Anyway, the cool thing is that this package is possible to build, and
does let one replicate Parisse's benchmarks with something built from
source (not a binary). Cool.

William

--

parisse

unread,
Jul 17, 2010, 11:49:16 AM7/17/10
to sage-devel
> After installing enough latex packages so that the docs seemed to
> build, my build still fails here:
>
> ln: creating symbolic link
> `/scratch/wstein/build/sage-4.4.4/local/share/giac/doc/aide_cas': File
> exists

Indeed, the command in the Makefile should be ln -sf instead of ln -s,
I will fix that.

> How do I surpress printing when using giac?  E.g., how do I keep giac
> from printing massive amounts when I do:
>
> 8>> f1:=(x+y+z+1)^20+1; p1:=normal(f1*(f1+1));
> 9>> factor(p1);
>

add :; at the end of a command to disable printing (in maple
compatibility mode : will also work). Here for example
f1:=(x+y+z+1)^20:;
I will try to install sage on a server somewhere to try your spkg (I'm
on a 2.6K/s internet connection), but perhaps the best is that I log
on your server and see what changes needs to be done on the giac side
so that I can generate a tar.bz2 identical to the spkg.

parisse

unread,
Jul 17, 2010, 12:02:50 PM7/17/10
to sage-devel
Hi,

It seems the original source is stored in a src subdirectory and you
added the shell script spkg-install and the readme SPKG.txt file, is
there any other changes?

William Stein

unread,
Jul 17, 2010, 12:05:23 PM7/17/10
to sage-...@googlegroups.com

Nope, that's all it is. spkg's are really simple.

-- William

parisse

unread,
Jul 19, 2010, 10:06:19 AM7/19/10
to sage-devel
Ok, I have added a link that should be updated when I make changes in
giac:
http://www-fourier.ujf-grenoble.fr/~parisse/giac/giac.spkg
(I also updated the benchmarks page
http://www-fourier.ujf-grenoble.fr/~parisse/giac/benchmarks/benchmarks.html).

William Stein

unread,
Jul 19, 2010, 11:50:46 AM7/19/10
to sage-...@googlegroups.com
On Mon, Jul 19, 2010 at 4:06 PM, parisse
<bernard...@ujf-grenoble.fr> wrote:
> Ok, I have added a link that should be updated when I make changes in
> giac:
> http://www-fourier.ujf-grenoble.fr/~parisse/giac/giac.spkg

Hi,

Please remark the spkg so when you type

tar jxvf giac.spkg

you get a directory named "giac". This is a requirement of the spkg
format, which
I forgot to mention. It's also a good idea to name the package with
a version, e.g.,
giac-20100719.spkg, in which case it would extract to giac-20100719.

I didn't expect your spkg to just work for me after the above rename,
so I wrote some stuff about debugging below. However, it turns out
your package *did* just 100% work for me! Thanks!

Could you post a very, very simple example C program that uses
libgiac.so to create a polynomial and factor it? That would be
helpful, in evaluating how giac could be used from Python most
efficiently.

--------------------------------------------------------------------------

It would be very helpful if you could:

(1) install Sage, then

(2) type
./sage -i giac.spkg

and fix the mistakes until it works.

Some hints:

- In order to make an spkg, the best thing to do is put the spkg in
a directory, say giac-20100719,
then do

sage -pkg giac-20100719

This is important. For example, your giac.spkg extracts to a
directory giac_spkg, which unfortunately won't work. The name of the
directory *must* be everything before spkg.

After fixing that, the next thing that breaks is:


- In order to try to fix this, I would type "sage -sh", then navigate to

SAGE_ROOT/spkg/build/giac-20100719/src

Doing "sage -sh" sets up the environment to be just as it was when
spkg-install was run.

I hope this helps. Don't hesitate to ask further questions, and also
note this page:

http://sagemath.org/doc/developer/producing_spkgs.html

entitled "Producing New Sage Packages".

-- William

-

> --
> To post to this group, send an email to sage-...@googlegroups.com
> To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/sage-devel
> URL: http://www.sagemath.org
>

--

parisse

unread,
Jul 19, 2010, 2:29:11 PM7/19/10
to sage-devel
> Please remark the spkg so when you type
>
>    tar jxvf giac.spkg
>
> you get a directory named "giac".    This is a requirement of the spkg
> format, which
> I forgot to mention.   It's also a good idea to name the package with
> a version, e.g.,
> giac-20100719.spkg, in which case it would extract to giac-20100719.
>
> I didn't expect your spkg to just work for me after the above rename,
> so I wrote some stuff about debugging below.  However, it turns out
> your package *did* just 100% work for me!  Thanks!
>

I renamed the package giac_0_9_0.spkg, to have a version number that
will not change too often (because I don't want to make errors when
updating the package and the HTML pages pointing to it). I'm glad it
worked for you, because it would not be easy to test it myself
currently.
Please note that the CXXFLAGS I gave are correct for x64 linux only,
is there a way to specify it in the spkg-install file (something like
#ifdef __x86_64__ in C++?

> Could you post a very, very simple example C program that uses
> libgiac.so to create a polynomial and factor it?  That would be
> helpful, in evaluating how giac could be used from Python most
> efficiently.
>
I don't have a C interface, here is one in C++ using symbolic format.
One could also directly use giac polynomial, but probably the symbolic
format is better for a Python interface. Most giac user commands have
an equivalent in C++ with a _ as prefix, a first argument that is the
argument of the giac user function or a vector of these arguments (if
multiargs) and a second argument which is a context pointer (this
makes possible to have independant sessions with one running kernel,
each session with it's own context pointer).

// -*- compile-command: "g++ -g example_factor.cc -lgiac -lgmp" -*-
#include <giac/giac.h>

using namespace std;
using namespace giac;

int main(){
cout << "Enter expression " ;
string s;
cin >> s; // or define e.g. string s("x^4-y^4")
cout << s << endl;
giac::context ct;
gen c(s,&ct);
c=eval(c,1,&ct);
cout << "Factor :" << _factor(c,&ct) << endl;
return 0;
}

parisse

unread,
Jul 19, 2010, 2:31:10 PM7/19/10
to sage-devel
I forgot to mention how to link it
g++ example_factor.cc -lgiac -lgmp -lpng

William Stein

unread,
Jul 19, 2010, 2:33:14 PM7/19/10
to sage-...@googlegroups.com
On Mon, Jul 19, 2010 at 8:29 PM, parisse
<bernard...@ujf-grenoble.fr> wrote:
>> Please remark the spkg so when you type
>>
>>    tar jxvf giac.spkg
>>
>> you get a directory named "giac".    This is a requirement of the spkg
>> format, which
>> I forgot to mention.   It's also a good idea to name the package with
>> a version, e.g.,
>> giac-20100719.spkg, in which case it would extract to giac-20100719.
>>
>> I didn't expect your spkg to just work for me after the above rename,
>> so I wrote some stuff about debugging below.  However, it turns out
>> your package *did* just 100% work for me!  Thanks!
>>
>
> I renamed the package giac_0_9_0.spkg, to have a version number that
> will not change too often (because I don't want to make errors when
> updating the package and the HTML pages pointing to it). I'm glad it
> worked for you, because it would not be easy to test it myself
> currently.

Why? You just have to:

(1) install sage
(2) type "./sage -i giac-0.9.0.spkg

Note that the package name must have the format:

package_name-version_number.spkg

Note that "-'.

If you want me to be giac user, I hope it isn't too much to ask that
you be a Sage user :-).

> Please note that the CXXFLAGS I gave are correct for x64 linux only,
> is there a way to specify it in the spkg-install file (something like
> #ifdef __x86_64__ in C++?

Not easily. It's just shell script though, so...

>> Could you post a very, very simple example C program that uses
>> libgiac.so to create a polynomial and factor it?  That would be
>> helpful, in evaluating how giac could be used from Python most
>> efficiently.
>>
> I don't have a C interface, here is one in C++ using symbolic format.

No problem, C++ is great.

> One could also directly use giac polynomial, but probably the symbolic
> format is better for a Python interface. Most giac user commands have
> an equivalent in C++ with a _ as prefix, a first argument that is the
> argument of the giac user function or a vector of these arguments (if
> multiargs) and a second argument which is a context pointer (this
> makes possible to have independant sessions with one running kernel,
> each session with it's own context pointer).
>
> // -*- compile-command: "g++ -g example_factor.cc -lgiac -lgmp" -*-
> #include <giac/giac.h>
>
> using namespace std;
> using namespace giac;
>
> int main(){
>  cout << "Enter expression " ;
>  string s;
>  cin >> s; // or define e.g. string s("x^4-y^4")
>  cout << s << endl;
>  giac::context ct;
>  gen c(s,&ct);
>  c=eval(c,1,&ct);
>  cout << "Factor :" << _factor(c,&ct) << endl;
>  return 0;
> }
>

Thanks! That will be very helpful.

parisse

unread,
Jul 19, 2010, 3:37:36 PM7/19/10
to sage-devel
> Why?  You just have to:
>
>   (1) install sage

This is not easy because I have a low Internet connection (max 2.6K/s)
and sage is huge.

> If you want me to be giac user, I hope it isn't too much to ask that
> you be a Sage user :-).
>

Of course, I will certainly install it next time I have a good
Internet connection.

> > Please note that the CXXFLAGS I gave are correct for x64 linux only,
> > is there a way to specify it in the spkg-install file (something like
> > #ifdef __x86_64__ in C++?
>
> Not easily.  It's just shell script though, so...
>

Then, until someone finds a trick to do so, I will disable the flag
for 32 bits (disabling SMARTPTR64 will cause sizeof(gen) to be 16
instead of 8 that will slow some operations because more memory will
be needed). I have also renamed giac_0...spkg to giac-0...

William Stein

unread,
Jul 19, 2010, 3:43:22 PM7/19/10
to sage-...@googlegroups.com
On Mon, Jul 19, 2010 at 9:37 PM, parisse
<bernard...@ujf-grenoble.fr> wrote:
>> Why?  You just have to:
>>
>>   (1) install sage
>
> This is not easy because I have a low Internet connection (max 2.6K/s)
> and sage is huge.

I understand. I think this explains why it would also be hard for
you to work on another remote computer, e.g.,
sage.math.washington.edu. I assume that even if you have account
there, typing is too slow to be useful.

>
>> If you want me to be giac user, I hope it isn't too much to ask that
>> you be a Sage user :-).
>>
>
> Of course, I will certainly install it next time I have a good
> Internet connection.
>
>> > Please note that the CXXFLAGS I gave are correct for x64 linux only,
>> > is there a way to specify it in the spkg-install file (something like
>> > #ifdef __x86_64__ in C++?
>>
>> Not easily.  It's just shell script though, so...
>>
>
> Then, until someone finds a trick to do so, I will disable the flag
> for 32 bits (disabling SMARTPTR64 will cause sizeof(gen) to be 16
> instead of 8 that will slow some operations because more memory will
> be needed). I have also renamed giac_0...spkg to giac-0...

Thanks!

Is there any chance you could write an email explaining the basic
ideas behind your new factorization implementation? I'm just curious
what goes into it. Even something along the lines of this page would
be nice:

http://magma.maths.usyd.edu.au/magma/htmlhelp/text315.htm

-- William

parisse

unread,
Jul 20, 2010, 2:36:46 AM7/20/10
to sage-devel

> Is there any chance you could write an email explaining the basic
> ideas behind your new factorization implementation?   I'm just curious
> what goes into it.  Even something along the lines of this page would
> be nice:
>
>    http://magma.maths.usyd.edu.au/magma/htmlhelp/text315.htm
>
>  -- William
>

It's not very different from what other systems do. First of course
square free factorization, then try a few (2) random values for all
indeterminates except the first one, for a fast check of
irreducibility. If not, lift the equality in one variable
P(x,0,...0)=product P_i(x,0,...,0)
more precisely, one must take care of the leading coefficient, hence
lift
P*lcoeff(P)^(#nfactors-1)=product P_i
where in P_i the leading coefficient is replaced by lcoeff(P).
In order to avoid densification (if lcoeff(P) has many coefficients),
I make a bivariate factorization, this way instead of using lcoeff(P)
for every P_i, I'm using a divisor of lcoeff(P). It's a little
different from what is describe in the thesis of Bernardin (I don't
make polynomial rational reconstructions for example).
There is also a try to do sparse factorization before.
If Hensel lift does not work (for example P(x,0...0) loose degree or
has non square-free factors), I'm using "heuristic factorization",
i.e. evaluate P at x>=2 linfnorm(P)+2, factor this polynomial,
reconstruct factors (using x as basis and symmetric remainder) and
check division, if remainder is 0 an irreducible factor has been
found, otherwise try with a larger value of x.
The corresponding code slices are in ezgcd.cc try_sparse_factor and
try_hensel_lift_factor, and do_factor in gausspol.cc.
Reply all
Reply to author
Forward
0 new messages