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
);

As you probably realize at this point, the inverse of a in Zp can be computed simply by

b:= 1/a mod p;

and this is not done by the "brute force" method of checking every possibility. It can be done as a byproduct of computing the gcd of a and p by the Euclidean algorithm. Of course, this gcd is 1; but with slight additions to the standard algorithm, the Bezout coefficients can also be computed. The Bezout coefficient of a, reduced mod p, is the inverse.

The command igcdex will tell you the Bezout coefficients, but unfortunately the method that it uses is totally fake and inefficient. I'll post a decent algorithm based on first principles in a little while.

@vv Yes, of course they aren't independent for fixed p.  If they were, there'd be little point in doing these computations. I just lacked the mathematical vocabulary to state my asymptotic intuition precisely. Do you have any intuition about it, or about how difficult it would be to prove? Do you have any intuition about the residual bounds?

@acer 

I accept your point about the obsfuscation clouding my judgment, but I'll point out that what Christian did could only be considered obsfuscation because of the simplicity of the example. Otherwise, what he did would just be considered normal procedure coding without raising an eyebrow.

I very often (daily) write procedures that take procedures as arguments. In almost all cases, I'd like that code to continue to work if the argument is replaced by an expression that can be treated like a procedure. I don't like the idea of handling the distinction via a local variable in the parent procedure. I'd rather place the onus on the one passing the procedure-like expression to do something to ensure correct behavior. That something is usually making it a procedure. This can usually be done ad hoc in the same line of code. Here's an example that I posted last night, in response to another Question:

MetaOp:= (gen, n::nonnegint, pred, EA)-> ()-> EA(select(pred, ['gen()'$n])): 

#The original, erroneous, code:
Op1:= MetaOp(rand(1..100), 20, x-> x::even, [add, add/nops]):

#The corrected code:
Op2:= MetaOp(rand(1..100), 20, x-> x::even, x-> (S-> (S, S/nops(x)))(add(x))):

The only difference between the two being the fourth argument, passed to MetaOp's parameter EA. As you can see, I made it into a double procedure, which both corrects the issue being discussed in this thread and avoids duplicating the computation done by add.

You can see in MetaOp that the obsfuscation (of the argument-distribution issue) is unintentional, yet it's several levels deeper than Christian's example: the random-numerator generation survives as an unevaluated expression all the way until it's passed to add (twice) and nops.

@Nata Red Grand 

I realize that the example of an external accountant procedure that I gave in my Answer---

x-> (S-> (S, S/nops(x)))(add(x))

---is a bit over-the-top, but that's because my goal was exploring the possibilities rather than directly doing your homework. You, on the other hand, should still be capable of writing a simpler external accountant procedure and passing it in. Here's an example:

SumAndAverage:= proc(Data::list(algebraic))
local S:= add(Data), n:= nops(Data);
   'Sum' = S, 'Average'= `if`(n=0, undefined, S/n)
end proc:

In addition to what Joe said: Unless the package is huge (like more than, say, 100,000 lines of code), there's no noticeable benefit to "updating" as opposed to "entirely recreating from the source code", because the Maple parser/"compiler" is so fast.

@Kitonum I'm not exactly sure what the OP's teacher means by "external accountant", but my guess is that they mean that there should be a separate procedure that does the sum and average and that that procedure should be passed into the final procedure.

One very small potential problem with your code is that it may fail if L1 is never created as a table because none of the numbers are even. You can get around that by explcitly declaring L1:= table() at the start.

@Christian Wolinski 

I like your example because I can look at it and immediately say Surely Op and Op2 ought^[*1] to produce equivalent output for any input for which they both produce output (not counting error messages as output). If there were cases where Op returned an error message or warning while Op2 returned regular output, I wouldn't be as upset.

^[*1]I mean ought in the moral sense, as I always mean it.

@Joe Riel I see what you're saying, Joe, and perhaps if all examples were like Acer's example---the direct application of an expression as if it were a procedure---then I wouldn't be so concerned about it. But in my indirect example, I want the second argument of Op to be a procedure or anything that can be used as if it were a procedure.