src/HOL/Library/Sum_of_Squares/sum_of_squares.ML
author wenzelm
Sun Nov 26 21:08:32 2017 +0100 (16 months ago)
changeset 67091 1393c2340eec
parent 63523 54e932f0c30a
child 67271 48ef58c6cf4c
permissions -rw-r--r--
more symbols;
     1 (*  Title:      HOL/Library/Sum_of_Squares/sum_of_squares.ML
     2     Author:     Amine Chaieb, University of Cambridge
     3     Author:     Philipp Meyer, TU Muenchen
     4 
     5 A tactic for proving nonlinear inequalities.
     6 *)
     7 
     8 signature SUM_OF_SQUARES =
     9 sig
    10   datatype proof_method = Certificate of RealArith.pss_tree | Prover of string -> string
    11   val sos_tac: (RealArith.pss_tree -> unit) -> proof_method -> Proof.context -> int -> tactic
    12   val trace: bool Config.T
    13   val debug: bool Config.T
    14   val trace_message: Proof.context -> (unit -> string) -> unit
    15   val debug_message: Proof.context -> (unit -> string) -> unit
    16   exception Failure of string;
    17 end
    18 
    19 structure Sum_of_Squares: SUM_OF_SQUARES =
    20 struct
    21 
    22 val max = Integer.max;
    23 
    24 val denominator_rat = Rat.dest #> snd #> Rat.of_int;
    25 
    26 fun int_of_rat a =
    27   (case Rat.dest a of
    28     (i, 1) => i
    29   | _ => error "int_of_rat: not an int");
    30 
    31 fun lcm_rat x y =
    32   Rat.of_int (Integer.lcm (int_of_rat x) (int_of_rat y));
    33 
    34 fun rat_pow r i =
    35  let fun pow r i =
    36    if i = 0 then @1 else
    37    let val d = pow r (i div 2)
    38    in d * d * (if i mod 2 = 0 then @1 else r)
    39    end
    40  in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end;
    41 
    42 fun round_rat r =
    43   let
    44     val (a,b) = Rat.dest (abs r)
    45     val d = a div b
    46     val s = if r < @0 then ~ o Rat.of_int else Rat.of_int
    47     val x2 = 2 * (a - (b * d))
    48   in s (if x2 >= b then d + 1 else d) end
    49 
    50 
    51 val trace = Attrib.setup_config_bool @{binding sos_trace} (K false);
    52 val debug = Attrib.setup_config_bool @{binding sos_debug} (K false);
    53 
    54 fun trace_message ctxt msg =
    55   if Config.get ctxt trace orelse Config.get ctxt debug then tracing (msg ()) else ();
    56 fun debug_message ctxt msg = if Config.get ctxt debug then tracing (msg ()) else ();
    57 
    58 exception Sanity;
    59 
    60 exception Unsolvable;
    61 
    62 exception Failure of string;
    63 
    64 datatype proof_method =
    65     Certificate of RealArith.pss_tree
    66   | Prover of (string -> string)
    67 
    68 (* Turn a rational into a decimal string with d sig digits.                  *)
    69 
    70 local
    71 
    72 fun normalize y =
    73   if abs y < @1/10 then normalize (@10 * y) - 1
    74   else if abs y >= @1 then normalize (y / @10) + 1
    75   else 0
    76 
    77 in
    78 
    79 fun decimalize d x =
    80   if x = @0 then "0.0"
    81   else
    82     let
    83       val y = abs x
    84       val e = normalize y
    85       val z = rat_pow @10 (~ e) * y + @1
    86       val k = int_of_rat (round_rat (rat_pow @10 d * z))
    87     in
    88       (if x < @0 then "-0." else "0.") ^
    89       implode (tl (raw_explode(string_of_int k))) ^
    90       (if e = 0 then "" else "e" ^ string_of_int e)
    91     end
    92 
    93 end;
    94 
    95 (* Iterations over numbers, and lists indexed by numbers.                    *)
    96 
    97 fun itern k l f a =
    98   (case l of
    99     [] => a
   100   | h::t => itern (k + 1) t f (f h k a));
   101 
   102 fun iter (m,n) f a =
   103   if n < m then a
   104   else iter (m + 1, n) f (f m a);
   105 
   106 (* The main types.                                                           *)
   107 
   108 type vector = int * Rat.rat FuncUtil.Intfunc.table;
   109 
   110 type matrix = (int * int) * Rat.rat FuncUtil.Intpairfunc.table;
   111 
   112 fun iszero (_, r) = r = @0;
   113 
   114 
   115 (* Vectors. Conventionally indexed 1..n.                                     *)
   116 
   117 fun vector_0 n = (n, FuncUtil.Intfunc.empty): vector;
   118 
   119 fun dim (v: vector) = fst v;
   120 
   121 fun vector_cmul c (v: vector) =
   122   let val n = dim v in
   123     if c = @0 then vector_0 n
   124     else (n,FuncUtil.Intfunc.map (fn _ => fn x => c * x) (snd v))
   125   end;
   126 
   127 fun vector_of_list l =
   128   let val n = length l in
   129     (n, fold_rev FuncUtil.Intfunc.update (1 upto n ~~ l) FuncUtil.Intfunc.empty): vector
   130   end;
   131 
   132 (* Matrices; again rows and columns indexed from 1.                          *)
   133 
   134 fun dimensions (m: matrix) = fst m;
   135 
   136 fun row k (m: matrix) : vector =
   137   let val (_, j) = dimensions m in
   138     (j,
   139       FuncUtil.Intpairfunc.fold (fn ((i, j), c) => fn a =>
   140         if i = k then FuncUtil.Intfunc.update (j, c) a else a) (snd m) FuncUtil.Intfunc.empty)
   141   end;
   142 
   143 (* Monomials.                                                                *)
   144 
   145 fun monomial_eval assig m =
   146   FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => a * rat_pow (FuncUtil.Ctermfunc.apply assig x) k)
   147     m @1;
   148 
   149 val monomial_1 = FuncUtil.Ctermfunc.empty;
   150 
   151 fun monomial_var x = FuncUtil.Ctermfunc.onefunc (x, 1);
   152 
   153 val monomial_mul =
   154   FuncUtil.Ctermfunc.combine Integer.add (K false);
   155 
   156 fun monomial_multidegree m =
   157   FuncUtil.Ctermfunc.fold (fn (_, k) => fn a => k + a) m 0;
   158 
   159 fun monomial_variables m = FuncUtil.Ctermfunc.dom m;
   160 
   161 (* Polynomials.                                                              *)
   162 
   163 fun eval assig p =
   164   FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => a + c * monomial_eval assig m) p @0;
   165 
   166 val poly_0 = FuncUtil.Monomialfunc.empty;
   167 
   168 fun poly_isconst p =
   169   FuncUtil.Monomialfunc.fold (fn (m, _) => fn a => FuncUtil.Ctermfunc.is_empty m andalso a)
   170     p true;
   171 
   172 fun poly_var x = FuncUtil.Monomialfunc.onefunc (monomial_var x, @1);
   173 
   174 fun poly_const c =
   175   if c = @0 then poly_0 else FuncUtil.Monomialfunc.onefunc (monomial_1, c);
   176 
   177 fun poly_cmul c p =
   178   if c = @0 then poly_0
   179   else FuncUtil.Monomialfunc.map (fn _ => fn x => c * x) p;
   180 
   181 fun poly_neg p = FuncUtil.Monomialfunc.map (K ~) p;
   182 
   183 
   184 fun poly_add p1 p2 =
   185   FuncUtil.Monomialfunc.combine (curry op +) (fn x => x = @0) p1 p2;
   186 
   187 fun poly_sub p1 p2 = poly_add p1 (poly_neg p2);
   188 
   189 fun poly_cmmul (c,m) p =
   190   if c = @0 then poly_0
   191   else
   192     if FuncUtil.Ctermfunc.is_empty m
   193     then FuncUtil.Monomialfunc.map (fn _ => fn d => c * d) p
   194     else
   195       FuncUtil.Monomialfunc.fold (fn (m', d) => fn a =>
   196           (FuncUtil.Monomialfunc.update (monomial_mul m m', c * d) a)) p poly_0;
   197 
   198 fun poly_mul p1 p2 =
   199   FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0;
   200 
   201 fun poly_square p = poly_mul p p;
   202 
   203 fun poly_pow p k =
   204   if k = 0 then poly_const @1
   205   else if k = 1 then p
   206   else
   207     let val q = poly_square(poly_pow p (k div 2))
   208     in if k mod 2 = 1 then poly_mul p q else q end;
   209 
   210 fun multidegree p =
   211   FuncUtil.Monomialfunc.fold (fn (m, _) => fn a => max (monomial_multidegree m) a) p 0;
   212 
   213 fun poly_variables p =
   214   sort FuncUtil.cterm_ord
   215     (FuncUtil.Monomialfunc.fold_rev
   216       (fn (m, _) => union (is_equal o FuncUtil.cterm_ord) (monomial_variables m)) p []);
   217 
   218 (* Conversion from HOL term.                                                 *)
   219 
   220 local
   221   val neg_tm = @{cterm "uminus :: real => _"}
   222   val add_tm = @{cterm "op + :: real => _"}
   223   val sub_tm = @{cterm "op - :: real => _"}
   224   val mul_tm = @{cterm "op * :: real => _"}
   225   val inv_tm = @{cterm "inverse :: real => _"}
   226   val div_tm = @{cterm "op / :: real => _"}
   227   val pow_tm = @{cterm "op ^ :: real => _"}
   228   val zero_tm = @{cterm "0:: real"}
   229   val is_numeral = can (HOLogic.dest_number o Thm.term_of)
   230   fun poly_of_term tm =
   231     if tm aconvc zero_tm then poly_0
   232     else
   233       if RealArith.is_ratconst tm
   234       then poly_const(RealArith.dest_ratconst tm)
   235       else
   236        (let
   237           val (lop, r) = Thm.dest_comb tm
   238         in
   239           if lop aconvc neg_tm then poly_neg(poly_of_term r)
   240           else if lop aconvc inv_tm then
   241             let val p = poly_of_term r in
   242               if poly_isconst p
   243               then poly_const(Rat.inv (eval FuncUtil.Ctermfunc.empty p))
   244               else error "poly_of_term: inverse of non-constant polyomial"
   245             end
   246           else
   247            (let
   248               val (opr,l) = Thm.dest_comb lop
   249             in
   250               if opr aconvc pow_tm andalso is_numeral r
   251               then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o Thm.term_of) r)
   252               else if opr aconvc add_tm
   253               then poly_add (poly_of_term l) (poly_of_term r)
   254               else if opr aconvc sub_tm
   255               then poly_sub (poly_of_term l) (poly_of_term r)
   256               else if opr aconvc mul_tm
   257               then poly_mul (poly_of_term l) (poly_of_term r)
   258               else if opr aconvc div_tm
   259               then
   260                 let
   261                   val p = poly_of_term l
   262                   val q = poly_of_term r
   263                 in
   264                   if poly_isconst q
   265                   then poly_cmul (Rat.inv (eval FuncUtil.Ctermfunc.empty q)) p
   266                   else error "poly_of_term: division by non-constant polynomial"
   267                 end
   268               else poly_var tm
   269             end handle CTERM ("dest_comb",_) => poly_var tm)
   270         end handle CTERM ("dest_comb",_) => poly_var tm)
   271 in
   272   val poly_of_term = fn tm =>
   273     if type_of (Thm.term_of tm) = @{typ real}
   274     then poly_of_term tm
   275     else error "poly_of_term: term does not have real type"
   276 end;
   277 
   278 
   279 (* String of vector (just a list of space-separated numbers). *)
   280 
   281 fun sdpa_of_vector (v: vector) =
   282   let
   283     val n = dim v
   284     val strs =
   285       map (decimalize 20 o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i @0)) (1 upto n)
   286   in space_implode " " strs ^ "\n" end;
   287 
   288 fun triple_int_ord ((a, b, c), (a', b', c')) =
   289   prod_ord int_ord (prod_ord int_ord int_ord) ((a, (b, c)), (a', (b', c')));
   290 structure Inttriplefunc = FuncFun(type key = int * int * int val ord = triple_int_ord);
   291 
   292 
   293 (* Parse back csdp output. *)
   294 
   295 local
   296 
   297 val decimal_digits = Scan.many1 Symbol.is_ascii_digit
   298 val decimal_nat = decimal_digits >> (#1 o Library.read_int);
   299 val decimal_int = decimal_nat >> Rat.of_int;
   300 
   301 val decimal_sig =
   302   decimal_int -- Scan.option (Scan.$$ "." |-- decimal_digits) >>
   303   (fn (a, NONE) => a
   304     | (a, SOME bs) => a + Rat.of_int (#1 (Library.read_int bs)) / rat_pow @10 (length bs));
   305 
   306 fun signed neg parse = $$ "-" |-- parse >> neg || $$ "+" |-- parse || parse;
   307 val exponent = ($$ "e" || $$ "E") |-- signed ~ decimal_nat;
   308 
   309 val decimal =
   310   signed ~ decimal_sig -- Scan.optional exponent 0
   311     >> (fn (a, b) => a * rat_pow @10 b);
   312 
   313 val csdp_output =
   314   decimal -- Scan.repeat (Scan.$$ " " |-- Scan.option decimal) --| Scan.many Symbol.not_eof
   315     >> (fn (a, bs) => vector_of_list (a :: map_filter I bs));
   316 
   317 in
   318 
   319 fun parse_csdpoutput s =
   320   Symbol.scanner "Malformed CSDP output" csdp_output (raw_explode s);
   321 
   322 end;
   323 
   324 
   325 (* Try some apparently sensible scaling first. Note that this is purely to   *)
   326 (* get a cleaner translation to floating-point, and doesn't affect any of    *)
   327 (* the results, in principle. In practice it seems a lot better when there   *)
   328 (* are extreme numbers in the original problem.                              *)
   329 
   330 (* Version for (int*int*int) keys *)
   331 local
   332   fun max_rat x y = if x < y then y else x
   333   fun common_denominator fld amat acc =
   334     fld (fn (_,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
   335   fun maximal_element fld amat acc =
   336     fld (fn (_,c) => fn maxa => max_rat maxa (abs c)) amat acc
   337   fun float_of_rat x =
   338     let val (a,b) = Rat.dest x
   339     in Real.fromInt a / Real.fromInt b end;
   340   fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0)
   341 in
   342 
   343 fun tri_scale_then solver (obj:vector) mats =
   344   let
   345     val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats @1
   346     val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj) @1
   347     val mats' = map (Inttriplefunc.map (fn _ => fn x => cd1 * x)) mats
   348     val obj' = vector_cmul cd2 obj
   349     val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' @0
   350     val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') @0
   351     val scal1 = rat_pow @2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0))
   352     val scal2 = rat_pow @2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0))
   353     val mats'' = map (Inttriplefunc.map (fn _ => fn x => x * scal1)) mats'
   354     val obj'' = vector_cmul scal2 obj'
   355   in solver obj'' mats'' end
   356 end;
   357 
   358 (* Round a vector to "nice" rationals.                                       *)
   359 
   360 fun nice_rational n x = round_rat (n * x) / n;
   361 fun nice_vector n ((d,v) : vector) =
   362   (d, FuncUtil.Intfunc.fold (fn (i,c) => fn a =>
   363       let val y = nice_rational n c in
   364         if c = @0 then a
   365         else FuncUtil.Intfunc.update (i,y) a
   366       end) v FuncUtil.Intfunc.empty): vector
   367 
   368 fun dest_ord f x = is_equal (f x);
   369 
   370 (* Stuff for "equations" ((int*int*int)->num functions).                         *)
   371 
   372 fun tri_equation_cmul c eq =
   373   if c = @0 then Inttriplefunc.empty
   374   else Inttriplefunc.map (fn _ => fn d => c * d) eq;
   375 
   376 fun tri_equation_add eq1 eq2 =
   377   Inttriplefunc.combine (curry op +) (fn x => x = @0) eq1 eq2;
   378 
   379 fun tri_equation_eval assig eq =
   380   let
   381     fun value v = Inttriplefunc.apply assig v
   382   in Inttriplefunc.fold (fn (v, c) => fn a => a + value v * c) eq @0 end;
   383 
   384 (* Eliminate all variables, in an essentially arbitrary order.               *)
   385 
   386 fun tri_eliminate_all_equations one =
   387   let
   388     fun choose_variable eq =
   389       let val (v,_) = Inttriplefunc.choose eq
   390       in
   391         if is_equal (triple_int_ord(v,one)) then
   392           let
   393             val eq' = Inttriplefunc.delete_safe v eq
   394           in
   395             if Inttriplefunc.is_empty eq' then error "choose_variable"
   396             else fst (Inttriplefunc.choose eq')
   397           end
   398         else v
   399       end
   400 
   401     fun eliminate dun eqs =
   402       (case eqs of
   403         [] => dun
   404       | eq :: oeqs =>
   405           if Inttriplefunc.is_empty eq then eliminate dun oeqs
   406           else
   407             let
   408               val v = choose_variable eq
   409               val a = Inttriplefunc.apply eq v
   410               val eq' =
   411                 tri_equation_cmul ((Rat.of_int ~1) / a) (Inttriplefunc.delete_safe v eq)
   412               fun elim e =
   413                 let val b = Inttriplefunc.tryapplyd e v @0 in
   414                   if b = @0 then e
   415                   else tri_equation_add e (tri_equation_cmul (~ b / a) eq)
   416                 end
   417             in
   418               eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun))
   419                 (map elim oeqs)
   420             end)
   421   in
   422     fn eqs =>
   423       let
   424         val assig = eliminate Inttriplefunc.empty eqs
   425         val vs = Inttriplefunc.fold (fn (_, f) => fn a =>
   426           remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
   427       in (distinct (dest_ord triple_int_ord) vs,assig) end
   428   end;
   429 
   430 (* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
   431 
   432 fun tri_epoly_pmul p q acc =
   433   FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a =>
   434     FuncUtil.Monomialfunc.fold (fn (m2, e) => fn b =>
   435       let
   436         val m =  monomial_mul m1 m2
   437         val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty
   438       in
   439         FuncUtil.Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b
   440       end) q a) p acc;
   441 
   442 (* Hence produce the "relevant" monomials: those whose squares lie in the    *)
   443 (* Newton polytope of the monomials in the input. (This is enough according  *)
   444 (* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal,       *)
   445 (* vol 45, pp. 363--374, 1978.                                               *)
   446 (*                                                                           *)
   447 (* These are ordered in sort of decreasing degree. In particular the         *)
   448 (* constant monomial is last; this gives an order in diagonalization of the  *)
   449 (* quadratic form that will tend to display constants.                       *)
   450 
   451 (* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form.  *)
   452 
   453 local
   454   fun diagonalize n i m =
   455     if FuncUtil.Intpairfunc.is_empty (snd m) then []
   456     else
   457       let
   458         val a11 = FuncUtil.Intpairfunc.tryapplyd (snd m) (i,i) @0
   459       in
   460         if a11 < @0 then raise Failure "diagonalize: not PSD"
   461         else if a11 = @0 then
   462           if FuncUtil.Intfunc.is_empty (snd (row i m))
   463           then diagonalize n (i + 1) m
   464           else raise Failure "diagonalize: not PSD ___ "
   465         else
   466           let
   467             val v = row i m
   468             val v' =
   469               (fst v, FuncUtil.Intfunc.fold (fn (i, c) => fn a =>
   470                 let val y = c / a11
   471                 in if y = @0 then a else FuncUtil.Intfunc.update (i,y) a
   472                 end) (snd v) FuncUtil.Intfunc.empty)
   473             fun upt0 x y a =
   474               if y = @0 then a
   475               else FuncUtil.Intpairfunc.update (x,y) a
   476             val m' =
   477               ((n, n),
   478                 iter (i + 1, n) (fn j =>
   479                   iter (i + 1, n) (fn k =>
   480                     (upt0 (j, k)
   481                       (FuncUtil.Intpairfunc.tryapplyd (snd m) (j, k) @0 -
   482                         FuncUtil.Intfunc.tryapplyd (snd v) j @0 *
   483                         FuncUtil.Intfunc.tryapplyd (snd v') k @0))))
   484                     FuncUtil.Intpairfunc.empty)
   485           in (a11, v') :: diagonalize n (i + 1) m' end
   486       end
   487 in
   488   fun diag m =
   489     let
   490       val nn = dimensions m
   491       val n = fst nn
   492     in
   493       if snd nn <> n then error "diagonalize: non-square matrix"
   494       else diagonalize n 1 m
   495     end
   496 end;
   497 
   498 (* Enumeration of monomials with given multidegree bound.                    *)
   499 
   500 fun enumerate_monomials d vars =
   501   if d < 0 then []
   502   else if d = 0 then [FuncUtil.Ctermfunc.empty]
   503   else if null vars then [monomial_1]
   504   else
   505     let val alts =
   506       map_range (fn k =>
   507         let
   508           val oths = enumerate_monomials (d - k) (tl vars)
   509         in map (fn ks => if k = 0 then ks else FuncUtil.Ctermfunc.update (hd vars, k) ks) oths end)
   510         (d + 1)
   511   in flat alts end;
   512 
   513 (* Enumerate products of distinct input polys with degree <= d.              *)
   514 (* We ignore any constant input polynomials.                                 *)
   515 (* Give the output polynomial and a record of how it was derived.            *)
   516 
   517 fun enumerate_products d pols =
   518   if d = 0 then [(poly_const @1, RealArith.Rational_lt @1)]
   519   else if d < 0 then []
   520   else
   521     (case pols of
   522       [] => [(poly_const @1, RealArith.Rational_lt @1)]
   523     | (p, b) :: ps =>
   524         let val e = multidegree p in
   525           if e = 0 then enumerate_products d ps
   526           else
   527             enumerate_products d ps @
   528             map (fn (q, c) => (poly_mul p q, RealArith.Product (b, c)))
   529               (enumerate_products (d - e) ps)
   530         end)
   531 
   532 (* Convert regular polynomial. Note that we treat (0,0,0) as -1.             *)
   533 
   534 fun epoly_of_poly p =
   535   FuncUtil.Monomialfunc.fold (fn (m, c) => fn a =>
   536       FuncUtil.Monomialfunc.update (m, Inttriplefunc.onefunc ((0, 0, 0), ~ c)) a)
   537     p FuncUtil.Monomialfunc.empty;
   538 
   539 (* String for block diagonal matrix numbered k.                              *)
   540 
   541 fun sdpa_of_blockdiagonal k m =
   542   let
   543     val pfx = string_of_int k ^" "
   544     val ents =
   545       Inttriplefunc.fold
   546         (fn ((b, i, j), c) => fn a => if i > j then a else ((b, i, j), c) :: a)
   547         m []
   548     val entss = sort (triple_int_ord o apply2 fst) ents
   549   in
   550     fold_rev (fn ((b,i,j),c) => fn a =>
   551       pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
   552       " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
   553   end;
   554 
   555 (* SDPA for problem using block diagonal (i.e. multiple SDPs)                *)
   556 
   557 fun sdpa_of_blockproblem nblocks blocksizes obj mats =
   558   let val m = length mats - 1
   559   in
   560     string_of_int m ^ "\n" ^
   561     string_of_int nblocks ^ "\n" ^
   562     (space_implode " " (map string_of_int blocksizes)) ^
   563     "\n" ^
   564     sdpa_of_vector obj ^
   565     fold_rev (fn (k, m) => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a)
   566       (1 upto length mats ~~ mats) ""
   567   end;
   568 
   569 (* Run prover on a problem in block diagonal form.                       *)
   570 
   571 fun run_blockproblem prover nblocks blocksizes obj mats =
   572   parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats))
   573 
   574 (* 3D versions of matrix operations to consider blocks separately.           *)
   575 
   576 val bmatrix_add = Inttriplefunc.combine (curry op +) (fn x => x = @0);
   577 fun bmatrix_cmul c bm =
   578   if c = @0 then Inttriplefunc.empty
   579   else Inttriplefunc.map (fn _ => fn x => c * x) bm;
   580 
   581 val bmatrix_neg = bmatrix_cmul (Rat.of_int ~1);
   582 
   583 (* Smash a block matrix into components.                                     *)
   584 
   585 fun blocks blocksizes bm =
   586   map (fn (bs, b0) =>
   587     let
   588       val m =
   589         Inttriplefunc.fold
   590           (fn ((b, i, j), c) => fn a =>
   591             if b = b0 then FuncUtil.Intpairfunc.update ((i, j), c) a else a)
   592         bm FuncUtil.Intpairfunc.empty
   593       val _ = FuncUtil.Intpairfunc.fold (fn ((i, j), _) => fn a => max a (max i j)) m 0
   594     in (((bs, bs), m): matrix) end)
   595   (blocksizes ~~ (1 upto length blocksizes));
   596 
   597 (* FIXME : Get rid of this !!!*)
   598 local
   599   fun tryfind_with msg _ [] = raise Failure msg
   600     | tryfind_with _ f (x::xs) = (f x handle Failure s => tryfind_with s f xs);
   601 in
   602   fun tryfind f = tryfind_with "tryfind" f
   603 end
   604 
   605 (* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
   606 
   607 fun real_positivnullstellensatz_general ctxt prover linf d eqs leqs pol =
   608   let
   609     val vars =
   610       fold_rev (union (op aconvc) o poly_variables)
   611         (pol :: eqs @ map fst leqs) []
   612     val monoid =
   613       if linf then
   614         (poly_const @1, RealArith.Rational_lt @1)::
   615         (filter (fn (p,_) => multidegree p <= d) leqs)
   616       else enumerate_products d leqs
   617     val nblocks = length monoid
   618     fun mk_idmultiplier k p =
   619       let
   620         val e = d - multidegree p
   621         val mons = enumerate_monomials e vars
   622         val nons = mons ~~ (1 upto length mons)
   623       in
   624         (mons,
   625           fold_rev (fn (m, n) =>
   626             FuncUtil.Monomialfunc.update (m, Inttriplefunc.onefunc ((~k, ~n, n), @1)))
   627           nons FuncUtil.Monomialfunc.empty)
   628       end
   629 
   630     fun mk_sqmultiplier k (p,_) =
   631       let
   632         val e = (d - multidegree p) div 2
   633         val mons = enumerate_monomials e vars
   634         val nons = mons ~~ (1 upto length mons)
   635       in
   636         (mons,
   637           fold_rev (fn (m1, n1) =>
   638             fold_rev (fn (m2, n2) => fn a =>
   639               let val m = monomial_mul m1 m2 in
   640                 if n1 > n2 then a
   641                 else
   642                   let
   643                     val c = if n1 = n2 then @1 else @2
   644                     val e = FuncUtil.Monomialfunc.tryapplyd a m Inttriplefunc.empty
   645                   in
   646                     FuncUtil.Monomialfunc.update
   647                       (m, tri_equation_add (Inttriplefunc.onefunc ((k, n1, n2), c)) e) a
   648                   end
   649               end) nons) nons FuncUtil.Monomialfunc.empty)
   650       end
   651 
   652     val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid)
   653     val (_(*idmonlist*),ids) =  split_list (map2 mk_idmultiplier (1 upto length eqs) eqs)
   654     val blocksizes = map length sqmonlist
   655     val bigsum =
   656       fold_rev (fn (p, q) => fn a => tri_epoly_pmul p q a) (eqs ~~ ids)
   657         (fold_rev (fn ((p, _), s) => fn a => tri_epoly_pmul p s a) (monoid ~~ sqs)
   658           (epoly_of_poly (poly_neg pol)))
   659     val eqns = FuncUtil.Monomialfunc.fold (fn (_, e) => fn a => e :: a) bigsum []
   660     val (pvs, assig) = tri_eliminate_all_equations (0, 0, 0) eqns
   661     val qvars = (0, 0, 0) :: pvs
   662     val allassig =
   663       fold_rev (fn v => Inttriplefunc.update (v, (Inttriplefunc.onefunc (v, @1)))) pvs assig
   664     fun mk_matrix v =
   665       Inttriplefunc.fold (fn ((b, i, j), ass) => fn m =>
   666           if b < 0 then m
   667           else
   668             let val c = Inttriplefunc.tryapplyd ass v @0 in
   669               if c = @0 then m
   670               else Inttriplefunc.update ((b, j, i), c) (Inttriplefunc.update ((b, i, j), c) m)
   671             end)
   672         allassig Inttriplefunc.empty
   673     val diagents =
   674       Inttriplefunc.fold
   675         (fn ((b, i, j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a)
   676         allassig Inttriplefunc.empty
   677 
   678     val mats = map mk_matrix qvars
   679     val obj =
   680       (length pvs,
   681         itern 1 pvs (fn v => fn i =>
   682           FuncUtil.Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v @0))
   683           FuncUtil.Intfunc.empty)
   684     val raw_vec =
   685       if null pvs then vector_0 0
   686       else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats
   687     fun int_element (_, v) i = FuncUtil.Intfunc.tryapplyd v i @0
   688 
   689     fun find_rounding d =
   690       let
   691         val _ =
   692           debug_message ctxt (fn () => "Trying rounding with limit "^Rat.string_of_rat d ^ "\n")
   693         val vec = nice_vector d raw_vec
   694         val blockmat =
   695           iter (1, dim vec)
   696             (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a)
   697             (bmatrix_neg (nth mats 0))
   698         val allmats = blocks blocksizes blockmat
   699       in (vec, map diag allmats) end
   700     val (vec, ratdias) =
   701       if null pvs then find_rounding @1
   702       else tryfind find_rounding (map Rat.of_int (1 upto 31) @ map (rat_pow @2) (5 upto 66))
   703     val newassigs =
   704       fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k))
   705         (1 upto dim vec) (Inttriplefunc.onefunc ((0, 0, 0), Rat.of_int ~1))
   706     val finalassigs =
   707       Inttriplefunc.fold (fn (v, e) => fn a =>
   708         Inttriplefunc.update (v, tri_equation_eval newassigs e) a) allassig newassigs
   709     fun poly_of_epoly p =
   710       FuncUtil.Monomialfunc.fold (fn (v, e) => fn a =>
   711           FuncUtil.Monomialfunc.updatep iszero (v, tri_equation_eval finalassigs e) a)
   712         p FuncUtil.Monomialfunc.empty
   713     fun mk_sos mons =
   714       let
   715         fun mk_sq (c, m) =
   716           (c, fold_rev (fn k => fn a =>
   717               FuncUtil.Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a)
   718             (1 upto length mons) FuncUtil.Monomialfunc.empty)
   719       in map mk_sq end
   720     val sqs = map2 mk_sos sqmonlist ratdias
   721     val cfs = map poly_of_epoly ids
   722     val msq = filter (fn (_, b) => not (null b)) (map2 pair monoid sqs)
   723     fun eval_sq sqs = fold_rev (fn (c, q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0
   724     val sanity =
   725       fold_rev (fn ((p, _), s) => poly_add (poly_mul p (eval_sq s))) msq
   726         (fold_rev (fn (p, q) => poly_add (poly_mul p q)) (cfs ~~ eqs) (poly_neg pol))
   727   in
   728     if not(FuncUtil.Monomialfunc.is_empty sanity) then raise Sanity
   729     else (cfs, map (fn (a, b) => (snd a, b)) msq)
   730   end
   731 
   732 
   733 (* Iterative deepening.                                                      *)
   734 
   735 fun deepen ctxt f n =
   736   (trace_message ctxt (fn () => "Searching with depth limit " ^ string_of_int n);
   737     (f n handle Failure s =>
   738       (trace_message ctxt (fn () => "failed with message: " ^ s); deepen ctxt f (n + 1))));
   739 
   740 
   741 (* Map back polynomials and their composites to a positivstellensatz.        *)
   742 
   743 fun cterm_of_sqterm (c, p) = RealArith.Product (RealArith.Rational_lt c, RealArith.Square p);
   744 
   745 fun cterm_of_sos (pr,sqs) =
   746   if null sqs then pr
   747   else RealArith.Product (pr, foldr1 RealArith.Sum (map cterm_of_sqterm sqs));
   748 
   749 (* Interface to HOL.                                                         *)
   750 local
   751   open Conv
   752   val concl = Thm.dest_arg o Thm.cprop_of
   753   fun simple_cterm_ord t u = Term_Ord.fast_term_ord (Thm.term_of t, Thm.term_of u) = LESS
   754 in
   755 (* FIXME: Replace tryfind by get_first !! *)
   756 fun real_nonlinear_prover proof_method ctxt =
   757   let
   758     val {add = _, mul = _, neg = _, pow = _, sub = _, main = real_poly_conv} =
   759       Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt
   760         (the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"}))
   761         simple_cterm_ord
   762     fun mainf cert_choice translator (eqs, les, lts) =
   763       let
   764         val eq0 = map (poly_of_term o Thm.dest_arg1 o concl) eqs
   765         val le0 = map (poly_of_term o Thm.dest_arg o concl) les
   766         val lt0 = map (poly_of_term o Thm.dest_arg o concl) lts
   767         val eqp0 = map_index (fn (i, t) => (t,RealArith.Axiom_eq i)) eq0
   768         val lep0 = map_index (fn (i, t) => (t,RealArith.Axiom_le i)) le0
   769         val ltp0 = map_index (fn (i, t) => (t,RealArith.Axiom_lt i)) lt0
   770         val (keq,eq) = List.partition (fn (p, _) => multidegree p = 0) eqp0
   771         val (klep,lep) = List.partition (fn (p, _) => multidegree p = 0) lep0
   772         val (kltp,ltp) = List.partition (fn (p, _) => multidegree p = 0) ltp0
   773         fun trivial_axiom (p, ax) =
   774           (case ax of
   775             RealArith.Axiom_eq n =>
   776               if eval FuncUtil.Ctermfunc.empty p <> @0 then nth eqs n
   777               else raise Failure "trivial_axiom: Not a trivial axiom"
   778           | RealArith.Axiom_le n =>
   779               if eval FuncUtil.Ctermfunc.empty p < @0 then nth les n
   780               else raise Failure "trivial_axiom: Not a trivial axiom"
   781           | RealArith.Axiom_lt n =>
   782               if eval FuncUtil.Ctermfunc.empty p <= @0 then nth lts n
   783               else raise Failure "trivial_axiom: Not a trivial axiom"
   784           | _ => error "trivial_axiom: Not a trivial axiom")
   785       in
   786         let val th = tryfind trivial_axiom (keq @ klep @ kltp) in
   787           (fconv_rule (arg_conv (arg1_conv (real_poly_conv ctxt))
   788             then_conv Numeral_Simprocs.field_comp_conv ctxt) th,
   789             RealArith.Trivial)
   790         end handle Failure _ =>
   791           let
   792             val proof =
   793               (case proof_method of
   794                 Certificate certs =>
   795                   (* choose certificate *)
   796                   let
   797                     fun chose_cert [] (RealArith.Cert c) = c
   798                       | chose_cert (RealArith.Left::s) (RealArith.Branch (l, _)) = chose_cert s l
   799                       | chose_cert (RealArith.Right::s) (RealArith.Branch (_, r)) = chose_cert s r
   800                       | chose_cert _ _ = error "certificate tree in invalid form"
   801                   in
   802                     chose_cert cert_choice certs
   803                   end
   804               | Prover prover =>
   805                   (* call prover *)
   806                   let
   807                     val pol = fold_rev poly_mul (map fst ltp) (poly_const @1)
   808                     val leq = lep @ ltp
   809                     fun tryall d =
   810                       let
   811                         val e = multidegree pol
   812                         val k = if e = 0 then 0 else d div e
   813                         val eq' = map fst eq
   814                       in
   815                         tryfind (fn i =>
   816                             (d, i, real_positivnullstellensatz_general ctxt prover false d eq' leq
   817                               (poly_neg(poly_pow pol i))))
   818                           (0 upto k)
   819                       end
   820                     val (_,i,(cert_ideal,cert_cone)) = deepen ctxt tryall 0
   821                     val proofs_ideal =
   822                       map2 (fn q => fn (_,ax) => RealArith.Eqmul(q,ax)) cert_ideal eq
   823                     val proofs_cone = map cterm_of_sos cert_cone
   824                     val proof_ne =
   825                       if null ltp then RealArith.Rational_lt @1
   826                       else
   827                         let val p = foldr1 RealArith.Product (map snd ltp) in
   828                           funpow i (fn q => RealArith.Product (p, q))
   829                             (RealArith.Rational_lt @1)
   830                         end
   831                   in
   832                     foldr1 RealArith.Sum (proof_ne :: proofs_ideal @ proofs_cone)
   833                   end)
   834           in
   835             (translator (eqs,les,lts) proof, RealArith.Cert proof)
   836           end
   837       end
   838   in mainf end
   839 end
   840 
   841 (* FIXME : This is very bad!!!*)
   842 fun subst_conv eqs t =
   843   let
   844     val t' = fold (Thm.lambda o Thm.lhs_of) eqs t
   845   in
   846     Conv.fconv_rule (Thm.beta_conversion true)
   847       (fold (fn a => fn b => Thm.combination b a) eqs (Thm.reflexive t'))
   848   end
   849 
   850 (* A wrapper that tries to substitute away variables first.                  *)
   851 
   852 local
   853   open Conv
   854   fun simple_cterm_ord t u = Term_Ord.fast_term_ord (Thm.term_of t, Thm.term_of u) = LESS
   855   val concl = Thm.dest_arg o Thm.cprop_of
   856   val shuffle1 =
   857     fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))"
   858       by (atomize (full)) (simp add: field_simps)})
   859   val shuffle2 =
   860     fconv_rule (rewr_conv @{lemma "(x + a == y) ==  (x == y - (a::real))"
   861       by (atomize (full)) (simp add: field_simps)})
   862   fun substitutable_monomial fvs tm =
   863     (case Thm.term_of tm of
   864       Free (_, @{typ real}) =>
   865         if not (member (op aconvc) fvs tm) then (@1, tm)
   866         else raise Failure "substitutable_monomial"
   867     | @{term "op * :: real => _"} $ _ $ (Free _) =>
   868         if RealArith.is_ratconst (Thm.dest_arg1 tm) andalso
   869           not (member (op aconvc) fvs (Thm.dest_arg tm))
   870         then (RealArith.dest_ratconst (Thm.dest_arg1 tm), Thm.dest_arg tm)
   871         else raise Failure "substitutable_monomial"
   872     | @{term "op + :: real => _"}$_$_ =>
   873          (substitutable_monomial (Drule.cterm_add_frees (Thm.dest_arg tm) fvs) (Thm.dest_arg1 tm)
   874            handle Failure _ =>
   875             substitutable_monomial (Drule.cterm_add_frees (Thm.dest_arg1 tm) fvs) (Thm.dest_arg tm))
   876     | _ => raise Failure "substitutable_monomial")
   877 
   878   fun isolate_variable v th =
   879     let
   880       val w = Thm.dest_arg1 (Thm.cprop_of th)
   881     in
   882       if v aconvc w then th
   883       else
   884         (case Thm.term_of w of
   885           @{term "op + :: real => _"} $ _ $ _ =>
   886             if Thm.dest_arg1 w aconvc v then shuffle2 th
   887             else isolate_variable v (shuffle1 th)
   888         | _ => error "isolate variable : This should not happen?")
   889    end
   890 in
   891 
   892 fun real_nonlinear_subst_prover prover ctxt =
   893   let
   894     val {add = _, mul = real_poly_mul_conv, neg = _, pow = _, sub = _, main = real_poly_conv} =
   895       Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt
   896         (the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"}))
   897         simple_cterm_ord
   898 
   899     fun make_substitution th =
   900       let
   901         val (c,v) = substitutable_monomial [] (Thm.dest_arg1(concl th))
   902         val th1 =
   903           Drule.arg_cong_rule
   904             (Thm.apply @{cterm "op * :: real => _"} (RealArith.cterm_of_rat (Rat.inv c)))
   905             (mk_meta_eq th)
   906         val th2 = fconv_rule (binop_conv (real_poly_mul_conv ctxt)) th1
   907       in fconv_rule (arg_conv (real_poly_conv ctxt)) (isolate_variable v th2) end
   908 
   909     fun oprconv cv ct =
   910       let val g = Thm.dest_fun2 ct in
   911         if g aconvc @{cterm "op <= :: real => _"} orelse g aconvc @{cterm "op < :: real => _"}
   912         then arg_conv cv ct else arg1_conv cv ct
   913       end
   914     fun mainf cert_choice translator =
   915       let
   916         fun substfirst (eqs, les, lts) =
   917           (let
   918               val eth = tryfind make_substitution eqs
   919               val modify =
   920                 fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv (real_poly_conv ctxt))))
   921             in
   922               substfirst
   923                 (filter_out
   924                   (fn t => (Thm.dest_arg1 o Thm.dest_arg o Thm.cprop_of) t aconvc @{cterm "0::real"})
   925                   (map modify eqs),
   926                   map modify les,
   927                   map modify lts)
   928             end handle Failure  _ =>
   929               real_nonlinear_prover prover ctxt cert_choice translator (rev eqs, rev les, rev lts))
   930       in substfirst end
   931   in mainf end
   932 
   933 (* Overall function. *)
   934 
   935 fun real_sos prover ctxt =
   936   RealArith.gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt)
   937 
   938 end;
   939 
   940 val known_sos_constants =
   941   [@{term "op ==>"}, @{term "Trueprop"},
   942    @{term HOL.False}, @{term HOL.implies}, @{term HOL.conj}, @{term HOL.disj},
   943    @{term "Not"}, @{term "op = :: bool => _"},
   944    @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"},
   945    @{term "op = :: real => _"}, @{term "op < :: real => _"},
   946    @{term "op <= :: real => _"},
   947    @{term "op + :: real => _"}, @{term "op - :: real => _"},
   948    @{term "op * :: real => _"}, @{term "uminus :: real => _"},
   949    @{term "op / :: real => _"}, @{term "inverse :: real => _"},
   950    @{term "op ^ :: real => _"}, @{term "abs :: real => _"},
   951    @{term "min :: real => _"}, @{term "max :: real => _"},
   952    @{term "0::real"}, @{term "1::real"},
   953    @{term "numeral :: num => nat"},
   954    @{term "numeral :: num => real"},
   955    @{term "Num.Bit0"}, @{term "Num.Bit1"}, @{term "Num.One"}];
   956 
   957 fun check_sos kcts ct =
   958   let
   959     val t = Thm.term_of ct
   960     val _ =
   961       if not (null (Term.add_tfrees t []) andalso null (Term.add_tvars t []))
   962       then error "SOS: not sos. Additional type varables"
   963       else ()
   964     val fs = Term.add_frees t []
   965     val _ =
   966       if exists (fn ((_,T)) => not (T = @{typ "real"})) fs
   967       then error "SOS: not sos. Variables with type not real"
   968       else ()
   969     val vs = Term.add_vars t []
   970     val _ =
   971       if exists (fn ((_,T)) => not (T = @{typ "real"})) vs
   972       then error "SOS: not sos. Variables with type not real"
   973       else ()
   974     val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
   975     val _ =
   976       if null ukcs then ()
   977       else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
   978   in () end
   979 
   980 fun core_sos_tac print_cert prover = SUBPROOF (fn {concl, context = ctxt, ...} =>
   981   let
   982     val _ = check_sos known_sos_constants concl
   983     val (th, certificates) = real_sos prover ctxt (Thm.dest_arg concl)
   984     val _ = print_cert certificates
   985   in resolve_tac ctxt [th] 1 end);
   986 
   987 fun default_SOME _ NONE v = SOME v
   988   | default_SOME _ (SOME v) _ = SOME v;
   989 
   990 fun lift_SOME f NONE a = f a
   991   | lift_SOME _ (SOME a) _ = SOME a;
   992 
   993 
   994 local
   995   val is_numeral = can (HOLogic.dest_number o Thm.term_of)
   996 in
   997   fun get_denom b ct =
   998     (case Thm.term_of ct of
   999       @{term "op / :: real => _"} $ _ $ _ =>
  1000         if is_numeral (Thm.dest_arg ct)
  1001         then get_denom b (Thm.dest_arg1 ct)
  1002         else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct)) (Thm.dest_arg ct, b)
  1003     | @{term "op < :: real => _"} $ _ $ _ =>
  1004         lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
  1005     | @{term "op <= :: real => _"} $ _ $ _ =>
  1006         lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
  1007     | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct)
  1008     | _ => NONE)
  1009 end;
  1010 
  1011 fun elim_one_denom_tac ctxt = CSUBGOAL (fn (P, i) =>
  1012   (case get_denom false P of
  1013     NONE => no_tac
  1014   | SOME (d, ord) =>
  1015       let
  1016         val simp_ctxt =
  1017           ctxt addsimps @{thms field_simps}
  1018           addsimps [@{thm power_divide}]
  1019         val th =
  1020           Thm.instantiate' [] [SOME d, SOME (Thm.dest_arg P)]
  1021             (if ord then @{lemma "(d=0 \<longrightarrow> P) \<and> (d>0 \<longrightarrow> P) \<and> (d<(0::real) \<longrightarrow> P) \<Longrightarrow> P" by auto}
  1022              else @{lemma "(d=0 \<longrightarrow> P) \<and> (d \<noteq> (0::real) \<longrightarrow> P) \<Longrightarrow> P" by blast})
  1023       in resolve_tac ctxt [th] i THEN Simplifier.asm_full_simp_tac simp_ctxt i end));
  1024 
  1025 fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
  1026 
  1027 fun sos_tac print_cert prover ctxt =
  1028   (* The SOS prover breaks if mult_nonneg_nonneg is in the simpset *)
  1029   let val ctxt' = Context_Position.set_visible false ctxt delsimps @{thms mult_nonneg_nonneg}
  1030   in Object_Logic.full_atomize_tac ctxt' THEN'
  1031      elim_denom_tac ctxt' THEN'
  1032      core_sos_tac print_cert prover ctxt'
  1033   end;
  1034 
  1035 end;