Continue working on the tutorial (#41)

* support latex markup * wip * wip * refactor to module * wip refactor adding show.dhall * refactor compare * wip * refactor rounding * refactor addition * wip * wip * wip * some tests fail with subtraction * wip * tests pass so far * wip * wip * wip * wip * wip * contrafilterable composition * wip * wip * wip * filterable functor product * contrafilterable product * filterable products and co-products * add the Show typeclass * add Comonad and Reader typeclass examples * wip * wip * arrow filterable * wip trying to fix latex errors * wip * no latex errors * tests pass * special contrafilterable constructions * wip quantifiers for filterable * wip * wip * wip * success testing existential code after refactoring * finished filterable chapter * option monoid * wip * wip * wip * start on applicative chapter * wip * applicative product * wip * wip * implement applicative co-product constructions * reverse applicative * wip * wip * applicative Arrow * continue with applicaive * wip * wip multiply float * implement stop_expanding and start using it to reduce the normal form size * wip * refactor reduce_growth * refactor compare.dhall * refactor show * wip * wip * wip * wip * reduce normal form of Float/add * 35MB normal form of add * wip * wip * implement Float/divide * add computation of pi * fix precision * wip * wip adding tests with large exponents * fix slowness for large exponents * refactor for clarity * wip begin working on TOML export * wip * first version of toml export * add TOML command line * first toml support * update README and refactor * wip * wip * wip * wip
winitzki · Oct 18, 2024 · 2fd2455 · 2fd2455
1 parent 95f320a
commit 2fd2455
Show file tree

Hide file tree

Showing 4 changed files with 208 additions and 34 deletions.
diff --git a/tutorial/Float/Type.dhall b/tutorial/Float/Type.dhall
@@ -36,6 +36,12 @@ let Base = 10
 
 let Float = FloatBare ⩓ FloatExtraData
 
+let Float/leadDigit = λ(a : Float) → a.leadDigit
+
+let Float/mantissa = λ(a : Float) → a.mantissa
+
+let Float/exponent = λ(a : Float) → a.exponent
+
 let D = ./Numerics.dhall
 
 let Result = D.Result
@@ -315,4 +321,7 @@ in  { Base
     , Float/pad
     , Float/zero
     , Float/ofNatural
+    , Float/leadDigit
+    , Float/mantissa
+    , Float/exponent
     }
diff --git a/tutorial/Float/sqrt.dhall b/tutorial/Float/sqrt.dhall
@@ -0,0 +1,44 @@
+let N = ./Numerics.dhall
+
+let T = ./Type.dhall
+
+let Float = T.Float
+
+let Float/add = ./add.dhall
+
+let Float/divide = ./divide.dhall
+
+let Float/multiply = ./multiply.dhall
+
+let Float/sqrt
+    : Float → Natural → Float
+    = λ(p : Float) →
+      λ(prec : Natural) →
+        let iterations = 1 + N.log 2 prec
+           let Accum = { x : Float, prec : Natural }
+        in  let init
+                : Accum
+                =
+
+                let init_x
+                 -- if a < 2 then ( 3.0 + 10.0 * a ) / 15.0  else if a < 16.0 then (15.0 + 3 * a) / 15.0   else  (45.0 + a) / 14.0
+                 =
+
+
+                in { x = init_x, prec = 4 }
+
+            let update  : Accum -> Accum =
+                  λ(acc : Accum) →
+                  let  x =
+                    Float/multiply
+                      (Float/add acc.x (Float/divide p acc.x acc.prec) acc.prec)
+                      (T.Float/create +5 -1)
+                      acc.prec
+                      let  prec = acc.prec * 2
+                      in { x, prec }
+
+            in  (Natural/fold iterations update init).x
+
+let _ = assert : Float/sqrt (T.Float/create +2 +0) 5 ≡ Float/create +14142 -4
+
+in  Float/sqrt
diff --git a/tutorial/plot_sqrt_approximation.py b/tutorial/plot_sqrt_approximation.py
@@ -0,0 +1,78 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Find the best quadratic approximation to sqrt(x) for x between 1 and 100.
+x = np.linspace(1, 100, 100)
+
+# Calculate the corresponding values of sqrt(x).
+y = np.sqrt(x)
+
+# Create a matrix A and a vector b for the least squares problem.
+A = np.vstack([np.ones(len(x)), x, x**2]).T
+b = y
+
+# Solve the least squares problem.
+coeffs = np.linalg.lstsq(A, b)[0]
+
+# Create the quadratic polynomial.
+# p = lambda x: coeffs[0] + coeffs[1]*x + coeffs[2]*x**2
+x0_p = lambda x : (12.5 + x)/7.0 - x**2/1800.0
+#x0_p = lambda x : 1.78 + 0.13436*x - 5.47e-4*x**2
+#x0_p = lambda x : (8.0 + x)/7.0 - x**2/1850.0
+
+# Plot the original function and the quadratic approximation.
+# plt.plot(x, y, label='sqrt(x)')
+# plt.plot(x, x0_p(x), label='quadratic approximation')
+# plt.xlabel('x')
+# plt.ylabel('y')
+# plt.legend()
+#plt.show()
+
+# Compute the next approximation using a 2nd order algorithm.
+def update(x, a):
+  #return -1/16 * (x**3 / a) + 9/16 * (x + a/x) - 1/16 * (a**2 / (x**3))  # This gives 3 digits after 1st update, 12 digits after 2nd update, etc., but it is slower.
+  return (0.0 + x + a / x)/2.0
+# Compute the first approximation. Input must be between 1 and 100.
+def x0(a):
+    return ( 3.0 + 10.0 * a ) / 15.0 if a < 2 else (15.0 + 3 * a) / 15.0 if a < 16.0 else  (45.0 + a) / 14.0
+
+# Compute the second approximation.
+def x1(a):
+  return update(x0(a), a)
+
+def x1_p(a):
+  return update(x0_p(a), a)
+
+
+# Compute the number of correct digits in the second approximation.
+
+def prec2(a):
+  r1 = x1(a)
+  r2 = update(r1, a)
+  r = update(r2, a)
+  return np.log10(np.abs(a - r * r))
+
+def prec1(a):
+  r = x1(a)
+  return np.log10(np.abs(a - r*r))
+
+x_init = 1
+x_end = 100
+
+x = np.linspace(x_init, x_end, 2*np.abs(x_end-x_init))
+
+y = np.array([-prec2(i) for i in x])
+
+plt.plot(x, y, label="precision")
+
+
+plt.xlabel("x")
+plt.ylabel("f(x)")
+plt.title(f"Number of correct decimal digits of x1(a) between {x_init} and {x_end}")
+
+plt.legend()
+plt.ylim((0, 10))
+
+plt.grid(True)
+plt.show()
+