@ecterrab Wow, that was shortest time between a bug report and a published patch for it that I've ever seen from Maple. Good job!

@kencom1 I thought that maybe they were intended to be functions. In that case, they should be defined as "arrow expressions":

L[1]:= f-> diff(f, eta$3) - diff(f, eta)
L[2]:= theta-> diff(theta, eta$2) - theta

Now, it doesn't matter that f and theta are re-used, because any variables on the left side of an arrow are local to that arrow expression.

Note that this WON'T work: L[1]:= f-> diff(f(eta), eta$3)  - diff(f(eta), eta). I mention it because it is closer to what you had.

Remember that now you'll need to take away the extra spaces that you may have added because of the previous Answer.

When you see unexpected tables appearing in error messages, the cause is usually attempting to use a "parent" name after an assignment has been made to one of its indexed "children". The assignment turns the parent into a table.
 

What is the range of your parameter A? Are complex solutions for x3 and x4 acceptable?

@torabi If V(t) is constant, then why are you interested in plotting it? Or did you mean that you want to plot Z(t) vs. t?

@David Sycamore 

Maple's Li is the integral from 0, so Li(x) - Li(2) is the integral from 2. Neither of these functions has any singularities for x > 2, which is the only useful domain for our purposes.

The issue with .../2/p vs. .../2*p has nothing to do with generating the data. It only affects the asymptotic function, so, only the plot. 

"Decoupling" the code that counts the pairs from the code that generates the plot is utterly, utterly trivial, so I encourage you to just do it, and if you do it incorrectly, I''ll let you know.

The thing with the r-tuples is vaguely interesting. Since it's simple to code, I'll do it. It will be slower than the code for pairs, but not horribly slow. It is the "independence" that makes these things less than fully interesting.

@David Sycamore 

There is a certain aesthetic to asymptotic expansions. I automatically think of the terms being "nice" functions. Li(x), x > 2, is a very nice function---strictly monotonic, analytic, and capable of being quickly evaluated to arbitrary precision for arbitrarily large arguments; whereas pi(x) is a nasty function with none of those properties. But we know pi(x) ~ Li(x). (*1)

Maple's Li is the Li used in analysis. The Li used in number theory requires a tiny adjustment: It's Li(x) - Li(2). This doesn't change the asymptotic relationship (*1).

Now here's the easy rest of that derivation that I left for the interested reader:

Let's suppose p > 3 is prime. The only two inverse pairs with a repeated number are (1,1) and (p-1,p-1), neither of which contain a prime. There remain (p-3)/2 pairs to check, each being a pair of distinct values. The expected number of prime pairs is

"number of pairs" * "probability that an arbitrary pair contains two primes". (*2)

Assuming the "asymptotic independence", (*2) becomes

(p-3)/2 * "prob. 1st member of pair is prime" * "prob 2nd member is prime" =

(p-3)/2 * pi(p-2)/(p-3) * (pi(p-2) - 1)/(p-4) =

pi(p-2)*(pi(p-2) - 1) / (2*(p-4)) ~

(Li(p-2) - Li(2))*(Li(p-2) - Li(2) - 1)/2/(p-4)

This is a slightly sharper estimate than I had before. Note that the singularity at p = 4 is removable. I can use pi(p-2) because it's guaranteed that p-1 is composite.

@Josci95 

The command randomize(42) does the same thing globally. So, the seed option to GenerateGaussian is superfluous. Using randomize() with no argument sets the seed to a value based on the system clock.

When the GenerateGaussian command was introduced, I immediately noticed that its code contained a randomize() call. If this had stayed, it would've prevented the possibility of doing exactly the thing that you want to do---generate reproducible random numbers. I complained, and MapleSoft immediately fixed the command and someone threw in that seed option, I guess not realizing that it was superfluous. That's why the seed option is not mentioned in the top section of the help page ?GenerateGaussian---because it was added at the last minute.

What is the size and datatype of your matrices? Even for dense matrices of hardware floats with millions of entries, getting the generalized eigenvectors on my computer takes less than 10 seconds:

(FF1,FF2):= 'LinearAlgebra:-RandomMatrix(1600,1300,datatype= float[8])'$2:
CodeTools:-Usage(
   LinearAlgebra:-Eigenvectors(LinearAlgebra:-MatrixInverse(FF2, method= pseudo).FF1)
);

memory used=260.99MiB, alloc change=247.71MiB, cpu time=8.89s, real time=4.47s, gc time=0ns

The fact that the real time is significantly less than the cpu time indicates to me that there's already some parallel programming being used.

I promised you an algorithm for the inverse-mod-p based on the Bezout coeffcients using only simple arithmetic. Here it is:

`mod/Inv`:= proc(a::posint, p::posint)
description "computes 1/a mod p";
local r0:= p, r1:= modp(a,p), r, t0:= 0, t:= 1;
   while r1 > 1 do
      (t0, t, r0):= (t, t0 - t*iquo(r0, r1, 'r'), r1);
      r1:= r
   od;
   `if`(r1 = 0, FAIL, modp(t,p)) #FAIL if gcd(a,p) <> 1
end proc:

Even if p isn't prime, there's still a multiplicative group mod p. The above code works in that case also.

