wip

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
     }

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

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()
+

tutorial/programming_dhall.md

@@ -1692,10 +1692,34 @@ let gcd : Natural → Natural → Natural = λ(x : Natural) → λ(y : Natural)
 
 The built-in Dhall type `Double` does not support any numerical operations.
 However, one can use values of type `Natural` to implement floating-point arithmetic.
-The `scall` repository contains [some proof-of-concept code](https://github.com/winitzki/scall/blob/master/tutorial/Float/) that implements a number of floating-point operations in arbitrary precision.
+The `scall` repository contains [some proof-of-concept code](https://github.com/winitzki/scall/blob/master/tutorial/Float/) that implements a number of floating-point operations: `Float/create`, `Float/show`, `Float/compare`, `Float/add`, `Float/subtract`, `Float/multiply`, `Float/divide` and so on.
+Floating-point numbers are represented by a decimal mantissa and a decimal exponent, and may have arbitrarily high precision (in both mantissa and exponent).
 
 An example of an arbitrary-precision numerical algorithm is the computation of a floating-point square root.
 
+We will use the following algorithm that computes successive approximations for $x = \sqrt p$, where $p$ is a given non-negative number:
+
+1. Compute the initial approximation $x_0$ that is close to $\sqrt p$.
+2. Estimate the total number of iterations $n$, where $n \ge 1$.
+3. Apply $n$ times the function `update`. 
+
+The result is the Dhall code `Natural/fold n update x0`.
+
+The initial approximation is defined as follows:
+
+1. Find the largest integer number $k$ such that $p = q * 10^{2k}$ and $q \ge 1$. Then we will have $1 \le q \lt 100$.
+2. If $q \lt 2$ then the initial value is $x0 = (3 + 10 * q) / 15$. If $2 \le q \lt 16$ then $x0 = (15 + 3 * q) / 15$. If $q 16 \le q \lt 100$ then $x0 = (45 + q) / 14$. The divisions here may be performed in very low precision (2-3 digits).
+3. The update function is computed as $u(x) = \frac{1}{2}(x+p/x) $.
+4. The number of correct decimal digits doubles after each update. The total number of iterations is estimated as $n = 1 + \log_2 N$. The first iteration gives 2 correct digits, the second 4 digits, the third 8 digits, etc.
+
+```dhall
+let Float/sqrt = λ(p : Float) → λ(prec : Natural) →
+  let iterations = 1 + (./Numerics.dhall).log 2 prec
+  let init : Float = ??? -- Code omitted for brevity.
+  let update = λ(x : Float) → Float/multiply (Float/add x (Float/divide p x prec) prec) (T.Float/create +5 -1) prec
+  in Natural/fold iterations update init
+```
+
 TODO
 
 ## Programming with functions
@@ -4133,7 +4157,7 @@ When we apply `example2` to arguments `r`, `leaf`, and `branch`, the code of `ex
 
 This explains how the Church encoding replaces iterative computations by non-recursive code.
 A data structure that contains, say, 1000 data values is Church-encoded into a certain higher-order function.
-That function will be hard-coded to call its arguments 1000 times.
+That function is hard-coded to call its arguments 1000 times.
 
 In this way, it is guaranteed that all recursive structures will be finite and all operations on those structures will terminate.
 That's why Dhall is able to accept Church encodings of recursive types and perform iterative and recursive operations on Church-encoded data without compromising any safety guarantees.