Re: How to do quickest

*To*: mathgroup at smc.vnet.net*Subject*: [mg116594] Re: How to do quickest*From*: Artur <grafix at csl.pl>*Date*: Mon, 21 Feb 2011 04:21:02 -0500 (EST)

Dear Bobby and the Rest, You was write that my procedure is black box. I'm explained that my procedure is one of the few know (published) algorhitms used to determination of proper Galois transitive group for polynomial of deg degree This is called p-adic factorization frequency or Lenstra-Czebotaryev method. This procedure is not sufficient to exact determination of Galois Group from order 8, up to order 15 is sufficient second precedure called n-set resolvent factorization pattern procedure, starting from 16 are next procedures need which work up to order 23 From 24 to 31 procedures are not yet published (according my best knowledge) but John Cannon know these and use in MAGMA programme, but not available on online calculator because exceed 1 minute. Not so far John Cannon was count exact total number of different Galois group in order 32. Best wishes Artur The recent question W dniu 2011-02-20 18:57, DrMajorBob pisze: > OK. Mathgroup wasn't listed among addresses on the e-mail I got, so > you must have sent to the group separately or BCC. > > Bobby > > On Sun, 20 Feb 2011 11:49:52 -0600, Artur<grafix at csl.pl> wrote: > >> I was do! >> Artur >> >> W dniu 2011-02-20 17:05, DrMajorBob pisze: >>> Send questions to the group, not just Daniel and myself. >>> >>> Bobby >>> >>> On Sun, 20 Feb 2011 09:04:47 -0600, Artur<grafix at csl.pl> wrote: >>> >>>> Dear Daniel and Bobby, >>>> >>>> To working Galois prcedure up to 16 degree is necessary dynamical >>>> procedure (independent from order of polynomials, n-th roots >>>> resolvend and operation + or *. I'm too poor programist yet to >>>> write these. Could You help me? >>>> >>>> In particular case sample procedure is >>>> order dynamic, 2 root resolvents fixed, operation + fixed (I need >>>> both three dynamic) >>>> (*2 root resolvent +*) >>>> p = x^9 - 9*x^7 - 21*x^6 + 72*x^5 + 99*x^4 - 99*x^3 - 585*x^2 + >>>> 549*x + >>>> 166;; deg = Length[CoefficientList[p, x]] - 1; k = >>>> N[x /. Solve[p == 0, x], 300]; pol = 1; aa = {}; Do[ >>>> Do[pol = pol (x - (k[[m]] + k[[n]])), {m, n + 1, deg}], {n, 1, >>>> deg - 1}]; kk = Round[CoefficientList[Expand[pol], x]]; pl = 0; Do[ >>>> pl = pl + kk[[n]] x^(n - 1), {n, 1, Length[kk]}]; Factor[pl] >>>> >>>> (*3 roots resolventi*) >>>> p = x^9 - 9*x^7 - 21*x^6 + 72*x^5 + 99*x^4 - 99*x^3 - 585*x^2 + >>>> 549*x + 166; deg = Length[CoefficientList[p, x]] - 1; k = >>>> N[x /. Solve[p == 0, x], 300]; pol = 1; aa = {}; Do[ >>>> Do[Do[pol = pol (x - (k[[m]] + k[[n]] + k[[r]])), {r, m + 1, >>>> deg}], {m, n + 1, deg - 1}], {n, 1, deg - 2}]; kk = >>>> Round[CoefficientList[Expand[pol], x]]; pl = 0; Do[ >>>> pl = pl + kk[[n]] x^(n - 1), {n, 1, Length[kk]}]; Factor[pl] >>>> >>>> and one question to Daniel or another Wolfram staff programist, how >>>> is upper limitation for order of polynomial to function Factor or >>>> FactorList >>>> on integer coefficients polynomial in Mathematica (or is unlimited >>>> with only one limitation memory of computer)? >>>> >>>> My sample code of new function of Mathematica Galois is following >>>> (this sample counts Galois groups of polynomials up to order 5) >>>> I don't expand them on higher orders yet because in case more than >>>> one variable polynomials this algorhitm will be useless because not >>>> possibility count disriminant in such case and is necessary rebuild >>>> them on algorhitm which will be don't use information if >>>> Sqrt[Discriminant]] is perfect square or not. >>>> >>>> I will be greatfull for help and rebuild my statical procedure on >>>> dynamical which is absolutely necessary start from order 8 which is >>>> smallest order where algorithm of Lenstra-Czebotaryev (my previosly >>>> ask about help in message How to do quickest) is not sufficient to >>>> exact determination of Galois groups. >>>> >>>> Follow finish algorhitm determination of Galois grups will be >>>> available to write next to solving polynomials higher degree as 4 >>>> for radicals. >>>> >>>> Best wishes >>>> Artur >>>> >>>> W dniu 2011-02-20 03:59, DrMajorBob pisze: >>>>> Factorizations modulo prime p are NOT square-free, hence the >>>>> factor degrees modulo p do not partition the polynomial order, in >>>>> the following cases: >>>>> >>>>> poly = x^8 - x - 1; >>>>> order = Exponent[poly, x]; >>>>> limit = Prime[order!]; >>>>> Cases[FactorInteger[ >>>>> Discriminant[poly, x]], {p_, _} :> {p, FactorList[poly, Modulus >>>>> -> p]} /; >>>>> 0< p< limit] >>>>> >>>>> {{11, {{1, 1}, {9 + x, 2}, {8 + 5 x + 3 x^2 + 10 x^3 + x^4 + 4 x^5 >>>>> + x^6, >>>>> 1}}}} >>>>> >>>>> That is, (9 + x)^2 is a factor mod 11. >>>>> >>>>> That being taken care of, we don't need to print exceptions. I >>>>> also removed Rest, Sort, Reverse and an instance of Set from the >>>>> loop. >>>>> >>>>> All that made surprisingly little difference: >>>>> >>>>> Timing[ >>>>> poly = x^8 - x - 1; >>>>> order = Exponent[poly, x]; >>>>> partitions = >>>>> Replace[IntegerPartitions@order, >>>>> x_List :> Flatten@{0, Reverse@x}, {1}]; >>>>> np = Length@partitions; >>>>> Clear[counts, index]; >>>>> counts[_] = 0; >>>>> index[_] = 0; >>>>> Do[index@partitions[[j]] = j, {j, np}]; >>>>> n = 0; >>>>> While[n< order!, >>>>> p = Prime[++n]; >>>>> ndx = index@Exponent[FactorList[poly, Modulus -> p][[All, 1]], x]; >>>>> counts[ndx]++ >>>>> ]; >>>>> Array[counts, np]] >>>>> >>>>> {9.86778, {4996, 5781, 3361, 3449, 2653, 4055, 1360, 1249, 3360, >>>>> 1321, >>>>> 2470, 412, 1103, 1114, 1651, 1129, 105, 102, 416, 206, 25, 1}} >>>>> >>>>> Bobby >>>>> >>>>> On Sat, 19 Feb 2011 19:02:18 -0600, Daniel Lichtblau >>>>> <danl at wolfram.com> wrote: >>>>> >>>>>> The idea, which got removed from my code, is that the factor >>>>>> pattern will be a partition of the degree EXCEPT when the prime >>>>>> in question divides the discriminant of the polynomial. Reason >>>>>> being, those are the cases where the factorization is not square >>>>>> free. >>>>>> >>>>>> I'm not sure exactly how Artur's rendition of my code even works >>>>>> without spewing a couple of messages. My hash table ("htab") was >>>>>> pattern-free, that is to say, it never should return 0. Bobby's >>>>>> code fixed that issue. I had worked around it simpl by keeping >>>>>> Artur's original discriminant divisor check. >>>>>> >>>>>> Daniel >>>>>> >>>>>> >>>>>> ----- Original Message ----- >>>>>>> From: "DrMajorBob"<btreat1 at austin.rr.com> >>>>>>> To: drmajorbob at yahoo.com, "Artur"<grafix at csl.pl> >>>>>>> Cc: mathgroup at smc.vnet.net, "Sjoerd C. de Vries" >>>>>>> <sjoerd.c.devries at gmail.com>, "Daniel Lichtblau"<danl at wolfram.com> >>>>>>> Sent: Saturday, February 19, 2011 4:20:10 PM >>>>>>> Subject: Re: [mg116541] Re: How to do quickest >>>>>>> Here's a more readable code (IMHO). I like variable names that mean >>>>>>> something. >>>>>>> >>>>>>> Timing[ >>>>>>> poly = x^8 - x - 1; >>>>>>> order = Exponent[poly, x]; >>>>>>> partitions = IntegerPartitions@order; >>>>>>> np = Length@partitions; >>>>>>> index[_] = 0; >>>>>>> Do[index[partitions[[j]]] = j, {j, np}]; >>>>>>> aa = ConstantArray[0, np]; >>>>>>> n = total = 0; >>>>>>> While[total< order!, >>>>>>> p = Prime[++n]; >>>>>>> factors = FactorList[poly, Modulus -> p][[All, 1]]; >>>>>>> ndx = index@Reverse@Sort@Rest@Exponent[factors, x]; >>>>>>> Positive@ndx&& (total++; aa[[ndx]]++) >>>>>>> ]; >>>>>>> aa] >>>>>>> >>>>>>> {9.84093, {4996, 5781, 3361, 3449, 2653, 4055, 1360, 1249, 3360, >>>>>>> 1321, >>>>>>> 2470, 412, 1103, 1114, 1652, 1129, 105, 102, 416, 206, 25, 1}} >>>>>>> >>>>>>> The following yields the tidbit of information that ndx == 0 >>>>>>> EXACTLY >>>>>>> ONCE: >>>>>>> >>>>>>> Timing[ >>>>>>> poly = x^8 - x - 1; >>>>>>> order = Exponent[poly, x]; >>>>>>> partitions = IntegerPartitions@order; >>>>>>> np = Length@partitions; >>>>>>> Clear[counts, index]; >>>>>>> counts[_] = 0; >>>>>>> index[_] = 0; >>>>>>> Do[index@partitions[[j]] = j, {j, np}]; >>>>>>> n = total = 0; >>>>>>> While[total< order!, >>>>>>> p = Prime[++n]; >>>>>>> factors = FactorList[poly, Modulus -> p][[All, 1]]; >>>>>>> ndx = index@Reverse@Sort@Rest@Exponent[factors, x]; >>>>>>> counts[ndx]++; >>>>>>> Positive@ndx&& total++ >>>>>>> ]; >>>>>>> Table[counts@i, {i, 0, np}]] >>>>>>> >>>>>>> {10.2563, {1, 4996, 5781, 3361, 3449, 2653, 4055, 1360, 1249, 3360, >>>>>>> 1321, 2470, 412, 1103, 1114, 1652, 1129, 105, 102, 416, 206, 25, >>>>>>> 1}} >>>>>>> >>>>>>> If we can COUNT on that being so for some reason, the code becomes >>>>>>> >>>>>>> Timing[ >>>>>>> poly = x^8 - x - 1; >>>>>>> order = Exponent[poly, x]; >>>>>>> partitions = IntegerPartitions@order; >>>>>>> np = Length@partitions; >>>>>>> Clear[counts, index]; >>>>>>> counts[_] = 0; >>>>>>> index[_] = 0; >>>>>>> Do[index@partitions[[j]] = j, {j, np}]; >>>>>>> n = 0; >>>>>>> While[n< order!, >>>>>>> p = Prime[++n]; >>>>>>> factors = FactorList[poly, Modulus -> p][[All, 1]]; >>>>>>> ndx = index@Reverse@Sort@Rest@Exponent[factors, x]; >>>>>>> counts[ndx]++; >>>>>>> ndx == 0&& >>>>>>> Print@{n, p, factors, Reverse@Sort@Rest@Exponent[factors, x]}; >>>>>>> ]; >>>>>>> Array[counts, np]] >>>>>>> >>>>>>> {5,11,{1,9+x,8+5 x+3 x^2+10 x^3+x^4+4 x^5+x^6},{6,1}} >>>>>>> >>>>>>> {9.95797, {4996, 5781, 3361, 3449, 2653, 4055, 1360, 1249, 3360, >>>>>>> 1321, >>>>>>> 2470, 412, 1103, 1114, 1651, 1129, 105, 102, 416, 206, 25, 1}} >>>>>>> >>>>>>> where I have shown the exception. >>>>>>> >>>>>>> The problem is a black box to me, but finding just ONE exception is >>>>>>> not >>>>>>> what I would have expected. >>>>>>> >>>>>>> Bobby >>>>>>> >>>>>>> On Sat, 19 Feb 2011 12:02:22 -0600, Artur<grafix at csl.pl> wrote: >>>>>>> >>>>>>> > Thank You for procedure! On my computer >>>>>>> > {12.312, {4996, 5781, 3361, 3449, 2653, 4055, 1360, 1249, 3360, >>>>>>> > 1321, >>>>>>> > 2470, 412, 1103, 1114, 1652, 1129, 105, 102, 416, 206, 25,1}} >>>>>>> > >>>>>>> > Bob try Your quickest steps combined with following which is >>>>>>> still >>>>>>> > little quickest >>>>>>> > >>>>>>> > {10.313, {4996, 5781, 3361, 3449, 2653, 4055, 1360, 1249, 3360, >>>>>>> > 1321, >>>>>>> > 2470, 412, 1103, 1114, 1652, 1129, 105, 102, 416, 206, 25, 1}} >>>>>>> > >>>>>>> > (*Daniel Lichtblau modified by Artur Jasinski*) >>>>>>> > >>>>>>> > Timing[cc = {}; pol = x^8 - x - 1; >>>>>>> > nn = Length[CoefficientList[pol, x]] - 1; >>>>>>> > pp = IntegerPartitions[nn]; >>>>>>> > Do[htab[pp[[j]]] = j, {j, Length[pp]}]; >>>>>>> > aa = Table[0, {Length[pp]}]; >>>>>>> > n = 1; cn = 0; >>>>>>> > While[cn< nn!, p = Prime[n]; >>>>>>> > n++; >>>>>>> > kk = FactorList[pol, Modulus -> p]; >>>>>>> > ww = Rest[Exponent[kk[[All, 1]], x]]; >>>>>>> > ww = Reverse[Sort[ww]]; >>>>>>> > pos = htab[ww]; >>>>>>> > If[pos == 0, , cn++; aa[[pos]] = aa[[pos]] + 1]]; >>>>>>> > aa] >>>>>>> > >>>>>>> > {10.313, {4996, 5781, 3361, 3449, 2653, 4055, 1360, 1249, 3360, >>>>>>> > 1321, >>>>>>> > 2470, 412, 1103, 1114, 1652, 1129, 105, 102, 416, 206, 25, 1}} >>>>>>> > >>>>>>> > W dniu 2011-02-19 17:37, DrMajorBob pisze: >>>>>>> >> This is easier to read, if no faster. >>>>>>> >> >>>>>>> >> Timing[ >>>>>>> >> Clear[a, c]; >>>>>>> >> pol = x^8 - x - 1; >>>>>>> >> nn = Length@CoefficientList[pol, x] - 1; >>>>>>> >> If[ >>>>>>> >> IrreduciblePolynomialQ[pol], >>>>>>> >> a[i_] = {}; >>>>>>> >> c[i_] := Length@Flatten[a@i]; >>>>>>> >> pp = IntegerPartitions@nn; >>>>>>> >> b = FactorInteger[Discriminant[pol, x]][[All, 1]]; >>>>>>> >> n = 1; >>>>>>> >> cn = 0; >>>>>>> >> While[cn< nn!, p = Prime@n; >>>>>>> >> If[! MemberQ[b, p], >>>>>>> >> cn++; >>>>>>> >> k = Reverse@Rest@FactorList[pol, Modulus -> p][[All, 1]]; >>>>>>> >> w = Length@CoefficientList[#, x] - 1& /@ k; >>>>>>> >> pos = Position[pp, w, 1, 1][[1, 1]]; >>>>>>> >> a[pos] = {a[pos], p}]; >>>>>>> >> n++]]; >>>>>>> >> Array[c, Length@pp] >>>>>>> >> ] >>>>>>> >> >>>>>>> >> {10.8518, {4996, 5781, 3361, 3449, 2653, 4055, 1360, 1249, 3360, >>>>>>> >> 1321, >>>>>>> >> 2470, 412, 1103, 1114, 1652, 1129, 105, 102, 416, 206, 25, >>>>>>> 1}} >>>>>>> >> >>>>>>> >> Bobby >>>>>>> >> >>>>>>> >> On Sat, 19 Feb 2011 04:12:14 -0600, Sjoerd C. de Vries >>>>>>> >> <sjoerd.c.devries at gmail.com> wrote: >>>>>>> >> >>>>>>> >>> Without changing the basic operation of your algorithm I've >>>>>>> >>> changed a >>>>>>> >>> couple of details. The difference is not huge, but about 20% of >>>>>>> >>> speed >>>>>>> >>> gain is still nice. >>>>>>> >>> >>>>>>> >>> pol = x^8 - x - 1; >>>>>>> >>> nn = Length[CoefficientList[pol, x]] - 1; >>>>>>> >>> If[IrreduciblePolynomialQ[pol], >>>>>>> >>> pp = IntegerPartitions[nn]; >>>>>>> >>> aa = Table[{}, {n, 1, Length[pp]}]; Print[aa]; >>>>>>> >>> ff = FactorInteger[Discriminant[pol, x]]; >>>>>>> >>> bb = Table[ff[[n, 1]], {n, 1, Length[ff]}]; >>>>>>> >>> n = 1; >>>>>>> >>> cn = 0; >>>>>>> >>> While[cn< nn!, >>>>>>> >>> p = Prime[n]; >>>>>>> >>> If[MemberQ[bb, p], >>>>>>> >>> (*True*), >>>>>>> >>> cn++; >>>>>>> >>> kk = FactorList[pol, Modulus -> p]; >>>>>>> >>> ww = Table[ >>>>>>> >>> Length[CoefficientList[kk[[m, 1]], x]] - 1, >>>>>>> >>> {m, Length[kk], 2, -1} >>>>>>> >>> ]; >>>>>>> >>> pos = Position[pp, ww, 1, 1][[1, 1]]; >>>>>>> >>> aa[[pos]] = {aa[[pos]], p}; >>>>>>> >>> ]; >>>>>>> >>> n++ >>>>>>> >>> ] >>>>>>> >>> ]; aa = Map[Flatten, aa, {1}]; >>>>>>> >>> Table[Length[aa[[m]]], {m, 1, Length[aa]}] >>>>>>> >>> ] >>>>>>> >>> >>>>>>> >>> >>>>>>> >>> Cheers -- Sjoerd >>>>>>> >>> >>>>>>> >>> >>>>>>> >>> On Feb 15, 12:33 pm, Artur<gra... at csl.pl> wrote: >>>>>>> >>>> Dear Mathematica Gurus, >>>>>>> >>>> How to do following procedure quickest? >>>>>>> >>>> (*start*) >>>>>>> >>>> pol = x^8 - x - 1; nn = Length[CoefficientList[pol, x]] - >>>>>>> 1; If[ >>>>>>> >>>> IrreduciblePolynomialQ[pol], pp = IntegerPartitions[nn]; aa = >>>>>>> >>>> {}; >>>>>>> >>>> Do[AppendTo[aa, {}], {n, 1, Length[pp]}]; Print[aa]; >>>>>>> >>>> ff = FactorInteger[Discriminant[pol, x]]; bb = {}; >>>>>>> >>>> Do[AppendTo[bb, ff[[n]][[1]]], {n, 1, Length[ff]}]; n = 1; >>>>>>> cn = >>>>>>> >>>> 0; >>>>>>> >>>> While[cn< nn!, p = Prime[n]; >>>>>>> >>>> If[MemberQ[bb, p], , cn = cn + 1; >>>>>>> >>>> kk = FactorList[pol, Modulus -> p]; ww = {}; >>>>>>> >>>> Do[cc = Length[CoefficientList[kk[[m]][[1]], x]]; >>>>>>> >>>> AppendTo[ww, cc - 1], {m, 2, Length[kk]}]; ww = >>>>>>> Reverse[ww]; >>>>>>> >>>> pos = Position[pp, ww][[1]][[1]]; AppendTo[aa[[pos]], >>>>>>> >>>> Prime[n]]]= >>>>>>> >>> ; >>>>>>> >>>> n++]]; Table[Length[aa[[m]]], {m, 1, Length[aa]}] >>>>>>> >>>> (*end*) >>>>>>> >>>> Best wishes >>>>>>> >>>> Artur >>>>>>> >>> >>>>>>> >>> >>>>>>> >> >>>>>>> >> >>>>>>> >>>>>>> >>>>>>> -- DrMajorBob at yahoo.com >>>>> >>>>> >>> >>> > >

**Re: How to do quickest**

**Re: Factor/Simplify complex expressions**

**Re: How to do quickest**

**(MathLink) Correct handling of Messages generated by slave kernel**