@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.

@Carl Love Here is a Maple worksheet with all the computations needed, and the final plots, as I outlined above. The plots match those from the paper. You should put a proper citation for the paper in here, and fill in the descriptions and units of the parameters. The plots could be dressed up a bit, but they clearly show that the computations are accurate.
 

Download Chemoimmunotherapy.mw

@tomleslie s(t) is the input function (or forcing function) of immunotherapy. Typically, it's a square wave, but for the purpose of generating plots (a) and (b), it's constant with value s_infinity. (Note: s0 <> s(0), unless by some extraordinary coincidence.) Likewise for q(t), which is the forcing function of chemotherapy.

@acer Vote up. So, I was wrong: It is possible (to some extent)^[*1] to control the order of multiplicands within a `*` term. I did know that it was possible to use sort to control the order of terms in a `+` expression, but I didn't know that it also worked within the `*`s. Please let me know if I was wrong about anything else.

Would you please show some dismantle examples, as you hinted at, with some exposition? To start with, the distinction between SUM (aka `+`) and POLY seems capricious: 3*a*b*c + a*b*d is a POLY and 3/4*a*b*c + a*b*d is a SUM. Also, I am worried about increasing divergence between the type system and the low-level storage. And by your aside about dismantle above, did you mean that certain data structures in the kernel exist primarily to affect the GUI display at higher levels of prettyprint (or perhaps typesetting)? 

[*1] I say "to some extent" because I still don't think that it's possible to put numeric coefficients any place other than the beginning or to have two expressions with equal `*` terms in different orders.

I'm just restating your Question for those who might find your rambling style confusing:

Can DocumentTools:-Components:-VideoPlayer (or any similar command) be used to stream live video?

Your plots (a) and (b) have their vertical axes labelled D[chemo]. How is D[chemo] determined in terms of the quantities that you've already defined, such as Q(t)?

@Annonymouse Apparently, axis[1]= [mode= log] doesn't work very well when there are bars whose lengths---before taking the logarithm---are so short that they wouldn't appear on a non-log plot. So, we need to take the logs before plotting and then rescale the horizontal axis after. So, here's a better, hand-rolled logarithmic bar chart. Use this after you've assigned the list difference.

difs:= remove(x-> rhs(x)=0, (lhs=abs@rhs)~(difference));
n:= nops(difs);
M:= floor(min((log[10]@rhs)~(difs)));

plot(
   [seq([[0,k], [log[10](rhs(difs[k]))-M, k]], k= 1..n)],
   axis[1]= [tickmarks= [seq(k= typeset(10^nprintf("%d", k+M)), k= 1..-M, 2)], gridlines],
   axis[2]= [tickmarks= [seq(k= typeset(lhs(difs[k])), k= 1..n)]],
   thickness= 9, axes= boxed, view= [default, 0..n+1]
);

Let me know if you need any adjustments.

@Annonymouse As you've no doubt figured out on your own (since apparently you found my Answer useful), axis[2] refers to the vertical axis and axis[1] refers to the horizontal axis. I'm sorry that I didn't catch at first that you were using Statistics:-BarChart, which reverses (IMO) the usual roles of the two axes, so my Answer should've used axis[1].

Still, I'm not satisfied with the results. I'll come up with something that uses regular plot instead.