share/contrib/sarag/multiCertificateOfPositivity.mac

   1 /*-------------------------------------
   2
   3 Package : Certificates of positivity in the simplicial Bernstein basis
   4
   5 Author: Richard Leroy
   6
   7 First version : 07/01/07
   8
   9 Updated : 09/01/11
  10
  11 Updated again: 01/02/2017 (fix on standard_mid global m in standard_mid)
  12 --------------------------------------*/
  13
  14 /*-------------------------------------
  15
  16 Procedure: nonneg
  17
  18 Input: List L of the Bernstein expansion of a polynomial P of degree d in k variables
  19
  20 Output: True if all the Bernstein coefficients are nonnegative, else false.
  21
  22 --------------------------------------*/
  23
  24 nonneg(L,d,k):=
  25   block([i,bool],
  26   bool:true,
  27   for i:1 thru length(L) do
  28   (
  29   if L[i]<0 then bool:false
  30   ),
  31   return(bool)
  32   );
  33
  34 /*-------------------------------------
  35
  36 Procedure : pos
  37
  38 Input: List L of the Bernstein expansion of a polynomial P of degree d in k variables
  39
  40 Output: True if the Bernstein coefficients satisfy the certificates of positivity, else false.
  41
  42 Reminder: the certificate of positivy says that all the Bernstein coefficients are nonnegative,
  43 and the coefficients corresponding to the vertices of the simplex are positive.
  44
  45 ---------------------------------------*/
  46
  47 pos(L,d,k):=
  48   block([i,som,bool],
  49   som:sommets(d,k),
  50   bool:true,
  51   for i:1 thru length(L) do
  52   (
  53   if L[i]<0 then bool:false
  54   ),
  55   for i:1 thru length(som) do
  56   (
  57   if L[som[i]]<=0 then bool:false
  58   ),
  59   return(bool)
  60   );
  61
  62 /*-------------------------------------
  63
  64 Procedure: monmin2
  65
  66 Input: list L
  67
  68 Output: the minimum of L, together with a corresponding index.
  69
  70 --------------------------------------*/
  71
  72 monmin2(L):=
  73   block([m,i,j],
  74   m:L[1],
  75   j:1,
  76   for i:2 thru length(L) do
  77   (
  78   if L[i]<m then
  79         (m:L[i],
  80         j:i)
  81   ),
  82   return([m,j])
  83   );
  84
  85 /*-------------------------------------
  86
  87 Procedure: compositions_leq
  88
  89 Input: multiindex gamm
  90
  91 Output: list of all the (weak) compositions of the same size as alpha
  92 that are lower or equal to gamm componentwise (wrt lexicographic order).
  93
  94 --------------------------------------*/
  95
  96 compositions_leq(gamm):=
  97   block([k,l],
  98   k:length(gamm),
  99   if k = 1 then return(makelist([i],i,0,gamm[1])) else
 100         (
 101         l:compositions_leq(rest(gamm,1)),
 102         return(create_list(cons(i,l[j]),i,0,gamm[1],j,1,length(l)))
 103         )
 104   );
 105
 106 /*-------------------------------------
 107
 108 Procedure: compositions
 109
 110 Input: nonnegative integers n,k
 111
 112 Output: list of all the (weak) compositions of n in k nonnegative integers (wrt lexicographic order).
 113
 114 --------------------------------------*/
 115
 116 compositions(n,k):=
 117   block([i,res,l,par],
 118   if k=1 then return([[n]]) else
 119   res:[],
 120     for i : 0 thru n do
 121     (
 122     par:makelist(cons(i,l),l,compositions(n-i,k-1)),
 123     res:append(par,res)
 124     ),
 125   return(res)
 126   );
 127
 128 /*-------------------------------------
 129
 130 Procedure: mulcoeff
 131
 132 Input: a polynomial pol in the variables vars, and a multiindex alpha.
 133
 134 Output: coefficient of pol corresponding to the monomial of degree alpha
 135
 136 --------------------------------------*/
 137
 138 mulcoeff(pol,vars,alpha):=
 139   block([par,res,coe],
 140   if length(vars)=1 then return(coeff(pol,vars[1],alpha[1])) else
 141         (
 142         return(mulcoeff(coeff(pol,vars[1],alpha[1]),rest(vars,1),rest(alpha,1)))
 143         )
 144   );
 145
 146 /*-------------------------------------
 147
 148 Procedure: simplex2bar
 149
 150 Input: a simplex simpl, given as an ordered list of vertices, and a list a variables vars
 151
 152 Output: list of the barycentric coordinates associated to simpl (expressed in the variables vars).
 153
 154 --------------------------------------*/
 155
 156 simplex2bar(simpl,vars):=
 157   block([i,j,len,res,l,eq],
 158   res:[],
 159   len: length(simpl[1]),
 160   for i : 1 thru len do
 161     (
 162     eq:sum(l[j]*simpl[j+1][i],j,0,len)-vars[i],
 163     res:endcons(eq,res)
 164     ),
 165   res:endcons(sum(l[j],j,0,len)-1,res),
 166   res:linsolve(res,makelist(l[j],j,0,len)),
 167   return(makelist(second(res[j]),j,1,len+1))
 168   );
 169
 170 /*-------------------------------------
 171
 172 Procedure: dominant
 173
 174 Input: polynomial P in the variables vars
 175
 176 Output: degree (given as a multiindex) of the dominant coefficient of P wrt lexicographic order on the variables vars
 177
 178 --------------------------------------*/
 179
 180 dominant(P,vars):=
 181   block([i,l,expo,c,d],
 182   l:length(vars),
 183   expo:[],
 184   c:P,
 185   for i : 1 thru l do
 186         (
 187         d:hipow(c,vars[i]),
 188         c:coeff(c,vars[i],d),
 189         expo:endcons(d,expo)
 190         ),
 191   return(expo)
 192   );
 193
 194 /*-------------------------------------
 195
 196 Procedure: sparse_pol
 197
 198 Input: a polynomial P in the variables vars
 199
 200 Output: list of all the monomials that appear in P with a nonzero coefficient
 201
 202 --------------------------------------*/
 203
 204 sparse_pol(P,vars):=
 205   block([i,k,reste,d,mon,c,expo],
 206   k:length(vars),
 207   reste:P,
 208   expo:[],
 209   while not(reste=0) do
 210         (
 211         d:dominant(reste,vars),
 212         c:mulcoeff(reste,vars,d),
 213         reste:reste-c*product(vars[i]**d[i],i,1,k),
 214         expo:endcons(d,expo)
 215         ),
 216   return(expo)
 217   );
 218
 219 /*-------------------------------------
 220
 221 Procedure: sparse_standard2bern
 222
 223 Input: a polynomial P in the variables vars, a nonnegative integer d
 224
 225 Output: Bernstein expansion of degree d of P wrt the standard simplex (uses the previous procedure)
 226
 227 --------------------------------------*/
 228
 229 sparse_standard2bern(P,vars,d):=
 230   block([alph,bet,gamm,i,deg,coeff,res,k,mon,multg,multa,expo],
 231   k:length(vars),
 232   res:[],
 233   expo:sparse_pol(P,vars),
 234   for gamm in compositions(d,k+1) do
 235         (
 236         multg:apply(multinomial_coeff,gamm),
 237         coeff:0,
 238                 for alph in expo do
 239                         (
 240                         bet:gamm-cons(0,alph),
 241                         if apply(min,bet) >= 0 then
 242                                 (
 243                                 deg:sum(alph[a],a,1,length(alph)),
 244                                 multa:apply(multinomial_coeff,bet),
 245                                 coeff:coeff+mulcoeff(P,vars,alph)*multa/multg
 246                                 )
 247                         ),
 248         res:endcons(coeff,res)
 249         ),
 250   return(res)
 251   );
 252
 253 /*-------------------------------------
 254
 255 Procedure: monom2bern
 256
 257 Input: a polynomial P in the variables vars, a nonnegative integer d and a simplex V (given as the list of its vertices)
 258
 259 Output: Bernstein expansion of degree d of P wrt the simplex V (uses the previous procedure)
 260
 261 --------------------------------------*/
 262
 263 monom2bern(P,V,vars,d):=
 264         block([i,j,k,aff,g,varbis,vc],
 265         k:length(vars),
 266         varbis:makelist(concat(vc,i),i,1,k),
 267         aff:transpose(apply(matrix,rest(V,1)))-transpose(apply(matrix,makelist(V[1],i,1,k))),
 268         aff:aff.varbis+V[1],
 269         g:expand(subst(makelist(vars[i]=aff[i,1],i,1,k),P)),
 270         return(sparse_standard2bern(g,varbis,d))
 271         );
 272
 273 /*-------------------------------------
 274
 275 Procedure: bern2monom
 276
 277 Input: the list L of the Bernstein expansion of degree d of a polynomial P (expressed in the variables vars) wrt a simplex V
 278
 279 Output: expression of P in the monomial basis
 280
 281 --------------------------------------*/
 282
 283 bern2monom(L,V,vars,d):=
 284         block([k,i,j,D,P],
 285         k:length(vars),
 286         /*D:concat(makelist(0,i,0,k),makelist(flatten([makelist(0,i,0,j-1),[1],makelist(0,i,j+1,k)]),j,0,k)),*/
 287         P:deCasteljau(V,L,d,vars),
 288         P:expand(P[1][2][1])
 289         );
 290
 291 /*-------------------------------------
 292
 293 Procedure: deCasteljau
 294
 295 Input: the Bernstein expansion ber of degree deg of a polynomial P wrt to a simplex simpl, and an point M.
 296
 297 Output: list of all the items [V,B] where V is a subsimplex of simpl and B the corresponding Bernstein expansion of P,
 298 for all subsimplices V of simpl defined by the point M.
 299
 300 --------------------------------------*/
 301
 302 deCasteljau(simpl,ber,deg,M):=
 303   block([i,j,k,alpha,bet,var,l,cmp,res],local(bcastel),
 304   k:length(simpl)-1,
 305   var:varCasteljau,
 306   bet:subst( makelist( var[j]=M[j],j,1,k), simplex2bar(simpl,var)),
 307   j:1,
 308   for alpha in compositions(deg,k+1) do
 309   (
 310   bcastel[alpha]:ber[j],
 311   j:j+1
 312   ),
 313   for l:1 thru deg do
 314     (
 315     cmp:compositions(deg-l,k+1),
 316     for alpha in cmp do
 317        (
 318        bcastel[alpha]:sum(bet[i+1]*
 319                 bcastel[append(rest(alpha,-length(alpha)+i),[alpha[i+1]+1],
 320                          rest(alpha,i+1))],i,0,k)
 321        )
 322     ),
 323    res:makelist([subst([simpl[j+1]=M],simpl),makelist(bcastel[append(makelist(alpha[i],i,1,j),[0],makelist(alpha[i],i,j+2,k+1))],alpha,compositions(deg,k+1))],j,0,k),
 324   return(res)
 325   );
 326
 327 /*-------------------------------------
 328
 329 Procedure: deCasteljauarete
 330
 331 Input: the Bernstein expansion ber of degree deg of a polynomial P wrt a simplex simpl = [V_0,...,V_k], indices i,j corresponding to two vertices M_i, M_j of simpl, and a real 0<x<1.
 332
 333 Output: list of the Bernstein expansion of degree deg of P wrt the two simplices resulted from the subdivision of simpl at the point M = (1-x)*M_i + x*M_j.
 334
 335 --------------------------------------*/
 336
 337 deCasteljauarete(simpl,ber,deg,i,j,x):=
 338   block([k,alpha,bet,gamm,l,cmp,res,M,p,L,dprime,indi,resi,resj,C],local(bcastel),
 339   k:length(simpl)-1,
 340   M:(1-x)*simpl[i]+x*simpl[j],
 341   if k<2 then return(deCasteljau(simpl,ber,deg,M)),
 342   p:1,
 343   cmp:compositions(deg,k+1),
 344   for alpha in cmp do
 345   (
 346   bcastel[alpha]:ber[p],
 347   p:p+1
 348   ),
 349   cmp:compositions(deg,k),
 350   for l:1 thru deg do
 351     cmp:append(cmp,compositions(deg-l,k-1)),
 352   for bet in cmp do
 353     (
 354     dprime:deg-sum(bet[l],l,1,length(bet)),
 355     L:[],
 356     for gamm:0 thru dprime do
 357       (
 358       indi:append(makelist(bet[p],p,1,i-1),[dprime-gamm],makelist(bet[p],p,i,j-2),[gamm],makelist(bet[p],p,j-1,k-1)),
 359       L:append(L,[bcastel[indi]])
 360       ),
 361     C:deCasteljau([[0],[1]],L,dprime,[x]),
 362     for gamm:0 thru dprime do
 363       (
 364       indi:append(makelist(bet[p],p,1,i-1),[dprime-gamm],makelist(bet[p],p,i,j-2),[gamm],makelist(bet[p],p,j-1,k-1)),
 365       resi[indi]:C[1][2][gamm+1],
 366       resj[indi]:C[2][2][gamm+1]
 367       )
 368     ),
 369   res:[[subst([simpl[i]=M],simpl),makelist(resi[alpha],alpha,compositions(deg,k+1))],[subst([simpl[j]=M],simpl),makelist(resj[alpha],alpha,compositions(deg,k+1))]],
 370   return(res)
 371   );
 372
 373 /*-------------------------------------
 374
 375 Procedure: sommets
 376
 377 Input: two nonnegative integers d and k
 378
 379 Output: list of the multiindices labelling the Bernstein coefficients corresponding to the vertices of a simplex
 380 among the Bernstein expansion of degree d of a polynomial in k variables (wrt the lexicographic order)
 381
 382 --------------------------------------*/
 383
 384 sommets(d,k):=
 385   block([i,j,som,indii],
 386   som:[1],
 387   indii:1,
 388   for j:0 thru k-1 do
 389   (
 390   indii:indii+sum(binomial(d-i+k-1-j,k-1-j),i,1,d-1)+1,
 391   som:endcons(indii,som)
 392   ),
 393   return(som)
 394   );
 395
 396 /*-------------------------------------
 397
 398 Procedure: norm
 399
 400 Input: vector u of IR^k
 401
 402 Output: the square of the Euclidean norm of u: ||u||^2
 403
 404 --------------------------------------*/
 405
 406 norm(u):=
 407   block([i],
 408   return(sum(u[i]^2,i,1,length(u)))
 409   );
 410
 411 /*-------------------------------------
 412
 413 Procedure: diam
 414
 415 Input: simplex V
 416
 417 Output: indices of two vertices of V that realize the diameter of V
 418
 419 --------------------------------------*/
 420
 421 diam(V):=
 422   block([i,j,k,m,L,indi,r,q],
 423   k:length(V)-1,
 424   L:create_list(norm(V[i]-V[j]), i, 1, k+1, j, 1, k+1),
 425   indi:1,
 426   m:L[1],
 427   for i:2 thru length(L) do
 428         (
 429         if L[i]>m then
 430         (
 431         indi:i,
 432         m:L[i]
 433         )
 434         ),
 435   r:mod(indi,k+1),
 436   q:quotient(indi,k+1),
 437   if r = 0 then return([min(k+1,q),max(k+1,q)]) else return([min(r,q+1),max(r,q+1)])
 438   );
 439
 440 /*-------------------------------------
 441
 442 Procedure: bound
 443
 444 Input: a polynomial P of degree d in the variables var
 445
 446 Output: theoritical lower bound on the minimum of P over the standard simplex, if P is positive on it.
 447
 448 --------------------------------------*/
 449
 450 bound(P,var,d):=
 451         block([i,j,k,c,tau],
 452         k:length(var),
 453         c:coeffs(P,var,d),
 454         tau:length(num2list(max(apply(max,c),-apply(min,c)),2)),
 455         return(2^(-(tau+1)*d^(k+1))*d^(-(k+1)*d^k)*binomial(d+k,k+1)^(-d^k*(d-1)))
 456         );
 457
 458 /*-------------------------------------
 459
 460 Procedure: coeffs
 461
 462 Input: a polynomial P of degree D in variables var
 463
 464 Output: list of all the coefficients of P (wrt the lexicographic order on the variables).
 465
 466 --------------------------------------*/
 467
 468 coeffs(P,var,D):=
 469         block([i,j,a,k,v],
 470         k:length(var),
 471         apply(append,makelist(makelist(subst(makelist(v=0,v,var),ratcoeff(P,product(var[i]^(a[i]),i,1,k))),a,compositions(D-j,k)),j,0,D)
 472         ));
 473
 474 /*-------------------------------------
 475
 476 Procedure: multiCertificate
 477
 478 Input: a polynomial P in the variables vars, a simplex V, a nonnegative integer d, a subdivision method sub (e.g. bisection, bisection_castel, standard_triang2),
 479 and a choice for the certificate (pos for a certificate of positivity, nonneg for a certificate of nonnegativity).
 480
 481 Output:
 482
 483 - if P is positive on the standard simplex :a certificate of positivity for P expressed as a list of items [V,B],
 484 where the collection of subsimplices V covers the standard simplex,
 485 and B is the corresponding Berstein expansion of degree d of P on each W
 486
 487 - if not, an item [M,v] where M is a point of the standard simplex and v is a the value of P at M that is lower than the theoritical bound
 488 on the minimum of a positive polynomial
 489
 490 Warning: if P is nonnegative but not positive, the computation can take a very long time, due to the theoritical bound, which can be very small.
 491
 492 --------------------------------------*/
 493
 494 multiCertificate(P,V,vars,d,sub,cert):=
 495   block([k,bool,actif,passif,alpha,a,v,arete,w,Delta,m,N0,L,s,res,value,value_new,N,W,l],
 496   k:length(V)-1,
 497   res:[],
 498   L:monom2bern(P,V,vars,d),
 499   s:sommets(d,k),
 500   actif:[[V,L]],
 501   passif:[],
 502   if cert=pos then m:bound(P,vars,d) else m:0,
 503   value:[0,m+1],
 504   while length(actif)>0 and value[2]>=m do
 505   (
 506   a:[],
 507   for v in actif do
 508         (
 509         if apply(cert,[v[2],d,k]) then passif:endcons(v,passif)
 510         else
 511         (
 512         for i:1 thru k+1 do
 513                 (
 514                 l:v[2][s[i]],
 515                 if not(apply(cert,[[l],0,0])) then
 516                         (
 517                         res:([v[1][i],l]),
 518                         return()
 519                         )
 520                 ),
 521         if length(res)>0 then
 522                 (
 523                 return()
 524                 ),
 525         W:apply(sub,[v[1],v[2],d]),
 526
 527         for w in W do
 528                 (
 529                 alpha:monmin2(w[2])[2],
 530                 alpha:compositions(d,k+1)[alpha],
 531                 N:sum(alpha[i]*w[1][i],i,1,k+1)/d,
 532                 value_new:subst(makelist(vars[i]=N[i],i,1,k),P),
 533                 if value_new<value[2] then value:[N,value_new]
 534                 ),
 535
 536         a:append(a,W)
 537         )
 538         ),
 539         actif:a
 540  ),
 541  if length(res)>0 then
 542     return(res),
 543  if length(actif)=0 then
 544     return(passif)
 545  else
 546     return(value)
 547  );
 548
 549 /*-------------------------------------
 550
 551 Procedure: bisection_castel
 552
 553 Input: the Bernstein expansion L of degree d of a polynomial P wrt a simplex V
 554
 555 Output: list of the items [W,B] where B is the Bernstein expansion of degree d of P wrt the simplex W,
 556 for all the subsimplices W of V generated by the bisection of V at the control point corresponding to
 557 the multiindex associated to a minimal Berstein coefficient among L.
 558
 559 --------------------------------------*/
 560
 561 bisection_castel(V,L,d):=
 562         block([arete,w,k,c,j,alpha,N,W,w,res,C],
 563         res:[],
 564         k:length(V)-1,
 565         c:compositions(d,k+1),
 566         i:monmin2(L)[2],
 567         alpha:c[i],
 568         j:[],
 569         for l:1 thru length(alpha) do
 570                 (
 571                 if not(alpha[l]=0) then j:endcons([l,alpha[l]],j)
 572                 ),
 573         if length(j)=2 then return(deCasteljauarete(V,L,d,j[1][1],j[2][1],j[1][2]/d))
 574         else
 575         (
 576         N:sum(alpha[i]*V[i],i,1,k+1)/d,
 577         return(deCasteljau(V,L,d,N))
 578         )
 579         );
 580
 581 /*-------------------------------------
 582
 583 Procedure: bisection
 584
 585 Input: the Bernstein expansion L of degree d of a polynomial P wrt a simplex V
 586
 587 Output: list of the items [W,B] where B is the Bernstein expansion of degree d of P wrt the simplex W,
 588 for the two subsimplices W of V generated by the bisection of V at the midpoint of a diameter of V.
 589
 590 --------------------------------------*/
 591
 592 bisection(V,L,d):=
 593         block([arete],
 594         arete:diam(V),
 595         return(deCasteljauarete(V,L,d,min(arete[1],arete[2]),max(arete[1],arete[2]),1/2))
 596         );
 597
 598 /*-------------------------------------
 599
 600 Procedure: bisection_max
 601
 602 Input: the Bernstein expansion L of degree d of a polynomial P wrt a simplex V
 603
 604 Output: list of the items [W,B] where B is the Bernstein expansion of degree d of P wrt the simplex W,
 605 for the two subsimplices W of V generated by the bisection of V at a point of a diameter of V
 606 "maximizing" the volume of one subsimplex satisfying the certificate of positivity (cert= pos) or nonnegativity (cert=nonneg)
 607
 608 --------------------------------------*/
 609
 610 bisection_max(V,L,d,cert):=
 611         block([arete,W,i,n,k,bool,res],
 612         arete:diam(V),
 613         i:1,
 614         n:20,
 615         k:length(V)-1,
 616         res:[V,L],
 617         bool:apply(cert,[L,d,k]),
 618         while ((not bool) and i<n/2) do
 619                 (
 620                 W:deCasteljauarete(V,L,d,min(arete[1],arete[2]),max(arete[1],arete[2]),i/n),
 621                 if (apply(cert,[W[1][2],d,k])) then
 622                         (bool:true,
 623                         res:W
 624                         )
 625                 else
 626                         (
 627                         W:deCasteljauarete(V,L,d,min(arete[1],arete[2]),max(arete[1],arete[2]),(n-i)/n),
 628                         if (apply(cert,[W[2][2],d,k])) then
 629                                 (bool:true,
 630                                 res:W
 631                                 )
 632                         ),
 633                 i:i+1
 634                 ),
 635         if i>=n/2 then return(deCasteljauarete(V,L,d,min(arete[1],arete[2]),max(arete[1],arete[2]),1/2)) else return(res)
 636         );
 637
 638 bisection_max_pos(V,L,d):=
 639         block(
 640         return(bisection_max(V,L,d,pos))
 641         );
 642
 643
 644 bisection_max_nonneg(V,L,d):=
 645         block(
 646         return(bisection_max(V,L,d,nonneg))
 647         );
 648
 649
 650 /*-------------------------------------
 651
 652 Procedure: num2list
 653
 654 Input: two nonnegative integers n and b
 655
 656 Output: the expression of n in base b.
 657
 658 --------------------------------------*/
 659
 660 num2list(n,b):=block([lst:[],a:n],
 661  for i:1 while a#0 do (
 662    lst:cons(?mod(a,b),lst),
 663    a:(a-lst[1])/b
 664    ),
 665  return(lst)
 666  );
 667
 668 /*-------------------------------------
 669
 670 Procedure: standard_triang2
 671
 672 Input: list L of the Bernstein expansion of degree d of a polynomial P wrt to a simplex V
 673
 674 Output: list of all the items [W,B] where W is a subsimplex of V and B the corresponding Bernstein expansion of degree d of P,
 675 for all the subsimplices W of the standard triangulation of degree 2 of V.
 676
 677 --------------------------------------*/
 678
 679 standard_mid(k,d):=
 680         block([M,m:make_array('any,d,k+1),n,y,i,j,color,t,s],
 681         t:[],
 682         s:matrix(makelist(1/d,i,1,d)),
 683         for n:0 thru d^k-1 do
 684         (
 685         y:num2list(n,d),
 686         y:append(makelist(0,i,1,k-length(y)),y),
 687         color:0,
 688         for i:0 thru d-1 do
 689                 (
 690                 m[i,0]:color,
 691                 for j:1 thru k do
 692                         (
 693                         (if y[j]=i then color:color+1),
 694                         m[i,j]:color
 695                         )
 696                 ),
 697         t:cons(apply(matrix,makelist(makelist(m[i,j],j,0,k),i,0,d-1)),t)
 698         ),
 699         return(t)
 700         );
 701
 702 standard(V,d):=
 703         block([t,k,T,i,j,l],
 704         k:length(V)-1,
 705         T:standard_mid(k,d),
 706         return(makelist(makelist(sum(V[t[i,j]+1],i,1,d)/2,j,1,k+1),t,T))
 707         );
 708
 709 standard_triang(V,P,vars,d,l):=
 710         block([i,j,k,T,L,res,t],
 711         k:length(V)-1,
 712         res:[],
 713         T:standard(V,l),
 714         for t in T do
 715                 (
 716                 i:monom2bern(P,t,vars,d),
 717                 res:cons([t,i],res)
 718                 ),
 719         return(res)
 720         );
 721
 722 standard_triang2(V,L,d):=
 723         block([P,k,varbis,vc],
 724         k:length(V)-1,
 725         varbis::makelist(concat(vc,i),i,1,k),
 726         P:bern2monom(L,V,varbis,d),
 727         return(standard_triang(V,P,varbis,d,2))
 728         );
 729
 730 /*-------------------------------------
 731
 732 Procedure: sec_diff
 733
 734 Input: Bernstein expansion B of degree d of a polynomial P wrt a simplex V
 735
 736 Output: List of all the second differences of P.
 737
 738 --------------------------------------*/
 739
 740 sec_diff(V,B,d):=
 741         block([i,j,k,alpha,bet,gam,var,l,cmp,res,e],local(bcastel),
 742         res:[],
 743         k:length(V)-1,
 744         e:makelist(0,i,0,k),
 745         j:1,
 746         for i:0 thru k do
 747         (
 748         e[i+1]:append(makelist(0,j,0,i-1),[1],makelist(0,j,i+1,k))
 749         ),
 750         for alpha in compositions(d,k+1) do
 751         (
 752         bcastel[alpha]:B[j],
 753         j:j+1
 754         ),
 755         for gam in compositions(d-2,k+1) do
 756         (
 757         for j:1 thru k do
 758                         (
 759                         res:endcons(bcastel[gam+e[1]+e[j]]+bcastel[gam+e[k+1]+e[j+1]]-bcastel[gam+e[k+1]+e[j]]-bcastel[gam+e[1]+e[j+1]],res)
 760                         ),
 761                 for i:1 thru k do
 762                 (
 763                         for j:i+1 thru k do
 764                         (
 765                         res:endcons(bcastel[gam+e[i+1]+e[j]]+bcastel[gam+e[i]+e[j+1]]-bcastel[gam+e[i]+e[j]]-bcastel[gam+e[i+1]+e[j+1]],res)
 766                         )
 767                 )
 768         ),
 769         res:max(apply(max,res),-apply(min,res)),
 770         return(res)
 771         );
 772
 773 /*-------------------------------------
 774
 775 Procedure: random_color
 776
 777 Input: none
 778
 779 Output: code of a random color
 780
 781 --------------------------------------*/
 782
 783 random_color():=
 784   block([randcol],
 785   randcol:['1,'2,'3,'4,'5,'6,'7,'8,'9,'0,'A,'B,'C,'D,'E,'F],
 786     concat("#",randcol[random(16)+1],randcol[random(16)+1],randcol[random(16)+1],randcol[random(16)+1],randcol[random(16)+1],randcol[random(16)+1])
 787   );
 788
 789 /*-------------------------------------
 790
 791 Procedure: drawMultiCerticificate
 792
 793 Input: Output L of multiCertificate in dimension 2 when a polynomial P has been proven to be positive, and an integer s
 794
 795 Output: Vizualization of the subsimplices generated by multiCertificate.
 796 The size of the image is defined by s (usually 500).
 797
 798 --------------------------------------*/
 799
 800 load(draw);
 801
 802 drawMultiCertificate(L,s):=
 803   block([gone,t],
 804   t:concat(length(L)," subsimplices"),
 805   if length(L[1]) = 2 then
 806   (
 807   gone:append([dimensions=[s,s],xtics=false,ytics=false,axis_bottom=false,axis_left=false,axis_top=false,axis_right=false,colorbox=false,title=t],flatten(makelist([fill_color=random_color(),polygon(L[i][1])],i,1,length(L))))
 808   )
 809   else
 810   (
 811   gone:append([dimensions=[s,s],xtics=false,ytics=false,axis_bottom=false,axis_left=false,axis_top=false,axis_right=false,colorbox=false,fill_color=red, title=t],flatten(makelist([polygon(L[i])],i,1,length(L))))
 812   ),
 813   apply(draw2d,gone)
 814   );
 815
 816 /*-------------------------------------
 817
 818 Procedure: wxdrawMultiCerticificate
 819
 820 Same procedure as drawMultiCertificate, except that the vizualization is generated inside wxMaxima
 821
 822 --------------------------------------*/
 823
 824 wxdrawMultiCertificate(L,s):=
 825   block([gone,t],
 826   t:concat(length(L)," subsimplices"),
 827   if length(L[1]) = 2 then
 828   (
 829   gone:append([dimensions=[s,s],xtics=false,ytics=false,axis_bottom=false,axis_left=false,axis_top=false,axis_right=false,colorbox=false,color=black],flatten(makelist([fill_color=white,polygon(L[i][1])],i,1,length(L))))
 830   )
 831   else
 832   (
 833   gone:append([dimensions=[s,s],xtics=false,ytics=false,axis_bottom=false,axis_left=false,axis_top=false,axis_right=false,colorbox=false,fill_color=white,color=black],flatten(makelist([polygon(L[i])],i,1,length(L))))
 834   ),
 835   apply(wxdraw2d,gone),
 836   return(t)
 837   );
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856