Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 2008

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

Search the Archive

Re: Implementation rest modulo 2 needed

Dear Daniel,
One question more. If we will be have strictly matematical formula on 
adding numbers in binary system (not algorhitm but explicite formula) 
what will be changed in our difficult situation ?
Best wishes
Artur

Daniel Lichtblau pisze:
> 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
>
> __________ Information from ESET NOD32 Antivirus, version of virus 
> signature database 3643 (20081126) __________
>
> The message was checked by ESET NOD32 Antivirus.
>
> http://www.eset.com
>
>
>
  • Prev by Date: Re: Compiled function
  • Next by Date: Version 7 Cursor for Locators
  • Previous by thread: Implementation rest modulo 2 needed
  • Next by thread: Re: Implementation rest modulo 2 needed