Re: Implementation rest modulo 2 needed
- To: mathgroup at smc.vnet.net
- Subject: [mg93859] Re: [mg93789] Implementation rest modulo 2 needed
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Thu, 27 Nov 2008 05:31:13 -0500 (EST)
- References: <200811251216.HAA03013@smc.vnet.net> <200811261010.FAA19492@smc.vnet.net>
Artur wrote: > Dear Mathematica Gurus, > > I want to construct algorhitm finding two divisors of numbers (without > FactorIneteger and with uses only Reduce, GroebnerBasis etc.) > > e.g. for divisors of number 15 in binary system these will be > > Reduce[{x a == 1 && y a + x b == 1 && y b + x c == 1 && y c == 1 && > x (x - 1) == 0 && y (y - 1) == 0 && c (c - 1) == 0 && > b (b - 1) == 0 && a (a - 1) == 0}, {x, y, a, b, c}, Integers] > > but to find divisors of 15=1111 wasn't necessary to use rest modulo 2 > > How do that for e.g. 735 = 1011011111 where we have to uses rest modulo 2 > > Any idea ? > > Best wishes > Artur > When you go this route the factoring becomes, computationally, a very difficult problem in 0-1 quadratic programming (QP). That said, here we go. We'll assume we know the number of bits in each factor (one can try all possibilities, after all). We'll also assume we've removed factors of 2, so the trailing bits of the factors are both 1. The short answer to your question is to use one equality constraint, that in effect multiplies two factors as polynomials in powers of 2, and subtracts the desired number. (I seem to recall there are other ways that separately keep track of each bit in the product under construction, but I think the degree of the constraints then becomes larger than 2. Which is bad.) Here are a couple of ways to set up the problem. One is direct and uses Groebnerbasis behind the scenes. iFactorQP[val_, n1_, n2_] := Module[{x1, x2, v1, v2, vars, poly1, poly2, poly, polys, rule, soln}, v1 = Array[x1, n1]; v2 = Array[x2, n2]; rule = {x1[1] -> 1, x1[n1] -> 1, x2[1] -> 1, x2[n2] -> 1}; poly1 = (v1. 2^Range[0, n1 - 1]) /. rule; poly2 = (v2. 2^Range[0, n2 - 1]) /. rule; poly = Expand[poly1*poly2 - val]; vars = Variables[poly]; polys = Join[Map[(#^2 - #) &, vars], {poly}]; soln = NSolve[polys, vars]; Round[{poly1, poly2} /. soln /. rule] ] In[23]:= Timing[iFactorQP[735, 6, 5]] Out[23]= {0.107983, {{35, 21}}} It hangs attempting ifactor[Prime[55]*Prime[79], 9, 9]. Another approach is to use what is called a lineasr reformulation of the QP. Solve the relaxed LPs, then branch on noninteger values in the result. Order so we treat high bits as more important than low bits (so it has a chance of completing before your disk dies). Here is some code to do this. You can change pPrint to Print in places, to see what it is doing while it runs. ifactor[val_, n1_, n2_] := Module[{x1, x2, x, v1, v2, prodvars, allvars, vars, c1, c2, c3, c4, c5, rule1, rule2, rule3, constraints, obj, program, min, stack, tmp, soln, fax = {}, isol, badpos, p2, posn, iter = 0, lmax, eps = 10^(-6), vv1, rvars}, v1 = Array[x1, n1]; v2 = Array[x2, n2]; vv1 = PadRight[v1, n2]; vv1 = Flatten[Transpose[{vv1, v2}]] /. 0 -> Sequence[]; prodvars = Flatten[Array[x, {n1, n2}]]; allvars = Flatten[{vv1, prodvars}]; pPrint[allvars]; c1 = Map[0 <= # <= 1 &, allvars]; c2 = Map[# <= x1[#[[1]]] &, prodvars]; c3 = Map[# <= x2[#[[2]]] &, prodvars]; c4 = Map[# >= x1[#[[1]]] + x2[#[[2]]] - 1 &, prodvars]; rule1 = {x[1, 1] -> 1, x[1, n2] -> 1, x[n1, 1] -> 1, x[n1, n2] -> 1}; rule2 = {x[1, a_] :> x2[a], x[n1, a_] :> x2[a], x[a_, 1] :> x1[a], x[a_, n2] :> x1[a]}; rule3 = {x1[1] -> 1, x1[n1] -> 1, x2[1] -> 1, x2[n2] -> 1}; constraints = Union[Flatten[{c1, c2, c3, c4}] /. rule1 /. rule2 /. rule3 /. {1 >= _ :> Sequence[], a_ <= a_ :> Sequence[], a_ >= a_ :> Sequence[], True :> Sequence[]}]; obj = Sum[2^(n1 + n2 - j - k)*x[j, k], {j, n1}, {k, n2}] /. rule1 /. rule2; constraints = Prepend[constraints, obj == val]; allvars = Union[allvars /. rule1 /. rule2 /. rule3 /. (1 :> Sequence[])]; vars = Union[Join[v1, v2] /. rule3 /. (1 :> Sequence[])]; vars = Sort[vars, OrderedQ[{#1[[1]], #2[[1]]}] &]; rvars = vars; pPrint["start loop"]; stack = {constraints, {}}; While[stack =!= {}, program = First[stack]; stack = Last[stack]; Internal`DeactivateMessages[ soln = NMinimize[{1, program}, allvars], NMinimize::"nsol"]; iter++; If[! FreeQ[soln, Indeterminate] || Head[soln] == NMinimize, Continue[]]; {tmp, soln} = soln; isol = rvars /. soln; badpos = Position[ isol, (a_ /; Abs[a - Round[a]] > eps && Head[a] =!= x), {1}, 1, Heads -> False]; If[badpos == {}, fax = soln; Break[]]; badpos = badpos[[1, 1]]; pPrint[rvars[[badpos]]]; tmp = Head[rvars[[badpos]]]; posn = First[rvars[[badpos]]]; If[tmp === x1, c1 = Cases[allvars, x[posn, _]], c1 = Cases[allvars, x[_, posn]]]; p2 = program /. {rvars[[badpos]] -> 0} /. Thread[c1 -> 0]; stack = {Union[ Flatten[{p2, rvars[[badpos]] == 0, Thread[c1 == 0]}]] /. True -> Sequence[], stack}; p2 = program /. {rvars[[badpos]] -> 1}; If[tmp === x1, posn = 2; tmp = x2, posn = 1; tmp = x1]; stack = {Union[ Flatten[{p2, rvars[[badpos]] == 1, Map[# == tmp[#[[posn]]] &, c1]}]] /. True -> Sequence[], stack};]; Round[{1 + 2^(n1 - 1) + Sum[2^(n1 - j)*x1[j], {j, 2, n1 - 1}], 1 + 2^(n2 - 1) + Sum[2^(n2 - j)*x2[j], {j, 2, n2 - 1}]} /. fax]] The following settings also seem to help. SetOptions[LinearProgramming, Method -> "InteriorPoint"]; SetSystemOptions[ "LinearProgrammingOptions" -> {"InteriorPointSize" -> 1, "Preprocessing" -> True}]; It's faster than the NSolve method, but still slow. In[37]:= Timing[ss = ifactor[Prime[55]*Prime[79], 9, 9]] Out[37]= {0.231965, {401, 257}} In[36]:= Timing[ss = ifactor[Prime[555]*Prime[779], 13, 12]] Out[36]= {151.839, {5927, 4019}} Let me address an unasked question: No, I did not just write this up. It's code I played with a couple of years ago. Believe me, this approach is not worth much. Daniel Lichtblau Wolfram Research
- References:
- Using Alt or Cmd Keys to Enter Cell Types
- From: "David Park" <djmpark@comcast.net>
- Implementation rest modulo 2 needed
- From: Artur <grafix@csl.pl>
- Using Alt or Cmd Keys to Enter Cell Types