Significant Digits

Purpose of exercise

• Practice writing equation-based programs
• Use engineering analysis to deal with numerical approximations
• Practice specifying correctness properties
• Use DrACuLa DoubleCheck for predicate-based, random testing
• Certify correctness properties using ACL2

Activities

• Define the following functions, according to the engineering analysis outlined below
• Package the functions you define as a book.
• DoubleCheck the properties stated below
• Use ACL2 to prove the properties
• Package the theorems as a book.

Requirements

• Mathematical details: StudentstEngineeringAnalysis.pdf

1. Define a function (log2-ceiling x) that delivers the least integer n such that x <= 2^n when x > 1.
   Note: The theorems in the arithmetic-3/floor-mod/floor-mod book will help admit log2-ceiling to the ACL2 logic.
2. Define a function (sqrt~ n x) that delivers an approximation to the square root of the rational number x >= 0 that is correct to n significant digits.
   Hints:
   (1) Range reduction: for non-zero values of x, find y and m such that y = x/(4^m) such that 1/4 < y <= 1
   (2) If y = x/(4^m), then sqrt(x) = sqrt(y)*(2^m)
   (3) Linear interpolation of sqrt(x) between x = 1/4 and x = 1 yields an approximation accurate to one digit.
   (4) Assume: the number of significant digits in the approximation doubles with each Newton iteration.
3. Use Machin's identity combined with Leibniz's arctangent series to define a function (pi~ n) that delivers an approximation to pi that is correct to n significant digits. Extract information from StudentstEngineeringAnalysis.pdf to justify your the accuracy of your approximation, put it in a document, and refer to the document at an appropriate point in your code's commentary. Construct the definitions so that the number of arithmetic operations required to carry out the computation is proportional to the number of terms in the series.
4. Use the Simpson's rule algorithm with automatic step-size halving to "convergence" and a minimal number of evaluations of the integrand to define a function that delivers approximations to Student's t distribution:
      (Students-t-dist~ n d x) = P[X < x], accurate to n significant digits
   where P[X < x] is the probability that a sample from a random variable X is less than x, where X conforms to Student's t distribution with d degrees of freedom.
5. Use Newton's method (the general one, not the sqrt one) to define a function that inverts Student's t distribution, accurate to n significant digits on the range 0.5 < p < 0.95:
     (Backwards-Students-t-dist~ n d p) = x, such that p = (Students-t-dist~ n d x)

Properties

• 1/2 < x/(2^(log2-ceiling x)) <= 1 when x > 1.
  Note: You will need definitions of log2-ceiling and 2^k that are easy to reason about. And, you will need some lemmas.

Engineering analysis

1. Arithmetic and significant digits. If x any y are approximations accurate to n significant digits, then so are x*y and x/y. Also, x + y is accurate to n significant digits if x and y are and both have the same sign.
2. Significant digits and relative error. An approximation with relative error is 10^(-n) is accurate to n significant digits.
3. Student's t symmetries. If s is a Student's t probability density function, then s(-x) = s(x). Therefore, for any positive value of x
       integral{s(u) | u < x} = 1/2 + integral{s(u) | 0 < u < x}, and
       integral{s(u) | u < -x} = 1/2 - integral{s(u) | 0 < u < x}
4. Simpson's error. Assume that the absolute error  an approximation given by Simpson's rule for a given step size, when integrating Student's t density, does not exceed twice the difference between that approximation and the one with twice its step size.
5. Pulling out the constant. When integrating a function of the form c*f(x), pull the constant c outside the integral. That is, compute the integral of f, then multiple by c. For example, if the variable of integration is x, then the formula Gamma((n+1)/2)/(sqrt(n*pi)*Gamma(n/2) is constant with respect to the variable of integration and can be pulled outside the integral.
6. Ratios of gamma values, not gamma values. To compute values of the density function for the Student's t distribution, you will need to compute ratios of Gamma-function values. You could do this by computing Gamma-function values, then dividing, but you can save some computation by planning ahead. When the number of degrees of freedom is even, the sqrt(pi) factors in the Gamma ratio will cancel. When the number of degrees of freedom is odd, a product of two sqrt(pi) factors emerges. So, you can go straight to pi and avoid a sqrt computation and a multiplication.
7. Student's t when x is near zero. When x is not near zero, you can compute the integral to n significant digits, then add 1/2, an exact value, without loosing any significant digits. When x is near zero, the integral is also near zero, and computing it to fixed, relative error is problematic, since the denominator in the relative error is the exact value of the integral, which is near zero. So, when x is near zero, you need to change your strategy. 
      The density function, s, that you are integrating is decreasing on the {u | 0 < u < x} and is convex near zero. So, if x is near zero, then s(x) < s(u) < s(0) < 1 whenever 0 < u < x. Therefore, the integral over that range is between x*s(x) and x*s(0). With a little algebraic manipulation, you can can conclude that the number of significant digits in the approximation to the integral will be at least log10(1/x), where log10 is the base-ten logarithm. So, if you want to make sure to have n significant digits of accuracy in the integral when x is near zero, you can just use the value 1/2 + x*s(0) in place of the integral when 1/x > 10^n.
      If you decide to use this fact (at least, I think it's a fact ... you need to confirm it), submit a document that justifies it, including the algebraic manipulations, and refer to the file at an appropriate point in the commentary of your program. If you decide not to use this fact, you will need to figure out another way to deal with computing an approximation accurate to n significant digits.
8. Initial approximation to start Newton's method to compute a value of the inverted Student's t distribution
   Some engineering analysis is needed to get an initial approximation with a known accuracy. Student's t distribution function converges toward the standard Gaussian distribution as the number of degrees of freedom increases. This might provide a way to guess a starting value for the iteration of Newton's method. One way to do this is to use a Gaussian statistics table to find the x-values that yield, say p1=0.5 (well, you don't have to look that one up ... it's x1=0), p2=0.8, and p3=0.9. Find the piecewise linear function through these points (p1, x1), (p2, x2), and (p3, x3). 
      Then, check to see of that curve approximates the inverse standard Guassian distribution to one significant digit of accuracy. You don't have to prove this by finding the place where there derivative is zero and so on. It's a smooth function. Just write a little program that displays the some closely space (p, x) values on the parabolic curve, and check the Guassian table to see if they're accurate to one digit. If they are, use that curve to generate your initial approximations. If not, adjust the p1, p2, and p3 points a little bit, and try again. If you find a curve that gets within 20% instead of 10% relative error, just go an extra iteration in Newton's method. If it's 50%, go two extra iterations. 
     Or maybe you should go a couple extra iterations anyway, to make up for the using the Gaussian distribution as an approximation to Student's t. (You could check this out by looking at the curves for Gaussian versus Student's t with various degrees of freedom, and adjust the number of iterations to make up for discrepancies.)