Re: How to do quickest
- To: mathgroup at smc.vnet.net
- Subject: [mg116566] Re: How to do quickest
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Sun, 20 Feb 2011 05:24:53 -0500 (EST)
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