Faa di Bruno's formula.

Here's how it's supposed to go.

> p:=2038074743;
  q:=1190494759;
  n:=p*q;
  x:=1414213562373095048;
  y:=Power(x,2) mod n;

y := 363278601055491705

> yp:= Power(y,(p+1)/4) mod p;

yp := 1701899961

> yq:= Power(y,(q+1)/4) mod q;

yq := 325656728

> b:= chrem([yp,yq],[p,q]); 
  -b mod n;

                       b := 1488466947368437834

                          937850352623334103
 

> a:= chrem([yp,-yq],[p,q]);
  -a mod n;

                       a := 1414213562373095048

                         1012103737618676889
 

Here's how it's supposed to go.

> p:=2038074743;
  q:=1190494759;
  n:=p*q;
  x:=1414213562373095048;
  y:=Power(x,2) mod n;

y := 363278601055491705

> yp:= Power(y,(p+1)/4) mod p;

yp := 1701899961

> yq:= Power(y,(q+1)/4) mod q;

yq := 325656728

> b:= chrem([yp,yq],[p,q]); 
  -b mod n;

                       b := 1488466947368437834

                          937850352623334103
 

> a:= chrem([yp,-yq],[p,q]);
  -a mod n;

                       a := 1414213562373095048

                         1012103737618676889
 

> sum(binomial(n,k)*2^(m*k), k=0..n);

(2^m+1)^n

What that gives me when the Array is displayed is two rectangles of the same size; the colour bar is narrow, but there's a lot of white space around it.

What that gives me when the Array is displayed is two rectangles of the same size; the colour bar is narrow, but there's a lot of white space around it.

Trigonometric equations are tricky to deal with.  It took some experimenting, but the following works.

> eq:= sin(3*x) - cos(5*x);
   expand(convert(eq, exp));

-1/2*I*exp(x*I)^3+1/2*I/exp(x*I)^3-1/2*exp(x*I)^5-1/2*1/(exp(x*I)^5)

> factor(subs(exp(x*I)=z, %));

1/2*(-z^8+I)*(z^2+I)/z^5

> solve(subs(z=exp(I*x),%),x,AllSolutions);

1/16*Pi*(1+4*_Z1), 1/4*Pi*(-1+4*_Z2)
  

Trigonometric equations are tricky to deal with.  It took some experimenting, but the following works.

> eq:= sin(3*x) - cos(5*x);
   expand(convert(eq, exp));

-1/2*I*exp(x*I)^3+1/2*I/exp(x*I)^3-1/2*exp(x*I)^5-1/2*1/(exp(x*I)^5)

> factor(subs(exp(x*I)=z, %));

1/2*(-z^8+I)*(z^2+I)/z^5

> solve(subs(z=exp(I*x),%),x,AllSolutions);

1/16*Pi*(1+4*_Z1), 1/4*Pi*(-1+4*_Z2)
  

So your real problem is to define the function f in such a way that f(x,y) will have a (real) numerical value when x and y are numbers.   It might help to use something like

 evalf(sum(alpha1(r)*exp(x*r), r = RootOf(_Z^3 + x*_Z + y)))

Note that the exponential function is exp(x), not exp^(x).

So your real problem is to define the function f in such a way that f(x,y) will have a (real) numerical value when x and y are numbers.   It might help to use something like

 evalf(sum(alpha1(r)*exp(x*r), r = RootOf(_Z^3 + x*_Z + y)))

Note that the exponential function is exp(x), not exp^(x).

A solution using RootOf is likely to be as "clean" as you could possibly get. 

A solution using RootOf is likely to be as "clean" as you could possibly get. 

You could try Laplace transform.  If U(x,s) is the Laplace transform of T(x,t) with respect to t,
then your PDE becomes

s*U(x,s) - T(x,0) = alpha^2 * diff(U(x,s),x$2)

which you want to solve with boundary conditions U(0, s) = F(s) and U(1,s) = P(s) (the Laplace transforms of f(t) and p(t).  If T(x,0) = constant the boundary value problem can be solved in closed form.  However, there's no guarantee that the inverse Laplace transform can be done in closed form.

You could try Laplace transform.  If U(x,s) is the Laplace transform of T(x,t) with respect to t,
then your PDE becomes

s*U(x,s) - T(x,0) = alpha^2 * diff(U(x,s),x$2)

which you want to solve with boundary conditions U(0, s) = F(s) and U(1,s) = P(s) (the Laplace transforms of f(t) and p(t).  If T(x,0) = constant the boundary value problem can be solved in closed form.  However, there's no guarantee that the inverse Laplace transform can be done in closed form.

en.wikipedia.org/wiki/Fourier_transform