share/calculus/qual.mac

   1 ttyoff:true $
   2 load("qualsp");
   3 matchdeclare([utrue,vtrue,wtrue],true)$
   4 tellsimp(qual(utrue), qual1(utrue, listofvars(utrue))) $
   5 tellsimp(qual(utrue,vtrue),qual1(utrue,qual_listify(vtrue))) $
   6 tellsimp(revelation(utrue), revelation1(utrue,200,300)) $
   7 tellsimp(revelation(utrue,vtrue), revelation1(utrue,vtrue,300))$
   8 tellsimp(revelation(utrue,vtrue,wtrue),
   9    revelation1(utrue,vtrue,wtrue)) $
  10 tellsimp(slopes(utrue),slopes1(utrue,listofvars(utrue)))$
  11 tellsimp(slopes(utrue,vtrue),slopes1(utrue,qual_listify(vtrue)))$
  12 tellsimp(symmetry(utrue),symmetry1(utrue,listofvars(utrue)))$
  13 tellsimp(symmetry(utrue,vtrue),symmetry1(utrue,qual_listify(vtrue)))$
  14 tellsimp(periods(utrue), periods1(utrue,listofvars(utrue))) $
  15 tellsimp(periods(utrue,vtrue),periods1(utrue,qual_listify(vtrue))) $
  16 tellsimp(limits(utrue),limits1(utrue,listofvars(utrue)))$
  17 tellsimp(limits(utrue,vtrue),limits1(utrue,qual_listify(vtrue)))$
  18 tellsimp(stationarypoints(utrue),stationarypoints1(utrue,
  19    listofvars(utrue)))$
  20 tellsimp(stationarypoints(utrue,vtrue),stationarypoints1(utrue,
  21    qual_listify(vtrue))) $
  22
  23 variablep(u) := is(atom(u) and not numberp(u) or subvarp(u)) $
  24
  25 qual_listify(u) :=
  26    if listp(u) then u else [u] $
  27
  28 qual1(u,v) := block(
  29    revelation1(u, 200, 300),
  30    return([first(ldisp('bounds=bounds(u))), slopes1(u,v),
  31       ldisp('curvature=curvature(u)), symmetry1(u:radcan(u),v),
  32       periods1(u,v), zerosandsingularities(u), limits1(u,v),
  33       stationarypoints1(u,v)])) $
  34
  35 revelation1(u,umin,revmax) := block(
  36    [rev, lold, lnew, lu],
  37    if (lu:?length(?makstring(u)))>umin then (lold:-1,
  38       for j:1 step 1 while (lnew:?length(?makstring(rev:reveal(u,j))))
  39             <=revmax and lnew#lold and lnew<lu do(
  40          disp('reveal("...", ''j) = rev),
  41          lold:lnew))) $
  42
  43 slopes1(u,v) := block(
  44    [ans, partswitch, prederror],
  45    partswitch:true,  prederror:false,  ans: [],
  46    for x in v do ans: cons(slopes2(u,x), ans),
  47    return(ans)) $
  48
  49 slopes2(u,x) := block(
  50    u: bounds1(diff(u,x)),
  51    return(first(ldisp(slope(x) =
  52       if posl(u[1]) then 'increasing
  53       else if negu(u[2]) then 'decreasing
  54       else if nonnegl(u[1]) then
  55          if nonposu(u[2]) then 'constant
  56          else 'nondecreasing
  57       else if nonposu(u[2]) then 'nonincreasing
  58       else 'unknown)))) $
  59
  60 curvature(u) := block(
  61    [v], v:listofvars(u),
  62    return(['strictconcave, 'concave, 'nonconvex, 'concaveandconvex,
  63       'nonconcave, 'convex, 'strictconvex,
  64       'neitherconcavenorconvex, 'unknown]
  65       [definitecode(qual_hessian(gradient(u,v),v))])) $
  66
  67 qual_hessian(g,v) := block(
  68    [ans],
  69    ans:[],
  70    for x in v do ans: endcons(diff(g,x), ans),
  71    funmake('matrix, ans))$
  72
  73 gradient(u,v) := block(
  74    [ans],
  75    ans: [],
  76    for x in v do ans: endcons(diff(u,x), ans),
  77    return(ans)) $
  78
  79
  80 symmetry1(u,v) := block(
  81    [ans],
  82    ans: [],
  83    if u=0 then return(['zero]),
  84    for x in v do ans: endcons(first(ldisp(symmetries(x)=symmetry2(u,
  85      x))), ans),
  86    return(ans)) $
  87
  88 symmetry2(u,x):= block(
  89    [umx, evn, od, temp, v],
  90    umx: subst(x=-x, u),
  91    temp: radcan(u-umx),
  92    if temp=0 then return('even),
  93    umx: radcan(u+umx),
  94    if umx=0 then return('odd),
  95    if numberp(temp) then evn:'no
  96    else if length(v:listofvars(umx))=1 then evn:zeroequiv(temp,v)
  97    else evn: 'unknown,
  98    if numberp(umx) then od: 'no
  99    else if length(v:listofvars(umx))=1 then od:zeroequiv(temp,v)
 100    else od: 'unknown,
 101    if evn=true then
 102       if od=true then
 103          if zeroequiv(u,v)=true then return('probablyzero)
 104          else return('unknown)
 105       else return('probablyeven),
 106    if od=true then return('probablyodd),
 107    if evn='no then
 108       if od='no then return('neither)
 109       else if od=false then return('nonevenandprobablynonodd)
 110       else return('noneven),
 111    if od='no then
 112       if evn=false then return('nonoddandprobablynoneven)
 113       else return('nonodd),
 114    if evn=false then
 115       if od=false then return('probablyneither)
 116       else return('probablynoneven),
 117    if od=false then return('probablynonodd),
 118    return('unknown)) $
 119
 120 periods1(u,v) := block(
 121    [ans, partswitch],
 122    partswitch: true,
 123    u: trigreduce(u),
 124    ans: [],
 125    for x in v do ans: endcons(first(ldisp(period(x)=period2(u,
 126       x))),ans),
 127    return(ans)) $
 128
 129 period2(u,x) := block(
 130    [ans],
 131    if numberp(u) then return(0),
 132    if variablep(u) then
 133       if u=x then return('inf)
 134       else return(0),
 135    if inpart(u,0)="*" or piece="+" then (
 136       ans: period2(inpart(u,1), x),
 137       for j:2 step 1 while ans # 'inf and inpart(u,j) # 'end do
 138          ans: lcmspec(ans,period2(piece,x)),
 139       return(ans)),
 140    if piece="^" then return(lcmspec(
 141       period2(inpart(u,1),x),period2(inpart(u,2),x))),
 142    if piece='sin or piece='cos or piece='sec or piece='csc then
 143       if freeof(x,inpart(u,1)) then return(0)
 144       else if freeof(x,ans:diff(piece,x)) then return(2*%pi/ans)
 145       else return('inf),
 146    if piece='tan or piece='cot then
 147       if freeof(x,inpart(u,1)) then return(0)
 148       else if freeof(x,ans:diff(piece,x)) then return(%pi/ans)
 149       else return('inf),
 150    return(period2(inpart(u,1),x))) $
 151
 152 lcmspec(u,v) :=
 153    if u=0 then v
 154    else if v=0 then u
 155    else if u='inf or v='inf then 'inf
 156    else num(u)*num(v)/gcd(num(u)*denom(v), num(v)*denom(u)) $
 157
 158
 159 limits1(u,v) := block(
 160    [ans, t, partswitch],
 161    ans: [],
 162    partswitch: true,
 163    for x in v do (t: lbatom(x),
 164       ans : endcons(first(ldisp(limitas(x,t) =
 165          if inpart(t,0)='strict then strict(limit(u,x,inpart(t,1),
 166             'plus))
 167          else limit(u,x,t,'plus))), ans),
 168       t: ubatom(x),
 169       ans: endcons(first(ldisp(limitas(x,t) =
 170          if inpart(t,0)='strict then strict(limit(u,x,inpart(t,1),
 171             'minus))
 172          else limit(u,x,t,'minus))), ans)),
 173    return(ans)) $
 174
 175 zerosandsingularities(u) := block(
 176    [partswitch, temp, prederror],
 177    prederror: false,
 178    partswitch:true,
 179    u: radcan(trigreduce(u)),
 180    temp: zp1(factor(ratdenom(u)), zp1(factor(ratnumer(u)),[[],[]])),
 181    return(ldisp('zeros = first(temp), 'singularities=temp[2])))$
 182
 183 zp1(n,zp) := block(
 184    [z,p],
 185    z:first(zp), p:zp[2],
 186    if not constantp(n) then
 187       if inpart(n,0)="*" then for j:1 step 1 while inpart(n,j)#'end
 188          do (if not constantp(piece) then (z:cons(piece=0,z),
 189             p:conssingularities(p,piece)))
 190       else(z:cons(n=0,z),
 191          p:conssingularities(p,n)),
 192    return([p,z])) $
 193 conssingularities(p,u) := block(
 194    [bas],
 195    if variablep(u) then return(p),
 196    if inpart(u,0)="+" or piece="*" then
 197       for j:1 step 1 while inpart(u,j)#'end do p:conssingularities(p,piece)
 198    else if piece="^" and not constantp(bas:inpart(u,1)) then
 199       if numberp(piece) then (
 200          if piece<0 then p:cons(bas=0, p))
 201       else piece: cons(bas=0 and piece<0, p)
 202    else if piece='log and not numberp(inpart(u,1)) then
 203       p:cons(piece=0,p)
 204    else if (piece='tan or piece='sec) and not numberp(inpart(u,1))
 205       then p: cons(piece-('integer+1/2)*%pi=0, p)
 206    else if (piece='cot or piece='csc) and not numberp(inpart(u,1))
 207       then p: cons(piece-'integer*%pi=0, p)
 208    else if piece='atanh and not numberp(inpart(u,1))
 209       then p: cons(piece-1=0, cons(piece+1=0, p)),
 210    return(p)) $
 211
 212 stationarypoints1(u,v) := block(
 213    [singsolve,grindswitch,dispflag,g,ans,uu,s],
 214    g:gradient(u,v),
 215    singsolve: grindswitch:  true,
 216    dispflag: false,
 217    s:errcatch(ev(solve(g,v),eval)),
 218    if s=[] or s=[[]] or s=[[false=0]]
 219       then return(ldisp("no stationary points found")),
 220    s:first(s),
 221    ans: ldisp("stationary points" = s),
 222    uu:[],
 223    for ss in s do uu: endcons(if length(v)>1 or first(v)=lhs(ss) and
 224       freeof(first(v),rhs(ss)) then subst(ss,u) else 'unknown ,uu),
 225    ans:endcons(first(ldisp("corresponding expression values" = uu)),
 226    ans),
 227    g: qual_hessian(g,v), uu: [],
 228    for ss in s do uu: endcons(type(definitecode(subst(ss,g))),uu),
 229    ans: endcons(first(ldisp("corresponding types" = uu)), ans),
 230    return(ans)) $
 231
 232 type(u) :=
 233    ['maximum, 'nonminimum, 'nonminimum, 'unknown, 'nonmaximum,
 234       'nonmaximum, 'minimum, 'saddle, 'unknown][u] $
 235
 236 bounds(w) := ev(bounds1(w),prederror:false,partswitch:true)$
 237
 238 bounds1(w) := block(/* W is an expression.  Returns list of its
 239       lower, then upper bounds.  (reference: file qual usage .  In
 240       comments below, "symbolic" means neither numerical, inf, minf,
 241       or strict with such an argument. */
 242   [u, v, t],
 243   if numberp(w) then return([w,w]),
 244   if variablep(w) then return([lbatom(w), ubatom(w)]),
 245   if inpart(w,0) = "+" then (u: bounds1(inpart(w,1)),
 246     for j:2 step 1 while u#['minf,'inf] and inpart(w,j) # 'end do
 247       (v: bounds1(piece),
 248       u: [addbnd(u[1],v[1]), addbnd(u[2],v[2])]),
 249     return(u)),
 250
 251   if piece = "*" then (u:bounds1(inpart(w,1)),
 252     for j:2 step 1 while inpart(w,j) # 'end do (
 253       v:bounds1(piece),
 254         /* Try standardizing lowerbound of 1st arg to nonnegative: */
 255       if nonnegl(u[1]) then u:bndnntimes(u,v)
 256       else if nonnegl(v[1]) then u:bndnntimes(v,u)
 257       else if nonposu(u[2]) then u:bndnntimes(bndminus(u),bndminus(v))
 258       else if nonposu(v[2]) then u:bndnntimes(bndminus(v),bndminus(u))
 259         /* Try standardizing lowerbound of 1st arg to negative: */
 260       else if negl(u[1]) then u:bndnegtimes(u,v)
 261       else if negl(v[1]) then u:bndnegtimes(v,u)
 262       else if posu(u[2]) then u:bndnegtimes(bndminus(u),bndminus(v))
 263       else if posu(v[2]) then u:bndnegtimes(bndminus(v),bndminus(u))
 264         /* Both bounds of both args are symbolic: */
 265       else (u:[u[1]*v[1], u[1]*v[2], u[2]*v[1], u[2]*v[2]],
 266         u: [apply('min,u), apply('max,u)])),
 267       return(u)),
 268
 269   if piece="^" then (u:bounds1(inpart(w,1)), v:bounds1(inpart(w,2)),
 270     if posl(u[1]) then
 271         /*Try standardizing lowerbound of 1st arg to >=1: */
 272       if ge1l(u[1]) then return(bndge1to(u,v))
 273       else if le1u(u[2]) then return(bndrecip(bndge1to(bndrecip(u),
 274         v)))
 275       else if ge1u(u[2]) then
 276           /* 0<=u[1]<1 and u[2]>1.  Try standardizing
 277              lower bound of 2nd arg to nonnegative: */
 278         if nonnegl(v[1]) then return(bndspan1tonn(u,v))
 279         else if nonposu(v[2]) then return(bndrecip(bndspan1tonn(u,
 280           bndminus(v))))
 281             /* v[1]<1 or symbolic & v[2]>1 or symbolic.  Standardize
 282                nonsymbolic args of ** to nonneg: */
 283         else return([min(nntonn(u[1],v[2]),recipl(nntonn(u[2],neg8(
 284           v[1])))), max(nntonn(u[2],v[2]),recipu(nntonn(u[1],neg8(
 285           v[1]))))])
 286             /* 0<=u[1]<1 & u[2] symbolic.  Try standardizing lower
 287                bound of 2nd arg to nonegative: */
 288       else if nonnegl(v[1]) then return(bndmayspan1tonn(u,v))
 289       else if nonposu(v[2]) then
 290         return(bndrecip(bndmayspan1tonn(u,bndminus(v))))
 291           /* u[1]<1 & u[2] symbolic: */
 292       else if posu(v[2]) then
 293         if negl(v[1]) then return([min(nntonn(u[1],v[2]),u[2]**v[1]),
 294           max(recipu(nntonn(u[1],neg8(v[1]))), u[2]**v[2])])
 295             /* v[1] symbolic too, so another possible upperbound:*/
 296         else return([min(nntonn(u[1],v[2]), u[2]**v[1]),
 297           max(u[1]**v[1], u[2]**v[2], u[2]**v[1])])
 298       else if negl(v[1]) then return([min(u[1]**v[2],u[2]**v[2],u[2]
 299         **v[1]),max(recipu(nntonn(u[1],neg8(v[1]))),u[2]**v[2])])
 300           /* v[1] & v[2] symbolic.  3 symbolic possibilities for
 301              each bound: */
 302       else return([min(u[1]**v[2], u[2]**v[2], u[2]**v[1]),
 303                    max(u[1]**v[1], u[2]**v[2], u[2]**v[1])])
 304       /* u[1]=0 or symbolic.  Negatives must not be raised to
 305          nonintegers: */
 306     else if integerp(v[1]) and integerp(v[2]) then
 307       if v[1]=v[2] then  /* interval ** integer: */
 308         if v[1]>=0 then
 309           if evnp(v[1]) then
 310             if nonposu(u[2]) then return([nntonn(neg8(u[2]),v[1]),
 311               nntonn(neg8(u[1]),v[1])])
 312                 /* interval spanning 0 ** nonnegative integer: */
 313             else if negl(u[1]) and posu(u[2]) then return([0,
 314               max(nntonn(u[2],v[1]), nntonn(neg8(u[1]),v[1]))])
 315                 /* u[1] or u[2] symbolic so that maybespan0 **
 316                    nonnegative even integer: */
 317             else return([if posu(u[2]) then 0 else u[2]**v[2],
 318                 max(nntonn(neg8(u[1]),v[2]), u[2]**v[2])])
 319           else return([neg8(nntonn(neg8(u[1]),v[1])),
 320                    /* Allow for symbolic or either-signed
 321                        upper bound of u: */
 322                  if negu(u[2]) then neg8(nntonn(neg8(u[2]),v[1]))
 323                  else nntonn(u[2],v[1])])
 324           /* u[1]<0: */
 325         else if nonposu(u[2]) then
 326           if evnp(v[1]) then return(bndrecip(bndge1tonn(bndminus(u),
 327             bndminus(v))))
 328           else return(bndminus(bndrecip(bndge1tonn(bndminus(u),
 329             bndminus(v)))))
 330         else return(['minf,'inf])
 331       else if negu(u[2]) then
 332           /* Try standardizing  lowerbound of 1st arg <=-1: */
 333         if lem1u(u[2]) then return(bndlem1to(u,v))
 334         else if gem1l(u[1]) then return(bndlem1to(bndrecip(u)
 335           ,bndminus(v)))
 336         else if lem1l(u[1]) then
 337             /* u[1]<-1 & u[2]>-1.  Try standardizing lower
 338               bound of v to nonnegative: */
 339           if nonnegl(v[1]) then return(bndspanm1tonn(u,v))
 340           else if nonposu(v[2]) then return(bndrecip(bndspanm1tonn(
 341             u, bndminus(v))))
 342           else (w: bndlem1tonn(u,v),
 343               /* v[1]<0 or symbolic & v[2]>0 or symbolic: */
 344             u: bndlem1tonn(bndrecip(u),bndminus(v)),
 345             return([min(u[1],w[1]), max(u[2],w[2])]))
 346               /* u[1] algebraic: */
 347         else return([lb(w),ub(w)])
 348       else if v[1]>=0 then (     /* 0<=v[1]<v[2]: */
 349         if lem1l(u[1]) then t: bndlem1tonn(u,v)
 350         else if gem1l(u[1]) then
 351             /* u[1] symbolic: */
 352           t: bndrecip(bndlem1tonn(bndrecip([1,u[1]]),v))
 353         else return([lb(w), ub(w)]),
 354         if ge1u(u[2]) then u: nntonn(u[2],v[2])
 355         else if le1u(u[2]) then u: nntonn(u[2],v[1])
 356           /* u[2] symbolic: */
 357         else return([lb(w), ub(w)]),
 358         return([t[1], max(t[2],u)]))
 359     else if v[2]<0 and negu(u[2]) and posl(u[1])
 360       then return(['minf,'inf])
 361     else return(['minf, 'inf])),
 362
 363   if piece='log or piece='atan or piece='erf or piece='sinh or
 364     piece='asinh or piece='acosh or
 365     piece='tanh then return(bndunary(piece, bounds(inpart(w,1)))),
 366   if piece = 'sin or piece = 'cos then return([-1,1]),
 367   if piece='acot or piece='asech then return(
 368     reverse(bndunary(piece, bounds(inpart(w,1))))),
 369   if piece = 'cosh then return([1, 'inf]),
 370   if piece='sech then return([0,1]),
 371   if piece='asec then return([0, 3.14159]),
 372   if piece='acsc then return([-1.57079, 1.57079]),
 373   if piece='asin or piece='atanh then return(bndrestrict(piece,w)),
 374   if piece='acos then return(reverse(bndrestrict(piece,w))),
 375   return(['minf, 'inf])) $
 376
 377 bndrestrict(p,w) := block(
 378    w:bounds(inpart(w,1)),
 379    if lem1l(w[1]) then w[1]:-1,
 380    if ge1u(w[2]) then w[2]:1,
 381    return(bndunary(p,w))) $
 382
 383 addbnd(b1,b2) := /* b1 and b2 are both lower or both upper
 384       bounds.  Returns their sum.  Assumes partswitch:true. */
 385    if b1='inf or b2='inf then 'inf
 386    else if b1='minf or b2='minf then 'minf
 387    else if inpart(b1,0)='strict then
 388       if inpart(b2,0)='strict then
 389          strict(addbnd(inpart(b1,1), inpart(b2,1)))
 390       else strict(addbnd(inpart(b1,1), b2))
 391    else if inpart(b2,0)='strict then strict(addbnd(b1,inpart(b2,1)))
 392    else b1+b2 $
 393
 394 bndge1to(u,v) := /* u & v are intervals, with u[1]>=1.  Returns
 395       interval of u**v.  First try standardizing to nonnegative
 396       lower bound of power: */
 397    if nonnegl(v[1]) then bndge1tonn(u,v)
 398    else if nonposu(v[2]) then bndrecip(bndge1tonn(u,bndminus(v)))
 399    else if negl(v[1]) then
 400       if posu(v[2]) then [recipl(nntonn(u[2], neg8(v[1]))),
 401          nntonn(u[2],v[2])]
 402         /* v[2] symbolic: */
 403       else [recipl(nntonn(u[2], neg8(v[1]))), max(u[1]**v[2],
 404          u[2]**v[2])]
 405       /* v[1] symbolic: */
 406    else if posu(v[2]) then [min(u[1]**v[1], u[2]**v[1]),
 407       nntonn(u[2],v[2])]
 408       /* v[1] and v[2] symbolic: */
 409    else [min(u[1]**v[1], u[2]**v[1]), max(u[1]**v[2], u[2]**v[2])] $
 410
 411 bndge1tonn(u,v) := /* u & v are intervals with u[2]>=1, v[1]>=0.
 412       Returns interval of u**v. */
 413    [nntonn(u[1],v[1]), nntonn(u[2],v[2])] $
 414
 415 bndlem1to(u,v) := /* u and v are intervals with u[2]<=-1 &
 416       v[1] & v[2] are unequal integers.  Returns interval of u**v.
 417       First, standardize to v[2]>0: */
 418    if v[2]>0 then bndlem1tonn(u,v)
 419    else if evnp(v[2]) then [recipl(neg8(nntonn(neg8(u[2]),1-v[2]))),
 420       recipu(nntonn(neg8(u[2]),-v[2]))]
 421    else [recipl(neg8(nntonn(neg8(u[2]),-v[2]))),
 422       recipu(nntonn(neg8(u[2]),1-v[2]))] $
 423
 424 bndlem1tonn(u,v) := /* u & v are intervals with u[1]>=1, v[2]>1.
 425       Returns interval for u**v. */
 426    if evnp(v[2]) then [neg8(nntonn(neg8(u[1]),v[2]-1)),
 427       nntonn(neg8(u[1]),v[2])]
 428    else [neg8(nntonn(neg8(u[1]),v[2])),nntonn(neg8(u[1]),v[2]-1)]$
 429
 430 bndmayspan1tonn(u,v) := /* u & v are intervals with 0<=u[1]<1 &
 431       u[2] symbolic & v[1]>=0.  Returns interval for u**v. */
 432    [nntonn(u[1],v[2]), max(u[2]**v[1], u[2]**v[2])] $
 433
 434 bndminus(u) := /* u is an interval.  returns interval for -u. */
 435    [neg8(u[2]), neg8(u[1])] $
 436
 437 bndnntimes(u,v) := /* u & v are intervals with u[1]>=0.  returns
 438       interval of u*v.  First, try to standardize lower bound
 439       of 2nd arg to nonnegative too: */
 440   if nonnegl(v[1]) then bndnntimnn(u,v)
 441   else if nonposu(v[2]) then bndminus(bndnntimnn(u,bndminus(v)))
 442   else if negl(v[1]) then
 443     if posu(v[2]) then [neg8(mgez(u[2],neg8(v[1]))),mgez(u[2],v[2])]
 444     else [neg8(mgez(u[2],neg8(v[1]))), max(u[1]*v[2], u[2]*v[2])]
 445   else if posu(v[2]) then [min(u[1]*v[1], u[2]*v[1]), mgez(u[2],v[2])]
 446   else [min(u[1]*v[1], u[2]*v[1]), max(u[1]*v[2], u[2]*v[2])] $
 447
 448 bndnegtimes(u,v) := /* u & v are intervals with u[1]<0.
 449       Returns interval of u*v. */
 450   if posu(u[2]) or posu(v[2]) and negl(v[1]) then
 451     [min(neg8(mgez(neg8(u[1]),v[2])), neg8(mgez(u[2],neg8(v[1])))),
 452     max(mgez(neg8(u[1]),neg8(v[1])), mgez(u[2],v[2]))]
 453   else if negl(v[1]) then [min(u[2]*v[2], u[2]*v[1], u[1]*v[2]),
 454     max(mgez(neg8(u[1]), neg8(v[1])), u[2]*v[2])]
 455   else if posu(v[2]) then [min(u[2]*v[1],neg8(mgez(neg8(u[1]),v[2]))),
 456     max(u[2]*v[2], u[2]*v[1], u[1]*v[1])]
 457   else [min(u[2]*v[2], u[2]*v[1], u[1]*v[2]),
 458     max(u[2]*v[2], u[2]*v[1], u[1]*v[1])] $
 459
 460 bndnntimnn(u,v) := /* u & v are intervals with u[1] & u[2]>=0.
 461       Returns interval for u*v. */
 462    [mgez(u[1],v[1]), mgez(u[2],v[2])] $
 463
 464 bndnptonnevn(u,v) := /* u & v are intervals with u[1]<=0 &
 465       v a nonnegative even integer.  Returns interval of u**v. */
 466    [nntonn(neg8(u[2]),v[1]), nntonn(neg8(u[1]), v[1])] $
 467
 468 bndrecip(u) := /* u is an interval not containing zzero in its
 469       interior.  Returns interval of 1/u. */
 470    [recipl(u[2]), recipu(u[1])] $
 471
 472 bndspan1tonn(u,v) := /* u & v are intervals with 0<=u[1]<1<u[2]
 473       & v[1]>=0.  Returns interval for u**v. */
 474    [nntonn(u[1],v[2]), nntonn(u[2],v[2])] $
 475
 476 bndunary(name,u) := /* Name is the name of a univariate
 477       nondecreasing function such as log, and u is the bounds of its
 478       argument.  Returns bounds1(name(argument)). */
 479    [unarybnd(name, u[1], 'plus), unarybnd(name, u[2], 'minus)] $
 480
 481 evnp(b) := /* b is integer. Returns true if it is even & false
 482       otherwise. */
 483    if integerp(b/2) then true else false $
 484
 485 gem1l(lb) := /* lb is a lowerbound.  Returns true if it is >=1,
 486       false otherwise. */
 487    if numberp(lb) and lb>=-1 or inpart(lb,0)='strict and
 488       numberp(inpart(lb,1)) and piece>=-1 then true
 489    else false $
 490
 491 ge1l(lb) := /* lb is a lowerbound.  Returns true if it is >=1,
 492       false otherwise. */
 493    if numberp(lb) and lb>=1 or inpart(lb,0)='strict and numberp(
 494       inpart(lb,1)) and piece>=1 then true
 495    else false $
 496
 497 ge1u(ub) := /* ub is an upperbound.  Returns true if it is >=1,
 498       false otherwise. */
 499    if ub='inf or numberp(ub) and ub>=1 or inpart(ub,0)='strict and
 500     (numberp(bounds1(inpart(ub,1))) and piece>1 or piece='inf)then true
 501    else false $
 502
 503 /*lbatom(w) := block(/* w is an indeterminate.  Returns its
 504       lowerbound, printing a message and establishing it as minf if
 505       none existed. */
 506    [ans],
 507    ans: get(w, lowerbound),
 508    if ans=false then (print("doing  put(", w, ", minf, lowerbound)"),
 509       put(w, 'minf, lowerbound),
 510       ans:'minf),
 511    return (ans)) $*/
 512 lbatom(w) := block(
 513    [ans],
 514    if w=%e then return(2.718281),
 515    if w=%pi then return(3.141592),
 516    ans: greaters(w),
 517    if ans=[] then (ans:geqs(w),
 518       if ans=[] then ans:'minf
 519       else ans: first(ans))
 520    else ans: strict(first(ans)),
 521    return(ans)) $
 522
 523 lem1l(lb) := /* lb is a lowerbound.  Returns true if it's <=-1,
 524       false otherwise. */
 525    if numberp(lb) and lb<=-1 or lb='minf or inpart(lb,0)='strict and
 526       (inpart(lb,1)='minf or numberp(piece) and piece<1) then true
 527    else false $
 528
 529 lem1u(ub) := /* ub is an upperbound.  Returns true if it's <=-1,
 530       false otherwise. */
 531    if numberp(ub) and ub<=-1 or inpart(ub,0)='strict and
 532       numberp(inpart(ub,1)) and piece<=-1 then true
 533    else false $
 534
 535 le1u(ub) := /* ub is an upperbound.  Returns true if it is <=1,
 536       false otherwise. */
 537    if numberp(ub) and ub<=1 or inpart(ub,0)='strict and
 538       numberp(inpart(ub,1)) and piece<=1 then true
 539    else false $
 540
 541 mgez(x,y) := /* x & y are bounds.  Returns x*y. */
 542    if x=0 or y=0 then 0
 543    else if x='inf or y='inf then 'inf
 544    else if inpart(x,0)='strict then
 545       if inpart(y,0)='strict then
 546          strict(mgez(inpart(x,1),inpart(y,1)))
 547       else strict(mgez(inpart(x,1),y))
 548    else if inpart(y,0)='strict then strict(mgez(x,inpart(y,1)))
 549    else x*y $
 550
 551 negl(lb) := /* lb is a lowerbound.  Returns true if it is <0,
 552       false otherwise. */
 553    if lb='minf or numberp(lb) and lb<0 or inpart(lb,0)='strict and
 554       (inpart(lb,1)='minf or numberp(piece) and piece<0) then true
 555    else false $
 556
 557 negu(ub) := /* ub is an upperbound.  Returns true if it is <0
 558       false otherwise. */
 559    if numberp(ub) and  ub<0 or inpart(ub,0)='strict and
 560       numberp(inpart(ub,1)) and piece<=0 then true
 561    else false $
 562
 563 neg8(b) := /* b is a bound.  Returns its negative. */
 564    if variablep(b) then
 565       if b='inf then 'minf
 566       else if b='minf then 'inf
 567       else -b
 568    else if inpart(b,0)='strict then strict(neg8(inpart(b,1)))
 569    else -b $
 570
 571 nntonn(x,y) := /* x & y are nonnegative bounds. Returns x**y. */
 572    if y=0 then 1
 573    else if x=0 then 0
 574    else if x='inf then 'inf
 575    else if x=1 then 1
 576    else if y='inf then
 577       if numberp(x) and x<1 or inpart(x,0)='strict and
 578          numberp(inpart(x,1)) and piece<1 then 0
 579       else 'inf
 580    else if inpart(x,0)='strict then
 581       if inpart(y,0)='strict then
 582          strict(nntonn(inpart(x,1),inpart(y,1)))
 583       else strict(nntonn(inpart(x,1),y))
 584    else if inpart(y,0)='strict then strict(nntonn(x,inpart(y,1)))
 585    else ev(x**y,numer) $
 586
 587 nonnegl(lb) := /* lb is a lower bound.  Returns true if it is
 588       nonnegative, false otherwise. */
 589    if lb=0 or posl(lb) then true else false $
 590
 591 nonposu(ub) := /* ub is an upperbound.  Returns true if it is
 592       positive, false otherwise. */
 593    if ub=0 or negu(ub) then true else false $
 594
 595 posl(lb) := /* lb is a lowerbound.  Returns true if it is >0,
 596       false otherwise. */
 597    if numberp(lb) and lb>0 or inpart(lb,0)='strict and
 598       numberp(inpart(lb,1)) and piece>=0 then true
 599    else false $
 600
 601 posu(ub) := /* ub is an upperbound.  Returns true if >0,
 602       false otherwise. */
 603    if ub='inf or numberp(ub) and ub>0 or inpart(ub,0)='strict and
 604       (inpart(ub,1)='inf or numberp(piece) and piece>=0) then true
 605    else false $
 606
 607
 608 recipl(ub) := /* ub is an upperbound.  Returns its 1/ub. */
 609    if ub = 'inf then 0
 610    else if ub=0 then 'minf
 611    else if inpart(ub,0)='strict then strict(recipl(inpart(ub,1)))
 612    else 1/ub $
 613
 614 recipu(lb) := /* lb is a lowerbound.  Returns its 1/lb. */
 615    if lb = 'minf then 0
 616    else if lb=0 then 'inf
 617    else if inpart(lb,0)='strict then
 618       strict(recipu(inpart(lb,1)))
 619    else 1/lb $
 620
 621 /*ubatom(w) := block(/* w is an indeterminate.
 622       Returns its upperbound, printing a message & establishing it as
 623       inf if none existed. */
 624    [ans],
 625    ans: get(w, upperbound),
 626    if ans=false then (print("doing  put(", w, ", inf, upperbound)"),
 627       put(w,'inf,upperbound),
 628       ans: 'inf),
 629    return(ans)) $*/
 630 ubatom(w) := block(
 631    [ans],
 632    if w=%e then return(2.718282),
 633    if w=%pi then return(3.141593),
 634    ans: lesses(w),
 635    if ans=[] then (ans:leqs(w),
 636       if ans=[] then ans:'inf
 637       else ans: first(ans))
 638    else ans: strict(first(ans)),
 639    return(ans)) $
 640
 641
 642 unarybnd(name, b, d) := block(/* Name is name of a univariate
 643       nondecreasing function like log, b is a bound of its argument,
 644       and d is plus for a lower bound or minus for an upperbound.
 645       Returns the corresponding bound of name(argument). */
 646    [arg],
 647    if inpart(b,0) = 'strict then
 648       arg: strict(limit(apply(name,[arg]), arg, inpart(b,1), d))
 649    else arg: limit(apply(name,[arg]), arg, b, d),
 650    return(ev(arg,numer))) $
 651
 652 definitecode(a) := block( /*lagrange's */
 653    [n, perm, b, ii, jj, kk, npos, nneg, nzero, nnpos, nnneg, nunkn,
 654       partswitch, prederror],
 655    prederror:false,  partswitch: true,  n: length(a),  perm: [],
 656    npos: nneg: nzero: nnpos: nnneg: nunkn: 0,
 657    for i:n step -1 thru 1 do perm: cons(i, perm),
 658    for i:1 thru n while (npos=0 or nneg=0) do(
 659       jj: i,
 660       while jj<=n and a[ii:perm[jj],ii]=0 do jj: jj+1,
 661       if jj>n then (nzero: n+1-i,
 662          for j:i thru n while npos=0 or nneg=0 do(ii: perm[j],
 663             for k:i thru n do if a[ii,perm[k]]#0 then npos:nneg:1))
 664       else (perm[jj]:perm[i],  perm[i]:ii,
 665          b: bounds1(a[ii,ii]),
 666          if posl(b[1]) then npos: npos+1
 667          else if negu(b[2]) then nneg: nneg+1
 668          else if b[1]=0 then
 669             if b[2]=0 then nzero:nzero+1
 670             else nnneg: nnneg+1
 671          else if b[2]=0 then nnpos: nnpos+1
 672          else nunkn: nunkn+1,
 673          for j:i+1 thru n do (jj: perm[j],
 674             b: -a[jj,ii]/a[ii,ii],
 675             for k:i+1 thru n do (kk: perm[k],
 676                a[jj,kk]: a[jj,kk] + b*a[ii,kk])))),
 677    if npos>0 then
 678       if nneg>0 then return(/*indefinite*/ 8)
 679       else if nnpos>0 then return(/*pos semi or indef*/ 5)
 680       else if nunkn=0 then
 681          if nzero=0 and nnneg=0 then return(/*pos def*/ 7)
 682          else return(/*pos semi*/ 6)
 683       else return(/*pos def, pos semi, or indef*/ 5)
 684    else if nneg>0 then
 685       if nnneg>0 then return(/*neg semi or indef*/ 3)
 686       else if nunkn=0 then
 687          if nzero=0 and nnpos=0  then return(/*neg def*/ 1)
 688          else return(/*neg semi*/ 2)
 689       else return(/*neg def, neg semi, or indef*/ 3)
 690    else if nunkn=0 then
 691       if nnpos=0 then
 692          if nnneg=0 then return(/*rank 0*/ 4)
 693          else return(/*pos semi*/ 6)
 694       else if nnneg=0 and nzero=0 then return(/*neg def or semi*/ 2)
 695       else return(/*unknown*/ 9)
 696    else return(/*unknown*/ 9)) $
 697
 698 ttyoff:false $
 699
 700