search

sagemath / sage

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

add computation of swinnerton-dyer polynomials to sage #9492

Open williamstein opened 14 years ago

williamstein commented 14 years ago

Magma has 'em, so we should to:

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

Sympy has them: http://github.com/mattpap/sympy-polys/commit/dc42fd1995c48ad426b6279d3d1914f74e0c3296

This page has a Mathematica notebook with a function that computes them:

http://mathworld.wolfram.com/Swinnerton-DyerPolynomial.html

I tried it in Mathematica 7, and it is massively, dramatically SLOW compared to Magma. For comparison, the 5th one takes 55 seconds in Mathematica, and in Magma it takes... 0.02 seconds. 

  n   time in seconds with Magma 2.15.11 on my macbook air
----------------
  5  | 0.02
  6  | 0.18
  7  | 6.68
  8  | 99.99

CC: @mstreng

Component: misc

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

williamstein commented 14 years ago

Description changed:

--- 
+++ 
@@ -3,7 +3,7 @@
    http://magma.maths.usyd.edu.au/magma/htmlhelp/text306.htm

 Sympy has them:
-   
+   http://github.com/mattpap/sympy-polys/commit/dc42fd1995c48ad426b6279d3d1914f74e0c3296

 This page has a Mathematica notebook with a function that computes them:
williamstein commented 14 years ago
comment:1

OK, I just wrote some Sage code from scratch that provably computes the swinnerton-dyer polynomials vastly more quickly than Magma or Mathematica.

williamstein commented 14 years ago
comment:2

See also http://wiki.sagemath.org/days23/CodingProjects

In the attached file sd.sage, I implemented a few functions, and also copied in Jeroen Demeyer's p-adic function. I think the interval arithmetic one (sdpoly3) is considerably faster than the p-adic one, and moreover sdpoly3 is provably correct. It simply computes the poly using intervals until there is a unique integer in each interval. So I think it should just be SAge-library-ified and put into sage... somewhere.

williamstein commented 14 years ago
comment:3

Attachment: sd.sage.gz

I think these are interesting because they are used in benchmarks for factoring and irreducibility testing. See http://www.shoup.net/ntl/doc/tour-time.html

The sdpoly3 function computes the 10th Swinnerton-Dyer poly of degree 2^10=1024 in < 20 seconds. 

sage: time c = sdpoly3(10)
CPU times: user 17.02 s, sys: 0.05 s, total: 17.07 s
Wall time: 17.50 s
mstreng commented 14 years ago
comment:4

Attachment: sd2.sage.gz

I just adapted sdpoly3 to use a binary tree. The result is called sdpoly5, see file sd2.sage. I found sdpoly5 to be slightly faster than sdpoly3 in the tests that I have run.

When using naive polynomial multiplication, the algorithms sdpoly3 and sdpoly5 are asymptotically equivalent. As soon as FFT quasi-linear polynomial multiplication is implemented for interval arithmetic and examples become large, the algorithm sdpoly5 should be the faster one (quasi-linear).

sage: time sdpoly5(12)
CPU times: user 845.82 s, sys: 0.10 s, total: 845.92 s
Wall time: 846.13 s

sage: time sdpoly3(12)
CPU times: user 861.84 s, sys: 0.01 s, total: 861.85 s
Wall time: 861.98 s
williamstein commented 14 years ago
comment:5

I just adapted sdpoly3 to use a binary tree.

Thanks. But just a remark -- the Sage prod command already uses a binary tree, at least if the input is a list or of length < 1000. See the file

devel/sage/sage/misc/misc_c.pyx

which Robert Bradshaw wrote. So putting 

  prod(list( ... ))

instead of 

  prod( ...)

in sdpoly3 might be another approach.

williamstein commented 14 years ago
comment:6

I'm at a Singular conference, so I decided to try this problem using Singular polynomial quotient rings. It's pretty good, though it doesn't beat interval arithmetic speed-wise.

def sdpoly6(n):
    R = PolynomialRing(QQ,n+1,names='x')
    x = R.gens()
    v = primes_first_n(n)
    I = R.ideal([ x[i]^2-v[i] for i in range(len(v)) ])
    S = R.quotient(I)
    x = S.gens()
    C = cartesian_product_iterator([[-1,1]]*n)
    f = prod([ x[-1] + sum(s[i]*x[i] for i in range(n)) for s in C])
    return f

Some timings:

sage: time a = sdpoly6(8)
Time: CPU 0.71 s, Wall: 0.71 s
sage: time a = sdpoly6(9)
Time: CPU 3.44 s, Wall: 3.47 s
sage: time a10 = sdpoly6(10)
Time: CPU 29.03 s, Wall: 29.19 s

Very impressive for something non-numerical, IMHO...

fredrik-johansson commented 11 years ago
comment:7

FLINT now includes a function that computes the nth Swinnerton-Dyer polynomial. It could be wrapped trivially.

fchapoton commented 10 years ago
comment:10

I have tried to wrap the flint function, but was not able to. Could someone more versed into interfaces do that ?