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