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 ---
BlockRandom[SeedRandom[12345678, Method -> {"Congruential", "Multiplier" -> 314159269, Increment -> 453806245, Modulus -> 2^31}]; RandomReal[1, 4, WorkingPrecision -> 20]]$\endgroup$