File: [local] / OpenXM / src / asir-contrib / packages / src / sm1.rr (download)
Revision 1.18, Wed Mar 30 04:45:14 2022 UTC (2 years, 2 months ago) by takayama
Branch: MAIN
CVS Tags: HEAD
Changes since 1.17: +6 -1
lines
Added an explanation on the behavior of update=4.
The option update=4 installs cygwin binary distribution (cygwin executables, ox servers and libraries including ox_pari.exe).
It rollbacks the asir-contrib library to the cygwin binary distribution. Then, you need to do update=1 if you want the latest asir-contrib.
OpenXM-win may be moved under c:/ProgramData/asir/OpenXM by hand.
|
/*$OpenXM: OpenXM/src/asir-contrib/packages/src/sm1.rr,v 1.18 2022/03/30 04:45:14 takayama Exp $ */
/* Old: sm1, see Attic */
module sm1;
/* ------------ list of local functions ---------- */
localf get_Sm1_proc$
localf set_Sm1_proc$
localf find_proc$
localf check_server$
localf start$
localf start_unix$
localf start_windows$
localf sm1flush$
localf sm1push$
localf sm1$
localf call_sm1$
localf sm1pop$
localf to_asir_form$
localf toAsirForm$
localf toOrdered$
localf toOrdered_dehomo$
localf push_poly0_R$
localf push_poly0$
localf pop_poly0$
localf pop_poly0_0_R$
localf pop_poly0_0$
localf push_int0_R$
localf push_int0$
localf push_0_R$
localf push_0$
localf push$
localf pop$
localf pop2$
localf check_arg_gb$
localf isListOfPoly$
localf isListOfVar$
localf gb$
localf gb_d$
localf pgb$
localf deRham$
localf vlist$
localf ringD$
localf expand_d$
localf mul_d$
localf dehomogenize_d$
localf homogenize_d$
localf groebner_d$
localf reduction_d$
localf reduction_noH_d$
localf hilbert$
localf genericAnn$
localf tensor0$
localf wTensor0$
localf reduction$
localf reduction_noH$
localf xml_tree_to_prefix_string$
localf wbf$
localf wbfRoots$
localf res_div$
localf syz$
localf res_solv$
localf res_solv_h$
localf mul$
localf mul_h$
localf adjoint$
localf transpose$
localf resol1$
localf gcd_aux$
localf lcm_aux$
localf mul_v$
localf div_v$
localf rat_to_p_aux$
localf denom_aux0$
localf num_aux0$
localf rat_to_p$
localf distraction$
localf sm1l_ntoint32$
localf to_ascii_array$
localf from_ascii_array$
localf gkz$
localf mgkz$
localf appell1$
localf sm1aux_dx$
localf sm1aux_x$
localf appell4$
localf rank$
localf rrank$
localf auto_reduce$
localf slope$
localf restriction$
localf integration$
localf ahg$
localf bfunction$
localf generalized_bfunction$
localf saturation$
localf ecartd_gb$
localf ecartd_reduction$
localf gb_reduction$
localf gb_reduction_noh$
localf gb_oxRingStructure$
localf ecartd_isSameIdeal_h$
localf ecart_homogenize01Ideal$
localf aux_parse_gb_option_list$
localf ecartd_gb_oxRingStructure$
localf isSameIdeal_in_Dalg$
localf ecartd_reduction_noh$
localf sm1vstr_to_vlist $
localf sm1vlist_to_asirvlist $
localf gb_ptod2 $
localf gb_ptod $
localf gb_dtop $
localf is_in_dp $
localf gb_dtosm1p $
localf gb_dehomo2 $
localf gb_dehomo $
localf ecartd_syz $
localf init_Sm1_proc $
localf laplace $
localf add_d$
localf regular_triangulation$
localf mod_reduction$
localf reduction_verbose$
/* ------------ static variables ---------------- */
static Sm1_proc$
def init_Sm1_proc() {
Sm1_proc = -1$
}
init_Sm1_proc() $
static V_sm1_pop$
/* Functions to get and set static variables out of the module */
def get_Sm1_proc() {
return Sm1_proc;
}
def set_Sm1_proc(A) {
Sm1_proc = A;
}
#define SM1L_FIND_PROC(P) P = getopt(proc);\
if (type(P) == -1) {\
P = find_proc();\
}
def find_proc() {
/*! extern Sm1_proc; */
if (Sm1_proc == -1) {
A = ox_get_serverinfo();
/* Look for ox_sm1. Not yet written */
/* Start sm1 automatically if there is not ox_sm1 */
Sm1_proc = start();
}
return(Sm1_proc);
}
/* Search : oxPrintMessage, see cmoDebugCMO, too. */
/************** end of configure *******************************/
def check_server(P) {
M=ox_get_serverinfo(P);
if (M == []) {
return(start());
}
if (M[0][1] != "Ox_system=ox_sm1_ox_sm1_forAsir") {
print("Warning: the server number ",0)$
print(P,0)$
print(" is not ox_sm1_forAsir server.")$
print("Starting ox_sm1_forAsir server on the localhost.")$
return(start());
}
return(P);
}
def start() {
extern Xm_unix;
if (ox_ostype()[0] == "windows" && Xm_unix == 0)
return start_windows(0);
else
return start_unix();
}
def start_unix() {
extern Sm1_lib;
extern Xm_noX;
/*! extern Sm1_proc; */
if (Xm_noX) {
P = ox_launch_nox(0,Sm1_lib+"/bin/ox_sm1_forAsir");
}else{
P = ox_launch(0,Sm1_lib+"/bin/ox_sm1_forAsir");
}
if (Xm_noX) {
sm1(P," oxNoX ");
}
ox_check_errors(P);
Sm1_proc = P;
return(P);
}
def start_windows(U) {
if (U) error("start_windows is called.");
Monitor = " " ;
// Monitor = " -monitor ";
/* if (sysinfo()[0] != "windows") error("This is not Windows"); */
if (type(getopt(openxm_home_win)) > 0) {
OpenXM_HOME_WIN=rtostr(getopt(openxm_home_win));
}else if (access(getenv("APPDATA")+"/OpenXM/OpenXM-win")) {
OpenXM_HOME_WIN=getenv("APPDATA")+"/OpenXM/OpenXM-win";
/* NOTE!! No space or damemoji are allowed in the user name. */
}else if (access("c:/ProgramData/asir/OpenXM")) {
OpenXM_HOME_WIN="c:/ProgramData/asir/OpenXM/OpenXM-win";
}else{
OpenXM_HOME_WIN=getenv("APPDATA")+"/OpenXM/OpenXM-win";
/* NOTE!! No space or damemoji are allowed in the user name. */
}
printf("OpenXM_HOME_WIN=%a\n",OpenXM_HOME_WIN);
printf(" NOTE!! No space or damemoji are allowed in the OpenXM_HOME_WIN\n");
printf("To set manually call pari.start_win(|openxm_home_win=\"...\");\n");
CPort = generate_port();
SPort = generate_port();
// Ox_sm1_forAsir = OpenXM_HOME_WIN+"/bin/ox_sm1" ;
Ox_sm1_forAsir = OpenXM_HOME_WIN+"/bin/ox_sm1_forAsir" ;
OpenXM_START_EXE="start ";
Xm_use_timer_to_start_server=1;
if (version() > 20160828) {
Com = OpenXM_HOME_WIN+"/bin/ox -ox ";
}else{
Com = OpenXM_HOME_WIN+"/bin/ox -protocol_1999 -ox ";
}
Com += Ox_sm1_forAsir + " " + Monitor +
" -data "+ rtostr(SPort) +" -control "+ rtostr(CPort);
Com = OpenXM_START_EXE+" /MIN "+Com ;
print(Com +"\n");
// 2016.08.26 generate start-up batch file.
BCom = "set OpenXM_HOME="+OpenXM_HOME_WIN+"\n";
BCom += "echo %OpenXM_HOME%\n";
BCom += Com;
BCom += "\nexit\n";
BatchFile=OpenXM_HOME_WIN+"/bin/tmp-oxe";
Fp = open_file(BatchFile+".bat","w");
fprintf(Fp,BCom);
close_file(Fp);
shell("start /MIN "+BatchFile);
// Todo, remove the batch file.
//shell(Com);
if (Xm_use_timer_to_start_server) {
print("Waiting for 3 seconds."); sleep(3000);
}else{
purge_stdin(); print("Type in Return to connect to the server.");
get_line();
}
print("Trying to connect to the server...",0)$
CSocket = try_connect("localhost",CPort);
SSocket = try_connect("localhost",SPort);
P = register_server(CSocket,CPort,SSocket,SPort);
print(" Done.");
if (Xm_noX) {
sm1(P," oxNoX ");
}
ox_check_errors(P);
Sm1_proc = P;
return(P);
}
/* ox_sm1 */
/* P is the process number */
def sm1flush(P) {
ox_execute_string(P,"[(flush)] extension pop");
}
def sm1push(P,F) {
G = ox_ptod(F);
ox_push_cmo(P,G);
}
def sm1(P,F) {
ox_execute_string(P,F);
sm1flush(P);
}
/*&usage begin: sm1.call_sm1(F)
It executes {F} on the sm1 server.
See also sm1.
end: */
def call_sm1(F) {
SM1L_FIND_PROC(P);
P = check_server(P);
ox_execute_string(P,F);
sm1flush(P);
}
def sm1pop(P) {
return(ox_pop_cmo(P));
}
def to_asir_form(V) { return(toAsirForm(V)); }
def toAsirForm(V) {
extern ToAsirForm_V; /* for debug */
if (type(V) == 4) { /* list */
if((length(V) == 3) && (V[0] == "sm1_dp")) {
/* For debugging. */
if (ToAsir_Debug != 0) {
ToAsirForm_V = V;
print(map(type,V[1]));
print(V);
}
/* */
Vlist = map(strtov,V[1]);
return(dp_dtop(V[2],Vlist));
} else {
return(map(toAsirForm,V));
}
}else{
return(V);
}
}
def toOrdered(V) {
Dehomo = getopt(dehomogenize);
if (type(Dehomo) == -1) {
Dehomo = 0;
}else{
Dehomo = 1;
}
if (type(V) == 4) { /* list */
if((length(V) == 3) && (V[0] == "sm1_dp")) {
Vlist = map(strtov,V[1]);
Ans = "";
F = V[2];
while (F != 0) {
G = dp_hm(F);
F = dp_rest(F);
if (dp_hc(G)>0) {
Ans += "+";
}
if (Dehomo) {
Ans += rtostr(subst(dp_dtop(G,Vlist),h,1));
}else{
Ans += rtostr(dp_dtop(G,Vlist));
}
}
return Ans;
} else {
if (Dehomo) {
return(map(toOrdered_dehomo,V));
}else{
return(map(toOrdered,V));
}
}
}else{
return(V);
}
}
def toOrdered_dehomo(F) {
return(toOrdered(F|dehomogenize=1));
}
def push_poly0_R(A,P,Vlist) {
return(push_poly0(P,A,Vlist));
}
def push_poly0(P,A,Vlist) {
if (type(Vlist[0]) == 4) {
Vlist = Vlist[2];
}
/* if Vlist=[[e,x,y,H,E,Dx,Dy,h],[e,x,y,hH,eE,dx,dy,h],[e,x,y,hH,eE,dx,dy,h]]
list of str (sm1) list of str (asir) list of var (asir)
then we execute the code above.
*/
if (type(A) == 2 || type(A) == 1) { /* recursive poly or number*/
A = dp_ptod(A,Vlist);
ox_push_cmo(P,A);
return;
}
if (type(A) == 0) { /* zero */
sm1(P," (0). ");
return;
}
if (type(A) == 4) { /* list */
ox_execute_string(P," [ ");
map(push_poly0_R,A,P,Vlist);
ox_execute_string(P," ] ");
return;
}
ox_push_cmo(P,A);
ox_check_errors2(P);
return;
}
/* push_poly0(0,[0,1,x+y,["Hello",y^3]],[x,y]); */
def pop_poly0(P,Vlist) {
if (type(Vlist[0]) == 4) {
Vlist = Vlist[2];
}
A = ox_pop_cmo(P);
return(pop_poly0_0(P,A,Vlist));
}
def pop_poly0_0_R(A,P,Vlist) {
return(pop_poly0_0(P,A,Vlist));
}
def pop_poly0_0(P,A,Vlist) {
if (type(A) == 4) {
return(map(pop_poly0_0_R,A,P,Vlist));
}
if (type(A)== 9) {return(dp_dtop(A,Vlist));}
return(A);
}
def push_int0_R(A,P) {
return(push_int0(P,A));
}
def push_int0(P,A) {
if (type(A) == 1 || type(A) == 0) {
/* recursive poly or number or 0*/
A = rtostr(A);
ox_push_cmo(P,A);
sm1(P," . (integer) dc ");
return;
}
if (type(A) == 2) {
A = rtostr(A); ox_push_cmo(P,A);
return;
}
if (type(A) == 4) { /* list */
ox_execute_string(P," [ ");
map(push_int0_R,A,P);
ox_execute_string(P," ] ");
return;
}
ox_push_cmo(P,A);
return;
}
def push_0_R(A,P) {
return(push_0(P,A));
}
def push_0(P,A) {
if (type(A) == 0) {
/* 0 */
A = rtostr(A);
ox_push_cmo(P,A);
sm1(P," .. ");
return;
}
if (type(A) == 2) {
/* Vlist = vars(A); One should check Vlist is a subset of Vlist3. */
Vlist2 = vlist(P);
Vlist3 = map(strtov,Vlist2[1]);
B = dp_ptod(A,Vlist3);
ox_push_cmo(P,B);
return;
}
if (type(A) == 4) { /* list */
ox_execute_string(P," [ ");
map(push_0_R,A,P);
ox_execute_string(P," ] ");
return;
}
ox_push_cmo(P,A);
return;
}
def push(P,A) {
push_0(P,A);
}
def pop(P) {
/*! extern V_sm1_pop; */
sm1(P," toAsirForm ");
V_sm1_pop = ox_pop_cmo(P);
return(toAsirForm(V_sm1_pop));
}
def pop2(P) {
/*! extern V_sm1_pop; */
sm1(P," toAsirForm ");
V_sm1_pop = ox_pop_cmo(P);
return([toAsirForm(V_sm1_pop),V_sm1_pop]);
}
def check_arg_gb(A,Fname) {
/* A = [[x^2+y^2-1,x*y],[x,y],[[x,-1,y,-1]]] */
if (type(A) != 4) {
error(Fname+" : argument should be a list.");
}
if (length(A) < 2) {
error(Fname+" : argument should be a list of 2 or 3 elements.");
}
if (type(A[0]) != 4) {
error(Fname+" : example: [[dx^2+dy^2-4,dx*dy-1]<== it should be a list,[x,y]]");
}
if (!isListOfPoly(A[0])) {
error(Fname+" : example: [[dx^2+dy^2-4,dx*dy-1]<== it should be a list of polynomials or strings,[x,y]]");
}
if (!isListOfVar(A[1])) {
error(Fname+" : example: [[dx^2+dy^2-4,dx*dy-1],[x,y]<== list of variables or \"x,y\"]");
}
if (length(A) >= 3) {
if (type(A[2]) != 4) {
error(Fname+" : example:[[dx^2+dy^2-4,dx*dy-1],[x,y],[[x,-1,dx,1]]<== a list of weights]");
}
if (type(A[2][0]) != 4) {
error(Fname+" : example:[[dx^2+dy^2-4,dx*dy-1],[x,y],[[x,-1,dx,1],[dy,1]]<== a list of lists of weight]");
}
}
return(1);
}
def isListOfPoly(A) {
if (type(A) !=4 ) return(0);
N = length(A);
for (I=0; I<N; I++) {
if (type(A[I]) == 4) {
if (!isListOfPoly(A[I])) return(0);
} else if (!(type(A[I]) == 0 || type(A[I]) == 1 || type(A[I]) == 2 ||
type(A[I]) == 7 || type(A[I]) == 9)) {
return(0);
}
}
return(1);
}
def isListOfVar(A) {
if (type(A) == 7) return(1); /* "x,y" */
if (type(A) != 4) return(0);
N = length(A);
for (I=0; I<N; I++) {
if (!(type(A[I]) == 2 || type(A[I]) == 7 )) {
return(0);
}
}
return(1);
}
def gb(A) {
SM1L_FIND_PROC(P);
Dp = is_in_dp(getopt(),A);
A=Dp[1]; Dp=Dp[0];
Sorted = getopt(sorted);
if (type(getopt(ringvar))>0) {
RingVar=rtostr(getopt(ringvar));
}else{
RingVar="ring_var_for_asir";
}
if (type(Sorted) == -1) {
Sorted = 0;
}else{
Sorted = 1;
}
Dehomo = getopt(dehomogenize);
if (type(Dehomo) == -1) {
Dehomo = 0;
}else{
Dehomo = 1;
}
P = check_server(P);
if ((length(A[0]) >= 1) && (type(A[0][0]) != type(""))) A = cons(rat_to_p(A[0])[0],cdr(A));
check_arg_gb(A,"Error in sm1.gb");
push_int0(P,A);
if (type(getopt(needBack))>0) {
sm1(P," [(needBack) 1] setAttributeList ");
}
sm1(P," gb /tmpAsirSm1Gb set ");
sm1(P,sprintf(" tmpAsirSm1Gb getRing /%a set ",RingVar));
sm1(P," tmpAsirSm1Gb ");
if (Sorted) {
T = pop2(P);
if (Dehomo) {
return(append(T[0],[toOrdered(T[1]|dehomogenize=1)]));
}else{
return(append(T[0],[toOrdered(T[1])]));
}
}else{
Ans = pop(P);
if (Dp) {
Opt = gb_oxRingStructure();
Ans = gb_ptod(Ans,Opt);
}
if (type(getopt(needBack))>0) {
sm1(P," tmpAsirSm1Gb getAttributeList ");
All = pop(P);
return append(Ans,All);
}
return Ans;
}
}
def gb_d(A) {
SM1L_FIND_PROC(P);
P = check_server(P);
check_arg_gb(A,"Error in sm1.gb_d");
push_int0(P,A);
sm1(P," gb /gb.tmp1 set ");
sm1(P," gb.tmp1 getOrderMatrix {{(universalNumber) dc} map } map /gb.tmp2 set ");
sm1(P," gb.tmp1 0 get 0 get getvNamesCR { [(class) (indeterminate)] dc } map /gb.tmp3 set ");
sm1(P," gb.tmp1 getRing ring_def "); /* Change the current ring! */
sm1(P,"[[ gb.tmp3 gb.tmp2] gb.tmp1] ");
return(ox_pop_cmo(P));
}
def pgb(A) {
SM1L_FIND_PROC(P);
P = check_server(P);
check_arg_gb(A,"Error in pgb");
sm1(P," set_timer ");
push_int0(P,A);
sm1(P," pgb ");
B = pop(P);
sm1(P," set_timer ");
return(B);
}
def deRham(A) {
SM1L_FIND_PROC(P);
P = check_server(P);
sm1(P," set_timer ");
push_int0(P,A);
sm1(P," deRham ");
B = pop(P);
sm1(P," set_timer ");
ox_check_errors2(P);
return(B);
}
def vlist(P) {
sm1(P," getvNamesC ");
B=ox_pop_cmo(P);
sm1(P," getvNamesC toAsirVar ");
C=ox_pop_cmo(P);
return([B,C,map(strtov,C)]);
}
/* [ sm1 names(string), asir names(string), asir names(var)] */
/* Vlist = vlist(P);
push_poly0( x + 20*x, Vlist[2]);
pop_poly0(Vlist[2]);
*/
/* ring of Differential operators */
def ringD(V,W) {
SM1L_FIND_PROC(P);
sm1(P," [ ");
if (type(V) == 7) { /* string */
ox_push_cmo(P,V);
}else if (type(V) == 4) {/* list */
V = map(rtostr,V);
ox_push_cmo(P,V);
sm1(P," from_records ");
}else { printf("Error: ringD"); return(-1); }
sm1(P," ring_of_differential_operators ");
if (type(W) != 0) {
push_int0(P,W); sm1(P," weight_vector ");
}
sm1(P," pstack ");
sm1(P," 0 ] define_ring getOrderMatrix {{(universalNumber) dc}map}map ");
ox_check_errors2(P);
M = ox_pop_cmo(P);
return([vlist(P)[2],M]);
}
def expand_d(F) {
SM1L_FIND_PROC(P);
ox_push_cmo(P,F);
sm1(P, " expand ");
return(ox_pop_cmo(P));
}
def mul_d(A,B) {
SM1L_FIND_PROC(P);
ox_push_cmo(P,A);
ox_push_cmo(P,B);
sm1(P," mul ");
return(ox_pop_cmo(P));
}
def dehomogenize_d(A) {
SM1L_FIND_PROC(P);
ox_push_cmo(P,A);
sm1(P," dehomogenize ");
return(ox_pop_cmo(P));
}
def homogenize_d(A) {
SM1L_FIND_PROC(P);
ox_push_cmo(P,A);
sm1(P," homogenize ");
return(ox_pop_cmo(P));
}
def groebner_d(A) {
SM1L_FIND_PROC(P);
ox_push_cmo(P,A);
sm1(P," groebner ");
return(ox_pop_cmo(P));
}
def reduction_d(F,G) {
SM1L_FIND_PROC(P);
ox_push_cmo(P,F);
ox_push_cmo(P,G);
sm1(P," reduction ");
Ans = ox_pop_cmo(P);
return([Ans[0],Ans[1],Ans[2],G]);
}
def reduction_noH_d(F,G) {
SM1L_FIND_PROC(P);
ox_push_cmo(P,F);
ox_push_cmo(P,G);
sm1(P," reduction-noH ");
Ans = ox_pop_cmo(P);
return([Ans[0],Ans[1],Ans[2],G]);
}
def hilbert(A) {
G=nd_gr(A[0],A[1],0,0); /* bug fix for duplicated generators. misc-2021/11/misc/sm1-hilb-bug.rr */
A=[G,A[1]];
SM1L_FIND_PROC(P);
P = check_server(P);
sm1(P,"[ ");
push_int0(P,A[0]);
push_int0(P,A[1]);
sm1(P," ] pgb /sm1_hilbert.gb set ");
sm1(P," sm1_hilbert.gb 0 get { init toString } map ");
push_int0(P,A[1]);
sm1(P, " hilbert ");
B = pop(P);
return(B[1]/fac(B[0]));
}
def genericAnn(F) {
SM1L_FIND_PROC(P);
push_int0(P,F[0]);
push_int0(P,F[1]);
sm1(P, " genericAnn ");
B = pop(P);
return(B);
}
def tensor0(F) {
SM1L_FIND_PROC(P);
push_int0(P,F);
sm1(P, " tensor0 ");
B = pop(P);
return(B);
}
def wTensor0(F) {
SM1L_FIND_PROC(P);
push_int0(P,F);
sm1(P, " wTensor0 ");
B = pop(P);
return(B);
}
def reduction(A) {
/* Example: reduction(A|proc=10) */
SM1L_FIND_PROC(P);
/* check the arguments */
if (type(A) != 4) {
error("sm1. reduction(A|proc=p): A must be a list.");
}
if (length(A) == 2) return mod_reduction(A | option_list=getopt());
AA = [rtostr(A[0])];
AA = append(AA,[ map(rtostr,A[1]) ]);
AA = append(AA, cdr(cdr(A)));
sm1(P," /reduction*.noH 0 def ");
push_int0(P,AA);
sm1(P," reduction* ");
ox_check_errors2(P);
Ans = pop(P); /* Throw away Ans[3] */
return([Ans[0],Ans[1],Ans[2],A[1]]);
}
def reduction_noH(A) {
/* Example: reduction(A|proc=10) */
SM1L_FIND_PROC(P);
/* check the arguments */
if (type(A) != 4) {
error("sm1. reduction_noH(A|proc=p): A must be a list.");
}
AA = [rtostr(A[0])];
AA = append(AA,[ map(rtostr,A[1]) ]);
AA = append(AA, cdr(cdr(A)));
sm1(P," /reduction*.noH 1 def ");
push_int0(P,AA);
sm1(P," reduction* ");
ox_check_errors2(P);
Ans = pop(P); /* Throw away Ans[3] */
return([Ans[0],Ans[1],Ans[2],A[1]]);
}
def xml_tree_to_prefix_string(A) {
SM1L_FIND_PROC(P);
/* check the arguments */
if (type(A) != 7) {
error("sm1. xml_tree_to_prefix_string(A|proc=p): A must be a string.");
}
ox_push_cmo(P,A);
sm1(P," xml_tree_to_prefix_string ");
ox_check_errors2(P);
return(ox_pop_cmo(P));
}
def wbf(A) {
SM1L_FIND_PROC(P);
/* check the arguments */
if (type(A) != 4) {
error("sm1. wbf(A): A must be a list.");
}
if (length(A) != 3) {
error("sm1. wbf(A): A must be a list of the length 3.");
}
if (type(A[0]) != 4 || type(A[1]) != 4 || type(A[2]) != 4) {
error("sm1. wbf([A,B,C]): A, B, C must be a list.");
}
if (! (type(A[2][0]) == 7 || type(A[2][0]) == 2)) {
error("sm1. wbf([A,B,C]): C must be of a form [v-name, v-weight, ...]");
}
push_int0(P,A);
sm1(P," wbf ");
ox_check_errors2(P);
return(pop(P));
}
def wbfRoots(A) {
SM1L_FIND_PROC(P);
/* check the arguments */
if (type(A) != 4) {
error("sm1. wbfRoots(A): A must be a list.");
}
if (length(A) != 3) {
error("sm1. wbfRoots(A): A must be a list of the length 3.");
}
if (type(A[0]) != 4 || type(A[1]) != 4 || type(A[2]) != 4) {
error("sm1. wbfRoots([A,B,C]): A, B, C must be a list.");
}
if (! (type(A[2][0]) == 7 || type(A[2][0]) == 2)) {
error("sm1. wbfRoots([A,B,C]): C must be of a form [v-name, v-weight, ...]");
}
push_int0(P,A);
sm1(P," wbfRoots ");
ox_check_errors2(P);
return(pop(P));
}
def res_div(A) {
SM1L_FIND_PROC(P);
push_int0(P,[[A[0],A[1]],A[2]]);
sm1(P," res*div ");
ox_check_errors2(P);
return(pop(P));
}
def syz(A) {
SM1L_FIND_PROC(P);
Dp = is_in_dp(getopt(),A);
A=Dp[1]; Dp=Dp[0];
push_int0(P,A);
sm1(P," syz ");
ox_check_errors2(P);
Ans = pop(P);
if (Dp) {
Opt = gb_oxRingStructure();
Ans = gb_ptod(Ans,Opt);
}
return Ans;
}
def res_solv(A) {
SM1L_FIND_PROC(P);
push_int0(P,[[A[0],A[1]],A[2]]);
sm1(P," res*solv ");
ox_check_errors2(P);
return(pop(P));
}
def res_solv_h(A) {
SM1L_FIND_PROC(P);
push_int0(P,[[A[0],A[1]],A[2]]);
sm1(P," res*solv*h ");
ox_check_errors2(P);
return(pop(P));
}
def mul(A,B,V) {
SM1L_FIND_PROC(P);
/* For rational coeff test: sm1.mul(dx/2,x/3,[x]);
sm1.mul(dx/a,x,[x,a]);
*/
AA = nm(A); BB= nm(B);
if (dn(B) != 1) error("sm1.mul: not implemented yet.");
D = 1/(dn(A)*dn(B));
AA = ptozp(AA | factor=1);
BB = ptozp(BB | factor=1);
D = AA[1]*BB[1]*D;
push_int0(P,[[AA[0],BB[0]],V]);
sm1(P," res*mul ");
ox_check_errors2(P);
return(D*pop(P));
}
def mul_h(A,B,V) {
SM1L_FIND_PROC(P);
/* For rational coeff test: sm1.mul(dx/2,x/3,[x]);
sm1.mul(dx/a,x,[x,a]);
*/
AA = nm(A); BB= nm(B);
if (dn(B) != 1) error("sm1.mul: not implemented yet.");
D = 1/(dn(A)*dn(B));
AA = ptozp(AA | factor=1);
BB = ptozp(BB | factor=1);
D = AA[1]*BB[1]*D;
push_int0(P,[[A[0],B[0]],V]);
sm1(P," res*mul*h ");
ox_check_errors2(P);
return(D*pop(P));
}
def adjoint(A,V) {
SM1L_FIND_PROC(P);
push_int0(P,[A,V]);
sm1(P," res*adjoint ");
ox_check_errors2(P);
return(pop(P));
}
def transpose(A) {
if (type(A) == 4) {
N = length(A); M = length(A[0]);
B = newmat(N,M,A);
C = newmat(M,N);
for (I=0; I<N; I++) {
for (J=0; J<M; J++) {
C[J][I] = B[I][J];
}
}
D = newvect(M);
for (J=0; J<M; J++) {
D[J] = C[J];
}
return(map(vtol,vtol(D)));
}else{
print(A)$
error("sm1. tranpose: traspose for this argument has not been implemented.");
}
}
def resol1(A) {
SM1L_FIND_PROC(P);
push_int0(P,A);
sm1(P," res*resol1 ");
ox_check_errors2(P);
return(pop(P));
}
def gcd_aux(A,B) {
if (type(A) == 1 && type(B) == 1) return(igcd(A,B));
else return(gcd(A,B));
}
def lcm_aux(V) { /* lcm_aux([3,5,6]); */
N = length(V);
if (N == 0) return(0);
if (N == 1) return(V[0]);
L = V[0];
for (I=1; I<N; I++) {
L = red(L*V[I]/gcd_aux(L,V[I]));
}
return(L);
}
def mul_v(V,S) {
if (type(V) == 4) {
return(map(mul_v,V,S));
} else {
return(V*S);
}
}
def div_v(V,S) {
if (type(V) == 4) {
return(map(div_v,V,S));
} else {
return(V/S);
}
}
def rat_to_p_aux(T) { /* cf. rat2plist2 */
if (T == 0) return([0,1]);
T = red(T);
T1 = nm(T); T1a = ptozp(T1);
T1b = red(T1a/T1);
T2 = dn(T);
return([T1a*dn(T1b),T2*nm(T1b)]);
}
def denom_aux0(A) {
return(A[1]);
}
def num_aux0(P) {
return(P[0]);
}
/*&usage begin: sm1.rat_to_p(F)
It returns the denominator of F and the numerator of F.
They are returned in a list.
All elements of the denominator and numerator are polynomial objects with integer coefficients. Note that dn and nm do not regard rational numbers as a factional object
and this function is necessary to send data to sm1, which accept only integers and
does not accept rational numbers.
example:
sm1.rat_to_p(1/2*x+1);
[x+2,2]
sm1.rat_to_p([1/2*x,1/3*x]);
[[3*x,2*x],6]
end: */
def rat_to_p(T) {
if (type(T) == 4) {
A = map(rat_to_p,T);
D = map(denom_aux0,A);
N = map(num_aux0,A);
L = lcm_aux(D);
B = newvect(length(N));
for (I=0; I<length(N); I++) {
B[I] = mul_v(N[I],L/D[I]);
}
return([vtol(B),L]);
}else{
return(rat_to_p_aux(T));
}
}
/* ---------------------------------------------- */
def distraction(A) {
SM1L_FIND_PROC(P);
push_int0(P,A);
sm1(P," distraction2* ");
ox_check_errors2(P);
return(pop(P));
}
/* Temporary functions */
/* Use this function for a while to wait a fix of asir. */
def sm1l_ntoint32(I) { /* Fixed */
SM1L_FIND_PROC(P);
if (I >= 0) return(ntoint32(I));
sm1(P," "+rtostr(I)+" ");
return(ox_pop_cmo(P));
}
def to_ascii_array(S) { /* Use strtoascii */
SM1L_FIND_PROC(P);
ox_push_cmo(P,S);
sm1(P," (array) dc { (universalNumber) dc } map ");
return(ox_pop_cmo(P));
}
def from_ascii_array(S) { /* Use asciitostr */
SM1L_FIND_PROC(P);
ox_push_cmo(P,S);
sm1(P," { (integer) dc (string) dc } map cat ");
return(ox_pop_cmo(P));
}
/*
[288] to_ascii_array("Hello");
[72,101,108,108,111]
[289] from_ascii_array(@@);
Hello
*/
/* end of temporary functions */
def gkz(S) {
SM1L_FIND_PROC(P);
A = S[0];
if (matrix_rank(A) != length(A)) error("Rank of A must be equal to the row numbers.\n");
B = S[1];
AA = [ ];
BB = [ ];
for (I=0; I<length(A); I++) {
AA = append(AA,[map(ntoint32,A[I])]);
BB = append(BB,[ntoint32(0)]);
}
sm1(P,"[ ");
push_int0(P,AA);
push_int0(P,BB);
sm1(P," ] gkz ");
ox_check_errors2(P);
R = pop(P);
RR0 = map(eval_str,R[0]);
RR1 = map(eval_str,R[1]);
RR3 = [ ];
for (I=0; I<length(B); I++) {
RR3 = append(RR3,[ rat_to_p(RR0[I]-B[I])[0] ]);
}
for (I=length(B); I<length(RR0); I++) {
RR3 = append(RR3,[RR0[I]]);
}
return([RR3,RR1]);
}
def mgkz(S) {
SM1L_FIND_PROC(P);
A = S[0];
W = S[1];
B = S[2];
AA = [ ];
WW = [ ];
BB = [ ];
for (I=0; I<length(A); I++) {
AA = append(AA,[map(ntoint32,A[I])]);
BB = append(BB,[ntoint32(0)]);
}
for (I=0; I<length(W); I++) {
WW = append(WW,[map(ntoint32,W[I])]);
}
sm1(P,"[ ");
push_int0(P,AA);
push_int0(P,WW);
push_int0(P,BB);
sm1(P," ] mgkz ");
ox_check_errors2(P);
R = pop(P);
RR0 = map(eval_str,R[0]);
RR1 = map(eval_str,R[1]);
RR3 = [ ];
for (I=0; I<length(B); I++) {
RR3 = append(RR3,[ rat_to_p(RR0[I]-B[I])[0] ]);
}
for (I=length(B); I<length(RR0); I++) {
RR3 = append(RR3,[RR0[I]]);
}
return([RR3,RR1]);
}
def appell1(S) {
N = length(S)-2;
B = cdr(cdr(S));
A = S[0];
C = S[1];
V = [ ];
for (I=0; I<N; I++) {
V = append(V,[sm1aux_x(I+1)]);
}
/* get vars of parameter vector */
VV=0;
for (I=0; I<length(S); I++) {
VV = VV+sm1_tttt^I*S[I];
}
VV = vars(VV); VVV = [];
for (I=0; I<length(VV); I++) {
if (VV[I] != sm1_tttt) {
VVV=cons(VV[I],VVV);
}
}
V0 = V;
V = append(V,reverse(VVV));
Ans = [ ];
Euler = 0;
for (I=0; I<N; I++) {
Euler = sm1aux_x(I+1)*sm1aux_dx(I+1) + Euler;
}
for (I=0; I<N; I++) {
T = mul(sm1aux_dx(I+1), Euler+C-1,V)-
mul(Euler+A, sm1aux_x(I+1)*sm1aux_dx(I+1)+B[I],V);
/* Tmp=rat_to_p(T);
print(Tmp[0]/Tmp[1]-T)$ */
T = rat_to_p(T)[0];
Ans = append(Ans,[T]);
}
for (I=0; I<N; I++) {
for (J=I+1; J<N; J++) {
T = (sm1aux_x(I+1)-sm1aux_x(J+1))*sm1aux_dx(I+1)*sm1aux_dx(J+1)
- B[J]*sm1aux_dx(I+1) + B[I]*sm1aux_dx(J+1);
/* Tmp=rat_to_p(T);
print(Tmp[0]/Tmp[1]-T)$ */
T = rat_to_p(T)[0];
Ans = append(Ans,[T]);
}
}
return([Ans,V0]);
}
def sm1aux_dx(I) {
return(strtov("dx"+rtostr(I)));
}
def sm1aux_x(I) {
return(strtov("x"+rtostr(I)));
}
def appell4(S) {
N = length(S)-2;
B = cdr(cdr(S));
A = S[0];
C = S[1];
V = [ ];
for (I=0; I<N; I++) {
V = append(V,[sm1aux_x(I+1)]);
}
/* get vars of parameter vector */
VV=0;
for (I=0; I<length(S); I++) {
VV = VV+sm1_tttt^I*S[I];
}
VV = vars(VV); VVV = [];
for (I=0; I<length(VV); I++) {
if (VV[I] != sm1_tttt) {
VVV=cons(VV[I],VVV);
}
}
V0 = V;
V = append(V,reverse(VVV));
Ans = [ ];
Euler = 0;
for (I=0; I<N; I++) {
Euler = sm1aux_x(I+1)*sm1aux_dx(I+1) + Euler;
}
for (I=0; I<N; I++) {
T = mul(sm1aux_dx(I+1), sm1aux_x(I+1)*sm1aux_dx(I+1)+B[I]-1,V)-
mul(Euler+A,Euler+C,V);
/* Tmp=rat_to_p(T);
print(Tmp[0]/Tmp[1]-T)$ */
T = rat_to_p(T)[0];
Ans = append(Ans,[T]);
}
return([Ans,V0]);
}
def rank(A) {
SM1L_FIND_PROC(P);
push_int0(P,A);
sm1(P," rank toString .. ");
ox_check_errors2(P);
R = pop(P);
return(R);
}
def rrank(A) {
SM1L_FIND_PROC(P);
push_int0(P,A);
sm1(P," rrank toString .. ");
ox_check_errors2(P);
R = pop(P);
return(R);
}
def auto_reduce(T) {
SM1L_FIND_PROC(P);
sm1(P,"[(AutoReduce) "+rtostr(T)+" ] system_variable ");
ox_check_errors2(P);
R = pop(P);
return(R);
}
def slope(II,V,FF,VF) {
SM1L_FIND_PROC(P);
A =[II,V,FF,VF];
push_int0(P,A);
sm1(P," slope toString ");
ox_check_errors2(P);
R = eval_str(pop(P));
return(R);
}
/*&usage begin: sm1.restriction(I,V,R|degree)
It computes the restriction of {I} as a D-module to the set defined by {R}.
{V} is the list of variables.
When the optional variable {degree}=d is given, only the restrictions
from 0 to d are computed.
Note that, in case of vector input, RESTRICTION VARIABLES MUST APPEAR FIRST
in the list of variable {V}. We are using wbfRoots to get the roots of
b-functions, so we can use only generic weight vector for now.
example: sm1.restriction([dx^2-x,dy^2-1],[x,y],[y]);
The output [[n0,F0],[n1,F1],...] means that H^0=D^n0/F0, H^(-1)=D^n1/F1, ...
The free basis of the vector space D^n is denoted by e0, e1, ...
algorithm: T.Oaku and N.Takayama, math.AG/9805006, http://xxx.lanl.gov
end: */
def restriction(II,V,R) {
SM1L_FIND_PROC(P);
Opt = getopt();
Deg = -1;
for (I=0; I < length(Opt); I++) {
if (Opt[I][0] == "degree") {
Deg = eval_str(rtostr(Opt[I][1]));
}else{
print(Opt[I][0],0);
error("sm1. : restriction: unknown option.");
}
}
if (Deg < 0) {
A =[II,R,[V,[ ]]];
}else{
A = [II,R,[V,[ ]],Deg];
}
push_int0(P,A);
sm1(P," restriction ");
ox_check_errors2(P);
R = pop(P);
return(R);
}
/*&usage begin: sm1.integration(I,V,R|degree)
It computes the integration of {I} as a D-module to the set defined by {R}.
{V} is the list of variables.
When the optional variable {degree}=d is given, only the integrations
from 0 to d are computed.
Note that, in case of vector input, INTEGRATION VARIABLES MUST APPEAR FIRST
in the list of variable {V}. We are using wbfRoots to get the roots of
b-functions, so we can use only generic weight vector for now.
example: sm1.integration([dt - (3*t^2-x), dx + t],[t,x],[t]);
The output [[n0,F0],[n1,F1],...] means that H^0=D^n0/F0, H^(-1)=D^n1/F1, ...
The free basis of the vector space D^n is denoted by e0, e1, ...
algorithm: T.Oaku and N.Takayama, math.AG/9805006, http://www.arxiv.org
end: */
def integration(II,V,R) {
SM1L_FIND_PROC(P);
Opt = getopt();
Deg = -1;
for (I=0; I < length(Opt); I++) {
if (Opt[I][0] == "degree") {
Deg = eval_str(rtostr(Opt[I][1]));
}else{
print(Opt[I][0],0);
error("sm1. : restriction: unknown option.");
}
}
if (Deg < 0) {
A =[II,R,[V,[ ]]];
}else{
A = [II,R,[V,[ ]],Deg];
}
push_int0(P,A);
sm1(P," integration ");
ox_check_errors2(P);
R = pop(P);
return(R);
}
/*&usage begin: sm1.ahg(A)
It idential with sm1.gkz({A}).
end: */
def ahg(A) {
return gkz(A);
}
/*&usage begin: sm1.bfunction(F)
It computes the global b-function of {F}.
description:
It no longer calls sm1's original bfunction. Instead, it calls
asir "bfct".
example: sm1.bfunction(x^2-y^3);
algorithm: M.Noro, Mathematical Software, icms 2002, pp.147--157.
end: */
def bfunction(F) {
return bfct(F);
}
/*&usage begin: sm1.generalized_bfunction(I,V,VD,W)
It computes the generalized b-function (indicial equation) of {I}
with respect to the weight {W}.
description:
It no longer calls sm1's original function. Instead, it calls
asir "generic_bfct".
example: sm1.generalized_bfunction([x^2*dx^2-1/2,dy^2],[x,y],[dx,dy],[-1,0,1,0]);
end: */
def generalized_bfunction(I,V,VD,W) {
return generic_bfct(I,V,VD,W);
}
/*&usage begin: sm1.saturation(T)
{T} = [{I},{J},{V}].
It returns saturation of {I} with respect to {J}^infty.
{V} is a list of variables.
example: sm1.saturation([[x2^2,x2*x4, x2, x4^2], [x2,x4], [x2,x4]]);
end: */
def saturation(A) {
SM1L_FIND_PROC(P);
push_int0(P,A);
sm1(P," saturation ");
ox_check_errors2(P);
R = pop(P);
return(R);
}
/*&usage begin: sm1.ecartd_gb(A)
It returns a standard basis of {A} by using
a tangent cone algorithm. h[0,1](D)-homogenization is used.
If the option rv="dp" (return_value="dp") is given,
the answer is returned in distributed polynomials.
Note. Functions in the category ecart changes the global environment
in the sm1 server. If you interrupted these functions,
run sm1.ecartd_gb with a small input and terminate it normally.
Then, the global environment is reset to the normal state.
Reference. G. Granger, T. Oaku, N. Takayama, Tangent cone algorithm for homogeized differential operators, 2005.
example:
input1
F=[2*(1-x-y)*dx+1,2*(1-x-y)*dy+1]$
FF=[F,"x,y",[[dx,1,dy,1],[x,-1,y,-1]]]$
sm1.ecartd_gb(FF);
output1
[[(-2*x-2*y+2)*dx+h,(-2*x-2*y+2)*dy+h],[(-2*x-2*y+2)*dx,(-2*x-2*y+2)*dy]]
input2
F=[2*(1-x-y)*dx+h,2*(1-x-y)*dy+h]$
FF=[F,"x,y",[[dx,1,dy,1],[x,-1,y,-1,dx,1,dy,1]],["noAutoHomogenize",1]]$
sm1.ecartd_gb(FF);
input3
F=[[x^2,y+x],[x+y,y^3], [2*x^2+x*y,y+x+x*y^3]]$
FF=[F,"x,y",[[dx,1,dy,1],[x, -1, y, -1,dx, 1, dy, 1]],
["degreeShift",[[0,1],[-3,1]]]]$
sm1.ecartd_gb(FF);
end: */
def ecartd_gb(A) {
SM1L_FIND_PROC(P);
Dp = is_in_dp(getopt(),A);
A=Dp[1]; Dp=Dp[0];
Sorted = getopt(sorted);
if (type(Sorted) == -1) {
Sorted = 0;
}else{
Sorted = 1;
}
Dehomo = getopt(dehomogenize);
if (type(Dehomo) == -1) {
Dehomo = 0;
}else{
Dehomo = 1;
}
P = check_server(P);
if ((length(A[0]) >= 1) && (type(A[0][0]) != type(""))) A = cons(rat_to_p(A[0])[0],cdr(A));
check_arg_gb(A,"Error in ecartd_gb");
push_int0(P,A);
sm1(P," ecartd.gb ");
if (Sorted) {
T = pop2(P);
if (Dehomo) {
return(append(T[0],[toOrdered(T[1]|dehomogenize=1)]));
}else{
return(append(T[0],[toOrdered(T[1])]));
}
}else{
Ans = pop(P);
if (Dp) {
Opt = ecartd_gb_oxRingStructure();
Ans = gb_ptod(Ans,Opt);
}
return Ans;
}
}
/*&usage begin: sm1.ecartd_reduction(F,A)
It returns a reduced form of {F} in terms of {A} by using
a tangent cone algorithm. h[0,1](D)-homogenization is used.
When the output is {G}, {G}[3] is {F} and
G[0]-(G[1]*A-sum(k,G[2][k]*G[3][k]))=0 holds.
{F} must be (0,1)-hohomogenized (see sm1.ecart_homogenize01Ideal).
This function does not check if the given order is admissible for the ecart
reduction. To do this check, use sm1.ecartd_gb.
example:
input
F=[2*(1-x-y)*dx+h,2*(1-x-y)*dy+h]$
FF=[F,"x,y",[[dx,1,dy,1],[x,-1,y,-1]]]$
G=sm1.ecartd_reduction(dx+dy,FF);
G[0]-(G[1]*(dx+dy)+G[2][0]*F[0]+G[2][1]*F[1]);
end: */
def ecartd_reduction(F,A) {
SM1L_FIND_PROC(P);
Dp = is_in_dp(getopt(),A);
A=Dp[1]; Dp=Dp[0];
F = gb_dtosm1p(F,A);
Sorted = getopt(sorted);
if (type(Sorted) == -1) {
Sorted = 0;
}else{
Sorted = 1;
}
Dehomo = getopt(dehomogenize);
if (type(Dehomo) == -1) {
Dehomo = 0;
}else{
Dehomo = 1;
}
P = check_server(P);
if ((length(A[0]) >= 1) && (type(A[0][0]) != type(""))) A = cons(rat_to_p(A[0])[0],cdr(A));
check_arg_gb(A,"Error in ecartd_reduction");
push_int0(P,F);
push_int0(P,A);
sm1(P," ecartd.reduction ");
if (Sorted) {
T = pop2(P);
if (Dehomo) {
return(append(T[0],[toOrdered(T[1]|dehomogenize=1)]));
}else{
return(append(T[0],[toOrdered(T[1])]));
}
}else{
Ans = pop(P);
if (Dp) {
Opt = ecartd_gb_oxRingStructure();
Ans = gb_ptod(Ans,Opt);
}
return Ans;
}
}
/*&usage begin: sm1.gb_reduction(F,A)
It returns a reduced form of {F} in terms of {A} by using
a normal form algorithm. h[1,1](D)-homogenization is used.
example:
input
F=[2*(h-x-y)*dx+h^2,2*(h-x-y)*dy+h^2]$
FF=[F,"x,y",[[dx,1,dy,1],[x,-1,y,-1,dx,1,dy,1]]]$
sm1.gb_reduction((h-x-y)^2*dx*dy,FF);
end: */
def gb_reduction(F,A) {
SM1L_FIND_PROC(P);
Dp = is_in_dp(getopt(),A);
A=Dp[1]; Dp=Dp[0];
F = gb_dtosm1p(F,A);
Sorted = getopt(sorted);
if (type(Sorted) == -1) {
Sorted = 0;
}else{
Sorted = 1;
}
Dehomo = getopt(dehomogenize);
if (type(Dehomo) == -1) {
Dehomo = 0;
}else{
Dehomo = 1;
}
P = check_server(P);
if ((length(A[0]) >= 1) && (type(A[0][0]) != type(""))) A = cons(rat_to_p(A[0])[0],cdr(A));
check_arg_gb(A,"Error in gb_reduction");
push_int0(P,F);
push_int0(P,A);
sm1(P," gb.reduction ");
if (Sorted) {
T = pop2(P);
if (Dehomo) {
return(append(T[0],[toOrdered(T[1]|dehomogenize=1)]));
}else{
return(append(T[0],[toOrdered(T[1])]));
}
}else{
Ans = pop(P);
if (Dp) {
Opt = gb_oxRingStructure();
Ans = gb_ptod(Ans,Opt);
}
return Ans;
}
}
/*&usage begin: sm1.gb_reduction_noh(F,A)
It returns a reduced form of {F} in terms of {A} by using
a normal form algorithm.
example:
input
F=[2*dx+1,2*dy+1]$
FF=[F,"x,y",[[dx,1,dy,1]]]$
sm1.gb_reduction_noh((1-x-y)^2*dx*dy,FF);
end: */
def gb_reduction_noh(F,A) {
SM1L_FIND_PROC(P);
Dp = is_in_dp(getopt(),A);
A=Dp[1]; Dp=Dp[0];
F = gb_dtosm1p(F,A);
P = check_server(P);
if ((length(A[0]) >= 1) && (type(A[0][0]) != type(""))) A = cons(rat_to_p(A[0])[0],cdr(A));
check_arg_gb(A,"Error in gb_reduction_noh");
push_int0(P,F);
push_int0(P,A);
sm1(P," gb.reduction_noh ");
Ans = pop(P);
if (Dp) {
Opt = gb_oxRingStructure();
Ans = gb_ptod(Ans,Opt);
Ans = gb_dehomo(Ans);
}
return Ans;
}
/*&usage begin: sm1.gb_oxRingStructure()
It returns the oxRingStructure of the most recent gb computation.
end: */
def gb_oxRingStructure() {
SM1L_FIND_PROC(P);
P = check_server(P);
sm1(P," gb.oxRingStructure ");
Option_list = pop(P);
return aux_parse_gb_option_list(Option_list[0]);
}
/*&usage begin: sm1.ecartd_isSameIdeal_h(F)
Here, {F}=[{II},{JJ},{V}].
It compares two ideals {II} and {JJ} in h[0,1](D)_alg.
example:
input
II=[(1-x)^2*dx+h*(1-x)]$ JJ = [(1-x)*dx+h]$
V=[x]$
sm1.ecartd_isSameIdeal_h([II,JJ,V]);
end: */
def ecartd_isSameIdeal_h(F) {
SM1L_FIND_PROC(P);
Sorted = getopt(sorted);
if (type(Sorted) == -1) {
Sorted = 0;
}else{
Sorted = 1;
}
Dehomo = getopt(dehomogenize);
if (type(Dehomo) == -1) {
Dehomo = 0;
}else{
Dehomo = 1;
}
P = check_server(P);
II = F[0]; JJ = F[1]; V=F[2];
if (type(V) != 7) V=map(rtostr,V);
push_int0(P,[II,JJ,V]);
sm1(P," ecartd.isSameIdeal_h ");
if (Sorted) {
T = pop2(P);
if (Dehomo) {
return(append(T[0],[toOrdered(T[1]|dehomogenize=1)]));
}else{
return(append(T[0],[toOrdered(T[1])]));
}
}else{
return(pop(P));
}
}
/*&usage begin: sm1.ecart_homogenize01Ideal(A)
It (0,1)-homogenizes the ideal {A}[0]. Note that it is not an elementwise
homogenization.
example:
input1
F=[(1-x)*dx+1]$ FF=[F,"x,y"]$
sm1.ecart_homogenize01Ideal(FF);
intput2
F=sm1.appell1([1,2,3,4]);
sm1.ecart_homogenize01Ideal(F);
end: */
def ecart_homogenize01Ideal(A) {
SM1L_FIND_PROC(P);
P = check_server(P);
check_arg_gb(A,"Error in ecart_homogenize01Ideal");
push_int0(P,A);
sm1(P," ecart.homogenize01Ideal ");
return(pop(P));
}
def aux_parse_gb_option_list(Ol) {
Ans = [];
while (length(Ol) != 0) {
L = car(Ol); Ol = cdr(Ol);
Ans = cons([L[0],map(eval_str,map(rtostr,L[1]))], Ans);
}
return reverse(Ans);
}
/*&usage begin: sm1.ecartd_gb_oxRingStructure()
It returns the oxRingStructure of the most recent ecartd_gb computation.
end: */
def ecartd_gb_oxRingStructure() {
SM1L_FIND_PROC(P);
P = check_server(P);
sm1(P," ecartd.gb.oxRingStructure ");
Option_list = pop(P);
return aux_parse_gb_option_list(Option_list[0]);
}
/*&usage begin: sm1.isSameIdeal_in_Dalg(I,J,V)
It compares two ideals {I} and {J} in D_alg (algebraic D
with variables {V}, no homogenization).
example:
Input1
II=[(1-x)^2*dx+(1-x)]$ JJ = [(1-x)*dx+1]$ V=[x]$
sm1.isSameIdeal_in_Dalg(II,JJ,V);
end: */
def isSameIdeal_in_Dalg(I,J,V) {
SM1L_FIND_PROC(P);
P = check_server(P);
Ih = ecart_homogenize01Ideal([I,V]);
Jh = ecart_homogenize01Ideal([J,V]);
/* print([Ih,Jh,V]); */
R = ecartd_isSameIdeal_h([Ih,Jh,V]);
return R;
}
/*&usage begin: sm1.ecartd_reduction_noh(F,A)
It returns a reduced form of {F} in terms of {A} by using a
tangent cone algorithm. h[0,1](D)-homogenization is NOT used.
{A}[0] must not contain the variable h.
example:
F=[2*(1-x-y)*dx+1,2*(1-x-y)*dy+1]$
FF=[F,"x,y",[[dx,1,dy,1],[x,-1,y,-1]]]$
sm1.ecartd_reduction_noh(dx+dy,FF);
end: */
def ecartd_reduction_noh(F,A) {
SM1L_FIND_PROC(P);
Dp = is_in_dp(getopt(),A);
A=Dp[1]; Dp=Dp[0];
F = gb_dtosm1p(F,A);
P = check_server(P);
if ((length(A[0]) >= 1) && (type(A[0][0]) != type(""))) A = cons(rat_to_p(A[0])[0],cdr(A));
check_arg_gb(A,"Error in ecartd_reduction_noh");
push_int0(P,F);
push_int0(P,A);
sm1(P," ecartd.reduction_noh ");
Ans = pop(P);
if (Dp) {
Opt = ecartd_gb_oxRingStructure();
Ans = gb_ptod(Ans,Opt);
Ans = gb_dehomo(Ans);
}
return Ans;
}
def sm1vstr_to_vlist(V) {
S = strtoascii(V);
/* , = 44 */
Ans = [ ];
while (S != []) {
P = [];
while (length(S) != 0 && S[0] != 44) {
P = cons(car(S),P); S = cdr(S);
}
if (P != []) {
P = reverse(P);
Ans = cons(asciitostr(P),Ans);
}
S = cdr(S);
}
return reverse(Ans);
}
def sm1vlist_to_asirvlist(V) {
if (type(V) == 7) {
V = sm1vstr_to_vlist(V);
}else{
V = map(rtostr,V);
}
DV = [];
for (I=0; I<length(V); I++) {
DV = cons("d"+V[I],DV);
}
DV = reverse(DV);
V = append(V,DV);
V = map(strtov,V);
return V;
}
def gb_ptod2(F,Opt) {
if (type(F) == 4) return gb_ptod(F,Opt);
else {
G = dp_ptod( F | option_list=Opt);
return G;
}
}
def gb_ptod(F,Opt) {
if (type(F) == 4) map(gb_ptod2,F,Opt);
else gb_ptod2(F,Opt);
}
def gb_dtop(F,V) {
if (type(F) == 4) {
return map(dp_dtop,F,V);
}else return dp_dtop(F,V);
}
/* A is in the format of gb. */
def is_in_dp(Opt,A) {
Dp = 0;
while (Opt != []) {
if (Opt[0][0] == "return_value") {
if (Opt[0][1] == "dp") {
Dp = 1;
}
}else if (Opt[0][0] == "rv") {
if (Opt[0][1] == "dp") {
Dp = 1;
}
}
Opt = cdr(Opt);
}
if (type(A[0][0]) == 9) { /* argument is dp */
Dp = 1;
VV = sm1vlist_to_asirvlist(A[1]);
N = length(dp_etov(A[0][0]));
if (N % 2 != 0) VV = append(VV,[h]);
AA = gb_dtop(A[0],VV);
return [Dp,cons(AA,cdr(A))];
}
return [Dp,A];
}
def gb_dtosm1p(F,A) {
if (type(F) == 9) {
N = length(dp_etov(F));
V = sm1vlist_to_asirvlist(A[1]);
if (N % 2 == 1) V = append(V,[h]);
return dp_dtop(F,V);
}else return F;
}
def gb_dehomo2(F) {
if (type(F) == 4) {
return gb_dehomo(F);
}else{
return dp_dehomo(F);
}
}
def gb_dehomo(F) {
if (type(F) == 4) {
return map(gb_dehomo2,F);
}else{
return dp_dehomo(F);
}
}
/*&testdata for rv="dp"
F=[2*(1-x-y)*dx+h,2*(1-x-y)*dy+h]$
FF=[F,"x,y",[[dx,1,dy,1],[x,-1,y,-1]]]$
sm1.ecartd_reduction(dx+dy,FF | rv="dp");
F=[(2)*<<0,0,1,0,0>>+(-2)*<<1,0,1,0,0>>+(-2)*<<0,1,1,0,0>>+(1)*<<0,0,0,0,1>>,(2)*<<0,0,0,1,0>>+(-2)*<<1,0,0,1,0>>+(-2)*<<0,1,0,1,0>>+(1)*<<0,0,0,0,1>>];
FF=[F,"x,y",[[dx,1,dy,1],[x,-1,y,-1]]]$
sm1.ecartd_reduction(<<0,0,1,0,0>>,FF);
F=[(2)*<<0,0,1,0>>+(-2)*<<1,0,1,0>>+(-2)*<<0,1,1,0>>+(1)*<<0,0,0,0>>,(2)*<<0,0,0,1>>+(-2)*<<1,0,0,1>>+(-2)*<<0,1,0,1>>+(1)*<<0,0,0,0>>];
FF=[F,"x,y",[[dx,1,dy,1],[x,-1,y,-1]]]$
sm1.ecartd_reduction_noh(<<0,0,1,0>>,FF);
sm1.gb([sm1.appell1([1/2,1/3,1/5,-1/7])[0],[x1,x2],[[dx1,1,dx2,1]]] | rv="dp");
*/
/*&usage begin: sm1.ecartd_syz(A)
It returns a syzygy of {A} by using
a tangent cone algorithm. h[0,1](D)-homogenization is used.
If the option rv="dp" (return_value="dp") is given,
the answer is returned in distributed polynomials.
The return value is in the format [s,[g,m,t]].
s is the generator of the syzygies, g is the Grobner basis,
m is the translation matrix from the generators to g.
t is the syzygy of g.
example:
input1
F=[2*(1-x-y)*dx+1,2*(1-x-y)*dy+1]$
FF=[F,"x,y",[[dx,1,dy,1],[x,-1,y,-1]]]$
sm1.ecartd_syz(FF);
input2
F=[2*(1-x-y)*dx+h,2*(1-x-y)*dy+h]$
FF=[F,"x,y",[[dx,1,dy,1],[x,-1,y,-1,dx,1,dy,1]],["noAutoHomogenize",1]]$
sm1.ecartd_syz(FF);
end: */
def ecartd_syz(A) {
SM1L_FIND_PROC(P);
Dp = is_in_dp(getopt(),A);
A=Dp[1]; Dp=Dp[0];
P = check_server(P);
if ((length(A[0]) >= 1) && (type(A[0][0]) != type(""))) A = cons(rat_to_p(A[0])[0],cdr(A));
check_arg_gb(A,"Error in ecartd_syz");
push_int0(P,A);
sm1(P," ecart.syz ");
Ans = pop(P);
if (Dp) {
Opt = ecartd_gb_oxRingStructure();
Ans = gb_ptod(Ans,Opt);
}
return Ans;
}
/*&usage begin: sm1.laplace(L,V,VL)
It returns the Laplace transformation of L for {VL}.
{V} is the list of space variables.
The numbers in coefficients must not be rational with a non-1 denominator.
cf. ptozp
example:
L1=sm1.laplace(dt-(3*t^2-x),[x,t],[t,dt]);
L2=sm1.laplace(dx+t,[x,t],[t,dt]);
sm1.restriction([L1,L2],[t,x],[t] | degree=0);
end: */
def laplace(L,V,VL) {
SM1L_FIND_PROC(P);
ringD(V,[[]]);
push_int0(P,L);
sm1(P," . ");
push_int0(P,VL);
sm1(P," laplace0 ");
Ans = pop(P);
return Ans;
}
/*&usage begin sm1.add_d(V)
It adds d to V.
example: sm1.add_d([x1,x2,x3]);
end: */
def add_d(V) {
if (type(V) >= 4) return map(add_d,V);
return(eval_str("d"+rtostr(V)));
}
/*&usage begin sm1.regular_triangulation(AW)
It returns the regular triangulation for the point configuration {AW[0]}
with respect to the weight {AW[1]}.
Ans[0] is the list of simpleces. Ans[1] is the initial ideal. If it is
not a monomial ideal, it is not a regular triangulation.
cf. sm1.gkz
example: sm1.regular_triangulation([ [[1,1,1,1],[0,0,1,1],[1,0,1,0]], [1,0,0,1]]);
sm1.regular_triangulation([ [[1,1,1,1,1,1],[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0]], [2,0,0,1,0,5]]);
sm1.regular_triangulation([ [[1,1,1,1,1],[3,2,2,2,1],[0,1,-1,0,0]], [1,1,1,0,1]]);
end: */
def regular_triangulation(AW) {
A=AW[0];
W=AW[1];
Beta=vtol(newvect(length(A)));
F=sm1.gkz([A,Beta]);
T=F[0]; V=F[1];
for (I=0; I<length(A); I++) T=cdr(T);
DV=add_d(V);
Wvec=newvect(2*length(DV));
for (I=0; I<length(DV); I++) {
Wvec[2*I]=DV[I]; Wvec[2*I+1]=W[I];
}
Wvec=vtol(Wvec);
print([T,[Wvec]]);
G=gb([T,DV,[Wvec]]);
P=primadec(G[1],DV);
Ans=[];
for (I=0; I<length(P); I++) {
Q=P[I][1];
Ans=cons(base_set_minus(DV,Q),Ans);
}
return [Ans,G[1],G[0],DV];
}
def mod_reduction(A) {
SM1L_FIND_PROC(P);
A=map(matrix_matrix_to_list,A);
/* check the arguments */
if (type(A) != 4) {
error("sm1.reduction_mod(A|ringvar=p): A must be a list.");
}
// auto_reduce(1);
// sm1(P," [(parse) [(getenv) (HOME)] extension (/this19/misc-2019/09/hgs/sred.sm1) 2 cat_n pushfile] extension ");
if (type(getopt(ring_var)) > 0) {
Ring_var = rtostr(getopt(ring_var));
}else{
Ring_var = "ring_var_for_asir";
}
sm1(P," [(CurrentRingp)] pushEnv /env_for_asir set ");
sm1(P,sprintf(" %a ring_def ",Ring_var));
push_int0(P,A);
sm1(P," mod_reduction* ");
ox_check_errors2(P);
Ans = pop(P); /* Throw away Ans[3] */
sm1(P," env_for_asir popEnv ");
return([Ans[0],Ans[1],Ans[2],A[1]]);
}
def reduction_verbose(A)
"reduction_verbose([F,G,V,Weight]) returns [R,R0,S,G,LT,Ord] where R*R0+S*G=0 and LT is the initial of R with the Order"
{
/* Example: reduction(A|proc=10) */
SM1L_FIND_PROC(P);
/* check the arguments */
if (type(A) != 4) {
error("sm1. reduction(A|proc=p): A must be a list.");
}
if (length(A) == 2) return mod_reduction(A | option_list=getopt());
if (type(A[0]) != 4) {
AA = [rtostr(A[0])];
AA = append(AA,[ map(rtostr,A[1]) ]);
}else{ // module case
AA = [map(rtostr,A[0])];
B=[];
for (I=0; I<length(A[1]); I++) {
B=cons(map(rtostr,A[1][I]),B);
}
AA = append(AA,[reverse(B)]);
}
AA = append(AA, cdr(cdr(A)));
/* printf("AA=%a\n",AA);
print(map(type,AA[0]));
print(map(type,AA[2]));
print(map(type,AA[3][0])); */
sm1(P," /reduction*.noH 0 def ");
push_int0(P,AA);
sm1(P," reduction* /tmp-reduction set tmp-reduction getRing :: tmp-reduction ");
ox_check_errors2(P);
Ans = pop(P); /* Throw away Ans[3] */
sm1(P," tmp-reduction 0 get init /tmp-reduction-init set tmp-reduction-init ");
Ans2 = pop(P);
sm1(P," tmp-reduction getRing (oxRingStructure) data_conversion ");
Ans3 = pop(P);
return([Ans[0],Ans[1],Ans[2],A[1],Ans2,Ans3]);
}
endmodule;
end$
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to sm1.rr CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / asir-contrib / packages / src