Soluzione

Definiamo un vettore ortogonale al piano

> n := vector([a,b,c]);

Calcoliamo un punto P_0 del piano, risolvendo il sistema lineare e dando dei valori arbitrari alle variabili libere

> P_0 := linsolve(matrix(1,3,[a,b,c]),vector([-d]));

> P_0 := subs(_t[1]=0,_t[2]=0,evalm(P_0));

Determiniamo la proiezione del vettore v = P-P_0 (ottenuto traslando P di -P_0 ) sul sottospazio di dimensione 1 generato da n (l'ortogonale del piano)

> v := P-P_0:
v_1 := dotprod(v,n)*n/dotprod(n,n);

> evalm(v_1);

Per differenza, calcoliamo la proiezione del vettore v sul sottospazio parallelo al piano (ottenuto traslando il piano di -P_0 )

> v_2 := evalm(v - v_1);

Traslando il risultato di P_0 si ottiene il punto cercato

> evalm(v_2 + P_0);

> Q := evalm(P-v_1);