@Annonymouse 

subsindets(P, specindex(integer,x), _XX)

replaces each x[n]   having n::integer  by _XX(x[n]);
Then,
remove(has, %, _XX);
removes the polynomials containing _XX

@nm 

It seems that Maple ignores the singular points. You can simply exclude them from Carl's solution:

sing := remove(has,[singular(f)],infinity);

@Adam Ledger 

The function is discontinuous; you may add option discont.

Choosing these coords=...  appears pointless, unless you have a good reason for it.

AFAIK there is no user-defined storage functions for Arrays.  A possibility is to use storage=sparse and write the index function; it would be interesting to see the efficiency in this case.

@Markiyan Hirnyk 

It would be nice to tell us your ideas. E.g.
1. A multi-dimensional Van Aubel theorem.
2. A complex program for animated flexible polyhedra.
3. A project to generalize Gaifullin's results.
4. Something else?
 

@mmcdara 

0<x<1, 0<y<1 is not necessary because we are going to integrate the product of densities.

E.g. for a=4/3 Maple gives;

So, it's not unreasonable to expect from PDF (when a=4/3) the expression

piecewise(z<=0, 0, F(z))
but it seems that you are ok with the present status.

@mmcdara 

I expect from PDF to return at least the (double) integral (not necessarily a closed form, of course, even if such a closed form exists). This can be easily done for the toy example because Maple is able to solve Z < z:

SolveTools:-SemiAlgebraic({x^2 + y^2 + a*x*y<z},[x,y], 'parameters'=[a,z]);

(not to mention the case when a is fixed, e.g. a=1/2 or a=3).

So, I don't agree about "Not really a severe limitation".

@Alger 

It is not clear what you are trying to do. If you want to find an orthogonal system in a function space, the simplest way is to start with a known one e.g. cos(nx), sin(nx), in L²([0,2*Pi], 1)  and use the Fact above for some h.

@Rouben Rostamian  

The fact applied here is simple:

If u_n(x) are orthogonal in L²([x1,x2], w(x))  (here w is the weight) and h : [t1, t2] --> [x1,x2] is a diffeomorphism
then   u_n(h(t))  are orthogonal in L²([t1,t2], w(h(t))*|h'(t)|).

@minhthien2016 

Thank you for the new recurrence.
Have you tried to replace the new coefficients in the provided answer?

@Mariusz Iwaniuk 

On my systems (Maple 2017 or 2018, Win 7, 64 bit) it works in all situations: 1D/2D, Document/Worksheet.

@Carl Love 

FAIL/FAIL and even FAIL/0 or 0/FAIL returns FAIL, so the previous version seems to be ok (with the typo corrected).

If d=0 why have you decided to return 1? It should be FAIL (or remain 0/0, hence error) because in this case the definition of the conditional probability is different, so not handled by the formula.

@tomleslie 

It was just a remark. I agree with you of course that the results of the two parsers should be the same.

Using complex numbers the proof is even shorter:

P[5]:=P[1]:
for i to 4 do Q[i] := P[i] + (P[i+1]-P[i])*exp(-I*Pi/4)/sqrt(2) od:
is( Q[3]-Q[1] = I* (Q[2]-Q[4]) );

        true

n dots (n>=2) surrounded by spaces (which is probably the preferable syntax) are interpreted as  ..  in both 1D and 2D.
seq(-3 ... 3);  seq(-3 .... 3);