# On Solutions of Diophantine Equation

By Mathematica, I try to obtain solutions of Diophantine equations such as:

$$F_{n_1}F_{n_2}\cdots F_{n_k}+1=F_t$$

where the sequence $\{F_n\}$ is the Fibonacci sequence,$n_1 < n_2 < \cdots < n_k$ are positive integers, and $1 < t < 21$.

For example, one of the solutions is $F_2\times F_3+1=F_4$.

How can I find the solutions by Mathematica?

• Do you want to count all combinations with f[1] and f[2] separately? That is, for any solution with 2<n1 there are three more solutions that also use f[1], f[2] or f[1]*f[2]. Mar 4 '16 at 13:13
• Yes, I want to obtain all solution of the diophantine equation including f[1] and f[2] separately Mar 4 '16 at 13:45

n       = 20;
ff      = Fibonacci /@ Range@n;
sus     = Subsets[Fibonacci /@ Range@n, {1, n}];
susf    = Subsets[f /@ Range@n, {1, n}];
results = 1 + Times @@@ sus;
e = Extract[susf, p = Position[results, x_?(MemberQ[ff, #] &)]]


These are the valid products:

(*{
{f[1]},
{f[2]},
{f[3]},
{f[1], f[2]},
{f[1], f[3]},
{f[2], f[3]},
{f[1], f[2], f[3]}}
*)


To get the full equations:

res = f /@ Flatten[Position[ff, #] & /@ Extract[results, p]];
Equal @@@ Transpose[{1 + Times @@@ e, res}] // Column

(*
1 + f[1] == f[3]
1 + f[2] == f[3]
1 + f[3] == f[4]
1 + f[1] f[2] == f[3]
1 + f[1] f[3] == f[4]
1 + f[2] f[3] == f[4]
1 + f[1] f[2] f[3] == f[4]
*)

• Ah nice, I like the simplicity, but this is probably not feasible if you want to go much beyond the first 20 Fibonacci numbers. Mar 4 '16 at 15:51
• @MartinBüttner The number of SetPartitions grows pretty fast too f = (Join @@ ConstantArray @@@ FactorInteger[Fibonacci@# - 1]) & /@ Range@100; Max[BellB[Length /@ f]] Mar 4 '16 at 16:50
• Sure it's still a combinatorics problem. But the set partitions grow a lot slower than $2^n$ (mainly because the number of factors doesn't grow very fast). My solution solves up to $n = 50$ in a few seconds (and is dominated by a few numbers with many prime factors), whereas this approach will stop being feasible around $n = 30$. Mar 4 '16 at 17:19
• @MartinBüttner Yup. My try was to find a simple way to solve the problem with n = 20, as stated. I don't pretend this approach to be scalable. Mar 4 '16 at 17:33

Here is a subset sum approach using lattice reduction. For this we work with logs, scaled up and rounded to get an integer lattice. It's not guaranteed to find a solution but in practice it seems quite good at finding near misses wherein a quotient of such products "almost" works.

fibProd[n_Integer] := Catch[Module[
{fibs = Fibonacci[Range[3, n - 1]], fiblogs,
fm1log = Log[Fibonacci[n] - 1], rlogs, lat, redlat, row},
fiblogs = Log[fibs];
rlogs = Round[2^20*Append[fiblogs, fm1log]];
lat = Join[Transpose[{rlogs}], IdentityMatrix[n - 2], 2];
lat[[-1, -1]] = 2^10;
lat = LatticeReduce[lat];
If[Length[Select[lat, Abs[Last[#]] == 2^10 &]] > 1,
Print[{"found more", n}]];
row = SelectFirst[lat, Abs[Last[#]] == 2^10 &];
If[! FreeQ[row, Missing], Throw[$Failed]]; row = row[[2 ;; -2]]; If[Max[Abs[row]] > 1, Throw[$Failed]];
1/fibs^row
]]


Running it to 220 produced only the trivial result fib(2)+1=fib(3).

This answer currently ignores solutions involving f[1] and f[2], because they are trivial to find by hand:

f[1] + 1 = f[3]
f[2] + 1 = f[3]
f[1]*f[2] + 1 = f[3]


Additionally, by prepending these three factors to any other solution you can get three more solutions. E.g. if you have f[3] + 1 = f[4], then 3 more solutions are:

f[1]*f[3] + 1 = f[4]
f[2]*f[3] + 1 = f[4]
f[1]*f[2]*f[3] + 1 = f[4]


In fact, it appears that these 6 are the only solutions up to $F_{50}$ (unless I made a mistake, of course).

I first tried putting a system of equations into FindInstance, but it's incredibly slow. Instead, it's much easier here to look for possible solutions manually, by checking in how many ways we can express starting from $F_t-1$ as a product of distinct Fibonacci numbers:

Needs["Combinatorica"];
maxN = 21;
(* Lookup table of admissible Fibonacci numbers. *)
fibs = Fibonacci@Range@maxN
solutions = {};
For[n = 4, n < maxN, ++n,
(* Get the list of (F_n-1)'s prime factors. *)
factors = Join @@ ConstantArray @@@ FactorInteger[Fibonacci@n - 1];
(* Now we get all, not necessarily contiguous partitions of the factors from
which we reconstruct all possible decompositions into factors of F_n-1. *)
decompositions = Union[Sort /@ (Times @@@ # & /@ SetPartitions[factors])];
(* We retain only those decompositions that contain only distinct numbers, all
of which have to be Fibonacci numbers, and append them to the solutions. *)
solutions = Join[
solutions, {#, Fibonacci@n} & /@
Select[decompositions,
AllTrue[#, MemberQ[fibs, #] &] && Unequal @@ # &]]
];
solutions
(* {{{2}, 3}} *)
`

I admit there's a lot of syntactic sugar in there, so let me know if anything is unclear.