1 ;; Author Barton Willis
2 ;; University of Nebraska at Kearney
3 ;; Copyright (C) 2011,2021 Barton Willis
5 ;; This program is free software; you can redistribute it and/or modify
6 ;; it under the terms of the GNU General Public License as published by
7 ;; the Free Software Foundation; either version 2 of the License, or
8 ;; (at your option) any later version.
10 ;; This program is distributed in the hope that it will be useful,
11 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
12 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 ;; GNU General Public License for more details.
17 ($load
"bernstein_utilities.mac")
19 ;; When bernstein_explicit is non-nil, bernstein_poly(k,n,x) simplifies to
20 ;; binomial(n,k) * x^k (1-x)^(n-k) regardless of the values of k or n;
22 (defmvar $bernstein_explicit nil
)
24 ;; numerical (complex rational, float, or big float) evaluation of bernstein polynomials
25 (in-package #:bigfloat
)
27 (defun bernstein-poly (k n x
)
28 (* (to (maxima::opcons
'maxima
::%binomial n k
)) (expt x k
) (expt (- 1 x
) (- n k
))))
32 (defun $bernstein_poly
(k n x
) (simplify (list '(%bernstein_poly
) k n x
)))
33 (defprop $bernstein_poly %bernstein_poly alias
)
34 (defprop $bernstein_poly %bernstein_poly verb
)
35 (defprop %bernstein_poly $bernstein_poly reversealias
)
36 (defprop %bernstein_poly $bernstein_poly noun
)
38 (setf (get '%bernstein_poly
'conjugate-function
)
43 (if (and ($featurep k
'$integer
) ($featurep n
'$integer
))
44 (opcons '%bernstein_poly k n
(opcons '$conjugate x
))
45 (list (list '$conjugate
'simp
) (take '(%bernstein_poly
) k n x
))))))
47 ;; integrate(bernstein_poly(k,n,x),x) = hypergeometric([k+1,k-n],[k+2],x)*x^(k+1)/(k+1)
49 (defun bernstein-integral (k n x
)
52 (opcons '%binomial n k
)
53 (opcons 'mexpt x
(add 1 k
))
54 (opcons '%hypergeometric
55 (opcons 'mlist
(add 1 k
) (sub k n
))
56 (opcons 'mlist
(add 2 k
))
60 (putprop '%bernstein_poly
`((k n x
) nil nil
,'bernstein-integral
) 'integral
)
61 (putprop '$bernstein_poly
`((k n x
) nil nil
,'bernstein-integral
) 'integral
)
63 (defun bernstein-poly-simp (e y z
)
65 (let* ((fn (car (pop e
)))
66 (k (if (consp e
) (simpcheck (pop e
) z
) (wna-err fn
)))
67 (n (if (consp e
) (simpcheck (pop e
) z
) (wna-err fn
)))
68 (x (if (consp e
) (simpcheck (pop e
) z
) (wna-err fn
))))
69 (if (consp e
) (wna-err fn
))
70 (cond ((and (integerp k
) (integerp n
) (>= k
0) (>= n k
)
71 (complex-number-p x
#'(lambda (s) (or (integerp s
) ($ratnump s
) (floatp s
) ($bfloatp s
)))))
72 (maxima::to
(bigfloat::bernstein-poly
(bigfloat::to k
) (bigfloat::to n
) (bigfloat::to x
))))
74 ((zerop1 x
) (opcons '%kron_delta k
0))
76 ((onep1 x
) (opcons '%kron_delta k n
))
78 ((or $bernstein_explicit
(and (integerp k
) (integerp n
)))
79 (if (and (integerp k
) (integerp n
) (or (< k
0) (> k n
))) (mul 0 x
)
80 (mul (opcons '%binomial n k
) (opcons 'mexpt x k
) (opcons 'mexpt
(sub 1 x
) (sub n k
)))))
82 (t (list (list fn
'simp
) k n x
)))))
84 (setf (get '%bernstein_poly
'operators
) 'bernstein-poly-simp
)
86 (defprop %bernstein_poly
89 ((mtimes) ((%bernstein_poly
) k n x
)
92 ((mqapply) (($psi array
) 0) ((mplus) 1 k
)))
93 ((mqapply) (($psi array
) 0)
94 ((mplus) 1 ((mtimes) -
1 k
) n
))
96 ((%log
) ((mplus) 1 ((mtimes) -
1 x
))))
99 ((mtimes)((%bernstein_poly
) k n x
)
101 ((mqapply) (($psi array
) 0) ((mplus) 1 n
))
103 ((mqapply) (($psi array
) 0)
104 ((mplus) 1 ((mtimes) -
1 k
) n
)))
105 ((%log
) ((mplus) 1 ((mtimes) -
1 x
)))))
109 ((%bernstein_poly
) ((mplus) -
1 k
) ((mplus) -
1 n
) x
)
111 ((%bernstein_poly
) k
((mplus) -
1 n
) x
))) n
))
114 (defun $bernstein_approx
(e vars n
)
115 (if (not ($listp vars
)) (merror "The second argument to bernstein_approx must be a list"))
117 (setq vars
(margs vars
))
118 (if (some #'(lambda (s) (not ($mapatom s
))) vars
)
119 (merror "The second argument to bernstein_approx must be a list of atoms"))
121 (if (or (not (integerp n
)) (< n
1)) (merror "The third argument to bernstein_approx must be a positive integer"))
123 (let* ((k (length vars
))
124 (d (make-list k
:initial-element
0))
125 (nn (make-list k
:initial-element n
))
126 (carry) (acc 0) (m) (x))
128 (setq m
(expt (+ n
1) k
))
129 (setq nn
(cons '(mlist) nn
))
133 (opcons '%multibernstein_poly
(cons '(mlist) d
) nn
(cons '(mlist) vars
))
134 ($substitute
(cons '(mlist) (mapcar #'(lambda (s x
) (opcons 'mequal x
(div s n
))) d vars
)) e
))))
137 (setq x
(+ carry
(nth j d
)))
138 (if (> x n
) (setq x
0 carry
1) (setq carry
0))
142 (defun $multibernstein_poly
(k n x
) (simplify (list '(%multibernstein_poly
) k n x
)))
143 (defprop $multibernstein_poly %multibernstein_poly alias
)
144 (defprop $multibernstein_poly %multibernstein_poly verb
)
145 (defprop %multibernstein_poly $multibernstein_poly reversealias
)
146 (defprop %multibernstein_poly $multibernstein_poly noun
)
148 (defun multi-bernstein-poly-simp (e y z
)
152 (k (if (consp e
) (simpcheck (pop e
) z
) (wna-err fn
)))
153 (n (if (consp e
) (simpcheck (pop e
) z
) (wna-err fn
)))
154 (x (if (consp e
) (simpcheck (pop e
) z
) (wna-err fn
))))
155 (if (consp e
) (wna-err fn
))
157 (if (or (not (and ($listp k
) ($listp n
) ($listp x
))) (/= ($length k
) ($length n
)) (/= ($length n
) ($length x
)))
158 (merror "Each argument to multibernstein_poly must be an equal length list"))
160 (muln (mapcar #'(lambda (a b z
) (opcons '%bernstein_poly a b z
)) (margs k
) (margs n
) (margs x
)) t
)))
162 (setf (get '%multibernstein_poly
'operators
) 'multi-bernstein-poly-simp
)
