Another step done: I’ve recoded the poincaré disk algorithm from scratch, so I can verify I’ve indeed understood the different geometric steps behind the poincaré disk.
Now, contrary to the original code of B. Horne, I’m also able to generate poincaré disks where the smallest primitive is not the (hyperbolic) triangle, but the whole (hyperbolic) polygon.
The latter will be better to have space to map image data into it.
Also, now I’ve understood the concept, the G’MIC code is also shorter 
Code:
foo :
1400,1400,1,1,":
begin(
const p = 5;
const q = 4;
const nb_iter = 20;
const is_polygon = 1;
const pi2 = 2*pi;
const pip = pi/p;
const pi2p = 2*pip;
const cosq2 = cos(pi/q)^2;
const sinp2 = sin(pip)^2;
const delta = cosq2 - sinp2;
const r0 = sqrt(sinp2/delta);
const x0 = sqrt(cosq2/delta);
const y0 = 0;
);
# Cartesian coordinates (x,y) and polar coordinates (r,t) (assumed to be always synchronized).
x = lerp(-1.1,1.1,x/(w-1));
y = lerp(-1.1,1.1,y/(h-1));
r = norm(x,y);
t = (atan2(y,x) + pi2)%pi2;
r>1?0:(
parity = 0;
repeat ((is_polygon?1:3)*nb_iter,k,
# Test if point (x,y) is in fundamental region.
is_polygon?(
trm = (round(t - pip,pi2p) + pip)%pi2;
P0 = rot(trm)*[ x0,y0 ];
is_inside_region = norm([x,y] - P0)>=r0;
):(
P0 = [ x0,y0 ];
is_inside_region = t<pip && norm(x - x0, y - y0)>=r0;
);
is_inside_region?break();
# Otherwise, transform (x,y) coordinates.
++parity;
is_polygon || (k3 = k%3)==2?(
u = x - P0[0];
v = y - P0[1];
nr = norm(u,v);
x = P0[0] + u*r0^2/nr^2;
y = P0[1] + v*r0^2/nr^2;
r = norm(x,y);
t = atan2(y,x)%pi2):
k3==1?(
y*=-1;
t = -t + pi2;
):
(
t = (((t + pip)%pi2p) - pip)%pi2;
x = r*cos(t); y = r*sin(t)
);
);
parity%3;
)"
r2dx 50%
Next step: map image pixels !
