Fix pow(x, 1/2) -> sqrt(x) miscompilation for integer division

When both operands in the exponent division are UInt, Stan performs
integer division (1/2 = 0), so pow(x, 1/2) is actually pow(x, 0) = 1.
The partial evaluator was incorrectly rewriting this to sqrt(x).

Add a type check to the guard so the simplification only fires when
at least one operand is real (i.e., actual floating-point division).
Also removes the comment that was flagging this as a known bug.

Add test case for pow(theta, 1/2) with pure integer literals.

Author: ishaan-arora-1

src/analysis_and_optimization/Partial_evaluator.ml

@@ -704,11 +704,11 @@ let rec eval_expr ?(preserve_stability = false) (e : Expr.Typed.t) =
                     ; { pattern=
                           FunApp (StanLib ("Divide__", FnPlain, mem), [y; z])
                       ; _ } ] )
-                  when is_int 1 y && is_int 2 z ->
+                  when is_int 1 y && is_int 2 z
+                       && not
+                            (y.meta.type_ = UInt && z.meta.type_ = UInt) ->
                     let lub_mem = lub_mem_pat [mem] in
                     FunApp (StanLib ("sqrt", suffix, lub_mem), [x])
-                    (* This is wrong; if both are type UInt the exponent is
-                       rounds down to zero. *)
                 | ( "square"
                   , [{pattern= FunApp (StanLib ("sd", FnPlain, mem), [x]); _}] )
                   ->

test/unit/Optimize.ml

@@ -2171,6 +2171,7 @@ model {
     target += pow(theta, 1./2.);
     target += pow(theta, 1/2.);
     target += pow(theta, 1./2);
+    target += pow(theta, 1/2);
     target += square(sd(x_vector));
     target += sqrt(2);
     target += sum(square(x_vector - y_vector));
@@ -2338,6 +2339,7 @@ model {
           target += sqrt(34.);
           target += sqrt(34.);
           target += sqrt(34.);
+          target += pow(34., 0);
           target += variance(x_vector);
           target += sqrt2();
           target += squared_distance(x_vector, y_vector);