Data commit

Task/Numerical-integration/ATS/numerical-integration.ats

@@ -0,0 +1,188 @@
+#include "share/atspre_staload.hats"
+
+%{^
+#include <math.h>
+%}
+
+typedef FILEstar = $extype"FILE *"
+extern castfn FILEref2star : FILEref -<> FILEstar
+
+(* This type declarations is for composite quadrature functions for
+   all the different g0float typekinds. The function must either prove
+   termination or mask the requirement. (All of ours will prove
+   termination.) The function to be integrated will not be passed as
+   an argument, but inlined via the template mechanism. (This design
+   is more general. It can easily be used to write a quadrature
+   function that takes the argument, but also can be used for faster
+   code that requires no function call.) *)
+typedef composite_quadrature (tk : tkind) =
+  (g0float tk, g0float tk, intGte 2) -<> g0float tk
+
+extern fn {tk : tkind}
+composite_quadrature$func : g0float tk -<> g0float tk
+
+extern fn {tk : tkind} left_rule : composite_quadrature tk
+extern fn {tk : tkind} right_rule : composite_quadrature tk
+extern fn {tk : tkind} midpoint_rule : composite_quadrature tk
+extern fn {tk : tkind} trapezium_rule : composite_quadrature tk
+extern fn {tk : tkind} simpson_rule : composite_quadrature tk
+
+extern fn {tk : tkind}
+_one_point_rule$init_x :
+  g0float tk -<> g0float tk
+
+fn {tk : tkind}
+_one_point_rule : composite_quadrature tk =
+  lam (a, b, n) =>
+    let
+      prval [n : int] EQINT () = eqint_make_gint n
+      macdef f = composite_quadrature$func
+      val h = (b - a) / g0i2f n
+      val x0 = _one_point_rule$init_x<tk> h
+      fun
+      loop {i   : nat | i <= n} .<n - i>.
+           (i   : int i,
+            sum : g0float tk) :<> g0float tk =
+        if i = n then
+          sum
+        else
+          loop (succ i, sum + f(x0 + (g0i2f i * h)))
+    in
+      loop (0, g0i2f 0) * h
+    end
+
+(* The left rule, for any floating point type. *)
+implement {tk}
+left_rule (a, b, n) =
+  let
+    implement _one_point_rule$init_x<tk> _ = a
+  in
+    _one_point_rule<tk> (a, b, n)
+  end
+
+(* The right rule, for any floating point type. *)
+implement {tk}
+right_rule (a, b, n) =
+  let
+    implement _one_point_rule$init_x<tk> h = a + h
+  in
+    _one_point_rule<tk> (a, b, n)
+  end
+
+(* The midpoint rule, for any floating point type. *)
+implement {tk}
+midpoint_rule (a, b, n) =
+  let
+    implement _one_point_rule$init_x<tk> h = a + (h / g0i2f 2)
+  in
+    _one_point_rule<tk> (a, b, n)
+  end
+
+implement {tk}
+trapezium_rule : composite_quadrature tk =
+  lam (a, b, n) =>
+    let
+      prval [n : int] EQINT () = eqint_make_gint n
+      macdef f = composite_quadrature$func
+      val h = (b - a) / g0i2f n
+      fun
+      loop {i   : pos | i <= n} .<n - i>.
+           (i   : int i,
+            sum : g0float tk) :<> g0float tk =
+        if i = n then
+          sum
+        else
+          loop (succ i, sum + f(a + (g0i2f i * h)))
+      val sum = loop (1, g0i2f 0)
+    in
+      ((f(a) + sum + sum + f(b)) * h) / g0i2f 2
+    end
+
+(* Simpson’s 1/3 rule, for any floating point type. *)
+implement {tk}
+simpson_rule : composite_quadrature tk =
+  lam (a, b, n) =>
+    let
+      (* I have noticed that the Simpson rule is a weighted average of
+         the trapezium and midpoint rules, which themselves evaluate
+         the function at different points. Therefore, the following
+         should be efficient and produce good results. *)
+      val estimate1 = trapezium_rule<tk> (a, b, n)
+      val estimate2 = midpoint_rule<tk> (a, b, n)
+    in
+      (estimate1 + estimate2 + estimate2) / (g0i2f 3)
+    end
+
+extern fn {tk : tkind}
+fprint_result$rule : composite_quadrature tk
+
+extern fn {tk : tkind}
+fprint_result (outf    : FILEref,
+               message : string,
+               a       : g0float tk,
+               b       : g0float tk,
+               n       : intGte 2,
+               nominal : g0float tk) : void
+
+implement
+fprint_result<dblknd> (outf, message, a, b, n, nominal) =
+  let
+    val integral = fprint_result$rule<dblknd> (a, b, n)
+  in
+    fprint! (outf, "  ", message, "  ");
+    ignoret ($extfcall (int, "fprintf", FILEref2star outf,
+                        "%18.15le", integral));
+    fprint! (outf, "   (nominal + ");
+    ignoret ($extfcall (int, "fprintf", FILEref2star outf,
+                        "% .6le", integral - nominal));
+    fprint! (outf, ")\n")
+  end
+
+fn {tk : tkind}
+fprint_rule_results (outf    : FILEref,
+                     a       : g0float tk,
+                     b       : g0float tk,
+                     n       : intGte 2,
+                     nominal : g0float tk) : void =
+  let
+    implement fprint_result$rule<tk> (a, b, n) = left_rule<tk> (a, b, n)
+    val () = fprint_result (outf, "left rule      ", a, b, n, nominal)
+    implement fprint_result$rule<tk> (a, b, n) = right_rule<tk> (a, b, n)
+    val () = fprint_result (outf, "right rule     ", a, b, n, nominal)
+    implement fprint_result$rule<tk> (a, b, n) = midpoint_rule<tk> (a, b, n)
+    val () = fprint_result (outf, "midpoint rule  ", a, b, n, nominal)
+    implement fprint_result$rule<tk> (a, b, n) = trapezium_rule<tk> (a, b, n)
+    val () = fprint_result (outf, "trapezium rule ", a, b, n, nominal)
+    implement fprint_result$rule<tk> (a, b, n) = simpson_rule<tk> (a, b, n)
+    val () = fprint_result (outf, "Simpson rule   ", a, b, n, nominal)
+  in
+  end
+
+implement
+main () =
+  let
+    val outf = stdout_ref
+
+    val () = fprint! (outf, "\nx³ in [0,1] with n = 100\n")
+    implement composite_quadrature$func<dblknd> x = x * x * x
+    val () = fprint_rule_results<dblknd> (outf, 0.0, 1.0, 100, 0.25)
+
+    val () = fprint! (outf, "\n1/x in [1,100] with n = 1000\n")
+    implement composite_quadrature$func<dblknd> x = g0i2f 1 / x
+    val () = fprint_rule_results<dblknd> (outf, 1.0, 100.0, 1000,
+                                          $extfcall (double, "log", 100.0))
+
+    val () = fprint! (outf, "\nx in [0,5000] with n = 5000000\n")
+    implement composite_quadrature$func<dblknd> x = x
+    val () = fprint_rule_results<dblknd> (outf, 0.0, 5000.0, 5000000,
+                                          12500000.0)
+
+    val () = fprint! (outf, "\nx in [0,6000] with n = 6000000\n")
+    implement composite_quadrature$func<dblknd> x = x
+    val () = fprint_rule_results<dblknd> (outf, 0.0, 6000.0, 6000000,
+                                          18000000.0)
+
+    val () = fprint! (outf, "\n")
+  in
+    0
+  end