# How to solve an equation system whose equations are based on Mod[x,y]

I'm trying to find the parameters of a specific linear congruential generator (pseudo random number generator) based on a list of the first five numbers generated by it. The list is this one:

{0.786953, 0.689557, 0.432612, 0.303753, 0.730747}


I tried to use NSolve for this, but can't seem to solve it (I wrote the complete numbers estimated by Mathematica in here):

NSolve[{N[Mod[z1*a1 + c1, m1]]/(m1 - 1) == 0.7869529629996759,
0.6895574464879732 ==
N[Mod[a1*0.7869529629996759*(m1 - 1) + c1, m1]]/(m1 - 1),
4326116672868895 ==
N[Mod[a1*0.6895574464879732*(m1 - 1) + c1, m1]]/(m1 - 1),
0.30375286671507773 ==
N[Mod[a1*4326116672868895*(m1 - 1) + c1, m1]]/(m1 - 1)
}, {z1, a1, c1, m1}]


Is there any other way to solve an equation system whose equations are based on Mathematica functions, such as these ones with Mod[x,y]? I'm pretty new to Mathematica, so I'm sorry if my code is too messy :c

P.S: the actual parameters are a=314 159 269, c=453 806 245, m=2^31 and the first seed z1=12345678, in case an approximation is needed.

• This generates your example numbers pretty quickly: BlockRandom[SeedRandom[12345678, Method -> {"Congruential", "Multiplier" -> 314159269, Increment -> 453806245, Modulus -> 2^31}]; RandomReal[1, 4, WorkingPrecision -> 20]] – J. M.'s discontentment Dec 18 '15 at 3:43
• That's an interesting function! But I wasn't asking how to generate those numbers, I need to know how to figure out the parameters of that generator based on those first generated numbers. – Sizigia Dec 18 '15 at 7:56
• I wrote the actual parameters because I was trying to achieve this with a generator I made, and so that they could be used if an approximation to the solution was needed. But I need to be able to estimate the unknown parameters of any other LCG. – Sizigia Dec 18 '15 at 8:01
• That's why it's a comment and not an answer. Anyway, have you seen this? – J. M.'s discontentment Dec 18 '15 at 8:05
• Nevertheless, that this can even be contemplated is a reason why LCGs are worse than useless in cryptographic applications. – J. M.'s discontentment Dec 18 '15 at 10:40

This almost scares me, but FindInstance is able to solve this under the assumption that we know the modulus m1. I had to make a correction or two along the way.

m1 = 2^31;
vals0 = {0.7869529629996759, 0.6895574464879732, .4326116672868895,
0.30375286671507773};
vals = Join[{v1}, Round[(m1 - 1)*vals0]]

(* Out[327]= {v1, 1689968619, 1480813340, 929026481, 652304314} *)


We want an integer system of equations. Notice that I have already converted the given values to integers, under the assumption that the decimals are approximated by corresponding integers divided by m1-1. I finesse the Mod part by adding an (unknown) integer multiple of m1 to get equalities. These unknowns simply get added to the variable list.

polys =
Table[vals[[j]]*a1 + c1 + k[j]*m1 - vals[[j+1]], {j,Length[vals]-1}]

(* Out[328]= {-1689968619 + c1 + a1 v1 + 2147483648 k[1], -1480813340 +
1689968619 a1 + c1 + 2147483648 k[2], -929026481 + 1480813340 a1 +
c1 + 2147483648 k[3], -652304314 + 929026481 a1 + c1 +
2147483648 k[4]} *)


Now we find a solution.

FindInstance[polys == 0,
Flatten[{a1, c1, v1, Array[k, 4]}], Integers]

(* Out[331]= {{a1 -> 314159269, c1 -> 453806245, v1 -> 12345678,
k[1] -> -1806071, k[2] -> -247228567, k[3] -> -216630863,
k[4] -> -135908965}} *)


--- edit ---

Okay, I realized later how to solve first for m1 and the integers from which the decimal values are constructed (by dividing the integers by m1-1 and using the nearest machine double). I should note that this might be somewhat uncommon insofar as it will turn out we divide not by a power of 2, but rather one less. Dividing by an exact power of 2 means we can fill in the mantissa bits of our random double directly from the integer bits of the numbers we generate. But I digress.

The method I will show is taken from this recent response to a different MSE post. The idea is to use simultaneous rational approximation of the given decimal values, multiplied by some "large" exponent. How large can take fiddling but if we suspect the decimal values came from integers divided by whatever, then since they fit in machine doubles I go with something modestly larger than the mantissa scale (which is 2^53). We form an appropriate lattice, reduce it, then find a usable row to recover our desired values.

vals0 = {0.7869529629996759, 0.6895574464879732, .4326116672868895,
0.30375286671507773};
deg = 60;
mult = 2^deg;
vals1 = Round[mult*vals0/Min[vals0]];
ivec = Prepend[vals1, 1];
lat = Join[{ivec}, Rest[mult*IdentityMatrix[Length[ivec]]]];
redlat = LatticeReduce[lat];
rows = Select[redlat, #[[1]] =!= 0 &];
rnum = Ordering[Map[N[#].# &, rows], 1][[1]];
xvals = Abs[(rows[[rnum, 1]]*Rest[ivec] - Rest[rows[[rnum]]])/mult]

(* Out[215]= {1689968619, 1480813340, 929026481, 652304314} *)


And we have recovered the integers that gave rise to the decimal values. Now it is trivial to recover m1. I take an average below but actually every entry in the list gave the same (correct) value.

m1 = Round[Mean[xvals/vals0]] + 1

(* Out[184]= 2147483648 *)


I also realized how FindInstance might be handling this. All but the first equation is linear, so possibly it is discarding that one initially, solving, and then using back substitution into the first equation to finish the job. This is just speculation though.

Also note that our handling of the modulus inflated the number of variables so that there were more than there were equations. That, it turns out, is not so big a deal. There is an example with some explanation in section 4, "Computation and use of matrix Hermite normal forms", of the paper found at this link (and please pardon the self reference).

--- end edit ---

• Wow, that was impressive, thanks! So I assume m1 needs to be known beforehand because if not, there would be 'unevaluable' Round[val*m1] expressions in polys? (I tried your code without converting the values to integers and yes, FindInstance just returns the trivial zero solution). – Sizigia Dec 17 '15 at 20:31
• I was not able to get this to work without knowing the modulus. But then again, I was surprised something worked even knowing the modulus. – Daniel Lichtblau Dec 17 '15 at 21:05
• As for guessing the modulus, this previous answer of yours seems to be useful. – J. M.'s discontentment Dec 18 '15 at 15:51
• @J.M. Umm, this is eerie timing. I was writing up exactly that method when you wrote this comment. I confess I am also happy to see someone else made that connection. – Daniel Lichtblau Dec 18 '15 at 16:11
• I was staring at your answer for a while, and it suddenly clicked that you had written an answer before on the old AffineRationalize[] function. :) In any event, did you see the paper I linked to in a comment in the OP? – J. M.'s discontentment Dec 18 '15 at 16:14