On Dec 12, 2007 9:00 AM, Martin Albrecht <[EMAIL PROTECTED]> wrote:
>
> On Wednesday 12 December 2007, William Stein wrote:
> > On Dec 12, 2007 8:47 AM, Joel B. Mohler <[EMAIL PROTECTED]> wrote:
> > > On Wed, Dec 12, 2007 at 10:12:07AM -0600, Daniel R. Grayson wrote:
> > > > The Macaulay2 code for factorization just calls the Singular-Factory
> > > > library that we link against. So I don't know why it should be faster
> > > > than Singular!
> > >
> > > I did wonder about that. However, my testing indicates that M2
> > > consistently takes about 1.5 seconds for this factorization, but singular
> > > takes much longer (and the time isn't consistent at all). I'd love to
> > > know if I'm just being tremendously stupid and missing something -- so it
> > > might be good if somebody else wants to confirm or deny my timings for
> > > the example at
> > > http://sage.math.washington.edu/home/jbmohler/singular/factorize01.singul
> > >ar
> > >
> > > How do you link against singular? I thought the singular library mode
> > > was new in the last 6 months and that M2 was around longer than that?
>
> They are linking against two components of singular (libfac and libcf) what we
> shared libraries before. We used to link against those too but they only
> provide factorization and stuff. My first guess for the different timings: M2
> uses an older version of libcf/libfac so maybe the last optimizations for us
> broke something?
I have to strongly reiterate that the above sort of work-arounds are totally the
wrong approach to this problem. The problem is that like 50% of CAS's
-- including
everything Sage uses -- use completely the wrong algorithm, and instead we
should find the right algorithm to use, then implement it ourselves.
I think I just realized what the right algorithm is and that it should do Joel's
challenge problem in about 0.03 seconds on my laptop instead of "1 to 3 minutes"
like it currently takes.
Anyway, Yi just made this planet sage blog thing:
http://sage.math.washington.edu/home/yqiang/psage/output/
So I looked at it, and saw Martin Albrecht's blog posts, which I
hadn't looked at before. Those led me to the Giac site, where
I read a talk about Giac and realized it does multivariate factorization,
and has the *ugliest* screenshots known to man...
Anyway, at the Giac site I found this paper by the author of Giac:
"A correct proof of the heuristic GCD algorithm"
For some reason I clicked on it and read the first 2 pages.
It says that to find the GCD of P(x1,...,xn) and Q(x1,...,xn),
the standard trick used by the closed source programs -- which
they don't prove (!) -- is to choose a prime p (say), and compute
GCD( P(x1,...,x_{n-1}, p), Q(x1,...,x_{n-1},p) )
recursively, then view the result as a polynomial in x_1,...,x_{n-1}
with coefficients in the p-adic numbers Z_p. Then take those
p-adics expansions (the coefficients) and turn them back into
polynomials in x_n. Duh -- that's a completely obvious trick,
that you want to use if you have fast integer arithmetic at your disposal.
He proves that if p is bigger than an explicit bound, etc., then the
GCD one gets this way is correct. For Joel's sort of problems,
where his polynomial is in 2 or 3 variables or so, the above is
_possibly better_ than anything Singular does, as far as I can tell.
Or maybe it is the same thing. I don't know.
Anyway, I checked what Joel really wanted, and he wants to
*factor* rather than compute GCD's, at least in this thread.
It takes about 3 seconds to realize that the same ideas as
in the above GCD algorithm also works for polynomial factorization.
And Joel has a very simple polynomial in *three* variables, for
which Singular takes minutes to factor it, but Magma is very
fast (< 1 second). So I took his polynomial, set one of the
variables equal to 97, factored the polynomial that is now
in 2 variables -- which took 0.01 seconds, then did
the above p-adic reconstruction trick, and got the correct
factorization.
I'm too busy with other things to implement this algorithm. It seems
as though Joel or somebody should implement this in Sage, so that
it is available as an option, since in some cases it is vastly better than
whatever Singular does.
Here's the Sage worksheet I used when playing around as above:
{{{id=171|
R.<p,X1,X2> = QQ[]
f =
-p^170*X1^10*X2^10+p^130*X1^10*X2^5+p^130*X1^5*X2^10-p^90*X1^5*X2^5+p^80*X1^5*X2^5-p^40*X1^5-p^40*X2^5+1
}}}
{{{id=172|
S.<X1,X2> = QQ[]
}}}
{{{id=178|
}}}
{{{id=170|
g = R.hom([97, X1, X2])(f); g
///
-56389056138356198216968526048767543091860801642367637717432942146382994680778455554250110793949802034049057298447017689263742993445880143642591140366642673142083758893391683798980349738411985421075076048312108908602231196435266097495925028306158172875468837901220347199073347808740182578622261748619920748084199390951601879629728968102849*X1^10*X2^10
+
1906889180504224276009953368161870579887188681832408838389939140200952007749109183799089214892653951940639044911734149798645662684504564746496022865286076242690037067023391046184651580740857887451370800736757437318895147621635439661232557442347754978260243649*X1^10*X2^5
+
1906889180504224276009953368161870579887188681832408838389939140200952007749109183799089214892653951940639044911734149798645662684504564746496022865286076242690037067023391046184651580740857887451370800736757437318895147621635439661232557442347754978260243649*X1^5*X2^10
-
64484610946532362921199276257300627575251780870138908996004274535399763240274543593565033337620873721598512321686420643087999739313741352977226275372397370643506329868496594982848*X1^5*X2^5
-
29571228739913279498948663112125511966004786122340140758535526194999034006060801*X1^5
-
29571228739913279498948663112125511966004786122340140758535526194999034006060801*X2^5
+ 1
}}}
{{{id=174|
time h = g.factor(); h
///
(-1) * (7837433594376961*X2 - 1) * (7837433594376961*X1 - 1) *
(577951262543040979328274795306257089*X1*X2 - 1) *
(3773075507922571775273744352333525624409237957166570157859261441*X2^4
+ 481417221911723703707680699904910046003049191681*X2^3 +
61425365346268570446197767595521*X2^2 + 7837433594376961*X2 + 1) *
(3773075507922571775273744352333525624409237957166570157859261441*X1^4
+ 481417221911723703707680699904910046003049191681*X1^3 +
61425365346268570446197767595521*X1^2 + 7837433594376961*X1 + 1) *
(111574478897742849177136197511285442432661298716331706733204593731410300884557578113685150976499394580506142057779342816269560527771563600874241*X1^4*X2^4
+
193051708905011198180170922822981440543185324376355040083133132815345376288225583305326871203951918668795969*X1^3*X2^3
+
334027661875095082935878165689897695724025886105865253311336914562753921*X1^2*X2^2
+ 577951262543040979328274795306257089*X1*X2 + 1)
CPU time: 0.02 s, Wall time: 0.02 s
}}}
{{{id=169|
h[1][0]
///
7837433594376961*X2 - 1
}}}
{{{id=168|
Qp(97)(7837433594376961)
///
97^8 + O(97^28)
}}}
{{{id=175|
h[4][0]
///
3773075507922571775273744352333525624409237957166570157859261441*X2^4
+ 481417221911723703707680699904910046003049191681*X2^3 +
61425365346268570446197767595521*X2^2 + 7837433594376961*X2 + 1
}}}
{{{id=177|
Qp(97)(3773075507922571775273744352333525624409237957166570157859261441)
///
97^32 + O(97^52)
}}}
{{{id=176|
>>>
correct answer, from REDUCE:
2:
factorize(-p^170*X1^10*X2^10+p^130*X1^10*X2^5+p^130*X1^5*X2^10-p^90*X1^5*X2^5+p^80*X1^5*X2^5-p^40*X1^5-p^40*X2^5+1)
2: ;
72 4 4 54 3 3 36 2 2 18
{{ - (p *x1 *x2 + p *x1 *x2 + p *x1 *x2 + p *x1*x2 + 1),
1},
32 4 24 3 16 2 8
{p *x1 + p *x1 + p *x1 + p *x1 + 1,1},
32 4 24 3 16 2 8
{p *x2 + p *x2 + p *x2 + p *x2 + 1,1},
18
{p *x1*x2 - 1,1},
8
{p *x1 - 1,1},
8
{p *x2 - 1,1}}
}}}
--~--~---------~--~----~------------~-------~--~----~
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-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~----------~----~----~----~------~----~------~--~---
