# Find a single max value

w[n_] := Expand[Sum[Binomial[n - k - 1, k]*(-1)^k*A^(n - 2*k - 1), {k, 0, n - 1}]]
f[x_, y_, z_] :=PolynomialRemainder[(w[z] - 1)*(w[y] - 1), (w[x] - 1), A]

For[i = 3, i < 450, i++,For[j = 3, j < 450, j++,For[k = 3, k < 450, k++,If[i < j < k,
Print[{i, j, k},
N[Max[FindInstance[{Abs[f[i, j, k]] - Abs[w[i] - 1]} == 0 &&
3 <= A, {A}]]] ]]] ]]


Above code gives max A for all each cases. My aim is to find max A for all cases.

Q2: How can i reduce the time for above code? Thank you.

-
sorry, i forgot to add the functions. – MATIRMAK Apr 6 '13 at 11:20
Now it makes more sense, however {Abs[f[i, j, k]] - Abs[w[i] - 1]} == 0 cannot be satisfied, it means a_List == b_Integer, moreover I can't understand how FindInstance can help you here. – Artes Apr 6 '13 at 11:31
@Nurettinırmak please add those functions into the question by editing it. Don't leave them hidden in comments. – Verbeia Apr 6 '13 at 11:37
For example suppose that i=4, j=10, k=12. Therefore, the functions f and w depend on 'A'. Using the 'FindInstance', i can find a solutions of the equation '{Abs[f[i, j, k]] - Abs[w[i] - 1]} == 0' for i=4, j=10, k=12. it runs. i checked it. But, it gives the greatest A for each cases and it seems that it will take a long. – MATIRMAK Apr 6 '13 at 11:46
@Nurettinırmak Does "Max A in all cases" mean max A for all i between 3 and 450 ? Since A is a symbol do you want real or integer A? If it works, could you demonstrate for this simple example i=4, j=10, k=12 how does your code work ? – Artes Apr 6 '13 at 11:56

Artes’ comment suggests there is a problem with the concept of optimum you are trying to find. Let me make the more general point about how your code may be improved.

• First, if you have a triply nested For loop, you are probably doing it wrong. Consider Table instead, or some of the other constructs mentioned here.

• Second, you don't need to go through every combination of {i, j, k} since you only want the cases where i<j<k. I think you can get this set of combinations using Subsets[Range[3,450],{3}] which gives the ordered subsets of exactly length 3.

• Third, functional programming is usually more efficient. Create the list of indices and then Map the desired function to the elements of that list.

• Fourth, Print doesn't store anything. You need to save the results for each case as part of a list somewhere and then find the maximum value of that list. The mapping- onto-a-list approach will actually store the results, if only temporarily.

Putting all of this together, something like this should work:

Max[yourfunction /@ Subsets[Range[3, 450], {3}]]

-
@Verbeia, Thank you very muck for your attention and your answer. But, Unfortunately, i'm a new user. Therefore, your method is ver complicated for me. – MATIRMAK Apr 6 '13 at 13:11

I tried to make a compilable version by converting the polynomial remainder function to a linear map (matrix multiplication). But I didn't think about how large the coefficients were (up to about 90 or so digits). So the compiled function eventually detects an integer overflow, bails out and proceeds with the uncompiled version.

Finding the remainder involves creating a matrix that maps each power of A to its remainder. This can be done recursively for each w[n] (see polyrem0). The matrix actually operates on a polynomial's CoefficientList. The thought was that for most w[n], this will be used many times, and, in addition to computing it only once for each n, some gain in speed might be realized from internally optimized matrix multiplication. The larger than machine size integers raises doubts about that. Perhaps another basis for the polynomials would keep the coefficients in polyrem0 small. There is some reason to think this is possible, as randomly tested remainders do not seem to have extremely large coefficients. The w[n] themselves seem worth testing as a basis. I may try it myself, when I have more time.

