For integer exponents this is called Hermite interpolation.
One can also interpret this task as an application of the chinese remainder theorem in polynomial rings
f(x) == 0 mod x^p
f(x) == 1 mod (x-1)^p
The solution has the form
f(x) = a(x)*x^p+b(x)*(x-1)^p with deg a(x) < p, deg b(x) < p.
Inserting into the second equations gives
a(x)*x^p == 1 mod (x-1)^p <=> a(y+1)*(y+1)^p == 1 mod y^p
which can be solved via binomial series, i.e., direct power series arithmetics without solving linear systems of equations,
a(y+1) = 1 - p*y + p*(p+1)/2*y^2 - p*(p+1)*(p+2)/2*y^3 +-... ...*y^(p-1)
For non-integer exponents, why would you want to do that?
Also look up "sigmoid functions".
Update
If continuity of derivatives is not so much a concern then for any positive a
you can use c*x^a
for 0<=x<=0.5
and 1-c*(1-x)^a
for 0.5<=x<=1
. To close the gap at x=0.5
, the constant has to be chosen such that
1 = 2*c*0.5^a <=> c = 2^(a-1)
which can be implemented for 0<=x<=1
as
y = (0.5>x) ? 0.5*pow(2*x,a) : 1-0.5*pow(2*(1-x), a);