%------------------------------------------------------------------------------
% Reciprocal
%
%     Author: David Lester, Manchester University
%
%     Version 1.0            18/2/09   Initial Release Version
%------------------------------------------------------------------------------

inv: THEORY

BEGIN

  IMPORTING prelude_aux, appendix, cauchy, unique, neg, min_nat[nat]

  cx:   VAR cauchy_real
  x:    VAR real
  s,p:  VAR nat
  nzcx: VAR cauchy_nzreal
  nzx:  VAR nzreal
  pcx:  VAR cauchy_posreal
  px:   VAR posreal
  pn:   VAR posnat
  n,r:  VAR int
  nn:   VAR negint
  nzn:  VAR nzint
  nx:   VAR negreal

  expt_x2:     LEMMA x^2 = x*x
  expt_times2: LEMMA nzx^(2*s) = nzx^s*nzx^s

  minimum_inv_exists: LEMMA nonempty?({s | 3 <= abs(nzcx(s))})

  minimum_inv(nzcx):nat = min_nat.min({s | 3 <= abs(nzcx(s))})

  minimum_inv_prop0: LEMMA s = minimum_inv(nzcx) => 3 <= abs(nzcx(s))

  minimum_inv_prop1: LEMMA s = minimum_inv(nzcx) AND cauchy_prop(x,nzcx)
                     => 2 <= abs(x)*2^s
  minimum_inv_prop2:
     LEMMA 2 <= abs(nzx)*2^s AND n-1 < nzx*2^(p+2*s+2) AND nzx*2^(p+2*s+2)< n+1
           => n /= 0

  minimum_inv_aux(nzcx:cauchy_nzreal,s:{n:nat | n <=  minimum_inv(nzcx)}):
   RECURSIVE nat
    = IF 3 <= abs(nzcx(s)) THEN s ELSE minimum_inv_aux(nzcx,s+1) ENDIF
      MEASURE (LAMBDA (nzcx:cauchy_nzreal,s:{n:nat | n <=  minimum_inv(nzcx)}):
               IF 3 <= abs(nzcx(s)) THEN 0 
                              ELSE minimum_inv(nzcx)-s ENDIF)

  minimum_inv_impl(nzcx):nat = minimum_inv_aux(nzcx,0)

  minimum_inv_impl_def: LEMMA minimum_inv(nzcx) = minimum_inv_impl(nzcx)

  inv_p0: LEMMA 2^(p+s+3) -1 < pn => 6 <= pn
  inv_p1: LEMMA 2^(p+s+3) -1 < pn => 2^(2*p+2*s+3) < pn*(pn-1)
  inv_p2: LEMMA 2^(p+s+3) -1 < pn => 2^(2*p+2*s+3) < pn*(pn+1)
  inv_p3: LEMMA 2^(p+s+3) -1 < pn => 2^(2*p+2*s+2)/pn -(pn-1)/2 < 0
  inv_p4: LEMMA 2^(p+s+3) -1 < pn => 2^(2*p+2*s+2)/pn -(pn+1)/2 < 0
  inv_p5: LEMMA r = round(2^(2*p+2*s+2)/pn) AND 2^(p+s+3) -1 < pn
                                  => (r-1)*(pn+1) < 2^(2*p+2*s+2)
  inv_p6: LEMMA r = round(2^(2*p+2*s+2)/pn) AND 2^(p+s+3) -1 < pn
                                  => 2^(2*p+2*s+2) < (r+1)*(pn-1)
  inv_p7: LEMMA r = round(2^(2*p+2*s+2)/pn) AND 2^(p+s+3) -1 < pn
              => (r-1)/2^p < 2^(p+2*s+2)/(pn+1) AND
                      2^(p+2*s+2)/(pn-1) < (r+1)/2^p

  inv_p8: LEMMA 2 <= px*2^s AND pn-1 < px*2^(p+2*s+2) AND px*2^(p+2*s+2) < pn+1
              => 2^(p+s+3) -1 < pn

  inv_p9: LEMMA 2 <= px*2^s AND r = round(2^(2*p+2*s+2)/pn) AND
                pn-1 < px*2^(p+2*s+2) AND px*2^(p+2*s+2) < pn+1
              => (r-1)*px < 2^p AND 2^p < (r+1)*px

  inv_n5: LEMMA r = round(2^(2*p+2*s+2)/nn) AND 2^(p+s+3) -1 < -1*nn
                                  => (r-1) < 2^(2*p+2*s+2)/(nn+1)
  inv_n6: LEMMA r = round(2^(2*p+2*s+2)/nn) AND 2^(p+s+3) -1 < -1*nn
                                  => 2^(2*p+2*s+2)/(nn-1) < (r+1)
  inv_n7: LEMMA r = round(2^(2*p+2*s+2)/nn) AND 2^(p+s+3) -1 < -nn
              => (r-1)/2^p < 2^(p+2*s+2)/(nn+1) AND
                      2^(p+2*s+2)/(nn-1) < (r+1)/2^p

  inv_n8: LEMMA 2 <= -nx*2^s AND nn-1 < nx*2^(p+2*s+2) AND nx*2^(p+2*s+2) < nn+1
              => 2^(p+s+3) -1 < -nn

  inv_n9: LEMMA 2 <= -nx*2^s AND r = round(2^(2*p+2*s+2)/nn) AND
                nn-1 < nx*2^(p+2*s+2) AND nx*2^(p+2*s+2) < nn+1
              => (r-1) < 2^p/nx AND 2^p/nx < (r+1)

  inv_p: LEMMA 2 <= abs(nzx)*2^s AND r = round(2^(2*p+2*s+2)/nzn) AND
                nzn-1 < nzx*2^(p+2*s+2) AND nzx*2^(p+2*s+2) < nzn+1
              => (r-1) < 2^p/nzx AND 2^p/nzx < (r+1)


  minimum_inv_prop3: LEMMA s = minimum_inv(nzcx) AND cauchy_prop(nzx,nzcx) =>
          cauchy_prop(1/nzx,LAMBDA p: round(2^(2*p+2*s+2)/nzcx(p+2*s+2)))

  cauchy_nz_inv(nzcx:cauchy_nzreal):cauchy_nzreal
    = (LAMBDA p: LET s = minimum_inv_impl(nzcx)
                 IN round(2^(2*p+2*s+2)/nzcx(p+2*s+2)))

  inv_nz_lemma: LEMMA cauchy_prop(nzx,nzcx)
         => cauchy_prop(1/nzx, cauchy_nz_inv(nzcx))

  cauchy_inv(nzcx): cauchy_nzreal = cauchy_nz_inv(nzcx)

  inv(nzx:nonzero_real):nzreal = 1/nzx

  inv_lemma: LEMMA cauchy_prop(nzx,nzcx) => cauchy_prop(1/nzx,cauchy_inv(nzcx))

END inv