In terms of efficiency, Inv(a) mod p isn't bad, but 1/a mod p is faster. I think that's simply because it uses builtin code.

@vv I think that you meant "For k = p - 2 one obtains the inverses mod p".

@David Sycamore 

Your error message suggests to me (although I'm not sure about it) that you've entered the code into Maple's 2D Input mode and that you've entered extra spaces. If you copy-and-paste it directly from the post to 2D Input (without adding any spaces!), I think that it'll work, but I don't know for sure because I never use 2D Input. Your best bet is to simply never use 2D Input. Instead, use Maple Input (aka 1D Input). Anyway, I doubt that it'd be possible to transcribe 2D code into OEIS.

a/b/c is equivalent to a/(b*c). I use the former because it's fewer characters. However, a/b/c is NOT equivalent to a/b*c. So, if you've changed my .../2/p to .../2*p, then you'll need to correct that.

The n in procedure CountPairs that appears stray is the procedure's return value. An expression that appears immediately before the end of a procedure is usually (although not 100% always) its return value.

So, I guess that for the purposes of OEIS you want code that lists the first several values of the sequence? Then immediately after the definition of CountPairs (and before try), put

seq(CountPairs(ithprime(k)), k= 1..23);

You can change 23 to whatever positive integer you want.

@David Sycamore The matrix PP has what you want. It has two columns: the first is the primes, and the second is the corresponding count of pairs. This data is also in the remember table of CountPairs. So, if you enter CountPairs(nextprime(10000)), the number of pairs will be returned instantaneously, because the procedure remembers having computed that already.

What does it matter whether you look at it before or after the plot?

The matrix PP is a "hardware float" matrix simply because that plots faster. That means that its entries have decimal points even though they're meant to represent integers. If you don't like that, you can convert it to a matrix of integers by

IntPP:= map(trunc, PP);

@waseem Ah, yes, that paper, which we've discussed at length in another thread. That one (which I'll call Gupta after the lead author's surname) is the most error-filled of all the papers in this field that I've tried to work through. Even the symbolic presentation of the BVP system itself (which occurs twice in the paper) contains an error (the same error, twice).

This Post may not be the appropriate place to discuss this further, but here goes anyway....

I do not understand the "variational FEM" algorithm presented in the paper. You've shown me other papers by the same author describing the algorithm, but unfortunately they all have a nearly identical presentation of it. The only additional information that I've learned is that evenly spaced mesh nodes have been used.

Unlike the BVP from the paper that I analyzed in this Post, Maple's BVP solver does have significant trouble solving Gupta's BVP (after Preben and I corrected it) for some configuations of parameters. So there might be some value in his "variational FEM" algorithm. However, I'll note that for some other configurations of parameters where Maple's solver converged, it only matched Gupta's numbers to 1 or 2 significant digits. I'll guarantee you that Maple's numbers are right and Gupta's are wrong! For one thing, his method has no error control or convergence check, and he doesn't even claim that it does. (Yet he mysteriously claims that his results are accurate to four decimal places.) And he's affiliated with a department of mathematics; shame on him!

Why won't Gupta provide his code? Until he does, his papers are pretty much worthless, and their only purpose is to inflate the CVs of him and his cronies. (Don't papers get refereed these days? How can a referee work without the code?)

That's not to say that Gupta's BVP isn't worth solving. It is. Without access to his code, I think the best hope is to use a "bootstrapping" approach with Maple's solver. That is, you build up a series of successively better approximate solutions to pruned-down BVPs which you feed back into the solver until you have one for which the original unpruned BVP converges. Once you have convergence, the numeric results will be accurate; there'll be no worries about that.

Like I said, this is not the appropriate place to discuss this further. If you have something to contribute towards a solution, please post it in the thread where we were discussing the Gupta paper.

@ecterrab Thank you, Edgardo; your praise means a lot to me.

@David Sycamore Here's the retrofitted code. As I said, the retrofit was trivial. In the future, please put your Maple version in Question headers using the pull-down list provided for that purpose.

#This code should run in Maple 2017. Let me know if you have any trouble with it.

CountPairs:= proc(p::prime)
option remember;
local b, a:= 2, n:= 0;
   while a < p do
      b:= 1/a mod p; 
      if a <= b and isprime(b) then n:= n+1 fi;
      a:= nextprime(a)
   od;
   n
end proc:

#Compute as many as we can within a certain time limit:
try
   timelimit(
      999, #seconds
      #This infinite loop is intentional:
      proc() local p:= 1; do p:= nextprime(p); CountPairs(p) od end proc() 
   )
catch:
end try:

#Extract remember table and convert to a Matrix:
PP:= Matrix(
   [lhs,rhs]~([entries(op(4, eval(CountPairs)), pairs, indexorder)]),
    datatype= float[8]
):
N:= max(PP):
Asy:= p-> Li(p)^2/2/p: #1st asymptotic term
Res:= p-> 2*sqrt(p)/Pi^2: #estimated residual bound
plot(
   [PP, Asy(n), Asy(n)+Res(n), Asy(n)-Res(n)], n= 2..N, 
   style= [point,line$3], color= [grey,red,blue$2], thickness= [0,1,0,0], 
   symbolsize= 1, symbol= point
);