MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Implementation rest modulo 2 needed


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


  • Prev by Date: Re: Solve this differential equation with periodic boundary conditions:
  • Next by Date: Mathematica 7: Histogram Y-axis range?
  • Previous by thread: Re: Implementation rest modulo 2 needed
  • Next by thread: Re: Re: Implementation rest modulo 2 needed