(* Produces number of seeds with n of this form. *) squareSeedCount[d_, prime_] :=   If[IntegerQ[Sqrt[d]] || (Mod[d, prime] == 0 && IntegerQ[Sqrt[d/prime]]),    Module[{i, w, y},     If[IntegerQ[Sqrt[d]], w = prime*Sqrt[d]*i; y = prime*i^2,      w = Sqrt[d/prime]*i; y = i^2];     Sum[If[w < d && y < d && GCD[w, y, d] == 1, 1, 0], {i, 0,       Sqrt[d/prime]}]], 0]; (* Produces number of seeds with n of this form. *) multipleOfPSeedCount[d_, prime_] :=   Module[{factorization, squarelist, finallist, i},    factorization = FactorInteger[d];    If[Length[factorization] > 1 && Mod[d, prime] == 0,     squarelist = MapAt[2 # &, 2] /@ factorization;     squarelist[[First@Flatten@Position[squarelist, {prime, _}]]][[       2]] = squarelist[[         First@Flatten@Position[squarelist, {prime, _}]]][[2]] - 1;     finallist = Times @@@ Subsets[Power @@@ squarelist];     Sum[If[(finallist[[i]] <          d ) && (finallist[[Length@finallist - i + 1]] < d), 1,       0], {i, 1, (Length@finallist)/2}], 0]]; (* Produces total number of new solutions at divisor d of n *) SolAtDivisorNCount[d_, prime_] := (SolAtDivisorNCount[d, prime] =     8 (multipleOfPSeedCount[d, prime] + squareSeedCount[d, prime])); (* Produces first 100 terms with prime 11 by calculating new solutions at each natural number <= 100. *) With[{n = 100, prime = 11},   Accumulate[Join[{1}, Table[Total[Function[d, SolAtDivisorNCount[d, prime]] /@ Divisors[i]] - 4, {i, 1, n}]]] ]