1 /*
2 * optimization.cc
3 *
4 * Copyright 2017 Luka Marohnić
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 *
20 * __________________________________________________________________________
21 * |Example of using 'extrema', 'minimize' and 'maximize' functions to solve|
22 * |the set of exercises found in: |
23 * |https://math.feld.cvut.cz/habala/teaching/veci-ma2/ema2r3.pdf |
24 * ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
25 * 1) Input:
26 * extrema(2x^3+9x*y^2+15x^2+27y^2,[x,y])
27 * Result: we get local minimum at origin, local maximum at (-5,0)
28 * and (-3,2),(-3,-2) as saddle points.
29 *
30 * 2) Input:
31 * extrema(x^3-2x^2+y^2+z^2-2x*y+x*z-y*z+3z,[x,y,z])
32 * Result: the given function has local minimum at (2,1,-2).
33 *
34 * 3) Input:
35 * minimize(x^2+2y^2,x^2-2x+2y^2+4y=0,[x,y]);
36 * maximize(x^2+2y^2,x^2-2x+2y^2+4y=0,[x,y])
37 * Result: the minimal value of x^2+2y^2 is 0 and the maximal is 12.
38 *
39 * 4) We need to minimize f(x,y,z)=x^2+(y+3)^2+(z-2)^2 for points (x,y,z)
40 * lying in plane x+y-z=1. Since the feasible area is not bounded, we're
41 * using the function 'extrema' because obviously the function has single
42 * local minimum.
43 * Input:
44 * extrema(x^2+(y+3)^2+(z-2)^2,x+y-z=1,[x,y,z])
45 * Result: the point closest to P in x+y-z=1 is (2,-1,0), and the distance
46 * is equal to sqrt(f(2,-1,0))=2*sqrt(3).
47 *
48 * 5) We're using the same method as in exercise 4.
49 * Input:
50 * extrema((x-1)^2+(y-2)^2+(z+1)^2,[x+y+z=1,2x-y+z=3],[x,y,z])
51 * Result: the closest point is (2,0,-1) and the corresponding distance
52 * equals to sqrt(5).
53 *
54 * 6) First we need to determine the feasible area. Plot its bounds with:
55 * implicitplot(x^2+(y+1)^2-4,x,y);line(y=-1);line(y=x+1)
56 * Now we see that the feasible area is given by set of inequalities:
57 * cond:=[x^2+(y+1)^2<=4,y>=-1,y<=x+1]
58 * Draw this area with:
59 * plotinequation(cond,[x,y],xstep=0.05,ystep=0.05)
60 * Now calculate global minimum and maximum of f(x,y)=x^2+4y^2 on that area:
61 * f(x,y):=x^2+4y^2;
62 * minimize(f(x,y),cond,[x,y]);maximize(f(x,y),cond,[x,y])
63 * Result: the minimum is 0 and the maximum is 8.
64 *
65 * 7) Input:
66 * minimize(x^2+y^2-6x+6y,x^2+y^2<=4,[x,y]);
67 * maximize(x^2+y^2-6x+6y,x^2+y^2<=4,[x,y])
68 * Result: the minimum is 4-12*sqrt(2) and the maximum is 4+12*sqrt(2).
69 *
70 * 8) Input:
71 * extrema(y,y^2+2x*y=2x-4x^2,[x,y])
72 * Result: we obtain (1/2,-1) as local minimum and (1/6,1/3) as local
73 * maximum of f(x,y)=y. Therefore, the maximal value is y(1/2)=-1
74 * and the maximal value is y(1/6)=1/3.
75 *
76 * The above set of exercises could be turned into an example Xcas worksheet.
77 */
78
79 #include "giacPCH.h"
80 #include "giac.h"
81 #include "optimization.h"
82 #include "signalprocessing.h"
83 #include <sstream>
84 #include <bitset>
85
86 using namespace std;
87
88 #ifndef NO_NAMESPACE_GIAC
89 namespace giac {
90 #endif // ndef NO_NAMESPACE_GIAC
91
92 #define GOLDEN_RATIO 1.61803398875
93 typedef unsigned long ulong;
94
95 /* simplify the expression */
simp(const gen & g,GIAC_CONTEXT)96 gen simp(const gen &g,GIAC_CONTEXT) {
97 if (g.type==_VECT) {
98 vecteur res;
99 for (const_iterateur it=g._VECTptr->begin();it!=g._VECTptr->end();++it) {
100 res.push_back(simp(*it,contextptr));
101 }
102 return res;
103 }
104 if (_evalf(g,contextptr).type==_DOUBLE_)
105 return simplify(g,contextptr);
106 return recursive_normal(g,contextptr);
107 }
108
109 /*
110 * Return true iff the expression 'e' is rational with respect to
111 * variables in 'vars'.
112 */
is_rational_wrt_vars(const gen & e,const vecteur & vars,GIAC_CONTEXT)113 bool is_rational_wrt_vars(const gen &e,const vecteur &vars,GIAC_CONTEXT) {
114 for (const_iterateur it=vars.begin();it!=vars.end();++it) {
115 vecteur l(rlvarx(e,*it));
116 if (l.size()>1)
117 return false;
118 }
119 return true;
120 }
121
122 /* simplify the critical points */
cpt_simp(vecteur & cv,const vecteur & vars,const gen & f,GIAC_CONTEXT)123 void cpt_simp(vecteur &cv,const vecteur &vars,const gen &f,GIAC_CONTEXT) {
124 for(int j=cv.size();j-->0;) {
125 gen sc=simp(_epsilon2zero(cv[j],contextptr),contextptr),img;
126 if (is_undef(sc) || !is_constant_wrt_vars(sc,vars,contextptr) ||
127 ((img=_evalf(im(sc,contextptr),contextptr)).type==_DOUBLE_ && !is_zero(img))) {
128 cv.erase(cv.begin()+j);
129 continue;
130 }
131 if (cv[j].type==_VECT && cv[j]._VECTptr->size()==vars.size()) {
132 gen val=simp(subst(f,vars,*cv[j]._VECTptr,false,contextptr),contextptr);
133 if (is_inf(val) || is_undef(val) || _evalf(val,contextptr).type==_CPLX) {
134 // cv[j] is not in the domain of f
135 cv.erase(cv.begin()+j);
136 continue;
137 }
138 }
139 cv[j]=simp(cv[j],contextptr);
140 }
141 }
142
solve_vect(const vecteur & e,const vecteur & v,GIAC_CONTEXT)143 vecteur solve_vect(const vecteur &e,const vecteur &v,GIAC_CONTEXT) {
144 gen tmp=_solve(makesequence(e,v),contextptr);
145 if (tmp.type!=_VECT)
146 return vecteur(0);
147 return *tmp._VECTptr;
148 }
149
150 int var_index=0;
151
is_greater_than_zero(const gen & g,const vecteur & vars,GIAC_CONTEXT)152 bool is_greater_than_zero(const gen &g,const vecteur &vars,GIAC_CONTEXT) {
153 vecteur terms(0);
154 if (g.is_symb_of_sommet(at_plus) && g._SYMBptr->feuille.type==_VECT)
155 terms=*g._SYMBptr->feuille._VECTptr;
156 else terms=makevecteur(g);
157 bool has_exp=false;
158 gen rest(0);
159 for (const_iterateur it=terms.begin();it!=terms.end();++it) {
160 if (_lin(*it,contextptr).is_symb_of_sommet(at_exp))
161 has_exp=true;
162 else rest+=*it;
163 }
164 if (!has_exp)
165 return false;
166 return is_positive(rest,contextptr);
167 }
168
remove_strictly_positive_factors(const gen & g,const vecteur & vars,GIAC_CONTEXT)169 gen remove_strictly_positive_factors(const gen &g,const vecteur &vars,GIAC_CONTEXT) {
170 gen f(g);
171 if (f.is_symb_of_sommet(at_neg))
172 f=f._SYMBptr->feuille;
173 if (f.is_symb_of_sommet(at_prod) && f._SYMBptr->feuille.type==_VECT) {
174 const vecteur &fv=*f._SYMBptr->feuille._VECTptr;
175 gen p(1);
176 for (const_iterateur jt=fv.begin();jt!=fv.end();++jt) {
177 if (is_greater_than_zero(*jt,vars,contextptr))
178 continue;
179 else p=*jt*p;
180 }
181 f=p;
182 }
183 return f;
184 }
185
186 /*
187 * Solves a system of equations.
188 * This function is based on _solve but handles cases where a variable
189 * is found inside trigonometric, hyperbolic or exponential functions.
190 */
solve2(const vecteur & e_orig,const vecteur & vars_orig,GIAC_CONTEXT)191 vecteur solve2(const vecteur &e_orig,const vecteur &vars_orig,GIAC_CONTEXT) {
192 int m=e_orig.size(),n=vars_orig.size(),i=0;
193 vecteur e_orig_simp=*expexpand(expand(_pow2exp(e_orig,contextptr),contextptr),contextptr)._VECTptr;
194 for (iterateur it=e_orig_simp.begin();it!=e_orig_simp.end();++it) {
195 if (it->is_symb_of_sommet(at_equal))
196 *it=equal2diff(*it);
197 gen f=(it->type==_SYMB?_factor(*it,contextptr):*it);
198 gen num=remove_strictly_positive_factors(_numer(f,contextptr),vars_orig,contextptr);
199 gen den=remove_strictly_positive_factors(_denom(f,contextptr),vars_orig,contextptr);
200 *it=num/den;
201 }
202 for (;i<m;++i) {
203 if (!is_rational_wrt_vars(e_orig_simp[i],vars_orig,contextptr))
204 break;
205 }
206 if (n==1 || i==m)
207 return solve_vect(e_orig_simp,vars_orig,contextptr);
208 vecteur e(*halftan(_texpand(hyp2exp(e_orig_simp,contextptr),contextptr),contextptr)._VECTptr);
209 vecteur lv(*exact(lvar(_evalf(lvar(e),contextptr)),contextptr)._VECTptr);
210 vecteur deps(n),depvars(n,gen(0));
211 vecteur vars(vars_orig);
212 const_iterateur it=lv.begin();
213 for (;it!=lv.end();++it) {
214 i=0;
215 for (;i<n;++i) {
216 if (is_undef(vars[i]))
217 continue;
218 if (*it==(deps[i]=vars[i]) ||
219 *it==(deps[i]=exp(vars[i],contextptr)) ||
220 is_zero(simp(*it-(deps[i]=tan(vars[i]/gen(2),contextptr)),contextptr))) {
221 vars[i]=undef;
222 depvars[i]=identificateur(" depvar"+print_INT_(i));
223 break;
224 }
225 }
226 if (i==n)
227 break;
228 }
229 if (it!=lv.end() || find(depvars.begin(),depvars.end(),gen(0))!=depvars.end())
230 return solve_vect(e_orig_simp,vars_orig,contextptr);
231 vecteur e_subs=subst(e,deps,depvars,false,contextptr);
232 vecteur sol=solve_vect(e_subs,depvars,contextptr);
233 vecteur ret;
234 for (const_iterateur it=sol.begin();it!=sol.end();++it) {
235 vecteur r(n);
236 i=0;
237 for (;i<n;++i) {
238 gen c(it->_VECTptr->at(i));
239 if (deps[i].type==_IDNT)
240 r[i]=c;
241 else if (deps[i].is_symb_of_sommet(at_exp) && is_strictly_positive(c,contextptr))
242 r[i]=simp(ln(c,contextptr),contextptr);
243 else if (deps[i].is_symb_of_sommet(at_tan))
244 r[i]=simp(2*atan(c,contextptr),contextptr);
245 else
246 break;
247 }
248 if (i==n)
249 ret.push_back(r);
250 }
251 return ret;
252 }
253
is_ineq_x_a(const gen & g,const gen & var,gen & a,GIAC_CONTEXT)254 bool is_ineq_x_a(const gen &g,const gen &var,gen &a,GIAC_CONTEXT) {
255 if ((g.is_symb_of_sommet(at_inferieur_egal) ||
256 g.is_symb_of_sommet(at_inferieur_strict) ||
257 g.is_symb_of_sommet(at_superieur_egal) ||
258 g.is_symb_of_sommet(at_superieur_strict)) &&
259 g._SYMBptr->feuille.type==_VECT &&
260 g._SYMBptr->feuille._VECTptr->front()==var &&
261 is_constant_wrt(g._SYMBptr->feuille._VECTptr->back(),var,contextptr))
262 {
263 a=g._SYMBptr->feuille._VECTptr->back();
264 return true;
265 }
266 return false;
267 }
268
269 /*
270 * Traverse the tree of symbolic expression 'e' and collect all points of
271 * transition in piecewise subexpressions.
272 */
collect_transition_points(const gen & e,const gen & var,vecteur & cv,GIAC_CONTEXT)273 void collect_transition_points(const gen &e,const gen &var,vecteur &cv,GIAC_CONTEXT) {
274 if (e.type==_VECT) {
275 for (const_iterateur it=e._VECTptr->begin();it!=e._VECTptr->end();++it) {
276 collect_transition_points(*it,var,cv,contextptr);
277 }
278 }
279 else if ((e.is_symb_of_sommet(at_piecewise) || e.is_symb_of_sommet(at_when)) &&
280 e._SYMBptr->feuille.type==_VECT) {
281 const vecteur &f=*e._SYMBptr->feuille._VECTptr;
282 int sz=f.size();
283 for (int i=0;i<sz/2;++i) {
284 gen g=_solve(makesequence(f[2*i],var),contextptr);
285 if (g.type==_VECT) {
286 gen a,b;
287 for (const_iterateur it=g._VECTptr->begin();it!=g._VECTptr->end();++it) {
288 if (is_ineq_x_a(*it,var,a,contextptr))
289 cv.push_back(a);
290 else if (it->is_symb_of_sommet(at_and) &&
291 it->_SYMBptr->feuille.type==_VECT &&
292 it->_SYMBptr->feuille._VECTptr->size()==2 &&
293 is_ineq_x_a(it->_SYMBptr->feuille._VECTptr->front(),var,a,contextptr) &&
294 is_ineq_x_a(it->_SYMBptr->feuille._VECTptr->back(),var,b,contextptr)) {
295 cv.push_back(a);
296 cv.push_back(b);
297 }
298 }
299 }
300 }
301 } else if (e.type==_SYMB)
302 collect_transition_points(e._SYMBptr->feuille,var,cv,contextptr);
303 }
304
next_binary_perm(vector<bool> & perm,int to_end=0)305 bool next_binary_perm(vector<bool> &perm,int to_end=0) {
306 if (to_end==int(perm.size()))
307 return false;
308 int end=int(perm.size())-1-to_end;
309 perm[end]=!perm[end];
310 return perm[end]?true:next_binary_perm(perm,to_end+1);
311 }
312
make_temp_vars(const vecteur & vars,const vecteur & ineq,bool open,GIAC_CONTEXT)313 vecteur make_temp_vars(const vecteur &vars,const vecteur &ineq,bool open,GIAC_CONTEXT) {
314 gen t,a,b,vmin,vmax;
315 vecteur tmpvars;
316 for (const_iterateur it=vars.begin();it!=vars.end();++it) {
317 vecteur as;
318 vmin=vmax=undef;
319 for (const_iterateur jt=ineq.begin();jt!=ineq.end();++jt) {
320 if (is_linear_wrt(*jt,*it,a,b,contextptr) &&
321 is_constant_wrt_vars(a,vars,contextptr) &&
322 is_constant_wrt_vars(b,vars,contextptr))
323 as.push_back(symb_inferieur_egal(*jt,0));
324 }
325 if (!as.empty()) {
326 gen s=_solve(makesequence(as,*it),contextptr);
327 if (s.type==_VECT)
328 as=*s._VECTptr;
329 else {
330 *logptr(contextptr) << "Warning: failed to set bounds for variable " << *it << "\n";
331 as.clear();
332 }
333 }
334 if (as.size()==1) {
335 const gen &s = as.front();
336 if (s.is_symb_of_sommet(at_inferieur_egal) &&
337 s._SYMBptr->feuille._VECTptr->front()==*it)
338 vmax=s._SYMBptr->feuille._VECTptr->back();
339 else if (s.is_symb_of_sommet(at_superieur_egal) &&
340 s._SYMBptr->feuille._VECTptr->front()==*it)
341 vmin=s._SYMBptr->feuille._VECTptr->back();
342 else if (s.is_symb_of_sommet(at_and) &&
343 s._SYMBptr->feuille._VECTptr->size()==2 &&
344 s._SYMBptr->feuille._VECTptr->front().is_symb_of_sommet(at_superieur_egal) &&
345 s._SYMBptr->feuille._VECTptr->front()._SYMBptr->feuille._VECTptr->front()==*it &&
346 s._SYMBptr->feuille._VECTptr->back().is_symb_of_sommet(at_inferieur_egal) &&
347 s._SYMBptr->feuille._VECTptr->back()._SYMBptr->feuille._VECTptr->front()==*it) {
348 vmin=s._SYMBptr->feuille._VECTptr->front()._SYMBptr->feuille._VECTptr->back();
349 vmax=s._SYMBptr->feuille._VECTptr->back()._SYMBptr->feuille._VECTptr->back();
350 } else *logptr(contextptr) << "Warning: failed to set bounds for variable " << *it << "\n";
351 }
352 gen v=identificateur(" "+it->print(contextptr));
353 if (!is_undef(vmax) && !is_undef(vmin))
354 assume_t_in_ab(v,vmin,vmax,open,open,contextptr);
355 else if (!is_undef(vmin))
356 giac_assume(open?symb_superieur_strict(v,vmin):symb_superieur_egal(v,vmin),contextptr);
357 else if (!is_undef(vmax))
358 giac_assume(open?symb_inferieur_strict(v,vmax):symb_inferieur_egal(v,vmax),contextptr);
359 tmpvars.push_back(v);
360 }
361 return tmpvars;
362 }
363
364 /*
365 * Determine critical points of function f under constraints g<=0 and h=0 using
366 * Karush-Kuhn-Tucker conditions.
367 */
solve_kkt(const gen & f,const vecteur & g,const vecteur & h,const vecteur & vars_orig,GIAC_CONTEXT)368 vecteur solve_kkt(const gen &f,const vecteur &g,const vecteur &h,const vecteur &vars_orig,GIAC_CONTEXT) {
369 int n=vars_orig.size(),m=g.size(),l=h.size();
370 vecteur vars(vars_orig),mug;
371 matrice gr_g,gr_h;
372 vars.resize(n+m+l);
373 gen gr_f_tmp=_grad(makesequence(f,vars_orig),contextptr);
374 if (gr_f_tmp.type!=_VECT || gr_f_tmp._VECTptr->size()!=vars_orig.size()) {
375 *logptr(contextptr) << "Error: failed to compute gradient of " << f << "\n";
376 return vecteur(0);
377 }
378 vecteur &gr_f=*gr_f_tmp._VECTptr;
379 for (int i=0;i<m;++i) {
380 vars[n+i]=identificateur(" mu"+print_INT_(++var_index));
381 giac_assume(symb_superieur_strict(vars[n+i],gen(0)),contextptr); // dual feasibility
382 gr_g.push_back(*_grad(makesequence(g[i],vars_orig),contextptr)._VECTptr);
383 }
384 for (int i=0;i<l;++i) {
385 vars[n+m+i]=identificateur(" lambda"+print_INT_(++var_index));
386 gr_h.push_back(*_grad(makesequence(h[i],vars_orig),contextptr)._VECTptr);
387 }
388 vecteur eqv;
389 for (int i=0;i<n;++i) {
390 gen eq(gr_f[i]);
391 for (int j=0;j<m;++j) {
392 eq+=vars[n+j]*gr_g[j][i];
393 }
394 for (int j=0;j<l;++j) {
395 eq+=vars[n+m+j]*gr_h[j][i];
396 }
397 eqv.push_back(eq);
398 }
399 eqv=mergevecteur(eqv,h); // primal feasibility
400 vector<bool> is_mu_zero(m,false);
401 matrice cv;
402 do {
403 vecteur e(eqv);
404 vecteur v(vars);
405 for (int i=m-1;i>=0;--i) {
406 if (is_mu_zero[i]) {
407 e=subst(e,v[n+i],gen(0),false,contextptr);
408 v.erase(v.begin()+n+i);
409 }
410 else
411 e.push_back(g[i]); // complementary slackness
412 }
413 gen res=solve2(e,v,contextptr);
414 if (res.type==_VECT)
415 cv=mergevecteur(cv,*res._VECTptr);
416 } while(next_binary_perm(is_mu_zero));
417 for (const_iterateur it=vars.begin()+n;it!=vars.end();++it) {
418 _purge(*it,contextptr);
419 }
420 vars.resize(n);
421 for (int i=cv.size();i-->0;) {
422 cv[i]._VECTptr->resize(n);
423 for (int j=0;j<m;++j) {
424 // check primal feasibility
425 if (is_strictly_positive(simp(subst(g[j],vars,cv[i],false,contextptr),contextptr),contextptr)) {
426 cv.erase(cv.begin()+i);
427 break;
428 }
429 }
430 }
431 cpt_simp(cv,vars,f,contextptr);
432 return cv;
433 }
434
435 /*
436 * Determine critical points of an univariate function f(x). Points where it is
437 * not differentiable are considered critical as well as zeros of the first
438 * derivative. Also, bounds of the range of x are critical points.
439 */
critical_univariate(const gen & f,const gen & x,GIAC_CONTEXT)440 matrice critical_univariate(const gen &f,const gen &x,GIAC_CONTEXT) {
441 gen df(derive(f,x,contextptr));
442 matrice cv;
443 gen z=_zeros(makesequence(df,x),contextptr);
444 if (z.type==_VECT)
445 cv=*z._VECTptr;
446 else *logptr(contextptr) << "Error: failed to compute zeros of " << df << "\n";
447 gen den(_denom(df,contextptr));
448 if (!is_constant_wrt(den,x,contextptr)) {
449 z=_zeros(makesequence(den,x),contextptr);
450 if (z.type==_VECT)
451 cv=mergevecteur(cv,*z._VECTptr);
452 }
453 collect_transition_points(f,x,cv,contextptr);
454 for (int i=cv.size();i-->0;) {
455 if (cv[i].is_symb_of_sommet(at_and))
456 cv.erase(cv.begin()+i);
457 else
458 cv[i]=vecteur(1,cv[i]);
459 }
460 return *_epsilon2zero(cv,contextptr)._VECTptr;
461 }
462
463 /*
464 * Compute global minimum mn and global maximum mx of function f(vars) under
465 * conditions g<=0 and h=0. The list of points where global minimum is achieved
466 * is returned.
467 */
global_extrema(const gen & f,const vecteur & g,const vecteur & h,const vecteur & vars,gen & mn,gen & mx,GIAC_CONTEXT)468 vecteur global_extrema(const gen &f,const vecteur &g,const vecteur &h,const vecteur &vars,gen &mn,gen &mx,GIAC_CONTEXT) {
469 int n=vars.size();
470 matrice cv;
471 vecteur tmpvars=make_temp_vars(vars,g,false,contextptr);
472 gen ff=subst(f,vars,tmpvars,false,contextptr);
473 if (n==1) {
474 cv=critical_univariate(ff,tmpvars[0],contextptr);
475 for (const_iterateur it=g.begin();it!=g.end();++it) {
476 gen a,b;
477 if (!is_linear_wrt(*it,vars[0],a,b,contextptr) || is_zero(a)) {
478 *logptr(contextptr) << "Warning: expected a linear function in " << vars[0]
479 << ", got " << *it << "\n";
480 continue;
481 }
482 cv.push_back(makevecteur(-b/a));
483 }
484 } else {
485 vecteur gg=subst(g,vars,tmpvars,false,contextptr);
486 vecteur hh=subst(h,vars,tmpvars,false,contextptr);
487 cv=solve_kkt(ff,gg,hh,tmpvars,contextptr);
488 }
489 for (const_iterateur it=tmpvars.begin();it!=tmpvars.end();++it) {
490 if (find(vars.begin(),vars.end(),*it)==vars.end())
491 _purge(*it,contextptr);
492 }
493 if (cv.empty())
494 return vecteur(0);
495 bool min_set=false,max_set=false;
496 matrice min_locations;
497 for (const_iterateur it=cv.begin();it!=cv.end();++it) {
498 gen val=_eval(subst(f,vars,*it,false,contextptr),contextptr);
499 if (is_inf(val)) {
500 *logptr(contextptr) << "Warning: the function is not bounded\n";
501 return vecteur(0);
502 }
503 if (is_undef(val)) {
504 *logptr(contextptr) << "Warning: failed to compute the function value at critical point "
505 << (it->_VECTptr->size()==1?it->_VECTptr->front():*it) << "\n";
506 return vecteur(0);
507 }
508 if (min_set && is_exactly_zero(simp(val-mn,contextptr))) {
509 if (find(min_locations.begin(),min_locations.end(),*it)==min_locations.end())
510 min_locations.push_back(*it);
511 }
512 else if (!min_set || is_strictly_greater(mn,val,contextptr)) {
513 mn=val;
514 min_set=true;
515 min_locations=vecteur(1,*it);
516 }
517 if (!max_set || is_strictly_greater(val,mx,contextptr)) {
518 mx=val;
519 max_set=true;
520 }
521 }
522 if (n==1) {
523 for (int i=0;i<int(min_locations.size());++i) {
524 min_locations[i]=min_locations[i][0];
525 }
526 }
527 return min_locations;
528 }
529
parse_varlist(const gen & g,vecteur & vars,vecteur & ineq,vecteur & initial,GIAC_CONTEXT)530 int parse_varlist(const gen &g,vecteur &vars,vecteur &ineq,vecteur &initial,GIAC_CONTEXT) {
531 vecteur varlist(g.type==_VECT ? *g._VECTptr : vecteur(1,g));
532 int n=0;
533 for (const_iterateur it=varlist.begin();it!=varlist.end();++it) {
534 if (it->is_symb_of_sommet(at_equal)) {
535 vecteur &ops=*it->_SYMBptr->feuille._VECTptr;
536 gen &v=ops.front(), &rh=ops.back();
537 if (v.type!=_IDNT)
538 return 0;
539 vars.push_back(v);
540 if (rh.is_symb_of_sommet(at_interval)) {
541 vecteur &range=*rh._SYMBptr->feuille._VECTptr;
542 if (!is_inf(range.front()))
543 ineq.push_back(range.front()-v);
544 if (!is_inf(range.back()))
545 ineq.push_back(v-range.back());
546 }
547 else
548 initial.push_back(rh);
549 }
550 else if (it->type!=_IDNT)
551 return 0;
552 else
553 vars.push_back(*it);
554 ++n;
555 }
556 return n;
557 }
558
559 /*
560 * Function 'minimize' minimizes a multivariate continuous real function on a
561 * closed and bounded region using the method of Lagrange multipliers. The
562 * feasible region is specified by bounding variables or by adding one or more
563 * (in)equality constraints.
564 *
565 * Usage
566 * ^^^^^
567 * minimize(obj,[constr],vars,[opt])
568 *
569 * Parameters
570 * ^^^^^^^^^^
571 * - obj : objective function to minimize
572 * - constr (optional) : single equality or inequality constraint or
573 * a list of constraints, if constraint is given as
574 * expression it is assumed that it is equal to zero
575 * - vars : single variable or a list of problem variables, where
576 * optional bounds of a variable may be set by appending '=a..b'
577 * - location (optional) : the option keyword 'locus' or 'coordinates' or 'point'
578 *
579 * Objective function must be continuous in all points of the feasible region,
580 * which is assumed to be closed and bounded. If one of these condinitions is
581 * not met, the final result may be incorrect.
582 *
583 * When the fourth argument is specified, point(s) at which the objective
584 * function attains its minimum value are also returned as a list of vector(s).
585 * The keywords 'locus', 'coordinates' and 'point' all have the same effect.
586 * For univariate problems, a vector of numbers (x values) is returned, while
587 * for multivariate problems it is a vector of vectors, i.e. a matrix.
588 *
589 * The function attempts to obtain the critical points in exact form, if the
590 * parameters of the problem are all exact. It works best for problems in which
591 * the lagrangian function gradient and the constraints are rational expressions.
592 * Points at which the function is not differentiable are also considered critical.
593 * This function handles univariate piecewise functions.
594 *
595 * If no critical points were obtained, the return value is undefined.
596 *
597 * Examples
598 * ^^^^^^^^
599 * minimize(sin(x),[x=0..4])
600 * >> sin(4)
601 * minimize(asin(x),x=-1..1)
602 * >> -pi/2
603 * minimize(x^2+cos(x),x=0..3)
604 * >> 1
605 * minimize(x^4-x^2,x=-3..3,locus)
606 * >> -1/4,[-sqrt(2)/2]
607 * minimize(abs(x),x=-1..1)
608 * >> 0
609 * minimize(x-abs(x),x=-1..1)
610 * >> -2
611 * minimize(abs(exp(-x^2)-1/2),x=-4..4)
612 * >> 0
613 * minimize(piecewise(x<=-2,x+6,x<=1,x^2,3/2-x/2),x=-3..2)
614 * >> 0
615 * minimize(x^2-3x+y^2+3y+3,[x=2..4,y=-4..-2],point)
616 * >> -1,[[2,-2]]
617 * minimize(2x^2+y^2,x+y=1,[x,y])
618 * >> 2/3
619 * minimize(2x^2-y^2+6y,x^2+y^2<=16,[x,y])
620 * >> -40
621 * minimize(x*y+9-x^2-y^2,x^2+y^2<=9,[x,y])
622 * >> -9/2
623 * minimize(sqrt(x^2+y^2)-z,[x^2+y^2<=16,x+y+z=10],[x,y,z])
624 * >> -4*sqrt(2)-6
625 * minimize(x*y*z,x^2+y^2+z^2=1,[x,y,z])
626 * >> -sqrt(3)/9
627 * minimize(sin(x)+cos(x),x=0..20,coordinates)
628 * >> -sqrt(2),[5*pi/4,13*pi/4,21*pi/4]
629 * minimize((1+x^2+3y+5x-4*x*y)/(1+x^2+y^2),x^2/4+y^2/3=9,[x,y])
630 * >> -2.44662490691
631 * minimize(x^2-2x+y^2+1,[x+y<=0,x^2<=4],[x,y],locus)
632 * >> 1/2,[[1/2,-1/2]]
633 * minimize(x^2*(y+1)-2y,[y<=2,sqrt(1+x^2)<=y],[x,y])
634 * >> -4
635 * minimize(4x^2+y^2-2x-4y+1,4x^2+y^2=1,[x,y])
636 * >> -sqrt(17)+2
637 * minimize(cos(x)^2+cos(y)^2,x+y=pi/4,[x,y],locus)
638 * >> (-sqrt(2)+2)/2,[[5*pi/8,-3*pi/8]]
639 * minimize(x^2+y^2,x^4+y^4=2,[x,y])
640 * >> 1.41421356237
641 * minimize(z*x*exp(y),z^2+x^2+exp(2y)=1,[x,y,z])
642 * >> -sqrt(3)/9
643 */
_minimize(const gen & args,GIAC_CONTEXT)644 gen _minimize(const gen &args,GIAC_CONTEXT) {
645 if (args.type==_STRNG && args.subtype==-1) return args;
646 if (args.type!=_VECT || args.subtype!=_SEQ__VECT || args._VECTptr->size()>4)
647 return gentypeerr(contextptr);
648 vecteur &argv=*args._VECTptr,g,h;
649 bool location=false;
650 int nargs=argv.size();
651 if (argv.back()==at_coordonnees || argv.back()==at_lieu || argv.back()==at_point) {
652 location=true;
653 --nargs;
654 }
655 if (nargs==3) {
656 vecteur constr(argv[1].type==_VECT ? *argv[1]._VECTptr : vecteur(1,argv[1]));
657 for (const_iterateur it=constr.begin();it!=constr.end();++it) {
658 if (it->is_symb_of_sommet(at_equal))
659 h.push_back(equal2diff(*it));
660 else if (it->is_symb_of_sommet(at_superieur_egal) ||
661 it->is_symb_of_sommet(at_inferieur_egal)) {
662 vecteur &s=*it->_SYMBptr->feuille._VECTptr;
663 g.push_back(it->is_symb_of_sommet(at_inferieur_egal)?s[0]-s[1]:s[1]-s[0]);
664 } else if (it->type==_IDNT || it->type==_SYMB)
665 h.push_back(*it);
666 else *logptr(contextptr) << "Warning: ignoring constraint " << *it << "\n";
667 }
668 }
669 vecteur vars,initial;
670 int n; // number of variables
671 if ((n=parse_varlist(argv[nargs-1],vars,g,initial,contextptr))==0 || !initial.empty())
672 return gensizeerr(contextptr);
673 gen &f=argv[0];
674 gen mn,mx;
675 vecteur loc;
676 bool has_symb=false;
677 gen simb=_lname(makevecteur(f,g,h),contextptr);
678 if (simb.type==_VECT) {
679 for (const_iterateur it=simb._VECTptr->begin();it!=simb._VECTptr->end();++it) {
680 if (find(vars.begin(),vars.end(),*it)==vars.end()) {
681 has_symb=true;
682 break;
683 }
684 }
685 }
686 try {
687 loc=global_extrema(f,g,h,vars,mn,mx,contextptr);
688 } catch (const std::exception &e) {
689 loc.clear();
690 }
691 if (loc.empty()) {
692 if (has_symb)
693 return undef;
694 //*logptr(contextptr) << "Warning: switching to approx mode\n";
695 gen asol=_nlpsolve(makesequence(f,mergevecteur(
696 *_zip(makesequence(at_inferieur_egal,g,vecteur(g.size(),0)),contextptr)._VECTptr,
697 *_zip(makesequence(at_equal,h,vecteur(h.size(),0)),contextptr)._VECTptr),vars),contextptr);
698 if (asol.type==_VECT && asol._VECTptr->size()>1) {
699 mn=asol._VECTptr->front();
700 if (location) {
701 loc.resize(n);
702 gen pos=asol._VECTptr->at(1),v;
703 if (pos.type!=_VECT || pos._VECTptr->size()!=n) return undef;
704 for (const_iterateur it=pos._VECTptr->begin();it!=pos._VECTptr->end();++it) {
705 if (!it->is_symb_of_sommet(at_equal)) return undef;
706 const_iterateur jt=find(vars.begin(),vars.end(),it->_SYMBptr->feuille._VECTptr->front());
707 if (jt==vars.end()) return undef;
708 loc[jt-vars.begin()]=it->_SYMBptr->feuille._VECTptr->back();
709 }
710 }
711 } else return undef;
712 } else mn=simp(mn,contextptr);
713 if (location)
714 return makevecteur(mn,loc);
715 return mn;
716 }
717 static const char _minimize_s []="minimize";
718 static define_unary_function_eval (__minimize,&_minimize,_minimize_s);
719 define_unary_function_ptr5(at_minimize,alias_at_minimize,&__minimize,0,true)
720
721 /*
722 * 'maximize' takes the same arguments as the function 'minimize', but
723 * maximizes the objective function. See 'minimize' for details.
724 *
725 * Examples
726 * ^^^^^^^^
727 * maximize(cos(x),x=1..3)
728 * >> cos(1)
729 * maximize(piecewise(x<=-2,x+6,x<=1,x^2,3/2-x/2),x=-3..2)
730 * >> 4
731 * minimize(x-abs(x),x=-1..1)
732 * >> 0
733 * maximize(x^2-3x+y^2+3y+3,[x=2..4,y=-4..-2])
734 * >> 11
735 * maximize(x*y*z,x^2+2*y^2+3*z^2<=1,[x,y,z],point)
736 * >> sqrt(2)/18,[[-sqrt(3)/3,sqrt(6)/6,-1/3],[sqrt(3)/3,-sqrt(6)/6,-1/3],
737 * [-sqrt(3)/3,-sqrt(6)/6,1/3],[sqrt(3)/3,sqrt(6)/6,1/3]]
738 * maximize(x^2-x*y+2*y^2,[x=-1..0,y=-1/2..1/2],coordinates)
739 * >> 2,[[-1,1/2]]
740 * maximize(x*y,[x+y^2<=2,x>=0,y>=0],[x,y],locus)
741 * >> 4*sqrt(6)/9,[[4/3,sqrt(6)/3]]
742 * maximize(y^2-x^2*y,y<=x,[x=0..2,y=0..2])
743 * >> 4/27
744 * maximize(2x+y,4x^2+y^2=8,[x,y])
745 * >> 4
746 * maximize(x^2*(y+1)-2y,[y<=2,sqrt(1+x^2)<=y],[x,y])
747 * >> 5
748 * maximize(4x^2+y^2-2x-4y+1,4x^2+y^2=1,[x,y])
749 * >> sqrt(17)+2
750 * maximize(3x+2y,2x^2+3y^2<=3,[x,y])
751 * >> sqrt(70)/2
752 * maximize(x*y,[2x+3y<=10,x>=0,y>=0],[x,y])
753 * >> 25/6
754 * maximize(x^2+y^2+z^2,[x^2/16+y^2+z^2=1,x+y+z=0],[x,y,z])
755 * >> 8/3
756 * assume(a>0);maximize(x^2*y^2*z^2,x^2+y^2+z^2=a^2,[x,y,z])
757 * >> a^6/27
758 */
_maximize(const gen & g,GIAC_CONTEXT)759 gen _maximize(const gen &g,GIAC_CONTEXT) {
760 if (g.type==_STRNG && g.subtype==-1) return g;
761 if (g.type!=_VECT || g.subtype!=_SEQ__VECT || g._VECTptr->size()<2)
762 return gentypeerr(contextptr);
763 vecteur gv(*g._VECTptr);
764 gv[0]=-gv[0];
765 gen res=_minimize(_feuille(gv,contextptr),contextptr);
766 if (res.type==_VECT && res._VECTptr->size()>0) {
767 res._VECTptr->front()=-res._VECTptr->front();
768 }
769 else if (res.type!=_VECT)
770 res=-res;
771 return ratnormal(res,contextptr);
772 }
773 static const char _maximize_s []="maximize";
774 static define_unary_function_eval (__maximize,&_maximize,_maximize_s);
775 define_unary_function_ptr5(at_maximize,alias_at_maximize,&__maximize,0,true)
776
sum_ivector(const ivector & v,bool drop_last)777 int ipdiff::sum_ivector(const ivector &v,bool drop_last) {
778 int res=0;
779 for (ivector_iter it=v.begin();it!=v.end()-drop_last?1:0;++it) {
780 res+=*it;
781 }
782 return res;
783 }
784
785 /*
786 * IPDIFF CLASS IMPLEMENTATION
787 */
788
ipdiff(const gen & f_orig,const vecteur & g_orig,const vecteur & vars_orig,GIAC_CONTEXT)789 ipdiff::ipdiff(const gen &f_orig,const vecteur &g_orig,const vecteur &vars_orig,GIAC_CONTEXT) {
790 ctx=contextptr;
791 f=f_orig;
792 g=g_orig;
793 vars=vars_orig;
794 ord=0;
795 nconstr=g.size();
796 nvars=vars.size()-nconstr;
797 assert(nvars>0);
798 pdv[ivector(nvars,0)]=f; // make the zeroth order derivative initially available
799 }
800
ipartition(int m,int n,ivectors & c,const ivector & p)801 void ipdiff::ipartition(int m,int n,ivectors &c,const ivector &p) {
802 for (int i=0;i<n;++i) {
803 if (!p.empty() && p[i]!=0)
804 continue;
805 ivector r;
806 if (p.empty())
807 r.resize(n,0);
808 else r=p;
809 for (int j=0;j<m;++j) {
810 ++r[i];
811 int s=sum_ivector(r);
812 if (s==m && find(c.begin(),c.end(),r)==c.end())
813 c.push_back(r);
814 else if (s<m)
815 ipartition(m,n,c,r);
816 else break;
817 }
818 }
819 }
820
derive_diffterms(const diffterms & terms,ivector & sig)821 ipdiff::diffterms ipdiff::derive_diffterms(const diffterms &terms,ivector &sig) {
822 while (!sig.empty() && sig.back()==0) {
823 sig.pop_back();
824 }
825 if (sig.empty())
826 return terms;
827 int k=sig.size()-1,p;
828 diffterms tv;
829 ivector u(nvars+1,0);
830 for (diffterms::const_iterator it=terms.begin();it!=terms.end();++it) {
831 int c=it->second;
832 diffterm t(it->first);
833 const ivector_map &h_orig=it->first.second;
834 ++t.first.at(k);
835 tv[t]+=c;
836 --t.first.at(k);
837 ivector_map h(h_orig);
838 for (ivector_map::const_iterator jt=h_orig.begin();jt!=h_orig.end();++jt) {
839 ivector v=jt->first;
840 if ((p=jt->second)==0)
841 continue;
842 if (p==1)
843 h.erase(h.find(v));
844 else
845 --h[v];
846 ++v[k];
847 ++h[v];
848 t.second=h;
849 tv[t]+=c*p;
850 --h[v];
851 --v[k];
852 ++h[v];
853 }
854 t.second=h_orig;
855 for (int i=0;i<nconstr;++i) {
856 ++t.first.at(nvars+i);
857 u[k]=1;
858 u.back()=i;
859 ++t.second[u];
860 tv[t]+=c;
861 --t.first.at(nvars+i);
862 --t.second[u];
863 u[k]=0;
864 }
865 }
866 --sig.back();
867 return derive_diffterms(tv,sig);
868 }
869
get_pd(const pd_map & pds,const ivector & sig) const870 const gen &ipdiff::get_pd(const pd_map &pds,const ivector &sig) const {
871 try {
872 return pds.at(sig);
873 }
874 catch (out_of_range &e) {
875 return undef;
876 }
877 }
878
differentiate(const gen & e,pd_map & pds,const ivector & sig)879 const gen &ipdiff::differentiate(const gen &e,pd_map &pds,const ivector &sig) {
880 const gen &pd=get_pd(pds,sig);
881 if (!is_undef(pd))
882 return pd;
883 vecteur v(1,e);
884 bool do_derive=false;
885 assert(vars.size()<=sig.size());
886 for (int i=0;i<int(vars.size());++i) {
887 if (sig[i]>0) {
888 v=mergevecteur(v,vecteur(sig[i],vars[i]));
889 do_derive=true;
890 }
891 }
892 if (do_derive)
893 return pds[sig]=_derive(_feuille(v,ctx),ctx);
894 return e;
895 }
896
compute_h(const vector<diffterms> & grv,int order)897 void ipdiff::compute_h(const vector<diffterms> &grv,int order) {
898 if (g.empty())
899 return;
900 ivectors hsigv;
901 matrice A;
902 vecteur b(g.size()*grv.size(),gen(0));
903 gen t;
904 int grv_sz=grv.size();
905 for (int i=0;i<nconstr;++i) {
906 for (int j=0;j<grv_sz;++j) {
907 vecteur eq(g.size()*grv_sz,gen(0));
908 const diffterms &grvj=grv[j];
909 for (diffterms::const_iterator it=grvj.begin();it!=grvj.end();++it) {
910 ivector sig(it->first.first),hsig;
911 sig.push_back(i);
912 t=gen(it->second)*differentiate(g[i],pdg,sig);
913 for (ivector_map::const_iterator ht=it->first.second.begin();ht!=it->first.second.end();++ht) {
914 if (ht->second==0)
915 continue;
916 const ivector &sigh=ht->first;
917 if (sum_ivector(sigh,true)<order) {
918 gen h(get_pd(pdh,sigh));
919 assert(!is_undef(h));
920 t=t*pow(h,ht->second);
921 }
922 else {
923 assert(ht->second==1);
924 hsig=sigh;
925 }
926 }
927 if (hsig.empty())
928 b[grv_sz*i+j]-=t;
929 else {
930 int k=0,hsigv_sz=hsigv.size();
931 for (;k<hsigv_sz;++k) {
932 if (hsigv[k]==hsig)
933 break;
934 }
935 eq[k]+=t;
936 if (k==hsigv_sz)
937 hsigv.push_back(hsig);
938 }
939 }
940 A.push_back(*simp(eq,ctx)._VECTptr);
941 }
942 }
943 matrice B;
944 B.push_back(*simp(b,ctx)._VECTptr);
945 matrice invA=*_inv(A,ctx)._VECTptr;
946 vecteur sol(*mtran(mmult(invA,mtran(B))).front()._VECTptr);
947 for (int i=0;i<int(sol.size());++i) {
948 pdh[hsigv[i]]=simp(sol[i],ctx);
949 }
950 }
951
find_nearest_terms(const ivector & sig,diffterms & match,ivector & excess)952 void ipdiff::find_nearest_terms(const ivector &sig,diffterms &match,ivector &excess) {
953 excess=sig;
954 int i;
955 for (map<ivector,diffterms>::const_iterator it=cterms.begin();it!=cterms.end();++it) {
956 ivector ex(nvars,0);
957 for (i=0;i<nvars;++i) {
958 if ((ex[i]=sig[i]-it->first.at(i))<0)
959 break;
960 }
961 if (i<nvars)
962 continue;
963 if (sum_ivector(ex)<sum_ivector(excess)) {
964 excess=ex;
965 match=it->second;
966 }
967 }
968 }
969
raise_order(int order)970 void ipdiff::raise_order(int order) {
971 if (g.empty())
972 return;
973 ivectors c;
974 ivector excess,init_f(nvars+nconstr,0);
975 diffterm init_term;
976 init_term.first=init_f;
977 diffterms init_terms;
978 init_terms[init_term]=1;
979 vector<diffterms> grv;
980 for (int k=ord+1;k<=order;++k) {
981 grv.clear();
982 c.clear();
983 ipartition(k,nvars,c);
984 for (ivectors::const_iterator it=c.begin();it!=c.end();++it) {
985 diffterms terms=init_terms;
986 find_nearest_terms(*it,terms,excess);
987 if (sum_ivector(excess)>0) {
988 terms=derive_diffterms(terms,excess);
989 cterms[*it]=terms;
990 }
991 grv.push_back(terms);
992 }
993 compute_h(grv,k);
994 }
995 ord=order;
996 }
997
compute_pd(int order,const ivector & sig)998 void ipdiff::compute_pd(int order,const ivector &sig) {
999 gen pd;
1000 ivectors c;
1001 ipartition(order,nvars,c);
1002 for (ivectors::const_iterator ct=c.begin();ct!=c.end();++ct) {
1003 if (!sig.empty() && sig!=*ct)
1004 continue;
1005 if (g.empty()) {
1006 differentiate(f,pdv,sig);
1007 continue;
1008 }
1009 diffterms &terms=cterms[*ct];
1010 pd=gen(0);
1011 for (diffterms::const_iterator it=terms.begin();it!=terms.end();++it) {
1012 ivector sig(it->first.first);
1013 gen t(gen(it->second)*differentiate(f,pdf,sig));
1014 if (!is_zero(t)) {
1015 for (ivector_map::const_iterator jt=it->first.second.begin();jt!=it->first.second.end();++jt) {
1016 if (jt->second==0)
1017 continue;
1018 gen h(get_pd(pdh,jt->first));
1019 assert(!is_undef(h));
1020 t=t*pow(h,jt->second);
1021 }
1022 pd+=t;
1023 }
1024 }
1025 pdv[*ct]=simp(pd,ctx);
1026 }
1027 }
1028
gradient(vecteur & res)1029 void ipdiff::gradient(vecteur &res) {
1030 if (nconstr==0)
1031 res=*_grad(makesequence(f,vars),ctx)._VECTptr;
1032 else {
1033 res.resize(nvars);
1034 ivector sig(nvars,0);
1035 if (ord<1) {
1036 raise_order(1);
1037 compute_pd(1);
1038 }
1039 for (int i=0;i<nvars;++i) {
1040 sig[i]=1;
1041 res[i]=derivative(sig);
1042 sig[i]=0;
1043 }
1044 }
1045 }
1046
hessian(matrice & res)1047 void ipdiff::hessian(matrice &res) {
1048 if (nconstr==0)
1049 res=*_hessian(makesequence(f,vars),ctx)._VECTptr;
1050 else {
1051 res.clear();
1052 ivector sig(nvars,0);
1053 if (ord<2) {
1054 raise_order(2);
1055 compute_pd(2);
1056 }
1057 for (int i=0;i<nvars;++i) {
1058 vecteur r(nvars);
1059 ++sig[i];
1060 for (int j=0;j<nvars;++j) {
1061 ++sig[j];
1062 r[j]=derivative(sig);
1063 --sig[j];
1064 }
1065 res.push_back(r);
1066 --sig[i];
1067 }
1068 }
1069 }
1070
derivative(const ivector & sig)1071 const gen &ipdiff::derivative(const ivector &sig) {
1072 if (nconstr==0)
1073 return differentiate(f,pdf,sig);
1074 int k=sum_ivector(sig); // the order of the derivative
1075 if (k>ord) {
1076 raise_order(k);
1077 compute_pd(k,sig);
1078 }
1079 return get_pd(pdv,sig);
1080 }
1081
derivative(const vecteur & dvars)1082 const gen &ipdiff::derivative(const vecteur &dvars) {
1083 ivector sig(nvars,0);
1084 const_iterateur jt;
1085 for (const_iterateur it=dvars.begin();it!=dvars.end();++it) {
1086 if ((jt=find(vars.begin(),vars.end(),*it))==vars.end())
1087 return undef;
1088 ++sig[jt-vars.begin()];
1089 }
1090 return derivative(sig);
1091 }
1092
partial_derivatives(int order,pd_map & pdmap)1093 void ipdiff::partial_derivatives(int order,pd_map &pdmap) {
1094 if (nconstr>0 && ord<order) {
1095 raise_order(order);
1096 compute_pd(order);
1097 }
1098 ivectors c;
1099 ipartition(order,nvars,c);
1100 for (ivectors::const_iterator it=c.begin();it!=c.end();++it) {
1101 pdmap[*it]=derivative(*it);
1102 }
1103 }
1104
taylor_term(const vecteur & a,int k,bool shift)1105 gen ipdiff::taylor_term(const vecteur &a,int k,bool shift) {
1106 assert(k>=0);
1107 if (k==0)
1108 return subst(f,vars,a,false,ctx);
1109 ivectors sigv;
1110 ipartition(k,nvars,sigv);
1111 gen term(0);
1112 if (nconstr>0) while (k>ord) {
1113 raise_order(ord+1);
1114 compute_pd(ord);
1115 }
1116 for (ivectors::const_iterator it=sigv.begin();it!=sigv.end();++it) {
1117 gen pd;
1118 if (g.empty()) {
1119 vecteur args(1,f);
1120 for (int i=0;i<nvars;++i) {
1121 for (int j=0;j<it->at(i);++j) {
1122 args.push_back(vars[i]);
1123 }
1124 }
1125 pd=_derive(_feuille(args,ctx),ctx);
1126 }
1127 else
1128 pd=derivative(*it);
1129 pd=subst(pd,vars,a,false,ctx);
1130 for (int i=0;i<nvars;++i) {
1131 int ki=it->at(i);
1132 if (ki==0)
1133 continue;
1134 pd=pd*(shift?pow(vars[i]-a[i],ki):pow(vars[i],ki))/factorial(ki);
1135 }
1136 term+=pd;
1137 }
1138 return term;
1139 }
1140
taylor(const vecteur & a,int order)1141 gen ipdiff::taylor(const vecteur &a,int order) {
1142 assert(order>=0);
1143 gen T(0);
1144 for (int k=0;k<=order;++k) {
1145 T+=taylor_term(a,k);
1146 }
1147 return T;
1148 }
1149
1150 /*
1151 * END OF IPDIFF CLASS
1152 */
1153
vars_arrangements(matrice J,ipdiff::ivectors & arrs,GIAC_CONTEXT)1154 void vars_arrangements(matrice J,ipdiff::ivectors &arrs,GIAC_CONTEXT) {
1155 int m=J.size(),n=J.front()._VECTptr->size();
1156 assert(n<=32 && m<n);
1157 matrice tJ(mtran(J));
1158 ulong N=std::pow(2,n);
1159 vector<ulong> sets(comb(n,m).val);
1160 int i=0;
1161 for (ulong k=1;k<N;++k) {
1162 bitset<32> b(k);
1163 if (b.count()==(size_t)m)
1164 sets[i++]=k;
1165 }
1166 matrice S;
1167 ipdiff::ivector arr(n);
1168 for (vector<ulong>::const_iterator it=sets.begin();it!=sets.end();++it) {
1169 for (i=0;i<n;++i) arr[i]=i;
1170 N=std::pow(2,n);
1171 for (i=n;i-->0;) {
1172 N/=2;
1173 if ((*it & N)!=0) {
1174 arr.erase(arr.begin()+i);
1175 arr.push_back(i);
1176 }
1177 }
1178 S.clear();
1179 for (ipdiff::ivector::const_iterator it=arr.end()-m;it!=arr.end();++it) {
1180 S.push_back(tJ[*it]);
1181 }
1182 if (!is_zero(_det(S,contextptr)))
1183 arrs.push_back(arr);
1184 }
1185 }
1186
jacobian(vecteur & g,vecteur & vars,GIAC_CONTEXT)1187 matrice jacobian(vecteur &g,vecteur &vars,GIAC_CONTEXT) {
1188 matrice J;
1189 for (int i=0;i<int(g.size());++i) {
1190 gen gr=_grad(makesequence(g[i],vars),contextptr);
1191 if (gr.type==_VECT && gr._VECTptr->size()==vars.size())
1192 J.push_back(*gr._VECTptr);
1193 else {
1194 *logptr(contextptr) << "Error: failed to compute gradient of " << g[i] << "\n";
1195 return vecteur(0);
1196 }
1197 }
1198 return J;
1199 }
1200
ck_jacobian(vecteur & g,vecteur & vars,GIAC_CONTEXT)1201 bool ck_jacobian(vecteur &g,vecteur &vars,GIAC_CONTEXT) {
1202 matrice J(jacobian(g,vars,contextptr));
1203 if (!g.empty() && J.empty())
1204 return false;
1205 int m=g.size();
1206 int n=vars.size()-m;
1207 if (_rank(J,contextptr).val<m)
1208 return false;
1209 J=mtran(J);
1210 J.erase(J.begin(),J.begin()+n);
1211 return !is_zero(_det(J,contextptr));
1212 }
1213
1214 /*
1215 * 'implicitdiff' differentiates function(s) defined by equation(s) or a
1216 * function f(x1,x2,...,xn,y1,y2,...,ym) where y1,...,ym are functions of
1217 * x1,x2,...,xn defined by m equality constraints.
1218 *
1219 * Usage
1220 * ^^^^^
1221 * implicitdiff(f,constr,depvars,diffvars)
1222 * implicitdiff(f,constr,vars,order_size=<posint>,[P])
1223 * implicitdiff(constr,[depvars],y,diffvars)
1224 *
1225 * Parameters
1226 * ^^^^^^^^^^
1227 * - f : expression
1228 * - constr : (list of) equation(s)
1229 * - depvars : (list of) dependent variable(s), each of them given
1230 * either as a symbol, e.g. y, or a function, e.g. y(x,z)
1231 * - diffvars : sequence of variables w.r.t. which the differentiation
1232 * will be carried out
1233 * - vars : list all variables on which f depends such that
1234 * : dependent variables come after independent ones
1235 * - P : (list of) coordinate(s) to compute derivatives at
1236 * - y : (list of) dependent variable(s) to differentiate w.r.t.
1237 * diffvars, each of them given as a symbol
1238 *
1239 * The return value is partial derivative specified by diffvars. If
1240 * 'order=m' is given as the fourth argument, all partial derivatives of
1241 * order m will be computed and returned as vector for m=1, matrix for m=2 or
1242 * table for m>2. The first two cases produce gradient and hessian of f,
1243 * respectively. For m>2, the partial derivative
1244 * pd=d^m(f)/(d^k1(x1)*d^k2(x2)*...*d^kn(xn)) is saved under key [k1,k2,...kn].
1245 * If P is specified, pd(P) is saved.
1246 *
1247 * Examples
1248 * ^^^^^^^^
1249 * implicitdiff(x^2*y+y^2=1,y,x)
1250 * >> -2*x*y/(x^2+2*y)
1251 * implicitdiff(R=P*V/T,P,T)
1252 * >> P/T
1253 * implicitdiff([x^2+y=z,x+y*z=1],[y(x),z(x)],y,x)
1254 * >> (-2*x*y-1)/(y+z)
1255 * implicitdiff([x^2+y=z,x+y*z=1],[y(x),z(x)],[y,z],x)
1256 * >> [(-2*x*y-1)/(y+z),(2*x*z-1)/(y+z)]
1257 * implicitdiff(y=x^2/z,y,x)
1258 * >> 2x/z
1259 * implicitdiff(y=x^2/z,y,z)
1260 * >> -x^2/z^2
1261 * implicitdiff(y^3+x^2=1,y,x)
1262 * >> -2*x/(3*y^2)
1263 * implicitdiff(y^3+x^2=1,y,x,x)
1264 * >> (-8*x^2-6*y^3)/(9*y^5)x+3y-z,2x^2+y^2=z,[x,y,z]
1265 * implicitdiff(a*x^3*y-2y/z=z^2,y(x,z),x)
1266 * >> -3*a*x^2*y*z/(a*x^3*z-2)
1267 * implicitdiff(a*x^3*y-2y/z=z^2,y(x,z),x,z)
1268 * >> (12*a*x^2*y-6*a*x^2*z^3)/(a^2*x^6*z^2-4*a*x^3*z+4)
1269 * implicitdiff([-2x*z+y^2=1,x^2-exp(x*z)=y],[y(x),z(x)],y,x)
1270 * >> 2*x/(y*exp(x*z)+1)
1271 * implicitdiff([-2x*z+y^2=1,x^2-exp(x*z)=y],[y(x),z(x)],[y,z],x)
1272 * >> [2*x/(y*exp(x*z)+1),(2*x*y-y*z*exp(x*z)-z)/(x*y*exp(x*z)+x)]
1273 * implicitdiff([a*sin(u*v)+b*cos(w*x)=c,u+v+w+x=z,u*v+w*x=z],[u(x,z),v(x,z),w(x,z)],u,z)
1274 * >> (a*u*x*cos(u*v)-a*u*cos(u*v)+b*u*x*sin(w*x)-b*x*sin(w*x))/
1275 * (a*u*x*cos(u*v)-a*v*x*cos(u*v)+b*u*x*sin(w*x)-b*v*x*sin(w*x))
1276 * implicitdiff(x*y,-2x^3+15x^2*y+11y^3-24y=0,y(x),x$2)
1277 * >> (162*x^5*y+1320*x^4*y^2-320*x^4-3300*x^3*y^3+800*x^3*y+968*x^2*y^4-1408*x^2*y^2+
1278 * 512*x^2-3630*x*y^5+5280*x*y^3-1920*x*y)/(125*x^6+825*x^4*y^2-600*x^4+
1279 * 1815*x^2*y^4-2640*x^2*y^2+960*x^2+1331*y^6-2904*y^4+2112*y^2-512)
1280 * implicitdiff((x-u)^2+(y-v)^2,[x^2/4+y^2/9=1,(u-3)^2+(v+5)^2=1],[v(u,x),y(u,x)],u,x)
1281 * >> (-9*u*x-4*v*y+27*x-20*y)/(2*v*y+10*y)
1282 * implicitdiff(x*y*z,-2x^3+15x^2*y+11y^3-24y=0,[x,z,y],order_size=1)
1283 * >> [(2*x^3*z-5*x^2*y*z+11*y^3*z-8*y*z)/(5*x^2+11*y^2-8),x*y]
1284 * implicitdiff(x*y*z,-2x^3+15x^2*y+11y^3-24y=0,[x,z,y],order_size=2,[1,-1,0])
1285 * >> [[64/9,-2/3],[-2/3,0]]
1286 * pd:=implicitdiff(x*y*z,-2x^3+15x^2*y+11y^3-24y=0,[x,z,y],order_size=4,[0,z,0]);pd[4,0,0]
1287 * >> -2*z
1288 */
_implicitdiff(const gen & g,GIAC_CONTEXT)1289 gen _implicitdiff(const gen &g,GIAC_CONTEXT) {
1290 if (g.type==_STRNG && g.subtype==-1) return g;
1291 if (g.type!=_VECT || g.subtype!=_SEQ__VECT || g._VECTptr->size()<2)
1292 return gentypeerr(contextptr);
1293 vecteur &gv=*g._VECTptr;
1294 gen &f=gv[0];
1295 if (int(gv.size())<3)
1296 return gensizeerr(contextptr);
1297 int ci=gv[0].type!=_VECT && !gv[0].is_symb_of_sommet(at_equal)?1:0;
1298 vecteur freevars,depvars,diffdepvars;
1299 gen_map diffvars;
1300 // get the constraints as a list of vanishing expressions
1301 vecteur constr(gv[ci].type==_VECT?*gv[ci]._VECTptr:vecteur(1,gv[ci]));
1302 for (int i=0;i<int(constr.size());++i) {
1303 if (constr[i].is_symb_of_sommet(at_equal))
1304 constr[i]=equal2diff(constr[i]);
1305 }
1306 int m=constr.size();
1307 int dvi=3;
1308 if (ci==0) {
1309 if (gv[ci+1].type==_VECT)
1310 diffdepvars=gv[ci+2].type==_VECT?*gv[ci+2]._VECTptr:vecteur(1,gv[ci+2]);
1311 else
1312 dvi=2;
1313 }
1314 bool compute_all=false;
1315 int order=0;
1316 if (ci==1 && gv[dvi].is_symb_of_sommet(at_equal)) {
1317 vecteur &v=*gv[dvi]._SYMBptr->feuille._VECTptr;
1318 if (v.front()!=at_order || !v.back().is_integer())
1319 return gentypeerr(contextptr);
1320 order=v.back().val;
1321 if (order<=0)
1322 return gendimerr(contextptr);
1323 compute_all=true;
1324 }
1325 // get dependency specification
1326 vecteur deplist(gv[ci+1].type==_VECT?*gv[ci+1]._VECTptr:vecteur(1,gv[ci+1]));
1327 if (compute_all) {
1328 // vars must be specified as x1,x2,...,xn,y1,y2,...,ym
1329 int nd=deplist.size();
1330 if (nd<=m)
1331 return gensizeerr(contextptr);
1332 for (int i=0;i<nd;++i) {
1333 if (i<nd-m)
1334 freevars.push_back(deplist[i]);
1335 else
1336 depvars.push_back(deplist[i]);
1337 }
1338 }
1339 else {
1340 // get (in)dependent variables
1341 for (const_iterateur it=deplist.begin();it!=deplist.end();++it) {
1342 if (it->type==_IDNT)
1343 depvars.push_back(*it);
1344 else if (it->is_symb_of_sommet(at_of)) {
1345 vecteur fe(*it->_SYMBptr->feuille._VECTptr);
1346 depvars.push_back(fe.front());
1347 if (fe.back().type==_VECT) {
1348 for (int i=0;i<int(fe.back()._VECTptr->size());++i) {
1349 gen &x=fe.back()._VECTptr->at(i);
1350 if (find(freevars.begin(),freevars.end(),x)==freevars.end())
1351 freevars.push_back(x);
1352 }
1353 }
1354 else
1355 freevars.push_back(fe.back());
1356 }
1357 else
1358 return gentypeerr(contextptr);
1359 }
1360 // get diffvars
1361 for (const_iterateur it=gv.begin()+dvi;it!=gv.end();++it) {
1362 gen v(eval(*it,contextptr));
1363 gen x;
1364 if (v.type==_IDNT)
1365 diffvars[(x=v)]+=1;
1366 else if (v.type==_VECT && v.subtype==_SEQ__VECT)
1367 diffvars[(x=v._VECTptr->front())]+=v._VECTptr->size();
1368 else
1369 return gentypeerr(contextptr);
1370 if (find(freevars.begin(),freevars.end(),x)==freevars.end())
1371 freevars.push_back(x);
1372 }
1373 }
1374 int n=freevars.size(); // number of independent variables
1375 if (m!=int(depvars.size()))
1376 return gensizeerr(contextptr);
1377 vecteur vars(mergevecteur(freevars,depvars)); // list of all variables
1378 // check whether the conditions of implicit function theorem hold
1379 if (!ck_jacobian(constr,vars,contextptr))
1380 return gensizeerr(contextptr);
1381 // build partial derivative specification 'sig'
1382 ipdiff::ivector sig(n,0); // sig[i]=k means: derive k times with respect to ith independent variable
1383 ipdiff ipd(f,constr,vars,contextptr);
1384 if (compute_all) {
1385 vecteur pt(0);
1386 if (int(gv.size())>4) {
1387 pt=gv[4].type==_VECT?*gv[4]._VECTptr:vecteur(1,gv[4]);
1388 if (int(pt.size())!=n+m)
1389 return gensizeerr(contextptr);
1390 }
1391 ipdiff::pd_map pdv;
1392 ipd.partial_derivatives(order,pdv);
1393 if (order==1) {
1394 vecteur gr;
1395 ipd.gradient(gr);
1396 return pt.empty()?gr:simp(subst(gr,vars,pt,false,contextptr),contextptr);
1397 }
1398 else if (order==2) {
1399 matrice hess;
1400 ipd.hessian(hess);
1401 return pt.empty()?hess:simp(subst(hess,vars,pt,false,contextptr),contextptr);
1402 }
1403 else {
1404 ipdiff::ivectors c;
1405 ipdiff::ipartition(order,n,c);
1406 gen_map ret_pdv;
1407 for (ipdiff::ivectors::const_iterator it=c.begin();it!=c.end();++it) {
1408 vecteur v;
1409 for (int i=0;i<n;++i) {
1410 v.push_back(gen(it->at(i)));
1411 }
1412 ret_pdv[v]=pt.empty()?pdv[*it]:simp(subst(pdv[*it],vars,pt,false,contextptr),contextptr);
1413 }
1414 return ret_pdv;
1415 }
1416 }
1417 for (gen_map::const_iterator it=diffvars.begin();it!=diffvars.end();++it) {
1418 int i=0;
1419 for (;i<n;++i) {
1420 if (it->first==freevars[i]) {
1421 sig[i]=it->second.val;
1422 break;
1423 }
1424 }
1425 assert(i<n);
1426 }
1427 // compute the partial derivative specified by 'sig'
1428 order=ipdiff::sum_ivector(sig);
1429 if (ci==1)
1430 return simp(ipd.derivative(sig),contextptr);
1431 vecteur ret;
1432 if (diffdepvars.empty()) {
1433 assert(m==1);
1434 diffdepvars=vecteur(1,depvars.front());
1435 }
1436 for (const_iterateur it=diffdepvars.begin();it!=diffdepvars.end();++it) {
1437 if (find(depvars.begin(),depvars.end(),*it)==depvars.end()) {
1438 // variable *it is not in depvars, so it's treated as a constant
1439 ret.push_back(gen(0));
1440 continue;
1441 }
1442 ipdiff tmp(*it,constr,vars,contextptr);
1443 ret.push_back(simp(tmp.derivative(sig),contextptr));
1444 }
1445 return ret.size()==1?ret.front():ret;
1446 }
1447 static const char _implicitdiff_s []="implicitdiff";
1448 static define_unary_function_eval (__implicitdiff,&_implicitdiff,_implicitdiff_s);
1449 define_unary_function_ptr5(at_implicitdiff,alias_at_implicitdiff,&__implicitdiff,0,true)
1450
find_cpt(vecteur & cpts,const vecteur & cand,GIAC_CONTEXT)1451 iterateur find_cpt(vecteur &cpts,const vecteur &cand,GIAC_CONTEXT) {
1452 iterateur ret;
1453 for (ret=cpts.begin();ret!=cpts.end();++ret) {
1454 if (ret->type!=_VECT || ret->_VECTptr->size()!=2)
1455 continue;
1456 if (ret->_VECTptr->front().type==_VECT) {
1457 const vecteur &c=*ret->_VECTptr->front()._VECTptr;
1458 if (c.size()!=cand.size())
1459 continue;
1460 const_iterateur it;
1461 for (it=c.begin();it!=c.end();++it) {
1462 if (!is_zero(simp(*it-cand[it-c.begin()],contextptr)))
1463 break;
1464 }
1465 if (it==c.end())
1466 break;
1467 }
1468 }
1469 return ret;
1470 }
1471
1472 /* classification using the bordered Hessian and Theorem 1 of David Spring (1985) */
critical_point_class(const matrice & hess,int n,int m,GIAC_CONTEXT)1473 int critical_point_class(const matrice &hess,int n,int m,GIAC_CONTEXT) {
1474 vecteur s;
1475 int i,j,k;
1476 for (k=1;k<=n;++k) {
1477 matrice M;
1478 for (i=0;i<2*m+k;++i) {
1479 const vecteur &row=*hess[i]._VECTptr;
1480 M.push_back(vecteur(row.begin(),row.begin()+2*m+k));
1481 }
1482 s.push_back(simp(pow(gen(-1),m)*_det(M,contextptr),contextptr));
1483 }
1484 if (is_one(_contains(makesequence(s,0),contextptr)))
1485 return _CPCLASS_UNDECIDED; // paranoid check, a counterexample exists for Theorem 1 of D.S.
1486 for (i=n-1;i>=0 && is_zero(s[i]);--i);
1487 if (i<0) return _CPCLASS_UNDECIDED; // the sequence is trivial
1488 for (j=0;j<=i && is_strictly_positive(s[j],contextptr);++j);
1489 if (j>i) return i==n-1?_CPCLASS_MIN:_CPCLASS_UNDECIDED;
1490 for (j=0;j<=i && is_strictly_positive(pow(gen(-1),j+1)*s[j],contextptr);++j);
1491 if (j>i) return i==n-1?_CPCLASS_MAX:_CPCLASS_UNDECIDED;
1492 return _CPCLASS_SADDLE;
1493 }
1494
find_vars(const gen & e,const vecteur & vars,GIAC_CONTEXT)1495 vecteur find_vars(const gen &e,const vecteur &vars,GIAC_CONTEXT) {
1496 vecteur ret;
1497 for (const_iterateur it=vars.begin();it!=vars.end();++it) {
1498 if (!is_constant_wrt(e,*it,contextptr))
1499 ret.push_back(*it);
1500 }
1501 return ret;
1502 }
1503
test_parameters(const vecteur & cpt,const vecteur & vars,const vecteur & ineq,GIAC_CONTEXT)1504 bool test_parameters(const vecteur &cpt,const vecteur &vars,const vecteur &ineq,GIAC_CONTEXT) {
1505 vecteur csts=*_lname(cpt,contextptr)._VECTptr,rest_conds;
1506 gen_map conds;
1507 if (csts.empty())
1508 return true;
1509 for (const_iterateur it=cpt.begin();it!=cpt.end();++it) {
1510 vecteur fv=find_vars(*it,csts,contextptr);
1511 if (fv.empty()) continue;
1512 const gen &vr=vars[it-cpt.begin()];
1513 for (const_iterateur jt=ineq.begin();jt!=ineq.end();++jt) {
1514 if (is_constant_wrt(*jt,vr,contextptr)) continue;
1515 gen inq=symb_inferieur_strict(subst(*jt,vr,*it,false,contextptr),0);
1516 if (fv.size()==1) {
1517 gen &iqv=conds[fv.front()];
1518 if (iqv.type==_VECT)
1519 iqv._VECTptr->push_back(inq);
1520 else iqv=vecteur(1,inq);
1521 } else rest_conds.push_back(inq);
1522 }
1523 }
1524 for (gen_map::const_iterator it=conds.begin();it!=conds.end();++it) {
1525 if (it->second.type!=_VECT) continue;
1526 gen inqsol=_solve(makesequence(it->second,it->first),contextptr);
1527 if (inqsol.type==_VECT) {
1528 if (inqsol._VECTptr->empty())
1529 return false;
1530 *logptr(contextptr) << "Warning: assuming ";
1531 for (const_iterateur jt=inqsol._VECTptr->begin();jt!=inqsol._VECTptr->end();++jt) {
1532 *logptr(contextptr) << *jt;
1533 if (jt+1!=inqsol._VECTptr->end())
1534 *logptr(contextptr) << " or ";
1535 }
1536 *logptr(contextptr) << "\n";
1537 }
1538 }
1539 if (!rest_conds.empty())
1540 *logptr(contextptr) << "Warning: assuming " << rest_conds << "\n";
1541 return true;
1542 }
1543
find_local_extrema(vecteur & cpts,const gen & f,const vecteur & g,const vecteur & vars,const ipdiff::ivector & arr,const vecteur & ineq,const vecteur & initial,int order_size,bool approx_hompol,GIAC_CONTEXT)1544 void find_local_extrema(vecteur &cpts,const gen &f,const vecteur &g,const vecteur &vars,const ipdiff::ivector &arr,const vecteur &ineq,const vecteur &initial,int order_size,bool approx_hompol,GIAC_CONTEXT) {
1545 assert(order_size>=0);
1546 int nv=vars.size(),m=g.size(),n=nv-m,cls;
1547 vecteur tmpvars=make_temp_vars(vars,ineq,true,contextptr);
1548 if (order_size==0 && m>0) { // apply the method of Lagrange
1549 gen L(f);
1550 vecteur multipliers(m),allinitial;
1551 if (!initial.empty())
1552 allinitial=mergevecteur(vecteur(m,0),initial);
1553 for (int i=m;i-->0;) {
1554 L+=-(multipliers[i]=identificateur(" lambda"+print_INT_(++var_index)))*g[i];
1555 }
1556 L=subst(L,vars,tmpvars,false,contextptr);
1557 vecteur allvars=mergevecteur(tmpvars,multipliers);
1558 matrice cv,bhess;
1559 gen tmpgr=_grad(makesequence(L,allvars),contextptr);
1560 if (tmpgr.type==_VECT) {
1561 vecteur &gr=*tmpgr._VECTptr;
1562 if (allinitial.empty()) {
1563 cv=solve2(gr,allvars,contextptr);
1564 cpt_simp(cv,tmpvars,subst(f,vars,tmpvars,false,contextptr),contextptr);
1565 } else {
1566 gen tmpfsol=_fsolve(makesequence(gr,allvars,allinitial),contextptr);
1567 if (tmpfsol.type==_VECT) {
1568 vecteur &fsol=*tmpfsol._VECTptr;
1569 if (!fsol.empty())
1570 cv.push_back(fsol);
1571 }
1572 }
1573 for (iterateur it=cv.begin();it!=cv.end();++it) {
1574 *it=mergevecteur(vecteur(it->_VECTptr->begin()+nv,it->_VECTptr->end()),
1575 vecteur(it->_VECTptr->begin(),it->_VECTptr->begin()+nv));
1576 }
1577 allvars=mergevecteur(vecteur(allvars.begin()+nv,allvars.end()),
1578 vecteur(allvars.begin(),allvars.begin()+nv));
1579 if (!cv.empty()) {
1580 gen tmphess=_hessian(makesequence(L,allvars),contextptr);
1581 if (ckmatrix(tmphess))
1582 bhess=*tmphess._VECTptr; // bordered Hessian
1583 }
1584 gen s;
1585 for (const_iterateur it=cv.begin();it!=cv.end();++it) {
1586 if (!test_parameters(*(it->_VECTptr),vars,ineq,contextptr))
1587 continue;
1588 cls=_CPCLASS_UNDECIDED;
1589 if (!bhess.empty())
1590 cls=critical_point_class(subst(bhess,allvars,*it,false,contextptr),n,m,contextptr);
1591 vecteur cpt=vecteur(it->_VECTptr->begin()+m,it->_VECTptr->end());
1592 iterateur cit=find_cpt(cpts,cpt,contextptr);
1593 if (cit!=cpts.end()) {
1594 if (cls!=_CPCLASS_UNDECIDED)
1595 cit->_VECTptr->back()=cls;
1596 } else cpts.push_back(makevecteur(cpt,cls));
1597 }
1598 }
1599 } else if (order_size>0) { // use implicit differentiation instead of Lagrange multipliers
1600 vecteur gr,taylor_terms,a(nv),cpt_arr(nv);
1601 ipdiff ipd(f,g,vars,contextptr);
1602 ipd.gradient(gr);
1603 matrice cv,hess;
1604 vecteur eqv=subst(mergevecteur(gr,g),vars,tmpvars,false,contextptr);
1605 if (initial.empty()) {
1606 cv=solve2(subst(eqv,vars,tmpvars,false,contextptr),tmpvars,contextptr);
1607 cpt_simp(cv,tmpvars,subst(f,vars,tmpvars,false,contextptr),contextptr);
1608 } else {
1609 gen tmpfsol=_fsolve(makesequence(eqv,tmpvars,initial),contextptr);
1610 if (tmpfsol.type==_VECT) {
1611 vecteur &fsol=*tmpfsol._VECTptr;
1612 if (!fsol.empty())
1613 cv.push_back(fsol);
1614 }
1615 }
1616 if (!cv.empty()) {
1617 if (nv==1) {
1618 gen d;
1619 const gen &x=vars.front();
1620 for (const_iterateur it=cv.begin();it!=cv.end();++it) {
1621 gen &x0=it->_VECTptr->front();
1622 cls=_CPCLASS_UNDECIDED;
1623 for (int k=2;k<=order_size;++k) {
1624 d=simp(subst(derive(f,x,k,contextptr),x,x0,false,contextptr),contextptr);
1625 if (is_zero(d))
1626 continue;
1627 if ((k%2)!=0)
1628 cls=_CPCLASS_SADDLE;
1629 else cls=is_strictly_positive(d,contextptr)?_CPCLASS_MIN:_CPCLASS_MAX;
1630 break;
1631 }
1632 cpts.push_back(makevecteur(x0,cls));
1633 }
1634 } else {
1635 vecteur fvars(vars),ip;
1636 fvars.resize(n);
1637 identificateur l(" lambda"+print_INT_(++var_index));
1638 bool approx_hp;
1639 gen pmin,pmax,sp(-1);
1640 for (int j=0;j<n;++j) {
1641 sp+=pow(vars[j],2);
1642 ip.push_back(sqrt(fraction(1,n),contextptr));
1643 }
1644 if (order_size>1)
1645 ipd.hessian(hess);
1646 for (int i=0;i<nv;++i) {
1647 a[i]=identificateur(" a"+print_INT_(i));
1648 }
1649 for (const_iterateur it=cv.begin();it!=cv.end();++it) {
1650 if (!test_parameters(*(it->_VECTptr),vars,ineq,contextptr))
1651 continue;
1652 for (int j=0;j<nv;++j) {
1653 cpt_arr[arr[j]]=it->_VECTptr->at(j);
1654 }
1655 iterateur ct=find_cpt(cpts,cpt_arr,contextptr);
1656 if (ct==cpts.end()) {
1657 cpts.push_back(makevecteur(cpt_arr,0));
1658 ct=cpts.begin()+cpts.size()-1;
1659 }
1660 gen &cpt_class=ct->_VECTptr->back();
1661 if (order_size==1 || !is_zero(cpt_class))
1662 continue;
1663 cls=_CPCLASS_UNDECIDED;
1664 if (order_size>=2) {
1665 for (int k=2;k<=order_size;++k) {
1666 if (int(taylor_terms.size())<k-1)
1667 taylor_terms.push_back(ipd.taylor_term(a,k,false));
1668 if (is_zero(expand(taylor_terms[k-2],contextptr)))
1669 break;
1670 vecteur csts=*_lname(taylor_terms[k-2],contextptr)._VECTptr;
1671 for (int j=csts.size();j-->0;) {
1672 if (find(vars.begin(),vars.end(),csts[j])!=vars.end() ||
1673 find(a.begin(),a.end(),csts[j])!=a.end())
1674 csts.erase(csts.begin()+j);
1675 }
1676 if (!csts.empty()) {
1677 gen csol=_solve(makesequence(symb_equal(taylor_terms[k-2],0),csts),
1678 contextptr);
1679 gen sub=subst(csol,a,*it,false,contextptr);
1680 if (!is_constant_wrt_vars(sub,csts,contextptr))
1681 sub=_solve(makesequence(sub,csts),contextptr);
1682 if (sub.type==_VECT && !sub._VECTptr->empty()) {
1683 *logptr(contextptr)
1684 << "Warning: assuming " << csts << " not in "
1685 << simplify(sub,contextptr) << "\n";
1686 }
1687 }
1688 gen p=simp(subst(taylor_terms[k-2],a,*it,false,contextptr),contextptr);
1689 if (is_zero(p))
1690 continue;
1691 else if (k%2)
1692 cls=_CPCLASS_SADDLE;
1693 else {
1694 gen gL,tmp(undef);
1695 if ((approx_hp=approx_hompol)) {
1696 vecteur hpvars=*_lname(p,contextptr)._VECTptr;
1697 const_iterateur hit=hpvars.begin();
1698 for (;hit!=hpvars.end() &&
1699 find(fvars.begin(),fvars.end(),*hit)!=fvars.end();++hit);
1700 if (hit==hpvars.end()) {
1701 gL=_fMin(makesequence(p,vecteur(1,symb_equal(sp,0)),fvars,ip),
1702 contextptr);
1703 if (gL.type==_VECT && int(gL._VECTptr->size())==n) {
1704 tmp=vecteur(1,_epsilon2zero(subst(p,fvars,*gL._VECTptr,
1705 false,contextptr),
1706 contextptr));
1707 gL=_fMax(makesequence(p,vecteur(1,symb_equal(sp,0)),
1708 fvars,ip),contextptr);
1709 if (gL.type==_VECT && int(gL._VECTptr->size())==n)
1710 tmp._VECTptr->push_back(_epsilon2zero(
1711 subst(p,fvars,*gL._VECTptr,false,contextptr),
1712 contextptr));
1713 else tmp=undef;
1714 }
1715 }
1716 }
1717 if (is_undef(tmp)) {
1718 approx_hp=false;
1719 fvars.push_back(l);
1720 gL=_grad(makesequence(p-l*sp,fvars),contextptr);
1721 if (gL.type==_VECT)
1722 tmp=_solve(makesequence(gL,fvars),contextptr);
1723 }
1724 if (tmp.type!=_VECT || tmp._VECTptr->empty() ||
1725 (!approx_hp && !ckmatrix(tmp))) {
1726 if (k>2) break;
1727 // apply the second order derivative test
1728 matrice chess=*simp(subst(hess,vars,*it,false,contextptr),
1729 contextptr)._VECTptr;
1730 if (is_undef(chess))
1731 break;
1732 if (_lname(chess,contextptr)._VECTptr->empty()) // no symbols here
1733 chess=*_evalf(chess,contextptr)._VECTptr;
1734 gen res=_eigenvals(chess,contextptr);
1735 if (res.type==_VECT) {
1736 vecteur eigvals=*res._VECTptr;
1737 gen e=undef;
1738 for (const_iterateur et=eigvals.begin();et!=eigvals.end();++et) {
1739 if (is_zero(*et,contextptr)) {
1740 cls=_CPCLASS_UNDECIDED;
1741 break;
1742 } else if (is_undef(e)) {
1743 e=*et;
1744 cls=is_positive(e,contextptr)?_CPCLASS_MIN:_CPCLASS_MAX;
1745 } else if (is_strictly_positive(-e*(*et),contextptr))
1746 cls=_CPCLASS_SADDLE;
1747 }
1748 }
1749 break;
1750 }
1751 if (approx_hp) {
1752 pmin=tmp._VECTptr->front();
1753 pmax=tmp._VECTptr->back();
1754 } else {
1755 vecteur lst;
1756 const_iterateur mt=tmp._VECTptr->begin();
1757 for (;mt!=tmp._VECTptr->end();++mt) {
1758 lst.push_back(simp(subst(p,fvars,*mt->_VECTptr,
1759 false,contextptr),
1760 contextptr));
1761 }
1762 pmin=_min(lst,contextptr);
1763 pmax=_max(lst,contextptr);
1764 fvars.pop_back();
1765 }
1766 if (is_zero(pmin,contextptr))
1767 cls=_CPCLASS_POSSIBLE_MIN;
1768 else if (is_zero(pmax,contextptr))
1769 cls=_CPCLASS_POSSIBLE_MAX;
1770 else if (is_strictly_positive(pmin,contextptr))
1771 cls=_CPCLASS_MIN;
1772 else if (is_strictly_positive(pmax,contextptr))
1773 cls=_CPCLASS_SADDLE;
1774 else
1775 cls=_CPCLASS_MAX;
1776 }
1777 break;
1778 }
1779 }
1780 cpt_class=gen(cls);
1781 }
1782 }
1783 }
1784 }
1785 for (const_iterateur it=tmpvars.begin();it!=tmpvars.end();++it) {
1786 if (find(vars.begin(),vars.end(),*it)==vars.end())
1787 _purge(*it,contextptr);
1788 }
1789 }
1790
1791 /*
1792 * 'extrema' attempts to find all points of strict local minima/maxima of a
1793 * smooth (uni/multi)variate function (possibly subject to equality constraints).
1794 *
1795 * Usage
1796 * ^^^^^
1797 * extrema(expr,[constr],vars,[order=n||lagrange])
1798 *
1799 * Parameters
1800 * ^^^^^^^^^^
1801 * - expr : differentiable expression
1802 * - constr (optional) : (list of) equality constraint(s)
1803 * - vars : (list of) problem variable(s)
1804 * - order_size (optional) : specify 'order_size=<nonnegative integer>' to
1805 * bound the order of the derivatives being
1806 * inspected when classifying the critical points
1807 *
1808 * The number of constraints must be less than the number of variables. When
1809 * there are more than one constraint/variable, they must be specified in
1810 * form of list.
1811 *
1812 * Variables may be specified with symbol, e.g. 'var', or by using syntax
1813 * 'var=a..b', which restricts the variable 'var' to the open interval (a,b),
1814 * where a and b are real numbers or +/-infinity. If variable list includes a
1815 * specification of initial point, such as, for example, [x=1,y=0,z=2], then
1816 * numeric solver is activated to find critical point in the vicinity of the
1817 * given point. In this case, the single critical point, if found, is examined.
1818 *
1819 * The function attempts to find the critical points in exact form, if the
1820 * parameters of the problem are all exact. It works best for problems in which
1821 * the gradient of lagrangian function consists of rational expressions. The
1822 * result may be inexact, however, if exact solutions could not be obtained.
1823 *
1824 * For classifying critical points, the bordered hessian method is used first.
1825 * It is only a second order test, so it may be inconclusive in some cases. In
1826 * these cases function looks at higher-order derivatives, up to order
1827 * specified by 'order_size' option (the extremum test). Set 'order_size' to 1
1828 * to use only the bordered hessian test or 0 to output critical points without
1829 * attempting to classify them. Setting 'order_size' to 2 or higher will
1830 * activate checking for saddle points and inspecting higher derivatives (up to
1831 * 'order_size') to determine the nature of some or all unclassified critical
1832 * points. By default 'order_size' equals to 5.
1833 *
1834 * The return value is a sequence with two elements: list of strict local
1835 * minima and list of strict local maxima. If only critical points are
1836 * requested (by setting 'order_size' to 0), the output consists of a single
1837 * list, as no classification was attempted. For univariate problems the points
1838 * are real numbers, while for multivariate problems they are specified as
1839 * lists of coordinates, so "lists of points" are in fact matrices with rows
1840 * corresponding to points in multivariate cases, i.e. vectors in univariate
1841 * cases.
1842 *
1843 * The function prints out information about saddle/inflection points and also
1844 * about critical points for which no decision could be made, so that the user
1845 * can inspect candidates for local extrema by plotting the graph, for example.
1846 *
1847 * Examples
1848 * ^^^^^^^^
1849 * extrema(-2*cos(x)-cos(x)^2,x)
1850 * >> [0],[pi]
1851 * extrema((x^3-1)^4/(2x^3+1)^4,x=0..inf)
1852 * >> [1],[]
1853 * extrema(x/2-2*sin(x/2),x=-12..12)
1854 * >> [2*pi/3,-10*pi/3],[10*pi/3,-2*pi/3]
1855 * extrema(x-ln(abs(x)),x)
1856 * >> [1],[]
1857 * assume(a>=0);extrema(x^2+a*x,x)
1858 * >> [-a/2],[]
1859 * extrema(x^7+3x^6+3x^5+x^4+2x^2-x,x)
1860 * >> [0.225847362349],[-1.53862319761]
1861 * extrema((x^2+x+1)/(x^4+1),x)
1862 * >> [],[0.697247087784]
1863 * extrema(x^2+exp(-x),x)
1864 * >> [LambertW(1/2)],[]
1865 * extrema(exp(-x)*ln(x),x)
1866 * >> [],[exp(LambertW(1))]
1867 * extrema(tan(x)*(x^3-5x^2+1),x=-0.5)
1868 * >> [-0.253519032024],[]
1869 * extrema(tan(x)*(x^3-5x^2+1),x=0.5)
1870 * >> [],[0.272551772027]
1871 * extrema(exp(x^2-2x)*ln(x)*ln(1-x),x=0.5)
1872 * >> [],[0.277769149124]
1873 * extrema(ln(2+x-sin(x)^2),x=0..2*pi)
1874 * >> [],[] (needed to compute third derivative to drop critical points pi/4 and 5pi/4)
1875 * extrema(x^3-2x*y+3y^4,[x,y])
1876 * >> [[12^(1/5)/3,(12^(1/5))^2/6]],[]
1877 * extrema((2x^2-y)*(y-x^2),[x,y]) //Peano surface
1878 * >> [],[] (saddle point at origin)
1879 * extrema(5x^2+3y^2+x*z^2-z*y^2,[x,y,z])
1880 * >> [],[] (possible local minimum at origin, in fact saddle)
1881 * extrema(3*atan(x)-2*ln(x^2+y^2+1),[x,y])
1882 * >> [],[[3/4,0]]
1883 * extrema(x*y,x+y=1,[x,y])
1884 * >> [],[[1/2,1/2]]
1885 * extrema(sqrt(x*y),x+y=2,[x,y])
1886 * >> [],[[1,1]]
1887 * extrema(x*y,x^3+y^3=16,[x,y])
1888 * >> [],[[2,2]]
1889 * extrema(x^2+y^2,x*y=1,[x=0..inf,y=0..inf])
1890 * >> [[1,1]],[]
1891 * extrema(ln(x*y^2),2x^2+3y^2=8,[x,y])
1892 * >> [],[[2*sqrt(3)/3,-4/3],[2*sqrt(3)/3,4/3]]
1893 * extrema(y^2+4y+2x-x^2,x+2y=2,[x,y])
1894 * >> [],[[-2/3,4/3]]
1895 * assume(a>0);extrema(x/a^2+a*y^2,x+y=a,[x,y])
1896 * >> [[(2*a^4-1)/(2*a^3),1/(2*a^3)]],[]
1897 * extrema(6x+3y+2z,4x^2+2y^2+z^2=70,[x,y,z])
1898 * >> [[-3,-3,-4]],[[3,3,4]]
1899 * extrema(x*y*z,x+y+z=1,[x,y,z])
1900 * >> [],[[1/3,1/3,1/3]]
1901 * extrema(x*y^2*z^2,x+y+z=5,[x,y,z])
1902 * >> [],[[1,2,2]]
1903 * extrema(4y-2z,[2x-y-z=2,x^2+y^2=1],[x,y,z])
1904 * >> [[2*sqrt(13)/13,-3*sqrt(13)/13,(7*sqrt(13)-26)/13]],
1905 * [[-2*sqrt(13)/13,3*sqrt(13)/13,(-7*sqrt(13)-26)/13]]
1906 * extrema((x-3)^2+(y-1)^2+(z-1)^2,x^2+y^2+z^2=4,[x,y,z])
1907 * >> [[6*sqrt(11)/11,2*sqrt(11)/11,2*sqrt(11)/11]],
1908 * [[-6*sqrt(11)/11,-2*sqrt(11)/11,-2*sqrt(11)/11]]
1909 * extrema(x+3y-z,2x^2+y^2=z,[x,y,z])
1910 * >> [],[[1/4,3/2,19/8]]
1911 * extrema(2x*y+2y*z+x*z,x*y*z=4,[x,y,z])
1912 * >> [[2,1,2]],[]
1913 * extrema(x+y+z,[x^2+y^2=1,2x+z=1],[x,y,z])
1914 * >> [[sqrt(2)/2,-sqrt(2)/2,-sqrt(2)+1]],[[-sqrt(2)/2,sqrt(2)/2,sqrt(2)+1]]
1915 * assume(a>0);extrema(x+y+z,[y^2-x^2=a,x+2z=1],[x,y,z])
1916 * >> [[-sqrt(3)*sqrt(a)/3,2*sqrt(3)*sqrt(a)/3,(sqrt(3)*sqrt(a)+3)/6]],
1917 * [[sqrt(3)*sqrt(a)/3,-2*sqrt(3)*sqrt(a)/3,(-sqrt(3)*sqrt(a)+3)/6]]
1918 * extrema((x-u)^2+(y-v)^2,[x^2/4+y^2/9=1,(u-3)^2+(v+5)^2=1],[u,v,x,y])
1919 * >> [[2.35433932354,-4.23637555425,0.982084902545,-2.61340692712]],
1920 * [[3.41406613851,-5.91024679428,-0.580422508346,2.87088778158]]
1921 * extrema(x2^6+x1^3+4x1+4x2,x1^5+x2^4+x1+x2=0,[x1,x2])
1922 * >> [[-0.787457596325,0.758772338924],[-0.784754836317,-1.23363062357]],
1923 * [[0.154340184382,-0.155005038065]]
1924 * extrema(x*y,-2x^3+15x^2*y+11y^3-24y=0,[x,y])
1925 * >> [[sqrt(17)*sqrt(3*sqrt(33)+29)/17,sqrt(-15*sqrt(33)+127)*sqrt(187)/187],
1926 * [-sqrt(17)*sqrt(3*sqrt(33)+29)/17,-sqrt(-15*sqrt(33)+127)*sqrt(187)/187],
1927 * [sqrt(-3*sqrt(33)+29)*sqrt(17)/17,-sqrt(15*sqrt(33)+127)*sqrt(187)/187],
1928 * [-sqrt(-3*sqrt(33)+29)*sqrt(17)/17,sqrt(15*sqrt(33)+127)*sqrt(187)/187]],
1929 * [[1,1],[-1,-1],[0,0]]
1930 * extrema(x2^4-x1^4-x2^8+x1^10,[x1,x2],order_size=1)
1931 * >> [[0,0],[0,-(1/2)^(1/4)],[0,(1/2)^(1/4)],[-(2/5)^(1/6),0],[-(2/5)^(1/6),-(1/2)^(1/4)],
1932 * [-(2/5)^(1/6),(1/2)^(1/4)],[(2/5)^(1/6),0],[(2/5)^(1/6),-(1/2)^(1/4)],[(2/5)^(1/6),(1/2)^(1/4)]]
1933 * extrema(x2^4-x1^4-x2^8+x1^10,[x1,x2])
1934 * >> [[],[]]
1935 * extrema(x2^6+x1^3+4x1+4x2,x1^5+x2^4+x1+x2=0,[x1,x2])
1936 * >> [[-0.787457596325,0.758772338924],[-0.784754836317,-1.23363062357]],
1937 * [[0.154340184382,-0.155005038065]]
1938 * extrema(x2^6+x1^3+2x1^2-x2^2+4x1+4x2,x1^5+x2^4+x1+x2=0,[x1,x2])
1939 * >> [[-0.662879934158,-1.18571027742],[0,0]],[[0.301887394815,-0.314132376868]]
1940 * extrema(3x^2-2x*y+y^2-8y,[x,y])
1941 * >> [[2,6]],[]
1942 * extrema(x^3+3x*y^2-15x-12y,[x,y])
1943 * >> [[2,1]],[[-2,-1]]
1944 * extrema(4x*y-x^4-y^4,[x,y])
1945 * >> [],[[1,1],[-1,-1]]
1946 * extrema(x*sin(y),[x,y])
1947 * >> [],[]
1948 * extrema(x^4+y^4,[x,y])
1949 * >> [[0,0]],[]
1950 * extrema(x^3*y-x*y^3,[x,y]) //dog saddle at origin
1951 * >> [],[]
1952 * extrema(x^2+y^2+z^2,x^4+y^4+z^4=1,[x,y,z])
1953 * >> [[0,0,1],[0,0,-1]],[]
1954 * extrema(3x+3y+8z,[x^2+z^2=1,y^2+z^2=1],[x,y,z])
1955 * >> [[-3/5,-3/5,-4/5]],[[3/5,3/5,4/5]]
1956 * extrema(2x^2+y^2,x^4-x^2+y^2=5,[x,y])
1957 * >> [[0,-sqrt(5)],[0,sqrt(5)]],
1958 * [[-1/2*sqrt(6),-1/2*sqrt(17)],[-1/2*sqrt(6),1/2*sqrt(17)],
1959 * [1/2*sqrt(6),-1/2*sqrt(17)],[1/2*sqrt(6),1/2*sqrt(17)]]
1960 * extrema((3x^4-4x^3-12x^2+18)/(12*(1+4y^2)),[x,y])
1961 * >> [[2,0]],[[0,0]]
1962 * extrema(x-y+z,[x^2+y^2+z^2=1,x+y+2z=1],[x,y,z])
1963 * >> [[(-2*sqrt(70)+7)/42,(4*sqrt(70)+7)/42,(-sqrt(70)+14)/42]],
1964 * [[(2*sqrt(70)+7)/42,(-4*sqrt(70)+7)/42,(sqrt(70)+14)/42]]
1965 * extrema(ln(x)+2*ln(y)+3*ln(z)+4*ln(u)+5*ln(v),x+y+z+u+v=1,[x,y,z,u,v])
1966 * >> [],[[1/15,2/15,1/5,4/15,1/3]]
1967 * extrema(x*y*z,-2x^3+15x^2*y+11y^3-24y=0,[x,y,z])
1968 * >> [],[]
1969 * extrema(x+y-exp(x)-exp(y)-exp(x+y),[x,y])
1970 * >> [],[[ln(-1/2*(-sqrt(5)+1)),ln(-1/2*(-sqrt(5)+1))]]
1971 * extrema(x^2*sin(y)-4*x,[x,y]) // has two saddle points
1972 * >> [],[]
1973 * extrema((1+y*sinh(x))/(1+y^2+tanh(x)^2),[x,y])
1974 * >> [],[[0,0]]
1975 * extrema((1+y*sinh(x))/(1+y^2+tanh(x)^2),y=x^2,[x,y])
1976 * >> [[1.42217627369,2.02258535346]],[[8.69443642205e-16,7.55932246971e-31]]
1977 * extrema(x^2*y^2,[x^3*y-2*x^2+3x-2y^2=0],[x=0..inf,y])
1978 * >> [[3/2,0]],[[0.893768095046,-0.5789326693514]]
1979 * extrema(alpha*x^2+beta*y^2+gamma*z^2,[a1*x+a2*y+a3*z=c,b1*x+b2*y+b3*z=d],[x,y,z])
1980 * assume(a>0):;extrema(sqrt(a+x+y)/(1+y+sqrt(a+x)),[x=0..inf,y=0..inf])
1981 */
_extrema(const gen & g,GIAC_CONTEXT)1982 gen _extrema(const gen &g,GIAC_CONTEXT) {
1983 if (g.type==_STRNG && g.subtype==-1) return g;
1984 if (g.type!=_VECT || g.subtype!=_SEQ__VECT)
1985 return gentypeerr(contextptr);
1986 vecteur &gv=*g._VECTptr,constr;
1987 int order_size=5; // will not compute the derivatives of order higher than 'order_size'
1988 int ngv=gv.size();
1989 bool approx_hp=true; // use nlpsolve to determine images of homogeneous polynomials on spheres
1990 if (gv.back()==at_lagrange) {
1991 order_size=0; // use Lagrange method
1992 --ngv;
1993 } else if (gv.back().is_symb_of_sommet(at_equal)) {
1994 vecteur &v=*gv.back()._SYMBptr->feuille._VECTptr;
1995 if (v[0]==at_order && is_integer(v[1])) {
1996 if ((order_size=v[1].val)<1)
1997 return gensizeerr("Expected a positive integer");
1998 --ngv;
1999 }
2000 }
2001 if (ngv<2 || ngv>3)
2002 return gensizeerr("Wrong number of input arguments");
2003 // get the variables
2004 int nv;
2005 vecteur vars,ineq,initial;
2006 // parse variables and their ranges, if given
2007 if ((nv=parse_varlist(gv[ngv-1],vars,ineq,initial,contextptr))==0)
2008 return gentypeerr("Failed to parse variables");
2009 if (nv>32)
2010 return gendimerr("Too many variables");
2011 if (!initial.empty() && int(initial.size())<nv)
2012 return gendimerr("Bad initial point assumption");
2013 if (ngv==3) {
2014 // get the constraints
2015 if (gv[1].type==_VECT)
2016 constr=*gv[1]._VECTptr;
2017 else
2018 constr=vecteur(1,gv[1]);
2019 }
2020 if (order_size==0 && constr.empty())
2021 return gensizeerr("At least one constraint is required for Lagrange method");
2022 for (iterateur it=constr.begin();it!=constr.end();++it) {
2023 if (it->is_symb_of_sommet(at_equal))
2024 *it=equal2diff(*it);
2025 }
2026 ipdiff::ivectors arrs;
2027 if (order_size>0 && !constr.empty()) {
2028 matrice J(jacobian(constr,vars,contextptr));
2029 if (J.empty())
2030 return gensizeerr(contextptr);
2031 if (constr.size()>=vars.size() || _rank(J,contextptr).val<int(constr.size()))
2032 return gendimerr("Too many constraints");
2033 vars_arrangements(J,arrs,contextptr);
2034 } else {
2035 ipdiff::ivector arr(nv);
2036 for (int i=0;i<nv;++i) {
2037 arr[i]=i;
2038 }
2039 arrs.push_back(arr);
2040 }
2041 vecteur tmp_vars(vars.size()),cpts;
2042 /* iterate through all possible variable arrangements */
2043 bool undetected=false;
2044 for (ipdiff::ivectors::const_iterator ait=arrs.begin();ait!=arrs.end();++ait) {
2045 const ipdiff::ivector &arr=*ait;
2046 for (ipdiff::ivector::const_iterator it=arr.begin();it!=arr.end();++it) {
2047 tmp_vars[it-arr.begin()]=vars[*it];
2048 }
2049 try {
2050 find_local_extrema(cpts,gv[0],constr,tmp_vars,arr,ineq,initial,order_size,approx_hp,contextptr);
2051 } catch (const std::exception &e) {
2052 undetected=true;
2053 }
2054 }
2055 if (undetected) {
2056 string msg="Warning: some critical points may have been undetected.";
2057 if (initial.empty())
2058 msg+=" Try setting an initial point for numerical approximation.";
2059 *logptr(contextptr) << msg << "\n";
2060 }
2061 if (order_size==1) { // return the list of critical points
2062 vecteur cv;
2063 for (const_iterateur it=cpts.begin();it!=cpts.end();++it) {
2064 cv.push_back(it->_VECTptr->front());
2065 }
2066 return cv;
2067 }
2068 // return sequence of minima and maxima in separate lists and report non- or possible extrema
2069 vecteur minv(0),maxv(0);
2070 for (const_iterateur it=cpts.begin();it!=cpts.end();++it) {
2071 gen dispt(nv==1?symb_equal(vars[0],it->_VECTptr->front()):_zip(makesequence(at_equal,vars,it->_VECTptr->front()),contextptr));
2072 switch(it->_VECTptr->back().val) {
2073 case _CPCLASS_MIN:
2074 minv.push_back(it->_VECTptr->front());
2075 break;
2076 case _CPCLASS_MAX:
2077 maxv.push_back(it->_VECTptr->front());
2078 break;
2079 case _CPCLASS_SADDLE:
2080 *logptr(contextptr) << dispt << (nv==1?": inflection point":": saddle point") << "\n";
2081 break;
2082 case _CPCLASS_POSSIBLE_MIN:
2083 *logptr(contextptr) << dispt << ": indeterminate critical point\n";
2084 break;
2085 case _CPCLASS_POSSIBLE_MAX:
2086 *logptr(contextptr) << dispt << ": indeterminate critical point\n";
2087 break;
2088 case _CPCLASS_UNDECIDED:
2089 *logptr(contextptr) << dispt << ": unclassified critical point\n";
2090 break;
2091 }
2092 }
2093 return makevecteur(minv,maxv);
2094 }
2095 static const char _extrema_s []="extrema";
2096 static define_unary_function_eval (__extrema,&_extrema,_extrema_s);
2097 define_unary_function_ptr5(at_extrema,alias_at_extrema,&__extrema,0,true)
2098
2099 /*
2100 * Compute the value of expression f(var) (or |f(var)| if 'absolute' is true)
2101 * for var=a.
2102 */
compf(const gen & f,identificateur & x,gen & a,bool absolute,GIAC_CONTEXT)2103 gen compf(const gen &f,identificateur &x,gen &a,bool absolute,GIAC_CONTEXT) {
2104 gen val(subst(f,x,a,false,contextptr));
2105 return _evalf(absolute?_abs(val,contextptr):val,contextptr);
2106 }
2107
2108 /*
2109 * find zero of expression f(x) for x in [a,b] using Brent solver
2110 */
find_zero(const gen & f,identificateur & x,gen & a,gen & b,GIAC_CONTEXT)2111 gen find_zero(const gen &f,identificateur &x,gen &a,gen &b,GIAC_CONTEXT) {
2112 gen I(symb_interval(a,b));
2113 gen var(symb_equal(x,I));
2114 gen tmpsol=_fsolve(makesequence(f,var,_BRENT_SOLVER),contextptr);
2115 if (tmpsol.type==_VECT) {
2116 vecteur &sol=*tmpsol._VECTptr;
2117 return sol.empty()?(a+b)/2:sol[0];
2118 } else return (a+b)/2;
2119 }
2120
2121 /*
2122 * Find point of maximum/minimum of unimodal expression f(x) in [a,b] using the
2123 * golden-section search.
2124 */
find_peak(const gen & f,identificateur & x,gen & a_orig,gen & b_orig,GIAC_CONTEXT)2125 gen find_peak(const gen &f,identificateur &x,gen &a_orig,gen &b_orig,GIAC_CONTEXT) {
2126 gen a(a_orig),b(b_orig);
2127 gen c(b-(b-a)/GOLDEN_RATIO),d(a+(b-a)/GOLDEN_RATIO);
2128 while (is_strictly_greater(_abs(c-d,contextptr),epsilon(contextptr),contextptr)) {
2129 gen fc(compf(f,x,c,true,contextptr)),fd(compf(f,x,d,true,contextptr));
2130 if (is_strictly_greater(fc,fd,contextptr))
2131 b=d;
2132 else
2133 a=c;
2134 c=b-(b-a)/GOLDEN_RATIO;
2135 d=a+(b-a)/GOLDEN_RATIO;
2136 }
2137 return (a+b)/2;
2138 }
2139
2140 /*
2141 * Compute n Chebyshev nodes in [a,b].
2142 */
chebyshev_nodes(gen & a,gen & b,int n,GIAC_CONTEXT)2143 vecteur chebyshev_nodes(gen &a,gen &b,int n,GIAC_CONTEXT) {
2144 vecteur nodes(1,a);
2145 for (int i=1;i<=n;++i) {
2146 nodes.push_back(_evalf((a+b)/2+(b-a)*symbolic(at_cos,((2*i-1)*cst_pi/(2*n)))/2,contextptr));
2147 }
2148 nodes.push_back(b);
2149 return *_sort(nodes,contextptr)._VECTptr;
2150 }
2151
2152 /*
2153 * Implementation of Remez method for minimax polynomial approximation of a
2154 * continuous bounded function, which is not necessary differentiable in all
2155 * points of (a,b).
2156 *
2157 * Source: http://homepages.rpi.edu/~tasisa/remez.pdf
2158 *
2159 * Usage
2160 * ^^^^^
2161 * minimax(expr,var=a..b,n,[opts])
2162 *
2163 * Parameters
2164 * ^^^^^^^^^^
2165 * - expr : expression to be approximated on [a,b]
2166 * - var : variable
2167 * - a,b : bounds of the function domain
2168 * - n : degree of the minimax approximation polynomial
2169 * - opts (optional) : sequence of options
2170 *
2171 * This function uses 'compf', 'find_zero' and 'find_peak'. It does not use
2172 * derivatives to determine points of local extrema of error function, but
2173 * instead implements the golden search algorithm to find these points in the
2174 * exchange phase of Remez method.
2175 *
2176 * The returned polynomial may have degree lower than n, because the latter is
2177 * decremented during execution of the algorithm if there is no unique solution
2178 * for coefficients of a nth degree polynomial. After each decrement, the
2179 * algorithm is restarted. If the degree of resulting polynomial is m<n, it
2180 * means that polynomial of degree between m and n cannot be obtained by using
2181 * this implementation.
2182 *
2183 * In 'opts' one may specify 'limit=<posint>' which limits the number of
2184 * iterations. By default, it is unlimited.
2185 *
2186 * Be aware that, in some cases, the result with high n may be unsatisfying,
2187 * producing larger error than the polynomials for smaller n. This happens
2188 * because of the rounding errors. Nevertheless, a good approximation of an
2189 * "almost" smooth function can usually be obtained with n less than 30. Highly
2190 * oscillating functions containing sharp cusps and spikes will probably be
2191 * approximated poorly.
2192 *
2193 * Examples
2194 * ^^^^^^^^
2195 * minimax(x*exp(-x),x=0..10,24)
2196 * minimax(x*sin(x),x=0..10,25)
2197 * minimax(ln(2+x-sin(x)^2),x=0..2*pi,20)
2198 * minimax(cos(x^2-x+1),x=-2..2,40)
2199 * minimax(atan(x),x=-5..5,25)
2200 * minimax(tanh(sin(9x)),x=-1/2..1/2,40)
2201 * minimax(abs(x),x=-1..1,20)
2202 * minimax(abs(x)*sqrt(abs(x)),x=-2..2,15)
2203 * minimax(min(1/cosh(3*sin(10x)),sin(9x)),x=-0.3..0.4,25)
2204 * minimax(when(x==0,0,exp(-1/x^2)),x=-1..1,30)
2205 */
_minimax(const gen & g,GIAC_CONTEXT)2206 gen _minimax(const gen &g,GIAC_CONTEXT) {
2207 if (g.type==_STRNG && g.subtype==-1) return g;
2208 if (g.type!=_VECT || g.subtype!=_SEQ__VECT)
2209 return gentypeerr(contextptr);
2210 vecteur &gv=*g._VECTptr;
2211 if (gv.size()<3)
2212 return gensizeerr(contextptr);
2213 if (!gv[1].is_symb_of_sommet(at_equal) || !is_integer(gv[2]))
2214 return gentypeerr(contextptr);
2215 // detect parameters
2216 vecteur s(*gv[1]._SYMBptr->feuille._VECTptr);
2217 if (s[0].type!=_IDNT || !s[1].is_symb_of_sommet(at_interval))
2218 return gentypeerr((contextptr));
2219 identificateur x(*s[0]._IDNTptr);
2220 s=*s[1]._SYMBptr->feuille._VECTptr;
2221 gen a(_evalf(s[0],contextptr)),b(_evalf(s[1],contextptr));
2222 if (!is_strictly_greater(b,a,contextptr))
2223 return gentypeerr(contextptr);
2224 gen &f=gv[0];
2225 int n=gv[2].val;
2226 gen threshold(1.02); // threshold for stopping criterion
2227 // detect options
2228 int limit=0;
2229 //bool poly=true;
2230 for (const_iterateur it=gv.begin()+3;it!=gv.end();++it) {
2231 if (it->is_symb_of_sommet(at_equal)) {
2232 vecteur &p=*it->_SYMBptr->feuille._VECTptr;
2233 if (p[0]==at_limit) {
2234 if (!is_integer(p[1]) || !is_strictly_positive(p[1],contextptr))
2235 return gentypeerr(contextptr);
2236 limit=p[1].val;
2237 }
2238 }
2239 else if (is_integer(*it)) {
2240 switch (it->val) {
2241 // case _FRAC:
2242 // poly=false;
2243 // break;
2244 }
2245 }
2246 }
2247 // create Chebyshev nodes to start with
2248 vecteur nodes(chebyshev_nodes(a,b,n,contextptr)),sol;
2249 gen p,best_p,best_emax,emax,emin,tmpsol;
2250 int iteration_count=0;
2251 while (true) { // iterate the algorithm
2252 iteration_count++;
2253 if (n<1 || (limit>0 && iteration_count>limit))
2254 break;
2255 // compute polynomial p
2256 matrice m;
2257 vecteur fv;
2258 for (int i=0;i<n+2;++i) {
2259 fv.push_back(_evalf(subst(f,x,nodes[i],false,contextptr),contextptr));
2260 vecteur r;
2261 for (int j=0;j<n+1;++j) {
2262 r.push_back(j==0?gen(1):pow(nodes[i],j));
2263 }
2264 r.push_back(pow(gen(-1),i));
2265 m.push_back(r);
2266 }
2267 tmpsol=_linsolve(makesequence(_epsilon2zero(m,contextptr),fv),contextptr);
2268 if (tmpsol.type==_VECT)
2269 sol=*tmpsol._VECTptr;
2270 else sol.clear();
2271 if (sol.empty() || !_lname(sol,contextptr)._VECTptr->empty()) {
2272 // Solution is not unique, it contains a symbol.
2273 // Decrease n and start over.
2274 nodes=chebyshev_nodes(a,b,--n,contextptr);
2275 continue;
2276 }
2277 p=gen(0);
2278 for (int i=0;i<n+1;++i) {
2279 p+=sol[i]*pow(x,i);
2280 }
2281 // compute the error function and its zeros
2282 gen e(f-p);
2283 vecteur zv(1,a);
2284 for (int i=0;i<n+1;++i) {
2285 zv.push_back(find_zero(e,x,nodes[i],nodes[i+1],contextptr));
2286 }
2287 zv.push_back(b);
2288 // remez exchange:
2289 // determine points of local extrema of error function e
2290 vecteur ev(n+2,0);
2291 for (int i=0;i<n+2;++i) {
2292 if (i>0 && i<n+1) {
2293 nodes[i]=find_peak(e,x,zv[i],zv[i+1],contextptr);
2294 ev[i]=compf(e,x,nodes[i],true,contextptr);
2295 continue;
2296 }
2297 gen e1(compf(e,x,zv[i],true,contextptr)),e2(compf(e,x,zv[i+1],true,contextptr));
2298 if (is_greater(e1,e2,contextptr)) {
2299 nodes[i]=zv[i];
2300 ev[i]=e1;
2301 }
2302 else {
2303 nodes[i]=zv[i+1];
2304 ev[i]=e2;
2305 }
2306 }
2307 // compute minimal and maximal absolute error
2308 emin=_min(ev,contextptr);
2309 emax=_max(ev,contextptr);
2310 if (is_exactly_zero(best_emax) || is_strictly_greater(best_emax,emax,contextptr)) {
2311 best_p=p;
2312 best_emax=emax;
2313 }
2314 // emin >= E is required to continue, also check
2315 // if the threshold is reached, i.e. the difference
2316 // between emax and emin is at least fifty times
2317 // smaller than emax
2318 if (is_strictly_greater(sol.back(),emin,contextptr) ||
2319 is_greater(threshold*emin,emax,contextptr)) {
2320 break;
2321 }
2322 }
2323 *logptr(contextptr) << "max. absolute error: " << best_emax << "\n";
2324 return best_p;
2325 }
2326 static const char _minimax_s []="minimax";
2327 static define_unary_function_eval (__minimax,&_minimax,_minimax_s);
2328 define_unary_function_ptr5(at_minimax,alias_at_minimax,&__minimax,0,true)
2329
2330 /*
2331 * TPROB CLASS IMPLEMENTATION
2332 */
2333
tprob(const vecteur & s,const vecteur & d,const gen & m,GIAC_CONTEXT)2334 tprob::tprob(const vecteur &s,const vecteur &d,const gen &m,GIAC_CONTEXT) {
2335 eps=exact(epsilon(contextptr)/2,contextptr);
2336 ctx=contextptr;
2337 supply=s;
2338 demand=d;
2339 M=m;
2340 }
2341
2342 /*
2343 * North-West-Corner method giving the initial feasible solution to the
2344 * transportatiom problem with given supply and demand vectors. It handles
2345 * degeneracy cases (assignment problems, for example, always have degenerate
2346 * solutions).
2347 */
north_west_corner(matrice & feas)2348 void tprob::north_west_corner(matrice &feas) {
2349 feas.clear();
2350 int m=supply.size(),n=demand.size();
2351 for (int k=0;k<m;++k) {
2352 feas.push_back(vecteur(n,0));
2353 }
2354 int i=0,j=0;
2355 while (i<m && j<n) {
2356 const gen &s=supply[i],&d=demand[j];
2357 gen u,v;
2358 for (int k=0;k<i;++k) {
2359 v+=_epsilon2zero(feas[k][j],ctx);
2360 }
2361 for (int k=0;k<j;++k) {
2362 u+=_epsilon2zero(feas[i][k],ctx);
2363 }
2364 gen a=min(s-u,d-v,ctx);
2365 feas[i]._VECTptr->at(j)=a;
2366 int k=i+j;
2367 if (u+a==s)
2368 ++i;
2369 if (v+a==d)
2370 ++j;
2371 if (i<m && j<n && i+j-k==2) // avoid degeneracy
2372 feas[i-1]._VECTptr->at(j)=eps;
2373 }
2374 }
2375
2376 /*
2377 * Stepping stone path method for determining a closed path "jumping" from one
2378 * positive element of X to another in the same row or column.
2379 */
stepping_stone_path(ipairs & path_orig,const matrice & X)2380 tprob::ipairs tprob::stepping_stone_path(ipairs &path_orig,const matrice &X) {
2381 ipairs path(path_orig);
2382 int I=path.back().first,J=path.back().second;
2383 int m=X.size(),n=X.front()._VECTptr->size();
2384 if (path.size()>1 && path.front().second==J)
2385 return path;
2386 bool hrz=path.size()%2==1;
2387 for (int i=0;i<(hrz?n:m);++i) {
2388 int cnt=0;
2389 for (ipairs::const_iterator it=path.begin();it!=path.end();++it) {
2390 if ((hrz && it->second==i) || (!hrz && it->first==i))
2391 ++cnt;
2392 }
2393 if (cnt<2 && !is_exactly_zero(X[hrz?I:i][hrz?i:J])) {
2394 path.push_back(make_pair(hrz?I:i,hrz?i:J));
2395 ipairs fullpath(stepping_stone_path(path,X));
2396 if (!fullpath.empty())
2397 return fullpath;
2398 path.pop_back();
2399 }
2400 }
2401 return ipairs(0);
2402 }
2403
2404 /*
2405 * Implementation of MODI (modified ditribution) method. It handles degenerate
2406 * solutions if they appear during the process.
2407 */
modi(const matrice & P_orig,matrice & X)2408 void tprob::modi(const matrice &P_orig,matrice &X) {
2409 matrice P(P_orig);
2410 int m=X.size(),n=X.front()._VECTptr->size();
2411 vecteur u(m),v(n);
2412 if (M.type==_IDNT) {
2413 gen largest(0);
2414 for (int i=0;i<m;++i) {
2415 for (int j=0;j<n;++j) {
2416 if (is_greater(X[i][j],largest,ctx))
2417 largest=X[i][j];
2418 }
2419 }
2420 P=subst(P,M,100*largest,false,ctx);
2421 }
2422 for (int i=0;i<m;++i) {
2423 u[i]=i==0?gen(0):identificateur(" u"+print_INT_(++var_index));
2424 }
2425 for (int j=0;j<n;++j) {
2426 v[j]=identificateur(" v"+print_INT_(++var_index));
2427 }
2428 vecteur vars(mergevecteur(vecteur(u.begin()+1,u.end()),v));
2429 while (true) {
2430 vecteur eqv;
2431 for (int i=0;i<m;++i) {
2432 for (int j=0;j<n;++j) {
2433 if (!is_exactly_zero(X[i][j]))
2434 eqv.push_back(u[i]+v[j]-P[i][j]);
2435 }
2436 }
2437 vecteur sol(*_linsolve(makesequence(eqv,vars),ctx)._VECTptr);
2438 vecteur U(1,0),V(sol.begin()+m-1,sol.end());
2439 U=mergevecteur(U,vecteur(sol.begin(),sol.begin()+m-1));
2440 gen cmin(0);
2441 bool optimal=true;
2442 int I,J;
2443 for (int i=0;i<m;++i) {
2444 for (int j=0;j<n;++j) {
2445 if (is_exactly_zero(X[i][j])) {
2446 gen c(P[i][j]-U[i]-V[j]);
2447 if (is_strictly_greater(cmin,c,ctx)) {
2448 cmin=c;
2449 optimal=false;
2450 I=i;
2451 J=j;
2452 }
2453 }
2454 }
2455 }
2456 if (optimal)
2457 break;
2458 ipairs path;
2459 path.push_back(make_pair(I,J));
2460 path=stepping_stone_path(path,X);
2461 gen d(X[path.at(1).first][path.at(1).second]);
2462 for (ipairs::const_iterator it=path.begin()+3;it<path.end();it+=2) {
2463 d=min(d,X[it->first][it->second],ctx);
2464 }
2465 for (int i=0;i<int(path.size());++i) {
2466 gen &Xij=X[path.at(i).first]._VECTptr->at(path.at(i).second);
2467 gen x(Xij+(i%2?-d:d));
2468 bool has_zero=false;
2469 for (ipairs::const_iterator it=path.begin();it!=path.end();++it) {
2470 if (is_exactly_zero(X[it->first][it->second])) {
2471 has_zero=true;
2472 break;
2473 }
2474 }
2475 if ((!is_exactly_zero(x) && is_strictly_greater(gen(1)/gen(2),x,ctx)) ||
2476 (is_exactly_zero(x) && has_zero))
2477 x=eps;
2478 Xij=x;
2479 }
2480 }
2481 X=*exact(_epsilon2zero(_evalf(X,ctx),ctx),ctx)._VECTptr;
2482 }
2483
solve(const matrice & cost_matrix,matrice & sol)2484 void tprob::solve(const matrice &cost_matrix,matrice &sol) {
2485 north_west_corner(sol);
2486 modi(cost_matrix,sol);
2487 }
2488
2489 /*
2490 * END OF TPROB CLASS
2491 */
2492
2493 /*
2494 * Function 'tpsolve' solves a transportation problem using MODI method.
2495 *
2496 * Usage
2497 * ^^^^^
2498 * tpsolve(supply,demand,cost_matrix)
2499 *
2500 * Parameters
2501 * ^^^^^^^^^^
2502 * - supply : source capacity (vector of m positive integers)
2503 * - demand : destination demand (vector of n positive integers)
2504 * - cost_matrix : real matrix C=[c_ij] of type mXn where c_ij is cost of
2505 * transporting an unit from ith source to jth destination
2506 * (a nonnegative number)
2507 *
2508 * Supply and demand vectors should contain only positive integers. Cost matrix
2509 * must be consisted of nonnegative real numbers, which do not have to be
2510 * integers. There is a possibility of adding a certain symbol to cost matrix,
2511 * usually M, to indicate the "infinite cost", effectively forbidding the
2512 * transportation on a certain route. The notation of the symbol may be chosen
2513 * arbitrarily, but must be used consistently within a single problem.
2514 *
2515 * The return value is a sequence of total (minimal) cost and matrix X=[x_ij]
2516 * of type mXn where x_ij is equal to number of units which have to be shipped
2517 * from ith source to jth destination, for all i=1,2,..,m and j=1,2,..,n.
2518 *
2519 * This function uses 'north_west_corner' to determine initial feasible
2520 * solution and then applies MODI method to optimize it (function 'modi', which
2521 * uses 'stepping_stone_path'). Also, it is capable of handling degeneracy of
2522 * the initial solution and during iterations of MODI method.
2523 *
2524 * If the given problem is not balanced, i.e. if supply exceeds demand or vice
2525 * versa, dummy supply/demand points will be automatically added to the
2526 * problem, augmenting the cost matrix with zeros. Resulting matrix will not
2527 * contain dummy point.
2528 *
2529 * Examples
2530 * ^^^^^^^^
2531 * Balanced transportation problem:
2532 * tpsolve([12,17,11],[10,10,10,10],[[500,750,300,450],[650,800,400,600],[400,700,500,550]])
2533 * >> 2020,[[0,0,2,10],[0,9,8,0],[10,1,0,0]]
2534 * Non-balanced transportation problem:
2535 * tpsolve([7,10,8,8,9,6],[9,6,12,8,10],[[36,40,32,43,29],[28,27,29,40,38],[34,35,41,29,31],[41,42,35,27,36],[25,28,40,34,38],[31,30,43,38,40]])
2536 * >> [[0,0,2,0,5],[0,0,10,0,0],[0,0,0,0,5],[0,0,0,8,0],[9,0,0,0,0],[0,6,0,0,0]]
2537 * Transportation problem with forbidden routes:
2538 * tpsolve([95,70,165,165],[195,150,30,45,75],[[15,M,45,M,0],[12,40,M,M,0],[0,15,25,25,0],[M,0,M,12,0]])
2539 * >> [[20,0,0,0,75],[70,0,0,0,0],[105,0,30,30,0],[0,150,0,15,0]]
2540 * Assignment problem:
2541 * tpsolve([1,1,1,1],[1,1,1,1],[[10,12,9,11],[5,10,7,8],[12,14,13,11],[8,15,11,9]])
2542 * >> [[0,0,1,0],[1,0,0,0],[0,1,0,0],[0,0,0,1]]
2543 */
_tpsolve(const gen & g,GIAC_CONTEXT)2544 gen _tpsolve(const gen &g,GIAC_CONTEXT) {
2545 if (g.type==_STRNG && g.subtype==-1) return g;
2546 if (g.type!=_VECT || g.subtype!=_SEQ__VECT)
2547 return gentypeerr(contextptr);
2548 vecteur &gv=*g._VECTptr;
2549 if (gv.size()<3)
2550 return gensizeerr(contextptr);
2551 if (gv[0].type!=_VECT || gv[1].type!=_VECT ||
2552 gv[2].type!=_VECT || !ckmatrix(*gv[2]._VECTptr))
2553 return gentypeerr(contextptr);
2554 vecteur supply(*gv[0]._VECTptr),demand(*gv[1]._VECTptr);
2555 matrice P(*gv[2]._VECTptr);
2556 vecteur sy(*_lname(P,contextptr)._VECTptr);
2557 int m=supply.size(),n=demand.size();
2558 if (sy.size()>1 || m!=int(P.size()) || n!=int(P.front()._VECTptr->size()))
2559 return gensizeerr(contextptr);
2560 gen M(sy.size()==1 && sy[0].type==_IDNT?sy[0]:0);
2561 gen ts(_sum(supply,contextptr)),td(_sum(demand,contextptr));
2562 if (ts!=td) {
2563 *logptr(contextptr) << "Warning: transportation problem is not balanced\n";
2564 if (is_greater(ts,td,contextptr)) {
2565 demand.push_back(ts-td);
2566 P=mtran(P);
2567 P.push_back(vecteur(m,0));
2568 P=mtran(P);
2569 }
2570 else {
2571 supply.push_back(td-ts);
2572 P.push_back(vecteur(n,0));
2573 }
2574 }
2575 matrice X;
2576 tprob tp(supply,demand,M,contextptr);
2577 tp.solve(P,X);
2578 if (is_strictly_greater(ts,td,contextptr)) {
2579 X=mtran(X);
2580 X.pop_back();
2581 X=mtran(X);
2582 }
2583 else if (is_strictly_greater(td,ts,contextptr))
2584 X.pop_back();
2585 gen cost(0);
2586 for (int i=0;i<m;++i) {
2587 for (int j=0;j<n;++j) {
2588 cost+=P[i][j]*X[i][j];
2589 }
2590 }
2591 return makesequence(cost,X);
2592 }
2593 static const char _tpsolve_s []="tpsolve";
2594 static define_unary_function_eval (__tpsolve,&_tpsolve,_tpsolve_s);
2595 define_unary_function_ptr5(at_tpsolve,alias_at_tpsolve,&__tpsolve,0,true)
2596
compute_invdiff(int n,int k,vecteur & xv,vecteur & yv,map<tprob::ipair,gen> & invdiff,GIAC_CONTEXT)2597 gen compute_invdiff(int n,int k,vecteur &xv,vecteur &yv,map<tprob::ipair,gen> &invdiff,GIAC_CONTEXT) {
2598 tprob::ipair I=make_pair(n,k);
2599 assert(n<=k);
2600 gen res(invdiff[I]);
2601 if (!is_zero(res))
2602 return res;
2603 if (n==0)
2604 return invdiff[I]=yv[k];
2605 if (n==1)
2606 return invdiff[I]=(xv[k]-xv[0])/(yv[k]-yv[0]);
2607 gen d1(compute_invdiff(n-1,n-1,xv,yv,invdiff,contextptr));
2608 gen d2(compute_invdiff(n-1,k,xv,yv,invdiff,contextptr));
2609 return invdiff[I]=(xv[k]-xv[n-1])/(d2-d1);
2610 }
2611
thiele(int k,vecteur & xv,vecteur & yv,identificateur & var,map<tprob::ipair,gen> & invdiff,GIAC_CONTEXT)2612 gen thiele(int k,vecteur &xv,vecteur &yv,identificateur &var,map<tprob::ipair,gen> &invdiff,GIAC_CONTEXT) {
2613 if (k==int(xv.size()))
2614 return gen(0);
2615 gen phi(compute_invdiff(k,k,xv,yv,invdiff,contextptr));
2616 return (var-xv[k-1])/(phi+thiele(k+1,xv,yv,var,invdiff,contextptr));
2617 }
2618
2619 /*
2620 * 'thiele' computes rational interpolation for the given list of points using
2621 * Thiele's method with continued fractions.
2622 *
2623 * Source: http://www.astro.ufl.edu/~kallrath/files/pade.pdf (see page 19)
2624 *
2625 * Usage
2626 * ^^^^^
2627 * thiele(data,v)
2628 * or thiele(data_x,data_y,v)
2629 *
2630 * Parameters
2631 * ^^^^^^^^^^
2632 * - data : list of points [[x1,y1],[x2,y2],...,[xn,yn]]
2633 * - v : identifier or a symbolic expression
2634 * - data_x : list of x coordinates [x1,x2,...,xn]
2635 * - data_y : list of y coordinates [y1,y2,...,yn]
2636 *
2637 * The return value is an expression R(v), where R is rational interpolant of
2638 * the given set of points.
2639 *
2640 * Note that the interpolant may have singularities in
2641 * [min(data_x),max(data_x)].
2642 *
2643 * Example
2644 * ^^^^^^^
2645 * Function f(x)=(1-x^4)*exp(1-x^3) is sampled on interval [-1,2] in 13
2646 * equidistant points:
2647 *
2648 * data_x:=[-1,-0.75,-0.5,-0.25,0,0.25,0.5,0.75,1,1.25,1.5,1.75,2],
2649 * data_y:=[0.0,2.83341735599,2.88770329586,2.75030303645,2.71828182846,
2650 * 2.66568510781,2.24894558809,1.21863761951,0.0,-0.555711613283,
2651 * -0.377871362418,-0.107135851128,-0.0136782294833]
2652 *
2653 * To obtain rational function passing through these points, input:
2654 * thiele(data_x,data_y,x)
2655 * Output:
2656 * (-1.55286115659*x^6+5.87298387514*x^5-5.4439152812*x^4+1.68655817708*x^3
2657 * -2.40784868317*x^2-7.55954205222*x+9.40462512097)/(x^6-1.24295718965*x^5
2658 * -1.33526268624*x^4+4.03629272425*x^3-0.885419321*x^2-2.77913222418*x+3.45976823393)
2659 */
_thiele(const gen & g,GIAC_CONTEXT)2660 gen _thiele(const gen &g,GIAC_CONTEXT) {
2661 if (g.type==_STRNG && g.subtype==-1) return g;
2662 if (g.type!=_VECT || g.subtype!=_SEQ__VECT)
2663 return gentypeerr(contextptr);
2664 vecteur &gv=*g._VECTptr;
2665 if (gv.size()<2)
2666 return gensizeerr(contextptr);
2667 vecteur xv,yv;
2668 gen x;
2669 if (gv[0].type!=_VECT)
2670 return gentypeerr(contextptr);
2671 if (ckmatrix(gv[0])) {
2672 matrice m(mtran(*gv[0]._VECTptr));
2673 if (m.size()!=2)
2674 return gensizeerr(contextptr);
2675 xv=*m[0]._VECTptr;
2676 yv=*m[1]._VECTptr;
2677 x=gv[1];
2678 }
2679 else {
2680 if (gv[1].type!=_VECT)
2681 return gentypeerr(contextptr);
2682 if (gv[0]._VECTptr->size()!=gv[1]._VECTptr->size())
2683 return gensizeerr(contextptr);
2684 xv=*gv[0]._VECTptr;
2685 yv=*gv[1]._VECTptr;
2686 x=gv[2];
2687 }
2688 gen var(x.type==_IDNT?x:identificateur(" x"));
2689 map<tprob::ipair,gen> invdiff;
2690 gen rat(yv[0]+thiele(1,xv,yv,*var._IDNTptr,invdiff,contextptr));
2691 if (x.type==_IDNT) {
2692 // detect singularities
2693 gen den(_denom(rat,contextptr));
2694 matrice sing;
2695 if (*_lname(den,contextptr)._VECTptr==vecteur(1,x)) {
2696 for (int i=0;i<int(xv.size())-1;++i) {
2697 gen y1(_evalf(subst(den,x,xv[i],false,contextptr),contextptr));
2698 gen y2(_evalf(subst(den,x,xv[i+1],false,contextptr),contextptr));
2699 if (is_positive(-y1*y2,contextptr))
2700 sing.push_back(makevecteur(xv[i],xv[i+1]));
2701 }
2702 }
2703 if (!sing.empty()) {
2704 *logptr(contextptr) << "Warning, interpolant has singularities in ";
2705 for (int i=0;i<int(sing.size());++i) {
2706 *logptr(contextptr) << "(" << sing[i][0] << "," << sing[i][1] << ")";
2707 if (i<int(sing.size())-1)
2708 *logptr(contextptr) << (i<int(sing.size())-2?", ":" and ");
2709 }
2710 *logptr(contextptr) << "\n";
2711 }
2712 }
2713 else
2714 rat=simp(subst(rat,var,x,false,contextptr),contextptr);
2715 return simp(rat,contextptr);
2716 }
2717 static const char _thiele_s []="thiele";
2718 static define_unary_function_eval (__thiele,&_thiele,_thiele_s);
2719 define_unary_function_ptr5(at_thiele,alias_at_thiele,&__thiele,0,true)
2720
add_identifiers(const gen & source,vecteur & dest,GIAC_CONTEXT)2721 void add_identifiers(const gen &source,vecteur &dest,GIAC_CONTEXT) {
2722 vecteur v(*_lname(source,contextptr)._VECTptr);
2723 for (const_iterateur it=v.begin();it!=v.end();++it) {
2724 if (!contains(dest,*it))
2725 dest.push_back(*it);
2726 }
2727 dest=*_sort(dest,contextptr)._VECTptr;
2728 }
2729
indexof(const gen & g,const vecteur & v)2730 int indexof(const gen &g,const vecteur &v) {
2731 int n=v.size();
2732 for (int i=0;i<n;++i) {
2733 if (v.at(i)==g)
2734 return i;
2735 }
2736 return -1;
2737 }
2738
2739 /*
2740 * 'nlpsolve' computes an optimum of a nonlinear objective function, subject to
2741 * nonlinear equality and inequality constraints, using the COBYLA algorithm.
2742 *
2743 * Syntax
2744 * ^^^^^^
2745 * nlpsolve(objective, [constr], [bd], [opts])
2746 *
2747 * - objective: objective function
2748 * - constr: list of constraints
2749 * - bd: sequence of arguments of type x=a..b, where x is a problem variable
2750 * and a,b are reals
2751 * - opts: one of:
2752 * assume=nlp_nonnegative
2753 * maximize[=true]
2754 * nlp_initialpoint=[x1=a,x2=b,...]
2755 * nlp_precision=real
2756 * nlp_iterationlimit=intg
2757 *
2758 * If initial point is not given, it will be automatically generated. The given
2759 * point does not need to be feasible. Note that choosing a good initial point
2760 * is needed for obtaining a correct solution in some cases.
2761 *
2762 * Examples
2763 * ^^^^^^^^
2764 * (problems taken from:
2765 * www.ai7.uni-bayreuth.de/test_problem_coll.pdf
2766 * https://www.maplesoft.com/support/help/maple/view.aspx?path=Optimization%2FNLPSolve)
2767 *
2768 * nlpsolve(
2769 * (x1-10)^3+(x2-20)^3,
2770 * [(x1-5)^2+(x2-5)^2>=100,(x2-5)^2+(x1-6)^2<=82.81],
2771 * nlp_initialpoint=[x1=20.1,x2=5.84]) // problem 19, using initial point
2772 * nlpsolve(sin(x1+x2)+(x1-x2)^2-1.5x1+2.5x2+1,x1=-1.5..4,x2=-3..3) // problem 5
2773 * nlpsolve(ln(1+x1^2)-x2,[(1+x1^2)^2+x2^2=4]) // problem 7
2774 * nlpsolve(
2775 * x1,[x2>=exp(x1),x3>=exp(x2)],maximize=true,
2776 * x1=0..100,x2=0..100,x3=0..10,nlp_initialpoint=[x1=1,x2=1.05,x3=2.9]) // problem 34 (modified)
2777 * nlpsolve(-x1*x2*x3,[72-x1-2x2-2x3>=0],x1=0..20,x2=0..11,x3=0..42) // problem 36
2778 * nlpsolve(2-1/120*x1*x2*x3*x4*x5,[x1<=1,x2<=2,x3<=3,x4<=4,x5<=5],assume=nlp_nonnegative) // problem 45
2779 * nlpsolve(sin(x)/x,x=1..30) // Maple computes wrong result for this example, at least on their web page
2780 * nlpsolve(x^3+2x*y-2y^2,x=-10..10,y=-10..10,nlp_initialpoint=[x=3,y=4],maximize) // Maple example
2781 * nlpsolve(w^3*(v-w)^2+(w-x-1)^2+(x-y-2)^2+(y-z-3)^2,[w+x+y+z<=5,3z+2v=3],assume=nlp_nonnegative) // Maple example
2782 * nlpsolve(sin(x)*Psi(x),x=1..20,nlp_initialpoint=[x=16]) // Maple example, needs an initial point
2783 */
_nlpsolve(const gen & g,GIAC_CONTEXT)2784 gen _nlpsolve(const gen &g,GIAC_CONTEXT) {
2785 if (g.type==_STRNG && g.subtype==-1) return g;
2786 if (g.type!=_VECT || g.subtype!=_SEQ__VECT || g._VECTptr->size() < 2)
2787 return gentypeerr(contextptr);
2788 vecteur &gv=*g._VECTptr;
2789 vecteur constr,vars,initp;
2790 gen &obj=gv.front();
2791 add_identifiers(obj,vars,contextptr);
2792 const_iterateur it=gv.begin();
2793 bool maximize=false;
2794 int maxiter=RAND_MAX;
2795 double eps=epsilon(contextptr);
2796 if (gv.at(1).type==_VECT) {
2797 constr=*gv.at(1)._VECTptr;
2798 add_identifiers(constr,vars,contextptr);
2799 ++it;
2800 }
2801 initp=vecteur(vars.size(),gen(1));
2802 while (++it!=gv.end()) {
2803 if (*it==at_maximize || (it->is_integer() && it->val==_NLP_MAXIMIZE))
2804 maximize=true;
2805 else if (it->is_symb_of_sommet(at_equal)) {
2806 gen &lh=it->_SYMBptr->feuille._VECTptr->front();
2807 gen &rh=it->_SYMBptr->feuille._VECTptr->back();
2808 if (lh==at_assume && rh.is_integer() && rh.val==_NLP_NONNEGATIVE) {
2809 for (const_iterateur jt=vars.begin();jt!=vars.end();++jt) {
2810 constr.push_back(symbolic(at_inferieur_egal,makevecteur(gen(0),*jt)));
2811 }
2812 } else if (lh==at_maximize && rh.is_integer())
2813 maximize=(bool)rh.val;
2814 else if (lh.is_integer() && lh.val==_NLP_INITIALPOINT && rh.type==_VECT) {
2815 vecteur &pnt=*rh._VECTptr;
2816 for (const_iterateur jt=pnt.begin();jt!=pnt.end();++jt) {
2817 gen var;
2818 if (jt->is_symb_of_sommet(at_equal) &&
2819 contains(vars,var=jt->_SYMBptr->feuille._VECTptr->front()))
2820 initp.at(indexof(var,vars))=jt->_SYMBptr->feuille._VECTptr->back();
2821 }
2822 } else if (lh.is_integer() && lh.val==_NLP_ITERATIONLIMIT && rh.is_integer())
2823 maxiter=rh.val;
2824 else if (lh.is_integer() && lh.val==_NLP_MAXIMIZE && rh.is_integer())
2825 maximize=(bool)rh.val;
2826 else if(lh.is_integer() && lh.val==_NLP_PRECISION && rh.type==_DOUBLE_)
2827 eps=rh.DOUBLE_val();
2828 else if (contains(vars,lh) && rh.is_symb_of_sommet(at_interval)) {
2829 gen &lb=rh._SYMBptr->feuille._VECTptr->front();
2830 gen &ub=rh._SYMBptr->feuille._VECTptr->back();
2831 if (!is_inf(lh))
2832 constr.push_back(symbolic(at_superieur_egal,makevecteur(lh,lb)));
2833 if (!is_inf(rh))
2834 constr.push_back(symbolic(at_inferieur_egal,makevecteur(lh,ub)));
2835 }
2836 }
2837 }
2838 if (constr.empty()) {
2839 *logptr(contextptr) << "Error: no contraints detected\n";
2840 return gensizeerr(contextptr);
2841 }
2842 bool feasible=true;
2843 for (it=constr.begin();it!=constr.end();++it) {
2844 if (it->is_symb_of_sommet(at_equal)) {
2845 gen expr=_equal2diff(*it,contextptr);
2846 if (!is_zero(_subs(makesequence(expr,vars,initp),contextptr))) {
2847 feasible=false;
2848 break;
2849 }
2850 } else if (it->is_symb_of_sommet(at_inferieur_egal) || it->is_symb_of_sommet(at_superieur_egal)) {
2851 if (_evalb(_subs(makesequence(*it,vars,initp),contextptr),contextptr).val==0) {
2852 feasible=false;
2853 break;
2854 }
2855 } else {
2856 *logptr(contextptr) << "Error: unrecognized constraint " << *it << "\n";
2857 return gentypeerr(contextptr);
2858 }
2859 }
2860 gen sol,optval;
2861 try {
2862 if (!feasible) {
2863 gen tmpinitp=_fMin(makesequence(gen(0),constr,vars,initp),contextptr);
2864 if (tmpinitp.type==_VECT)
2865 initp=*tmpinitp._VECTptr;
2866 else initp.clear();
2867 if (is_undef(initp) || initp.empty()) {
2868 *logptr(contextptr) << "Error: unable to generate a feasible initial point\n";
2869 return undef;
2870 }
2871 //*logptr(contextptr) << "Using a generated feasible initial point " << initp << "\n";
2872 }
2873 gen args=makesequence(obj,constr,vars,initp,gen(eps),gen(maxiter));
2874 if (maximize)
2875 sol=_fMax(args,contextptr);
2876 else
2877 sol=_fMin(args,contextptr);
2878 } catch (std::runtime_error &err) {
2879 *logptr(contextptr) << "Error: " << err.what() << "\n";
2880 return undef;
2881 }
2882 if (is_undef(sol))
2883 return undef;
2884 optval=_subs(makesequence(obj,vars,sol),contextptr);
2885 return gen(makevecteur(optval,_zip(makesequence(at_equal,vars,sol),contextptr)),_LIST__VECT);
2886 }
2887 static const char _nlpsolve_s []="nlpsolve";
2888 static define_unary_function_eval (__nlpsolve,&_nlpsolve,_nlpsolve_s);
2889 define_unary_function_ptr5(at_nlpsolve,alias_at_nlpsolve,&__nlpsolve,0,true)
2890
2891 /*
2892 * Returns the trigonometric polynomial in variable x passing through points
2893 * with ordinate componets in 'data' and the abscissa components equally spaced between
2894 * a and b (the first being equal a and the last being equal to b).
2895 */
triginterp(const vecteur & data,const gen & a,const gen & b,const identificateur & x,GIAC_CONTEXT)2896 gen triginterp(const vecteur &data,const gen &a,const gen &b,const identificateur &x,GIAC_CONTEXT) {
2897 int n=data.size();
2898 if (n<2)
2899 return gensizeerr(contextptr);
2900 int N=(n%2)==0?n/2:(n-1)/2;
2901 gen T=(b-a)*fraction(n,n-1),twopi=2*_IDNT_pi(),X;
2902 matrice cos_coeff=*_matrix(makesequence(N,n,0),contextptr)._VECTptr;
2903 matrice sin_coeff=*_matrix(makesequence(N,n,0),contextptr)._VECTptr;
2904 for (int k=0;k<n;++k) {
2905 X=twopi*(a/T+fraction(k,n));
2906 for (int j=1;j<=N;++j) {
2907 cos_coeff[j-1]._VECTptr->at(k)=cos(j*X,contextptr);
2908 sin_coeff[j-1]._VECTptr->at(k)=sin(j*X,contextptr);
2909 }
2910 }
2911 gen tp=_mean(data,contextptr);
2912 for (int j=0;j<N;++j) {
2913 gen c=fraction(((n%2)==0 && j==N-1)?1:2,n);
2914 gen ak=_evalc(trig2exp(scalarproduct(data,*cos_coeff[j]._VECTptr,contextptr),contextptr),contextptr);
2915 gen bk=_evalc(trig2exp(scalarproduct(data,*sin_coeff[j]._VECTptr,contextptr),contextptr),contextptr);
2916 tp+=simp(c*ak,contextptr)*cos(simp((j+1)*twopi/T,contextptr)*x,contextptr);
2917 tp+=simp(c*bk,contextptr)*sin(simp((j+1)*twopi/T,contextptr)*x,contextptr);
2918 }
2919 return tp;
2920 }
2921
_triginterp(const gen & g,GIAC_CONTEXT)2922 gen _triginterp(const gen &g,GIAC_CONTEXT) {
2923 if (g.type==_STRNG && g.subtype==-1) return g;
2924 if (g.type!=_VECT || g.subtype!=_SEQ__VECT)
2925 return gentypeerr(contextptr);
2926 vecteur &args=*g._VECTptr;
2927 if (args.size()<2)
2928 return gensizeerr(contextptr);
2929 if (args.front().type!=_VECT)
2930 return gentypeerr(contextptr);
2931 vecteur &data=*args.front()._VECTptr;
2932 gen x,ab,a,b,&vararg=args.at(1);
2933 if (vararg.is_symb_of_sommet(at_equal) &&
2934 (x=_lhs(vararg,contextptr)).type==_IDNT &&
2935 (ab=_rhs(vararg,contextptr)).is_symb_of_sommet(at_interval)) {
2936 a=_lhs(ab,contextptr);
2937 b=_rhs(ab,contextptr);
2938 } else if (args.size()==4 && (x=args.back()).type==_IDNT) {
2939 a=args.at(1);
2940 b=args.at(2);
2941 } else return gensizeerr(contextptr);
2942 gen tp=triginterp(data,a,b,*x._IDNTptr,contextptr);
2943 if (is_approx(data) || is_approx(a) || is_approx(b))
2944 tp=_evalf(tp,contextptr);
2945 return tp;
2946 }
2947 static const char _triginterp_s []="triginterp";
2948 static define_unary_function_eval (__triginterp,&_triginterp,_triginterp_s);
2949 define_unary_function_ptr5(at_triginterp,alias_at_triginterp,&__triginterp,0,true)
2950
2951 /* select a good bandwidth for kernel density estimation using a direct plug-in method (DPI),
2952 * Gaussian kernel is assumed */
select_bandwidth_dpi(const vector<double> & data,double sd)2953 double select_bandwidth_dpi(const vector<double> &data,double sd) {
2954 int n=data.size();
2955 double g6=1.23044723*sd,s=0,t,t2;
2956 for (vector<double>::const_iterator it=data.begin();it!=data.end();++it) {
2957 for (vector<double>::const_iterator jt=it+1;jt!=data.end();++jt) {
2958 t=(*it-*jt)/g6;
2959 t2=t*t;
2960 s+=(2*t2*(t2*(t2-15)+45)-30)*std::exp(-t2/2);
2961 }
2962 }
2963 s-=15*n;
2964 double g4=g6*std::pow(-(6.0*n)/s,1/7.0);
2965 s=0;
2966 for (vector<double>::const_iterator it=data.begin();it!=data.end();++it) {
2967 for (vector<double>::const_iterator jt=it+1;jt!=data.end();++jt) {
2968 t=(*it-*jt)/g4;
2969 t2=t*t;
2970 s+=(2*t2*(t2-6)+6)*std::exp(-t2/2);
2971 }
2972 }
2973 s+=3*n;
2974 return std::pow(double(n)/(M_SQRT2*s),0.2)*g4;
2975 }
2976
fft_sum(const vecteur & c,const vecteur & k,int M,GIAC_CONTEXT)2977 gen fft_sum(const vecteur &c,const vecteur &k,int M,GIAC_CONTEXT) {
2978 return _scalar_product(makesequence(c,_mid(makesequence(_convolution(makesequence(c,k),contextptr),M,M),contextptr)),contextptr);
2979 }
2980
2981 /* faster bandwidth DPI selector using binned data and FFT */
select_bandwidth_dpi_bins(int n,const vecteur & c,double d,double sd,GIAC_CONTEXT)2982 double select_bandwidth_dpi_bins(int n,const vecteur &c,double d,double sd,GIAC_CONTEXT) {
2983 int M=c.size();
2984 vecteur k(2*M+1);
2985 double g6=1.23044723*sd,s=0,t,t2;
2986 for (int i=0;i<=2*M;++i) {
2987 t=d*double(i-M)/g6;
2988 t2=t*t;
2989 k[i]=gen((2*t2*(t2*(t2-15)+45)-30)*std::exp(-t2/2));
2990 }
2991 s=_evalf(fft_sum(c,k,M,contextptr),contextptr).DOUBLE_val();
2992 double g4=g6*std::pow(-(6.0*n)/s,1/7.0);
2993 for (int i=0;i<=2*M;++i) {
2994 t=d*double(i-M)/g4;
2995 t2=t*t;
2996 k[i]=gen((2*t2*(t2-6)+6)*std::exp(-t2/2));
2997 }
2998 s=_evalf(fft_sum(c,k,M,contextptr),contextptr).DOUBLE_val();
2999 return std::pow(double(n)/(M_SQRT2*s),0.2)*g4;
3000 }
3001
3002 /* kernel density estimation with Gaussian kernel */
kernel_density(const vector<double> & data,double bw,double sd,int bins,double a,double b,int interp,const gen & x,GIAC_CONTEXT)3003 gen kernel_density(const vector<double> &data,double bw,double sd,int bins,double a,double b,int interp,const gen &x,GIAC_CONTEXT) {
3004 int n=data.size();
3005 double SQRT_2PI=std::sqrt(2.0*M_PI);
3006 if (bins<=0) { // return density as a sum of exponential functions, usable for up to few hundred samples
3007 if (bw<=0)
3008 bw=select_bandwidth_dpi(data,sd);
3009 double fac=bw*n*SQRT_2PI;
3010 gen res(0),h(2.0*bw*bw);
3011 for (vector<double>::const_iterator it=data.begin();it!=data.end();++it) {
3012 res+=exp(-pow(x-gen(*it),2)/h,contextptr);
3013 }
3014 return res/gen(fac);
3015 }
3016 /* FFT method, constructs an approximation on [a,b] with the specified number of bins.
3017 * If interp>0, interpolation of order interp is performed and the density is returned piecewise. */
3018 assert(b>a && bins>0);
3019 double d=(b-a)/(bins-1);
3020 vecteur c(bins,0);
3021 int index;
3022 for (vector<double>::const_iterator it=data.begin();it!=data.end();++it) {
3023 index=(int)((*it-a)/d+0.5);
3024 if (index>=0 && index<bins) c[index]+=1;
3025 }
3026 if (bw<=0) { // select bandwidth
3027 if (n<=1000)
3028 bw=select_bandwidth_dpi(data,sd);
3029 else bw=select_bandwidth_dpi_bins(n,c,d,sd,contextptr);
3030 *logptr(contextptr) << "selected bandwidth: " << bw << "\n";
3031 }
3032 int L=std::min(bins-1,(int)std::floor(1+4*bw/d));
3033 vecteur k(2*L+1);
3034 for (int i=0;i<=2*L;++i) {
3035 k[i]=gen(1.0/(n*bw*SQRT_2PI)*std::exp(-std::pow(d*double(i-L)/bw,2)/2.0));
3036 }
3037 gen res=_mid(makesequence(_convolution(makesequence(c,k),contextptr),L,bins),contextptr);
3038 if (interp>0) { // interpolate the obtained points
3039 int pos0=0;
3040 if (x.type!=_IDNT) {
3041 double xd=_evalf(x,contextptr).DOUBLE_val();
3042 if (xd<a || xd>=b || (pos0=std::floor((xd-a)/d))>bins-2)
3043 return 0;
3044 if (interp==1) {
3045 gen &y1=res._VECTptr->at(pos0),&y2=res._VECTptr->at(pos0+1),x1=a+pos0*d;
3046 return y1+(x-x1)*(y2-y1)/gen(d);
3047 }
3048 }
3049 vecteur pos(bins);
3050 for (int i=0;i<bins;++i) pos[i]=a+d*i;
3051 identificateur X=x.type==_IDNT?*x._IDNTptr:identificateur(" X");
3052 vecteur p=*_spline(makesequence(pos,res,X,interp),contextptr)._VECTptr;
3053 vecteur args(0);
3054 if (x.type==_IDNT)
3055 args.reserve(2*bins+1);
3056 for (int i=0;i<bins;++i) {
3057 if (x.type==_IDNT) {
3058 args.push_back(i+1<bins?symb_inferieur_strict(X,pos[i]):symb_inferieur_egal(X,pos[i]));
3059 args.push_back(i==0?gen(0):p[i-1]);
3060 } else if (i==pos0) res=simp(subst(p[i],X,x,false,contextptr),contextptr);
3061 if (i+1<bins && interp>1 && !_solve(makesequence(p[i],symb_equal(X,symb_interval(pos[i],pos[i+1]))),contextptr)._VECTptr->empty())
3062 *logptr(contextptr) << "Warning: interpolated density has negative values in ["
3063 << pos[i] << "," << pos[i+1] << "]\n";
3064 }
3065 if (x.type!=_IDNT) return res;
3066 args.push_back(0);
3067 res=symbolic(at_piecewise,change_subtype(args,_SEQ__VECT));
3068 return res;
3069 }
3070 return res;
3071 }
3072
parse_interval(const gen & feu,double & a,double & b,GIAC_CONTEXT)3073 bool parse_interval(const gen &feu,double &a,double &b,GIAC_CONTEXT) {
3074 vecteur &v=*feu._VECTptr;
3075 gen l=v.front(),r=v.back();
3076 if ((l=_evalf(l,contextptr)).type!=_DOUBLE_ || (r=_evalf(r,contextptr)).type!=_DOUBLE_ ||
3077 !is_strictly_greater(r,l,contextptr))
3078 return false;
3079 a=l.DOUBLE_val(); b=r.DOUBLE_val();
3080 return true;
3081 }
3082
_kernel_density(const gen & g,GIAC_CONTEXT)3083 gen _kernel_density(const gen &g,GIAC_CONTEXT) {
3084 if (g.type==_STRNG && g.subtype==-1) return g;
3085 if (g.type!=_VECT)
3086 return gentypeerr(contextptr);
3087 gen x=identificateur("x");
3088 double a=0,b=0,bw=0,sd,d,sx=0,sxsq=0;
3089 int bins=100,interp=1,method=_KDE_METHOD_LIST,bw_method=_KDE_BW_METHOD_DPI;
3090 if (g.subtype==_SEQ__VECT) {
3091 // parse options
3092 for (const_iterateur it=g._VECTptr->begin()+1;it!=g._VECTptr->end();++it) {
3093 if (it->is_symb_of_sommet(at_equal)) {
3094 gen &opt=it->_SYMBptr->feuille._VECTptr->front();
3095 gen &v=it->_SYMBptr->feuille._VECTptr->back();
3096 if (opt==_KDE_BANDWIDTH) {
3097 if (v==at_select)
3098 bw_method=_KDE_BW_METHOD_DPI;
3099 else if (v==at_gauss || v==at_normal || v==at_normald)
3100 bw_method=_KDE_BW_METHOD_ROT;
3101 else {
3102 gen ev=_evalf(v,contextptr);
3103 if (ev.type!=_DOUBLE_ || !is_strictly_positive(ev,contextptr))
3104 return gensizeerr(contextptr);
3105 bw=ev.DOUBLE_val();
3106 }
3107 } else if (opt==_KDE_BINS) {
3108 if (!v.is_integer() || !is_strictly_positive(v,contextptr))
3109 return gensizeerr(contextptr);
3110 bins=v.val;
3111 } else if (opt==at_range) {
3112 if (v.type==_VECT) {
3113 if (v._VECTptr->size()!=2 || !parse_interval(v,a,b,contextptr))
3114 return gensizeerr(contextptr);
3115 } else if (!v.is_symb_of_sommet(at_interval) ||
3116 !parse_interval(v._SYMBptr->feuille,a,b,contextptr))
3117 return gensizeerr(contextptr);
3118 } else if (opt==at_output || opt==at_Output) {
3119 if (v==at_exact)
3120 method=_KDE_METHOD_EXACT;
3121 else if (v==at_piecewise)
3122 method=_KDE_METHOD_PIECEWISE;
3123 else if (v==_MAPLE_LIST)
3124 method=_KDE_METHOD_LIST;
3125 else return gensizeerr(contextptr);
3126 } else if (opt==at_interp) {
3127 if (!v.is_integer() || (interp=v.val)<1)
3128 return gensizeerr(contextptr);
3129 } else if (opt==at_spline) {
3130 if (!v.is_integer() || (interp=v.val)<1)
3131 return gensizeerr(contextptr);
3132 method=_KDE_METHOD_PIECEWISE;
3133 } else if (opt.type==_IDNT) {
3134 x=opt;
3135 if (!v.is_symb_of_sommet(at_interval) || !parse_interval(v._SYMBptr->feuille,a,b,contextptr))
3136 return gensizeerr(contextptr);
3137 } else if (opt==at_eval) x=v;
3138 else return gensizeerr(contextptr);
3139 } else if (it->type==_IDNT) x=*it;
3140 else if (it->is_symb_of_sommet(at_interval)) {
3141 if (!parse_interval(it->_SYMBptr->feuille,a,b,contextptr))
3142 return gensizeerr(contextptr);
3143 } else if (*it==at_exact)
3144 method=_KDE_METHOD_EXACT;
3145 else if (*it==at_piecewise)
3146 method=_KDE_METHOD_PIECEWISE;
3147 else return gensizeerr(contextptr);
3148 }
3149 }
3150 if (x.type!=_IDNT && (_evalf(x,contextptr).type!=_DOUBLE_ || method==_KDE_METHOD_LIST))
3151 return gensizeerr(contextptr);
3152 vecteur &data=g.subtype==_SEQ__VECT?*g._VECTptr->front()._VECTptr:*g._VECTptr;
3153 int n=data.size();
3154 if (n<2)
3155 return gensizeerr(contextptr);
3156 vector<double> ddata(n);
3157 gen e;
3158 for (const_iterateur it=data.begin();it!=data.end();++it) {
3159 if ((e=_evalf(*it,contextptr)).type!=_DOUBLE_)
3160 return gensizeerr(contextptr);
3161 d=ddata[it-data.begin()]=e.DOUBLE_val();
3162 sx+=d;
3163 sxsq+=d*d;
3164 }
3165 sd=std::sqrt(1/double(n-1)*(sxsq-1/double(n)*sx*sx));
3166 if (bw_method==_KDE_BW_METHOD_ROT) { // Silverman's rule of thumb
3167 double iqr=_evalf(_quartile3(data,contextptr)-_quartile1(data,contextptr),contextptr).DOUBLE_val();
3168 bw=1.06*std::min(sd,iqr/1.34)*std::pow(double(data.size()),-0.2);
3169 *logptr(contextptr) << "selected bandwidth: " << bw << "\n";
3170 }
3171 if (bins>0 && a==0 && b==0) {
3172 a=_evalf(_min(data,contextptr),contextptr).DOUBLE_val()-3*bw;
3173 b=_evalf(_max(data,contextptr),contextptr).DOUBLE_val()+3*bw;
3174 }
3175 if (method==_KDE_METHOD_EXACT)
3176 bins=0;
3177 else if (method==_KDE_METHOD_LIST) {
3178 if (bins<1)
3179 return gensizeerr(contextptr);
3180 interp=0;
3181 } else if (method==_KDE_METHOD_PIECEWISE) {
3182 if (bins<1 || interp<1)
3183 return gensizeerr(contextptr);
3184 }
3185 return kernel_density(ddata,bw,sd,bins,a,b,interp,x,contextptr);
3186 }
3187 static const char _kernel_density_s []="kernel_density";
3188 static define_unary_function_eval (__kernel_density,&_kernel_density,_kernel_density_s);
3189 define_unary_function_ptr5(at_kernel_density,alias_at_kernel_density,&__kernel_density,0,true)
3190
3191 static const char _kde_s []="kde";
3192 static define_unary_function_eval (__kde,&_kernel_density,_kde_s);
3193 define_unary_function_ptr5(at_kde,alias_at_kde,&__kde,0,true)
3194
3195 /* maximum likelihood estimation for Weibull distribution */
weibull_mle(const vecteur & S,const gen & k0,const gen & eps,GIAC_CONTEXT)3196 gen weibull_mle(const vecteur &S,const gen &k0,const gen &eps,GIAC_CONTEXT) {
3197 gen s0(0),s1(0),s2(0),L(0),t,l,n(S.size());
3198 for (const_iterateur it=S.begin();it!=S.end();++it) {
3199 L+=(l=ln(*it,contextptr));
3200 t=exp(k0*l,contextptr);
3201 s0+=t; s1+=t*l; s2+=t*sq(l);
3202 }
3203 L=L/n;
3204 gen ik0=_inv(k0,contextptr),k=k0-(s1-s0*(ik0+L))/(s2+s0*sq(ik0)-s1*(ik0+L));
3205 if (!is_strictly_positive(k,contextptr))
3206 return undef;
3207 if (is_greater(_abs(k-k0,contextptr),eps,contextptr))
3208 return weibull_mle(S,k,eps,contextptr);
3209 gen lambda=exp(_inv(k,contextptr)*ln(s0/n,contextptr),contextptr);
3210 return symbolic(at_weibulld,makesequence(k,lambda));
3211 }
3212
3213 /* maximum likelihood estimation for Cauchy distribution */
cauchy_mle(const vecteur & S,const gen & x0_init,const gen & gama_init,const gen & eps,GIAC_CONTEXT)3214 gen cauchy_mle(const vecteur &S,const gen &x0_init,const gen &gama_init,const gen &eps,GIAC_CONTEXT) {
3215 matrice m=*_matrix(makesequence(2,2,0),contextptr)._VECTptr;
3216 gen n(S.size()),d,dsq,den,densq,gsq=sq(gama_init),F(0),G(-n/2);
3217 gen &Fx0=m.front()._VECTptr->front(),&Fgama=m.front()._VECTptr->back(),
3218 &Gx0=m.back()._VECTptr->front(),&Ggama=m.back()._VECTptr->back();
3219 for (const_iterateur it=S.begin();it!=S.end();++it) {
3220 d=*it-x0_init; dsq=sq(d); den=gsq+dsq; densq=sq(den);
3221 F+=d/den; G+=gsq/den;
3222 Fx0+=(dsq-gsq)/densq; Fgama+=d/densq; Ggama+=dsq/densq;
3223 }
3224 Gx0=2*gsq*Fgama;
3225 Fgama=-2*gama_init*Fgama;
3226 Ggama=2*gama_init*Ggama;
3227 gen tmpdelta=_linsolve(makesequence(m,makevecteur(-F,-G)),contextptr);
3228 if (tmpdelta.type==_VECT) {
3229 vecteur &delta=*tmpdelta._VECTptr;
3230 if (delta.empty())
3231 return undef;
3232 gen x0=x0_init+delta.front(),gama=gama_init+delta.back();
3233 if (is_greater(_l2norm(delta,contextptr),eps,contextptr))
3234 return cauchy_mle(S,x0,gama,eps,contextptr);
3235 return symbolic(at_cauchyd,makesequence(x0,gama));
3236 }
3237 return undef;
3238 }
3239
3240 /* fit distribution of the given type to the given data using the method of maximum likelihood */
_fitdistr(const gen & g,GIAC_CONTEXT)3241 gen _fitdistr(const gen &g,GIAC_CONTEXT) {
3242 if (g.type==_STRNG && g.subtype==-1) return g;
3243 if (g.type!=_VECT || g.subtype!=_SEQ__VECT)
3244 return gentypeerr(contextptr);
3245 if (g._VECTptr->size()!=2)
3246 return gensizeerr(contextptr);
3247 if (g._VECTptr->front().type!=_VECT)
3248 return gensizeerr(contextptr);
3249 vecteur &S=*g._VECTptr->front()._VECTptr;
3250 int n=S.size(); // number of samples
3251 if (n<2)
3252 return gensizeerr(contextptr);
3253 gen &dist=g._VECTptr->back(),N(n);
3254 /* compute sample mean and variance */
3255 gen mean(0),var(0),ev;
3256 for (const_iterateur it=S.begin();it!=S.end();++it) {
3257 ev=_evalf(*it,contextptr);
3258 if (ev.type!=_DOUBLE_)
3259 return gensizeerr(contextptr);
3260 mean+=ev;
3261 }
3262 mean=mean/(N-1);
3263 for (const_iterateur it=S.begin();it!=S.end();++it) {
3264 var+=sq(*it-mean);
3265 }
3266 var=var/(N-1);
3267 gen sdev=sqrt(var,contextptr);
3268 /* fit the parameters of the specified distribution by the method of maximum likelihood */
3269 if (dist==at_normal || dist==at_normald || dist==at_NORMALD) {
3270 return symbolic(at_normald,makesequence(mean,sdev));
3271 } else if (dist==at_poisson || dist==at_POISSON) {
3272 for (const_iterateur it=S.begin();it!=S.end();++it) {
3273 if (!it->is_integer() || it->val<0) return gensizeerr(contextptr);
3274 }
3275 return symbolic(at_poisson,mean);
3276 } else if (dist==at_exp || dist==at_EXP || dist==at_exponential || dist==at_exponentiald) {
3277 for (const_iterateur it=S.begin();it!=S.end();++it) {
3278 if (!is_strictly_positive(*it,contextptr)) return gensizeerr(contextptr);
3279 }
3280 return symbolic(at_exponentiald,_inv(mean,contextptr));
3281 } else if (dist==at_geometric) {
3282 for (const_iterateur it=S.begin();it!=S.end();++it) {
3283 if (!it->is_integer() || it->val<1) return gensizeerr(contextptr);
3284 }
3285 return symbolic(at_geometric,_inv(mean,contextptr));
3286 } else if (dist==at_gammad || dist==at_Gamma) {
3287 gen slog(0);
3288 for (const_iterateur it=S.begin();it!=S.end();++it) {
3289 if (!is_strictly_positive(*it,contextptr)) return gensizeerr(contextptr);
3290 slog+=ln(*it,contextptr);
3291 }
3292 gen a_init=sq(mean)/var,aidn=identificateur(" a");
3293 gen e=ln(aidn,contextptr)-Psi(aidn,contextptr)-ln(mean,contextptr)+slog/N;
3294 gen a=_solve(makesequence(e,symb_equal(aidn,a_init),_NEWTON_SOLVER),contextptr);
3295 return symbolic(at_gammad,makesequence(a,a/mean));
3296 } else if (dist==at_betad || dist==at_Beta) {
3297 gen slog(0),s1log(0),aidn=identificateur(" a"),bidn=identificateur(" b");
3298 for (const_iterateur it=S.begin();it!=S.end();++it) {
3299 if (!is_greater(*it,0,contextptr) || is_strictly_greater(*it,1,contextptr))
3300 return gensizeerr(contextptr);
3301 slog+=ln(*it,contextptr); s1log+=ln(1-*it,contextptr);
3302 }
3303 gen fac=(var+sq(mean)-mean)/var,a_init=-mean*fac,b_init=(mean-1)*fac;
3304 gen e1=Psi(aidn,contextptr)-Psi(aidn+bidn,contextptr)-slog/N;
3305 gen e2=Psi(bidn,contextptr)-Psi(aidn+bidn,contextptr)-s1log/N;
3306 gen tmpsol=_fsolve(makesequence(makevecteur(e1,e2),makevecteur(aidn,bidn),makevecteur(a_init,b_init),_NEWTONJ_SOLVER),contextptr);
3307 if (tmpsol.type!=_VECT || tmpsol._VECTptr->size()!=2)
3308 return gensizeerr(contextptr);
3309 vecteur &sol=*tmpsol._VECTptr;
3310 return symbolic(at_betad,change_subtype(sol,_SEQ__VECT));
3311 } else if (dist==at_cauchy || dist==at_cauchyd) {
3312 gen x0_init=_median(S,contextptr);
3313 gen gama_init=(_quartile3(S,contextptr)-_quartile1(S,contextptr))/2;
3314 return cauchy_mle(S,x0_init,gama_init,1e-5,contextptr);
3315 } else if (dist==at_weibull || dist==at_weibulld) {
3316 if (is_zero(var)) return gensizeerr(contextptr);
3317 gen kidn=identificateur(" k"),slog(0);
3318 for (const_iterateur it=S.begin();it!=S.end();++it) {
3319 if (!is_positive(*it,contextptr)) return gensizeerr(contextptr);
3320 slog+=ln(*it,contextptr);
3321 }
3322 gen e=_Gamma(1+gen(2)/kidn,contextptr)/_Gamma(1+gen(1)/kidn,contextptr)-var/sq(mean)-1;
3323 gen k_init=_fsolve(makesequence(e,symb_equal(kidn,max(1,_inv(var,contextptr),contextptr)),_NEWTON_SOLVER),contextptr);
3324 return weibull_mle(S,k_init,1e-5,contextptr);
3325 }
3326 return gensizeerr(contextptr); // the distribution is not recognized
3327 }
3328 static const char _fitdistr_s []="fitdistr";
3329 static define_unary_function_eval (__fitdistr,&_fitdistr,_fitdistr_s);
3330 define_unary_function_ptr5(at_fitdistr,alias_at_fitdistr,&__fitdistr,0,true)
3331
3332 /* evaluates the function f(x,y,y') in the specified point and returns the result as a floating point number. */
eval_func(const gen & f,const vecteur & vars,const gen & x,const gen & y,const gen & dy,bool & errflag,GIAC_CONTEXT)3333 double eval_func(const gen &f,const vecteur &vars,const gen &x,const gen &y,const gen &dy,bool &errflag,GIAC_CONTEXT) {
3334 gen e=_evalf(subst(f,vars,makevecteur(x,y,dy),false,contextptr),contextptr);
3335 if (e.type!=_DOUBLE_) {
3336 errflag=false;
3337 return 0;
3338 }
3339 return e.DOUBLE_val();
3340 }
3341
3342 /* approximate the solution of the following boundary-value problem:
3343 * y''=f(x,y,y'), a<=x<=b, y(a)=alpha, y(b)=beta.
3344 * The solution is stored in the lists X, w1 and w2 such that
3345 * X=[x0=a,x1,x2,..,xN=b] and w1[k]=y(xk), w2[k]=y'(xk) for k=0,1,...,N.
3346 * Return value: 0 on success, 1 if the maximum number of iterations M is exceeded and 2 on computation failure.
3347 */
shooting(const gen & f,const gen & x_idn,const gen & y_idn,const gen & dy_idn,const gen & TK_orig,const gen & x1,const gen & x2,const gen & y1,const gen & y2,int N,double tol,int M,vecteur & X,vecteur & Y,vecteur & dY,GIAC_CONTEXT)3348 int shooting(const gen &f,const gen &x_idn,const gen &y_idn,const gen &dy_idn,const gen &TK_orig,
3349 const gen &x1,const gen &x2,const gen &y1,const gen &y2,
3350 int N,double tol,int M,vecteur &X,vecteur &Y,vecteur &dY,GIAC_CONTEXT) {
3351 gen dfy=derive(f,y_idn,contextptr),dfdy=derive(f,dy_idn,contextptr);
3352 vecteur vars=makevecteur(x_idn,y_idn,dy_idn);
3353 double a=x1.DOUBLE_val(),b=x2.DOUBLE_val(),alpha=y1.DOUBLE_val(),beta=y2.DOUBLE_val();
3354 double h=(b-a)/N,x,k11,k12,k21,k22,k31,k32,k41,k42,dk11,dk12,dk21,dk22,dk31,dk32,dk41,dk42,u1,u2,fv1,fv2,w1i,w2i;
3355 double TK=is_undef(TK_orig)?(beta-alpha)/(b-a):TK_orig.DOUBLE_val();
3356 vector<double> w1(N+1),w2(N+1);
3357 int k=1;
3358 dfy=simp(dfy,contextptr);
3359 dfdy=simp(dfdy,contextptr);
3360 bool ef=true;
3361 while (k<=M) {
3362 w1[0]=alpha;
3363 w2[0]=TK;
3364 u1=0;
3365 u2=1;
3366 for (int i=0;i<N;++i) {
3367 x=a+i*h;
3368 w1i=w1[i];
3369 w2i=w2[i];
3370 k11=h*w2i;
3371 k12=h*eval_func(f,vars,x,w1i,w2i,ef,contextptr);
3372 k21=h*(w2i+k12/2);
3373 k22=h*eval_func(f,vars,x+h/2,w1i+k11/2,w2i+k12/2,ef,contextptr);
3374 k31=h*(w2i+k22/2);
3375 k32=h*eval_func(f,vars,x+h/2,w1i+k21/2,w2i+k22/2,ef,contextptr);
3376 k41=h*(w2i+k32);
3377 k42=h*eval_func(f,vars,x+h,w1i+k31,w2i+k32,ef,contextptr);
3378 w1[i+1]=w1i+(k11+2*k21+2*k31+k41)/6;
3379 w2[i+1]=w2i+(k12+2*k22+2*k32+k42)/6;
3380 dk11=h*u2;
3381 dk12=h*(eval_func(dfy,vars,x,w1i,w2i,ef,contextptr)*u1+
3382 eval_func(dfdy,vars,x,w1i,w2i,ef,contextptr)*u2);
3383 dk21=h*(u2+dk12/2);
3384 fv1=eval_func(dfy,vars,x+h/2,w1i,w2i,ef,contextptr);
3385 fv2=eval_func(dfdy,vars,x+h/2,w1i,w2i,ef,contextptr);
3386 dk22=h*(fv1*(u1+dk11/2)+fv2*(u2+dk12/2));
3387 dk31=h*(u2+dk22/2);
3388 dk32=h*(fv1*(u1+dk21/2)+fv2*(u2+dk22/2));
3389 dk41=h*(u2+dk32);
3390 dk42=h*(eval_func(dfy,vars,x+h,w1i,w2i,ef,contextptr)*(u1+dk31)+
3391 eval_func(dfdy,vars,x+h,w1i,w2i,ef,contextptr)*(u2+dk32));
3392 u1+=(dk11+2*dk21+2*dk31+dk41)/6;
3393 u2+=(dk12+2*dk22+2*dk32+dk42)/6;
3394 if (!ef) return 2;
3395 }
3396 if (std::abs(w1[N]-beta)<=tol) {
3397 X.resize(N+1);
3398 Y.resize(N+1);
3399 dY.resize(N+1);
3400 for (int i=0;i<=N;++i) {
3401 X[i]=a+i*h;
3402 Y[i]=w1[i];
3403 dY[i]=w2[i];
3404 }
3405 return 0; // success
3406 }
3407 TK-=(w1[N]-beta)/u1;
3408 ++k;
3409 }
3410 return 1; // max number of iterations is exceeded
3411 }
3412
3413 /* the finite-difference method as an slower but more stable alternative to shooting method */
finitediff(const gen & f,const gen & x_idn,const gen & y_idn,const gen & dy_idn,const gen & x1,const gen & x2,const gen & y1,const gen & y2,int N,double tol,int M,vecteur & X,vecteur & Y,GIAC_CONTEXT)3414 int finitediff(const gen &f,const gen &x_idn,const gen &y_idn,const gen &dy_idn,const gen &x1,const gen &x2,
3415 const gen &y1,const gen &y2,int N,double tol,int M,vecteur &X,vecteur &Y,GIAC_CONTEXT) {
3416 gen dfy=derive(f,y_idn,contextptr),dfdy=derive(f,dy_idn,contextptr);
3417 vecteur vars=makevecteur(x_idn,y_idn,dy_idn);
3418 double a=x1.DOUBLE_val(),b=x2.DOUBLE_val(),alpha=y1.DOUBLE_val(),beta=y2.DOUBLE_val();
3419 double h=(b-a)/(N+1),fac=(beta-alpha)/(b-a)*h,x,t;
3420 vector<double> W(N+2,alpha),A(N+1),B(N+1),C(N+1),D(N+1),U(N+1),L(N+1),Z(N+1);
3421 vecteur V(N+1);
3422 W[N+1]=beta;
3423 for (int i=1;i<=N;++i) W[i]+=i*fac;
3424 int k=1;
3425 bool ef=true;
3426 dfy=simp(dfy,contextptr);
3427 dfdy=simp(dfdy,contextptr);
3428 while (k<=M) {
3429 x=a+h;
3430 t=(W[2]-alpha)/(2*h);
3431 A[1]=2+h*h*eval_func(dfy,vars,x,W[1],t,ef,contextptr);
3432 B[1]=h*eval_func(dfdy,vars,x,W[1],t,ef,contextptr)/2-1;
3433 D[1]=W[2]+alpha-2*W[1]-h*h*eval_func(f,vars,x,W[1],t,ef,contextptr);
3434 for (int i=2;i<N;++i) {
3435 x=a+i*h;
3436 t=(W[i+1]-W[i-1])/(2*h);
3437 A[i]=2+h*h*eval_func(dfy,vars,x,W[i],t,ef,contextptr);
3438 B[i]=h*eval_func(dfdy,vars,x,W[i],t,ef,contextptr)/2-1;
3439 C[i]=-h*eval_func(dfdy,vars,x,W[i],t,ef,contextptr)/2-1;
3440 D[i]=W[i+1]+W[i-1]-2*W[i]-h*h*eval_func(f,vars,x,W[i],t,ef,contextptr);
3441 }
3442 x=b-h;
3443 t=(beta-W[N-1])/(2*h);
3444 A[N]=2+h*h*eval_func(dfy,vars,x,W[N],t,ef,contextptr);
3445 C[N]=-h*eval_func(dfdy,vars,x,W[N],t,ef,contextptr)/2-1;
3446 D[N]=W[N-1]+beta-2*W[N]-h*h*eval_func(f,vars,x,W[N],t,ef,contextptr);
3447 if (!ef) return 2;
3448 L[1]=A[1];
3449 U[1]=B[1]/A[1];
3450 Z[1]=D[1]/L[1];
3451 for (int i=2;i<N;++i) {
3452 L[i]=A[i]-C[i]*U[i-1];
3453 U[i]=B[i]/L[i];
3454 Z[i]=(D[i]-C[i]*Z[i-1])/L[i];
3455 }
3456 L[N]=A[N]-C[N]*U[N-1];
3457 Z[N]=(D[N]-C[N]*Z[N-1])/L[N];
3458 V[N]=Z[N];
3459 W[N]+=Z[N];
3460 for (int i=N;i-->1;) {
3461 V[i]=gen(Z[i])-gen(U[i])*V[i+1];
3462 W[i]+=Z[i]-U[i]*V[i+1].DOUBLE_val();
3463 }
3464 if (is_greater(tol,_l2norm(V,contextptr),contextptr)) {
3465 X.resize(N+2);
3466 Y.resize(N+2);
3467 for (int i=0;i<=N+1;++i) {
3468 X[i]=a+i*h;
3469 Y[i]=W[i];
3470 }
3471 return 0; // success
3472 }
3473 ++k;
3474 }
3475 return 1; // maximum number of iterations is exceeded
3476 }
3477
_bvpsolve(const gen & g,GIAC_CONTEXT)3478 gen _bvpsolve(const gen &g,GIAC_CONTEXT) {
3479 if (g.type==_STRNG && g.subtype==-1) return g;
3480 if (g.type!=_VECT || g.subtype!=_SEQ__VECT)
3481 return gentypeerr(contextptr);
3482 vecteur &gv=*g._VECTptr;
3483 gen tk(undef);
3484 int maxiter=RAND_MAX;
3485 if (gv.size()<3 || gv[1].type!=_VECT || gv[2].type!=_VECT ||
3486 gv[1]._VECTptr->size()!=2 || (gv[2]._VECTptr->size()!=2 && gv[2]._VECTptr->size()!=3))
3487 return gensizeerr(contextptr);
3488 gen &f=gv.front(),&t=gv[1]._VECTptr->front(),&y=gv[1]._VECTptr->back(),
3489 y1=_evalf(gv[2]._VECTptr->at(0),contextptr),y2=_evalf(gv[2]._VECTptr->at(1),contextptr);
3490 if (gv[2]._VECTptr->size()==3 && (tk=_evalf(gv[2]._VECTptr->at(2),contextptr)).type!=_DOUBLE_)
3491 return gensizeerr(contextptr);
3492 if (y.type!=_IDNT || !t.is_symb_of_sommet(at_equal) ||
3493 t._SYMBptr->feuille._VECTptr->front().type!=_IDNT ||
3494 !t._SYMBptr->feuille._VECTptr->back().is_symb_of_sommet(at_interval))
3495 return gensizeerr(contextptr);
3496 gen &x=t._SYMBptr->feuille._VECTptr->front();
3497 vecteur &rng=*t._SYMBptr->feuille._VECTptr->back()._SYMBptr->feuille._VECTptr;
3498 gen x1=_evalf(rng.front(),contextptr),x2=_evalf(rng.back(),contextptr);
3499 if (x1.type!=_DOUBLE_ || x2.type!=_DOUBLE_ || y1.type!=_DOUBLE_ || y2.type!=_DOUBLE_ ||
3500 !is_strictly_greater(x2,x1,contextptr))
3501 return gensizeerr(contextptr);
3502 int N=100;
3503 /* parse options */
3504 int output_type=_BVP_LIST;
3505 for (const_iterateur it=gv.begin()+3;it!=gv.end();++it) {
3506 if (it->is_symb_of_sommet(at_equal)) {
3507 gen &lh=it->_SYMBptr->feuille._VECTptr->front();
3508 gen &rh=it->_SYMBptr->feuille._VECTptr->back();
3509 if (lh==at_output || lh==at_Output) {
3510 if (rh==_MAPLE_LIST)
3511 output_type=_BVP_LIST;
3512 else if (rh==at_derive)
3513 output_type=_BVP_DIFF;
3514 else if (rh==at_piecewise)
3515 output_type=_BVP_PIECEWISE;
3516 else if (rh==at_spline)
3517 output_type=_BVP_SPLINE;
3518 else return gensizeerr(contextptr);
3519 } else if (lh==at_limit) {
3520 if (!rh.is_integer() || (maxiter=rh.val)<1)
3521 return gensizeerr(contextptr);
3522 } else return gensizeerr(contextptr);
3523 } else if (it->is_integer()) {
3524 if ((N=it->val)<2)
3525 return gensizeerr(contextptr);
3526 } else return gensizeerr(contextptr);
3527 }
3528 gen dy=identificateur(" dy");
3529 gen F=subst(f,derive(symb_of(y,x),x,contextptr),dy,false,contextptr);
3530 F=subst(F,symb_of(y,x),y,false,contextptr);
3531 vecteur X,Y,dY;
3532 double tol=_evalf(_epsilon(change_subtype(vecteur(0),_SEQ__VECT),contextptr),contextptr).DOUBLE_val();
3533 int ec=shooting(F,x,y,dy,tk,x1,x2,y1,y2,N,tol,maxiter,X,Y,dY,contextptr);
3534 if (ec==1) {
3535 *logptr(contextptr) << "Error: maximum number of iterations exceeded\n";
3536 return undef;
3537 }
3538 if (ec==2) {
3539 *logptr(contextptr) << "Error: the shooting method failed to converge";
3540 if (is_undef(tk))
3541 *logptr(contextptr) << ", try to set an initial guess for y'(a)";
3542 *logptr(contextptr) << "\n";
3543 if (N>=3 && (output_type==_BVP_LIST || output_type==_BVP_PIECEWISE)) {
3544 *logptr(contextptr) << "Trying the finite-difference method instead\n";
3545 ec=finitediff(F,x,y,dy,x1,x2,y1,y2,N-1,tol,maxiter,X,Y,contextptr);
3546 if (ec==2) {
3547 *logptr(contextptr) << "Error: failed to converge\n";
3548 return undef;
3549 }
3550 if (ec==1) {
3551 *logptr(contextptr) << "Error: maximum number of iterations exceeded\n";
3552 return undef;
3553 }
3554 } else return undef;
3555 }
3556 vecteur res,coeff;
3557 matrice m=*_matrix(makesequence(4,4,0),contextptr)._VECTptr;
3558 m[0]._VECTptr->at(3)=m[1]._VECTptr->at(3)=m[2]._VECTptr->at(2)=m[3]._VECTptr->at(2)=1;
3559 switch (output_type) {
3560 case _BVP_LIST:
3561 res.reserve(N+1);
3562 for (int i=0;i<=N;++i)
3563 res.push_back(makevecteur(X[i],Y[i]));
3564 break;
3565 case _BVP_DIFF:
3566 res.reserve(N+1);
3567 for (int i=0;i<=N;++i)
3568 res.push_back(makevecteur(X[i],Y[i],dY[i]));
3569 break;
3570 case _BVP_PIECEWISE:
3571 case _BVP_SPLINE:
3572 res.reserve(2*(N+1)+1);
3573 for (int i=0;i<=N;++i) {
3574 res.push_back(i<N?symb_inferieur_strict(x,X[i]):symb_inferieur_egal(x,X[i]));
3575 if (i==0)
3576 res.push_back(0);
3577 else if (output_type==_BVP_PIECEWISE)
3578 res.push_back(Y[i-1]+(x-X[i-1])*(Y[i]-Y[i-1])/(X[i]-X[i-1]));
3579 else {
3580 m[0]._VECTptr->at(0)=pow(X[i-1],3);
3581 m[0]._VECTptr->at(1)=pow(X[i-1],2);
3582 m[0]._VECTptr->at(2)=X[i-1];
3583 m[1]._VECTptr->at(0)=pow(X[i],3);
3584 m[1]._VECTptr->at(1)=pow(X[i],2);
3585 m[1]._VECTptr->at(2)=X[i];
3586 m[2]._VECTptr->at(0)=3*m[0][1];
3587 m[2]._VECTptr->at(1)=2*X[i-1];
3588 m[3]._VECTptr->at(0)=3*m[1][1];
3589 m[3]._VECTptr->at(1)=2*X[i];
3590 coeff=*_linsolve(makesequence(m,makevecteur(Y[i-1],Y[i],dY[i-1],dY[i])),contextptr)._VECTptr;
3591 res.push_back(pow(x,3)*coeff[0]+pow(x,2)*coeff[1]+x*coeff[2]+coeff[3]);
3592 }
3593 }
3594 res.push_back(0);
3595 return symbolic(at_piecewise,change_subtype(res,_SEQ__VECT));
3596 default:
3597 break;
3598 }
3599 return res;
3600 }
3601 static const char _bvpsolve_s []="bvpsolve";
3602 static define_unary_function_eval (__bvpsolve,&_bvpsolve,_bvpsolve_s);
3603 define_unary_function_ptr5(at_bvpsolve,alias_at_bvpsolve,&__bvpsolve,0,true)
3604
strip_constants(const gen & g,GIAC_CONTEXT)3605 gen strip_constants(const gen &g,GIAC_CONTEXT) {
3606 if (g.is_symb_of_sommet(at_neg))
3607 return g._SYMBptr->feuille;
3608 if (g.is_symb_of_sommet(at_prod) && g._SYMBptr->feuille.type==_VECT) {
3609 const vecteur &feu=*g._SYMBptr->feuille._VECTptr;
3610 gen ret(1);
3611 for (const_iterateur it=feu.begin();it!=feu.end();++it) {
3612 if (_evalf(*it,contextptr).type==_DOUBLE_)
3613 continue;
3614 ret=ret*strip_constants(*it,contextptr);
3615 }
3616 return ret;
3617 }
3618 if (g.is_symb_of_sommet(at_inv))
3619 return _inv(strip_constants(g._SYMBptr->feuille,contextptr),contextptr);
3620 return g;
3621 }
3622
3623 /* return the expression for conjugate points */
_conjugate_equation(const gen & g,GIAC_CONTEXT)3624 gen _conjugate_equation(const gen &g,GIAC_CONTEXT) {
3625 if (g.type==_STRNG && g.subtype==-1) return g;
3626 if (g.type!=_VECT || g.subtype!=_SEQ__VECT)
3627 return gentypeerr(contextptr);
3628 vecteur &gv=*g._VECTptr;
3629 if (gv.size()!=5)
3630 return gensizeerr(contextptr);
3631 gen &y0=gv[0],&parm=gv[1],&pvals=gv[2],&t=gv[3],&a=gv[4];
3632 if (y0.type!=_SYMB || t.type!=_IDNT || parm.type!=_VECT || parm._VECTptr->size()!=2 ||
3633 pvals.type!=_VECT || pvals._VECTptr->size()!=2 ||
3634 _evalf(a,contextptr).type!=_DOUBLE_)
3635 return gensizeerr(contextptr);
3636 gen &alpha=parm._VECTptr->front(),&beta=parm._VECTptr->back();
3637 if (alpha.type!=_IDNT || beta.type!=_IDNT)
3638 return gensizeerr(contextptr);
3639 gen y1=derive(y0,alpha,contextptr),y2=derive(y0,beta,contextptr);
3640 gen ret=_collect(simp(subst(y1*subst(y2,t,a,false,contextptr)-y2*subst(y1,t,a,false,contextptr),
3641 parm,pvals,false,contextptr),contextptr),contextptr);
3642 return strip_constants(ret,contextptr);
3643 }
3644 static const char _conjugate_equation_s []="conjugate_equation";
3645 static define_unary_function_eval (__conjugate_equation,&_conjugate_equation,_conjugate_equation_s);
3646 define_unary_function_ptr5(at_conjugate_equation,alias_at_conjugate_equation,&__conjugate_equation,0,true)
3647
3648 static int cnst_count=0;
3649
3650 /* return the (list of) Euler-Lagrange equation(s) for functional L(u,du,t) */
_euler_lagrange(const gen & g,GIAC_CONTEXT)3651 gen _euler_lagrange(const gen &g,GIAC_CONTEXT) {
3652 if (g.type==_STRNG && g.subtype==-1) return g;
3653 gen L,t=identificateur("x");
3654 vecteur u=makevecteur(identificateur("y"));
3655 if (g.type!=_VECT) {
3656 L=g;
3657 if (!contains(lidnt(g),t))
3658 t=t__IDNT_e;
3659 } else {
3660 if (g.subtype!=_SEQ__VECT)
3661 return gensizeerr(contextptr);
3662 vecteur &gv=*g._VECTptr;
3663 L=gv.front();
3664 if (gv.size()>1) {
3665 if (gv.size()==2 && gv[1].is_symb_of_sommet(at_of)) {
3666 u.front()=gv[1]._SYMBptr->feuille._VECTptr->front();
3667 t=gv[1]._SYMBptr->feuille._VECTptr->back();
3668 } else if (gv.size()==2 && gv[1].type==_VECT && !gv[1]._VECTptr->empty()) {
3669 u.clear();
3670 t=undef;
3671 for (const_iterateur it=gv[1]._VECTptr->begin();it!=gv[1]._VECTptr->end();++it) {
3672 if (!it->is_symb_of_sommet(at_of))
3673 return gensizeerr(contextptr);
3674 u.push_back(it->_SYMBptr->feuille._VECTptr->front());
3675 if (is_undef(t))
3676 t=it->_SYMBptr->feuille._VECTptr->back();
3677 else if (t!=it->_SYMBptr->feuille._VECTptr->back())
3678 return gensizeerr(contextptr);
3679 }
3680 } else t=gv[1];
3681 if (t.type!=_IDNT)
3682 return gensizeerr(contextptr);
3683 }
3684 if (gv.size()>2) {
3685 if (gv[2].type==_IDNT)
3686 u.front()=gv[2];
3687 else if (gv[2].type==_VECT)
3688 u=*gv[2]._VECTptr;
3689 else return gensizeerr(contextptr);
3690 }
3691 }
3692 L=idnteval(L,contextptr);
3693 int n=u.size();
3694 vecteur du(n),Du(n),Dut(n),DU(n),D2U(n),d2u(n),ut(n);
3695 for (int i=0;i<n;++i) {
3696 if (u[i].type!=_IDNT)
3697 return gensizeerr(contextptr);
3698 ut[i]=symb_of(u[i],t);
3699 du[i]=identificateur(" du"+print_INT_(i));
3700 d2u[i]=identificateur(" d2u"+print_INT_(i));
3701 Du[i]=symbolic(at_derive,u[i]);
3702 Dut[i]=symb_of(Du[i],t);
3703 DU[i]=derive(ut[i],t,contextptr);
3704 D2U[i]=derive(ut[i],t,2,contextptr);
3705 }
3706 L=subst(L,Dut,du,false,contextptr);
3707 L=subst(L,Du,du,false,contextptr);
3708 L=subst(L,DU,du,false,contextptr);
3709 L=subst(L,ut,u,false,contextptr);
3710 vecteur ret;
3711 if (n==1 && !depend(L,*t._IDNTptr)) {
3712 ret.push_back(symb_equal(simp(du[0]*derive(L,du[0],contextptr)-L,contextptr),
3713 identificateur("K_"+print_INT_(cnst_count++))));
3714 }
3715 for (int i=0;i<n;++i) {
3716 if (!depend(L,*u[i]._IDNTptr)) {
3717 ret.push_back(symb_equal(simp(derive(L,du[i],contextptr),contextptr),
3718 identificateur("K_"+print_INT_(cnst_count++))));
3719 } else {
3720 gen eq=derive(L,u[i],contextptr),sol;
3721 eq-=derive(subst(derive(L,du[i],contextptr),makevecteur(du[i],u[i]),
3722 makevecteur(DU[i],ut[i]),false,contextptr),t,contextptr);
3723 eq=subst(eq,makevecteur(DU[i],D2U[i]),makevecteur(du[i],d2u[i]),false,contextptr);
3724 eq=symb_equal(simp(subst(eq,ut[i],u[i],false,contextptr),contextptr),0);
3725 if (depend(eq,*d2u[i]._IDNTptr) && (sol=_solve(makesequence(eq,d2u[i]),contextptr)).type==_VECT)
3726 eq=symb_equal(d2u[i],simp(sol._VECTptr->front(),contextptr));
3727 ret.push_back(eq);
3728 }
3729 }
3730 ret=subst(ret,ut,u,false,contextptr);
3731 gen tmprs=radsimp(ret,contextptr);
3732 if (tmprs.type==_VECT && tmprs._VECTptr->size()==ret.size())
3733 ret=*tmprs._VECTptr;
3734 ret=subst(subst(ret,u,ut,false,contextptr),mergevecteur(du,d2u),
3735 mergevecteur(DU,D2U),false,contextptr);
3736 return ret.size()==1?ret.front():ret;
3737 }
3738 static const char _euler_lagrange_s []="euler_lagrange";
3739 static define_unary_function_eval_quoted (__euler_lagrange,&_euler_lagrange,_euler_lagrange_s);
3740 define_unary_function_ptr5(at_euler_lagrange,alias_at_euler_lagrange,&__euler_lagrange,_QUOTE_ARGUMENTS,true)
3741
parse_functional(const gen & L,const gen & t,const gen & y,const gen & dy,GIAC_CONTEXT)3742 gen parse_functional(const gen &L,const gen &t,const gen &y,const gen &dy,GIAC_CONTEXT) {
3743 assert(t.type==_IDNT && y.type==_IDNT && dy.type==_IDNT);
3744 gen ret=subst(L,symb_of(symbolic(at_derive,y),t),dy,false,contextptr);
3745 ret=subst(ret,symbolic(at_derive,y),dy,false,contextptr);
3746 ret=subst(ret,derive(symb_of(y,t),t,contextptr),dy,false,contextptr);
3747 ret=subst(ret,symb_of(y,t),y,false,contextptr);
3748 return ret;
3749 }
3750
3751 /* return the Jacobi equation(s) for the functional L(y,y',a<=t<=b) and stationary function Y,
3752 * use h as unknown with h(a)=0 */
_jacobi_equation(const gen & g,GIAC_CONTEXT)3753 gen _jacobi_equation(const gen &g,GIAC_CONTEXT) {
3754 if (g.type==_STRNG && g.subtype==-1) return g;
3755 gen L,t,y,Y,h,a;
3756 if (g.type!=_VECT || g.subtype!=_SEQ__VECT)
3757 return gensizeerr(contextptr);
3758 vecteur &gv=*g._VECTptr;
3759 L=gv.front();
3760 if (gv.size()==5) {
3761 if (!gv[1].is_symb_of_sommet(at_of))
3762 return gensizeerr(contextptr);
3763 y=gv[1]._SYMBptr->feuille._VECTptr->front();
3764 t=gv[1]._SYMBptr->feuille._VECTptr->back();
3765 Y=gv[2]; h=gv[3]; a=gv[4];
3766 } else if (gv.size()==6) {
3767 t=gv[1]; y=gv[2]; Y=gv[3]; h=gv[4]; a=gv[5];
3768 } else return gensizeerr(contextptr);
3769 if (t.type!=_IDNT || h.type!=_IDNT || y.type!=_IDNT || _evalf(a,contextptr).type!=_DOUBLE_)
3770 return gensizeerr(contextptr);
3771 L=idnteval(L,contextptr);
3772 gen dy=identificateur(" dy");
3773 L=parse_functional(L,t,y,dy,contextptr);
3774 gen dY=derive(Y,t,contextptr);
3775 gen Ldydy=subst(simp(derive(L,dy,dy,contextptr),contextptr),
3776 makevecteur(y,dy),makevecteur(Y,dY),false,contextptr);
3777 gen Lyy=subst(simp(derive(L,y,y,contextptr),contextptr),
3778 makevecteur(y,dy),makevecteur(Y,dY),false,contextptr);
3779 gen Lydy=subst(simp(derive(L,y,dy,contextptr),contextptr),
3780 makevecteur(y,dy),makevecteur(Y,dY),false,contextptr);
3781 gen d=simp(Lyy-derive(Lydy,t,contextptr),contextptr);
3782 if (is_zero(Ldydy))
3783 return is_zero(d)?change_subtype(gen(0),_INT_BOOLEAN):symb_equal(h,0);
3784 gen jeq=simp(derive(Ldydy,t,contextptr),contextptr)*derive(symb_of(h,t),t,contextptr);
3785 jeq+=Ldydy*derive(symb_of(h,t),t,2,contextptr)-d*symb_of(h,t);
3786 jeq=symb_equal(jeq,0);
3787 gen sol=_deSolve(makesequence(makevecteur(jeq,symb_equal(symb_of(h,a),0)),makevecteur(t,h)),contextptr);
3788 if (sol.type==_VECT && sol._VECTptr->size()==1 && sol._VECTptr->front().type==_STRNG)
3789 return jeq;
3790 return makesequence(jeq,sol.is_symb_of_sommet(at_prod) || sol.is_symb_of_sommet(at_neg)?simp(sol,contextptr):sol);
3791 }
3792 static const char _jacobi_equation_s []="jacobi_equation";
3793 static define_unary_function_eval_quoted (__jacobi_equation,&_jacobi_equation,_jacobi_equation_s);
3794 define_unary_function_ptr5(at_jacobi_equation,alias_at_jacobi_equation,&__jacobi_equation,_QUOTE_ARGUMENTS,true)
3795
makevars(const gen & e,const gen & t,const vecteur & depvars,const vecteur & diffvars,GIAC_CONTEXT)3796 gen makevars(const gen &e,const gen &t,const vecteur &depvars,const vecteur &diffvars,GIAC_CONTEXT) {
3797 if (e.is_symb_of_sommet(at_of) && e._SYMBptr->feuille._VECTptr->back()==t) {
3798 gen &u=e._SYMBptr->feuille._VECTptr->front();
3799 for (const_iterateur it=depvars.begin();it!=depvars.end();++it) {
3800 if (*it==u) return u;
3801 if (u==symbolic(at_derive,*it)) return diffvars[it-depvars.begin()];
3802 }
3803 } else if (e.is_symb_of_sommet(at_derive)) {
3804 gen &feu=e._SYMBptr->feuille;
3805 if (feu.type!=_VECT || (feu._VECTptr->size()==2 && feu._VECTptr->at(1)==t)) {
3806 gen f=makevars(feu.type==_VECT?feu._VECTptr->front():feu,t,depvars,diffvars,contextptr);
3807 for (const_iterateur it=depvars.begin();it!=depvars.end();++it) {
3808 if (*it==f) return diffvars[it-depvars.begin()];
3809 }
3810 }
3811 } else if (e.type==_SYMB) {
3812 gen &feu=e._SYMBptr->feuille;
3813 if (feu.type==_VECT) {
3814 vecteur nf;
3815 nf.reserve(feu._VECTptr->size());
3816 for (const_iterateur it=feu._VECTptr->begin();it!=feu._VECTptr->end();++it) {
3817 nf.push_back(makevars(*it,t,depvars,diffvars,contextptr));
3818 }
3819 return symbolic(e._SYMBptr->sommet,change_subtype(nf,_SEQ__VECT));
3820 }
3821 return symbolic(e._SYMBptr->sommet,makevars(feu,t,depvars,diffvars,contextptr));
3822 }
3823 return _eval(e,contextptr);
3824 }
3825
apply_sign(const gen & g,const gen & simp,GIAC_CONTEXT)3826 gen apply_sign(const gen &g,const gen &simp,GIAC_CONTEXT) {
3827 gen e=_sign(g,contextptr);
3828 if (!e.is_symb_of_sommet(at_sign))
3829 return e;
3830 gen arg=_apply(makesequence(simp,vecteur(1,e._SYMBptr->feuille)),contextptr)._VECTptr->front();
3831 arg=_factor(arg,contextptr);
3832 if (is_zero(simp(arg-g,contextptr)))
3833 return e;
3834 return apply_sign(arg,simp,contextptr);
3835 }
3836
3837 /* return the sign of the expression, simplified */
determine_sign(const gen & e_orig,const gen & simp,GIAC_CONTEXT)3838 gen determine_sign(const gen &e_orig,const gen &simp,GIAC_CONTEXT) {
3839 gen e=_apply(makesequence(simp,vecteur(1,e_orig)),contextptr)._VECTptr->front();
3840 if (e.type==_SYMB)
3841 return apply_sign(_factor(e,contextptr),simp,contextptr);
3842 return _sign(e,contextptr);
3843 }
3844
strip_sign(const gen & g)3845 gen strip_sign(const gen &g) {
3846 if (g.is_symb_of_sommet(at_neg))
3847 return -strip_sign(g._SYMBptr->feuille);
3848 if (g.is_symb_of_sommet(at_sign))
3849 return g._SYMBptr->feuille;
3850 if (g.is_symb_of_sommet(at_prod) && g._SYMBptr->feuille.type==_VECT) {
3851 gen ret(1);
3852 vecteur &v=*g._SYMBptr->feuille._VECTptr;
3853 for (const_iterateur it=v.begin();it!=v.end();++it)
3854 ret=ret*strip_sign(*it);
3855 return ret;
3856 }
3857 return g;
3858 }
3859
3860 /* return the condition(s) under which the given function is convex */
_convex(const gen & g,GIAC_CONTEXT)3861 gen _convex(const gen &g,GIAC_CONTEXT) {
3862 if (g.type==_STRNG && g.subtype==-1) return g;
3863 if (g.type!=_VECT || g.subtype!=_SEQ__VECT)
3864 return gentypeerr(contextptr);
3865 vecteur &gv=*g._VECTptr,vars;
3866 gen f=gv.front(),t(undef);
3867 if (gv.size()<2)
3868 return gensizeerr(contextptr);
3869 f=idnteval(f,contextptr);
3870 if (gv[1].type==_VECT)
3871 vars=*gv[1]._VECTptr;
3872 else vars.push_back(gv[1]);
3873 if (vars.empty())
3874 return gensizeerr(contextptr);
3875 gen simp_func=at_simplify;
3876 if (gv.size()>2) {
3877 gen lh,rh;
3878 if (!gv[2].is_symb_of_sommet(at_equal) || (lh=gv[2]._SYMBptr->feuille._VECTptr->front())!=at_simplify ||
3879 !(rh=gv[2]._SYMBptr->feuille._VECTptr->back()).is_integer() || rh.subtype!=_INT_BOOLEAN)
3880 return false;
3881 if (rh.val==0)
3882 simp_func=at_ratnormal;
3883 }
3884 vecteur fvars,depvars;
3885 for (iterateur it=vars.begin();it!=vars.end();++it) {
3886 if (it->is_symb_of_sommet(at_of)) {
3887 gen &u=it->_SYMBptr->feuille._VECTptr->front();
3888 gen &v=it->_SYMBptr->feuille._VECTptr->back();
3889 if (v.type!=_IDNT || u.type!=_IDNT)
3890 return gensizeerr(contextptr);
3891 if (is_undef(t))
3892 t=v;
3893 else if (t!=v)
3894 return gensizeerr(contextptr);
3895 if (is_zero(_contains(makesequence(depvars,u),contextptr)))
3896 depvars.push_back(u);
3897 } else if (it->type==_IDNT) {
3898 if (is_zero(_contains(makesequence(fvars,*it),contextptr)))
3899 fvars.push_back(*it);
3900 } else return gensizeerr(contextptr);
3901 }
3902 vecteur diffvars,diffs;
3903 stringstream ss;
3904 int cnt=0;
3905 for (const_iterateur it=depvars.begin();it!=depvars.end();++it) {
3906 ss.str(""); ss << " tmp" << ++cnt;
3907 diffvars.push_back(identificateur(ss.str().c_str()));
3908 diffs.push_back(derive(symb_of(*it,t),t,contextptr));
3909 }
3910 gen F=is_undef(t)?_eval(f,contextptr):makevars(f,t,depvars,diffvars,contextptr);
3911 vecteur allvars=mergevecteur(fvars,mergevecteur(depvars,diffvars));
3912 gen TRUE=change_subtype(1,_INT_BOOLEAN),FALSE=change_subtype(0,_INT_BOOLEAN);
3913 if (allvars.size()==1) { // univariate function
3914 if (is_undef(simp_func))
3915 return derive(F,allvars.front(),2,contextptr);
3916 gen e=determine_sign(derive(F,allvars.front(),2,contextptr),simp_func,contextptr);
3917 if (is_one(e))
3918 return TRUE;
3919 if (is_minus_one(e))
3920 return FALSE;
3921 e=_apply(makesequence(simp_func,vecteur(1,strip_sign(e))),contextptr)._VECTptr->front();
3922 return e.is_symb_of_sommet(at_neg)?symb_inferieur_egal(e._SYMBptr->feuille,0):symb_superieur_egal(e,0);
3923 }
3924 // multivariate case, compute the eigenvalues of the Hessian
3925 int n=allvars.size(),N=std::pow(2,n);
3926 matrice H=*_hessian(makesequence(F,allvars),contextptr)._VECTptr;
3927 if (is_undef(simp_func))
3928 return subst(H,diffvars,diffs,false,contextptr);
3929 vecteur cond;
3930 cond.reserve(N);
3931 for (int i=1;i<N;++i) {
3932 matrice m_tmp,m;
3933 for (int j=0;j<n;++j) {
3934 if (((int)std::pow(2,j) & i)==0) continue;
3935 m_tmp.push_back(H[j]);
3936 }
3937 m_tmp=mtran(m_tmp);
3938 for (int j=0;j<n;++j) {
3939 if (((int)std::pow(2,j) & i)==0) continue;
3940 m.push_back(m_tmp[j]);
3941 }
3942 cond.push_back(determine_sign(_det(m,contextptr),simp_func,contextptr));
3943 }
3944 vecteur res;
3945 for (const_iterateur it=cond.begin();it!=cond.end();++it) {
3946 if (is_minus_one(*it))
3947 return FALSE;
3948 if (is_one(*it) || is_zero(*it))
3949 continue;
3950 res.push_back(strip_sign(*it));
3951 }
3952 if (res.empty())
3953 return TRUE;
3954 for (int i=0;i<int(res.size());++i) {
3955 for (int j=res.size();j-->i+1;) {
3956 if (is_zero(simp(res[j]-res[i],contextptr)))
3957 res.erase(res.begin()+j);
3958 }
3959 }
3960 for (int i=res.size();i-->0;) {
3961 if (res[i].is_integer() && res[i].subtype==_INT_BOOLEAN) {
3962 if (res[i].val==1)
3963 res.erase(res.begin()+i);
3964 else if (res[i].val==0)
3965 return FALSE;
3966 }
3967 gen &r=res[i];
3968 r=_apply(makesequence(simp_func,vecteur(1,r)),contextptr)._VECTptr->front();
3969 r=symb_superieur_egal(subst(r,diffvars,diffs,false,contextptr),0);
3970 }
3971 if (res.empty())
3972 return TRUE;
3973 gen simb=_lname(res,contextptr);
3974 if (simb.type==_VECT && simb._VECTptr->size()==1) {
3975 gen s=_solve(makesequence(res,simb._VECTptr->front()),contextptr);
3976 if (s.type==_VECT && s._VECTptr->size()>0)
3977 res=*s._VECTptr;
3978 }
3979 return res;
3980 }
3981 static const char _convex_s []="convex";
3982 static define_unary_function_eval_quoted (__convex,&_convex,_convex_s);
3983 define_unary_function_ptr5(at_convex,alias_at_convex,&__convex,_QUOTE_ARGUMENTS,true)
3984
3985 #ifndef NO_NAMESPACE_GIAC
3986 }
3987 #endif // ndef NO_NAMESPACE_GIAC
3988
3989