share/calculus/asympa.mac

   1 ttyoff: true $
   2
   3 /* Set up pattern-match rules to standardize bases of logs and
   4 exponentials: */
   5 matchdeclare(xtrue,true, non%econst,non%econstp) $
   6 tellsimpafter(log2(xtrue), 1.44269504*log(xtrue)) $
   7 tellsimp(non%econst^xtrue, %e^(log(non%econst)*xtrue)) $
   8
   9 non%econstp(ex) := is(baseconvert=true and constantp(ex) and ex#%e) $
  10
  11 /* Initialize global variables: */
  12 baseconvert: false$
  13
  14 taylormax: 0 $
  15
  16
  17 asymp(fn) := block(/* Returns 'ASYMP(expr), where expr is a simplified
  18       asymptotic equivalent to FN.  Uses the global integer variable
  19       named MAXTAYLOR to control optional use of the TAYLOR-series
  20       technique. */
  21    [asymvar, asympnt, ans, ratexpand, ratdenomdivide],
  22    asympnt: asympdata(fn),   asymvar: asympnt[1],
  23    if asymvar=[] then return('asymp(fn)),
  24    asympnt: ans: asympnt[2],
  25    ans: for jj:1 thru taylormax while ans#[]  do (
  26       ans: errcatch(apply('taylor, [fn,asymvar,asympnt,jj])),
  27       if ans#[] and(ans:ratdisrep(first(ans)))#0 then return(ans)),
  28    if ans#'done then return(ans),
  29    ans: for jj:1 thru taylormax while ans#[] do (
  30       ans: errcatch(apply('taylor, [1/fn,asymvar,asympnt,jj])),
  31       if ans#[] and(ans:ratdisrep(first(ans)))#0 then return(1/ans)),
  32    if ans#'done then return(ans),
  33    ratexpand: true,   ratdenomdivide: false,
  34    apply('ratvars, asymvar),
  35    ans: ratsimp(leadterm(ratnumer(fn))/leadterm(ratdenom(fn))),
  36    ratvars() /* Warning:  might alter user environment */,
  37    'asymp(ans)) $
  38
  39 asympdata(ex) := block(/* Returns a list of 2 lists, the first
  40       being the asymptotic variables in expression EX and the
  41       second being their limits. */
  42    [allvar, asymvar, asympnt, ans],
  43    allvar: listofvars(ex),
  44    asymvar: asympnt: [],
  45    for var in allvar do
  46        if (ans:get(var,'limit))#false then
  47           (asymvar: cons(var, asymvar),
  48            asympnt: cons(ans, asympnt)),
  49    [asymvar, asympnt]) $
  50
  51 leadterm(ex) := block(/* Same as ASYMP, except EX is a numerator
  52       or denominator */
  53    [saved],
  54    saved: pass1(ex),
  55    if saved=[] then return(ex),
  56    saved: pass2(saved),
  57    if saved=[] then ex
  58    else saved) $
  59
  60 pass1(ex) := block(/* Returns a list of lists, with each of the
  61       latter being an asymptotically significant cofactor followed
  62       by its coefficient.  If a coefficient is 0 or contains
  63       'IND, then cancellation does or could take place.  [] is
  64       returned when EX is not a sum or when dominance relations
  65       between terms are uncertain. */
  66    [saved],
  67    if mapatom(ex) or inpart(ex,0)#"+" then return([]),
  68    ex: substinpart("[",ex,0),
  69    saved: [[first(ex), 1]],
  70    while (ex:rest(ex))#[] do saved: update(saved,first(ex)),
  71    saved) $
  72
  73 update(saved,term) := block(/* SAVED is as indicated for PASS1,
  74       and TERM is a nonsum.  Returns an appropriately updated version
  75       of SAVED, according to the dominance between TERM and each
  76       cofactor in SAVED. Also, sets global EX to [[]] when it is time to
  77       preemptively terminate the loop in PASS1.  */
  78    [sav, nusaved, lim],
  79    nusaved: [],
  80    while saved#[] do(
  81       sav: first(saved),
  82       lim: multilim(term/sav[1]),
  83       if lim='savdominates then(
  84          nusaved: append(saved,nusaved),
  85          saved: [])
  86       else if lim='termdominates then saved: rest(saved)
  87       else if lim='bothdominate then(
  88          nusaved: cons(sav, nusaved),
  89          saved: rest(saved))
  90       else if lim='invalid then ex: [saved: []]
  91       else (/* asymptotically proportional */
  92          nusaved: cons(mrge(sav, term, lim),append(nusaved,rest(saved))),
  93          saved: [])),
  94    if lim='termdominates or lim='bothdominate then
  95       nusaved: cons([term,1], nusaved),
  96    nusaved) $
  97
  98 multilim(ratio) := block(/* RATIO is term/savedterm.  Returns
  99       'SAVDOMINATES, 'TERMDOMINATES, 'BOTHDOMINATE, 'INVALID, or
 100       the finite asymptotic limit of RATIO.  Uses global variables
 101       ASYMPNT and ASYMVAR established by ASYMP. */
 102    [mlim, ulim, pnts],
 103    ratio: [ratio,1],
 104    pnts: asympnt,
 105    mlim: 1,
 106    for var in asymvar while mlim#'invalid do(
 107       ratio: mypartition(ratio[1],var),
 108       if freeof(delete(var,asymvar),ratio[2]) then (
 109          ulim: errcatch(limit(ratio[2], var, first(pnts))),
 110          if ulim=[] or (ulim:first(ulim))='und or
 111                not freeof(nounify('limit),ulim) then mlim:'invalid
 112          else if mlim#'bothdominate then (
 113             if ulim=0 then
 114                if mlim='termdominates then mlim:'bothdominate
 115                else mlim: savdominates
 116             else if member(ulim,'[inf,minf,infinity]) then
 117                if mlim='savdominates then mlim:'bothdominate
 118                else mlim: 'termdominates
 119             else if mlim#'termdominates and mlim#savdominates then
 120                mlim: ulim*mlim)),
 121       pnts: rest(pnts)),
 122    if not member(mlim,'[savdominates,termdominates,bothdominate,invalid])
 123       then return(ratio[1]*mlim),
 124    mlim) $
 125
 126 mypartition(nonsum,var) := /* Returns a list consisting of
 127       factors of NONSUM freeof VAR, then other factors. */
 128    if mapatom(nonsum) or inpart(nonsum,0)#"*" then
 129       if freeof(var, nonsum) then [nonsum,1]
 130       else [1,nonsum]
 131    else partition(nonsum,var) $
 132
 133 mrge(sav,term,lim) := block( /* Returns a list consisting of the
 134       cofactor and coefficient of the simplest asymptotic form for TERM
 135       + SAV[1]*SAV[2], where TERM/SAV[1] approaches LIM. */
 136    [told, tnew],
 137    told: ratsimp(sav[2]+lim), tnew: ratsimp((lim+sav[2])/lim),
 138    if complexity(sav[1]*told) < complexity(term*tnew)
 139       then [sav[1], told]
 140    else [term, tnew]) $
 141
 142 complexity(expr) := /* Returns a complexity measure of
 143       expression EXPR. */
 144    if mapatom(expr) then 1
 145       else 1 + argscomplexity(substinpart("[", expr, 0)) $
 146
 147 argscomplexity(args) := /* Returns a complexity measure of
 148       a list of expressions, ARGS. */
 149    if args=[] then 0
 150    else complexity(first(args)) + argscomplexity(rest(args)) $
 151
 152 pass2(saved) := block(/* Given the list of lists returned by
 153       PASS1, this function returns [] if any cancellations are
 154       possible.  Otherwise, this function returns the sum of the
 155       terms represented by SAVED */
 156    [res, sav],
 157    res: 0,
 158    while saved#[] do (
 159       sav: first(saved),
 160       if freeof('ind,sav[2]) and asksign(abs(sav[2]))#'zero then(
 161          res: res + sav[1]*sav[2],
 162          saved: rest(saved))
 163       else res: saved: []),
 164    res) $
 165
 166 theta(ex) :=  (/* Returns 'THETA(expr), where expr is a
 167       simplified expression having the same exact order as EX. */
 168    ex: first(asymp(ex)),
 169    'theta(ratsimp(theta1(num(ex))/theta1(denom(ex))))) $
 170
 171 theta1(ex) := block(/* Same as THETA, except EX is a numerator
 172       or denominator. */
 173    [asymvar, terms, gc],
 174    asymvar: asympdata(ex)[1],
 175    if mapatom(ex) or inpart(ex,0)#"+" then terms:[ex]
 176    else terms: substinpart("[",ex,0),
 177    gc: cfasymp(first(terms)),
 178    if is(gc=0) then return(gc),
 179    terms: rest(terms),
 180    for term in terms while is(gc#1) do gc: gcd(gc,cfasymp(term)),
 181    ex/gc) $
 182
 183 cfasymp(term) := block(/* Returns the coefficient of TERM with
 184       respect to factors containing asymptotic variables. */
 185    [cf],
 186    if mapatom(term) or inpart(term,0)#"*" then term: [term]
 187    else term: substinpart("[",term,0),
 188    cf: 1,
 189    for fctr in term do if apply('freeof, endcons(fctr,asymvar))
 190       then cf: cf*fctr,
 191    cf) $
 192
 193 omega(ex) := /* Returns 'OMEGA(expr), where expr is a simplified
 194       expression having at most the same order as EX. */
 195    'omega(ooromega(ex,false)) $
 196
 197 lomega(ex) := /* Returns 'LOMEGA(expr), where expr is a simplified
 198       expression which is of lesser order than every expression which
 199       EX is lesser order than. */
 200    'lomega(ooromega(ex,false)) $
 201
 202 o(ex) := /* Returns 'O(expr), where expr is a simplified expression
 203       having at least the same order as EX. */
 204    'o(ooromega(ex,true)) $
 205
 206 lo(ex) := /* Returns 'LO(expr), where expr is a simplified
 207       expression which is of greater order than every expression
 208       which EX is greater order than. */
 209     'lo(ooromega(ex,true)) $
 210
 211 ooromega(ex,flag) := /* If FLAG is TRUE, returns product of
 212       highest-order asymptotic factors from terms of EX.  Otherwise
 213       returns product of lowest-order factors. */
 214    (ex: first(theta(ex)),
 215    ratdisrep(radcan(oom(ratnumer(ex),flag)/oom(ratdenom(ex),not flag)))) $
 216
 217 oom(ex,bigo) := block(/* Same as OOROMEGA, except EX is a
 218       numerator or denominator. */
 219    [asymvar, asympnt, ratio, limratio, ans, pnts],
 220    asympnt: asympdata(ex),
 221    asymvar: asympnt[1], asympnt: asympnt[2],
 222    if mapatom(ex) or inpart(ex,0)#"+" then return(ex),
 223    ex: substinpart("[", ex, 0),
 224    ans: first(ex),   ex: rest(ex),
 225    for term in ex while ans#'invalid do (
 226       ratio: [term/ans, 1],
 227       pnts: asympnt,
 228       for var in asymvar while ans#'invalid do (
 229          ratio: mypartition(ratio[1],var),
 230          if freeof(delete(var,asymvar),ratio[2]) then (
 231             limratio: errcatch(limit(ratio[2], var, first(pnts))),
 232             if limratio=[] or (limratio:first(limratio))='und
 233                or not freeof(nounify('limit),limratio) then ans:'invalid
 234             else if limratio=0 and not bigo
 235                 or member(limratio,'[inf,minf,infinity]) and bigo
 236                 then ans: ans*ratio[2]),
 237          pnts: rest(pnts))),
 238    ans) $
 239
 240 asympseries(uu,nt) := block(/* Returns an NT term asymptotic
 241       expansion for UU. */
 242    [answer, term],
 243    answer: 0,
 244    for j: 1 thru nt do (
 245       term: inpart(asymp(uu), 1),
 246       answer: answer + term,
 247       uu: uu - term),
 248    answer) $
 249
 250 /* Establish table for combining different kinds of approximate
 251    expressions: */
 252 asym[1]:'asymp$  asym[2]:'theta$  asym[3]:'o$
 253 asym[4]:'lo$  asym[5]:'omega$  asym[6]:'lomega$
 254 put('asymp,1,'level)$  put('theta,2,'level)$  put('o,3,'level)$
 255 put('lo,4,'level)$  put('omega,5,'level)$  put('lomega,6,'level)$
 256
 257 asympsimp(ex) := /* Returns a simplified version of EX, if EX
 258       contains approximate subexpressions. */
 259    scanmap('asimp, ex, 'bottomup) $
 260
 261 asimp(ex) := block(/* EX has its subexpressions already simplified recursively.
 262       This function properly combines any top-level approximate
 263       subexpressions therein. */
 264    [inflag, maxlev, lev, asyms, nonasyms],
 265    if mapatom(ex) then return(ex),
 266    if inpart(ex,0)="^" then
 267       if not mapatom(inpart(ex,1))
 268          and get(verbify(inpart(piece,0)),'level)#false
 269          and is(sign(inpart(ex,2))='pos)
 270          then return(apply(inpart(ex,1,0),
 271             [inpart(ex,1,1)^inpart(ex,2)]))
 272       else return(ex),
 273    inflag:true,  maxlev:0,  asyms:nonasyms:[],
 274    if piece="*" then (
 275       for fctr in ex do
 276          if mapatom(fctr) or(lev:get(verbify(inpart(fctr,0)),'level))=false
 277             then asyms: cons(fctr,asyms)
 278          else if incomensurate(lev,maxlev) then
 279             nonasyms: cons(inpart(fctr,1), nonasyms)
 280          else (asyms: cons(inpart(fctr,1), asyms),
 281             maxlev: max(lev, maxlev)),
 282       if maxlev=0 then return(ex),
 283       asyms: apply(asym[maxlev], if rest(asyms)=[] then asyms
 284             else [apply("*", asyms)]),
 285       if nonasyms#[] then asyms: asyms*apply(
 286          if maxlev<=4 then 'omega else 'o, if rest(nonasyms)=[] then
 287          nonasyms else apply("*",nonasyms)),
 288       return(asyms)),
 289    if piece="+" then (
 290       for trm in ex do
 291          if mapatom(trm) or (lev:get(verbify(inpart(trm,0)),'level))=false or
 292             incomensurate(lev,maxlev) then
 293             nonasyms: cons(trm,nonasyms)
 294          else (asyms: cons(inpart(trm,1),asyms),
 295             maxlev: if maxlev<=3 and lev=4 or maxlev=4 and lev<=3 then 3
 296                else if lev=6 and (maxlev=5 or maxlevel<3) or
 297                   maxlev=6 and (lev=5 or lev<3) then 5
 298                else max(maxlev,lev)),
 299       if maxlev=0 then return(ex),
 300       asyms: apply(asym[maxlev], if rest(asyms)=[] then asyms
 301             else [apply("+",asyms)]),
 302       if nonasyms#[] then asyms: asyms+(if rest(nonasyms)=[]
 303             then first(nonasyms) else apply("+",nonasyms)),
 304       return(asyms)),
 305    asyms: piece,
 306    if (lev:get(verbify(piece),'level))#false and not mapatom(inpart(ex,1))
 307       and (maxlev:get(verbify(inpart(piece,0)),'level))#false then
 308       if incomensurate(lev,maxlev) then return(ex)
 309       else return(apply(asym[max(lev,maxlev)], [inpart(ex,1,1)])),
 310    ex) $
 311
 312 incomensurate(l1,l2) := /* Returns TRUE if approximation operators having
 313       numeric codes L1 and L2 cannot coalesce, returning FALSE
 314       otherwise. */
 315     (l1=3 or l1=4) and (l2=5 or l2=6) $
 316
 317 ttyoff: false $