This is the same as my last column graph but with a logarithmic color scale for the probabilities. With this, you can see the many orders of magnitude of variation among the vast "field" of values  (in front of the "mountain") that originally just appeared to be close to 0.

p1:= 1/5:  p2:= 1/2:  p3:= 1-p1-p2:  n:= 20:
Multinom:= table([seq](seq(
    (i,j)= n!/i!/j!/(n-i-j)!*p1^i*p2^j*p3^(n-i-j), i= 0..n-j), j= 0..n)
):
E:= [evalf({entries}(Multinom, nolist))[]]:
(m,M):= [min,max](ln~(E))[]:
plots:-display(
    seq(seq(
        plottools:-cuboid(
            [i-2/5, j-2/5, 0], [i+2/5, j+2/5, (p:= Multinom[i,j])], 
            style= surface, 
            color= COLOR(HSV, 0.85*((ln(p)-m)/(M-m)), 1, 1)
        ),
        i= 0..n-j), j= 0..n
    ),
    plot3d(
        0, (i,j)=~ -1/2..n+1/2, grid= [(n+2)$2],
        style= wireframe, color= black, thickness= 0.6
    ),
    labels= ["Cat. 1", "Cat. 2", "Prob."], labelfont= [times, bold, 12],
    axis[1,2]= [tickmarks= [seq](x= `if`(x::odd,``,x), x= 0..n)],
    axes= frame, glossiness= 1, 
    projection= 0.8, orientation= [-125, 65], lightmodel= light1,
    view= [(-1/2..n+1/2)$2, default]
);
#Color legend:
plots:-display(
    seq(
        plot([[p,0],[p,1]], color= COLOR(HSV, 0.85*((ln(p)-m)/(M-m)), 1, 1)),
        p= E
    ),
    axis[1]= [mode= log], axis[2]= [tickmarks= []], size= 100*[7,2], 
    thickness= 3, title= "Logarithmic Color Scale:\n", 
    titlefont= [times, bold, 12],
    labels= [`probabilities by color\n`, ``], labelfont= [times,10]
);

@sand15 You are not using the view option that I included:

view= [min(X)-W..max(X)+W, 0.1..max(Y)]

The purpose of its X component is to avoid losing a bar on the left or right. I believe that that will work in all cases. The purpose of its Y component is to avoid seeing any "0 height" bars between the true bars.

And how does this have any significant effect on the length of the plot structure?

Here is another way to make a column graph similar to the matrixplot from raw components. I like this better than the matrixplot because like the trianglar-base plots shown, it avoids the area i+j > n (i.e., the part outside the support of the distributuion).

p1:= 1/5:  p2:= 1/2:  p3:= 1-p1-p2:  n:= 20:
Multinom:= (i,j)-> n!/i!/j!/(n-i-j)!*p1^i*p2^j*p3^(n-i-j):
plots:-display(
    seq(seq(
        plottools:-cuboid(
            [i-2/5, j-2/5, 0], [i+2/5, j+2/5, Multinom(i,j)], 
            style= surface, transparency= 0.15
        ),
        i= 0..n-j), j= 0..n
    ),
    plot3d(
        0, (i,j)=~ -1/2..n+1/2, grid= [(n+2)$2],
        style= wireframe, color= black, thickness= 0.6
    ),
    labels= ["Cat. 1", "Cat. 2", "Prob."], labelfont= [times, bold, 12],
    axis[1,2]= [tickmarks= [seq](x= `if`(x::odd,``,x), x= 0..n)],
    axes= frame, glossiness= 1, 
    projection= .8, orientation= [-125, 65], lightmodel= light1,
    view= [(-1/2..n+1/2)$2, default]
);

@sand15 I totally agree that no fixed setting of binwidth will work adequately over a variety of cases. But likewise no fixed setting of thickness will work either. A further complication of thickness is that it's not measured in any units that are meaningful within the plot's coordinate system. For example (I found these by trial-and-error), in your first plot, thickness= 70 is equivalent to binwidth= 1; whereas in your 2nd plot (after commenting out randomize()), thickness= 39 is equivalent to binwidth= 1.

