1 // -*- mode:C++ ; compile-command: "g++-3.4 -I.. -I../include -g -c maple.cc -DIN_GIAC -DHAVE_CONFIG_H" -*-
2 #include "giacPCH.h"
3 /*
4 * Copyright (C) 2000,2014 B. Parisse, Institut Fourier, 38402 St Martin d'Heres
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 #ifdef HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 using namespace std;
23 #include <stdexcept>
24 #include <cmath>
25 #include <cstdlib>
26 #include <stdio.h>
27 #if !defined GIAC_HAS_STO_38 && !defined NSPIRE && !defined FXCG && !defined POCKETCAS
28 #include <fstream>
29 #endif
30 #if !defined HAVE_NO_SYS_TIMES_H && defined HAVE_SYS_TIME_H
31 #include <fcntl.h>
32 #include <sys/time.h>
33 #include <time.h>
34 #endif
35 // #include <sys/resource.h>
36 // #include <unistd.h>
37 #include "sym2poly.h"
38 #include "usual.h"
39 #include "moyal.h"
40 #include "solve.h"
41 #include "intg.h"
42 #include "prog.h"
43 #include "misc.h"
44 #include "maple.h"
45 #include "plot.h"
46 #include "ifactor.h"
47 #include "subst.h"
48 #include "input_parser.h"
49 #include "sym2poly.h"
50 #include "modpoly.h"
51 #include "lin.h"
52 #include "derive.h"
53 #include "ti89.h"
54 #include "giacintl.h"
55 #ifdef HAVE_LIBGSL
56 #include <gsl/gsl_errno.h>
57 #include <gsl/gsl_fft_complex.h>
58 #include <gsl/gsl_fft_real.h>
59 #endif
60 #ifndef HAVE_PNG_H
61 #undef HAVE_LIBPNG
62 #endif
63 #ifdef HAVE_LIBPNG
64 #include <png.h>
65 #endif
66 #ifdef HAVE_UNISTD_H
67 #include <unistd.h>
68 #endif
69
70 #ifdef EMCC
71 #include <emscripten.h>
72 #endif
73
74
75 #if 0 // def GIAC_HAS_STO_38
76 TMillisecs PrimeGetNow();
77 #endif
78
79 #ifdef KHICAS
80 int time_shift;
81 #ifdef NUMWORKS
82 namespace Ion {
83 namespace Timing {
84
85 void usleep(uint32_t us);
86 void msleep(uint32_t ms);
87
88 /* millis is the number of milliseconds ellapsed since a random epoch.
89 * On the device, epoch is the boot time. */
90 volatile uint64_t millis();
91
92 }
93 }
millis()94 double millis(){
95 return double(Ion::Timing::millis()); // RTC_GetTicks();
96 }
97 #else
98 extern "C" double millis();
99 #endif
100 #endif
101
102
103 #if defined(EMCC) && !defined(PNACL)
104 extern "C" double emcctime();
105 // definition of emcctime should be added in emscripten directory src/library.js
106 // search for _time definition, and return only Date.now() for _emcctime
107 // otherwise time() will not work
108 #endif
109
110 #ifndef NO_NAMESPACE_GIAC
111 namespace giac {
112 #endif // ndef NO_NAMESPACE_GIAC
113
_zip(const gen & g,const context * contextptr)114 gen _zip(const gen & g,const context * contextptr){
115 if ( g.type==_STRNG && g.subtype==-1) return g;
116 int s=-1;
117 if (g.type!=_VECT || (s=g._VECTptr->size())<2)
118 return symbolic(at_zip,g);
119 vecteur & v=*g._VECTptr;
120 gen & f=v[0];
121 if (s==2){
122 if (f.type!=_VECT || v[1].type!=_VECT || f._VECTptr->size()!=v[1]._VECTptr->size())
123 return gendimerr(contextptr);
124 return _tran(g,contextptr);
125 }
126 if (v[1].type!=_VECT || v[2].type!=_VECT)
127 return f(gen(makevecteur(v[1],v[2]),_SEQ__VECT),contextptr);
128 vecteur & w1=*v[1]._VECTptr;
129 vecteur & w2=*v[2]._VECTptr;
130 int s1=int(w1.size()),s2=int(w2.size());
131 vecteur res;
132 int ss=giacmin(s1,s2),i=0;
133 res.reserve(ss);
134 for (;i<ss;++i)
135 res.push_back(_zip(gen(makevecteur(f,w1[i],w2[i]),_SEQ__VECT),contextptr));
136 if (s==3)
137 return res;
138 gen & ff=v[3];
139 for (;i<s1;++i)
140 res.push_back(_zip(gen(makevecteur(f,w1[i],ff),_SEQ__VECT),contextptr));
141 for (;i<s2;++i)
142 res.push_back(_zip(gen(makevecteur(f,ff,w2[i]),_SEQ__VECT),contextptr));
143 return res;
144 }
145 static const char _zip_s []="zip";
146 static define_unary_function_eval (__zip,&_zip,_zip_s);
147 define_unary_function_ptr5( at_zip ,alias_at_zip,&__zip,0,true);
148
_accumulate_head_tail(const gen & g,GIAC_CONTEXT)149 gen _accumulate_head_tail(const gen & g,GIAC_CONTEXT){
150 if ( g.type==_STRNG && g.subtype==-1) return g;
151 if (g.type!=_VECT || g._VECTptr->size()!=3)
152 return symbolic(at_accumulate_head_tail,g);
153 vecteur & v=*g._VECTptr;
154 if (v[0].type!=_VECT || v[1].type!=_INT_ || v[2].type!=_INT_)
155 return gensizeerr();
156 vecteur & w=*v[0]._VECTptr;
157 int tail=v[1].val;
158 int head=v[2].val;
159 int s=int(w.size());
160 if (tail<1 || head<1 || tail+head>s)
161 return gensizeerr();
162 gen tmp;
163 int i=0;
164 for (;i<tail;++i)
165 tmp=tmp+w[i];
166 vecteur res;
167 res.push_back(tmp);
168 int t=s-head;
169 for (;i<t;++i)
170 res.push_back(w[i]);
171 tmp=zero;
172 for (;i<s;++i)
173 tmp=tmp+w[i];
174 res.push_back(tmp);
175 return gen(res,v[0].subtype);
176 }
177 static const char _accumulate_head_tail_s []="accumulate_head_tail";
178 static define_unary_function_eval (__accumulate_head_tail,&_accumulate_head_tail,_accumulate_head_tail_s);
179 define_unary_function_ptr5( at_accumulate_head_tail ,alias_at_accumulate_head_tail,&__accumulate_head_tail,0,true);
180
_divide(const gen & g,GIAC_CONTEXT)181 gen _divide(const gen & g,GIAC_CONTEXT){
182 if ( g.type==_STRNG && g.subtype==-1) return g;
183 if (g.type!=_VECT || g._VECTptr->size()<2)
184 return gensizeerr();
185 vecteur & v=*g._VECTptr;
186 int s=int(v.size());
187 if (s==2)
188 return _quorem(g,contextptr);
189 gen q("Quo",contextptr),r("Rem",contextptr),arg;
190 gen mode=v.back();
191 if (s==4)
192 arg=gen(makevecteur(v[0],v[1],v[2]),_SEQ__VECT);
193 else
194 arg=gen(makevecteur(v[0],v[1]),_SEQ__VECT);
195 if (mode==q)
196 return _quo(arg,contextptr);
197 if (mode==r)
198 return _rem(arg,contextptr);
199 if (s==4)
200 return _quorem(arg,contextptr);
201 if (s==3 && mode.type==_IDNT)
202 return _quorem(g,contextptr);
203 return gensizeerr();
204 // return 0;
205 }
206 static const char _divide_s []="divide";
207 static define_unary_function_eval (__divide,&_divide,_divide_s);
208 define_unary_function_ptr5( at_divide ,alias_at_divide,&__divide,0,true);
209
_ratnormal(const gen & g,GIAC_CONTEXT)210 gen _ratnormal(const gen & g,GIAC_CONTEXT){
211 if ( g.type==_STRNG && g.subtype==-1) return g;
212 return ratnormal(g,contextptr);
213 }
214 static const char _ratnormal_s []="ratnormal";
215 static define_unary_function_eval (__ratnormal,&_ratnormal,_ratnormal_s);
216 define_unary_function_ptr5( at_ratnormal ,alias_at_ratnormal,&__ratnormal,0,true);
217
_about(const gen & g,GIAC_CONTEXT)218 gen _about(const gen & g,GIAC_CONTEXT){
219 if ( g.type==_STRNG && g.subtype==-1) return g;
220 if (g.type==_VECT)
221 return apply(g,contextptr,_about);
222 if (g.type==_IDNT)
223 return g._IDNTptr->eval(1,g,contextptr);
224 return g;
225 }
226 static const char _about_s []="about";
227 static define_unary_function_eval (__about,&_about,_about_s);
228 define_unary_function_ptr5( at_about ,alias_at_about,&__about,0,true);
229
_inverse(const gen & a_orig,GIAC_CONTEXT)230 gen _inverse(const gen & a_orig,GIAC_CONTEXT){
231 if ( a_orig.type==_STRNG && a_orig.subtype==-1) return a_orig;
232 matrice a;
233 bool convert_internal,minor_det,keep_pivot;
234 int algorithm,last_col;
235 if (!read_reduction_options(a_orig,a,convert_internal,algorithm,minor_det,keep_pivot,last_col))
236 return inv(a_orig,contextptr);
237 if (keep_pivot)
238 return gensizeerr(gettext("Option keep_pivot not applicable"));
239 if (minor_det){ // not really minors, use Leverrier algorithm
240 vecteur b;
241 vecteur p(mpcar(a,b,true,contextptr));
242 gen res=b.back()/p.back();
243 // if (a.size()%2==0)
244 res=-res;
245 return res;
246 }
247 matrice res;
248 if (!minv(a,res,convert_internal,algorithm,contextptr))
249 return gendimerr(contextptr);
250 return res;
251 }
252 static const char _inverse_s []="inverse";
253 static define_unary_function_eval (__inverse,&_inverse,_inverse_s);
254 define_unary_function_ptr5( at_inverse ,alias_at_inverse,&__inverse,0,true);
255
_Inverse(const gen & g,GIAC_CONTEXT)256 gen _Inverse(const gen & g,GIAC_CONTEXT){
257 if ( g.type==_STRNG && g.subtype==-1) return g;
258 return symbolic(at_inverse,g);
259 }
260 static const char _Inverse_s []="Inverse";
261 static define_unary_function_eval (__Inverse,&_Inverse,_Inverse_s);
262 define_unary_function_ptr5( at_Inverse ,alias_at_Inverse,&__Inverse,0,true);
263
_inverser(const gen & g,GIAC_CONTEXT)264 gen _inverser(const gen & g,GIAC_CONTEXT){
265 if ( g.type==_STRNG && g.subtype==-1) return g;
266 return sto(inv(eval(g,1,contextptr),contextptr),g,contextptr);
267 }
268 static const char _inverser_s []="inverser";
269 static define_unary_function_eval2_quoted (__inverser,&_inverser,_inverser_s,&printastifunction);
270 define_unary_function_ptr5( at_inverser ,alias_at_inverser,&__inverser,_QUOTE_ARGUMENTS,T_LOGO);
271
maple_gcdigcd(const gen & a_orig,const unary_function_ptr * u,GIAC_CONTEXT)272 static gen maple_gcdigcd(const gen & a_orig,const unary_function_ptr * u,GIAC_CONTEXT){
273 if (a_orig.type!=_VECT || a_orig._VECTptr->size()<2)
274 return gensizeerr();
275 vecteur v=*a_orig._VECTptr;
276 int s=int(v.size());
277 for (int i=2;i<s;++i){
278 v[i]=unmodunprod(v[i]);
279 }
280 if (s<4)
281 return (*u)(gen(v,a_orig.subtype),contextptr);
282 gen res;
283 if (s==4)
284 res=(*u)(makevecteur(v[0],v[1]),contextptr);
285 else
286 res=(*u)(makevecteur(v[0],v[1],v[2]),contextptr);
287 if (res.type==_VECT && res._VECTptr->size()==3 && v[s-2].type==_IDNT && v[s-1].type==_IDNT){
288 gen tmpsto=sto(res[0],v[s-2],contextptr);
289 if (is_undef(tmpsto))
290 return tmpsto;
291 tmpsto=sto(res[1],v[s-1],contextptr);
292 if (is_undef(tmpsto))
293 return tmpsto;
294 return res[2];
295 }
296 else
297 return res;
298 return gensizeerr();
299 }
_igcdex(const gen & a_orig,GIAC_CONTEXT)300 gen _igcdex(const gen & a_orig,GIAC_CONTEXT){
301 if ( a_orig.type==_STRNG && a_orig.subtype==-1) return a_orig;
302 return maple_gcdigcd(a_orig,at_iegcd,contextptr);
303 }
304 static const char _igcdex_s []="igcdex";
305 static define_unary_function_eval (__igcdex,&_igcdex,_igcdex_s);
306 define_unary_function_ptr5( at_igcdex ,alias_at_igcdex,&__igcdex,0,true);
307
308 static const char _igcd_s []="igcd";
309 static define_unary_function_eval (__igcd,&_gcd,_igcd_s);
310 define_unary_function_ptr5( at_igcd ,alias_at_igcd,&__igcd,0,true);
311
_gcdex(const gen & a_orig,GIAC_CONTEXT)312 gen _gcdex(const gen & a_orig,GIAC_CONTEXT){
313 if ( a_orig.type==_STRNG && a_orig.subtype==-1) return a_orig;
314 return maple_gcdigcd(a_orig,at_egcd,contextptr);
315 }
316 static const char _gcdex_s []="gcdex";
317 static define_unary_function_eval (__gcdex,&_gcdex,_gcdex_s);
318 define_unary_function_ptr5( at_gcdex ,alias_at_gcdex,&__gcdex,0,true);
319
320 static const char _indets_s []="indets";
321 static define_unary_function_eval (__indets,&_lname,_indets_s);
322 define_unary_function_ptr5( at_indets ,alias_at_indets,&__indets,0,true);
323
_revlist(const gen & a,GIAC_CONTEXT)324 gen _revlist(const gen & a,GIAC_CONTEXT){
325 if ( a.type==_STRNG && a.subtype==-1) return a;
326 if (a.type==_VECT){
327 vecteur v=*a._VECTptr;
328 reverse(v.begin(),v.end());
329 return gen(v,a.subtype);
330 }
331 if (a.type==_STRNG){
332 string s=*a._STRNGptr;
333 int l=int(s.size());
334 for (int i=0;i<l/2;++i){
335 char c=s[i];
336 s[i]=s[l-1-i];
337 s[l-1-i]=c;
338 }
339 return string2gen(s,false);
340 }
341 return a;
342 }
343 static const char _revlist_s []="revlist";
344 static define_unary_function_eval (__revlist,&_revlist,_revlist_s);
345 define_unary_function_ptr5( at_revlist ,alias_at_revlist,&__revlist,0,true);
346
347 static const char _reverse_s []="reverse";
348 static define_unary_function_eval (__reverse,&_revlist,_reverse_s);
349 define_unary_function_ptr5( at_reverse ,alias_at_reverse,&__reverse,0,true);
350
_restart(const gen & args,GIAC_CONTEXT)351 gen _restart(const gen & args,GIAC_CONTEXT){
352 if ( args.type==_STRNG && args.subtype==-1) return args;
353 intvar_counter=0;
354 realvar_counter=0;
355 if (args==at_solve) return 1;
356 init_context((context *) ((void *) contextptr));
357 gen res= _rm_all_vars(args,contextptr);
358 *logptr(contextptr) << "============== restarted ===============" << '\n';
359 if (args.type==_VECT && args.subtype==_SEQ__VECT && args._VECTptr->empty())
360 _srand(_time(gen(vecteur(0),_SEQ__VECT),contextptr),contextptr);
361 return res;
362 }
363 static const char _restart_s []="restart";
364 static define_unary_function_eval (__restart,&_restart,_restart_s);
365 define_unary_function_ptr5( at_restart ,alias_at_restart,&__restart,0,T_RETURN);
366
367 static const char _restart_vars_s []="restart_vars";
368 static define_unary_function_eval (__restart_vars,&_rm_all_vars,_restart_vars_s);
369 define_unary_function_ptr5( at_restart_vars ,alias_at_restart_vars,&__restart_vars,0,T_RETURN);
370
_restart_modes(const gen & args,GIAC_CONTEXT)371 gen _restart_modes(const gen & args,GIAC_CONTEXT){
372 if ( args.type==_STRNG && args.subtype==-1) return args;
373 if (!contextptr)
374 return 0;
375 init_context((context *) ((void *) contextptr));
376 return 1;
377 }
378 static const char _restart_modes_s []="restart_modes";
379 static define_unary_function_eval (__restart_modes,&_restart_modes,_restart_modes_s);
380 define_unary_function_ptr5( at_restart_modes ,alias_at_restart_modes,&__restart_modes,0,T_RETURN);
381
382 #ifdef KHICAS
_time(const gen & a,GIAC_CONTEXT)383 gen _time(const gen & a,GIAC_CONTEXT){
384 if ( a.type==_STRNG && a.subtype==-1) return a;
385 if (a.type==_VECT && a.subtype==_SEQ__VECT){
386 if (a._VECTptr->size()==2 && a._VECTptr->front().type==_INT_ && a._VECTptr->back().type==_INT_){
387 int h=a._VECTptr->front().val;
388 h=h%24;
389 if (h<0)
390 h+=24;
391 int m=a._VECTptr->back().val;
392 m=m%60;
393 if (m<0)
394 m+=60;
395 #ifdef NSPIRE_NEWLIB
396 unsigned NSPIRE_RTC_WADDR=0x90090008;
397 * (volatile unsigned *) NSPIRE_RTC_WADDR = (h*60+m)*60;
398 #else
399 time_shift=h*60+m;
400 #endif
401 return 1;
402 }
403 return millis();
404 }
405 double delta;
406 int ntimes=1,i=0;
407 int level=eval_level(contextptr);
408 double t1= millis(); // RTC_GetTicks(); // 1 tick=1/128 s
409 // CERR << t1 << endl;
410 for (unsigned i=1;i<=1000;++i){
411 eval(a,level,contextptr);
412 double t2= millis(); // RTC_GetTicks();
413 #ifdef NSPIRE_NEWLIB
414 if (t2>=t1+5000)
415 return double(t2-t1)/double(i)/1000;
416 #else
417 if (t2>=t1+32)
418 return double(t2-t1)/double(i)/1000;
419 #endif
420 }
421 return 0.0;
422 }
423
424 #else // KHICAS
425
_time(const gen & a,GIAC_CONTEXT)426 gen _time(const gen & a,GIAC_CONTEXT){
427 if ( a.type==_STRNG && a.subtype==-1) return a;
428 if (a.type==_VECT && a.subtype==_SEQ__VECT){
429 #ifdef GIAC_HAS_STO_38
430 return PrimeGetNow()/1000.;
431 #endif
432 #if 0 && defined(EMCC) && !defined(GIAC_GGB)
433 double res;
434 res=EM_ASM_DOUBLE_V({
435 var hw=Date.now();
436 return hw;
437 });
438 return res*1e-6;
439 #endif // GIAC_GGB
440
441 #if defined(EMCC) && !defined(PNACL)
442 return emcctime()/1e6;
443 #endif
444 return total_time(contextptr);
445 }
446 double delta;
447 int ntimes=1,i=0;
448 int level=eval_level(contextptr);
449 #if defined NSPIRE || defined NSPIRE_NEWLIB
450 unsigned NSPIRE_RTC_ADDR=0x90090000;
451 unsigned t1= * (volatile unsigned *) NSPIRE_RTC_ADDR;
452 // CERR << t1 << '\n';
453 for (unsigned i=1;i<=100;++i){
454 eval(a,level,contextptr);
455 unsigned t2= * (volatile unsigned *) NSPIRE_RTC_ADDR;
456 if (t2>=t1+3)
457 return double(t2-t1)/double(i);
458 }
459 return 0.0;
460 #endif
461 #if 0 && defined(EMCC) && !defined(GIAC_GGB)
462 double T1=EM_ASM_DOUBLE_V({
463 var hw=Date.now();
464 return hw;
465 });
466 eval(a,level,contextptr);
467 double T2=EM_ASM_DOUBLE_V({
468 var hw=Date.now();
469 return hw;
470 });
471 return (T2-T1)*1e-6;
472 #endif // GIAC_GGB
473 #if defined(EMCC) && !defined(PNACL)
474 // time_t t1,t2;
475 // time(&t1);
476 // eval(a,level,contextptr);
477 // time(&t2);
478 // return difftime(t2,t1);
479 double t1=emcctime();
480 eval(a,level,contextptr);
481 return (emcctime()-t1)/1e6;
482 #endif
483 #if defined(__APPLE__) || defined(PNACL)
484 unsigned u1=CLOCK();
485 struct timezone tz;
486 struct timeval debut,fin;
487 gettimeofday(&debut,&tz);
488 eval(a,level,contextptr);
489 gettimeofday(&fin,&tz);
490 u1=CLOCK()-u1;
491 return makevecteur(double(u1)/CLOCKS_PER_SEC,fin.tv_sec-debut.tv_sec+(fin.tv_usec-debut.tv_usec)/1e6);
492 #endif
493 #ifdef GIAC_HAS_STO_38
494 int t1=PrimeGetNow(),t2;
495 #endif
496 #ifdef _RUSAGE
497 struct rusage tmp1,tmp2,tmpc1,tmpc2;
498 getrusage(RUSAGE_SELF,&tmp1);
499 getrusage(RUSAGE_CHILDREN,&tmpc1);
500 #else // _RUSAGE
501 #ifdef HAVE_SYS_TIMES_H
502 struct tms tmp1,tmp2;
503 times(&tmp1);
504 #ifdef HAVE_LIBRT
505 struct timespec real1,real2;
506 clock_gettime(CLOCK_REALTIME,&real1);
507 #endif
508 #endif
509 #endif // _RUSAGE
510 for (;;){ // do it 10 times more
511 for (;i<ntimes;++i){
512 eval(a,level,contextptr);
513 }
514 #ifdef GIAC_HAS_STO_38
515 t2=PrimeGetNow();
516 delta=(t2-t1)/1000.;
517 #else // GIAC_HAS_STO_38
518 #ifdef _RUSAGE
519 getrusage(RUSAGE_SELF,&tmp2);
520 getrusage(RUSAGE_CHILDREN,&tmpc2);
521 delta = tmp2.ru_utime.tv_sec+1e-6*tmp2.ru_utime.tv_usec+tmp2.ru_stime.tv_sec+1e-6*tmp2.ru_stime.tv_usec-(tmp1.ru_utime.tv_sec+1e-6*tmp1.ru_utime.tv_usec+tmp1.ru_stime.tv_sec+1e-6*tmp1.ru_stime.tv_usec);
522 delta += tmpc2.ru_utime.tv_sec+1e-6*tmpc2.ru_utime.tv_usec+tmpc2.ru_stime.tv_sec+1e-6*tmpc2.ru_stime.tv_usec-(tmpc1.ru_utime.tv_sec+1e-6*tmpc1.ru_utime.tv_usec+tmpc1.ru_stime.tv_sec+1e-6*tmpc1.ru_stime.tv_usec);
523 #else // RUSAGE
524 #ifdef HAVE_SYS_TIMES_H
525 #ifdef HAVE_LIBRT
526 clock_gettime(CLOCK_REALTIME,&real2);
527 #endif // HAVE_LIBRT
528 times(&tmp2);
529 delta=delta_tms(tmp1,tmp2);
530 #else // HAVE_SYS_TIMES_H
531 delta=0.0;
532 break;
533 #endif // HAVE_SYS_TIMES_H
534 #endif // RUSAGE
535 #endif // GIAC_HAS_STO_38
536 if (delta>0.1)
537 break;
538 if (delta>0.05) // max wait time will be 2 seconds
539 ntimes *= 2;
540 else {
541 if (delta>0.02)
542 ntimes *= 5;
543 else
544 ntimes *= 10;
545 }
546 }
547 #ifdef GIAC_HAS_STO_38
548 return (delta/ntimes);
549 #endif
550 #ifdef HAVE_SYS_TIME_H
551 #ifdef HAVE_LIBRT
552 return makevecteur(delta/ntimes,(real2.tv_sec-real1.tv_sec+(real2.tv_nsec-real1.tv_nsec)/1e9)/ntimes);
553 #else
554 return (delta/ntimes);
555 #endif // HAVE_LIBRT
556 #else // HAVE_SYS_TIME_H
557 return (delta/ntimes);
558 #endif
559 }
560
561 #endif // KHICAS
562 static const char _time_s []="time";
563 static define_unary_function_eval_quoted (__time,&_time,_time_s);
564 define_unary_function_ptr5( at_time ,alias_at_time,&__time,_QUOTE_ARGUMENTS,true);
565
566 static const char _Phi_s []="Phi";
567 static define_unary_function_eval (__Phi,&_euler,_Phi_s);
568 define_unary_function_ptr5( at_Phi ,alias_at_Phi,&__Phi,0,true);
569
570 static const char _powermod_s []="powermod";
571 static define_unary_function_eval (__powermod,&_powmod,_powermod_s);
572 define_unary_function_ptr5( at_powermod ,alias_at_powermod,&__powermod,0,true);
573
574 static const char _length_s []="length";
575 static define_unary_function_eval (__length,&_size,_length_s);
576 define_unary_function_ptr5( at_length ,alias_at_length,&__length,0,true);
577
578 static const char _len_s []="len";
579 static define_unary_function_eval (__len,&_size,_len_s);
580 define_unary_function_ptr5( at_len ,alias_at_len,&__len,0,true);
581
582 static const char _nops_s []="nops";
583 static define_unary_function_eval (__nops,&_size,_nops_s);
584 define_unary_function_ptr5( at_nops ,alias_at_nops,&__nops,0,true);
585
_cat(const gen & a_orig,GIAC_CONTEXT)586 gen _cat(const gen & a_orig,GIAC_CONTEXT){
587 if ( a_orig.type==_STRNG && a_orig.subtype==-1) return a_orig;
588 vecteur v(gen2vecteur(a_orig));
589 string s;
590 const_iterateur it=v.begin(),itend=v.end();
591 for (;it!=itend;++it){
592 if (it->type==_STRNG)
593 s = s + *it->_STRNGptr;
594 else
595 s = s+it->print(contextptr);
596 }
597 return string2gen(s,false);
598 }
599 static const char _cat_s []="cat";
600 static define_unary_function_eval (__cat,&_cat,_cat_s);
601 define_unary_function_ptr5( at_cat ,alias_at_cat,&__cat,0,true);
602
_pivot(const gen & a_orig,GIAC_CONTEXT)603 gen _pivot(const gen & a_orig,GIAC_CONTEXT){
604 if ( a_orig.type==_STRNG && a_orig.subtype==-1) return a_orig;
605 vecteur v(gen2vecteur(a_orig));
606 int s=int(v.size());
607 if (s!=3 && s!=4)
608 return gensizeerr();
609 if (!ckmatrix(v.front()) || v[1].type!=_INT_ || v[2].type!=_INT_)
610 return gentypeerr();
611 matrice m=*v.front()._VECTptr;
612 int ml,mc;
613 mdims(m,ml,mc);
614 for (int i=0;i<ml;++i)
615 m[i]=*m[i]._VECTptr; // create a copy of the matrix
616 int l=v[1].val,c=v[2].val;
617 int shift = array_start(contextptr); // xcas_mode(contextptr)!=0 || abs_calc_mode(contextptr)==38;
618 l -= shift ;
619 c -= shift ;
620 if (l<0 || l>=ml || c<0 || c>=mc)
621 return gensizeerr();
622 gen p=m[l][c];
623 int l1=0,l2=ml-1;
624 if (s==4 && v[3].type==_INT_){
625 int lmin=v[3].val;
626 if (lmin<0){
627 l1 = -lmin;
628 l1 -= shift;
629 }
630 else {
631 lmin -= shift;
632 l2=giacmax(0,giacmin(lmin,ml-1));
633 l1=l2;
634 }
635 }
636 for (;l1<=l2;++l1){
637 if (l1!=l && !is_zero(m[l1][c]))
638 linear_combination(p,*m[l1]._VECTptr,-m[l1][c],*m[l]._VECTptr,plus_one,1,*m[l1]._VECTptr,epsilon(contextptr));
639 }
640 return m;
641 }
642 static const char _pivot_s []="pivot";
643 static define_unary_function_eval (__pivot,&_pivot,_pivot_s);
644 define_unary_function_ptr5( at_pivot ,alias_at_pivot,&__pivot,0,true);
645
rowcolspace(const gen & g,bool transpose,GIAC_CONTEXT)646 static gen rowcolspace(const gen & g,bool transpose,GIAC_CONTEXT){
647 if (g.type==_VECT && g._VECTptr->size()==2 && g._VECTptr->back().type!=_VECT){
648 gen tmp=rowcolspace(g._VECTptr->front(),transpose,contextptr),tmpsto;
649 if (transpose)
650 tmpsto=sto(int(tmp._VECTptr->front()._VECTptr->size()),g._VECTptr->back(),contextptr);
651 else
652 tmpsto=sto(int(tmp._VECTptr->size()),g._VECTptr->back(),contextptr);
653 if (is_undef(tmpsto))
654 return tmpsto;
655 return tmp;
656 }
657 if (!ckmatrix(g))
658 return gensizeerr(contextptr);
659 vecteur v=*g._VECTptr;
660 if (transpose)
661 v=mtran(v);
662 v=mrref(v,contextptr);
663 vecteur newv;
664 int s=int(v.size());
665 vecteur cmp(v.front()._VECTptr->size());
666 for (int i=0;i<s;++i){
667 if (v[i]!=cmp)
668 newv.push_back(v[i]);
669 }
670 if (transpose)
671 newv=mtran(newv);
672 return newv; // gen(newv,_SET__VECT);
673 }
_rowspace(const gen & g,GIAC_CONTEXT)674 gen _rowspace(const gen & g,GIAC_CONTEXT){
675 if ( g.type==_STRNG && g.subtype==-1) return g;
676 return rowcolspace(g,false,contextptr);
677 }
678 static const char _rowspace_s []="rowspace";
679 static define_unary_function_eval (__rowspace,&_rowspace,_rowspace_s);
680 define_unary_function_ptr5( at_rowspace ,alias_at_rowspace,&__rowspace,0,true);
681
_colspace(const gen & g,GIAC_CONTEXT)682 gen _colspace(const gen & g,GIAC_CONTEXT){
683 if ( g.type==_STRNG && g.subtype==-1) return g;
684 return rowcolspace(g,true,contextptr);
685 }
686 static const char _colspace_s []="colspace";
687 static define_unary_function_eval (__colspace,&_colspace,_colspace_s);
688 define_unary_function_ptr5( at_colspace ,alias_at_colspace,&__colspace,0,true);
689
690 static const char _nullspace_s []="nullspace";
691 static define_unary_function_eval (__nullspace,&_ker,_nullspace_s);
692 define_unary_function_ptr5( at_nullspace ,alias_at_nullspace,&__nullspace,0,true);
693
_copy(const gen & g,GIAC_CONTEXT)694 gen _copy(const gen & g,GIAC_CONTEXT){
695 if ( g.type==_STRNG && g.subtype==-1) return g;
696 if (g.type==_VECT){
697 vecteur v(*g._VECTptr);
698 iterateur it=v.begin(),itend=v.end();
699 for (;it!=itend;++it)
700 *it=_copy(*it,contextptr);
701 return gen(v,g.subtype);
702 }
703 if (g.type==_MAP)
704 return gen(*g._MAPptr);
705 return g;
706 }
707 static const char _copy_s []="copy";
708 static define_unary_function_eval (___copy,&_copy,_copy_s);
709 define_unary_function_ptr5( at_copy ,alias_at_copy,&___copy,0,true);
710
_row(const gen & g,GIAC_CONTEXT)711 gen _row(const gen & g,GIAC_CONTEXT){
712 if ( g.type==_STRNG && g.subtype==-1) return g;
713 if (g.type!=_VECT || g._VECTptr->size()!=2)
714 return gensizeerr();
715 int shift = array_start(contextptr); //xcas_mode(contextptr)!=0 || abs_calc_mode(contextptr)==38;
716 gen indice=g._VECTptr->back();
717 if (indice.is_symb_of_sommet(at_interval) && indice._SYMBptr->feuille.type==_VECT)
718 indice=symbolic(at_interval,indice._SYMBptr->feuille-gen(shift)*vecteur(indice._SYMBptr->feuille._VECTptr->size(),1));
719 else
720 indice -= shift;
721 gen res=g._VECTptr->front().operator_at(indice,contextptr);
722 if (ckmatrix(res))
723 return gen(*res._VECTptr,_SEQ__VECT);
724 else
725 return res;
726 }
727 static const char _row_s []="row";
728 static define_unary_function_eval (__row,&_row,_row_s);
729 define_unary_function_ptr5( at_row ,alias_at_row,&__row,0,true);
730
_col(const gen & g,GIAC_CONTEXT)731 gen _col(const gen & g,GIAC_CONTEXT){
732 if ( g.type==_STRNG && g.subtype==-1) return g;
733 if (g.type!=_VECT || g._VECTptr->size()!=2)
734 return gensizeerr();
735 int shift = array_start(contextptr); //xcas_mode(contextptr)!=0 || abs_calc_mode(contextptr)==38;
736 gen indice=g._VECTptr->back();
737 if (indice.is_symb_of_sommet(at_interval) && indice._SYMBptr->feuille.type==_VECT)
738 indice=symbolic(at_interval,indice._SYMBptr->feuille-gen(shift)*vecteur(indice._SYMBptr->feuille._VECTptr->size(),1));
739 else
740 indice -= shift;
741 gen res=_tran(g._VECTptr->front(),contextptr)[indice];
742 if (ckmatrix(res))
743 return gen(*res._VECTptr,_SEQ__VECT);
744 else
745 return res;
746 }
747 static const char _col_s []="col";
748 static define_unary_function_eval (__col,&_col,_col_s);
749 define_unary_function_ptr5( at_col ,alias_at_col,&__col,0,true);
750
count(const gen & f,const gen & v,const context * contextptr,const gen & param)751 static gen count(const gen & f,const gen & v,const context * contextptr,const gen & param){
752 if (v.type!=_VECT)
753 return f(v,contextptr);
754 const_iterateur it=v._VECTptr->begin(),itend=v._VECTptr->end();
755 if (param==at_row){
756 vecteur res;
757 for (;it!=itend;++it){
758 res.push_back(count(f,*it,contextptr,0));
759 }
760 return res;
761 }
762 if (param==at_col){
763 if (!ckmatrix(v))
764 return gentypeerr();
765 return count(f,mtran(*v._VECTptr),contextptr,at_row);
766 }
767 gen res;
768 for (;it!=itend;++it){
769 res=res+count(f,*it,contextptr,0);
770 }
771 return res;
772 }
_count(const gen & args,const context * contextptr)773 gen _count(const gen & args,const context * contextptr){
774 if ( args.type==_STRNG && args.subtype==-1) return args;
775 if ( (args.type!=_VECT) || (args._VECTptr->size()<2))
776 return gensizeerr(contextptr);
777 if (args.subtype!=_SEQ__VECT) {
778 if (!is_integer_vecteur(*args._VECTptr))
779 return gensizeerr(contextptr);
780 // count elements in list of integers
781 vector<int> x=vecteur_2_vector_int(*args._VECTptr);
782 int m=giacmin(x),M=giacmax(x),s=int(x.size());
783 if (M-m<3*s){
784 vector<int> eff(M-m+1);
785 effectif(x,eff,m);
786 matrice res;
787 for (int i=m;i<=M;++i){
788 int e=eff[i-m];
789 if (e>0)
790 res.push_back(makevecteur(i,e));
791 }
792 return res;
793 }
794 sort(x.begin(),x.end());
795 int old=M,eff=0,cur;
796 matrice res;
797 for (int pos=0;pos<s;++pos){
798 cur=x[pos];
799 if (cur==old){
800 ++eff;
801 continue;
802 }
803 if (eff)
804 res.push_back(makevecteur(old,eff));
805 old=cur;
806 eff=1;
807 }
808 if (eff)
809 res.push_back(makevecteur(old,eff));
810 return res;
811 }
812 gen v((*args._VECTptr)[1]);
813 gen f(args._VECTptr->front());
814 if (f.type==_STRNG && v.type==_STRNG){
815 // (Python-like) count occurences of v in f
816 int count=0,pos=-1,s=f._STRNGptr->size();
817 for (;pos<s;++count){
818 pos=f._STRNGptr->find(*v._STRNGptr,pos+1);
819 if (pos<0 || pos>=s)
820 break;
821 }
822 return count;
823 }
824 gen param;
825 if (args._VECTptr->size()>2)
826 param=(*args._VECTptr)[2];
827 else {
828 if (v.type!=_VECT)
829 return _count_eq(makesequence(v,f),contextptr);
830 }
831 return count(f,v,contextptr,param);
832 }
833 static const char _count_s []="count";
834 static define_unary_function_eval (___count,&_count,_count_s);
835 define_unary_function_ptr5( at_count ,alias_at_count,&___count,0,true);
836
count_eq0(const gen & f,const gen & v,GIAC_CONTEXT)837 static longlong count_eq0(const gen & f,const gen & v,GIAC_CONTEXT){
838 if (v.type!=_VECT)
839 return v==f;
840 const_iterateur it=v._VECTptr->begin(),itend=v._VECTptr->end();
841 longlong res=0;
842 for (;it!=itend;++it){
843 if (it->type!=_VECT){
844 if (*it==f)
845 ++res;
846 }
847 else
848 res += count_eq0(f,*it,contextptr);
849 }
850 return res;
851 }
852
count_eq1(const gen & f,const gen & v,GIAC_CONTEXT)853 static gen count_eq1(const gen & f,const gen & v,GIAC_CONTEXT){
854 if (v.type!=_VECT)
855 return is_strictly_greater(v,f,contextptr);
856 const_iterateur it=v._VECTptr->begin(),itend=v._VECTptr->end();
857 gen res;
858 for (;it!=itend;++it){
859 res=res+count_eq1(f,*it,contextptr);
860 }
861 return res;
862 }
863
count_eq2(const gen & f,const gen & v,GIAC_CONTEXT)864 static gen count_eq2(const gen & f,const gen & v,GIAC_CONTEXT){
865 if (v.type!=_VECT)
866 return is_strictly_greater(f,v,contextptr);
867 const_iterateur it=v._VECTptr->begin(),itend=v._VECTptr->end();
868 gen res;
869 for (;it!=itend;++it){
870 res=res+count_eq2(f,*it,contextptr);
871 }
872 return res;
873 }
874
count_eq(const gen & f,const gen & v,GIAC_CONTEXT,int type,const gen & rowcol)875 static gen count_eq(const gen & f,const gen & v,GIAC_CONTEXT,int type,const gen & rowcol){
876 if (rowcol==at_row || rowcol==at_col){
877 if (!ckmatrix(v))
878 return gentypeerr();
879 if (rowcol==at_row){
880 const_iterateur it=v._VECTptr->begin(),itend=v._VECTptr->end();
881 vecteur res;
882 for (;it!=itend;++it){
883 res.push_back(count_eq(f,*it,contextptr,type,0));
884 }
885 return res;
886 }
887 if (rowcol==at_col)
888 return count_eq(f,mtran(*v._VECTptr),contextptr,type,at_row);
889 }
890 if (type==0)
891 return count_eq0(f,v,contextptr);
892 if (type==1)
893 return count_eq1(f,v,contextptr);
894 if (type==2)
895 return count_eq2(f,v,contextptr);
896 return gentypeerr();
897 #if 0
898 if (v.type!=_VECT){
899 switch (type){
900 case 0:
901 return v==f;
902 case 1:
903 return is_strictly_greater(v,f,contextptr);
904 case 2:
905 return is_strictly_greater(f,v,contextptr);
906 default:
907 return gentypeerr();
908 }
909 }
910 const_iterateur it=v._VECTptr->begin(),itend=v._VECTptr->end();
911 gen res;
912 for (;it!=itend;++it){
913 res=res+count_eq(f,*it,contextptr,type,0);
914 }
915 return res;
916 #endif
917 }
918
_count_eq(const gen & args,const context * contextptr)919 gen _count_eq(const gen & args,const context * contextptr){
920 if ( args.type==_STRNG && args.subtype==-1) return args;
921 if ( args.type!=_VECT || args.subtype!=_SEQ__VECT || args._VECTptr->size()<2)
922 return gensizeerr(contextptr);
923 gen v((*args._VECTptr)[1]);
924 gen f(args._VECTptr->front());
925 gen param;
926 if (args._VECTptr->size()>2)
927 param=(*args._VECTptr)[2];
928 return count_eq(f,v,contextptr,0,param);
929 }
930 static const char _count_eq_s []="count_eq";
931 static define_unary_function_eval (___count_eq,&_count_eq,_count_eq_s);
932 define_unary_function_ptr5( at_count_eq ,alias_at_count_eq,&___count_eq,0,true);
933
_count_sup(const gen & args,GIAC_CONTEXT)934 gen _count_sup(const gen & args,GIAC_CONTEXT){
935 if ( args.type==_STRNG && args.subtype==-1) return args;
936 if ( args.type!=_VECT || args.subtype!=_SEQ__VECT || args._VECTptr->size()<2)
937 return gensizeerr(contextptr);
938 gen v((*args._VECTptr)[1]);
939 gen f(args._VECTptr->front());
940 gen param;
941 if (args._VECTptr->size()>2)
942 param=(*args._VECTptr)[2];
943 return count_eq(f,v,contextptr,1,param);
944 }
945 static const char _count_sup_s []="count_sup";
946 static define_unary_function_eval (___count_sup,&_count_sup,_count_sup_s);
947 define_unary_function_ptr5( at_count_sup ,alias_at_count_sup,&___count_sup,0,true);
948
_count_inf(const gen & args,GIAC_CONTEXT)949 gen _count_inf(const gen & args,GIAC_CONTEXT){
950 if ( args.type==_STRNG && args.subtype==-1) return args;
951 if ( args.type!=_VECT || args.subtype!=_SEQ__VECT || args._VECTptr->size()<2)
952 return gensizeerr(contextptr);
953 gen v((*args._VECTptr)[1]);
954 gen f(args._VECTptr->front());
955 gen param;
956 if (args._VECTptr->size()>2)
957 param=(*args._VECTptr)[2];
958 return count_eq(f,v,contextptr,2,param);
959 }
960 static const char _count_inf_s []="count_inf";
961 static define_unary_function_eval (___count_inf,&_count_inf,_count_inf_s);
962 define_unary_function_ptr5( at_count_inf ,alias_at_count_inf,&___count_inf,0,true);
963
_trunc(const gen & args,GIAC_CONTEXT)964 gen _trunc(const gen & args,GIAC_CONTEXT){
965 if ( args.type==_STRNG && args.subtype==-1) return args;
966 if (is_equal(args))
967 return apply_to_equal(args,_trunc,contextptr);
968 if (args.is_symb_of_sommet(at_unit))
969 return apply_unit(args,_trunc,contextptr);
970 if (args.type==_VECT) {
971 if (args.subtype==_SEQ__VECT && args._VECTptr->size()==2 && args._VECTptr->back().type==_INT_){
972 #ifdef BCD
973 if (args._VECTptr->front().type==_FLOAT_)
974 return ftrunc(args._VECTptr->front()._FLOAT_val,args._VECTptr->back().val);
975 #endif
976 gen b=args._VECTptr->back();
977 if (b.val<0){
978 gen gf=_floor(log10(abs(args._VECTptr->front(),contextptr),contextptr),contextptr);
979 if (gf.type!=_INT_ && gf.type!=_FLOAT_)
980 return gensizeerr(contextptr);
981 b=-1-b-gf;
982 }
983 #ifdef _SOFTMATH_H
984 double d=std::giac_gnuwince_pow(10.0,double(b.val));
985 #else
986 double d=std::pow(10.0,double(b.val));
987 #endif
988 return _floor(d*args._VECTptr->front(),contextptr)/d;
989 }
990 if (args.subtype==_SEQ__VECT && args._VECTptr->size()==2 && args._VECTptr->front().type==_INT_){
991 if (args._VECTptr->front().val)
992 return zero;
993 return _trunc(args._VECTptr->back(),contextptr);
994 }
995 return apply(args,contextptr,_trunc);
996 }
997 if (is_strictly_positive(-args,contextptr))
998 return -_floor(-args,contextptr);
999 return _floor(args,contextptr);
1000 }
1001 static const char _trunc_s []="trunc";
1002 static define_unary_function_eval (___trunc,&_trunc,_trunc_s);
1003 define_unary_function_ptr5( at_trunc ,alias_at_trunc,&___trunc,0,true);
1004
1005 static const char _giac_s []="giac";
1006 static define_unary_function_eval (__giac,&_version,_giac_s);
1007 define_unary_function_ptr5( at_giac ,alias_at_giac,&__giac,0,true);
1008
1009 static const char _div_s []="div";
1010 static define_unary_function_eval4 (__div,&_iquo,_div_s,&printsommetasoperator,&texprintsommetasoperator);
1011 define_unary_function_ptr5( at_div ,alias_at_div,&__div,0,T_MOD);
1012
_evalc(const gen & g,GIAC_CONTEXT)1013 gen _evalc(const gen & g,GIAC_CONTEXT){
1014 if ( g.type==_STRNG && g.subtype==-1) return g;
1015 if (g.type==_VECT)
1016 return apply(g,_evalc,contextptr);
1017 gen tmp(_exp2pow(_lin(recursive_normal(g,contextptr),contextptr),contextptr));
1018 vecteur l(lop(tmp,at_arg));
1019 if (!l.empty()){
1020 vecteur lp=*apply(l,gen_feuille)._VECTptr;
1021 lp=*apply(lp,contextptr,arg_CPLX)._VECTptr;
1022 tmp=subst(tmp,l,lp,false,contextptr);
1023 }
1024 tmp=recursive_normal(tmp,contextptr);
1025 vecteur vtmp(lvar(tmp));
1026 if (vtmp.size()==1 && vtmp[0].is_symb_of_sommet(at_exp)){
1027 tmp=ratnormal(_halftan(_exp2trig(tmp,contextptr),contextptr),contextptr);
1028 }
1029 gen r,i;
1030 reim(tmp,r,i,contextptr);
1031 gen tmp2=_lin(inv(tmp,contextptr),contextptr);
1032 gen re2,im2;
1033 reim(tmp2,re2,im2,contextptr);
1034 if (lvar(makevecteur(re2,im2)).size()<lvar(makevecteur(r,i)).size())
1035 reim(inv(ratnormal(re2,contextptr)+cst_i*ratnormal(im2,contextptr),contextptr),r,i,contextptr);
1036 // tmp=simplify(i,contextptr);
1037 if (is_zero(i))
1038 return r;
1039 if (is_zero(r))
1040 return cst_i*i;
1041 return symbolic(at_plus,gen(makevecteur(r,cst_i*i),_SEQ__VECT));
1042 }
1043 static const char _evalc_s []="evalc";
1044 static define_unary_function_eval (__evalc,&_evalc,_evalc_s);
1045 define_unary_function_ptr5( at_evalc ,alias_at_evalc,&__evalc,0,true);
1046
1047 static const char _evala_s []="evala";
1048 static define_unary_function_eval (__evala,(const gen_op_context)recursive_normal,_evala_s);
1049 define_unary_function_ptr5( at_evala ,alias_at_evala,&__evala,0,true);
1050
1051 #ifndef S_IRUSR
1052 #define S_IRUSR 00400
1053 #endif
1054 #ifndef S_IWUSR
1055 #define S_IWUSR 00200
1056 #endif
1057
1058 // open a file, returns a FD
_open(const gen & g,GIAC_CONTEXT)1059 gen _open(const gen & g,GIAC_CONTEXT){
1060 if ( g.type==_STRNG && g.subtype==-1) return g;
1061 #if defined(VISUALC) || defined(__MINGW_H) || defined (FIR) || defined(FXCG) || defined(NSPIRE) || defined(__ANDROID__) || defined(NSPIRE_NEWLIB) || defined(EMCC) || defined(GIAC_GGB) || defined KHICAS
1062 return gensizeerr(gettext("not implemented"));
1063 #else
1064 gen tmp=check_secure();
1065 if (is_undef(tmp)) return tmp;
1066 gen filename(g);
1067 if (filename.type!=_STRNG)
1068 return gensizeerr();
1069 int access=O_RDWR | O_CREAT | O_TRUNC | O_APPEND ;
1070 int res=open(filename._STRNGptr->c_str(),access, S_IRUSR|S_IWUSR);
1071 if (res==-1)
1072 return gensizeerr(gettext("Could not open file"));
1073 gen r(res);
1074 r.subtype=_INT_FD;
1075 return r;
1076 #endif
1077 }
1078 static const char _open_s []="open";
1079 static define_unary_function_eval (__open,&_open,_open_s);
1080 define_unary_function_ptr5( at_open ,alias_at_open,&__open,0,true);
1081
1082 // open a file, returns a FILE *
_fopen(const gen & g,GIAC_CONTEXT)1083 gen _fopen(const gen & g,GIAC_CONTEXT){
1084 if ( g.type==_STRNG && g.subtype==-1) return g;
1085 gen tmp=check_secure();
1086 if (is_undef(tmp)) return tmp;
1087 gen filename(g);
1088 string mode="w+";
1089 if (g.type==_VECT && g.subtype==_SEQ__VECT && g._VECTptr->size()==2 && g._VECTptr->back().type==_STRNG){
1090 filename=g._VECTptr->front();
1091 mode=*g._VECTptr->back()._STRNGptr;
1092 }
1093 if (filename.type!=_STRNG || mode.size()>2)
1094 return gensizeerr();
1095 FILE * f=fopen(filename._STRNGptr->c_str(),mode.c_str());
1096 return gen((void *) f,_FILE_POINTER);
1097 }
1098 static const char _fopen_s []="fopen";
1099 static define_unary_function_eval (__fopen,&_fopen,_fopen_s);
1100 define_unary_function_ptr5( at_fopen ,alias_at_fopen,&__fopen,0,true);
1101
1102 // fprint first arg=FD, rest is printed
_fprint(const gen & g,GIAC_CONTEXT)1103 gen _fprint(const gen & g,GIAC_CONTEXT){
1104 if ( g.type==_STRNG && g.subtype==-1) return g;
1105 gen tmp=check_secure();
1106 if (is_undef(tmp)) return tmp;
1107 if (g.type!=_VECT || g._VECTptr->size()<1)
1108 return gensizeerr(gettext("1st arg=open result, then other args"));
1109 vecteur & v=*g._VECTptr;
1110 int s=int(v.size());
1111 FILE * f=0;
1112 #if !defined(BESTA_OS) && !defined(NSPIRE) && !defined(FXCG)
1113 if (v[0].type==_INT_ && v[0].subtype==_INT_FD)
1114 f= fdopen(v[0].val,"a");
1115 #endif
1116 if (v[0].type==_POINTER_ && v[0].subtype==_FILE_POINTER)
1117 f=(FILE *) v[0]._POINTER_val;
1118 if (f){
1119 if (s>1 && v[1]==gen("Unquoted",contextptr)){
1120 for (int i=2;i<s;++i)
1121 fprintf(f,"%s",v[i].type==_STRNG?v[i]._STRNGptr->c_str():unquote(v[i].print(contextptr)).c_str());
1122 }
1123 else {
1124 for (int i=1;i<s;++i)
1125 fprintf(f,"%s",v[i].print(contextptr).c_str());
1126 }
1127 // fclose(f);
1128 return plus_one;
1129 }
1130 return zero;
1131 }
1132 static const char _fprint_s []="fprint";
1133 static define_unary_function_eval (__fprint,&_fprint,_fprint_s);
1134 define_unary_function_ptr5( at_fprint ,alias_at_fprint,&__fprint,0,true);
1135
_close(const gen & g0,GIAC_CONTEXT)1136 gen _close(const gen & g0,GIAC_CONTEXT){
1137 gen g=eval(g0,1,contextptr);
1138 if ( g.type==_STRNG && g.subtype==-1) return g;
1139 #if !defined(VISUALC) && !defined(BESTA_OS) && !defined(__MINGW_H) && !defined(NSPIRE) && !defined(FXCG)
1140 if (g.type==_INT_ && g.subtype==_INT_FD){
1141 purgenoassume(g0,contextptr);
1142 close(g.val);
1143 return plus_one;
1144 }
1145 #endif
1146 if (g.type==_POINTER_){
1147 purgenoassume(g0,contextptr);
1148 fclose((FILE *)g._POINTER_val);
1149 return plus_one;
1150 }
1151 return zero;
1152 }
1153 static const char _close_s []="close";
1154 static define_unary_function_eval_quoted (__close,&_close,_close_s);
1155 define_unary_function_ptr5( at_close ,alias_at_close,&__close,_QUOTE_ARGUMENTS,true);
1156
1157 static const char _fclose_s []="fclose";
1158 static define_unary_function_eval_quoted (__fclose,&_close,_fclose_s);
1159 define_unary_function_ptr5( at_fclose ,alias_at_fclose,&__fclose,_QUOTE_ARGUMENTS,true);
1160
_blockmatrix(const gen & g,GIAC_CONTEXT)1161 gen _blockmatrix(const gen & g,GIAC_CONTEXT){
1162 if ( g.type==_STRNG && g.subtype==-1) return g;
1163 if (g.type!=_VECT || g._VECTptr->size()!=3)
1164 return gentypeerr();
1165 vecteur & v = *g._VECTptr;
1166 if (v[0].type!=_INT_ || v[1].type!=_INT_ || v[2].type!=_VECT)
1167 return gensizeerr();
1168 int n=v[0].val,m=v[1].val;
1169 vecteur l=*v[2]._VECTptr;
1170 int s=int(l.size());
1171 if (n<=0 || m<=0 || n*m!=s)
1172 return gendimerr();
1173 for (int k=0;k<s;++k){
1174 if (!ckmatrix(l[k])){
1175 vecteur vtmp(vecteur(1,l[k]));
1176 if (ckmatrix(vtmp))
1177 l[k]=mtran(vtmp);
1178 if (!ckmatrix(l[k]))
1179 return gentypeerr();
1180 }
1181 }
1182 matrice res;
1183 int pos=0;
1184 unsigned final=0;
1185 for (int i=0;i<n;++i){
1186 // glue m matrices from l vertically, these matrices must have same #rows
1187 pos=i*m;
1188 int nrows=int(l[pos]._VECTptr->size());
1189 for (int j=1;j<m;++j){
1190 if (l[pos+j]._VECTptr->size()!=unsigned(nrows))
1191 return gendimerr();
1192 }
1193 for (int k=0;k<nrows;++k){
1194 // glue row k of l[pos+j] for j=0..m
1195 vecteur tmp;
1196 for (int j=0;j<m;++j){
1197 gen cur_row = (*l[pos+j]._VECTptr)[k];
1198 const_iterateur it=cur_row._VECTptr->begin(),itend=cur_row._VECTptr->end();
1199 for (;it!=itend;++it)
1200 tmp.push_back(*it);
1201 }
1202 if (final && tmp.size()!=final)
1203 return gendimerr();
1204 else
1205 final=unsigned(tmp.size());
1206 res.push_back(tmp);
1207 }
1208 }
1209 return res;
1210 }
1211 static const char _blockmatrix_s []="blockmatrix";
1212 static define_unary_function_eval (__blockmatrix,&_blockmatrix,_blockmatrix_s);
1213 define_unary_function_ptr5( at_blockmatrix ,alias_at_blockmatrix,&__blockmatrix,0,true);
1214
delrowscols(const gen & g,bool isrow,GIAC_CONTEXT)1215 static gen delrowscols(const gen & g,bool isrow,GIAC_CONTEXT){
1216 gen gm,interval,f,fa,fb;
1217 if (g.type!=_VECT || g._VECTptr->size()!=2 )
1218 return gentypeerr();
1219 gm=g._VECTptr->front();
1220 if (is_Ans(gm))
1221 gm=eval(gm,1,contextptr);
1222 if (gm.type==_IDNT){
1223 return sto(delrowscols(eval(g,eval_level(contextptr),contextptr),isrow,contextptr),gm,contextptr);
1224 }
1225 interval=g._VECTptr->back();
1226 // if (interval.type==_IDNT || interval.type==_SYMB)
1227 interval=eval(interval,1,contextptr);
1228 if (!interval.is_symb_of_sommet(at_interval))
1229 interval=symb_interval(interval,interval);
1230 if (!ckmatrix(gm) || !interval.is_symb_of_sommet(at_interval) || (f=interval._SYMBptr->feuille).type!=_VECT || f._VECTptr->size()!=2 || !is_integral(fa=f._VECTptr->front()) || !is_integral(fb=f._VECTptr->back()) )
1231 return gentypeerr();
1232 int shift = array_start(contextptr); //xcas_mode(contextptr)!=0 || abs_calc_mode(contextptr)==38;
1233 int a=fa.val-shift,b=fb.val-shift,s;
1234 matrice m=*gm._VECTptr;
1235 if (!isrow)
1236 m=mtran(m);
1237 s=int(m.size());
1238 if (a>=s || b>=s || a<0 || b<0 || a>b)
1239 return gendimerr();
1240 m.erase(m.begin()+a,m.begin()+b+1);
1241 if (!isrow)
1242 m=mtran(m);
1243 return m;
1244 }
1245
_delrows(const gen & g,GIAC_CONTEXT)1246 gen _delrows(const gen & g,GIAC_CONTEXT){
1247 if ( g.type==_STRNG && g.subtype==-1) return g;
1248 return delrowscols(g,true,contextptr);
1249 }
1250 static const char _delrows_s []="delrows";
1251 static define_unary_function_eval (__delrows,&_delrows,_delrows_s);
1252 define_unary_function_ptr5( at_delrows ,alias_at_delrows,&__delrows,0,true);
1253
_delcols(const gen & g,GIAC_CONTEXT)1254 gen _delcols(const gen & g,GIAC_CONTEXT){
1255 if ( g.type==_STRNG && g.subtype==-1) return g;
1256 return delrowscols(g,false,contextptr);
1257 }
1258 static const char _delcols_s []="delcols";
1259 static define_unary_function_eval (__delcols,&_delcols,_delcols_s);
1260 define_unary_function_ptr5( at_delcols ,alias_at_delcols,&__delcols,0,true);
1261
1262 static const char _gaussjord_s []="gaussjord";
1263 static define_unary_function_eval (__gaussjord,&_rref,_gaussjord_s);
1264 define_unary_function_ptr5( at_gaussjord ,alias_at_gaussjord,&__gaussjord,0,true);
1265
_JordanBlock(const gen & g,GIAC_CONTEXT)1266 gen _JordanBlock(const gen & g,GIAC_CONTEXT){
1267 if ( g.type==_STRNG && g.subtype==-1) return g;
1268 int s;
1269 if (g.type!=_VECT || g._VECTptr->size()!=2 || g._VECTptr->back().type!=_INT_ )
1270 return gentypeerr();
1271 s=g._VECTptr->back().val;
1272 if (s<=0 || double(s)*s>LIST_SIZE_LIMIT)
1273 return gendimerr();
1274 --s;
1275 gen x=g._VECTptr->front();
1276 matrice m;
1277 m.reserve(s+1);
1278 for (int i=0;i<=s;++i){
1279 #ifdef TIMEOUT
1280 control_c();
1281 #endif
1282 if (ctrl_c || interrupted)
1283 return gensizeerr(gettext("Stopped by user interruption."));
1284 vecteur v(s+1);
1285 v[i]=x;
1286 if (i<s)
1287 v[i+1]=plus_one;
1288 m.push_back(v);
1289 }
1290 return m;
1291 }
1292 static const char _JordanBlock_s []="JordanBlock";
1293 static define_unary_function_eval (__JordanBlock,&_JordanBlock,_JordanBlock_s);
1294 define_unary_function_ptr5( at_JordanBlock ,alias_at_JordanBlock,&__JordanBlock,0,true);
1295
_companion(const gen & g,GIAC_CONTEXT)1296 gen _companion(const gen & g,GIAC_CONTEXT){
1297 if ( g.type==_STRNG && g.subtype==-1) return g;
1298 vecteur v;
1299 if (g.type!=_VECT)
1300 return _companion(makesequence(g,vx_var),contextptr); // gentypeerr();
1301 if (g.subtype==_SEQ__VECT && g._VECTptr->size()==2){
1302 gen P=g._VECTptr->front();
1303 gen x=g._VECTptr->back();
1304 gen Px=_e2r(makevecteur(P,x),contextptr);
1305 if (Px.type==_FRAC)
1306 Px=inv(Px._FRACptr->den,contextptr)*Px._FRACptr->num;
1307 if (Px.type!=_VECT)
1308 return gensizeerr();
1309 v=*Px._VECTptr;
1310 }
1311 else
1312 v=*g._VECTptr;
1313 return companion(v);
1314 }
1315 static const char _companion_s []="companion";
1316 static define_unary_function_eval (__companion,&_companion,_companion_s);
1317 define_unary_function_ptr5( at_companion ,alias_at_companion,&__companion,0,true);
1318
_border(const gen & g,GIAC_CONTEXT)1319 gen _border(const gen & g,GIAC_CONTEXT){
1320 if ( g.type==_STRNG && g.subtype==-1) return g;
1321 if (g.type!=_VECT || g._VECTptr->size()!=2 || !ckmatrix(g._VECTptr->front()) || g._VECTptr->back().type!=_VECT)
1322 return gensizeerr();
1323 matrice m=*g._VECTptr->front()._VECTptr,v=*g._VECTptr->back()._VECTptr;
1324 if (m.size()!=v.size())
1325 return gendimerr();
1326 m=mtran(m);
1327 if (ckmatrix(v))
1328 m=mergevecteur(m,mtran(v));
1329 else
1330 m.push_back(v);
1331 return mtran(m);
1332 }
1333 static const char _border_s []="border";
1334 static define_unary_function_eval (__border,&_border,_border_s);
1335 define_unary_function_ptr5( at_border ,alias_at_border,&__border,0,true);
1336
1337 // pade(f(x),x,n,p) or pade(f(x),x,N(x),p)
1338 // Find P/Q such that P/Q=f(x) mod N(x) [1st case N(x)=x^n], deg(P)<p
1339 // Assume deg(f)<n
_pade(const gen & g,GIAC_CONTEXT)1340 gen _pade(const gen & g,GIAC_CONTEXT){
1341 if ( g.type==_STRNG && g.subtype==-1) return g;
1342 if (g.type!=_VECT || g._VECTptr->size()!=4)
1343 return gensizeerr();
1344 vecteur v =*g._VECTptr;
1345 if (v[1].type!=_IDNT || v[3].type!=_INT_)
1346 return gensizeerr();
1347 int n,p=v[3].val,ns;
1348 gen f=v[0],x=v[1],N=v[2];
1349 if (v[2].type==_INT_){
1350 n=v[2].val+1;
1351 N=pow(v[1],n);
1352 ns=n;
1353 }
1354 else {
1355 n=p;
1356 ns=_degree(makevecteur(N,x),contextptr).val;
1357 }
1358 if (p<=0 || n<=0)
1359 return gensizeerr();
1360 f=_series(makevecteur(f,x,zero,ns-1),contextptr);
1361 vecteur l1(lop(f,at_order_size));
1362 vecteur l2(l1.size(),zero);
1363 f=subst(f,l1,l2,false,contextptr);
1364 // Convert f and N to poly1
1365 vecteur l(1,x);
1366 lvar(f,l);
1367 lvar(N,l);
1368 int ls=int(l.size());
1369 gen ff(sym2r(f,l,contextptr));
1370 gen fn,fd;
1371 fxnd(ff,fn,fd);
1372 gen Nf(sym2r(N,l,contextptr));
1373 gen Nn,Nd;
1374 fxnd(Nf,Nn,Nd); // Note nd has no meaning
1375 vecteur fp;
1376 if (fn.type==_POLY)
1377 fp=polynome2poly1(*fn._POLYptr,1);
1378 else
1379 return gensizeerr();
1380 vecteur Np;
1381 if (Nn.type==_POLY)
1382 Np=polynome2poly1(*Nn._POLYptr,1);
1383 else
1384 return gensizeerr();
1385 n=int(Np.size())-1;
1386 // Check that fp is a poly of degree less than n
1387 if (n<1 || p>n || signed(fp.size())>n)
1388 return gendimerr();
1389 vecteur a,b;
1390 if (!egcd_pade(Np,fp,p,a,b,0))
1391 *logptr(contextptr) << gettext("Solution may be wrong since a and b are not prime together: ")+gen(a).print(contextptr)+","+gen(b).print(contextptr) << '\n';
1392 gen res=poly12polynome(a,1,ls);
1393 res=res/(fd*gen(poly12polynome(b,1,ls)));
1394 res=r2sym(res,l,contextptr);
1395 return res;
1396 }
1397 static const char _pade_s []="pade";
1398 static define_unary_function_eval (__pade,&_pade,_pade_s);
1399 define_unary_function_ptr5( at_pade ,alias_at_pade,&__pade,0,true);
1400
1401 static const char _interp_s []="interp";
1402 static define_unary_function_eval (__interp,&_lagrange,_interp_s);
1403 define_unary_function_ptr5( at_interp ,alias_at_interp,&__interp,0,true);
1404
lhsrhs(const gen & g,int i)1405 static gen lhsrhs(const gen & g,int i){
1406 if (!is_equal(g) && !g.is_symb_of_sommet(at_interval))
1407 return gensizeerr();
1408 gen & f=g._SYMBptr->feuille;
1409 if (f.type!=_VECT || f._VECTptr->size()!=2)
1410 return gensizeerr();
1411 return (*f._VECTptr)[i];
1412 }
_lhs(const gen & g,GIAC_CONTEXT)1413 gen _lhs(const gen & g,GIAC_CONTEXT){
1414 if ( g.type==_STRNG && g.subtype==-1) return g;
1415 return lhsrhs(g,0);
1416 }
1417 static const char _lhs_s []="lhs";
1418 static define_unary_function_eval (__lhs,&_lhs,_lhs_s);
1419 define_unary_function_ptr5( at_lhs ,alias_at_lhs,&__lhs,0,true);
1420
_rhs(const gen & g,GIAC_CONTEXT)1421 gen _rhs(const gen & g,GIAC_CONTEXT){
1422 if ( g.type==_STRNG && g.subtype==-1) return g;
1423 return lhsrhs(g,1);
1424 }
1425 static const char _rhs_s []="rhs";
1426 static define_unary_function_eval (__rhs,&_rhs,_rhs_s);
1427 define_unary_function_ptr5( at_rhs ,alias_at_rhs,&__rhs,0,true);
1428
1429 // Given [v_0 ... v_(2n-1)] (begin of the recurrence sequence)
1430 // return [b_n...b_0] such that b_n*v_{n+k}+...+b_0*v_k=0
1431 // Example [1,-1,3,3] -> [1,-3,-6]
1432 // -> the recurrence relation is v_{n+2}=3v_{n+1}+6v_n
_reverse_rsolve(const gen & g,GIAC_CONTEXT)1433 gen _reverse_rsolve(const gen & g,GIAC_CONTEXT){
1434 if ( g.type==_STRNG && g.subtype==-1) return g;
1435 if (g.type!=_VECT)
1436 return gensizeerr();
1437 vecteur v=reverse_rsolve(*g._VECTptr);
1438 return v/v.front();
1439 }
1440 static const char _reverse_rsolve_s []="reverse_rsolve";
1441 static define_unary_function_eval (__reverse_rsolve,&_reverse_rsolve,_reverse_rsolve_s);
1442 define_unary_function_ptr5( at_reverse_rsolve ,alias_at_reverse_rsolve,&__reverse_rsolve,0,true);
1443
1444 // Approx fft or exact if args=poly1,omega,n
fft(const gen & g_orig,int direct,GIAC_CONTEXT)1445 gen fft(const gen & g_orig,int direct,GIAC_CONTEXT){
1446 if (g_orig.type==_VECT && g_orig.subtype==_SEQ__VECT && g_orig._VECTptr->size()>=3 && g_orig._VECTptr->front().type==_VECT){
1447 vecteur & v =*g_orig._VECTptr->front()._VECTptr;
1448 int n=int(v.size());
1449 if (n<2)
1450 return gendimerr();
1451 gen omega=(*g_orig._VECTptr)[1];
1452 gen modulo=(*g_orig._VECTptr)[2];
1453 if (direct==1)
1454 omega=invmod(omega,modulo);
1455 if (omega.type==_INT_ && modulo.type==_INT_ && n==(1<<(sizeinbase2(n)-1))){
1456 if (debug_infolevel)
1457 CERR << CLOCK()*1e-6 << " fft start" << '\n';
1458 vector<int> A; A.reserve(n);
1459 for (int i=0;i<n;++i){
1460 if (v[i].type==_INT_)
1461 A.push_back(v[i].val);
1462 else
1463 A.push_back(smod(v[i],modulo).val);
1464 }
1465 int p=modulo.val;
1466 fft2(&A.front(),n,omega.val,p,g_orig._VECTptr->size()==3);
1467 gen r=vecteur(0);
1468 vecteur & res=*r._VECTptr;
1469 res.reserve(n);
1470 if (direct){
1471 for (int i=0;i<n;++i)
1472 res.push_back(A[i]);
1473 }
1474 else {
1475 longlong ninv=invmod(n,p);
1476 if (ninv<0) ninv += p;
1477 for (int i=0;i<n;++i)
1478 res.push_back((ninv*A[i])%p);
1479 }
1480 if (debug_infolevel)
1481 CERR << CLOCK()*1e-6 << "fft end" << '\n';
1482 return r;
1483 }
1484 vecteur w(n),res;
1485 iterateur it=w.begin(),itend=w.end();
1486 if (omega.type==_INT_ && modulo.type==_INT_){
1487 int Omega=omega.val,m=modulo.val;
1488 longlong omegan=1;
1489 for (;it!=itend;++it,omegan=(Omega*omegan) %m){
1490 *it=int(omegan);
1491 }
1492 }
1493 else {
1494 gen omegan=1;
1495 for (;it!=itend;++it,omegan=smod(omega*omegan,modulo)){
1496 *it=omegan;
1497 }
1498 }
1499 environment * env = new environment;
1500 env->moduloon=true;
1501 env->modulo=modulo;
1502 fft(v,w,res,env);
1503 delete env;
1504 if (direct==1)
1505 return res;
1506 else
1507 return smod(invmod(n,modulo)*res,modulo);
1508 }
1509 gen g=g_orig;
1510 if (g.type!=_VECT)
1511 return gensizeerr();
1512 int n=int(g._VECTptr->size());
1513 if (n<2)
1514 return gendimerr();
1515 vector< complex<double> > vd;
1516 if (convert(*g._VECTptr,vd,true)){
1517 if (debug_infolevel)
1518 CERR << CLOCK()*1e-6 << " fft start" << '\n';
1519 double theta=2.0*M_PI/n;
1520 if (direct) theta=-theta;
1521 bool done=false;
1522 if (n==(1<<(sizeinbase2(n)-1))){
1523 fft2(&vd.front(),n,theta);
1524 done=true;
1525 }
1526 #if 1 //ndef HAVE_LIBGSL
1527 if (!done){
1528 complex<double> w(std::cos(theta),std::sin(theta));
1529 vector< complex<double> > W(n),res(n);
1530 for (int i=0;i<n;++i){
1531 if ( i % 64==0)
1532 W[i]=complex<double>(std::cos(i*theta),std::sin(i*theta));
1533 else
1534 W[i]=W[i-1]*w;
1535 }
1536 fft(&vd[0],n,&W[0],n,&res[0]);
1537 done=true;
1538 }
1539 #endif
1540 if (done){
1541 gen r=vecteur(0);
1542 vecteur & v=*r._VECTptr;
1543 v.reserve(n);
1544 if (direct){
1545 for (int i=0;i<n;++i){
1546 v.push_back(vd[i]);
1547 }
1548 }
1549 else {
1550 double invn=1.0/n;
1551 for (int i=0;i<n;++i){
1552 v.push_back(gen(invn*vd[i].real(),invn*vd[i].imag()));
1553 }
1554 }
1555 if (debug_infolevel)
1556 CERR << CLOCK()*1e-6 << " fft end" << '\n';
1557 return r;
1558 }
1559 }
1560 g=evalf_double(g_orig,1,contextptr);
1561 if (debug_infolevel)
1562 CERR << CLOCK()*1e-6 << " fft start" << '\n';
1563 vecteur v =*g._VECTptr;
1564 #ifdef HAVE_LIBGSL
1565 if (direct && is_zero(im(v,contextptr))){
1566 double * data=new double[n];
1567 for (int i=0;i<n;++i){
1568 if (v[i].type!=_DOUBLE_){
1569 delete [] data;
1570 return gensizeerr(contextptr);
1571 }
1572 data[i]=v[i]._DOUBLE_val;
1573 }
1574 if (n==(1<<(sizeinbase2(n)-1))){
1575 gsl_fft_real_radix2_transform(data,1,n);
1576 v[0]=data[0];
1577 v[n/2]=data[n/2];
1578 for (int i=1;i<n/2;++i){
1579 v[i]=gen(data[i],data[n-i]);
1580 v[n-i]=conj(v[i],contextptr);
1581 }
1582 }
1583 else {
1584 gsl_fft_real_wavetable * wavetable = gsl_fft_real_wavetable_alloc (n);
1585 gsl_fft_real_workspace * workspace = gsl_fft_real_workspace_alloc (n);
1586 gsl_fft_real_transform (data, 1, n,wavetable,workspace);
1587 gsl_fft_real_wavetable_free (wavetable);
1588 gsl_fft_real_workspace_free (workspace);
1589 v[0]=data[0];
1590 int n2;
1591 if (n%2==0){
1592 v[n/2]=data[n-1];
1593 n2=n/2;
1594 }
1595 else
1596 n2=n/2+1;
1597 for (int i=1;i<n2;++i){
1598 v[i]=gen(data[2*i-1],data[2*i]);
1599 v[n-i]=conj(v[i],contextptr);
1600 }
1601 }
1602 delete [] data;
1603 if (debug_infolevel)
1604 CERR << CLOCK()*1e-6 << " fft end" << '\n';
1605 return v;
1606 }
1607 // Could be improved by keeping the wavetable
1608 double * data=new double[2*n];
1609 gen gr,gi;
1610 for (int i=0;i<n;++i){
1611 reim(v[i],gr,gi,contextptr);
1612 gi=evalf_double(gi,1,contextptr);
1613 if (gr.type!=_DOUBLE_ || gi.type!=_DOUBLE_){
1614 delete [] data;
1615 return gensizeerr(contextptr);
1616 }
1617 data[2*i]=gr._DOUBLE_val;
1618 data[2*i+1]=gi._DOUBLE_val;
1619 }
1620 if (n==(1<<(sizeinbase2(n)-1))){
1621 if (direct)
1622 gsl_fft_complex_radix2_forward(data,1,n);
1623 else
1624 gsl_fft_complex_radix2_backward(data,1,n);
1625 }
1626 else {
1627 gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (n);
1628 gsl_fft_complex_workspace * workspace=gsl_fft_complex_workspace_alloc (n);
1629 if (direct)
1630 gsl_fft_complex_forward (data, 1, n,wavetable,workspace);
1631 else
1632 gsl_fft_complex_backward (data, 1, n,wavetable,workspace);
1633 gsl_fft_complex_wavetable_free (wavetable);
1634 gsl_fft_complex_workspace_free (workspace);
1635 }
1636 if (direct){
1637 for (int i=0;i<n;++i){
1638 v[i]=gen(data[2*i],data[2*i+1]);
1639 }
1640 }
1641 else {
1642 for (int i=0;i<n;++i){
1643 v[i]=gen(data[2*i],data[2*i+1])/n;
1644 }
1645 }
1646 delete [] data;
1647 if (debug_infolevel)
1648 CERR << CLOCK()*1e-6 << " fft end" << '\n';
1649 return v;
1650 #endif
1651 /*
1652 unsigned m=gen(n).bindigits()-1;
1653 if (n!=1<<m)
1654 return gensizeerr(gettext("Size is not a power of 2 ")+print_INT_(n));
1655 */
1656 vecteur res;
1657 double theta;
1658 vecteur w(n);
1659 iterateur it=w.begin(),itend=w.end();
1660 for (int i=0;it!=itend;++it,++i){
1661 theta = (2*M_PI*i)/n;
1662 if (direct==1)
1663 *it = gen(std::cos(theta),-std::sin(theta));
1664 else
1665 *it = gen(std::cos(theta),std::sin(theta));
1666 }
1667 fft(v,w,res,0);
1668 if (direct==1)
1669 return res;
1670 else
1671 return gen(res)/n;
1672 }
1673
_fft(const gen & g,GIAC_CONTEXT)1674 gen _fft(const gen & g,GIAC_CONTEXT){
1675 if ( g.type==_STRNG && g.subtype==-1) return g;
1676 return fft(g,1,contextptr);
1677 }
1678 static const char _fft_s []="fft";
1679 static define_unary_function_eval (__fft,&_fft,_fft_s);
1680 define_unary_function_ptr5( at_fft ,alias_at_fft,&__fft,0,true);
1681
_ifft(const gen & g,GIAC_CONTEXT)1682 gen _ifft(const gen & g,GIAC_CONTEXT){
1683 if ( g.type==_STRNG && g.subtype==-1) return g;
1684 return fft(g,0,contextptr);
1685 }
1686 static const char _ifft_s []="ifft";
1687 static define_unary_function_eval (__ifft,&_ifft,_ifft_s);
1688 define_unary_function_ptr5( at_ifft ,alias_at_ifft,&__ifft,0,true);
1689
_Resultant(const gen & g,GIAC_CONTEXT)1690 gen _Resultant(const gen & g,GIAC_CONTEXT){
1691 if ( g.type==_STRNG && g.subtype==-1) return g;
1692 return symbolic(at_resultant,g);
1693 }
1694 static const char _Resultant_s []="Resultant";
1695 static define_unary_function_eval (__Resultant,&_Resultant,_Resultant_s);
1696 define_unary_function_ptr5( at_Resultant ,alias_at_Resultant,&__Resultant,0,true);
1697
_Nullspace(const gen & g,GIAC_CONTEXT)1698 gen _Nullspace(const gen & g,GIAC_CONTEXT){
1699 if ( g.type==_STRNG && g.subtype==-1) return g;
1700 return symbolic(at_nullspace,g);
1701 }
1702 static const char _Nullspace_s []="Nullspace";
1703 static define_unary_function_eval (__Nullspace,&_Nullspace,_Nullspace_s);
1704 define_unary_function_ptr5( at_Nullspace ,alias_at_Nullspace,&__Nullspace,0,true);
1705
_assign(const gen & g,GIAC_CONTEXT)1706 gen _assign(const gen & g,GIAC_CONTEXT){
1707 if ( g.type==_STRNG && g.subtype==-1) return g;
1708 if (g.type==_VECT && g.subtype==_SEQ__VECT && g._VECTptr->size()==2){
1709 return sto(g._VECTptr->back(),g._VECTptr->front(),contextptr);
1710 }
1711 if (is_equal(g)){
1712 gen & f=g._SYMBptr->feuille;
1713 if (f.type==_VECT && f._VECTptr->size()==2)
1714 return sto(f._VECTptr->back(),f._VECTptr->front(),contextptr);
1715 }
1716 if (g.type!=_VECT)
1717 return gensizeerr();
1718 return apply(g,_assign,contextptr);
1719 }
1720 static const char _assign_s []="assign";
1721 static define_unary_function_eval (__assign,&_assign,_assign_s);
1722 define_unary_function_ptr5( at_assign ,alias_at_assign,&__assign,0,true);
1723
1724
_implicitplot3d(const gen & g,GIAC_CONTEXT)1725 gen _implicitplot3d(const gen & g,GIAC_CONTEXT){
1726 if ( g.type==_STRNG && g.subtype==-1) return g;
1727 return _plotimplicit(g,contextptr);
1728 }
1729 static const char _implicitplot3d_s []="implicitplot3d";
1730 static define_unary_function_eval (__implicitplot3d,&_implicitplot3d,_implicitplot3d_s);
1731 define_unary_function_ptr5( at_implicitplot3d ,alias_at_implicitplot3d,&__implicitplot3d,0,true);
1732
convert_double_int(vecteur & v)1733 void convert_double_int(vecteur & v){
1734 for (unsigned i=0;i<v.size();++i){
1735 if (v[i].type==_DOUBLE_)
1736 v[i]=int(v[i]._DOUBLE_val+.5);
1737 if (v[i].type==_VECT){
1738 vecteur w=*v[i]._VECTptr;
1739 convert_double_int(w);
1740 v[i]=w;
1741 }
1742 }
1743 }
1744
read_audio(vecteur & v,int & channels,int & sample_rate,int & bits_per_sample,unsigned int & data_size)1745 bool read_audio(vecteur & v,int & channels,int & sample_rate,int & bits_per_sample,unsigned int & data_size){
1746 convert_double_int(v);
1747 if (v.size()>1 && v[1].type!=_VECT)
1748 v=makevecteur(1,v);
1749 gen g=v[0];
1750 if (ckmatrix(g)){
1751 if (v.size()==2 && v[1].type==_INT_){
1752 sample_rate=v[1].val;
1753 v=*g._VECTptr;
1754 v.insert(v.begin(),makevecteur(1,16,sample_rate));
1755 g=v[0];
1756 }
1757 else {
1758 if (v.size()==3 && v[1].type==_INT_ && v[2].type==_INT_){
1759 sample_rate=v[1].val;
1760 bits_per_sample=v[2].val;
1761 v=*g._VECTptr;
1762 v.insert(v.begin(),makevecteur(1,bits_per_sample,sample_rate));
1763 g=v[0];
1764 }
1765 }
1766 }
1767 // g=channels,bits_per_sample,sample_rate,data_size,
1768 if (g.type==_INT_)
1769 g=makevecteur(g,16);
1770 if (g.type!=_VECT || g._VECTptr->empty())
1771 return false;
1772 vecteur w=*g._VECTptr;
1773 if (w.size()==1)
1774 w.push_back(16);
1775 int ws=int(w.size());
1776 if (w[0].type!=_INT_ || w[1].type!=_INT_)
1777 return false;
1778 channels=w[0].val;
1779 if (channels<=0 || channels>int(v.size())){
1780 channels=int(v.size());
1781 if (v.front().type!=_VECT || !is_integer_vecteur(*v.front()._VECTptr))
1782 return false;
1783 data_size=unsigned(v.front()._VECTptr->size());
1784 for (int i=1;i<channels;++i){
1785 if (v[i].type!=_VECT || !is_integer_vecteur(*v[i]._VECTptr))
1786 return false;
1787 if (data_size>v[i]._VECTptr->size())
1788 data_size=unsigned(v[i]._VECTptr->size());
1789 }
1790 w=makevecteur(channels,16,44100); ws=3;
1791 g=w;
1792 v.insert(v.begin(),g);
1793 }
1794 if (channels<=0 || channels>4 || int(v.size())<=channels)
1795 return false;
1796 bits_per_sample=(w[1].val/8)*8;
1797 if (bits_per_sample<=0)
1798 return false;
1799 if (ws>=3)
1800 sample_rate=w[2].val;
1801 else
1802 sample_rate=44100;
1803 if (sample_rate<1)
1804 return false;
1805 if (ws>=4 && w[3].type==_INT_)
1806 data_size=w[3].val;
1807 else
1808 data_size=RAND_MAX;
1809 for (int i=1;i<=channels;++i){
1810 if (v[i].type!=_VECT || !is_integer_vecteur(*v[i]._VECTptr))
1811 return false;
1812 if (data_size>v[i]._VECTptr->size())
1813 data_size=unsigned(v[i]._VECTptr->size());
1814 }
1815 return true;
1816 }
1817
1818 #ifdef HAVE_LIBAO
1819 #include <ao/ao.h>
1820 typedef unsigned short aou16;
1821 typedef unsigned int aou32;
_playsnd(const gen & args,GIAC_CONTEXT)1822 gen _playsnd(const gen & args,GIAC_CONTEXT){
1823 if (args.type==_STRNG){
1824 if (args.subtype==-1) return args;
1825 return _playsnd(_readwav(args,contextptr),contextptr);
1826 }
1827 ao_device *device=0;
1828 ao_sample_format format;
1829 int default_driver;
1830 ao_initialize();
1831 default_driver = ao_default_driver_id();
1832 memset(&format, 0, sizeof(format));
1833 format.bits = 16;
1834 format.channels = 2;
1835 format.rate = 44100;
1836 format.byte_format = AO_FMT_LITTLE;
1837 unsigned int data_size=0;
1838 vecteur v;
1839 if (args.type==_VECT && !args._VECTptr->empty()){
1840 // set format
1841 v=*args._VECTptr;
1842 if (!read_audio(v,format.channels,format.rate,format.bits,data_size))
1843 return gensizeerr(gettext("Invalid sound data"));
1844 }
1845 if (data_size){
1846 *logptr(contextptr) << gettext("Using sound parameters: channels, rate, bits, records ") << format.channels << "," << format.rate << "," << format.bits << "," << data_size << '\n';
1847 device = ao_open_live(default_driver, &format, NULL /* no options */);
1848 if (device == NULL)
1849 return gensizeerr(gettext("Error opening audio device."));
1850 unsigned n=data_size*format.channels*format.bits/8;
1851 char * buffer=(char *)malloc(n*sizeof(char));
1852 // aou16 * bufshort=(aou16 *) buffer;
1853 aou32 * bufint=(aou32 *) buffer;
1854 if (buffer){
1855 // copy data from v into buffer and play it
1856 unsigned c=format.channels,b=format.bits/8;
1857 for (unsigned i=0;i<data_size;++i){
1858 for (unsigned j=0;j<c;++j){
1859 unsigned u=(*v[j+1]._VECTptr)[i].val;
1860 if (b==1)
1861 buffer[i*c+j]=u;
1862 if (b==2){
1863 buffer[2*(i*c+j)]=u & 0xff;
1864 buffer[2*(i*c+j)+1]=(u>>8) & 0xff;
1865 }
1866 if (b==4)
1867 bufint[i*c+j]=u;
1868 }
1869 }
1870 ao_play(device, buffer, n);
1871 }
1872 }
1873 ao_close(device);
1874 ao_shutdown();
1875 return 1;
1876 }
1877 static const char _playsnd_s []="playsnd";
1878 static define_unary_function_eval (__playsnd,&_playsnd,_playsnd_s);
1879 define_unary_function_ptr5( at_playsnd ,alias_at_playsnd,&__playsnd,0,true);
1880 #else
1881 #if !defined GIAC_GGB && defined EMCC // must have EM_ASM code javascript inlined (emscripten 1.30.4 at least?)
1882 #include <emscripten.h>
_playsnd(const gen & args,GIAC_CONTEXT)1883 gen _playsnd(const gen & args,GIAC_CONTEXT){
1884 if (args.type==_STRNG){
1885 if (args.subtype==-1) return args;
1886 return _playsnd(_readwav(args,contextptr),contextptr);
1887 }
1888
1889 int nbits = 16;
1890 int nchannels = 2;
1891 int nrate = 44100;
1892 unsigned int data_size=0;
1893 vecteur v;
1894 if (args.type==_VECT && !args._VECTptr->empty()){
1895 // set format
1896 v=*args._VECTptr;
1897 if (!read_audio(v,nchannels,nrate,nbits,data_size))
1898 return gensizeerr(gettext("Invalid sound data"));
1899 }
1900 if (data_size){
1901 *logptr(contextptr) << gettext("Using sound parameters: channels, rate, bits, records ") << nchannels << "," << nrate << "," << data_size << '\n';
1902 unsigned nDataBytes=data_size*nchannels*sizeof(float);
1903 // copy data from v into buffer and play it
1904 unsigned b=nbits/8;
1905 float * ptr = (float *) malloc(nDataBytes);
1906 for (unsigned j=0;j<nchannels;++j){
1907 vecteur & w=(*v[j+1]._VECTptr);
1908 COUT << "channel " << j << '\n';
1909 for (unsigned i=0;i<data_size;++i){
1910 unsigned u=w[i].val;
1911 double ud=0;
1912 if (b==1)
1913 ud=u/128.0-1;
1914 if (b==2)
1915 ud=u/32768.0-1;
1916 if (b==4)
1917 ud=u/2147483648.0-1;
1918 ptr[j*data_size+i]=ud;
1919 }
1920 }
1921 COUT << "playing" << '\n';
1922 EM_ASM_ARGS({
1923 var nchannels;
1924 var nDataBytes;
1925 var nrate;
1926 var ptr;
1927 var data_size;
1928 nchannels=$0;nDataBytes=$1;nrate=$2;ptr=$3;
1929 data_size=nDataBytes/4/nchannels;
1930 var audioCtx = new (window.AudioContext || window.webkitAudioContext)();
1931 var SoundArrayBuffer = audioCtx.createBuffer(nchannels, nDataBytes, audioCtx.sampleRate);
1932 var dataHeap = new Uint8Array(Module.HEAPU8.buffer, ptr, nDataBytes);
1933 var result = new Float32Array(dataHeap.buffer, dataHeap.byteOffset, nDataBytes/4);
1934 var j;
1935 var i;
1936 for (j=0;j<nchannels;j++){
1937 var v=SoundArrayBuffer.getChannelData(j);
1938 for (i=0;i<data_size;++i)
1939 v[i]=result[j*data_size+i];
1940 }
1941 var source = audioCtx.createBufferSource();
1942 // set the buffer in the AudioBufferSourceNode
1943 source.buffer = SoundArrayBuffer;
1944 // connect the AudioBufferSourceNode to the
1945 // destination so we can hear the sound
1946 source.connect(audioCtx.destination);
1947 // start the source playing
1948 source.start();
1949 },nchannels,nDataBytes,nrate,ptr);
1950 free(ptr);
1951 }
1952 return 1;
1953 }
1954 #else
_playsnd(const gen & args,GIAC_CONTEXT)1955 gen _playsnd(const gen & args,GIAC_CONTEXT){
1956 return gensizeerr("Sorry! libao is not present on system");
1957 }
1958 #endif
1959 static const char _playsnd_s []="playsnd";
1960 static define_unary_function_eval (__playsnd,&_playsnd,_playsnd_s);
1961 define_unary_function_ptr5( at_playsnd ,alias_at_playsnd,&__playsnd,0,true);
1962 #endif
1963
_soundsec(const gen & args,GIAC_CONTEXT)1964 gen _soundsec(const gen & args,GIAC_CONTEXT){
1965 // nseconds [,rate]
1966 gen n,rate=44100;
1967 if (args.type==_VECT && args._VECTptr->size()==2){
1968 n=args._VECTptr->front();
1969 rate=args._VECTptr->back();
1970 }
1971 else
1972 n=args;
1973 n=evalf_double(n,1,contextptr);
1974 if (n.type!=_DOUBLE_ || n._DOUBLE_val<=0 || rate.type!=_INT_ || rate.val < 1 )
1975 return gensizeerr(gettext("Invalid sound parameters"));
1976 double r=evalf_double(rate,1,contextptr)._DOUBLE_val;
1977 double nr=r*n._DOUBLE_val;
1978 if (nr>LIST_SIZE_LIMIT)
1979 return gensizeerr("Too many records");
1980 vecteur v;
1981 v.reserve(int(nr));
1982 for (int i=0;i<nr;++i){
1983 v.push_back(double(i)/r);
1984 }
1985 return v;
1986 }
1987 static const char _soundsec_s []="soundsec";
1988 static define_unary_function_eval (__soundsec,&_soundsec,_soundsec_s);
1989 define_unary_function_ptr5( at_soundsec ,alias_at_soundsec,&__soundsec,0,true);
1990
1991 /*
1992 gen _beep(const gen & args,GIAC_CONTEXT){
1993 }
1994 static const char _beep_s []="beep";
1995 static define_unary_function_eval (__beep,&_beep,_beep_s);
1996 define_unary_function_ptr5( at_beep ,alias_at_beep,&__beep,0,true);
1997 */
1998
1999 #ifdef RTOS_THREADX
_readwav(const gen & args,GIAC_CONTEXT)2000 gen _readwav(const gen & args,GIAC_CONTEXT){
2001 return undef;
2002 }
2003 static const char _readwav_s []="readwav";
2004 static define_unary_function_eval (__readwav,&_readwav,_readwav_s);
2005 define_unary_function_ptr5( at_readwav ,alias_at_readwav,&__readwav,0,true);
2006
2007 static const char _writewav_s []="writewav";
2008 static define_unary_function_eval (__writewav,&_readwav,_writewav_s);
2009 define_unary_function_ptr5( at_writewav ,alias_at_writewav,&__writewav,0,true);
2010
2011 #else // RTOS_THREADX
2012
2013 // http://ccrma.stanford.edu/courses/422/projects/WaveFormat/
in_readwav(FILE * f,gen & g)2014 static bool in_readwav(FILE * f,gen & g){
2015 unsigned char c,channels;
2016 unsigned int u,s,sample_rate,byte_rate,block_align=0,bits_per_sample=0,data_size=0;
2017 // Header
2018 if (fread(&u,4,1,f)!=1 || u!=0x46464952) // "RIFF"
2019 return false;
2020 if (fread(&s,4,1,f)!=1)
2021 return false;
2022 if (fread(&u,4,1,f)!=1 || u!=0x45564157) // "WAVE"
2023 return false;
2024 if (fread(&u,4,1,f)!=1 || u!=0x20746d66) // "fmt "
2025 return false;
2026 if (fread(&u,4,1,f)!=1 || u!=0x10) // 16 for PCM
2027 return false;
2028 c=fgetc(f);
2029 if (c!=1)
2030 return false;
2031 c=fgetc(f);
2032 if (c!=0)
2033 return false;
2034 channels=fgetc(f);
2035 c=fgetc(f);
2036 if (c!=0)
2037 return false;
2038 if (fread(&sample_rate,4,1,f)!=1)
2039 return false;
2040 if (fread(&byte_rate,4,1,f)!=1)
2041 return false;
2042 block_align=fgetc(f);
2043 block_align=block_align+(fgetc(f)<<8);
2044 bits_per_sample=fgetc(f); // 8 or 16
2045 bits_per_sample=bits_per_sample+(fgetc(f)<<8);
2046 bits_per_sample /= 8;
2047 if (fread(&u,4,1,f)!=1 || u!=0x61746164) // "data"
2048 return false;
2049 if (fread(&data_size,4,1,f)!=1)
2050 return false;
2051 int n=data_size;
2052 // data_size=bits_per_sample/8 * num_samples * channels
2053 vecteur v(channels+1);
2054 v[0]=makevecteur(int(channels),int(bits_per_sample*8),int(sample_rate),int(data_size));
2055 g=v;
2056 vecteur & w=*g._VECTptr;
2057 for (int i=1;i<=channels;++i){
2058 w[i]=vecteur(0);
2059 w[i]._VECTptr->reserve(n/(channels*bits_per_sample));
2060 }
2061 while (n>0 && !feof(f)){
2062 for (int i=1;i<=channels;++i){
2063 u=0;
2064 if (fread(&u,bits_per_sample,1,f)!=1)
2065 return false;
2066 // for (int j=0;j<bits_per_sample;++j)
2067 // u += fgetc(f) << (j*8) ;
2068 n -= bits_per_sample;
2069 if (feof(f))
2070 break;
2071 w[i]._VECTptr->push_back(int(u));
2072 if (n<=0)
2073 break;
2074 }
2075 }
2076 return true;
2077 }
2078
readwav(const string & s,gen & g)2079 static bool readwav(const string & s,gen & g){
2080 FILE * f = fopen(s.c_str(),"r");
2081 if (!f)
2082 return false;
2083 bool res=in_readwav(f,g);
2084 fclose(f);
2085 return res;
2086 }
2087
_readwav(const gen & g,GIAC_CONTEXT)2088 gen _readwav(const gen & g,GIAC_CONTEXT){
2089 if ( g.type==_STRNG && g.subtype==-1) return g;
2090 if (g.type!=_STRNG)
2091 return gensizeerr();
2092 gen res;
2093 bool ok= readwav(*g._STRNGptr,res);
2094 if (!ok)
2095 return gensizeerr(gettext("File not found or unrecognized wav file format"));
2096 return res;
2097 }
2098 static const char _readwav_s []="readwav";
2099 static define_unary_function_eval (__readwav,&_readwav,_readwav_s);
2100 define_unary_function_ptr5( at_readwav ,alias_at_readwav,&__readwav,0,true);
2101
in_writewav(FILE * f,const vecteur & v_)2102 static bool in_writewav(FILE * f,const vecteur & v_){
2103 if (v_.empty())
2104 return false;
2105 vecteur v(v_);
2106 int channels,sample_rate=44100,bits_per_sample=0;
2107 unsigned int u,byte_rate,block_align=0,data_size=1U<<31;
2108 if (!read_audio(v,channels,sample_rate,bits_per_sample,data_size))
2109 return false;
2110 u=0x46464952;
2111 if (fwrite(&u,4,1,f)!=1)
2112 return false;
2113 u=36+data_size*bits_per_sample/8*channels;
2114 if (fwrite(&u,4,1,f)!=1)
2115 return false;
2116 u=0x45564157;
2117 if (fwrite(&u,4,1,f)!=1)
2118 return false;
2119 u=0x20746d66;
2120 if (fwrite(&u,4,1,f)!=1)
2121 return false;
2122 u=0x10;
2123 if (fwrite(&u,4,1,f)!=1)
2124 return false;
2125 fputc(1,f);
2126 fputc(0,f);
2127 fputc(channels,f);
2128 fputc(0,f);
2129 if (fwrite(&sample_rate,4,1,f)!=1)
2130 return false;
2131 byte_rate = sample_rate * channels * bits_per_sample/8;
2132 if (fwrite(&byte_rate,4,1,f)!=1)
2133 return false;
2134 block_align=channels * bits_per_sample/8;
2135 if (fwrite(&block_align,2,1,f)!=1)
2136 return false;
2137 if (fwrite(&bits_per_sample,2,1,f)!=1)
2138 return false;
2139 u=0x61746164; // "data"
2140 if (fwrite(&u,4,1,f)!=1)
2141 return false;
2142 u= bits_per_sample/8* data_size*channels; // should be data_size
2143 if (fwrite(&u,4,1,f)!=1)
2144 return false;
2145 // write data
2146 u /= channels;
2147 bits_per_sample /= 8;
2148 unsigned n = u/ bits_per_sample;
2149 for (unsigned i=0;i<n;++i){
2150 for (int j=1;j<=channels;++j){
2151 u=(*v[j]._VECTptr)[i].val;
2152 if (fwrite(&u,bits_per_sample,1,f)!=1)
2153 return false;
2154 }
2155 }
2156 return true;
2157 }
2158
writewav(const string & s,const vecteur & v)2159 static bool writewav(const string & s,const vecteur & v){
2160 FILE * f = fopen(s.c_str(),"w");
2161 if (!f)
2162 return false;
2163 bool res=in_writewav(f,v);
2164 fclose(f);
2165 return res;
2166 }
2167
_writewav(const gen & g,GIAC_CONTEXT)2168 gen _writewav(const gen & g,GIAC_CONTEXT){
2169 if ( g.type==_STRNG && g.subtype==-1) return g;
2170 if (g.type!=_VECT || g._VECTptr->size()!=2 || g._VECTptr->front().type!=_STRNG || g._VECTptr->back().type!=_VECT)
2171 return gensizeerr();
2172 bool ok= writewav(*g._VECTptr->front()._STRNGptr,*g._VECTptr->back()._VECTptr);
2173 if (!ok)
2174 return gensizeerr(gettext("Unable to open file or unable to code data"));
2175 return 1;
2176 }
2177 static const char _writewav_s []="writewav";
2178 static define_unary_function_eval (__writewav,&_writewav,_writewav_s);
2179 define_unary_function_ptr5( at_writewav ,alias_at_writewav,&__writewav,0,true);
2180
2181 #endif // RTOS_THREADX
2182
animate2d3d(const gen & g,bool dim3,GIAC_CONTEXT)2183 static gen animate2d3d(const gen & g,bool dim3,GIAC_CONTEXT){
2184 int s=0,frames=10;
2185 if (g.type!=_VECT || (s=int(g._VECTptr->size()))<3)
2186 return gensizeerr();
2187 vecteur v = *g._VECTptr;
2188 gen t;
2189 double tmin,tmax;
2190 if (!readrange(v[2],gnuplot_tmin,gnuplot_tmax,t,tmin,tmax,contextptr))
2191 return gensizeerr();
2192 // find a frame argument, otherwise use 10
2193 for (int i=3;i<s;++i){
2194 if (is_equal(v[i])){
2195 gen & f = v[i]._SYMBptr->feuille;
2196 if (f.type==_VECT && f._VECTptr->size()==2 && f._VECTptr->front().type==_INT_ && f._VECTptr->front().val==_FRAMES && f._VECTptr->back().type==_INT_)
2197 frames=f._VECTptr->back().val;
2198 }
2199 }
2200 if (frames==0)
2201 return gensizeerr();
2202 v.erase(v.begin()+2);
2203 gen h=symbolic(dim3?at_plotfunc:at_plot,gen(v,g.subtype));
2204 v=makevecteur(h,t,tmin,tmax,(tmax-tmin)/frames);
2205 gen tmpseq=seqprod(v,0,contextptr);
2206 if (is_undef(tmpseq)) return tmpseq;
2207 return symbolic(at_animation,tmpseq);
2208 }
2209
_animate(const gen & g,GIAC_CONTEXT)2210 gen _animate(const gen & g,GIAC_CONTEXT){
2211 if ( g.type==_STRNG && g.subtype==-1) return g;
2212 return animate2d3d(g,false,contextptr);
2213 }
2214 static const char _animate_s []="animate";
2215 static define_unary_function_eval (__animate,&_animate,_animate_s);
2216 define_unary_function_ptr5( at_animate ,alias_at_animate,&__animate,0,true);
2217
_animate3d(const gen & g,GIAC_CONTEXT)2218 gen _animate3d(const gen & g,GIAC_CONTEXT){
2219 if ( g.type==_STRNG && g.subtype==-1) return g;
2220 return animate2d3d(g,true,contextptr);
2221 }
2222 static const char _animate3d_s []="animate3d";
2223 static define_unary_function_eval (__animate3d,&_animate3d,_animate3d_s);
2224 define_unary_function_ptr5( at_animate3d ,alias_at_animate3d,&__animate3d,0,true);
2225
_even(const gen & g_,GIAC_CONTEXT)2226 gen _even(const gen & g_,GIAC_CONTEXT){
2227 gen g(g_);
2228 if ( g.type==_STRNG && g.subtype==-1) return g;
2229 if (!is_integral(g)) return gentypeerr(contextptr);
2230 return is_zero(smod(g,2));
2231 }
2232 static const char _even_s []="even";
2233 static define_unary_function_eval (__even,&_even,_even_s);
2234 define_unary_function_ptr5( at_even ,alias_at_even,&__even,0,true);
2235
_odd(const gen & g_,GIAC_CONTEXT)2236 gen _odd(const gen & g_,GIAC_CONTEXT){
2237 gen g(g_);
2238 if ( g.type==_STRNG && g.subtype==-1) return g;
2239 if (!is_integral(g)) return gentypeerr(contextptr);
2240 return !is_zero(smod(g,2));
2241 }
2242 static const char _odd_s []="odd";
2243 static define_unary_function_eval (__odd,&_odd,_odd_s);
2244 define_unary_function_ptr5( at_odd ,alias_at_odd,&__odd,0,true);
2245
2246 #ifdef HAVE_LIBPNG
write_png(const char * file_name,void * rows_,int w,int h,int colortype,int bitdepth)2247 int write_png(const char *file_name, void *rows_, int w, int h, int colortype, int bitdepth){
2248 png_bytep * rows=(png_bytep *) rows_;
2249 png_structp png_ptr;
2250 png_infop info_ptr;
2251 FILE *fp = fopen(file_name, "wb");
2252 const char *doing = "open for writing";
2253 if (!(fp = fopen(file_name, "wb"))) goto fail;
2254 doing = "create png write struct";
2255 if (!(png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL,
2256 NULL, NULL))) goto fail;
2257 doing = "create png info struct";
2258 if (!(info_ptr = png_create_info_struct(png_ptr))) goto fail;
2259 if (setjmp(png_jmpbuf(png_ptr))) goto fail;
2260 doing = "init IO";
2261 png_init_io(png_ptr, fp);
2262 doing = "write header";
2263 png_set_IHDR(png_ptr, info_ptr, w, h, bitdepth, colortype,
2264 PNG_INTERLACE_NONE,
2265 PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
2266 doing = "write info";
2267 png_write_info(png_ptr, info_ptr);
2268 doing = "write image";
2269 png_write_image(png_ptr, rows);
2270 doing = "write end";
2271 png_write_end(png_ptr, NULL);
2272 fclose(fp);
2273 return 0;
2274 fail:
2275 printf("Write_png: could not %s\n", doing);
2276 return -1;
2277 }
2278
writergb(const string & filename,const vecteur & v)2279 static bool writergb(const string & filename,const vecteur & v){
2280 if (v.size()==2 && ckmatrix(v[1])){
2281 int w,h;
2282 mdims(*v[1]._VECTptr,h,w);
2283 unsigned i=0,taille=w*h;
2284 unsigned char * screenbuf=new unsigned char[w*h];
2285 unsigned rowbytes = w;
2286 // fill screenbuf with data from v[] order 1,2,4,3
2287 for (;i<taille;++i){
2288 vecteur & wnv =*v[1]._VECTptr;
2289 int nr=i/rowbytes;
2290 vecteur & row =*wnv[nr]._VECTptr;
2291 screenbuf[i]=row[i%w].val;
2292 }
2293 unsigned char *rows[h];
2294 for (i = 0; int(i) < h; i++) {
2295 rows[i] = &screenbuf[i*w]; // &screenbuf[(h - i - 1)*4*w];
2296 }
2297 int res=write_png(filename.c_str(), rows, w, h, PNG_COLOR_TYPE_GRAY, 8);
2298 delete screenbuf;
2299 return res+1;
2300 }
2301 if (v.size()!=5 || !ckmatrix(v[1]) || !ckmatrix(v[2])
2302 || !ckmatrix(v[3]) || !ckmatrix(v[4]))
2303 return false;
2304 int w,h,w2,w3,w4,h2,h3,h4;
2305 mdims(*v[1]._VECTptr,h,w);
2306 mdims(*v[2]._VECTptr,h2,w2);
2307 mdims(*v[3]._VECTptr,h3,w3);
2308 mdims(*v[4]._VECTptr,h4,w4);
2309 if (w!=w2 || w!=w3 || w!=w4 || h!=h2 || h!=h3 || h!=h4)
2310 return false;
2311 unsigned i=0,taille=w*h*4;
2312 unsigned char * screenbuf=new unsigned char[w*h*4];
2313 unsigned rowbytes = w*4;
2314 // fill screenbuf with data from v[] order 1,2,4,3
2315 for (;i<taille;++i){
2316 int nv,nr;
2317 switch (i%4){
2318 case 0: nv=1; break;
2319 case 1: nv=2; break;
2320 case 2: nv=4; break;
2321 case 3: nv=3; break;
2322 }
2323 vecteur & wnv =*v[nv]._VECTptr;
2324 nr=i/rowbytes;
2325 vecteur & row =*wnv[nr]._VECTptr;
2326 screenbuf[i]=row[(i/4)%w].val;
2327 }
2328 unsigned char *rows[h];
2329 for (i = 0; int(i) < h; i++) {
2330 rows[i] = &screenbuf[i*4*w]; // &screenbuf[(h - i - 1)*4*w];
2331 }
2332 int res=write_png(filename.c_str(), rows, w, h, PNG_COLOR_TYPE_RGBA, 8);
2333 delete screenbuf;
2334 return res+1;
2335 }
2336 #else // LIBPNG
write_png(const char * file_name,void * rows,int w,int h,int colortype,int bitdepth)2337 int write_png(const char *file_name, void *rows, int w, int h, int colortype, int bitdepth){
2338 return -1;
2339 }
2340
writergb(const string & s,const vecteur & w)2341 static bool writergb(const string & s,const vecteur & w){
2342 return false;
2343 }
2344 #endif // LIBPNG
2345
_writergb(const gen & g,GIAC_CONTEXT)2346 gen _writergb(const gen & g,GIAC_CONTEXT){
2347 if ( g.type==_STRNG && g.subtype==-1) return g;
2348 vecteur v(gen2vecteur(g));
2349 if (ckmatrix(v[1])){
2350 int l,c;
2351 mdims(*v[1]._VECTptr,l,c);
2352 vecteur w(1,makevecteur(v.size()==2?1:4,c,l));
2353 for (unsigned i=1;i<v.size();++i){
2354 if (i==3 && v.size()==4)
2355 w.push_back(vecteur(l,vecteur(c,255)));
2356 w.push_back(v[i]);
2357 }
2358 v=makevecteur(v[0],w);
2359 }
2360 if (v.size()!=2 || v[0].type!=_STRNG || v[1].type!=_VECT)
2361 return gensizeerr();
2362 vecteur w=*v[1]._VECTptr;
2363 // w[0]==[d,w,h], w[1..4]=data
2364 bool ok= writergb(*v[0]._STRNGptr,w);
2365 if (!ok)
2366 return gensizeerr(gettext("File not found or unrecognized image file format"));
2367 gen tmp(_GL_TEXTURE);
2368 tmp.subtype=_INT_PLOT;
2369 vecteur arg=makevecteur(w[0][2],0,w[0][1]/w[0][2],
2370 symb_equal(tmp,makevecteur(v[0],undef)));
2371 tmp=gen(arg,_SEQ__VECT);
2372 return _rectangle(tmp,contextptr);
2373 }
2374 static const char _writergb_s []="writergb";
2375 static define_unary_function_eval (__writergb,&_writergb,_writergb_s);
2376 define_unary_function_ptr5( at_writergb ,alias_at_writergb,&__writergb,0,true);
2377
2378 bool (* readrgb_ptr)(const std::string & s,int W,int H,gen & res)=0;
2379
_readrgb(const gen & g,GIAC_CONTEXT)2380 gen _readrgb(const gen & g,GIAC_CONTEXT){
2381 if ( g.type==_STRNG && g.subtype==-1) return g;
2382 vecteur v(gen2vecteur(g));
2383 if (v.empty() || v[0].type!=_STRNG)
2384 return gensizeerr();
2385 int s=int(v.size());
2386 gen res;
2387 bool ok= false;
2388 if (readrgb_ptr)
2389 ok = readrgb_ptr(*v[0]._STRNGptr,
2390 ( (s>2 && v[1].type==_INT_) ?v[1].val:0),
2391 ( (s>2 && v[2].type==_INT_) ?v[2].val:0),
2392 res);
2393 if (!ok)
2394 return gensizeerr(gettext("File not found or unrecognized image file format"));
2395 return res;
2396 }
2397 static const char _readrgb_s []="readrgb";
2398 static define_unary_function_eval (__readrgb,&_readrgb,_readrgb_s);
2399 define_unary_function_ptr5( at_readrgb ,alias_at_readrgb,&__readrgb,0,true);
2400
linear_apply(const gen & e,const gen & x,const gen & l,gen & remains,GIAC_CONTEXT,gen (* f)(const gen &,const gen &,const gen &,gen &,const context *))2401 gen linear_apply(const gen & e,const gen & x,const gen & l,gen & remains, GIAC_CONTEXT, gen (* f)(const gen &,const gen &,const gen &,gen &,const context *)){
2402 if (is_constant_wrt(e,x,contextptr) || (e==x) )
2403 return f(e,x,l,remains,contextptr);
2404 // e must be of type _SYMB
2405 if (e.type!=_SYMB) return gensizeerr(gettext("in linear_apply"));
2406 unary_function_ptr u(e._SYMBptr->sommet);
2407 gen arg(e._SYMBptr->feuille);
2408 gen res;
2409 if (u==at_neg){
2410 res=-linear_apply(arg,x,l,remains,contextptr,f);
2411 remains=-remains;
2412 return res;
2413 } // end at_neg
2414 if (u==at_plus){
2415 if (arg.type!=_VECT)
2416 return linear_apply(arg,x,l,remains,contextptr,f);
2417 const_iterateur it=arg._VECTptr->begin(),itend=arg._VECTptr->end();
2418 for (gen tmp;it!=itend;++it){
2419 res = res + linear_apply(*it,x,l,tmp,contextptr,f);
2420 remains =remains + tmp;
2421 }
2422 return res;
2423 } // end at_plus
2424 if (u==at_prod){
2425 if (arg.type!=_VECT)
2426 return linear_apply(arg,x,l,remains,contextptr,f);
2427 // find all constant terms in the product
2428 vecteur non_constant;
2429 gen prod_constant;
2430 decompose_prod(*arg._VECTptr,x,non_constant,prod_constant,true,contextptr);
2431 if (non_constant.empty()) return gensizeerr(gettext("in linear_apply 2")); // otherwise the product would be constant
2432 if (non_constant.size()==1)
2433 res = linear_apply(non_constant.front(),x,l,remains,contextptr,f);
2434 else
2435 res = f(symbolic(at_prod,non_constant),x,l,remains,contextptr);
2436 remains = prod_constant * remains;
2437 return prod_constant * res;
2438 } // end at_prod
2439 return f(e,x,l,remains,contextptr);
2440 }
2441
2442 // discrete product primitive of a polynomial P
product(const polynome & P,const vecteur & v,const gen & n,gen & remains,GIAC_CONTEXT)2443 gen product(const polynome & P,const vecteur & v,const gen & n,gen & remains,GIAC_CONTEXT){
2444 // factor P, for a linear factor a*n+b, multiply res by a^n*gamma(n+b/a)
2445 // for other factors multiply remains by the factor
2446 polynome Pcont;
2447 factorization f;
2448 gen divan=1,res,extra_div=1;
2449 if (!factor(P,Pcont,f,/* is_primitive*/false,/* with_sqrt*/true,/* complex */true,divan,extra_div) || extra_div!=1){
2450 remains=r2e(P,v,contextptr);
2451 return 1;
2452 }
2453 res = pow(divan,-n,contextptr);
2454 factorization::const_iterator it=f.begin(),itend=f.end();
2455 for (;it!=itend;++it){
2456 gen tmp=r2e(it->fact,v,contextptr);
2457 if (it->fact.lexsorted_degree()!=1){
2458 remains = remains * pow(tmp,it->mult);
2459 }
2460 else {
2461 gen a=derive(tmp,n,contextptr);
2462 if (is_undef(a))
2463 return a;
2464 gen b=normal(tmp-a*n,contextptr);
2465 res = res * pow(a,it->mult*n,contextptr) * pow(symbolic(at_factorial,n+b/a-1),it->mult,contextptr);
2466 }
2467 }
2468 return res*pow(r2e(Pcont,v,contextptr),n,contextptr);
2469 }
2470
2471 // product(P,n=a..b) where the first variable in v is n
product(const polynome & P,const vecteur & v,const gen & n,const gen & a,const gen & b,GIAC_CONTEXT)2472 gen product(const polynome & P,const vecteur & v,const gen & n,const gen & a,const gen & b,GIAC_CONTEXT){
2473 gen remains(1),res=product(P,v,n,remains,contextptr);
2474 res=subst(res,n,b+1,false,contextptr)/subst(res,n,a,false,contextptr);
2475 if (is_one(remains))
2476 return res;
2477 else
2478 return res*symbolic(at_product,gen(makevecteur(remains,n,a,b),_SEQ__VECT));
2479 }
2480
2481 // solve vnext*sol[n+1]+v*sol[n]+vcst=0
2482 // example: u_{n+1}=u_n/n+1 -> n*u_{n+1}-u_n-n=0
2483 // vnext=[1,0], v=[1], vcst=[-1,0] -> n solution
rsolve_particular(const modpoly & vnext,const modpoly & v,const modpoly & vcst,modpoly & sol,GIAC_CONTEXT)2484 static bool rsolve_particular(const modpoly & vnext,const modpoly & v,const modpoly & vcst,modpoly & sol,GIAC_CONTEXT){
2485 // find majoration of degree for solution
2486 int vnextdeg=int(vnext.size())-1,vdeg=int(v.size())-1,vcstdeg=int(vcst.size())-1,soldeg;
2487 soldeg=vcstdeg-giacmax(vnextdeg,vdeg);
2488 if (vnextdeg==vdeg && is_zero(vnext.front()+v.front()))
2489 soldeg++;
2490 if (soldeg<0)
2491 return false;
2492 modpoly vars(soldeg+1);
2493 for (int i=0;i<=soldeg;++i){
2494 vars[i]=identificateur("rsolve_x"+print_INT_(i));
2495 }
2496 modpoly equations=operator_plus(operator_times(vnext,taylor(vars,1,0),0),
2497 operator_times(v,vars,0),0);
2498 equations=operator_plus(equations,vcst,0);
2499 sol=linsolve(equations,vars,contextptr);
2500 return !sol.empty() && !is_undef(sol);
2501 }
2502
2503 // solve a recurrence relation x_{n+1}=l0*x_n+e
rsolve(const gen & e,const gen & n,const gen & l0,gen & remains,GIAC_CONTEXT)2504 static gen rsolve(const gen & e,const gen & n,const gen & l0,gen & remains,GIAC_CONTEXT){
2505 if (is_zero(e)){
2506 remains=0;
2507 return e;
2508 }
2509 vecteur v;
2510 polynome P,Q,R;
2511 if (!is_hypergeometric(e,*n._IDNTptr,v,P,Q,R,contextptr)){
2512 *logptr(contextptr) << gettext("Cst part must be hypergeometric") << '\n';
2513 remains=e;
2514 return 0;
2515 }
2516 if (Q.lexsorted_degree() || R.lexsorted_degree()){
2517 *logptr(contextptr) << gettext("Cst part must be of type a^n*P(n)") << '\n';
2518 remains=e;
2519 return 0;
2520 }
2521 gen a=r2e(Q,v,contextptr)/r2e(R,v,contextptr);
2522 gen P0=r2e(P,v,contextptr);
2523 if (is_zero(P0))
2524 return 0;
2525 for (int i=0;;i++){
2526 gen tmp=normal(subst(e/P0,n,i,false,contextptr),contextptr);
2527 if (!is_undef(tmp)){
2528 P0=tmp*P0/pow(a,i);
2529 break;
2530 }
2531 }
2532 v.clear();
2533 v.push_back(n);
2534 lvar(P0,v);
2535 lvar(l0,v);
2536 lvar(a,v);
2537 vecteur v1(v.begin()+1,v.end());
2538 gen l=e2r(l0,v1,contextptr);
2539 P0=e2r(P0,v,contextptr);
2540 gen P0n,P0d;
2541 fxnd(P0,P0n,P0d);
2542 if (P0n.type!=_POLY)
2543 P=polynome(P0n,int(v.size()));
2544 else
2545 P=*P0n._POLYptr;
2546 a=e2r(a,v1,contextptr);
2547 remains=0;
2548 // solve u(n+1)=l*u(n)+a^n*P(n)
2549 if (a==l){ // let v(n)=u(n)/l^n, then v(n+1)=v(n)+P(n)/l
2550 return pow(l0,n,contextptr)*sum(r2e(P,v,contextptr),n,remains,contextptr)/r2e(l,v1,contextptr)/r2e(P0d,v,contextptr);
2551 }
2552 // search u(n)=a^n*Q(n) with Q a polynomial of same degree than P
2553 // we have u(n+1)=a^(n+1)*Q(n+1)=l*a^n*Q(n)+a^n*P(n)
2554 // hence a*Q(n+1)-l*Q(n)=P(n)
2555 vecteur p=polynome2poly1(P,1);
2556 int pdeg=int(p.size())-1;
2557 vecteur q(pdeg+1);
2558 reverse(p.begin(),p.end());
2559 for (int j=pdeg;j>=0;j--){
2560 gen tmp;
2561 for (int k=j+1;k<=pdeg;++k){
2562 tmp = tmp + comb(k,j)*q[k];
2563 }
2564 q[j]=(p[j]-a*tmp)/(a-l);
2565 }
2566 reverse(q.begin(),q.end());
2567 gen res=r2e(q,v1,contextptr);
2568 if (res.type!=_VECT)
2569 return gentypeerr();
2570 res=symb_horner(*res._VECTptr,n);
2571 res=res*pow(r2e(a,v1,contextptr),n,contextptr)/r2e(P0d,v,contextptr);
2572 return res;
2573 }
2574
2575 /*
2576 void gen2fracpoly1(const gen & g,vecteur & num,vecteur & den){
2577 if (g.type==_FRAC){
2578 num=gen2vecteur(g._FRACptr->num);
2579 den=gen2vecteur(g._FRACptr->den);
2580 }
2581 else {
2582 num=gen2vecteur(g);
2583 den=vecteur(1,1);
2584 }
2585 }
2586 */
2587
2588 // solveseq(expression,[var,value])
2589 // var is a vector of dim the number of terms in the recurrence
_seqsolve(const gen & args,GIAC_CONTEXT)2590 gen _seqsolve(const gen & args,GIAC_CONTEXT){
2591 if ( args.type==_STRNG && args.subtype==-1) return args;
2592 gen f,x=vx_var,uzero=zero,n;
2593 int dim=1,udim=1;
2594 if (args.type==_VECT){
2595 vecteur & v=*args._VECTptr;
2596 int s=int(v.size());
2597 if (!s)
2598 return gentoofewargs("");
2599 f=v[0];
2600 if (s>1)
2601 x=v[1];
2602 if (s>2)
2603 uzero=eval(v[2],eval_level(contextptr),contextptr);
2604 if (x.type==_VECT){
2605 dim=int(x._VECTptr->size());
2606 if (uzero.type==_VECT)
2607 udim=int(uzero._VECTptr->size());
2608 }
2609 }
2610 else
2611 f=args;
2612 x=gen2vecteur(x);
2613 vecteur fv=gen2vecteur(f);
2614 fv=quote_eval(fv,*x._VECTptr,contextptr);
2615 bool fv1=false;
2616 if (fv.size()==1 && fv.front().type==_VECT){
2617 fv1=true;
2618 vecteur tmp=*fv.front()._VECTptr;
2619 fv=tmp;
2620 }
2621 if (dim==udim+1){
2622 n=x._VECTptr->back();
2623 if (n.type!=_IDNT){
2624 identificateur nn(" rsolve_n");
2625 x._VECTptr->back()=nn;
2626 f=quotesubst(f,n,nn,contextptr);
2627 gen res=_seqsolve(gen(makevecteur(f,x,uzero),_SEQ__VECT),contextptr);
2628 return quotesubst(res,nn,n,contextptr);
2629 }
2630 x=vecteur(x._VECTptr->begin(),x._VECTptr->end()-1);
2631 --dim;
2632 }
2633 else {
2634 if (dim!=udim)
2635 return gendimerr();
2636 n=gen(identificateur("rsolve_n"));
2637 }
2638 int fs=int(fv.size());
2639 if (dim<fs) return gensizeerr(contextptr);
2640 int add=dim-fs;
2641 if (f.type!=_VECT && !fv1){
2642 f=vecteur(x._VECTptr->end()-add,x._VECTptr->end());
2643 f._VECTptr->push_back(fv[0]);
2644 }
2645 else {
2646 for (int i=0;i<add;++i)
2647 fv.push_back(x[i]);
2648 f=fv;
2649 }
2650 // check linear recurrence + admissible cst term wrt n
2651 gen difff=derive(f,x,contextptr);
2652 if (is_undef(difff) || difff.type!=_VECT)
2653 return difff;
2654 matrice m=mtran(*difff._VECTptr);
2655 if (!is_zero(derive(m,x,contextptr))){
2656 if (f._VECTptr->size()==1 && x._VECTptr->size()==1 && is_zero(derive(f,n,contextptr))){
2657 // homographic?
2658 gen tmp=ratnormal(f._VECTptr->front(),contextptr),var(x._VECTptr->front());
2659 gen tmpn=_getNum(tmp,contextptr),tmpd=_getDenom(tmp,contextptr),a,b,c,d;
2660 if (is_linear_wrt(tmpn,var,a,b,contextptr) && is_linear_wrt(tmpd,var,c,d,contextptr)&& !is_zero(c)){
2661 // u_{n+1}=(a*u_n+b)/(c*u_n+d)
2662 // find fixed points
2663 gen fixed=solve(a*n+b-(c*n+d)*n,n,1,contextptr);
2664 if (fixed.type==_VECT && fixed._VECTptr->size()==2){
2665 gen r1=fixed._VECTptr->front(),r2=fixed._VECTptr->back();
2666 // (u_n-r1)/(u_n-r2) is geometric, compute ratio
2667 gen un1=(a*n+b)/(c*n+d);
2668 gen r=normal((un1-r1)*(n-r2)/((un1-r2)*(n-r1)),contextptr);
2669 if (is_zero(derive(r,n,contextptr))){
2670 un1=pow(r,n,contextptr)*(var-r1)/(var-r2);
2671 return (r1-un1*r2)/(1-un1);
2672 }
2673 }
2674 if (fixed.type==_VECT && fixed._VECTptr->size()==1){
2675 gen r1=fixed._VECTptr->front();
2676 // 1/(u_n-r1) is arithmetic
2677 gen un1=(a*n+b)/(c*n+d);
2678 gen r=normal(inv(un1-r1,contextptr)-inv(n-r1,contextptr),contextptr);
2679 if (is_zero(derive(r,n,contextptr))){
2680 un1=r*n+inv(var-r1,contextptr);
2681 return inv(un1,contextptr)+r1;
2682 }
2683 }
2684 }
2685 }
2686 return symbolic(at_seqsolve,args);
2687 }
2688 vecteur cst=*normal(subvecteur(*f._VECTptr,multmatvecteur(m,*x._VECTptr)),contextptr)._VECTptr;
2689 if (m.size()==1 && cst.size()==1){
2690 gen l(m[0][0]),c(cst[0]),remains;
2691 gen u0=gen2vecteur(uzero)[0];
2692 if (n.type!=_IDNT)
2693 return gentypeerr();
2694 if (!is_zero(derive(l,n,contextptr))){
2695 vecteur v;
2696 polynome P,Q,R;
2697 if (!is_hypergeometric(l,*n._IDNTptr,v,P,Q,R,contextptr)){
2698 *logptr(contextptr) << gettext("Cst part must be hypergeometric") << '\n';
2699 return symbolic(at_seqsolve,args);
2700 }
2701 // l(n+1)/l(n) = P(n+1)/P(n)*Q(n)/R(n+1)
2702 std::vector< monomial<gen> >::const_iterator it=Q.coord.begin();
2703 polynome q0=Tnextcoeff<gen>(it,Q.coord.end()).untrunc1();
2704 it=R.coord.begin();
2705 polynome r0=Tnextcoeff<gen>(it,R.coord.end()).untrunc1();
2706 if (!is_zero(r0*Q-q0*R)){
2707 *logptr(contextptr) << gettext("Unable to handle coeff of homogeneous part") << '\n';
2708 return symbolic(at_seqsolve,args);
2709 }
2710 // -> l(n)=l(0)*P(n)/P(0) -> product(l(n))
2711 gen pn=r2e(P,v,contextptr)/r2e(Q,v,contextptr);
2712 gen q0r0=r2e(q0,v,contextptr)/r2e(r0,v,contextptr);
2713 gen res=pow(subst(normal(l/pn,contextptr),n,0,false,contextptr),n,contextptr)*simplify(product(P,v,n,0,n-1,contextptr)/product(Q,v,n,0,n-1,contextptr),contextptr)*pow(q0r0,n*(n-1)/2,contextptr);
2714 // then we might search for a polynomial particular solution to e
2715 if (is_zero(c))
2716 return u0/subst(res,n,0,false,contextptr)*res;
2717 if (!is_one(q0r0))
2718 return gensizeerr("Unable to find particular solution, general solution is "+res.print(contextptr));
2719 // u_{n+1}=l*u_{n}+c
2720 gen tmp=l*x[0]+c,tmpnum,tmpden;
2721 vecteur tmpv(1,n);
2722 lvar(tmp,tmpv);
2723 tmp=e2r(tmp,tmpv,contextptr);
2724 fxnd(tmp,tmpnum,tmpden);
2725 tmpnum=r2e(tmpnum,tmpv,contextptr);
2726 tmpden=r2e(tmpden,tmpv,contextptr);
2727 tmpnum=_e2r(gen(makevecteur(tmpnum,n),_SEQ__VECT),contextptr);
2728 tmpden=_e2r(gen(makevecteur(tmpden,n),_SEQ__VECT),contextptr);
2729 if (is_zero(derive(tmpnum,n,contextptr)) &&
2730 is_zero(derive(tmpnum,n,contextptr)) &&
2731 tmpnum.type==_VECT && tmpden.type==_VECT){
2732 vecteur tmpn,tmpd,ln,cn;
2733 tmpn=*tmpnum._VECTptr;
2734 tmpd=*tmpden._VECTptr;
2735 gen difftmpn=derive(tmpn,x[0],contextptr);
2736 if (is_undef(difftmpn) || difftmpn.type!=_VECT)
2737 return gensizeerr(contextptr);
2738 ln=trim(*difftmpn._VECTptr,0);
2739 cn=trim(subst(tmpn,x[0],0,false,contextptr),0);
2740 if (rsolve_particular(-tmpd,ln,cn,tmpn,contextptr)){
2741 gen sol=_r2e(gen(makevecteur(tmpn,n),_SEQ__VECT),contextptr);
2742 // general solution is sol+C*res, sol(0)+C*res(0)=u0
2743 gen C=subst((u0-sol)/res,n,0,false,contextptr);
2744 return sol+C*res;
2745 }
2746 }
2747 *logptr(contextptr) << gettext("Unable to find a particular solution for inhomogeneous part") << '\n';
2748 return symbolic(at_seqsolve,args);
2749 }
2750 gen d=linear_apply(c,n,l,remains,contextptr,rsolve);
2751 if (!is_zero(remains))
2752 *logptr(contextptr) << gettext("Unable to solve recurrence") << '\n';
2753 // d is a particular solution of u(n+1)=l*u(n)+c(n)
2754 // add a general solution d(n)+C*l^n
2755 // such that at n=0 we get u0 -> C+d(0)=u0
2756 gen C=normal(u0-quotesubst(d,n,0,contextptr),contextptr);
2757 return d+C*pow(l,n,contextptr);
2758 }
2759 // u(n+1)=M*u(n)+cst
2760 // Let M=PDP^-1, v=P^-1*u, v satisfies v(n+1)=Dv(n)+P^-1*cst
2761 // solve for v then for u
2762 if (!is_zero(derive(m,n,contextptr)))
2763 return symbolic(at_seqsolve,args);
2764 if (has_num_coeff(m))
2765 m=*evalf(m,1,contextptr)._VECTptr;
2766 matrice P,Pinv,D;
2767 bool b=complex_mode(contextptr);
2768 complex_mode(true,contextptr);
2769 egv(m,P,D,contextptr,true,false,false);
2770 complex_mode(b,contextptr);
2771 Pinv=minv(P,contextptr);
2772 if (is_undef(Pinv))
2773 return Pinv;
2774 // if (fs==1) uzero=_revlist(uzero);
2775 vecteur vzero=multmatvecteur(Pinv,*uzero._VECTptr);
2776 vecteur vcst=multmatvecteur(Pinv,cst);
2777 int taille=int(m.size());
2778 vecteur res(taille);
2779 for (int i=taille-1;i>=0;--i){
2780 // find cst coefficient
2781 gen c=vcst[i],l=D[i][i];
2782 for (int j=i+1;j<taille;j++){
2783 c += D[i][j]*res[j];
2784 }
2785 c=normal(c,contextptr);
2786 gen remains,d=linear_apply(c,n,l,remains,contextptr,rsolve);
2787 if (!is_zero(remains))
2788 *logptr(contextptr) << gettext("Unable to solve recurrence") << '\n';
2789 gen C=normal(vzero[i]-quotesubst(d,n,0,contextptr),contextptr);
2790 if (is_zero(l))
2791 res[i]=d+C*symbolic(at_same,gen(makevecteur(n,0),_SEQ__VECT));
2792 else
2793 res[i]=d+C*pow(l,n,contextptr);
2794 }
2795 // u=P*v
2796 res=multmatvecteur(P,res);
2797 if (fs==1)
2798 //return normal(subst(res[0],n,n-add,false,contextptr),contextptr);// normal(res[0]);
2799 return ratnormal(res[0],contextptr);
2800 else
2801 return ratnormal(res,contextptr);
2802 }
2803 static const char _seqsolve_s []="seqsolve";
2804 static define_unary_function_eval_quoted (__seqsolve,&_seqsolve,_seqsolve_s);
2805 define_unary_function_ptr5( at_seqsolve ,alias_at_seqsolve,&__seqsolve,_QUOTE_ARGUMENTS,true);
2806
rsolve_initcond(const vecteur & initcond,const gen & un,const gen & u,const gen & n,const vecteur & uinit,GIAC_CONTEXT)2807 static vecteur rsolve_initcond(const vecteur & initcond,const gen & un,const gen & u,const gen & n,const vecteur & uinit,GIAC_CONTEXT){
2808 if (initcond.empty())
2809 return makevecteur(un);
2810 gen initc=subst(initcond,u,_unapply(gen(makevecteur(un,n),_SEQ__VECT),contextptr),false,contextptr);
2811 initc=initc.eval(1,contextptr);
2812 if (initc.type!=_VECT)
2813 return vecteur(1,gensizeerr());
2814 vecteur valv=gsolve(*initc._VECTptr,uinit,/* complex mode */ true,0,contextptr);
2815 if (is_undef(valv))
2816 return valv;
2817 vecteur resv(valv.size());
2818 for (unsigned int i=0;i<valv.size();++i){
2819 resv[i]=normal(quotesubst(un,uinit,valv[i],contextptr),contextptr);
2820 }
2821 return resv;
2822 }
2823
2824 // example f : u(n+1)=2*u(n)+n, initcond [u(0)=1]
crsolve(const gen & f0,const gen & u,const gen & n,vecteur & initcond0,GIAC_CONTEXT)2825 static gen crsolve(const gen & f0,const gen &u,const gen & n,vecteur & initcond0,GIAC_CONTEXT){
2826 if (f0.type==_VECT){
2827 if (u.type!=_VECT || f0._VECTptr->size()!=u._VECTptr->size())
2828 return gendimerr();
2829 }
2830 else {
2831 if (u.type==_VECT)
2832 return gendimerr();
2833 }
2834 vecteur uv(gen2vecteur(u));
2835 int uvs=int(uv.size());
2836 vecteur initcond;
2837 aplatir(*apply(initcond0,equal2diff)._VECTptr,initcond);
2838 gen f=apply(f0,equal2diff);
2839 if (n.type!=_IDNT){
2840 identificateur N(" rsolve_N");
2841 gen F=quotesubst(f,n,N,contextptr);
2842 gen tmp=crsolve(F,u,N,initcond0,contextptr);
2843 return quotesubst(tmp,N,n,contextptr);
2844 }
2845 vecteur vof0(lop(f,at_of));
2846 // keep only those with u
2847 vecteur vofa,vofb,vopos,vofun;
2848 bool all_one=true,all_bzero=true;
2849 for (const_iterateur it=vof0.begin();it!=vof0.end();++it){
2850 gen & itf=it->_SYMBptr->feuille;
2851 int pos;
2852 if (itf.type==_VECT && itf._VECTptr->size()==2 && (pos=equalposcomp(uv,itf._VECTptr->front())) ){
2853 gen tmp=it->_SYMBptr->feuille._VECTptr->back(),a,b;
2854 if (is_linear_wrt(tmp,n,a,b,contextptr) && !is_zero(a)){
2855 if (!is_one(a))
2856 all_one=false;
2857 if (!is_zero(b))
2858 all_bzero=false;
2859 vofa.push_back(a);
2860 vofb.push_back(b);
2861 vopos.push_back(pos-1);
2862 vofun.push_back(*it);
2863 }
2864 else
2865 return gensizeerr(gettext("Unable to handle this kind of recurrence"));
2866 }
2867 }
2868 if (vofa.empty())
2869 return undef;
2870 if (all_one){
2871 gen bmin=vofb.front(),bmax=vofb.front();
2872 int nmax=0;
2873 for (const_iterateur it=vofb.begin();it!=vofb.end();++it){
2874 const gen & tmp = *it;
2875 if (is_strictly_greater(tmp,bmax,contextptr)){
2876 nmax=int(it-vofb.begin());
2877 bmax=tmp;
2878 }
2879 if (is_strictly_greater(bmin,tmp,contextptr))
2880 bmin=tmp;
2881 }
2882 if (bmin.type!=_INT_ || bmax.type!=_INT_)
2883 return gensizeerr(gettext("Bad indexes in recurrence"));
2884 int m=bmin.val,M=bmax.val;
2885 if (uvs>1){
2886 if (M!=m+1)
2887 return gendimerr(gettext("Multi-recurrences implemented only for 1 step"));
2888 // generate identifiers x0,...,x_{uvs-1} and y0,...,y_{uvs-1}
2889 vecteur idx,idX;
2890 for (int i=0;i<uvs;++i){
2891 idx.push_back(gen(identificateur("rsolve_x"+print_INT_(i))));
2892 idX.push_back(gen(identificateur("rsolve_X"+print_INT_(i))));
2893 }
2894 vecteur vofid;
2895 for (const_iterateur it=vofb.begin();it!=vofb.end();++it){
2896 if (it->type!=_INT_)
2897 return gensizeerr();
2898 bool isx=it->val==m;
2899 int pos=vopos[it-vofb.begin()].val;
2900 vofid.push_back(isx?idx[pos]:idX[pos]);
2901 }
2902 gen F=quotesubst(f,vofun,vofid,contextptr);
2903 if (F.type!=_VECT)
2904 return gensizeerr();
2905 // solves for idX in terms of idx
2906 vecteur Fv=gsolve(*F._VECTptr,idX,1/* complex */,0/* approx mode=no*/,contextptr),res;
2907 if (is_undef(Fv))
2908 return Fv;
2909 for (const_iterateur it=Fv.begin();it!=Fv.end();++it){
2910 vecteur idxn(idx);
2911 idxn.push_back(n);
2912 gen vn=_seqsolve(makevecteur(*it,idxn,idx),contextptr);
2913 gen un=quotesubst(vn,n,n+1-M,contextptr);
2914 if (un.type!=_VECT)
2915 return gensizeerr("Unable to solve this recurrence");
2916 if (int(un._VECTptr->size())!=uvs)
2917 return gendimerr(contextptr);
2918 if (initcond.empty())
2919 return vecteur(1,un);
2920 // solve initial conditions
2921 vecteur vout;
2922 for (int i=0;i<uvs;++i){
2923 vout.push_back(_unapply(gen(makevecteur(un[i],n),_SEQ__VECT),contextptr));
2924 }
2925 gen initc=subst(initcond,uv,vout,false,contextptr);
2926 initc=initc.eval(1,contextptr);
2927 if (initc.type!=_VECT)
2928 return gensizeerr();
2929 vecteur valv=gsolve(*initc._VECTptr,idx,/* complex mode */ true,/* approx=no */ 0,contextptr);
2930 if (is_undef(valv))
2931 return valv;
2932 for (unsigned int i=0;i<valv.size();++i){
2933 res.push_back(normal(quotesubst(un,idx,valv[i],contextptr),contextptr));
2934 }
2935 }
2936 return res;
2937 }
2938 // generate identifiers xm,...,xM
2939 vecteur ids;
2940 for (int i=m;i<=M;++i){
2941 ids.push_back(gen(identificateur("rsolve_x"+print_INT_(i-m))));
2942 }
2943 // replace vofun
2944 vecteur vofid;
2945 for (const_iterateur it=vofb.begin();it!=vofb.end();++it){
2946 if (it->type!=_INT_)
2947 return gensizeerr();
2948 vofid.push_back(ids[it->val-m]);
2949 }
2950 gen F=quotesubst(f,vofun,vofid,contextptr),a,b;
2951 if (!is_linear_wrt(F,ids.back(),a,b,contextptr))
2952 return gensizeerr();
2953 F=-b/a; // xM in terms of xm to xM-1
2954 ids.back()=n;
2955 vecteur uinit(ids);
2956 uinit.pop_back();
2957 // let v_n=u_{n+m} -> u_n=v_{n-m}
2958 // old code
2959 // gen vn=_seqsolve(makevecteur(F,ids,uinit),contextptr);
2960 // gen un=quotesubst(vn,n,n-m,contextptr);
2961 // new code that solves rsolve(u(n)=(n)*u(n-1),u(n),u(1)=3);
2962 gen Fm=quotesubst(F,n,n-m,contextptr);
2963 gen un=_seqsolve(makevecteur(Fm,ids,uinit),contextptr);
2964 // end new code
2965 return rsolve_initcond(initcond,un,u,n,uinit,contextptr);
2966 }
2967 // divide and conquier?
2968 if (uvs==1 && vofa.size()==2 && all_bzero){
2969 gen idy=identificateur("rsolve_y"),x(identificateur("rsolve_x"));
2970 vecteur ids(makevecteur(x,idy));
2971 // replace vofun
2972 vecteur vofid;
2973 gen b;
2974 if (is_one(vofa[0])){
2975 vofid=ids;
2976 b=vofa[1];
2977 }
2978 if (is_one(vofa[1])){
2979 vofid=makevecteur(ids[1],ids[0]);
2980 b=vofa[0];
2981 }
2982 if (vofid.empty() || ck_is_greater(1,b,contextptr))
2983 return gensizeerr();
2984 vofun.push_back(n);
2985 vofid.push_back(pow(b,n,contextptr));
2986 gen F=quotesubst(f,vofun,vofid,contextptr);
2987 // F(x,y,n)=0 where y=t(b*n) x=t(n)
2988 // auxiliary sequence u(n)=t(b^n), verifies F(u(n+1),u(n),b^n)=0
2989 gen A,B;
2990 if (is_linear_wrt(F,idy,A,B,contextptr)){
2991 F=-B/A; // u(n+1) in terms of u(n)=x
2992 gen vn=_seqsolve(makevecteur(F,makevecteur(x,n),x),contextptr);
2993 gen un=quotesubst(vn,n,ln(n,contextptr)/ln(b,contextptr),contextptr);
2994 // solve initial cond
2995 return rsolve_initcond(initcond,un,u,n,makevecteur(x),contextptr);
2996 }
2997 }
2998 return gensizeerr(gettext("Not yet implemented"));
2999 }
rsolve(const gen & f0,const gen & u,const gen & n,vecteur & initcond0,int st,GIAC_CONTEXT)3000 static gen rsolve(const gen & f0,const gen &u,const gen & n,vecteur & initcond0,int st,GIAC_CONTEXT){
3001 bool b=complex_mode(contextptr);
3002 complex_mode(true,contextptr);
3003 gen res=crsolve(f0,u,n,initcond0,contextptr);
3004 if (!b){
3005 complex_mode(b,contextptr);
3006 // FIXME take real part for res
3007 }
3008 return res;
3009 }
3010 // example args=u(n+1)=2*u(n)+n,u(n),u(0)=1
_rsolve(const gen & args,GIAC_CONTEXT)3011 gen _rsolve(const gen & args,GIAC_CONTEXT){
3012 if ( args.type==_STRNG && args.subtype==-1) return args;
3013 vecteur varg=gen2vecteur(args);
3014 if (debug_infolevel>20)
3015 varg.dbgprint();
3016 int s=int(varg.size());
3017 if (!s)
3018 return gendimerr();
3019 gen f,u,n;
3020 vecteur initcond;
3021 if (s>1){
3022 gen un=varg[1];
3023 if (un.is_symb_of_sommet(at_of)){
3024 gen & unf=un._SYMBptr->feuille;
3025 if (unf.type==_VECT && unf._VECTptr->size()==2){
3026 u=unf._VECTptr->front();
3027 n=unf._VECTptr->back();
3028 }
3029 }
3030 if (un.type==_VECT){
3031 vecteur & unv=*un._VECTptr;
3032 vecteur uv;
3033 for (const_iterateur it=unv.begin();it!=unv.end();++it){
3034 const gen & itg=*it;
3035 if (itg.is_symb_of_sommet(at_of)){
3036 gen & unf=itg._SYMBptr->feuille;
3037 if (unf.type==_VECT && unf._VECTptr->size()==2){
3038 uv.push_back(unf._VECTptr->front());
3039 if (is_zero(n))
3040 n=unf._VECTptr->back();
3041 else
3042 if (n!=unf._VECTptr->back())
3043 return gentypeerr();
3044 }
3045 }
3046 }
3047 u=uv;
3048 }
3049 if (is_zero(n))
3050 return gentypeerr();
3051 }
3052 else {
3053 u=gen(identificateur("rsolve_u"));
3054 n=gen(identificateur("rsolve_n"));
3055 }
3056 vecteur quoted=gen2vecteur(u);
3057 quoted.push_back(n);
3058 varg=quote_eval(varg,quoted,contextptr);
3059 f=varg[0];
3060 if (s>2){
3061 initcond=vecteur(varg.begin()+2,varg.end());
3062 if (initcond.size()==1 && initcond.front().type==_VECT){
3063 // use a temporary vector,
3064 // otherwise the source is destroyed when copying the first element of the vector
3065 // another way to do that could be to have a gen with a copy of initcond.front
3066 vecteur tmp=*initcond.front()._VECTptr;
3067 initcond=tmp;
3068 }
3069 }
3070 else {
3071 if (f.type==_VECT && !f._VECTptr->empty()){
3072 if (u.type!=_VECT || f._VECTptr->size()!=u._VECTptr->size()){
3073 initcond=*f._VECTptr;
3074 f=initcond.front();
3075 initcond.erase(initcond.begin());
3076 }
3077 }
3078 }
3079 int st=step_infolevel(contextptr);
3080 step_infolevel(0,contextptr);
3081 gen res=rsolve(f,u,n,initcond,st,contextptr);
3082 step_infolevel(st,contextptr);
3083 return res;
3084 }
3085 static const char _rsolve_s []="rsolve";
3086 static define_unary_function_eval_quoted (__rsolve,&_rsolve,_rsolve_s);
3087 define_unary_function_ptr5( at_rsolve ,alias_at_rsolve,&__rsolve,_QUOTE_ARGUMENTS,true);
3088
islessthanf(const gen & a,const gen & b)3089 static bool islessthanf(const gen & a,const gen & b){
3090 if (a.type!=_VECT || b.type!=_VECT)
3091 return is_strictly_greater(b,a,context0);
3092 vecteur & av =*a._VECTptr;
3093 vecteur & bv =*b._VECTptr;
3094 int avs=int(av.size()),bvs=int(bv.size());
3095 if (avs!=bvs)
3096 return avs<bvs;
3097 for (int i=0;i<avs;++i){
3098 if (av[i]==bv[i])
3099 continue;
3100 return islessthanf(av[i],bv[i]);
3101 }
3102 return false;
3103 }
3104
_array(const gen & g,GIAC_CONTEXT)3105 gen _array(const gen & g,GIAC_CONTEXT){
3106 if ( g.type==_STRNG && g.subtype==-1) return g;
3107 // create a table with specified index, may be initialized
3108 vecteur args;
3109 if (g.type==_VECT && g.subtype==_SEQ__VECT){
3110 args=*g._VECTptr;
3111 if (!args.empty()){
3112 gen f=args.back();
3113 if (f==at_float || (f.is_symb_of_sommet(at_equal) && f._SYMBptr->feuille[0]==at_dtype)){
3114 args.pop_back();
3115 return _convert(makesequence(_array(args.size()==1?args.front():gen(args,_SEQ__VECT),contextptr),f.type==_FUNC?f:f._SYMBptr->feuille[1]),contextptr);
3116 }
3117 }
3118 }
3119 else
3120 args=gen2vecteur(g);
3121 #if 1 // def NSPIRE
3122 gen_map m;
3123 #else
3124 gen_map m(ptr_fun(islessthanf));
3125 #endif
3126 int s=int(args.size());
3127 vector<int> indexbegin,indexsize;
3128 int nindexes=1;
3129 gen initv(vecteur(0));
3130 vecteur & init = *initv._VECTptr;
3131 int shift = array_start(contextptr); //xcas_mode(contextptr)!=0 || abs_calc_mode(contextptr)==38;
3132 for (int i=0;i<s;++i){
3133 if (args[i].is_symb_of_sommet(at_interval)){
3134 gen & f =args[i]._SYMBptr->feuille;
3135 if (f.type!=_VECT || f._VECTptr->size()!=2)
3136 return gendimerr(args[i].print(contextptr));
3137 vecteur & fv= *f._VECTptr;
3138 if (fv[0].type!=_INT_ || fv[1].type!=_INT_ || fv[1].val<fv[0].val)
3139 return gendimerr(args[i].print(contextptr));
3140 indexbegin.push_back(fv.front().val-shift);
3141 indexsize.push_back(fv.back().val-fv.front().val+1);
3142 nindexes *= indexsize.back();
3143 }
3144 if (args[i].type==_VECT)
3145 init=mergevecteur(init,*args[i]._VECTptr);
3146 }
3147 if (nindexes>>24)
3148 return gendimerr(gettext("Array too large")+print_INT_(nindexes));
3149 if (nindexes==1)
3150 return g;
3151 int is=int(indexsize.size());
3152 for (int i=0;i<nindexes;++i){
3153 // generate index by writing nindexes in bases indexsize
3154 vecteur curidx(is);
3155 int pos=i,posj;
3156 gen initval=initv;
3157 bool initialize=true;
3158 for (int j=0;j<is;++j){
3159 posj = (pos % indexsize[j]);
3160 if (initialize && initval.type==_VECT && (int)initval._VECTptr->size()>posj)
3161 initval=(*initval._VECTptr)[posj];
3162 else
3163 initialize=false;
3164 curidx[j] = posj + indexbegin[j];
3165 pos /= indexsize[j];
3166 }
3167 // initialize m[curidx]
3168 if (initialize)
3169 m[curidx]=initval;
3170 else
3171 m[curidx]=0;
3172 }
3173 gen res=m;
3174 res.subtype=1;
3175 return res;
3176 }
3177 static const char _array_s []="array";
3178 static define_unary_function_eval (__array,&_array,_array_s);
3179 define_unary_function_ptr5( at_array ,alias_at_array,&__array,0,true);
3180
3181 static const char _whattype_s []="whattype";
3182 static define_unary_function_eval (__whattype,&_type,_whattype_s);
3183 define_unary_function_ptr5( at_whattype ,alias_at_whattype,&__whattype,0,true);
3184
3185 static const char _modp_s []="modp";
3186 static define_unary_function_eval (__modp,&_irem,_modp_s);
3187 define_unary_function_ptr5( at_modp ,alias_at_modp,&__modp,0,true);
3188
_makemod(const gen & args,GIAC_CONTEXT)3189 gen _makemod(const gen & args,GIAC_CONTEXT){
3190 if ( args.type==_STRNG && args.subtype==-1) return args;
3191 if (args.type!=_VECT || args._VECTptr->size()!=2)
3192 return gentypeerr();
3193 gen a=args._VECTptr->front(),b=args._VECTptr->back();
3194 if (is_zero(b))
3195 return unmod(a);
3196 if (!is_integer(a) || !is_integer(b))
3197 return gentypeerr();
3198 return makemod(a,b);
3199 }
3200 static const char _makemod_s []="makemod";
3201 static define_unary_function_eval (__makemod,&_makemod,_makemod_s);
3202 define_unary_function_ptr5( at_makemod ,alias_at_makemod,&__makemod,0,true);
3203
cpp_convert_0(const gen & g,GIAC_CONTEXT)3204 gen cpp_convert_0(const gen &g,GIAC_CONTEXT){
3205 return g;
3206 }
3207
cpp_convert_2(const gen & g,GIAC_CONTEXT)3208 longlong cpp_convert_2(const gen & g,GIAC_CONTEXT){
3209 if (g.type==_INT_)
3210 return g.val;
3211 gen h=g;
3212 if (!is_integral(h)){
3213 gensizeerr(contextptr);
3214 return 0;
3215 }
3216 if (h.type==_INT_)
3217 return h.val;
3218 #ifdef USE_GMP_REPLACEMENTS
3219 gensizeerr(contextptr);
3220 return 0;
3221 #else
3222 if (mpz_sizeinbase(*h._ZINTptr,2)>62){
3223 gensizeerr(contextptr);
3224 return 0;
3225 }
3226 if (is_greater(0,g,context0))
3227 return -cpp_convert_2(-g,contextptr);
3228 unsigned int lo, hi;
3229 mpz_t tmp;
3230 mpz_init( tmp );
3231 mpz_mod_2exp( tmp, *h._ZINTptr, 64 ); /* tmp = (lower 64 bits of n) */
3232 lo = mpz_get_ui( tmp ); /* lo = tmp & 0xffffffff */
3233 mpz_div_2exp( tmp, tmp, 32 ); /* tmp >>= 32 */
3234 hi = mpz_get_ui( tmp ); /* hi = tmp & 0xffffffff */
3235 mpz_clear( tmp );
3236 longlong i=(((unsigned long long)hi) << 32) + lo;
3237 return i;
3238 #endif
3239 }
3240
cpp_convert_1(const gen & g,GIAC_CONTEXT)3241 double cpp_convert_1(const gen & g,GIAC_CONTEXT){
3242 if (g.type==_DOUBLE_) return g._DOUBLE_val;
3243 gen h=evalf_double(g,1,context0);
3244 if (h.type!=_DOUBLE_){
3245 gensizeerr(contextptr);
3246 return 0;
3247 }
3248 return h._DOUBLE_val;
3249 }
3250
cpp_convert_4(const gen & g,GIAC_CONTEXT)3251 complex<double> cpp_convert_4(const gen & g,GIAC_CONTEXT){
3252 gen h=evalf_double(g,1,context0);
3253 if (h.type==_DOUBLE_)
3254 return h._DOUBLE_val;
3255 if (h.type!=_CPLX || h.subtype!=3){
3256 gensizeerr(contextptr);
3257 return 0;
3258 }
3259 return complex<double>(h._CPLXptr->_DOUBLE_val,(h._CPLXptr+1)->_DOUBLE_val);
3260 }
3261
cpp_convert_7(const gen & g,GIAC_CONTEXT)3262 vecteur cpp_convert_7(const gen & g,GIAC_CONTEXT){
3263 if (g.type!=_VECT){
3264 gensizeerr(contextptr);
3265 return 0;
3266 }
3267 return *g._VECTptr;
3268 }
3269
cpp_convert_12(const gen & g,GIAC_CONTEXT)3270 std::string cpp_convert_12(const gen & g,GIAC_CONTEXT){
3271 if (g.type!=_STRNG){
3272 gensizeerr(contextptr);
3273 return "";
3274 }
3275 return *g._STRNGptr;
3276 }
3277
is_int_or_double(const gen & g)3278 int is_int_or_double(const gen & g){
3279 if (g.type==_INT_)
3280 return 2;
3281 if (g.type==_DOUBLE_)
3282 return _DOUBLE_;
3283 if (g.type==_IDNT){
3284 const char * ch =g._IDNTptr->id_name;
3285 int s=int(strlen(ch));
3286 if (s>=3 && ch[s-2]=='_') {
3287 switch (ch[s-1]){
3288 case 'i': case 'l':
3289 return 2;
3290 case 'd':
3291 return _DOUBLE_;
3292 }
3293 }
3294 return 0;
3295 }
3296 if (g.type==_SYMB){
3297 if (g._SYMBptr->sommet==at_floor || g._SYMBptr->sommet==at_irem || g._SYMBptr->sommet==at_iquo || g._SYMBptr->sommet==at_size)
3298 return 2;
3299 vecteur v=lvar(g._SYMBptr->feuille);
3300 return is_int_or_double(v); // this assumes that the function has the same type as the arguments
3301 }
3302 if (g.type!=_VECT)
3303 return 0;
3304 const vecteur & v=*g._VECTptr;
3305 for (int i=0;i<int(v.size());++i){
3306 if (!is_int_or_double(v[i]))
3307 return 0;
3308 }
3309 return 3;
3310 }
3311
cpp_vartype(const gen & g,const vecteur & allnames,const vecteur & allchecks)3312 int cpp_vartype(const gen & g,const vecteur & allnames,const vecteur & allchecks){
3313 switch (g.type){
3314 case _INT_: case _ZINT:
3315 return 2;
3316 case _DOUBLE_:
3317 return _DOUBLE_;
3318 case _CPLX:
3319 return _CPLX;
3320 case _VECT:
3321 return _VECT;
3322 case _STRNG:
3323 return _STRNG;
3324 }
3325 if (g.type==_SYMB){
3326 if (g._SYMBptr->sommet==at_at)
3327 return 0; // could be improved if we add support for vector of long/double
3328 if (g._SYMBptr->sommet==at_of){
3329 gen f=g._SYMBptr->feuille[0];
3330 int pos=0;
3331 if ( (pos=equalposcomp(allnames,f)) && pos<=allchecks.size()){
3332 gen g=allchecks[pos-1];
3333 if (g.is_symb_of_sommet(at_deuxpoints) && g._SYMBptr->feuille[0]==f)
3334 return g._SYMBptr->feuille[1].val;
3335 }
3336 }
3337 if (g._SYMBptr->sommet==at_floor || g._SYMBptr->sommet==at_irem || g._SYMBptr->sommet==at_iquo)
3338 return 2;
3339 if (g._SYMBptr->sommet==at_size)
3340 return 2;
3341 if (g._SYMBptr->sommet==at_abs && cpp_vartype(g._SYMBptr->feuille,allnames,allchecks)==3)
3342 return 2;
3343 vecteur v=lvar(g._SYMBptr->feuille);
3344 const_iterateur it=v.begin(),itend=v.end();
3345 int cur=-1;
3346 for (;it!=itend;++it){
3347 int i=cpp_vartype(*it,allnames,allchecks);
3348 if (i==0)
3349 return 0;
3350 if (cur==2 && i==_DOUBLE_)
3351 cur=_DOUBLE_;
3352 if (cur==i || (cur==_DOUBLE_ && i==2))
3353 continue;
3354 if (cur==-1)
3355 cur=i;
3356 else
3357 return 0;
3358 }
3359 return cur;
3360 }
3361 if (g.type!=_IDNT)
3362 return 0;
3363 const char * ch=g._IDNTptr->id_name;
3364 int cl=int(strlen(ch));
3365 if (cl>=3 && ch[cl-2]=='_'){
3366 switch (ch[cl-1]){
3367 case 'i':case 'l':
3368 return 2;
3369 case 'd':
3370 return _DOUBLE_;
3371 case 'c':
3372 return _CPLX;
3373 case 'v':
3374 return _VECT;
3375 case 's':
3376 return _STRNG;
3377 }
3378 }
3379 return 0;
3380 }
3381
cprintvars(const gen & f0,bool cst,bool ctx0,const string & sep,int * typeptr,const vecteur & allnames,const vecteur & allchecks,GIAC_CONTEXT)3382 std::string cprintvars(const gen & f0,bool cst,bool ctx0,const string & sep,int * typeptr,const vecteur & allnames,const vecteur & allchecks,GIAC_CONTEXT){
3383 vecteur fv=gen2vecteur(f0);
3384 if (fv.size()==2 && fv[0].type==_VECT)
3385 fv=*fv[0]._VECTptr;
3386 string res;
3387 for (int i=0;;){
3388 gen fvi=fv[i];
3389 bool eq=fvi.is_symb_of_sommet(at_equal);
3390 bool eqsto=fvi.is_symb_of_sommet(at_sto);
3391 if (eq)
3392 fvi=fvi[1];
3393 if (eqsto)
3394 fvi=fvi[2];
3395 if (fvi.type!=_IDNT)
3396 return "Invalid parameter "+fvi.print(contextptr);
3397 const char * ch=fvi._IDNTptr->id_name;
3398 int cl=int(strlen(ch));
3399 string vtype=cst?"const giac::gen & ":"giac::gen ";
3400 if (typeptr) *typeptr=0;
3401 if (cl>=3 && ch[cl-2]=='_'){
3402 switch (ch[cl-1]){
3403 case 'i': case 'l':
3404 vtype="long ";
3405 if (typeptr) *typeptr=2;
3406 break;
3407 case 'd':
3408 vtype="double ";
3409 if (typeptr) *typeptr=_DOUBLE_;
3410 break;
3411 case 'c':
3412 vtype="complex<double> ";
3413 if (typeptr) *typeptr=_CPLX;
3414 break;
3415 case 'v':
3416 vtype="giac::vecteur ";
3417 if (typeptr) *typeptr=_VECT;
3418 break;
3419 case 's':
3420 vtype="std::string ";
3421 if (typeptr) *typeptr=_STRNG;
3422 break;
3423 }
3424 }
3425 if (typeptr) ++typeptr;
3426 vtype = vtype+ch;
3427 if (eq || eqsto){
3428 gen g=(fv[i])[1];
3429 if (eq)
3430 g=(fv[i])[2];
3431 #if 1
3432 int i=cpp_vartype(fvi,allnames,allchecks);
3433 int gt=cpp_vartype(g,allnames,allchecks);
3434 if (!i || i==gt || (i==2 && gt==_DOUBLE_) || (i==_DOUBLE_ && gt==2))
3435 vtype = vtype + " = " +cprint(g,0,contextptr);
3436 else
3437 vtype = vtype + " = cpp_convert_"+print_INT_(i)+"("+cprint(g,0,contextptr)+(ctx0?",giac::context0)":",contextptr)");
3438 #else
3439 vtype = vtype+'=';
3440 int i=is_int_or_double(fvi);
3441 if (!i || is_int_or_double(g))
3442 vtype = vtype +g.print(contextptr);
3443 else {
3444 vtype = vtype + "("+g.print(contextptr);
3445 vtype += (i==_DOUBLE_)?")._DOUBLE_val":").val";
3446 }
3447 #endif
3448 }
3449 ++i;
3450 if (i==int(fv.size()))
3451 return res+vtype;
3452 res = res + vtype + sep;
3453 }
3454 }
3455
cpp_stoprint(const gen & g,const vecteur & allnames,const vecteur & allchecks,GIAC_CONTEXT)3456 std::string cpp_stoprint(const gen & g,const vecteur & allnames,const vecteur & allchecks,GIAC_CONTEXT){
3457 if (g.is_symb_of_sommet(at_at) && g._SYMBptr->feuille.type==_VECT && g._SYMBptr->feuille._VECTptr->size()==2){
3458 gen f=g._SYMBptr->feuille._VECTptr->front();
3459 gen i=g._SYMBptr->feuille._VECTptr->back();
3460 int t=cpp_vartype(i,allnames,allchecks);
3461 if (t!=2)
3462 return cprint(f,0,contextptr)+"[cpp_convert_2("+cprint(i,0,contextptr)+",contextptr)]";
3463 else
3464 return cprint(f,0,contextptr)+"["+cprint(i,0,contextptr)+"]";
3465 }
3466 return cprint(g,0,contextptr);
3467 }
3468
3469 // name is the program name if args==program(..),
3470 // name=1 for vectors of instructions or 0 otherwise or -1 for vector to put in a makesequence
3471 // allnames is a list of all variables names that will be converted
3472 // allchecks is the same list but maybe with type returned
cprint(const gen & args,const gen & name_,const vecteur & allnames,const vecteur & allchecks,GIAC_CONTEXT)3473 std::string cprint(const gen & args,const gen & name_,const vecteur & allnames,const vecteur & allchecks,GIAC_CONTEXT){
3474 gen name=name_.is_symb_of_sommet(at_deuxpoints)?name_._SYMBptr->feuille[0]:name_;
3475 if (args.type==_INT_ && args.subtype==_INT_BOOLEAN)
3476 return args.val?"true":"false";
3477 if (args.type==_IDNT){
3478 if (args==cst_pi)
3479 return "cst_pi";
3480 }
3481 if (args.type==_VECT){
3482 string sep=name==1?";\n":",";
3483 string s;
3484 bool ret=!args._VECTptr->empty() && args._VECTptr->back().is_symb_of_sommet(at_return);
3485 if (name==0 && args.subtype!=_SEQ__VECT && !ret){
3486 s="giac::makevecteur(";
3487 }
3488 if (name==-1)
3489 s="giac::makesequence(";
3490 vecteur::const_iterator it=args._VECTptr->begin(), itend=args._VECTptr->end();
3491 if (it==itend)
3492 return "gen(vecteur(0),"+print_INT_(args.subtype)+")";
3493 for(int i=0;;++i){
3494 s += cprint(*it,(name.type==_VECT && i<name._VECTptr->size())?name[i]:zero,allnames,allchecks,contextptr);
3495 ++it;
3496 if (it==itend){
3497 break;
3498 }
3499 s += sep;
3500 }
3501 if (name==0 && args.subtype!=_SEQ__VECT && !ret) s += ")";
3502 if (name==-1) s += ")";
3503 if (name==1 || ret) s += ";";
3504 return s;
3505 }
3506 if (args.type==_CPLX)
3507 return "gen("+args._CPLXptr->print(contextptr)+","+(args._CPLXptr+1)->print(contextptr)+")";
3508 if (args.type==_FUNC){
3509 if (args==at_break)
3510 return "break";
3511 if (args==at_continue)
3512 return "continue";
3513 string name=args.print(contextptr);
3514 if (name.size()>2 && name[0]=='\'' && name[name.size()-1]=='\'')
3515 name=name.substr(1,name.size()-2);
3516 return "at_"+name;
3517 }
3518 if (args.type!=_SYMB)
3519 return args.print(contextptr);
3520 unary_function_ptr & u=args._SYMBptr->sommet;
3521 if (u==at_break)
3522 return "break";
3523 if (u==at_continue)
3524 return "continue";
3525 gen f=args._SYMBptr->feuille;
3526 if ( (u==at_sto || u==at_array_sto) && f.type==_VECT && f._VECTptr->size()==2){
3527 #if 1
3528 if (f[1].type==_VECT){
3529 string res="{\n";
3530 const vecteur & v=*f[1]._VECTptr;
3531 int vs=v.size();
3532 res += "giac::vecteur __tmp_=giac::makevecteur("+cprint(f[0],0,allnames,allchecks,contextptr)+");\n";
3533 for (int j=0;j<vs;++j){
3534 gen var=f[1][j];
3535 int i=cpp_vartype(var,allnames,allchecks);
3536 res += cpp_stoprint(var,allnames,allchecks,contextptr)+" = cpp_convert_"+print_INT_(i)+"("+"__tmp_["+print_INT_(j)+"],contextptr);";
3537 }
3538 return res+"\n}";
3539 }
3540 int i=cpp_vartype(f[1],allnames,allchecks);
3541 int f0t=cpp_vartype(f[0],allnames,allchecks);
3542 if (!i || i==f0t || (i==_DOUBLE_ && f0t==2) || (i==2 && f0t==_DOUBLE_))
3543 return cpp_stoprint(f[1],allnames,allchecks,contextptr)+'='+cprint(f[0],0,allnames,allchecks,contextptr);;
3544 string res= cpp_stoprint(f[1],allnames,allchecks,contextptr)+" = cpp_convert_"+print_INT_(i)+"("+cprint(f[0],0,allnames,allchecks,contextptr)+",contextptr)";
3545 return res;
3546 #else
3547 int i=is_int_or_double(f[1]);
3548 if (!i || is_int_or_double(f[0]))
3549 return f[1].print(contextptr)+'='+cprint(f[0],0,allnames,allchecks,contextptr);;
3550 string res= f[1].print(contextptr)+"=("+cprint(f[0],0,allnames,allchecks,contextptr);
3551 res += (i==_DOUBLE_)?")._DOUBLE_val":").val";
3552 return res;
3553 #endif
3554 }
3555 if (u==at_equal && f.type==_VECT && f._VECTptr->size()==2){
3556 #if 1
3557 int i=cpp_vartype(f[0],allnames,allchecks);
3558 if (!i || i==cpp_vartype(f[1],allnames,allchecks))
3559 return cpp_stoprint(f[0],allnames,allchecks,contextptr)+'='+cprint(f[1],0,allnames,allchecks,contextptr);
3560 string res= cpp_stoprint(f[0],allnames,allchecks,contextptr)+" = cpp_convert_"+print_INT_(i)+"("+cprint(f[1],0,allnames,allchecks,contextptr)+",contextptr)";
3561 return res;
3562 #else
3563 int i=is_int_or_double(f[0]);
3564 if (!i || is_int_or_double(f[1]))
3565 return f[0].print(contextptr)+'='+cprint(f[1],0,allnames,allchecks,contextptr);;
3566 string res= f[0].print(contextptr)+"=("+cprint(f[1],0,allnames,allchecks,contextptr);
3567 res += (i==_DOUBLE_)?")._DOUBLE_val":").val";
3568 return res;
3569 #endif
3570 }
3571 if (u==at_size){
3572 int i=cpp_vartype(f,allnames,allchecks);
3573 if (i==_VECT)
3574 return cprint(f,0,allnames,allchecks,contextptr)+".size()";
3575 return "cpp_convert_7("+cprint(f,0,allnames,allchecks,contextptr)+",contextptr).size()";
3576 }
3577 if (f.type==_VECT && f._VECTptr->size()==2){
3578 if (u==at_at){
3579 int i=cpp_vartype(f._VECTptr->front(),allnames,allchecks);
3580 int j=cpp_vartype(f._VECTptr->back(),allnames,allchecks);
3581 if (i==0 || j==2)
3582 return f._VECTptr->front().print(contextptr)+"["+cprint(f._VECTptr->back(),0,allnames,allchecks,contextptr)+"]";
3583 return f._VECTptr->front().print(contextptr)+"[cpp_convert_2("+cprint(f._VECTptr->back(),0,allnames,allchecks,contextptr)+",contextptr)]";
3584 }
3585 if (u==at_of){
3586 gen fnc=f._VECTptr->front();
3587 if (equalposcomp(allnames,fnc)){
3588 gen evaled_fnc=eval(fnc,1,contextptr);
3589 if (evaled_fnc.is_symb_of_sommet(at_program)){
3590 gen F=evaled_fnc._SYMBptr->feuille;
3591 if (F.type==_VECT && F._VECTptr->size()==3){
3592 vector<int> argtype(gen2vecteur(F[0]).size());
3593 vecteur argF(gen2vecteur(f._VECTptr->back()));
3594 if (argtype.size()==argF.size()){
3595 cprintvars(F[0],false,true,",",&argtype.front(),allnames,allchecks,contextptr);
3596 string res(fnc.print(contextptr));
3597 res += '(';
3598 for (size_t pos=0;pos<argF.size();++pos){
3599 int i=argtype[pos];
3600 int j=cpp_vartype(argF[pos],allnames,allchecks);
3601 if (i==j || (i==0 || j==2))
3602 res += cprint(argF[pos],0,allnames,allchecks,contextptr);
3603 else
3604 res += "cpp_convert_"+print_INT_(i)+"("+cprint(argF[pos],0,allnames,allchecks,contextptr)+",contextptr)";
3605 if (pos!=argF.size()-1)
3606 res += ',';
3607 }
3608 res += ')';
3609 return res;
3610 }
3611 }
3612 }
3613 }
3614 }
3615 }
3616 if (u==at_program && f.type==_VECT && f._VECTptr->size()==3){
3617 // header
3618 bool head=name.type==_IDNT;
3619 string heads;
3620 int convert=0;
3621 if (head){
3622 // IMPROVE if name=="main" make it argv,argc->int
3623 string rtype="giac::gen ";
3624 convert=cpp_vartype(name,allnames,allchecks);
3625 if (name_.is_symb_of_sommet(at_deuxpoints))
3626 convert = name_._SYMBptr->feuille[1].val;
3627 if (convert){
3628 switch (convert){
3629 case _DOUBLE_:
3630 rtype="double ";
3631 break;
3632 case 2:
3633 rtype="long ";
3634 break;
3635 case _VECT:
3636 rtype="giac::vecteur ";
3637 break;
3638 case _STRNG:
3639 rtype="std::string ";
3640 break;
3641 }
3642 }
3643 heads=rtype+name._IDNTptr->id_name+"(";
3644 }
3645 else
3646 heads="giac::gen f(";
3647 vector<int> argtype(gen2vecteur(f[0]).size());
3648 heads += cprintvars(f[0],false,true,",",&argtype.front(),allnames,allchecks,contextptr)+",const giac::context * contextptr=0)";
3649 // core
3650 string core=cprint(f[2],0,allnames,allchecks,contextptr);
3651 if (!core.empty()){
3652 char ch=core[core.size()-1];
3653 if (ch!='}' && ch!=';')
3654 core += ';';
3655 }
3656 string res= heads+"{\n"+core+"\n}\n";
3657 heads="f";
3658 if (head)
3659 heads=name._IDNTptr->id_name;
3660 if (f[0].type!=_VECT || f[0]._VECTptr->size()<=1){
3661 if (f[0].type==_VECT && f[0]._VECTptr->empty())
3662 core ="return "+heads+"(contextptr);";
3663 else {
3664 #if 1
3665 if (0 && name_.is_symb_of_sommet(at_deuxpoints))
3666 core = "return cpp_convert_"+print_INT_(name_._SYMBptr->feuille[1].val)+"("+heads+"(cpp_convert_"+print_INT_(cpp_vartype(f[0]._VECTptr->front(),allnames,allchecks))+"(g,contextptr),contextptr));";
3667 else
3668 core = "return "+heads+"(cpp_convert_"+print_INT_(cpp_vartype(f[0]._VECTptr->front(),allnames,allchecks))+"(g,contextptr),contextptr);";
3669 #else
3670 switch (convert){
3671 case 0:
3672 core="return "+heads+"(g,contextptr);";
3673 break;
3674 case 2:
3675 core="if (g.type!=_INT_) return gensizeerr(contextptr); else return "+heads+"(g.val,contextptr);";
3676 break;
3677 case _DOUBLE_:
3678 core="if (g.type!=_DOUBLE_) return gensizeerr(contextptr); else return "+heads+"(g._DOUBLE_val,contextptr);";
3679 break;
3680 }
3681 #endif
3682 }
3683 }
3684 else {
3685 int nargs=int(f[0]._VECTptr->size());
3686 core = "if (g.type!=_VECT || g.subtype!=_SEQ__VECT || g._VECTptr->size()!="+print_INT_(nargs)+") return gendimerr(contextptr);\n";
3687 core += "vecteur v = *g._VECTptr;";
3688 for (int i=0;i<nargs;++i){
3689 string vi="v["+print_INT_(i)+"]";
3690 switch(argtype[i]){
3691 case _DOUBLE_:
3692 core += vi+"=evalf_double("+vi+",1,contextptr);\n";
3693 core += "if ("+vi+".type!=_DOUBLE_) return gensizeerr(contextptr);\n";
3694 break;
3695 case 2:
3696 core += "if (!is_integral("+vi+")) return gensizeerr(contextptr);\n";
3697 break;
3698 case _VECT:
3699 core += "if ("+vi+".type!=_VECT) return gensizeerr(contextptr);\n";
3700 break;
3701 case _STRNG:
3702 core += "if ("+vi+".type!=_STRNG) return gensizeerr(contextptr);\n";
3703 break;
3704 }
3705 }
3706 core += "return "+heads+"(" ;
3707 for (int i=0;;){
3708 string vi="v["+print_INT_(i)+"]";
3709 switch (argtype[i]){
3710 case _DOUBLE_:
3711 core += vi+"._DOUBLE_val";
3712 break;
3713 case 2:
3714 core += vi+".val";
3715 break;
3716 case _VECT:
3717 core += "*"+vi+"._VECTptr";
3718 break;
3719 case _STRNG:
3720 core += "*"+vi+"._STRNGptr";
3721 break;
3722 default:
3723 core += vi;
3724 }
3725 ++i;
3726 if (i==nargs) break;
3727 core += ",";
3728 }
3729 core += ",contextptr);";
3730 }
3731 heads="giac::gen _"+heads+"(const giac::gen & g,const giac::context * contextptr)";
3732 res += heads+"{\n"+core+"\n}\n";
3733 return res;
3734 }
3735 if (u==at_irem && f.type==_VECT && f._VECTptr->size()==2){
3736 int f0t=cpp_vartype(f[0],allnames,allchecks);
3737 int f1t=cpp_vartype(f[1],allnames,allchecks);
3738 if ( (f0t==2 || f0t==_DOUBLE_ ) && (f1t==2 || f1t==_DOUBLE_))
3739 return "long("+cprint(f[0],0,allnames,allchecks,contextptr)+") % long("+ cprint(f[1],0,allnames,allchecks,contextptr)+")" ;
3740 }
3741 if (u==at_floor){
3742 int ft=cpp_vartype(f,allnames,allchecks);
3743 if (ft==2 || ft==_DOUBLE_)
3744 return "long("+cprint(f,0,allnames,allchecks,contextptr)+")";
3745 }
3746 if (u==at_local){
3747 gen f0=f[0];
3748 string s(cprintvars(f0,false,false,";\n",0,allnames,allchecks,contextptr)+";\n");
3749 return s+cprint(f[1],1,allnames,allchecks,contextptr);
3750 }
3751 if ( (u==at_for || u==at_pour) && f.type==_VECT && f._VECTptr->size()==4){
3752 string res="for(";
3753 res = res +cprint(f[0],0,allnames,allchecks,contextptr);
3754 res = res +';';
3755 res = res +cprint(f[1],0,allnames,allchecks,contextptr);
3756 res = res +';';
3757 res = res +cprint(f[2],0,allnames,allchecks,contextptr);
3758 res = res +"){\n";
3759 res = res +cprint(f[3],0,allnames,allchecks,contextptr);
3760 res = res +";\n}\n";
3761 return res;
3762 }
3763 if (u==at_bloc)
3764 return cprint(f,1,allnames,allchecks,contextptr);
3765 if ( (u==at_ifte || u==at_si) && f.type==_VECT && f._VECTptr->size()==3){
3766 string res="if(";
3767 res = res +cprint(f[0],0,allnames,allchecks,contextptr);
3768 res = res +"){\n";
3769 res = res +cprint(f[1],0,allnames,allchecks,contextptr);
3770 res = res +";\n}";
3771 if (!is_zero(f[2])){
3772 res = res + " else {\n";
3773 res = res + cprint(f[2],0,allnames,allchecks,contextptr) + ";\n}";
3774 }
3775 return res;
3776 }
3777 if (u==at_when){
3778 string res("giac::_when(giac::makesequence(");
3779 res += cprint(f,0,allnames,allchecks,contextptr);
3780 res += "),contextptr)";
3781 return res;
3782 }
3783 if (u==at_try_catch && f.type==_VECT && f._VECTptr->size()>=3){
3784 string res="try {\n";
3785 res = res +cprint(f[0],0,allnames,allchecks,contextptr);
3786 res = res +";\n} catch(";
3787 res = res +cprint(f[1],0,allnames,allchecks,contextptr);
3788 res = res +"){\n";
3789 res = res +cprint(f[2],0,allnames,allchecks,contextptr);
3790 res = res +";\n}";
3791 return res;
3792 }
3793 if (u==at_return){
3794 if (f.type==_VECT)
3795 return "return giac::makevecteur("+cprint(f,0,allnames,allchecks,contextptr)+")";
3796 else
3797 return "return "+cprint(f,0,allnames,allchecks,contextptr);
3798 }
3799 if (u==at_struct_dot){
3800 if (f.type==_VECT && f._VECTptr->size()==2){
3801 gen f0=f._VECTptr->front(),f1=f._VECTptr->back();
3802 if (cpp_vartype(f0,allnames,allchecks)==_VECT){
3803 string f0s=f0.print(contextptr);
3804 gen f1f=f1._SYMBptr->feuille;
3805 if (f1.is_symb_of_sommet(at_append))
3806 return f0s+".push_back("+cprint(f1f,0,allnames,allchecks,contextptr)+")";
3807 if (f1.is_symb_of_sommet(at_prepend))
3808 return f0s+".insert("+f0s+".begin(),"+cprint(f1f,0,allnames,allchecks,contextptr)+")";
3809 if (f1.is_symb_of_sommet(at_suppress)){
3810 f0s=f0s+".erase("+f0s+".begin()+";
3811 if (cpp_vartype(f1f,allnames,allchecks)!=2)
3812 f0s += "cpp_convert_2("+cprint(f1f,0,allnames,allchecks,contextptr)+")";
3813 else
3814 f0s += cprint(f1f,0,allnames,allchecks,contextptr);
3815 return f0s+")";
3816 }
3817 if (f1f.type==_VECT && f1f._VECTptr->empty()){
3818 if (f1.is_symb_of_sommet(at_sort))
3819 return "sort("+f0s+".begin(),"+f0s+".end())";
3820 if (f1.is_symb_of_sommet(at_reverse) || f1.is_symb_of_sommet(at_revlist))
3821 return "reverse("+f0s+".begin(),"+f0s+".end())";
3822 }
3823 }
3824 }
3825 return "// Unable to convert "+args.print(contextptr)+"\n";
3826 }
3827 if (u.ptr()->cprint)
3828 return u.ptr()->cprint(args._SYMBptr->feuille,u.ptr()->s,contextptr);
3829 // operators: if all vars are double or int args.print()
3830 // otherwise or not an operator giac::opname(,contextptr)
3831 string res=u.ptr()->print(contextptr);
3832 int idf0=is_int_or_double(f);
3833 bool idf= idf0 || cpp_vartype(f,allnames,allchecks)==_CPLX;
3834 if (u==at_sin || u==at_cos || u==at_tan || u==at_asin || u==at_acos || u==at_atan || u==at_exp || u==at_ln || u==at_abs){
3835 if (idf)
3836 return "std::"+args.print(contextptr);
3837 }
3838 else
3839 res = '_'+res;
3840 if (u.ptr()->printsommet==printsommetasoperator || u==at_division ||
3841 (idf0 && (u==at_inferieur_strict || u==at_inferieur_egal || u==at_superieur_strict || u==at_superieur_egal))){
3842 int rs=int(res.size());
3843 if ( (idf && (rs!=2 || res[1]!='^')) || (rs==2 && (res[1]=='+' || res[1]=='*' || res[1]=='>' || res[1]=='<')) || (rs==3 && (res[1]=='>' || res[1]=='<') && res[2]=='=') ){
3844 res=u.ptr()->print(contextptr);
3845 string sres;
3846 vecteur v=gen2vecteur(f);
3847 for (int i=0;;){
3848 bool b=need_parenthesis(v[i]);
3849 if (b) sres += '(';
3850 sres += cprint(v[i],0,allnames,allchecks,contextptr);
3851 if (b) sres += ')';
3852 ++i;
3853 if (i>=v.size())
3854 break;
3855 sres += res;
3856 }
3857 return sres;
3858 }
3859 if (rs==2){
3860 switch (res[1]){
3861 case '+':
3862 res="_plus";
3863 break;
3864 case '*':
3865 res="_prod";
3866 break;
3867 case '/':
3868 res="_division";
3869 break;
3870 case '<':
3871 res="_inferieur_strict";
3872 break;
3873 case '>':
3874 res="_superieur_strict";
3875 break;
3876 case '=':
3877 res="_equal";
3878 break;
3879 case '^':
3880 res="_pow";
3881 break;
3882 }
3883 }
3884 }
3885 if (u==at_increment || u==at_decrement){
3886 if (f.type!=_VECT && is_int_or_double(f)){
3887 return cprint(f,0,allnames,allchecks,contextptr)+(u==at_increment?"++":"--");
3888 }
3889 if (f.type==_VECT && f._VECTptr->size()==2 && is_int_or_double(f._VECTptr->front()) && is_int_or_double(f._VECTptr->back()) ){
3890 return cprint(f._VECTptr->front(),0,allnames,allchecks,contextptr)+(u==at_increment?"+=":"-=")+cprint(f._VECTptr->back(),0,allnames,allchecks,contextptr);
3891 }
3892 }
3893 if (u==at_same || u==at_different || u==at_and || u==at_ou ){
3894 return cprint(args._SYMBptr->feuille[0],-1,allnames,allchecks,contextptr)+" "+u.ptr()->s+" "+cprint(args._SYMBptr->feuille[1],-1,allnames,allchecks,contextptr);
3895 }
3896 if (u==at_neg || u==at_not)
3897 return u.ptr()->s+cprint(args._SYMBptr->feuille,-1,allnames,allchecks,contextptr);
3898 if (u==at_inferieur_egal)
3899 return "giac::is_greater("+cprint(args._SYMBptr->feuille[1],-1,allnames,allchecks,contextptr)+","+cprint(args._SYMBptr->feuille[0],0,allnames,allchecks,contextptr)+",contextptr)";
3900 if (u==at_superieur_egal)
3901 return "giac::is_greater("+cprint(args._SYMBptr->feuille[0],-1,allnames,allchecks,contextptr)+","+cprint(args._SYMBptr->feuille[1],0,allnames,allchecks,contextptr)+",contextptr)";
3902 if (u==at_inferieur_strict)
3903 return "!giac::is_greater("+cprint(args._SYMBptr->feuille[0],-1,allnames,allchecks,contextptr)+","+cprint(args._SYMBptr->feuille[1],0,allnames,allchecks,contextptr)+",contextptr)";
3904 if (u==at_superieur_strict)
3905 return "!giac::is_greater("+cprint(args._SYMBptr->feuille[1],-1,allnames,allchecks,contextptr)+","+cprint(args._SYMBptr->feuille[0],0,allnames,allchecks,contextptr)+",contextptr)";
3906 res = "giac::"+res + '(';
3907 res = res + cprint(args._SYMBptr->feuille,-1,allnames,allchecks,contextptr);
3908 res = res +",contextptr)";
3909 return res;
3910 // ...
3911 }
3912
3913 // name is the program name if args==program(..) or 1 for vectors of instructions or 0 otherwise or -1 for vector to put in a makesequence
cprint(const gen & args,const gen & name,GIAC_CONTEXT)3914 std::string cprint(const gen & args,const gen & name,GIAC_CONTEXT){
3915 vecteur v,w;
3916 if (name.is_symb_of_sommet(at_deuxpoints)){
3917 gen f=name._SYMBptr->feuille;
3918 if (f.type==_IDNT){
3919 v.push_back(f);
3920 w.push_back(name);
3921 }
3922 }
3923 else {
3924 if (name.type==_IDNT){
3925 v.push_back(name);
3926 w.push_back(name);
3927 }
3928 }
3929 return cprint(args,name,v,w,contextptr);
3930 }
3931
cprint(const gen & args,GIAC_CONTEXT)3932 std::string cprint(const gen & args,GIAC_CONTEXT){
3933 return cprint(args,0,contextptr);
3934 }
3935
_cprint(const gen & args,GIAC_CONTEXT)3936 gen _cprint(const gen & args,GIAC_CONTEXT){
3937 if ( args.type==_STRNG && args.subtype==-1) return args;
3938 int x=xcas_mode(contextptr),l=language(contextptr);
3939 xcas_mode(0,contextptr);
3940 language(2,contextptr);
3941 gen args_evaled=eval(args,1,contextptr);
3942 string s=cprint(args_evaled,args,contextptr);
3943 xcas_mode(x,contextptr);
3944 language(l,contextptr);
3945 return string2gen(s,false);
3946 }
3947 static const char _cprint_s []="cprint";
3948 static define_unary_function_eval_quoted (__cprint,&_cprint,_cprint_s);
3949 define_unary_function_ptr5( at_cprint ,alias_at_cprint,&__cprint,_QUOTE_ARGUMENTS,true);
3950
3951 #if !defined GIAC_HAS_STO_38 && !defined NSPIRE && !defined FXCG && !defined POCKETCAS
cpp_write_compile(const string & filename,const string & funcname,const string & s,GIAC_CONTEXT)3952 int cpp_write_compile(const string & filename,const string & funcname,const string &s,GIAC_CONTEXT){
3953 ofstream of(filename.c_str());
3954 #ifdef __APPLE__
3955 of << "// -*- mode:c++; compile-command:\" c++ -Wall -fno-strict-aliasing -I/Applications/usr/include -I.. -I. -fPIC -DPIC -g -O2 -dynamiclib giac_" << funcname << ".cpp -o libgiac_" << funcname << ".dylib -L/Applications/usr/lib -L/Applications/usr/64/local/lib -lgiac \" -*-" << '\n';
3956 #else
3957 of << "// -*- mode:c++; compile-command:\" c++ -Wall -fno-strict-aliasing -I.. -I. -fPIC -DPIC -g -O2 -c giac_" << funcname << ".cpp -o giac_" << funcname << ".lo && cc -shared giac_"<<funcname<<".lo -lgiac -lc -Wl,-soname -Wl,libgiac_"<<funcname<<"so.0 -o libgiac_"<<funcname<<".so.0.0.0 && ln -sf libgiac_"<<funcname<<".so.0.0.0 libgiac_"<<funcname<<".so\" -*-" << '\n';
3958 #endif
3959 of << "#include <giac/config.h>" << '\n';
3960 of << "#include <giac/giac.h>" << '\n';
3961 of << "using namespace std;" << '\n';
3962 of << "namespace giac {" << '\n';
3963 of << "inline gen _sqrt(const gen & a,GIAC_CONTEXT){return sqrt(a,contextptr);}" << '\n';
3964 of << "inline gen _arg(const gen & a,GIAC_CONTEXT){return arg(a,contextptr);}" << '\n';
3965 of << "inline gen _complex(const gen & a,GIAC_CONTEXT){return _build_complex(a,contextptr);}" << '\n';
3966 of << "inline gen operator *(const gen & a,double b){return b*a;}" << '\n';
3967 of << "inline bool operator <= (const gen & a,const gen &b){ return is_greater(b,a,giac::context0); }\ninline bool operator >= (const gen & a,const gen &b){ return is_greater(a,b,giac::context0); }" <<'\n';
3968 of << s << '\n';
3969 of << "const string _"+funcname+"_s(\""+funcname+"\");" << '\n';
3970 of << "unary_function_eval __"+funcname+"(0,&_"+funcname+",_"+funcname+"_s);" << '\n';
3971 of << "unary_function_ptr at_"+funcname+"(&__"+funcname+",0,true);" << '\n';
3972 of << "}" << '\n';
3973 of.close();
3974 *logptr(contextptr) << "File " << filename << " created." << '\n';
3975 string cmd="indent -br -brf -l256 giac_"+funcname+".cpp"; int not_ok;
3976 #ifndef __APPLE__
3977 *logptr(contextptr) << "Running " << cmd << '\n';
3978 not_ok=system_no_deprecation(cmd.c_str());
3979 if (not_ok)
3980 *logptr(contextptr) << "Warning, indent not found, please install for nice output" << '\n';
3981 #endif
3982 #ifdef __APPLE__
3983 cmd="g++ -Wall -fno-strict-aliasing -dynamiclib -I/Applications/usr/include -I.. -I. -fPIC -DPIC -g -O2 -dy giac_"+funcname+".cpp -o libgiac_"+funcname+".dylib -L/Applications/usr/lib -L/Applications/usr/64/local/lib -lgiac";
3984 #else
3985 cmd="c++ -Wall -fno-strict-aliasing -I.. -I. -fPIC -DPIC -g -O2 -c giac_"+funcname+".cpp -o giac_"+funcname+".lo";
3986 #endif
3987 *logptr(contextptr) << "Running\n" << cmd << '\n';
3988 not_ok=system_no_deprecation(cmd.c_str());
3989 if (not_ok){
3990 *logptr(contextptr) << "Unable to compile, please fix cpp file" << '\n';
3991 return -1;
3992 }
3993 //cmd="ln -sf giac_"+funcname+".lo giac_"+funcname+".o";
3994 //*logptr(contextptr) << "Running " << cmd << '\n';
3995 //not_ok=system_no_deprecation(cmd.c_str());
3996 #ifndef __APPLE__
3997 cmd="cc -shared giac_"+funcname+".lo -lgiac -lc -Wl,-soname -Wl,libgiac_"+funcname+".so.0 -o libgiac_"+funcname+".so.0.0.0";
3998 *logptr(contextptr) << cmd << '\n';
3999 not_ok=system_no_deprecation(cmd.c_str());
4000 if (not_ok){
4001 *logptr(contextptr) << "Unable to create shared library, perhaps missing libraries?" << '\n';
4002 return -2;
4003 }
4004 //cmd="ln -sf libgiac_"+funcname+".so.0.0.0 libgiac_"+funcname+".so.0";
4005 //system_no_deprecation(cmd.c_str());
4006 cmd="ln -sf libgiac_"+funcname+".so.0.0.0 libgiac_"+funcname+".so";
4007 *logptr(contextptr) << cmd << '\n';
4008 not_ok=system_no_deprecation(cmd.c_str());
4009 #endif
4010 *logptr(contextptr) << "You can now run insmod(\"" << funcname << "\")" << '\n';
4011 return 1;
4012 }
4013
4014 // make a cpp translation file
_cpp(const gen & args,GIAC_CONTEXT)4015 gen _cpp(const gen & args,GIAC_CONTEXT){
4016 if ( args.type==_STRNG && args.subtype==-1) return args;
4017 gen tmp=check_secure();
4018 if (is_undef(tmp)) return tmp;
4019 int x=xcas_mode(contextptr),l=language(contextptr),l0=language(context0),p=python_compat(contextptr);
4020 if (args.type==_IDNT){
4021 xcas_mode(0,contextptr);
4022 language(2,contextptr);
4023 language(2,context0);
4024 python_compat(0,contextptr);
4025 gen argnock=args.is_symb_of_sommet(at_deuxpoints)?args._SYMBptr->feuille[0]:args;
4026 gen args_evaled=eval(argnock,1,contextptr);
4027 //gen args_evaled=eval(args,1,contextptr);
4028 string s=cprint(args_evaled,args,contextptr);
4029 python_compat(p,contextptr);
4030 language(l,contextptr);
4031 language(l0,context0);
4032 xcas_mode(x,contextptr);
4033 string funcname=argnock.print(contextptr);
4034 string filename="giac_"+funcname+".cpp";
4035 return cpp_write_compile(filename,funcname,s,contextptr);
4036 }
4037 if (args.type!=_VECT || args.subtype!=_SEQ__VECT || args._VECTptr->size()<2)
4038 return gensizeerr(contextptr);
4039 xcas_mode(0,contextptr);
4040 language(2,contextptr);
4041 language(2,context0);
4042 python_compat(0,contextptr);
4043 vecteur allnames(*args._VECTptr);
4044 for (int i=0;i<allnames.size();++i){
4045 if (allnames[i].is_symb_of_sommet(at_deuxpoints))
4046 allnames[i]=allnames[i]._SYMBptr->feuille[0];
4047 }
4048 gen arg,args_evaled,argnock;
4049 string s;
4050 for (int i=0;i<args._VECTptr->size();++i){
4051 arg=(*args._VECTptr)[i];
4052 argnock=allnames[i];
4053 args_evaled=eval(argnock,1,contextptr);
4054 if (args_evaled.is_symb_of_sommet(at_program))
4055 s += cprint(args_evaled,arg,allnames,*args._VECTptr,contextptr);
4056 else {
4057 if (args_evaled.type==_DOUBLE_)
4058 s += "const double "+arg.print(contextptr)+"="+args_evaled.print(contextptr)+";";
4059 else {
4060 if (args_evaled.type==_INT_)
4061 s += "const int "+arg.print(contextptr)+"="+args_evaled.print(contextptr)+";";
4062 else
4063 s += "const gen "+arg.print(contextptr)+"(\""+args_evaled.print(contextptr)+"\",context0);";
4064 }
4065 }
4066 }
4067 python_compat(p,contextptr);
4068 language(l,contextptr);
4069 language(l0,context0);
4070 xcas_mode(x,contextptr);
4071 string funcname=argnock.print(contextptr);
4072 string filename="giac_"+funcname+".cpp";
4073 return cpp_write_compile(filename,funcname,s,contextptr);
4074 return 1;
4075 }
4076 static const char _cpp_s []="cpp";
4077 static define_unary_function_eval_quoted (__cpp,&_cpp,_cpp_s);
4078 define_unary_function_ptr5( at_cpp ,alias_at_cpp,&__cpp,_QUOTE_ARGUMENTS,true);
4079 #endif
4080
_hexprint(const gen & g,GIAC_CONTEXT)4081 gen _hexprint(const gen & g,GIAC_CONTEXT){
4082 if ( g.type==_STRNG && g.subtype==-1) return g;
4083 if (g.type == _INT_)
4084 return string2gen(hexa_print_INT_(g.val), false);
4085 if (g.type == _ZINT)
4086 return string2gen(hexa_print_ZINT(*g._ZINTptr), false);
4087 return gentypeerr();
4088 }
4089 static const char _hexprint_s []="hexprint";
4090 static define_unary_function_eval (__hexprint,&_hexprint,_hexprint_s);
4091 define_unary_function_ptr5( at_hexprint ,alias_at_hexprint,&__hexprint,0,true);
4092 static const char _giac_hex_s []="hex";
4093 static define_unary_function_eval (__giac_hex,&_hexprint,_giac_hex_s);
4094 define_unary_function_ptr5( at_giac_hex ,alias_at_giac_hex,&__giac_hex,0,true)
_octprint(const gen & g,GIAC_CONTEXT)4095 gen _octprint(const gen & g,GIAC_CONTEXT){
4096 if ( g.type==_STRNG && g.subtype==-1) return g;
4097 if (g.type == _INT_)
4098 return string2gen(octal_print_INT_(g.val), false);
4099 if (g.type == _ZINT)
4100 return string2gen(octal_print_ZINT(*g._ZINTptr), false);
4101 return gentypeerr();
4102 }
4103 static const char _octprint_s []="octprint";
4104 static define_unary_function_eval (__octprint,&_octprint,_octprint_s);
4105 define_unary_function_ptr5( at_octprint ,alias_at_octprint,&__octprint,0,true);
4106
_binprint(const gen & g,GIAC_CONTEXT)4107 gen _binprint(const gen & g,GIAC_CONTEXT){
4108 if ( g.type==_STRNG && g.subtype==-1) return g;
4109 if (g.type == _INT_)
4110 return string2gen(binary_print_INT_(g.val), false);
4111 if (g.type == _ZINT)
4112 return string2gen(binary_print_ZINT(*g._ZINTptr), false);
4113 return gentypeerr();
4114 }
4115 static const char _binprint_s[]="binprint";
4116 static define_unary_function_eval(__binprint,&_binprint,_binprint_s);
4117 define_unary_function_ptr5(at_binprint,alias_at_binprint,&__binprint,0,true);
4118 // const unary_function_ptr at_binprint (&__binprint,0,true);
4119
4120 static const char _giac_bin_s []="bin";
4121 static define_unary_function_eval (__giac_bin,&_binprint,_giac_bin_s);
4122 define_unary_function_ptr5( at_giac_bin ,alias_at_giac_bin,&__giac_bin,0,true);
4123
4124 #ifndef NO_NAMESPACE_GIAC
4125 } // namespace giac
4126 #endif // ndef NO_NAMESPACE_GIAC
4127