@@ -191,7 +191,12 @@ local str, msolve_path, fname1, fname2, file_dir, verb, param, nthreads, output,
191191 fname2 := cat(file_dir, RandomTools[Generate](string(8,alpha)), ".ms");;
192192 end if;
193193
194- param:=0;
194+ str:=subs(opts,"param");
195+ if type(str, integer) and str in [0, 1, 2] then
196+ param:=str;
197+ else
198+ param:=0:
199+ end if;
195200 msolve_path, fname1, fname2, file_dir, verb, param, nthreads, output, gb, elim;
196201end proc;
197202
@@ -297,7 +302,7 @@ end proc:
297302# [ vars[1] = (a1+b1)/2, ..., vars[n] = (an+bn)/2 ]
298303MSolveRealRoots:=proc(F, vars, opts:={})
299304local results, dim, fname1, fname2, verb, param, msolve_path, file_dir,
300- lsols, nl, i, j, gb, output, nthreads, str, sols, prec, elim;
305+ lsols, nl, i, j, gb, output, nthreads, str, sols, prec, elim, b, ratpar ;
301306 if type(F, list(polynom(rational))) = false then
302307 printf("First argument is not a list of polynomials with rational coefficients\n");
303308 end if;
@@ -328,7 +333,6 @@ lsols, nl, i, j, gb, output, nthreads, str, sols, prec, elim;
328333 fi:
329334 str := cat(msolve_path, " -v ", verb, " -P ", param, " -p ", prec, " -t ", nthreads, " -f ", fname1, " -o ", fname2):
330335 gb:=0; #Needed to avoid the user stops GB comp once a prime computation is done
331- param:=0;
332336 try
333337 system(str):
334338 read(fname2):
@@ -345,6 +349,22 @@ lsols, nl, i, j, gb, output, nthreads, str, sols, prec, elim;
345349 return [];
346350 end if;
347351
352+ if param in [1, 2] then
353+ b, ratpar := GetRootsFromMSolve(results[2]);
354+ if b = -1 then
355+ printf("System has infinitely many complex solutions\n");
356+ return [];
357+ elif b = -2 then
358+ printf("System not in generic position. You may add to your system\n");
359+ printf("a random linear form of your variables and a new variable\n");
360+ return [];
361+ end if;
362+ end if;
363+
364+ if param = 2 then
365+ return ratpar;
366+ end if;
367+
348368 dim := results[1];
349369 if dim = -1 then
350370 if verb >= 1 then
@@ -359,7 +379,11 @@ lsols, nl, i, j, gb, output, nthreads, str, sols, prec, elim;
359379 return [1];
360380 end if;
361381 if dim = 0 then
362- lsols := results[2];
382+ if param = 1 then
383+ lsols := [op(results[3][1]), map(s -> subs(zip((k,v) -> k = v, results[2][4], s), vars), results[3][2])];
384+ else
385+ lsols := results[2];
386+ end if;
363387 nl := lsols[1]:
364388 sols:=[]:
365389 for i from 1 to nl do
@@ -368,6 +392,9 @@ lsols, nl, i, j, gb, output, nthreads, str, sols, prec, elim;
368392 if output=1 then
369393 sols := map(_p->map(_c->lhs(_c)=(rhs(_c)[1]+rhs(_c)[2])/2, _p), sols);
370394 end if;
395+ if param = 1 then
396+ return ratpar, [0, sols];
397+ end if;
371398 return [0, sols];
372399 end if;
373400 return results;
@@ -388,6 +415,6 @@ end proc:
388415# #Usage
389416
390417# #with rational parametrization
391- # param, sols:=MSolveRealRoots(F,vars," ../binary/msolve","/tmp/ in.ms","/tmp/ out.ms",1 );
418+ # param, sols:=MSolveRealRoots(F,vars,{"mspath"=" ../binary/msolve","file_in"=" in.ms","file_out"=" out.ms","param"=1} );
392419# #or with solutions only
393- # sols:=MSolveRealRoots(F,vars);
420+ # sols:=MSolveRealRoots(F,vars,{"mspath"="../binary/msolve"} );
0 commit comments