Merge pull request #1602 from ishaan-arora-1/fix/pow-int-division-sqrt-miscompilation

Fix pow(x, 1/2) miscompilation when exponent uses integer division

Author: WardBrian

src/analysis_and_optimization/Partial_evaluator.ml

@@ -704,11 +704,10 @@ 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

@@ -2257,6 +2257,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));
@@ -2424,6 +2425,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);