diff --git a/books/bookvol10.2.pamphlet b/books/bookvol10.2.pamphlet index 88c7428..24ea0fe 100644 --- a/books/bookvol10.2.pamphlet +++ b/books/bookvol10.2.pamphlet @@ -872,6 +872,8 @@ CombinatorialFunctionCategory(): Category == with ++ binomial(n,r) returns the \spad{(n,r)} binomial coefficient ++ (often denoted in the literature by \spad{C(n,r)}). ++ Note: \spad{C(n,r) = n!/(r!(n-r)!)} where \spad{n >= r >= 0}. + ++ + ++X [binomial(5,i) for i in 0..5] factorial : $ -> $ ++ factorial(n) computes the factorial of n ++ (denoted in the literature by \spad{n!}) diff --git a/books/bookvol10.4.pamphlet b/books/bookvol10.4.pamphlet index b603e5b..4537e98 100644 --- a/books/bookvol10.4.pamphlet +++ b/books/bookvol10.4.pamphlet @@ -6326,6 +6326,145 @@ CoerceVectorMatrixPackage(R : CommutativeRing): public == private where @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{package COMBF CombinatorialFunction} +<>= +)set break resume +)sys rm -f CombinatorialFunction.output +)spool CombinatorialFunction.output +)set message test on +)set message auto off +)clear all +--S 1 of 6 +f := operator 'f +--R +--R +--R (1) f +--R Type: BasicOperator +--E 1 + +--S 2 of 6 +D(product(f(i,x),i=1..m),x) +--R +--R +--R m m f (i,x) +--R ++-++ --+ ,2 +--R (2) | | f(i,x)> -------- +--R | | --+ f(i,x) +--R i= 1 i= 1 +--R Type: Expression Integer +--E 2 + +--S 3 of 6 +)set expose add constructor OutputForm +--R +--I OutputForm is already explicitly exposed in frame frame0 +--E 3 + +--S 4 of 6 +pascalRow(n) == [right(binomial(n,i),4) for i in 0..n] +--R +--R Type: Void +--E 4 + +--S 5 of 6 +displayRow(n)==output center blankSeparate pascalRow(n) +--R +--R Type: Void +--E 5 + +--S 6 of 6 +for i in 0..7 repeat displayRow i +--R +--R Compiling function pascalRow with type NonNegativeInteger -> List +--R OutputForm +--R Compiling function displayRow with type NonNegativeInteger -> Void +--R 1 +--R 1 1 +--R 1 2 1 +--R 1 3 3 1 +--R 1 4 6 4 1 +--R 1 5 10 10 5 1 +--R 1 6 15 20 15 6 1 +--R 1 7 21 35 35 21 7 1 +--R Type: Void +--E 6 + + +)spool +)lisp (bye) +@ +<>= +==================================================================== +CombinatorialFunction examples +==================================================================== + +f := operator 'f + + + (1) f + Type: BasicOperator + +D(product(f(i,x),i=1..m),x) + + + m m f (i,x) + ++-++ --+ ,2 + (2) | | f(i,x)> -------- + | | --+ f(i,x) + i= 1 i= 1 + Type: Expression Integer + + +The binomial(n, r) returns the number of subsets of r objects +taken among n objects, i.e. n!/(r! * (n-r)!) + +The binomial coefficients are the coefficients of the series expansion +of a power of a binomial, that is + + n + --+ / n \ k n + > | | x = (1 + x) + --+ \ k / + k= 0 + +This leads to the famous pascal triangle. First we expose the OutputForm +domain, which is normally hidden, so we can use it to format the lines. + + )set expose add constructor OutputForm + +Next we define a function that will output the list of binomial coefficients +right justified with proper spacing: + + pascalRow(n) == [right(binomial(n,i),4) for i in 0..n] + +and now we format the whole line so that it looks centered: + + displayRow(n)==output center blankSeparate pascalRow(n) + +and we compute the triangle + + for i in 0..7 repeat displayRow i + +giving the pretty result: + + Compiling function pascalRow with type NonNegativeInteger -> List + OutputForm + Compiling function displayRow with type NonNegativeInteger -> Void + 1 + 1 1 + 1 2 1 + 1 3 3 1 + 1 4 6 4 1 + 1 5 10 10 5 1 + 1 6 15 20 15 6 1 + 1 7 21 35 35 21 7 1 + +See Also: +o )show CombinatorialFunction +o )d op binomial +o )show OutputForm +o )help set + +@ \pagehead{CombinatorialFunction}{COMBF} \pagepic{ps/v104combinatorialfunction.ps}{COMBF}{1.00} @@ -6351,52 +6490,11 @@ CoerceVectorMatrixPackage(R : CommutativeRing): public == private where \cross{COMBF}{?**?} && \end{tabular} -<>= -)abbrev package COMBF CombinatorialFunction -++ Provides the usual combinatorial functions -++ Author: Manuel Bronstein, Martin Rubey -++ Date Created: 2 Aug 1988 -++ Date Last Updated: 30 October 2005 -++ Description: -++ Provides combinatorial functions over an integral domain. -++ Keywords: combinatorial, function, factorial. -++ Examples: )r COMBF INPUT - -CombinatorialFunction(R, F): Exports == Implementation where - R: Join(OrderedSet, IntegralDomain) - F: FunctionSpace R - - OP ==> BasicOperator - K ==> Kernel F - SE ==> Symbol - O ==> OutputForm - SMP ==> SparseMultivariatePolynomial(R, K) - Z ==> Integer - - POWER ==> "%power"::Symbol - OPEXP ==> "exp"::Symbol - SPECIALDIFF ==> "%specialDiff" - SPECIALDISP ==> "%specialDisp" - SPECIALEQUAL ==> "%specialEqual" - - Exports ==> with - belong? : OP -> Boolean - ++ belong?(op) is true if op is a combinatorial operator; - operator : OP -> OP - ++ operator(op) returns a copy of op with the domain-dependent - ++ properties appropriate for F; - ++ error if op is not a combinatorial operator; - "**" : (F, F) -> F - ++ a ** b is the formal exponential a**b; - binomial : (F, F) -> F - ++ binomial(n, r) returns the number of subsets of r objects - ++ taken among n objects, i.e. n!/(r! * (n-r)!); -@ - +\subsubsection{binomial} We currently simplify binomial coefficients only for non-negative integral second argument, using the formula $$ \binom{n}{k}=\frac{1}{k!}\prod_{i=0..k-1} (n-i),$$ -except if the second argument is symbolic: in this case [[binomial(n,n)]] is +except if the second argument is symbolic: in this case binomial(n,n) is simplified to one. Note that there are at least two different ways to define binomial coefficients @@ -6411,7 +6509,7 @@ second argument. This is, partially, also the approach taken in m = 0 => 1 \end{verbatim} -Of course, here [[n]] and [[m]] are integers. This definition agrees with the +Of course, here $n$ and $m$ are integers. This definition agrees with the recurrence $$\binom{n}{k}+\binom{n}{k+1}=\binom{n+1}{k+1}.$$ @@ -6433,11 +6531,116 @@ $\Gamma(n_0+1)$ is unbounded. However, since for $k\in {\bf Z}$, $n\in {\bf Z}$ and $0 < k < n$ both definitions agree, one could also combine them. This is what, for example, -Mathematica does. It seems that MuPAD sets [[binomial(n,n)=1]] for all -arguments [[n]], and returns [[binomial(-2, n)]] unevaluated. Provisos may help +Mathematica does. It seems that MuPAD sets binomial(n,n)=1 for all +arguments $n$, and returns binomial(-2, n) unevaluated. Provisos may help here. +\subsubsection{dvsum and dvdsum} +The dvsum and dvdsum operations implement differentiation of sums with +and without bounds. Note that the function +$$n\mapsto\sum_{k=1}^n f(k,n)$$ +is well defined only for integral values of $n$ greater than or equal to zero. +There is not even consensus how to define this function for $n<0$. Thus, it is +not differentiable. Therefore, we need to check whether we erroneously are +differentiating with respect to the upper bound or the lower bound, where the +same reasoning holds. + +Differentiating a sum with respect to its indexing variable correctly gives +zero. This is due to the introduction of dummy variables in the internal +representation of a sum: the operator [[%defsum]] takes 5 arguments, namely + +\begin{enumerate} +\item the summands, where each occurrence of the indexing variable is replaced + by +\item the dummy variable, +\item the indexing variable, +\item the lower bound, and +\item the upper bound. +\end{enumerate} + +\subsubsection{dvprod and dvdprod} +The dvprod and dvdprod operations implement differentiation of products with +and without bounds. Note again, that we cannot even properly define +products with bounds that are not integral. + +To differentiate the product, we use Leibniz rule: +$$\frac{d}{dx}\prod_{i=a}^b f(i,x) = + \sum_{i=a}^b \frac{\frac{d}{dx} f(i,x)}{f(i,x)}\prod_{i=a}^b f(i,x) +$$ + +There is one situation where this definition might produce wrong results, +namely when the product is zero, but axiom failed to recognize it: in this +case, +$$ + \frac{d}{dx} f(i,x)/f(i,x) +$$ +is undefined for some $i$. However, I was not able to come up with an +example. The alternative definition +$$ + \frac{d}{dx}\prod_{i=a}^b f(i,x) = + \sum_{i=a}^b \left(\frac{d}{dx} f(i,x)\right)\prod_{j=a,j\neq i}^b f(j,x) +$$ +has the slight (display) problem that we would have to come up with a new index +variable, which looks very ugly. Furthermore, it seems to me that more +simplifications will occur with the first definition. + +\begin{verbatim} + f := operator 'f + D(product(f(i,x),i=1..m),x) +\end{verbatim} + +\subsubsection{dvpow2} +The dvpow2 operation implements the differentiation of the power operator +\verb|%power| with respect to its second argument, i.e., the exponent. +It uses the formula +$$\frac{d}{dx} g(y)^x = \frac{d}{dx} e^{x\log g(y)} = \log g(y) g(y)^x.$$ + +If $g(y)$ equals zero, this formula is not valid, since the logarithm is not +defined there. Although strictly speaking $0^x$ is not differentiable at zero, +we return zero for convenience. + <>= +)abbrev package COMBF CombinatorialFunction +++ Provides the usual combinatorial functions +++ Author: Manuel Bronstein, Martin Rubey +++ Date Created: 2 Aug 1988 +++ Date Last Updated: 30 October 2005 +++ Description: +++ Provides combinatorial functions over an integral domain. +++ Keywords: combinatorial, function, factorial. +++ Examples: )r COMBF INPUT + +CombinatorialFunction(R, F): Exports == Implementation where + R: Join(OrderedSet, IntegralDomain) + F: FunctionSpace R + + OP ==> BasicOperator + K ==> Kernel F + SE ==> Symbol + O ==> OutputForm + SMP ==> SparseMultivariatePolynomial(R, K) + Z ==> Integer + + POWER ==> "%power"::Symbol + OPEXP ==> "exp"::Symbol + SPECIALDIFF ==> "%specialDiff" + SPECIALDISP ==> "%specialDisp" + SPECIALEQUAL ==> "%specialEqual" + + Exports ==> with + belong? : OP -> Boolean + ++ belong?(op) is true if op is a combinatorial operator; + operator : OP -> OP + ++ operator(op) returns a copy of op with the domain-dependent + ++ properties appropriate for F; + ++ error if op is not a combinatorial operator; + "**" : (F, F) -> F + ++ a ** b is the formal exponential a**b; + binomial : (F, F) -> F + ++ binomial(n, r) returns the number of subsets of r objects + ++ taken among n objects, i.e. n!/(r! * (n-r)!); + ++ + ++X [binomial(5,i) for i in 0..5] permutation: (F, F) -> F ++ permutation(n, r) returns the number of permutations of ++ n objects taken r at a time, i.e. n!/(n-r)!; @@ -6503,25 +6706,14 @@ here. K2fact : (K, List SE) -> F smpfact : (SMP, List SE) -> F - dummy == new()$SE :: F -@ -This macro will be used in [[product]] and [[summation]], both the $5$ and $3$ -argument forms. It is used to introduce a dummy variable in place of the -summation index within the summands. This in turn is necessary to keep the -indexing variable local, circumventing problems, for example, with -differentiation. +-- This macro will be used in product and summation, both the 5 and 3 +-- argument forms. It is used to introduce a dummy variable in place of the +-- summation index within the summands. This in turn is necessary to keep the +-- indexing variable local, circumventing problems, for example, with +-- differentiation. -This works if we don't accidently use such a symbol as a bound of summation or -product. - -Note that up to [[patch--25]] this used to read -\begin{verbatim} - dummy := new()$SE :: F -\end{verbatim} -thus introducing the same dummy variable for all products and summations, which -caused nested products and summations to fail. (Issue~\#72) + dummy == new()$SE :: F -<>= opfact := operator("factorial"::Symbol)$CommonOperators opperm := operator("permutation"::Symbol)$CommonOperators opbinom := operator("binomial"::Symbol)$CommonOperators @@ -6576,11 +6768,9 @@ caused nested products and summations to fail. (Issue~\#72) dm := dummy opsum [eval(x, k := kernel(i)$K, dm), dm, k::F] -@ -These two operations return the product or the sum as unevaluated operators. A -dummy variable is introduced to make the indexing variable \lq local\rq. +-- These two operations return the product or the sum as unevaluated operators +-- A dummy variable is introduced to make the indexing variable local. -<>= dvsum(l, x) == opsum [differentiate(first l, x), second l, third l] @@ -6592,51 +6782,6 @@ dummy variable is introduced to make the indexing variable \lq local\rq. else opdsum [differentiate(first l, x), second l, y, g, h] -@ -The above two operations implement differentiation of sums with and without -bounds. Note that the function -$$n\mapsto\sum_{k=1}^n f(k,n)$$ -is well defined only for integral values of $n$ greater than or equal to zero. -There is not even consensus how to define this function for $n<0$. Thus, it is -not differentiable. Therefore, we need to check whether we erroneously are -differentiating with respect to the upper bound or the lower bound, where the -same reasoning holds. - -Differentiating a sum with respect to its indexing variable correctly gives -zero. This is due to the introduction of dummy variables in the internal -representation of a sum: the operator [[%defsum]] takes 5 arguments, namely - -\begin{enumerate} -\item the summands, where each occurrence of the indexing variable is replaced - by -\item the dummy variable, -\item the indexing variable, -\item the lower bound, and -\item the upper bound. -\end{enumerate} - -Note that up to [[patch--40]] the following incorrect code was used, which -tried to parallel the known rules for integration: (Issue~\#180) - -\begin{verbatim} - dvdsum(l, x) == - x = retract(y := third l)@SE => 0 - k := retract(d := second l)@K - differentiate(h := third rest rest l,x) * eval(f := first l, k, h) - - differentiate(g := third rest l, x) * eval(f, k, g) - + opdsum [differentiate(f, x), d, y, g, h] -\end{verbatim} - -Up to [[patch--45]] a similar mistake could be found in the code for -differentiation of formal sums, which read -\begin{verbatim} - dvsum(l, x) == - k := retract(second l)@K - differentiate(third l, x) * summand l - + opsum [differentiate(first l, x), second l, third l] -\end{verbatim} - -<>= dvprod(l, x) == dm := retract(dummy)@SE f := eval(first l, retract(second l)@K, dm::F) @@ -6653,42 +6798,9 @@ differentiation of formal sums, which read else opdsum cons(differentiate(first l, x)/first l, rest l) * opdprod l -@ -The above two operations implement differentiation of products with and without -bounds. Note again, that we cannot even properly define products with bounds -that are not integral. - -To differentiate the product, we use Leibniz rule: -$$\frac{d}{dx}\prod_{i=a}^b f(i,x) = - \sum_{i=a}^b \frac{\frac{d}{dx} f(i,x)}{f(i,x)}\prod_{i=a}^b f(i,x) -$$ - -There is one situation where this definition might produce wrong results, -namely when the product is zero, but axiom failed to recognize it: in this -case, -$$ - \frac{d}{dx} f(i,x)/f(i,x) -$$ -is undefined for some $i$. However, I was not able to come up with an -example. The alternative definition -$$ - \frac{d}{dx}\prod_{i=a}^b f(i,x) = - \sum_{i=a}^b \left(\frac{d}{dx} f(i,x)\right)\prod_{j=a,j\neq i}^b f(j,x) -$$ -has the slight (display) problem that we would have to come up with a new index -variable, which looks very ugly. Furthermore, it seems to me that more -simplifications will occur with the first definition. - -<>= - f := operator 'f - D(product(f(i,x),i=1..m),x) -@ - -Note that up to [[patch--45]] these functions did not exist and products were -differentiated according to the usual chain rule, which gave incorrect -results. (Issue~\#211) +-- These four operations handle the conversion of sums and products to +-- OutputForm -<>= dprod l == prod(summand(l)::O, third(l)::O) @@ -6701,33 +6813,26 @@ results. (Issue~\#211) ddsum l == sum(summand(l)::O, third(l)::O = fourth(l)::O, fourth(rest l)::O) -@ -These four operations handle the conversion of sums and products to -[[OutputForm]]. Note that up to [[patch--45]] the definitions for sums and -products without bounds were missing and output was illegible. +-- The two operations handle the testing for equality of sums and products. +-- The corresponding property \verb|%specialEqual| set below is checked in +-- Kernel. Note that we can assume that the operators are equal, since this is +-- checked in Kernel itself. -<>= equalsumprod(s1, s2) == l1 := argument s1 l2 := argument s2 - (eval(first l1, retract(second l1)@K, second l2) = first l2) equaldsumprod(s1, s2) == l1 := argument s1 l2 := argument s2 - ((third rest l1 = third rest l2) and (third rest rest l1 = third rest rest l2) and (eval(first l1, retract(second l1)@K, second l2) = first l2)) -@ -The preceding two operations handle the testing for equality of sums and -products. This functionality was missing up to [[patch--45]]. (Issue~\#213) The -corresponding property [[%specialEqual]] set below is checked in -[[Kernel]]. Note that we can assume that the operators are equal, since this is -checked in [[Kernel]] itself. -<>= +-- These two operations return the product or the sum as unevaluated operators +-- A dummy variable is introduced to make the indexing variable local. + product(x:F, s:SegmentBinding F) == k := kernel(variable s)$K dm := dummy @@ -6738,11 +6843,6 @@ checked in [[Kernel]] itself. dm := dummy opdsum [eval(x,k,dm), dm, k::F, lo segment s, hi segment s] -@ -These two operations return the product or the sum as unevaluated operators. A -dummy variable is introduced to make the indexing variable \lq local\rq. - -<>= smpfact(p, l) == map(K2fact(#1, l), #1::F, p)$PolynomialCategoryLifting( IndexedExponents K, K, R, SMP, F) @@ -6889,15 +6989,9 @@ dummy variable is introduced to make the indexing variable \lq local\rq. => ibinom l binomial(r1::R, r2::R)::F -@ - -[[iibinom]] checks those cases in which the binomial coefficient may be -evaluated explicitly. Note that up to [[patch--51]], the case where the second -argument is a positive integer was not checked.(Issue~\#336) Currently, the -naive iterative algorithm is used to calculate the coefficient, there is room -for improvement here. - -<>= +-- iibinom checks those cases in which the binomial coefficient may +-- be evaluated explicitly. Currently, the naive iterative algorithm is +-- used to calculate the coefficient, there is room for improvement here. else iibinom l == @@ -6926,23 +7020,6 @@ for improvement here. else log(first l) * first(l) ** second(l) -@ -This operation implements the differentiation of the power operator [[%power]] -with respect to its second argument, i.e., the exponent. It uses the formula -$$\frac{d}{dx} g(y)^x = \frac{d}{dx} e^{x\log g(y)} = \log g(y) g(y)^x.$$ - -If $g(y)$ equals zero, this formula is not valid, since the logarithm is not -defined there. Although strictly speaking $0^x$ is not differentiable at zero, -we return zero for convenience. - -Note that up to [[patch--25]] this used to read -\begin{verbatim} - if F has ElementaryFunctionCategory then - dvpow2 l == log(first l) * first(l) ** second(l) -\end{verbatim} -which caused differentiating $0^x$ to fail. (Issue~\#19) - -<>= evaluate(opfact, iifact)$BasicOperatorFunctions1(F) evaluate(oppow, iipow) evaluate(opperm, iiperm) @@ -6952,17 +7029,18 @@ which caused differentiating $0^x$ to fail. (Issue~\#19) evaluate(opprod, iprod) evaluate(opdprod, iidprod) derivative(oppow, [dvpow1, dvpow2]) + +-- These four properties define special differentiation rules for sums and +-- products. + setProperty(opsum, SPECIALDIFF, dvsum@((List F, SE) -> F) pretend None) setProperty(opdsum, SPECIALDIFF, dvdsum@((List F, SE)->F) pretend None) setProperty(opprod, SPECIALDIFF, dvprod@((List F, SE)->F) pretend None) setProperty(opdprod, SPECIALDIFF, dvdprod@((List F, SE)->F) pretend None) -@ -The last four properties define special differentiation rules for sums and -products. Note that up to [[patch--45]] the rules for products were missing. -Thus products were differentiated according the usual chain-rule, which gave -incorrect results. -<>= +-- Set the properties for displaying sums and products and testing for +-- equality. + setProperty(opsum, SPECIALDISP, dsum@(List F -> O) pretend None) setProperty(opdsum, SPECIALDISP, ddsum@(List F -> O) pretend None) setProperty(opprod, SPECIALDISP, dprod@(List F -> O) pretend None) @@ -6973,10 +7051,6 @@ incorrect results. setProperty(opdprod, SPECIALEQUAL, equaldsumprod@((K,K) -> Boolean) pretend None) @ -Finally, we set the properties for displaying sums and products and testing for -equality. - - <>= "COMBF" [color="#FF4488",href="bookvol10.4.pdf#nameddest=COMBF"] "FS" [color="#4488FF",href="bookvol10.2.pdf#nameddest=FS"] @@ -43431,6 +43505,107 @@ IntegerBits: with @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{package COMBINAT IntegerCombinatoricFunctions} +<>= +)set break resume +)sys rm -f IntegerCombinatoricFunctions.output +)spool IntegerCombinatoricFunctions.output +)set message test on +)set message auto off +)clear all +--S 1 of 4 +)set expose add constructor OutputForm +--R +--I OutputForm is already explicitly exposed in frame frame0 +--E 1 + +--S 2 of 4 +pascalRow(n) == [right(binomial(n,i),4) for i in 0..n] +--R +--R Type: Void +--E 2 + +--S 3 of 4 +displayRow(n)==output center blankSeparate pascalRow(n) +--R +--R Type: Void +--E 3 + +--S 4 of 4 +for i in 0..7 repeat displayRow i +--R +--R Compiling function pascalRow with type NonNegativeInteger -> List +--R OutputForm +--R Compiling function displayRow with type NonNegativeInteger -> Void +--R 1 +--R 1 1 +--R 1 2 1 +--R 1 3 3 1 +--R 1 4 6 4 1 +--R 1 5 10 10 5 1 +--R 1 6 15 20 15 6 1 +--R 1 7 21 35 35 21 7 1 +--R Type: Void +--E 4 + +)spool +)lisp (bye) +@ +<>= +==================================================================== +IntegerCombinatoricFunctions examples +==================================================================== + +The binomial(n, r) returns the number of subsets of r objects +taken among n objects, i.e. n!/(r! * (n-r)!) + +The binomial coefficients are the coefficients of the series expansion +of a power of a binomial, that is + + n + --+ / n \ k n + > | | x = (1 + x) + --+ \ k / + k= 0 + +This leads to the famous pascal triangle. First we expose the OutputForm +domain, which is normally hidden, so we can use it to format the lines. + + )set expose add constructor OutputForm + +Next we define a function that will output the list of binomial coefficients +right justified with proper spacing: + + pascalRow(n) == [right(binomial(n,i),4) for i in 0..n] + +and now we format the whole line so that it looks centered: + + displayRow(n)==output center blankSeparate pascalRow(n) + +and we compute the triangle + + for i in 0..7 repeat displayRow i + +giving the pretty result: + + Compiling function pascalRow with type NonNegativeInteger -> List + OutputForm + Compiling function displayRow with type NonNegativeInteger -> Void + 1 + 1 1 + 1 2 1 + 1 3 3 1 + 1 4 6 4 1 + 1 5 10 10 5 1 + 1 6 15 20 15 6 1 + 1 7 21 35 35 21 7 1 + +See Also: +o )show IntegerCombinatoricFunctions +o )d op binomial +o )show OutputForm +o )help set + +@ \pagehead{IntegerCombinatoricFunctions}{COMBINAT} \pagepic{ps/v104integercombinatoricfunctions.ps}{COMBINAT}{1.00} @@ -43469,6 +43644,8 @@ IntegerCombinatoricFunctions(I:IntegerNumberSystem): with ++ \spad{binomial(n,r)} returns the binomial coefficient ++ \spad{C(n,r) = n!/(r! (n-r)!)}, where \spad{n >= r >= 0}. ++ This is the number of combinations of n objects taken r at a time. + ++ + ++X [binomial(5,i) for i in 0..5] factorial: I -> I ++ \spad{factorial(n)} returns \spad{n!}. this is the product of all ++ integers between 1 and n (inclusive). @@ -123791,7 +123968,7 @@ PolynomialNumberTheoryFunctions(): Exports == Implementation where Exports ==> with bernoulli : I -> SUP RN ++ bernoulli(n) returns the nth Bernoulli polynomial \spad{B[n](x)}. - ++ Note: Bernoulli polynomials denoted \spad{B(n,x)} computed by solving the + ++ Bernoulli polynomials denoted \spad{B(n,x)} computed by solving the ++ differential equation \spad{differentiate(B(n,x),x) = n B(n-1,x)} where ++ \spad{B(0,x) = 1} and initial condition comes from \spad{B(n) = B(n,0)}. chebyshevT: I -> SUP I @@ -129823,49 +130000,45 @@ RandomFloatDistributions(): Cat == Body where ++ Description: ++ This package exports integer distributions RandomIntegerDistributions(): with - uniform: Segment Integer -> (() -> Integer) - ++ uniform(s) \undocumented - binomial: (Integer, RationalNumber) -> (() -> Integer) - ++ binomial(n,f) \undocumented - poisson: RationalNumber -> (() -> Integer) - ++ poisson(f) \undocumented - geometric: RationalNumber -> (() -> Integer) - ++ geometric(f) \undocumented - - ridHack1: (Integer,Integer,Integer,Integer) -> Integer - ++ ridHack1(i,j,k,l) \undocumented - == add - import RandomNumberSource() - import IntegerBits() - - -- Compute uniform(a..b) as - -- - -- l + U0 + w*U1 + w**2*U2 +...+ w**(n-1)*U-1 + w**n*M - -- - -- where - -- l = min(a,b) - -- m = abs(b-a) + 1 - -- w**n < m < w**(n+1) - -- U0,...,Un-1 are uniform on 0..w-1 - -- M is uniform on 0..(m quo w**n)-1 - - uniform aTob == - a := lo aTob; b := hi aTob - l := min(a,b); m := abs(a-b) + 1 - - w := 2**(bitLength size() quo 2)::NonNegativeInteger - - n := 0 - mq := m -- m quo w**n - while (mqnext := mq quo w) > 0 repeat - n := n + 1 - mq := mqnext - ridHack1(mq, n, w, l) - - ridHack1(mq, n, w, l) == - r := randnum mq - for i in 1..n repeat r := r*w + randnum w - r + l + uniform: Segment Integer -> (() -> Integer) + ++ uniform(s) as + ++ l + u0 + w*u1 + w**2*u2 +...+ w**(n-1)*u-1 + w**n*m + ++ where + ++ s = a..b + ++ l = min(a,b) + ++ m = abs(b-a) + 1 + ++ w**n < m < w**(n+1) + ++ u0,...,un-1 are uniform on 0..w-1 + ++ m is uniform on 0..(m quo w**n)-1 + binomial: (Integer, RationalNumber) -> (() -> Integer) + ++ binomial(n,f) \undocumented + poisson: RationalNumber -> (() -> Integer) + ++ poisson(f) \undocumented + geometric: RationalNumber -> (() -> Integer) + ++ geometric(f) \undocumented + ridHack1: (Integer,Integer,Integer,Integer) -> Integer + ++ ridHack1(i,j,k,l) \undocumented + == add + import RandomNumberSource() + import IntegerBits() + + uniform aTob == + a := lo aTob; b := hi aTob + l := min(a,b); m := abs(a-b) + 1 + + w := 2**(bitLength size() quo 2)::NonNegativeInteger + + n := 0 + mq := m -- m quo w**n + while (mqnext := mq quo w) > 0 repeat + n := n + 1 + mq := mqnext + ridHack1(mq, n, w, l) + + ridHack1(mq, n, w, l) == + r := randnum mq + for i in 1..n repeat r := r*w + randnum w + r + l @ <>= diff --git a/changelog b/changelog index 7e8e2d2..fed1b33 100644 --- a/changelog +++ b/changelog @@ -1,3 +1,7 @@ +20090417 tpd src/axiom-website/patches.html 20090417.01.tpd.patch +20090417 tpd src/algebra/Makefile add help, regression tests +20090417 tpd books/bookvol10.4 document binomial +20090417 tpd books/bookvol10.2 document binomial 20090416 tpd src/axiom-website/patches.html 20090416.03.tpd.patch 20090416 tpd books/bookvol10.4 update bezier documentation 20090416 tpd books/ps/v104beziermove.eps graph the bezier curves diff --git a/src/algebra/Makefile.pamphlet b/src/algebra/Makefile.pamphlet index f454c55..8af306b 100644 --- a/src/algebra/Makefile.pamphlet +++ b/src/algebra/Makefile.pamphlet @@ -16445,6 +16445,7 @@ SPADHELP=\ ${HELP}/BinarySearchTree.help ${HELP}/CardinalNumber.help \ ${HELP}/CartesianTensor.help ${HELP}/Character.help \ ${HELP}/CharacterClass.help ${HELP}/CliffordAlgebra.help \ + ${HELP}/CombinatorialFunction.help \ ${HELP}/Complex.help ${HELP}/ContinuedFraction.help \ ${HELP}/CycleIndicators.help ${HELP}/DeRhamComplex.help \ ${HELP}/DecimalExpansion.help ${HELP}/Dequeue.help \ @@ -16463,7 +16464,9 @@ SPADHELP=\ ${HELP}/GroebnerPackage.help \ ${HELP}/Heap.help ${HELP}/HexadecimalExpansion.help \ ${HELP}/HomogeneousDistributedMultivariatePolynomial.help \ - ${HELP}/Integer.help ${HELP}/IntegerLinearDependence.help \ + ${HELP}/Integer.help \ + ${HELP}/IntegerCombinatoricFunctions.help \ + ${HELP}/IntegerLinearDependence.help \ ${HELP}/IntegerNumberTheoryFunctions.help \ ${HELP}/Kernel.help ${HELP}/KeyedAccessFile.help \ ${HELP}/LazardSetSolvingPackage.help \ @@ -16536,6 +16539,7 @@ REGRESS=\ BinarySearchTree.regress CardinalNumber.regress \ CartesianTensor.regress Character.regress \ CharacterClass.regress CliffordAlgebra.regress \ + CombinatorialFunction.regress \ Complex.regress ContinuedFraction.regress \ CycleIndicators.regress DeRhamComplex.regress \ DecimalExpansion.regress Dequeue.regress \ @@ -16553,7 +16557,8 @@ REGRESS=\ GroebnerPackage.regress \ Heap.regress HexadecimalExpansion.regress \ HomogeneousDistributedMultivariatePolynomial.regress \ - Integer.regress IntegerLinearDependence.regress \ + Integer.regress IntegerCombinatoricFunctions.regress \ + IntegerLinearDependence.regress \ IntegerNumberTheoryFunctions.regress \ Kernel.regress KeyedAccessFile.regress \ LazardSetSolvingPackage.regress \ @@ -16744,6 +16749,18 @@ ${HELP}/CliffordAlgebra.help: ${BOOKS}/bookvol10.3.pamphlet >${INPUT}/CliffordAlgebra.input @echo "CliffordAlgebra (CLIF)" >>${HELPFILE} +${HELP}/CombinatorialFunction.help: ${BOOKS}/bookvol10.4.pamphlet + @echo 7010 create CombinatorialFunction.help from \ + ${BOOKS}/bookvol10.4.pamphlet + @${TANGLE} -R"CombinatorialFunction.help" \ + ${BOOKS}/bookvol10.4.pamphlet \ + >${HELP}/CombinatorialFunction.help + @cp ${HELP}/CombinatorialFunction.help ${HELP}/COMBF.help + @${TANGLE} -R"CombinatorialFunction.input" \ + ${BOOKS}/bookvol10.4.pamphlet \ + >${INPUT}/CombinatorialFunction.input + @echo "CombinatorialFunction (COMBF)" >>${HELPFILE} + ${HELP}/Complex.help: ${BOOKS}/bookvol10.3.pamphlet @echo 7011 create Complex.help from ${BOOKS}/bookvol10.3.pamphlet @${TANGLE} -R"Complex.help" ${BOOKS}/bookvol10.3.pamphlet \ @@ -17028,6 +17045,18 @@ ${HELP}/Integer.help: ${BOOKS}/bookvol10.3.pamphlet >${INPUT}/Integer.input @echo "Integer (INT)" >>${HELPFILE} +${HELP}/IntegerCombinatoricFunctions.help: ${BOOKS}/bookvol10.4.pamphlet + @echo 7039 create IntegerCombinatoricFunctions.help from \ + ${BOOKS}/bookvol10.4.pamphlet + @${TANGLE} -R"IntegerCombinatoricFunctions.help" \ + ${BOOKS}/bookvol10.4.pamphlet \ + >${HELP}/IntegerCombinatoricFunctions.help + @cp ${HELP}/IntegerCombinatoricFunctions.help ${HELP}/COMBINAT.help + @${TANGLE} -R"IntegerCombinatoricFunctions.input" \ + ${BOOKS}/bookvol10.4.pamphlet \ + >${INPUT}/IntegerCombinatoricFunctions.input + @echo "IntegerCombinatoricFunctions (COMBINAT)" >>${HELPFILE} + ${HELP}/IntegerLinearDependence.help: ${BOOKS}/bookvol10.4.pamphlet @echo 7039 create IntegerLinearDependence.help from \ ${BOOKS}/bookvol10.4.pamphlet diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html index 758d24e..d8cd037 100644 --- a/src/axiom-website/patches.html +++ b/src/axiom-website/patches.html @@ -1094,5 +1094,7 @@ regress.lisp tighten checks on regression tests
bookvol5 move more interpreter code
20090416.03.tpd.patch bookvol10.4 update bezier documentation
+20090416.03.tpd.patch +bookvol10.4, 10.2 document binomial