@Carl Love 

Isn't this overkill? Why not simply

`&||` := (x,y) -> 1/(1/x+1/y);

if we use it only as operator (not as function `&||`(a,b,...)).

@Carl Love 

M:= f(output((J:= Iterator:-Combination(6,4))), 0, [2,2]):
(in J do M+~1 od);

f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1, f(Vector[row](4, [0, 1, 2, 3]), 0, [2, 2]) + 1

@Carl Love 

Do you mean that in your code Alias is invoked only once? It creates a 2x2 Matrix.

(Actually I am a little afraid of what Alias does and try to avoid it).

@Carl Love 

That's interesting, but I find the resulting code unintuitive.
Probably most users would prefer:

J:=Iterator:-Combination(6,4):
seq(ArrayTools:-Alias(x[]+~1, [2,2]), x in J);

@emendes 

The denominator of R[n] is  4*2^(2^n).

Actually, using rsolve ==>  R[n] = -1/2*cos(arccos(3/4)*2^n)+1/2

Using this, R can of course be formed even for n >> 1000.

@JanBSDenmark 

If you want an explanation. In mathematics sqrt and arcsinh are multivalued functions (defined on C).
All CASes have problems in choosing the wanted branches for such functions.
E.g. with your notation, sol[2] is a solution for another branch:  - x*sqrt(4*x^2 + 1) + arcsinh(2*x) - 5 = 0.

@David Sycamore 

P calls PP (from my first answer), so, PP must be present and executed in order to use P.

Just execute (use the (!!!) icon) the attached:

PP.mw

@David Sycamore 

You do not have the procudure PP in the worksheet (and executed, of course).

@David Sycamore 

The procedure PP can be easily used for any such questions:

PPrime:=proc(p::prime,n::posint)  
evalb(nops([PP(n,p,[p],1)])=1)
end:

P := proc(n)
local p, P:=select(isprime, [seq(1..n)]);
select(PPrime, P, n)
end:

P(50)
     [2, 3, 5, 7, 11, 13, 19]

@David Sycamore 

Just replace

Sol[++nSol] := XXX

with

nSol:=nSol+1; Sol[nSol] := XXX

(in 2 places)

@Carl Love 

It works in 2018 (even if ++x etc are not documented).

@Rouben Rostamian  

There are 2 solutions: take the symmetry  (x,a) -> (-x, -a).
(For external tangency there are also 2 solutions).

@Kitonum 

A continuous function e.g. on R has a continuous antiderivative but a CAS does not necessarily find one.

E.g.  int(1/(2+cos(x)), x) .  Even for a human it's not trivial.

@Kitonum 

The antiderivative is computed correctly, but it is not continuous, which techically is not a bug:

restart;
F:=int(sqrt(1+((k*Pi)/l*cos((Pi*x)/l))^2), x) assuming k>0,l>0:
F22:=eval(F,[l = 2, k = 2]):
limit(F22,x=2,left) - limit(F22,x=0, right); 
evalf(%);

      4*EllipticE(Pi*I)/Pi
      4.609785320 - 0.*I

@mehdi jafari 

I tried DirectSearch instead of fsolve.
It seems that the system has no solutions, see the attached worksheet.

Inflation-Inverse-1-vv.mw