@tomleslie You are a Moderator of MaplePrimes. You have the power to convert Posts to Questions or to edit them in any other way. Please use that power; it only takes 10-15 seconds. Chastising the user also is okay in addition to editing.

@Moh Huda 

It's not possible that the error message that you showed came directly from the worksheet that I posted. It looks like you (probably unintentionally) changed the line c1 = 5e-11 to c1, 5e-11. Anyway, here's a much-updated version of the worksheet. By using method= rosenbrock, I was able to duplicate the curves from the paper in every visual detail.

Let me know how it goes!
 

Download Chemoimmunotherapy.mw

@vv The Remez algorithm is implemented in Maple as numapprox:-minimax. I was going to post an example of it applied to this function, but I didn't have the patience to wait for it to finish.

@Kitonum I believe (although I'm not sure) that you're misinterpretting what the OP means by neighboring points. Under your interpretation, all points are neighboring points of all other points; whereas in my interpretation, they're neighboring iff they're in the same cluster. The term cluster here is not precisely defined, but I think that if you look at the Wikipedia page "Cluster analysis", then you'll see what I mean.

@mmcdara Upon a closer re-reading, I just saw that actually both of your "Not so fast" worksheets show that the Newton's method takes 9+ seconds while Quantile takes 16+ seconds. So either I misinterpretted what you were trying to say, or you switched the measurements in your mind. Either way, the Newton's method is significantly faster.

@Kitonum Note that your Answer applies only to what the OP did in the worksheet, i.e. find a single focal point. That's not the more-general question posed in the 3rd paragraph of the Question.

@mmcdara Ah, you were using Digits = 20, which is outside the domain of hardware floats (aka double precision). My quickly written code had been optimized for hardware floats, since that's quite easy to do (just use evalhf). I just made some changes so that it'll work efficiently regardless of Digits, so please try testing it again (code and worksheet below). And please try on the full interval 0..1, because I made some changes regarding that.

Note that the "magic coefficient" methods have an accuracy far, far less than the 15 digits of hardware floats, as your Post shows, and I suspect that they'd have far more trouble in the extreme tails than my Newton's method.

There's no need for your procedure f that wraps NormInv; it just wastes (a small amout of) time. Use NormInv~ just like Quantile~.
 

Download NormalInveseTest.mw

@Teep @Daniel Skoog

Daniel asked me if he could use my code in his package, and I said yes. It's very likely that he will see this and respond soon because I tagged him above, and he's a regular correspondent here. 

Wikipedia is good place to start reading about cluster analysis. My code used algorithms that I found on Wikipedia.

That's interesting. My test shows that using Newton's method is much faster and much more accurate than Statistics:-Quantile and works over the whole interval 0..1.  I wonder if this is a version difference. If you post another test, please show kernelopts(version), Digits, and a randomize seed so that your results can be duplicated.
 

Download NormalInveseTest.mw

@mmcdara I'm surprised that there's been such a change in the numeric computation model. Anyway, here are some minor adjustments, which I hope will also work in your Maple. I also corrected the issue for Digits > 15.

restart:
Digits:= 20: #Changing this doesn't change the code; it's just for testing.

#Symbolic steps are only needed to create the numeric procedures, not
#needed to run them:
Ncdf:= unapply(int(exp(-x^2/2)/sqrt(2*Pi), x= -infinity..X), X);
Newton:= unapply(simplify(x-(Ncdf(x)-Q)/D(Ncdf)(x)), (x,Q));

#Purely numeric procedures:
NormInv_inner:= proc(Q, eps)
local z:= 0, z1;
   do  z1:= evalf(Newton(z,Q));  if abs(z1-z) < eps then return z1 fi;  z:= z1  od
end proc:
      
NormInv:= (Q::And(numeric, realcons, satisfies(Q-> 0 <= Q and Q <= 1)))->
   if Q > 1/2 then -thisproc(1-Q)
   elif Q = 0 then -infinity
   elif Digits <= evalhf(Digits) then evalhf(NormInv_inner(Q, DBL_EPSILON))
   else NormInv_inner(evalf(Q), 10.^(1-Digits))
   fi
:

#Testing code:
P := [seq(1 - 10.^(-k), k=1..Digits)]:
printf("           p                  z\n");
printf("-----------------------------------\n");         
for n from 1 to numelems(P) do
   printf("%2d  %0.*f  %2.5f\n", n, Digits, P[n], NormInv(P[n]))
end do;
           p                  z
-----------------------------------
 1  0.90000000000000000000  1.28155
 2  0.99000000000000000000  2.32635
 3  0.99900000000000000000  3.09023
 4  0.99990000000000000000  3.71902
 5  0.99999000000000000000  4.26489
 6  0.99999900000000000000  4.75342
 7  0.99999990000000000000  5.19934
 8  0.99999999000000000000  5.61200
 9  0.99999999900000000000  5.99781
10  0.99999999990000000000  6.36134
11  0.99999999999000000000  6.70602
12  0.99999999999900000000  7.03448
13  0.99999999999990000000  7.34880
14  0.99999999999999000000  7.65063
15  0.99999999999999900000  7.94135
16  0.99999999999999990000  8.22208
17  0.99999999999999999000  8.49377
18  0.99999999999999999900  8.75680
19  0.99999999999999999990  9.01253
20  0.99999999999999999999  9.24939

I haven't researched or time- or accuracy-tested them, but I find it hard to believe that any of these "magic coefficients" methods are generally better than this simple Newton's method, which only took 10 minutes to write, can be understood immediately, and works for any setting of Digits or in evalhf mode.

(*>>>------------------------------------------------------
| A procedure to compute quantiles of the standard normal |
| (Gaussian) distribution in sfloat or hfloat arithmetic  |
|---------------------------------------------------------|
| Author: Carl Love <carl.j.love@gmail.com> 15-Aug-2019   |
------------------------------------------------------<<<*)
restart:
Digits:= 15: #Changing this doesn't change the code; it's just for testing.

#Symbolic steps are only needed to create the numeric procedures, not
#needed to run them:
Ncdf:= unapply(int(exp(-x^2/2)/sqrt(2*Pi), x= -infinity..X), X);
Newton:= unapply(simplify(x-(Ncdf(x)-Q)/D(Ncdf)(x)), (x,Q));

#Purely numeric procedures:
NormInv_inner:= proc(Q)
local z:= 0, z1;
   do  z1:= Newton(z,Q);  if z1=z then return z1 fi;  z:= z1  od
end proc:
      
NormInv:= (Q::And(realcons, satisfies(Q-> 0 <= Q and Q <= 1)))->
   if Q > 1/2 then -thisproc(1-Q)
   elif Q=0 then -infinity
   elif Digits <= evalhf(Digits) then evalhf(NormInv_inner(Q))
   else NormInv_inner(evalf(Q))
   fi
:

I'll admit that the rational-approximation methods may be a little faster for some small values of Digits.

@syhue The following very short and simple paper explains how and why to divide out the 2s when computing the decryption exponent d for a two-prime RSA via the Chinese Remainder Theorem (CRT). In particular, I draw your attention to step 4 at the very top of page 258 and the explanatory paragraph a little below that begins "To apply [CRT] in step 4...." The paper is "Faster RSA Algorithm for Decryption Using Chinese Remainder Theorem"^[*1]. This idea can be extended fairly easily to your three-prime RSA. Note that this special CRT is only needed once---for computing d; when you're decrypting, you'll use the three primes themselves as the moduli for CRT.

Let me know if you need further help.

[*1] G.N. Shinde and H.S. Fadewar, (2008), "Faster RSA Algorithm for Decryption Using Chinese Remainder Theorem", ICCES, vol.5, no.4, pp.255-261

@acer Thanks, your first interpretation of my question was correct: how the box was obtained within Maple's GUI.

How did you get that box shown in your Question?

@syhue The pairwise GCDs of p1-1, p2-1, p3-1 are all 2, as your code shows. Since this isn't 1, they aren't pairwise relatively prime, which is required for chrem. Perhaps you're supposed to divide out the 2s before doing chrem?

Everything after the error message is nonsense due to the error, so there's not much point talking about it until the error is fixed.