1ttyoff:true $
2load("qualsp");
3matchdeclare([utrue,vtrue,wtrue],true)$
4tellsimp(qual(utrue), qual1(utrue, listofvars(utrue))) $
5tellsimp(qual(utrue,vtrue),qual1(utrue,qual_listify(vtrue))) $
6tellsimp(revelation(utrue), revelation1(utrue,200,300)) $
7tellsimp(revelation(utrue,vtrue), revelation1(utrue,vtrue,300))$
8tellsimp(revelation(utrue,vtrue,wtrue),
9   revelation1(utrue,vtrue,wtrue)) $
10tellsimp(slopes(utrue),slopes1(utrue,listofvars(utrue)))$
11tellsimp(slopes(utrue,vtrue),slopes1(utrue,qual_listify(vtrue)))$
12tellsimp(symmetry(utrue),symmetry1(utrue,listofvars(utrue)))$
13tellsimp(symmetry(utrue,vtrue),symmetry1(utrue,qual_listify(vtrue)))$
14tellsimp(periods(utrue), periods1(utrue,listofvars(utrue))) $
15tellsimp(periods(utrue,vtrue),periods1(utrue,qual_listify(vtrue))) $
16tellsimp(limits(utrue),limits1(utrue,listofvars(utrue)))$
17tellsimp(limits(utrue,vtrue),limits1(utrue,qual_listify(vtrue)))$
18tellsimp(stationarypoints(utrue),stationarypoints1(utrue,
19   listofvars(utrue)))$
20tellsimp(stationarypoints(utrue,vtrue),stationarypoints1(utrue,
21   qual_listify(vtrue))) $
22
23variablep(u) := is(atom(u) and not numberp(u) or subvarp(u)) $
24
25qual_listify(u) :=
26   if listp(u) then u else [u] $
27
28qual1(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
35revelation1(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
43slopes1(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
49slopes2(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
60curvature(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
67qual_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
73gradient(u,v) := block(
74   [ans],
75   ans: [],
76   for x in v do ans: endcons(diff(u,x), ans),
77   return(ans)) $
78
79
80symmetry1(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
88symmetry2(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
120periods1(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
129period2(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
152lcmspec(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
159limits1(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
175zerosandsingularities(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
183zp1(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])) $
193conssingularities(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
212stationarypoints1(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
232type(u) :=
233   ['maximum, 'nonminimum, 'nonminimum, 'unknown, 'nonmaximum,
234      'nonmaximum, 'minimum, 'saddle, 'unknown][u] $
235
236bounds(w) := ev(bounds1(w),prederror:false,partswitch:true)$
237
238bounds1(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
377bndrestrict(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
383addbnd(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
394bndge1to(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
411bndge1tonn(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
415bndlem1to(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
424bndlem1tonn(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
430bndmayspan1tonn(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
434bndminus(u) := /* u is an interval.  returns interval for -u. */
435   [neg8(u[2]), neg8(u[1])] $
436
437bndnntimes(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
448bndnegtimes(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
460bndnntimnn(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
464bndnptonnevn(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
468bndrecip(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
472bndspan1tonn(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
476bndunary(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
481evnp(b) := /* b is integer. Returns true if it is even & false
482      otherwise. */
483   if integerp(b/2) then true else false $
484
485gem1l(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
491ge1l(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
497ge1u(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)) $*/
512lbatom(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
523lem1l(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
529lem1u(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
535le1u(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
541mgez(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
551negl(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
557negu(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
563neg8(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
571nntonn(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
587nonnegl(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
591nonposu(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
595posl(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
601posu(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
608recipl(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
614recipu(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)) $*/
630ubatom(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
642unarybnd(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
652definitecode(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
698ttyoff:false $
699
700
701