Here is, I think, a reasonable way to choose the binwidth, as you asked at the end of your worksheet. I assume that the data are given as equal-length lists X and Y with all X-values distinct.

# What binwidth to choose in Carl Love's code?
Xs:= sort(X); W:= min(Xs[2..] -~ Xs[..-2]);
Statistics:-Histogram(
    `$`~(X, Y), frequencyscale= absolute, binwidth= W, 
    axis[1]= [tickmarks= (Xs=~ typeset~(Xs))],
    view= [min(X)-W..max(X)+W, 0.1..max(Y)],
    axesfont= [times, bold, 9], axes= frame
);

You need to add frequencyscale= absolute to your earlier histograms to achieve what the OP wants.

While randomize() is essential for "production" code, it's not great when generating examples because they're not reproducible unless you display the "key" that is the value returned by the command randomize().

Regarding tickmarks: This option works the same way for all plotting commands. If a tickmark is specified as an "equation", the left side must be a number in the background coordinate system for its axis. The coordinate system for an n x n matrix in a matrixplot is (1..n) x (1..n). The right side of the "equation" can be anything, not necessarily a number, and is what should be printed at the position indicated by the left side. (I put "equation" in quotes to indicate that I simply mean a Maple expression with an =, not a mathematical equation.) So, for example, 1.5= 0 means to display a tickmark labeled 0 at the position that was originally 1.5.

Regarding small histograms: To make a small histogram such as you asked for: 

Statistics:-Histogram(
    (`$`@op)~([[1,4],[2,3],[4,5]]), frequencyscale= absolute, binwidth= 1
);

Actually, even the combine is not necessary. It might help a human see what's going on, but is can do the verification without it:

Download CosProductInduct.mw

@dharr I tried the is with assuming n::posint, x > 0, x < Pi/2, but I still got false. It's not at all surprising to me that is cannot prove the identity, and I wouldn't even expect it to be able to. But IMO it should never return false without a witness. (A witness is an assignment to the variables under universal quantification that fulfills their assumptions and for which the identity is definitively false.)

The fully symbolic induction can be done with a single combine command, like this:

Download CosProductInduct.mw

@dharr My comments regarding the is and coulditbe statements were strictly regarding those statements in the OP's followup entitled "p. s.". It looks to me as if he took term1 from the Question and substituted x= t/2^n, but also took term2 and substituted x= t/2^(n-1). Thus the is(term1 = term2) correctly returns false in "p. s.". The false return from is in the Question is a problem. If it can't prove the identity, but also finds no counterexample, it should return FAIL.

@vv Thank you; I stand corrected. I was previously aware of this distinction, but I forgot to apply it here. Nonetheless, the graphical evidence of the (nearly) space-filling curves suggests to me that in this case those are indeed the liminf and limsup.

Regarding the limit that you asked for at label (5) in the original Question: Clearly the limit proper doesn't exist. The expression is a (nearly) space-filling curve for most real x. By giving appropriate assumptions, limit will return the liminf and limsup (wrt x). This is a special case where this works. (As @vv pointed out, in general it doesn't work.) The cos(2^(n+1)*x) is essentially a uniform random-number generator on [-1,1]. So the numerator oscillates between -1 and 3; hence, the range answer given below.

restart:
T:= (n,x)-> (1+2*cos(2^(n+1)*x))/(1+2*cos(2*x)):
plot(T(n,Pi/6), n= 0..999); #space-filling curve
limit(T(n,x), n= infinity) assuming x>0, x<Pi/4;

Once again, I took the challenge of seeing how many such triples I could find with 1 second of CPU time. To make things more efficient, note that the first prime in any such triple must be of the form 6*k+1. I'll leave it to you to prove why that's true. I found 1694 such triples. Below, I only list the first 9 and last 9.

restart
:
p:= 1:
T:= 1:  Stop:= T + time():
R:= table():
while time() < Stop do
    p:= nextprime(p);
    if irem(p,6)=1 then
        if isprime(p+4) and isprime(p+6) then
            R[p, p+4, p+6]:= ()
        fi;
        p:= p+5
    fi
od:
L:= [indices](R, indexorder):
nops(L);
L[..9];
L[-9..]; 

                              1694

