@Markiyan Hirnyk 

You took another definition for the equidistance.

If you think that it is better, consider in your example (and notice also that implicitplot did not find all the points):

plots:-implicitplot(['F(s, t)[1] = 2', t = sin(s)], s = -Pi .. 3*Pi, t = -4 .. 4, color = [red, blue], gridrefine = 2)

If you use the first definition:

dist = 1/10 ==>

dist = 1 ==>

Why not, according to the definition considered: each point on the surface is moved by 1/2 units in the direction of the normal vector. I did not check whether the computations made by Maple are correct, have you?

@Markiyan Hirnyk 

And what is not so simple here?

@Markiyan Hirnyk

@one man

It would be nice to say something about the method used in the worksheet.
What "universal parametrization" have you used? Otherwise, the reader must use a "reverse engineering"  technique to find out.

If the surface is already parametrized, one may use something like this:

CP:=(v,w)->[v[2]*w[3]-v[3]*w[2], v[3]*w[1]-v[1]*w[3], v[1]*w[2]-v[2]*w[1]]: #Cross Product
IP:=(v,w)->v[1]*w[1]+v[2]*w[2]+v[3]*w[3]: # Inner Product
UV:=v -> v/~sqrt( IP(v,v)); #Unit Vector

r:=[s*cos(t), s*sin(t), s^2]:          # Example, paraboloid (parametric)
Nr:=UV( CP( diff(r,s), diff(r,t) ) ): #Unit Normal vector to r
Er:=simplify(  r + 1/2*Nr) ;           # Equidistant to r


p1:=plot3d( r, s=1 .. 2, t=0..2*Pi ):  p2:=plot3d( Er, s=1 .. 2, t=0..2*Pi ):
plots[display](p1,p2);






Edit. OP's animation using this method:

The problem is that Maple gives the result without assumptions, even if for a complex k, the integral may not exist (converge) e.g. for k = I.

@m3rcur1al 

You may use the Explore command and play with the parameters, as in the example:

Explore(fsolve(c1*exp(x)-c2*sin(x), x), parameters=[c1=1..15, c2=-3..13]);

You may prefer to replace fsolve with a plot.

@abcd 

You are right. Both our approaches are correct.
The reasons I suspected that you are wrong are:
1. I am a 'pure' mathematician. I am not used to see  r + dr  etc. For me it's easier to compute a Jacobian
(see a previous intervention).
2. The Triangular distribution was misleading. I would have used the (equivalent in this case) Power distribution.
3. I had the visual impression that the density is higher off the center.

All the best, and sorry for my mistake.

What are  f^n(x_0) and df^n(x_0)/dx ?

@Carl Love 

Probably a good idea (from an esthetical point of view) would be
`convert/If`:= ...

convert(... , 'If');

PS. Is such functional style considerably more efficient?

@Carl Love 

piecewise(x<1,10,x<2,20,30);

is equivalent to

`if`(x<1,10,`if`(x<2,20,30));

@Carl Love 

convert/`if`  seems to be not documened.
What type of expression does it obtain?
P:=piecewise(x<7,x+1,x);
PP:=convert(P,`if`);
lprint(PP);   #     Error, invalid expression for eval

Edit. OK, you answered the question while I was asking, thank you  :-),

@abcd 

A random variable X (with values in T subset R^2) is considered uniformly distributed in the sense
Prob(X in D) = Area(D)/Area(T)
for each measurable subset D of T.

My sample generated for the ellipse satisfies this, while yours - for the unit disc (obtained from uniform and triangular distributions) does not.
You can verify this by counting the number of points in a disc (with radius 1/2 e.g.) not centered at the origin.
[You did not apply the definition, but a kind of statistical intuition].

@Markiyan Hirnyk 

Nevermind, it was pure curiosity.

@Markiyan Hirnyk 

Yes, but what about the Mathematica timing?

@Markiyan Hirnyk 

We do not know how Mathematica measures the timing.
So, I thing that the best way to compare would be to execute in a loop the procedure (maybe with distind input data) and compute the difference of cpu or real time (or even use a physical clock).
Also, the  graphical part should be excluded if we want to compare only the algorithms.

I also think that Mathematica will win only if the random generator is better, because the algorithm is too simple and taylored exactly for the problem.