share/simplification/functs.mac

   1 /* eval_when([translate,load,loadfile,batch,demo],
   2  matchdeclare(a,nonzeroandfreeof(x),[b,c],freeof(x)))$ */
   3
   4 define_variable(takegcd,true,boolean,"used in gcdivide to decide the
   5 gcd choice");
   6
   7 rempart&&
   8   rempart (expr, n) :=
   9     if symbolp (expr) then
  10       'rempart (expr, n)
  11     else if atom (expr) then
  12       error ("rempart expects a compound object as its first argument")
  13     else if symbolp (n) then
  14       'rempart (expr, n)
  15     else
  16       block ([elts: args (expr)],
  17         local (elts),
  18         if numberp (n) then
  19           if not integerp (n) then
  20             error ("Non-integer n in rempart")
  21           else if is (n < 1 or n > length (elts)) then
  22             expr
  23           else
  24             apply (op (expr),
  25                    append (makelist (elts[i], i, 1, n-1),
  26                            makelist (elts[i], i, n+1, length (elts))))
  27         else if not (listp (n)) then
  28           'rempart (expr, n)
  29         else if not (length (n) = 2) then
  30           error ("If second argument to rempart is a list, it should be a pair")
  31         else
  32           block ([a: first (n), b: second (n)],
  33             if not (numberp (a) and numberp (b)) then
  34               'rempart (expr, n)
  35             else if not (integerp (a) and integerp (b)) then
  36               error ("Non-integer [a,b] to remove in rempart")
  37             else if is (b < a) then
  38               expr
  39             else (
  40               a: min (max (a, 0), length (elts) + 1),
  41               b: min (max (b, 0), length (elts) + 1),
  42               apply (op (expr),
  43                      append (makelist (elts[i], i, 1, a-1),
  44                              makelist (elts[i], i, b+1, length (elts)))))))$
  45
  46 wronskian&&  wronskian(functlist,var):=block([end],
  47     end:length(functlist)-1,
  48     functlist:[functlist],
  49     thru end do functlist:endcons(map(lambda([x],diff(x,var)),
  50         last(functlist)),functlist),
  51     apply('matrix,functlist))$
  52
  53 tracematrix&&  tracematrix(m):=block([sum,len],sum:0,len:length(m),
  54 for i:1 thru len do sum:sum+part(m,i,i),sum)$
  55
  56 rational&&  rational(z):=block([n,d,cd,ratfac],
  57     ratfac:false,
  58     n:ratdisrep(ratnumer(z)*(cd:conjugate(d:ratdenom(z)))),
  59     d:rat(n/ratdisrep(d*cd)),
  60     if ratp(z) then d else ratdisrep(d))$
  61
  62 nonzeroandfreeof&&  nonzeroandfreeof(x,e):=is(e#0 and freeof(x,e))$
  63
  64 lcm&& lcm([list]):=block([listconstvars:false],if listofvars(list)=[] then
  65 lcm1(list) else factor(lcm1(list)))$
  66
  67 /* Replaced by the following routine
  68 lcm1(list):=if list=[] then 1 else block([rlist:rest(list),flist:first(list),
  69 frlist,partswitch:true,inflag:true,piece], if rlist=[] then flist else
  70 lcm1(cons(flist*(frlist:first(rlist))/gcd(flist,frlist),rest(rlist))))$
  71 */
  72
  73 /* New implementation of the function lcm
  74  * Do not use a recursive, but an iterative algorithm.
  75  * Check more carefully the case division by zero.
  76  */
  77
  78 lcm1(list):=
  79    block
  80    (
  81     [flist, rlist, result, keepfloat:true, ratprint:false],
  82
  83     result : flatten(list),
  84     if result = [] then result : [1],
  85     unless rest(result)=[] do
  86     (
  87      flist : first(result),
  88      rlist : first(rest(result)),
  89      if rlist=0 then
  90         result : [0]
  91      else
  92         result : cons( (flist*rlist / gcd(flist, rlist)),
  93                        (rest(rest(result))))
  94     ),
  95     first(result)
  96 )$
  97
  98 gcdivide&&  gcdivide(poly1,poly2):=block([gcdlist],
  99                 gcdlist:if takegcd then ezgcd(poly1,poly2)
 100                         else [1,poly1,poly2],
 101                 gcdlist[2]/gcdlist[3])$
 102
 103 series&&  arithmetic(a,d,n):=a+(n-1)*d$
 104           geometric(a,r,n):=a*r^(n-1)$
 105           harmonic(a,b,c,n):=a/(b+(n-1)*c)$
 106           arithsum(a,d,n):=n*(a+(n-1)*d/2)$
 107           geosum (a,r,n) :=
 108             block([dummyvar: gensym()],
 109               if n = 'inf then
 110                 block ([coef_sign: sign (a)],
 111                   if member (sign (r-1), ['pos, 'zero, 'pz]) then
 112                     if member (coef_sign, ['pos, 'neg, 'zero]) then
 113                       limit (a * inf)
 114                     else
 115                       'geosum (a, r, n)
 116                   else if coef_sign = 'zero then 0
 117                   else
 118                     block ([sgn_abs: sign (abs (r) - 1)],
 119                       if member (sgn_abs, ['pos, 'pz, 'zero]) then
 120                         if not (member (coef_sign, ['pz, 'nz, 'pnz])) then
 121                           'und
 122                         else
 123                           'geosum (a, r, n)
 124                       else if sgn_abs # 'neg then
 125                         'geosum (a, r, n)
 126                       else
 127                         a/(1 - r)))
 128               else
 129                 a * (1 - r^n) / (1 - r))$
 130
 131 gauss&&  gaussprob(x):=1/sqrt(2*%pi)*%e^(-x^2/2)$
 132
 133 /* See, for example, http://en.wikipedia.org/wiki/Gudermannian_function */
 134 gd&&  gd(x):=2*atan(%e^x)-%pi/2$
 135         agd(x):=log(tan(%pi/4+x/2))$
 136
 137 trig&&  vers(x):=1-cos(x)$
 138         covers(x):=1-sin(x)$
 139         exsec(x):=sec(x)-1$
 140         hav(x):=(1-cos(x))/2$
 141
 142 combination&&  combination(n,r):=binomial(n,r)$
 143         permutation(n,r):=binomial(n,r)*r!$