$maxWIdx = 50; (* maximum index n to test for w[n]; can be changed in Block[] *) w[n_] := Expand[Sum[Binomial[n - m - 1, m]*(-1)^m*A^(n - 2*m - 1), {m, 0, (n - 1)/2}]]; f[{x_, y_, z_}] := f[x, y, z]; f[x_, y_, z_] := CoefficientList[(w[z] - 1)*(w[y] - 1), A]. Take[polyrem[x,$maxWIdx + $maxWIdx - 3], y + z - 1].A^Range[0, x - 2]; polyrem0 = Compile[{{coeff, _Integer, 1}, {maxDeg, _Integer}}, IdentityMatrix[Length[coeff] - 1]~Join~ NestList[Most[{0}~Join~# - #[[-1]] coeff] &, -Most@coeff, maxDeg - Length[coeff] + 1] ]; polyrem[n0_Integer] := polyrem[n0,$maxWIdx + $maxWIdx - 3]; polyrem[n0_Integer, maxDeg0_Integer] := (Print[{n0, maxDeg0}]; polyrem[n0, maxDeg0] = polyrem0[CoefficientList[w[n0] - 1, A], maxDeg0]);  Test up to 50 (3 <= x < y < z <=50): Max A is 16.9443. (222 sec.) DistributeDefinitions[w, f, polyrem]; Block[{$maxWIdx = 50},
polyrem /@ Range[3, $maxWIdx - 2]; (* initialize *) Max[ParallelMap[ Max@(A /. NSolve[{Abs[f[#]] - Abs[w[First[#]] - 1]} == 0, {A}, Reals]) &, Subsets[Range[3,$maxWIdx], {3}]]]
] // AbsoluteTiming
(* {222.149237, 16.9443} *)


The maximum seems to occur when the first argument is around 5. I don't know why. (Oops: I should have set $maxWIdx = 449.) Max A (x==5 < y < z <= 450) is 150., curiously nice number! (575 sec. Ouch.) DistributeDefinitions[w, f, polyrem]; Block[{$maxWIdx = 450},
polyrem[5]; (* initialize *)
Max[ParallelMap[
Max@(A /. NSolve[{Abs[f[5, Sequence @@ #]] - Abs[w[5] - 1]} == 0, {A}, Reals]) &,
Subsets[Range[6, \$maxWIdx], {2}]]]
] // AbsoluteTiming
(* {574.866569, 150.} *)


It is a rather large problem, with large coefficients. Nearly 15 million equations to solve, and it took 4 sec. just to generate all the indices:

Subsets[Range[3, 449], {3}] // Length // Timing
(* {4.024577, 14786015} *)


Using approximate Reals, which is often faster, has a drawback. When the type of coeff in the compliled polyrem0 was changed to _Real, the maximum A returned was 664.01. Computing the polynomial remainder involves subtracting large numbers, and I would guess signifcant error is introduced.

-
Thank you for response. But, unfortunately the function w was defined falsely. – MATIRMAK Apr 8 '13 at 13:27
Do you mean summing only up to (n-1)/2? Further terms are zero, aren't they? – Michael E2 Apr 8 '13 at 13:31
Sorry.yes, you are right. And do you say that max A is 150. Right? – MATIRMAK Apr 8 '13 at 13:36
Max A is 150. for x == 5 and 5 < y < z <= 450. I didn't check all cases for x, since that one case took 575 sec. on a quad-core i7. – Michael E2 Apr 8 '13 at 14:44

Up to 50:

r = Quiet@ Cases[Evaluate[({{##},
FindRoot[Abs[f[#1, #2, #3]] - Abs[w[#1] - 1], {A, 3}]} & @@@
Subsets[Range[3, 50], {3}])], {x_, (y_ /; (A /. y) >  3)} :> {x, (A /. y)}];
ListPointPlot3D[r[[All, 1]]]


-
.Thank you for your answer. still, i have no idea for how to find the largest A – MATIRMAK Apr 6 '13 at 14:21