[[7, 11, 13], [13, 17, 19], [37, 41, 43], [67, 71, 73], [97, 101, 103], 
 [103, 107, 109], [193, 197, 199], [223, 227, 229], [277, 281, 283]]

[[1207117, 1207121, 1207123], [1208017, 1208021, 1208023], [1210393, 1210397, 1210399],
 [1210873, 1210877, 1210879], [1211593, 1211597, 1211599], [1211653, 1211657, 1211659],
 [1212433, 1212437, 1212439], [1212847, 1212851, 1212853], [1213627, 1213631, 1213633]]

Your wording seems to imply that you think that using isprime is somewhat equivalent to using ifactor and seeing if there is one or more than one factors. If the number being checked is indeed prime, that might be true. But if it's a large composite, isprime will determine that in seconds, whereas ifactor could take years or centuries. If it weren't for this fundamental difference between primality testing and integer factorization, moderm cryptography would be impossible.

Some improvements to your algorithm: 

1. For any odd positive p (prime or not), 2^p+1 is divisible by 3, so you don't need to check whether it's divisible. Indeed, 

(2^p+1)/3 = sum((-2)^k, k= 0..p-1),

which proves the divisibility. However, that formula is not computationally efficient, so I don't use it in the code below.

2. Here's the efficient formula: For positive odd p, let W(p) = (2^p+1)/3. Let q>p be odd and positive. (The formula below doesn't depend on p and q being prime; however, they always are in my code below.) Then 
W(q) = 2^(q-p)*W(p) - floor(2^(q-p)/3).
So, each Wagstaff number can be obtained by updating the previous one. Since q-p is always the gap between consecutive small primes, the only large number in the computation is W.

3. You don't seem to know about the Maple command nextprime, which makes it easier to iterate over primes.

By those methods, I can find the first 21 Wagstaff primes using just 1 second of CPU time:

restart
:
#Procedure to elide middle digits from the display of long integers:
Elide:= proc(n::posint, w::posint:= interface(screenwidth), ht::posint:= 20)
local L:= length(n) - 2*ht, r;
    if L < w then sprintf("%d", n)
    else 
        sprintf(
            "%d...[%d digits elided]...%0*d", 
            iquo(iquo(n, 10^ht, r), 10^L), L, ht, r
        )
    fi
end proc
:
#T is the number of seconds of cpu time to allow.
T:= 1:  Stop:= T + time():  
p:= 1:  W:= (2^p + 1)/3:  nW:= 0:
printf(
    "\nExponent Wagstaff prime"
    "\n======== ==============\n"
);
while time() < Stop do
    q:= nextprime(p+1);
    M:= 2^(q-p);
    p:= q;
    W:= M*W - iquo(M,3);
    if isprime(W) then nW:= nW+1; printf("%7d  %s\n", p, Elide(W,64)) fi
od:

printf(
    "\nLargest prime exponent checked: %d"
    "\nNumber of prime exponents checked: %d"
    "\nNumber of Wagstaff primes found: %d\n",
    p, numtheory:-pi(p), nW
); 

Exponent Wagstaff prime
======== ==============
      3  3
      5  11
      7  43
     11  683
     13  2731
     17  43691
     19  174763
     23  2796203
     31  715827883
     43  2932031007403
     61  768614336404564651
     79  201487636602438195784363
    101  845100400152152934331135470251
    127  56713727820156410577229101238628035243
    167  62357403192785191176690552862561408838653121833643
    191  1046183622564446793972631570534611069350392574077339085483
    199  267823007376498379256993682056860433753700498963798805883563
    313  5562466239377370006237035693149875298444543026970449921737087520370363869220418099018130434731
    347  95562442332919646317...[64 digits elided]...17439712921903606443
    701  35067572676989156714...[171 digits elided]...78180862167823854251
   1709  96192527248701069005...[474 digits elided]...12036857299070528171

Largest prime exponent checked: 2351
Number of prime exponents checked: 349
Number of Wagstaff primes found: 21

@C_R I think so, but I'm not sure. I think that you would not be allowed to create another MaplePrimes account from the same email address.

@acer I didn't have the trouble using Android+Chrome (latest automatically updated versions) that you did.