Carl Love

Carl Love

28010 Reputation

25 Badges

12 years, 286 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

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:

term1:= n-> (2*cos(2^n*x)+1)/(2*cos(x)+1):

term2_inner:= k-> 2*cos(2^k*x)-1:

term2:= n-> product(term2_inner(k), k= 0..n-1)
:

#Verify base case:
term1(1) = term2(1); is(%);

(2*cos(2*x)+1)/(2*cos(x)+1) = 2*cos(x)-1

true

#Induction hypothesis:
IH:= (term2 = term1)(n);

product(2*cos(2^k*x)-1, k = 0 .. n-1) = (2*cos(2^n*x)+1)/(2*cos(x)+1)

term2(n+1) = term2(n)*term2_inner(n);

product(2*cos(2^k*x)-1, k = 0 .. n) = (product(2*cos(2^k*x)-1, k = 0 .. n-1))*(2*cos(2^n*x)-1)

subs(IH, rhs(%)) = term1(n+1); is(%);

(2*cos(2^n*x)+1)*(2*cos(2^n*x)-1)/(2*cos(x)+1) = (2*cos(2^(n+1)*x)+1)/(2*cos(x)+1)

true

 

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:

term1:= n-> (2*cos(2^n*x)+1)/(2*cos(x)+1):

term2_inner:= k-> 2*cos(2^k*x)-1:

term2:= n-> product(term2_inner(k), k= 0..n-1):

#Base case:
term1(1) = term2(1); is(%);

(2*cos(2*x)+1)/(2*cos(x)+1) = 2*cos(x)-1

true

#Induction hypothesis:
IH:= term2(n) = term1(n);

product(2*cos(2^k*x)-1, k = 0 .. n-1) = (2*cos(2^n*x)+1)/(2*cos(x)+1)

#Obviously,...
#term2(n+1) =
term2(n)*term2_inner(n);

(product(2*cos(2^k*x)-1, k = 0 .. n-1))*(2*cos(2^n*x)-1)

is(combine(subs(IH,%)) = term1(n+1));

true

 

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.

The Display settings are in the Options dialog/menu. In Maple 2023, Options is accessed from the Tools menu, but in Maple 2025 , Options is the last item on the File menu.

While your typing cursor is in an input field (usually, that's a line with a > on the left), go to the Edit menu. Is there an entry "Switch to Math Mode"? If so, select it. Then try typing something. Does that make any difference, and is it the difference that you want?

The moral of this story is that Maplesoft should test code both with and without the with command.

@dharr Your Answer has some good tips on indexing, but it has nothing to do with the cause of the OP's problem. The error message "invalid subscript selector" can only come from applying an invalid numeric index to a set, list, or expression sequence. In this particular case, in the 3rd iteration of the loop (i=2), dsolve returns NULL. Thus S[2] = NULL, which is an expression sequence, and thus S[2][1] has an invalid subscript selector.

In Maple 2025 (which the OP didn't use), you are getting Typesetting errors in the loop that sets up the conds. There are two easy workarounds for these. The first is to set interface(typesetting= standard). The second is to not use PDEtools:-declare.

In Maple 2025, dsolve works for S[2], but fails for S[3].

@acer The spam button is gone now!

I'm amazed that several hours of work was done on MaplePrimes this morning, yet this Post hasn't been deleted.

1 2 3 4 5 6 7 Last Page 1 of 708