A little benchmark

[Waldek Hebisch]

Sage folks frequently boast how fast their base routines are.
So I decided to compare our new modular routine with Sage.
I used Sage 4.1.1 (few months old but IIUC there were no
significant speedups for routines in question) and FriCAS
trunk (the routines in 1.0.8 are essentially the same).
Imporatant point is that FriCAS was compiled using sbcl
and '--with-algebra-optimization="((speed 3) (safety 0))"'
option, that is algebra was compiled at safety 0. For low
level routines safety 0 gives 2-4 times code then regular
build. Given that core Sage routines are in C/assembler
there is no extra safety checks, so I think that using safety
0 is the only fair way to do comparison. I compared
seed of mutiplcation and gcd for polynomials over Z_2500009.

For multiplication I computed square of (x+5)^n, that is
both arguments were the same (FriCAS does not try to take
advantage from this, Sage probably also not). For gcd
I used (((x+5)*(x+7))^5)^(n/10) and (((x+5)*(x+3))^5)^(n/10),
so that gcd has degree n/2. For Sage I set

R.<x> = PolynomialRing(GF(2500009))

then computed polynomials and used %timeit to get timings.
For FriCAS I used the following functions:

p_to_u32(p) ==
pu := p :: UnivariatePolynomial(x, Integer)
n := degree(pu)
res := new(n+1, 0)$U32Vector
for i in 0..n repeat
res(i) := positiveRemainder(coefficient(pu, i), 2500009)
res

fun2(p, n) ==
#[mul(p, p, 2500009)$U32VectorPolynomialOperations for i in 1..n]

fun3(p1, p2, n) ==
#[gcd(p1, p2, 2500009)$U32VectorPolynomialOperations for i in 1..n]

and timed apropriate number of iterations.

The results are as follows:

Sage 4.1.1 FriCAS trunk
safety 0
* gcd * gcd
deg 10 9.78 30.1 5.4 12.3
deg 20 10.4 44.1 8.3 19.4
deg 30 10.8 58.7 11.2 31.2
deg 40 11.4 74.7 16.4 48.3
deg 50 12 93 21.9 67.0
deg 60 12.7 112 27.5 91.9
deg 70 13.3 209 35.5 120
deg 80 14.1 265 43.7 153
deg 90 14.7 321 54.5 189
deg 100 15.9 351 63.5 230


So for small degree we are faster both at multiplication and
at gcd. For slightly larger degree Sage is faster at multiplication,
getting about 4 times faster at degree 100. For gcd we are
faster up to degree 100, but Sage scales better.

I did not try higher degrees for three reasons:
- low degrees are most important for symbolic computations
- our current routines are quadratic, so at higher degrees
we will be slower for obvious reasons
- at high degree it makes sense to switch to asymptitically
faster routine, possibly in external library. At low
degree one has to pay interfacing cost, and one point of
this exercise was to determine possible treshhold for
such switch.

I must admit that I a bit surprised, given the whole speed propaganda
surrounding Sage I expected it to perform better. I know that our
routines are far from "best possible" -- ATM we stay within Spad and
are limited by sbcl code generator. Also, there are some algorithmic
improvement which give better speed, but complicate code (I wanted
reasonably quick code which is reasonable fast to write, leaving
more to the future).

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Waldek Hebisch]

Another little benchmark.

p1 := (1 + x + y + 2*z^2 + 3*t^3 + 5*u^5)^12;
p2 := (1 + u + t + 2*z^2 + 3*y^3 + 5*x^5)^12;

res := p1*p1;

FriCAS (trunk, sbcl, safety 0): 34.72 sec
Sage 4.1.1: 44.28 sec

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Martin Rubey]

Waldek Hebisch <heb...@math.uni.wroc.pl> writes:

> For Sage I set
>
> R.<x> = PolynomialRing(GF(2500009))
>
> then computed polynomials and used %timeit to get timings.
> For FriCAS I used the following functions:
>
> p_to_u32(p) ==
> pu := p :: UnivariatePolynomial(x, Integer)
> n := degree(pu)
> res := new(n+1, 0)$U32Vector
> for i in 0..n repeat
> res(i) := positiveRemainder(coefficient(pu, i), 2500009)
> res
>
> fun2(p, n) ==
> #[mul(p, p, 2500009)$U32VectorPolynomialOperations for i in 1..n]
>
> fun3(p1, p2, n) ==
> #[gcd(p1, p2, 2500009)$U32VectorPolynomialOperations for i in 1..n]

Isn't this benchmark unfair to sage? I would think that the natural
counterpart to R.<x> = PolynomialRing(GF(2500009)) in fricas would be
SUP PF 2500009.

I have no idea, but I wouldn't be surprised if you could speed up things
in sage if you call appropriate library functions directly?

On the other hand, I guess it wouldn't be too hard to incorporate
p_to_u32 into SUP, i.e., special case small prime fields.

Martin

[Waldek Hebisch]

Matin Rubey wrote:
> Waldek Hebisch <heb...@math.uni.wroc.pl> writes:
>
> > For Sage I set
> >
> > R.<x> = PolynomialRing(GF(2500009))
> >
> > then computed polynomials and used %timeit to get timings.
> > For FriCAS I used the following functions:
> >
> > p_to_u32(p) ==
> > pu := p :: UnivariatePolynomial(x, Integer)
> > n := degree(pu)
> > res := new(n+1, 0)$U32Vector
> > for i in 0..n repeat
> > res(i) := positiveRemainder(coefficient(pu, i), 2500009)
> > res
> >
> > fun2(p, n) ==
> > #[mul(p, p, 2500009)$U32VectorPolynomialOperations for i in 1..n]
> >
> > fun3(p1, p2, n) ==
> > #[gcd(p1, p2, 2500009)$U32VectorPolynomialOperations for i in 1..n]
>
> Isn't this benchmark unfair to sage? I would think that the natural
> counterpart to R.<x> = PolynomialRing(GF(2500009)) in fricas would be
> SUP PF 2500009.
>
> I have no idea, but I wouldn't be surprised if you could speed up
> things
> in sage if you call appropriate library functions directly?

PolynomialRing(GF(2500009)) is supposed to use routines from FLINT
which are advertised as very fast (at least that is my reading of
documetation). My understanding is that PolynomialRing uses
dense representation, so comparing with dense look fair. I leave
to Sage folks more detailed investigation of Sage result. I do
not know if there are any faster routines available from command
line. Timings for multiplication suggest that there is significant
overhead (possibly coming from Python interpreter) beyond
multiplication. OTOH overhead for multiplication and for GCD
should be the same, but Sage gcd takes significantly more time
than multiplication, which suggests that at least for GCD core
routines are slower.

BTW: I wonder how many systems have polynomials as general as
SparseUnivariatePolynomial -- I think that in most systems
'x^(10^100)' would fail.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Martin Rubey]

Waldek Hebisch <heb...@math.uni.wroc.pl> writes:

> BTW: I wonder how many systems have polynomials as general as
> SparseUnivariatePolynomial -- I think that in most systems
> 'x^(10^100)' would fail.

I just noticed that in sage one has to use an option:

sage: R.<y> = PolynomialRing(QQ, sparse=True); R
Sparse Univariate Polynomial Ring in y over Rational Field
sage: y^(10^10)
y^10000000000

Martin