A358404 - OEIS

%I #12 Mar 20 2023 19:43:39

%S 2,3,5,8,13,21,23,34,41,61,85,89,144,233,255,264,377,383,397,443,610,

%T 762,875,987

%N Multipliers involving Fibonacci-like sequences and Pythagorean triples.

%C A positive integer m is an element of this sequence if and only if there exists a Pythagorean triple of the form (m*G_0, m*G_1, G_n), where (G_k) is a Fibonacci-like sequence, i.e., a sequence with arbitrary positive integer starting values G_0 and G_1 and satisfying the recurrence G_k = G_{k-1} + G_{k-2} for every index k > 1.

%F Each element of this list has a unique representation of the form m = m(c, j) = G_n / c, where j is an arbitrary nonnegative integer and c is "good", meaning that all of its prime divisors are of the form 4k + 1 and the Fibonacci entry point t of c is odd, in which case n = ((2j + 1)t + 1)/2 and (G_0, G_1, c) is the unique primitive Pythagorean triple such that G_0/G_1 is congruent to F_n/F_{n-1} modulo c.

%e m = m(5, 0) = 2, since the Fibonacci-like sequence (G_n) with G_0 = 4 and G_1 = 3 has G_3 = 10 and (m*G_0, m*G_1, G_3) = (8, 6, 10) is a Pythagorean triple. Since m = 2 is the smallest positive integer with this property, m(1) = 2.

%t (* Fibonacci entry point *)

%t T[m_] :=

%t   Module[{fi = FactorInteger[m], lenN, i, fi2, p, e, q, n1, divs,

%t     nDivs, d, found, preres, result = 1},

%t    If[m == 1, Return[1]];

%t    len = Length[fi];

%t    {p, e} = fi[[1]];

%t    q = p^e;

%t    If[len == 1,

%t     If[p == 5, Return[q]];

%t     If[e == 1,

%t      result = p - JacobiSymbol[p, 5];

%t      While[EvenQ[result] && Mod[Fibonacci[result], m] == 0,

%t       result /= 2];

%t      If[Mod[Fibonacci[result], m] != 0, result *= 2];

%t      fi2 = FactorInteger[result];

%t      If[EvenQ[result], Drop[fi2, 1]];

%t      n1 = Product[fi2[[i, 1]]^fi2[[i, 2]], {i, Length[fi2]}];

%t      divs = Divisors[n1];

%t      nDivs = Length[divs];

%t      found = False;

%t      For[i = 2, i <= nDivs && ! found, i++,

%t       d = divs[[i]];

%t       If[Mod[Fibonacci[d], m] == 0,

%t        found = True;

%t        result = d;

%t        Return[result];

%t        ];

%t       ],

%t      result = p^(e - 1 - If[p == 2 && e > 2, 1, 0])*T[p];

%t      Return[result];

%t      ],

%t     result = LCM[T[q], T[m/q]];

%t     ];

%t    result

%t    ];

%t (* Good moduli *)

%t GoodQ[m_] :=

%t   Module[{fi, len, i, p, t},

%t    If[m < 5, Return[False]];

%t    fi = FactorInteger[m];

%t    len = Length[fi];

%t    For[i = 1, i <= len, i++,

%t     p = fi[[i, 1]];

%t     If[Mod[p, 4] != 1, Return[False]];

%t     ];

%t    True

%t    ];

%t {* Great moduli *)

%t GreatQ[m_] := GoodQ[m] && OddQ[T[m]];

%t (* Fibonacci modular ratio *)

%t R[c_, k_] :=

%t   Module[{f0 = Fibonacci[k], f1},

%t    If[GCD[f0, c] > 1, Return[$Failed]];

%t    f1 = Fibonacci[k + 1];

%t    Mod[f1*PowerMod[f0, -1, c], c]

%t    ];

%t (* Starting pair for Fibonacci-like sequence *)

%t StartingPair[c_] :=

%t   Module[{pr, len, i, r0, t, n, r, u, v, g0, g1, preres},

%t    If[! GreatQ[c], Return[$Failed]];

%t    t = T[c];

%t    n = (t + 1)/2;

%t    r0 = R[c, n - 1];

%t    pr = PowersRepresentations[c, 2, 2];

%t    len = Length[pr];

%t    For[i = 1, i <= len, i++,

%t     {u, v} = pr[[i]];

%t     If[GCD[u, v] == 1,

%t      r = Mod[v*PowerMod[u, -1, c], c];

%t      preres = {Abs[u^2 - v^2], 2 u*v};

%t      If[r == c - r0, Return[preres]];

%t      If[r == r0, Return[Reverse[preres]]];

%t      ];

%t     ];

%t    $Failed

%t    ];

%t (* Great modulus multiplier *)

%t M[c_, j_] :=

%t   Module[{t, n0, n, g0, g1, result},

%t    If[! GreatQ[c], Return[$Failed]];

%t    {g0, g1} = StartingPair[c];

%t    t = T[c];

%t    n0 = (t + 1)/2;

%t    n = n0 + j*t;

%t    (g0*Fibonacci[n - 1] + g1*Fibonacci[n])/c

%t    ];

%t (* Master table *)

%t MasterTable[mMax_] :=

%t   Module[{c, j, m, g0, g1, t, n0, n, done, result = {}},

%t    For[c = 5, c <= GoldenRatio*mMax^2, c += 4,

%t     While[! GreatQ[c], c += 4];

%t     If[c <= GoldenRatio*mMax^2,

%t      {g0, g1} = StartingPair[c];

%t      t = T[c];

%t      n0 = (t + 1)/2;

%t      For[j = 0, j <= JMax[mMax, n0], j++,

%t       n = n0 + j*t;

%t       m = M[c, j];

%t       If[m <= mMax, AppendTo[result, {g0, g1, c, m, n}]];

%t       ];

%t      ];

%t     ];

%t    result

%t    ];

%t (* Multiplier list *)

%t MList[mMax_] := Union[MasterTable[mMax][[All, 4]]];

%Y Cf. A000045.

%K nonn

%O 1,1

%A _David Terr_, Nov 14 2022