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