[SCM] an open source computer algebra system branch, cleanedupstream, updated. 6125e540ca6d66c307958938a9d53b245507c323
Bernhard R. Link
brlink at debian.org
Tue Apr 24 15:54:04 UTC 2012
The following commit has been merged in the cleanedupstream branch:
commit a9f6dff7b653daf660bc8ef16d3b30f2c323e50b
Author: Hans Schoenemann <hannes at mathematik.uni-kl.de>
Date: Thu Mar 8 19:28:25 2012 +0100
add: ffsolve.lib; contribution by Gergo Gyula Borus
diff --git a/Singular/LIB/ffsolve.lib b/Singular/LIB/ffsolve.lib
new file mode 100644
index 0000000..b454724
--- /dev/null
+++ b/Singular/LIB/ffsolve.lib
@@ -0,0 +1,752 @@
+////////////////////////////////////////////////////////////////////
+version="$Id$";
+category="Symbolic-numerical solving";
+info="
+LIBRARY: ffsolve.lib multivariate equation solving over finite fields
+AUTHOR: Gergo Gyula Borus, borisz at borisz.net
+KEYWORDS: multivariate equations; finite field
+
+PROCEDURES:
+ffsolve(); finite field solving using heuristically chosen method
+speoff(); solve system of multivariate equations over finite field
+simpleSolver(); solver using modified brute-force search
+gbsolve(); multivariate solver using Groebner-basis
+";
+
+LIB "presolve.lib";
+LIB "general.lib";
+LIB "ring.lib";
+LIB "standard.lib";
+LIB "matrix.lib";
+
+////////////////////////////////////////////////////////////////////
+proc ffsolve(ideal equations, list #)
+"USAGE: ffsolve(I[, L]); I ideal, L list of strings
+RETURN: list L, the common roots of I as ideal
+ASSUME: basering is a finite field of type (p^n,a)
+"
+{
+ list solutions, possibleSolvers, tempsols;
+ int i,j, k, found;
+ ideal factors, linfacs;
+ poly lp;
+ // check assumptions
+ if(npars(basering)>1){
+ ERROR("Basering must have at most one parameter");
+ }
+ if(char(basering)==0){
+ ERROR("Basering must have finite characteristic");
+ }
+
+ // heuristic choice of solver
+ // not yet implemented
+ if(size(#)){
+ if(size(#)==1 and typeof(#[1])=="list"){
+ possibleSolvers = #[1];
+ }else{
+ possibleSolvers = #;
+ }
+ }else{
+ possibleSolvers = "simpleSolver", "speoff", "gbsolve";
+ }
+ string solver = possibleSolvers[random(1,size(possibleSolvers))];
+
+ if(nvars(basering)==1){
+ return(facstd(equations));
+ }
+
+ // search for the first univariate polynomial
+ found = 0;
+ for(i=1; i<=ncols(equations); i++){
+ if(univariate(equations[i])>0){
+ factors=factorize(equations[i],1);
+ for(j=1; j<=ncols(factors); j++){
+ if(deg(factors[j])==1){
+ linfacs[size(linfacs)+1] = factors[j];
+ }
+ }
+ if(deg(linfacs[1])>0){
+ found=1;
+ break;
+ }
+ }
+ }
+ // if there is, collect its the linear factors
+ if(found){
+ // substitute the root and call recursively
+ ideal neweqs, invmapideal, ti;
+ map invmap;
+ for(k=1; k<=ncols(linfacs); k++){
+ lp = linfacs[k];
+ neweqs = reduce(equations, lp);
+
+ intvec varexp = leadexp(lp);
+ def original_ring = basering;
+ def newRing = clonering(nvars(original_ring)-1);
+ setring newRing;
+ ideal mappingIdeal;
+ j=1;
+ for(i=1; i<=size(varexp); i++){
+ if(varexp[i]){
+ mappingIdeal[i] = 0;
+ }else{
+ mappingIdeal[i] = var(j);
+ j++;
+ }
+ }
+ map recmap = original_ring, mappingIdeal;
+ list tsols = ffsolve(recmap(neweqs), possibleSolvers);
+ if(size(tsols)==0){
+ tsols = list(ideal(1));
+ }
+ setring original_ring;
+ j=1;
+ for(i=1;i<=size(varexp);i++){
+ if(varexp[i]==0){
+ invmapideal[j] = var(i);
+ j++;
+ }
+ }
+ invmap = newRing, invmapideal;
+ tempsols = invmap(tsols);
+
+ // combine the solutions
+ for(j=1; j<=size(tempsols); j++){
+ ti = std(tempsols[j]+lp);
+ if(deg(ti)>0){
+ solutions = insert(solutions,ti);
+ }
+ }
+ }
+ }else{
+ execute("solutions="+solver+"(equations);") ;
+ }
+ return(solutions);
+}
+example
+{
+ "EXAMPLE:";echo=2;
+ ring R = (2,a),x(1..3),lp;
+ minpoly=a2+a+1;
+ ideal I;
+ I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
+ I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
+ I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
+ ffsolve(I);
+}
+////////////////////////////////////////////////////////////////////
+proc speoff (ideal L, list #)
+"USAGE: speoff(I[, i]); I ideal, i optional integer
+RETURN: list if optional parameter is not given or set to 2,
+integer if optional is set to 0
+ASSUME: basering is a finite field of type (p^n,a)
+NOTE: When the optional parameter is set to 0, speoff only
+checks if I has common roots, then return 1, otherwise return 0.
+
+"
+{
+ system("--ticks-per-sec", 1000);
+ int t = rtimer;
+ int mode, i,j;
+ list results, rs, start;
+ poly g;
+ // check assumptions
+ if(npars(basering)>1){
+ ERROR("Basering must have at most one parameter");
+ }
+ if(char(basering)==0){
+ ERROR("Basering must have finite characteristic");
+ }
+
+ if( size(#) > 0 ){
+ mode = #[1];
+ }else{
+ mode = 2;
+ }
+ L = simplify(L,15);
+ g = productOfEqs( L );
+
+ if(g == 0){
+ if(voice==2){printf("Runtime: %s ms", rtimer-t)};
+ if(mode==0){
+ return(0);
+ }
+ return( list() );
+ }
+ if(g == 1){
+ list vectors = every_vector();
+ for(j=1; j<=size(vectors); j++){
+ ideal res;
+ for(i=1; i<=nvars(basering); i++){
+ res[i] = var(i)-vectors[j][i];
+ }
+ results[size(results)+1] = std(res);
+ }
+ if(voice==2){printf("Runtime: %s ms", rtimer-t)};
+ return( results );
+ }
+
+ if( mode == 0 ){
+ if(voice==2){printf("Runtime: %s ms", rtimer-t)};
+ return( 1 );
+ }else{
+ for(i=1; i<=nvars(basering); i++){
+ start[i] = 0:order_of_extension();
+ }
+
+ if( mode == 1){
+ results[size(results)+1] = melyseg(g, start);
+ }else{
+ while(1){
+ start = melyseg(g, start);
+ if( size(start) > 0 ){
+ ideal res;
+ for(i=1; i<=nvars(basering); i++){
+ res[i] = var(i)-vec2elm(start[i]);
+ }
+ results[size(results)+1] = std(res);
+ start = increment(start);
+ }else{
+ break;
+ }
+ }
+ }
+ }
+
+ if(voice==2){printf("Runtime: %s ms", rtimer-t)};
+ return(results);
+}
+example
+{
+ "EXAMPLE:";echo=2;
+ ring R = (2,a),x(1..3),lp;
+ minpoly=a2+a+1;
+ ideal I;
+ I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
+ I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
+ I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
+ ffsolve(I);
+}
+////////////////////////////////////////////////////////////////////
+proc simpleSolver(ideal E)
+"USAGE: simpleSolver(I); I ideal
+RETURN: list L, the common roots of I as ideal
+ASSUME: basering is a finite field of type (p^n,a)
+"
+{
+ int SStime = rtimer;
+ int i,j,k,t, correct;
+ list solutions = list(std(ideal()));
+ list partial_solutions;
+ ideal partial_system, curr_sol, curr_sys, factors;
+ poly univar_poly;
+ E = E+defaultIdeal();
+ // check assumptions
+ if(npars(basering)>1){
+ ERROR("Basering must have at most one parameter");
+ }
+ if(char(basering)==0){
+ ERROR("Basering must have finite characteristic");
+ }
+
+ for(k=1; k<=nvars(basering); k++){
+ partial_solutions = list();
+ for(i=1; i<=size(solutions); i++){
+ partial_system = reduce(E, solutions[i]);
+ for(j=1; j<=ncols(partial_system); j++){
+ if(univariate(partial_system[j])>0){
+ univar_poly = partial_system[j];
+ break;
+ }
+ }
+ factors = factorize(univar_poly,1);
+ for(j=1; j<=ncols(factors); j++){
+ if(deg(factors[j])==1){
+ curr_sol = std(solutions[i]+ideal(factors[j]));
+ curr_sys = reduce(E, curr_sol);
+ correct = 1;
+ for(t=1; t<=ncols(curr_sys); t++){
+ if(deg(curr_sys[t])==0){
+ correct = 0;
+ break;
+ }
+ }
+ if(correct){
+ partial_solutions = insert(partial_solutions, curr_sol);
+ }
+ }
+ }
+ }
+ solutions = partial_solutions;
+ }
+ if(voice==2){printf("Runtime: %s ms", rtimer-SStime)};
+ return(solutions);
+}
+example
+{
+ "EXAMPLE:";echo=2;
+ ring R = (2,a),x(1..3),lp;
+ minpoly=a2+a+1;
+ ideal I;
+ I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
+ I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
+ I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
+ ffsolve(I);
+}
+////////////////////////////////////////////////////////////////////
+proc gbsolve(ideal equation_system)
+"USAGE: gbsolve(I); I ideal
+RETURN: list L, the common roots of I as ideal
+ASSUME: basering is a finite field of type (p^n,a)
+"
+{
+ option(redSB);
+ system("--ticks-per-sec",1000);
+ int i,j, prop, newelement, number_new_vars;
+ int t=rtimer;
+ ideal ls;
+ list results, slvbl, linsol, ctrl, new_sols, varinfo;
+ ideal I, linear_solution, unsolved_part, univar_part, multivar_part, unsolved_vars;
+ intvec unsolved_var_nums;
+ string new_vars;
+ // check assumptions
+ if(npars(basering)>1){
+ ERROR("Basering must have at most one parameter");
+ }
+ if(char(basering)==0){
+ ERROR("Basering must have finite characteristic");
+ }
+
+ def original_ring = basering;
+ if(npars(basering)==1){
+ int prime_coeff_field=0;
+ string minpolystr = "minpoly="+
+ get_minpoly_str(size(original_ring),parstr(original_ring,1))+";" ;
+ }else{
+ int prime_coeff_field=1;
+ }
+
+ equation_system = simplify(equation_system,15);
+
+ ideal standard_basis = std(equation_system);
+ printf("std(): %s ms", rtimer-t);
+ list basis_factors = facstd(standard_basis);
+ printf("facstd(): %s ms", rtimer-t);
+ if( basis_factors[1][1] == 1){
+ printf("Runtime: %s ms", rtimer-t);
+ return(results)
+ };
+
+ for(i=1; i<= size(basis_factors); i++){
+ prop = 0;
+ for(j=1; j<=size(basis_factors[i]); j++){
+ if( univariate(basis_factors[i][j])>0 and deg(basis_factors[i][j])>1){
+ prop =1;
+ break;
+ }
+ }
+
+ if(prop == 0){
+ ls = Presolve::solvelinearpart( basis_factors[i] );
+ if(ncols(ls) == nvars(basering) ){
+ ctrl, newelement = add_if_new(ctrl, ls);
+ if(newelement){
+ results = insert(results, ls);
+ }
+ }else{
+ slvbl = insert(slvbl, list(basis_factors[i],ls) );
+ }
+ }
+ }
+
+ if(size(slvbl)<>0){
+ for(int E = 1; E<= size(slvbl); E++){
+ I = slvbl[E][1];
+ linear_solution = slvbl[E][2];
+ unsolved_part = reduce(I,linear_solution);
+ univar_part = ideal();
+ multivar_part = ideal();
+ for(i=1; i<=ncols(I); i++){
+ if(univariate(I[i])>0){
+ univar_part = univar_part+I[i];
+ }else{
+ multivar_part = multivar_part+I[i];
+ }
+ }
+
+ // list varinfo = Presolve::findvars(linear_solution,1);
+ varinfo = Presolve::findvars(univar_part,1);
+ unsolved_vars = varinfo[3];
+ unsolved_var_nums = varinfo[4];
+ number_new_vars = ncols(unsolved_vars);
+
+ new_vars = "@y(1.."+string(number_new_vars)+")";
+ def R_new = Ring::changevar(new_vars, original_ring);
+ setring R_new;
+ if( !prime_coeff_field ){
+ execute(minpolystr);
+ }
+
+ ideal mapping_ideal;
+ for(i=1; i<=size(unsolved_var_nums); i++){
+ mapping_ideal[unsolved_var_nums[i]] = var(i);
+ }
+
+ map F = original_ring, mapping_ideal;
+ ideal I_new = F( multivar_part );
+
+ list sol_new;
+ int unsolvable = 0;
+ sol_new = simpleSolver(I_new);
+ if( size(sol_new) == 0){
+ unsolvable = 1;
+ }
+
+ setring original_ring;
+
+ if(unsolvable){
+ list sol_old = list();
+ }else{
+ map G = R_new, unsolved_vars;
+ new_sols = G(sol_new);
+ for(i=1; i<=size(new_sols); i++){
+ ideal sol = new_sols[i]+linear_solution;
+ sol = std(sol);
+ ctrl, newelement = add_if_new(ctrl, sol);
+ if(newelement){
+ results = insert(results, sol);
+ }
+ kill sol;
+ }
+ }
+ kill G;
+ kill R_new;
+ }
+ }
+
+ if(voice==2){printf("Runtime: %s ms", rtimer-t)};
+ return( results );
+}
+example
+{
+ "EXAMPLE:";echo=2;
+ ring R = (2,a),x(1..3),lp;
+ minpoly=a2+a+1;
+ ideal I;
+ I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
+ I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
+ I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
+ ffsolve(I);
+}
+
+
+////////////////////////////////////////////////////////////////////
+
+
+static proc melyseg(poly g, list start)
+{
+ list gsub = g;
+ int i = 1;
+
+ while( start[1][1] <> char(basering) ){
+ gsub[i+1] = subst( gsub[i], var(i), vec2elm(start[i]));
+ if( gsub[i+1] == 0 ){
+ list new = increment(start,i);
+ for(int l=1; l<=size(start); l++){
+ if(start[l]<>new[l]){
+ i = l;
+ break;
+ }
+ }
+ start = new;
+ }else{
+ if(i == nvars(basering)){
+ return(start);
+ }else{
+ i++;
+ }
+ }
+ }
+ return(list());
+}
+
+static proc productOfEqs(ideal I)
+{
+ system("--no-warn", 1);
+ ideal eqs = sort_ideal(I);
+ int i,q;
+ poly g = 1;
+ q = size(basering);
+ ideal I = defaultIdeal();
+
+ for(i=1; i<=size(eqs); i++){
+ if(g==0){return(g);}
+ g = reduce(g*(eqs[i]^(q-1)-1), I);
+ }
+ return( g );
+}
+
+static proc notSoQuickReduce(poly P)
+{
+ int q = ringlist( basering )[1];
+ int nov = nvars( basering );
+ poly reducedPoly = 0;
+ poly M;
+ int i, j, e;
+ intvec exps;
+
+ for(i=1; i<=size(P); i++){
+ if( deg(P[i]) > q-1 ){
+ M = leadcoef( P[i] );
+ exps = leadexp( P[i] );
+ for(j=1; j<=size(exps); j++){
+ if( exps[j] > q-1 ){
+ e = exps[j] % (q-1);
+ if( e == 0 ){ e = q-1; }
+ }else{
+ e = exps[j];
+ }
+ M = M*var(j)^e;
+ }
+ }else{
+ M = P[i];
+ }
+ reducedPoly = reducedPoly + M;
+ }
+ return( reducedPoly );
+}
+
+static proc polypowmod(poly P, int E, ideal I)
+{
+ system("--no-warn", 1);
+ list pot;
+ poly rs = 1;
+
+ while( E > 0 ){
+ pot[ size(pot)+1 ] = E % 2;
+ E = E / 2;
+ if( pot[size(pot)] ){
+ rs = reduce(rs*P, I);
+ }
+ P = reduce(P*P,I);
+ }
+ return( reduce(rs, I) );
+}
+
+static proc clonering(list #)
+{
+ def original_ring = basering;
+ int n = nvars(original_ring);
+ int prime_field=npars(basering);
+ if(prime_field){
+ string minpolystr = "minpoly="+
+ get_minpoly_str(size(original_ring),parstr(original_ring,1))+";" ;
+ }
+
+ if(size(#)){
+ int newvars = #[1];
+
+ }else{
+ int newvars = nvars(original_ring);
+ }
+ string newvarstr = "v(1.."+string(newvars)+")";
+ def newring = Ring::changevar(newvarstr, original_ring);
+ setring newring;
+ if( prime_field ){
+ execute(minpolystr);
+ }
+ return(newring);
+}
+
+static proc defaultIdeal()
+{
+ ideal I;
+ for(int i=1; i<=nvars(basering); i++){
+ I[i] = var(i)^size(basering)-var(i);
+ }
+ return( std(I) );
+}
+
+static proc order_of_extension()
+{
+ int oe=1;
+ list rl = ringlist(basering);
+ if( size(rl[1]) <> 1){
+ oe = deg( subst(minpoly,par(1),var(1)) );
+ }
+ return(oe);
+}
+
+
+static proc vec2elm(intvec v)
+{
+ number g = 1;
+ if(npars(basering) == 1){ g=par(1); }
+ number e=0;
+ int oe = size(v);
+ for(int i=1; i<=oe; i++){
+ e = e+v[i]*g^(oe-i);
+ }
+ return(e);
+}
+
+static proc random_vector()
+{
+ int c = char(basering);
+ intvec v = 0:order_of_extension();
+ for(int i=1; i<=size(v); i++){
+ v[i] = random(0, c-1);
+ }
+ return(v);
+}
+
+static proc increment(list l, list #)
+{
+ int c, i, j, oe;
+ oe = order_of_extension();
+ c = char(basering);
+
+ if( size(#) == 1 ){
+ i = #[1];
+ }else{
+ i = size(l);
+ }
+
+ l[i] = nextVec(l[i]);
+ while( l[i][1] == c && i>1 ){
+ l[i] = 0:oe;
+ i--;
+ l[i] = nextVec(l[i]);
+ }
+ if( i < size(l) ){
+ for(j=i+1; j<=size(l); j++){
+ l[j] = 0:oe;
+ }
+ }
+ return(l);
+}
+
+static proc nextVec(intvec l)
+{
+ int c, i, j;
+ i = size(l);
+ c = char(basering);
+ l[i] = l[i] + 1;
+ while( l[i] == c && i>1 ){
+ l[i] = 0;
+ i--;
+ l[i] = l[i] + 1;
+ }
+ return(l);
+}
+
+
+
+static proc every_vector()
+{
+ list element, list_of_elements;
+
+ for(int i=1; i<=nvars(basering); i++){
+ element[i] = 0:order_of_extension();
+ }
+
+ while(size(list_of_elements) < size(basering)^nvars(basering)){
+ list_of_elements = list_of_elements + list(element);
+ element = increment(element);
+ }
+ for(int i=1; i<=size(list_of_elements); i++){
+ for(int j=1; j<=size(list_of_elements[i]); j++){
+ list_of_elements[i][j] = vec2elm(list_of_elements[i][j]);
+ }
+ }
+ return(list_of_elements);
+}
+
+static proc num2int(number a)
+{
+ int N=0;
+ if(order_of_extension() == 1){
+ N = int(a);
+ if(N<0){
+ N = N + char(basering);
+ }
+ }else{
+ ideal C = coeffs(subst(a,par(1),var(1)),var(1));
+ for(int i=1; i<=ncols(C); i++){
+ int c = int(C[i]);
+ if(c<0){ c = c + char(basering); }
+ N = N + c*char(basering)^(i-1);
+ }
+ }
+ return(N);
+}
+
+static proc listnum2int(list L){
+ int N=0;
+ int q = size(basering);
+ for(int i=1; i<=nvars(basering); i++){
+ N = N + num2int(number(L[i]))*q^(nvars(basering)-i);
+ }
+ return(N);
+}
+
+static proc get_minpoly_str(int size_of_ring, string parname)
+{
+ def original_ring = basering;
+ ring new_ring = (size_of_ring, A),x,lp;
+ string S = string(minpoly);
+ string SMP;
+ if(S=="0"){
+ SMP = SMP+parname;
+ }else{
+ for(int i=1; i<=size(S); i++){
+ if(S[i]=="A"){
+ SMP = SMP+parname;
+ }else{
+ SMP = SMP+S[i];
+ }
+ }
+ }
+ return(SMP);
+}
+
+static proc sort_ideal(ideal I)
+{
+ ideal OI;
+ int i,j,M;
+ poly P;
+ M = ncols(I);
+ OI = I;
+ for(i=2; i<=M; i++){
+ j=i;
+ while(size(OI[j-1])>size(OI[j])){
+ P = OI[j-1];
+ OI[j-1] = OI[j];
+ OI[j] = P;
+ j--;
+ if(j==1){ break; }
+ }
+ }
+ return(OI);
+}
+
+static proc add_if_new(list L, ideal I)
+{
+ int i, newelement;
+ poly P;
+
+ for(i=1; i<=nvars(basering); i++){
+ P = P + reduce(var(i),I)*var(1)^(i-1);
+ }
+ newelement=1;
+ for(i=1; i<=size(L); i++){
+ if(L[i]==P){
+ newelement=0;
+ break;
+ }
+ }
+ if(newelement){
+ L = insert(L, P);
+ }
+ return(L,newelement);
+}
--
an open source computer algebra system
More information about the debian-science-commits
mailing list