果树问题讨论:这两个问题等价么？

mathe · 发表于 2008-8-21 09:59:19

Now let's summarize the progress on orchard planting problem with 4 trees per row. And the latest version of C++ code used to solve the problem is provided in the next two posts: 102# and 103# We know that optimal packing of $K_4$'s into a $K_n$ is an upbound for the orchard planting problem, and $u(n)=[{[{n-1}/3]n}/4]$ is upbound for both of the two problems. First, A C++ program is written to try to enumerate all non-equivalent method to packing at least u(n)-"MAX_FREEDOM" $K_4$ into $K_n$ in attachment file at 11#, the result for $12<=n<=14$ with MAX_FREEDOM=3 is first privided. Next in 33#, the result for n=11 with 6 rows is transformed into equations and WxMaxima used to solve the correspondent equations. We could find a solution for planting 6 rows with total 11 trees (4 trees per row) . Also WxMaxima helped to verify that for n=12, we could plant 6 rows of trees; but we could not plant 8 rows of trees. But WxMaxima complains that those equations for 7 rows are too complex so that it gives up. In 38# I found that if we project one row of trees into line at infinity (so that there'll be four infinity points and four groups of parallel straight lines in the result), it could greatly reduce the complexity of the equations. We could even manually solve all those equations for n=13 with at least 11 rows trees. With the help of WxMaxima, it is proved that we could plant at most 9 rows of trees (4 trees per row) when there're total 13 trees in 42# And at the sametime, we could find a solution to plant total 13 trees in 9 rows. (It is used to check the correctness of the C++ code used to find all method to packing $K_4$ into $K_n$ too). After projecting one rows of trees into infinity line and select another two lines as axises in the affine plane, the equations are 2 order multivariate polynomial (each variable is x or y coordinate of a tree) And we could know in advance that all of those variables are non-zero since the coordinate of all trees in X or Y axises is known in advance and they're not used as variables in the equations. So if we could find a linear combination of those equations and each non-zero term of the linear combination contains a given variable, we could divide the given variable from the linear combination to form a linear relationship among those variables. In 51# an algorithm is designed to use the idea to simplify/prove no solution/solve those nonlinear equations [But it may fail in some times]. and in 82# it is used to prove that we could plant at most 9 rows with 13 trees and at most 10 rows with 14 trees. And in 83# an example to plant10 rows with 14 trees is found by the program (used to verify correctness of the code). And in 84# it is proved that we could plant at most 12 rows with 15 trees. In 67# the C++ source code to enumerate all non-equivalent method to pack at leat u(n)-"MAX_FREEDOM" $K_4$ into $K_n$ is attached. A C++ template class node_edge_set is defined to generate a unique representation for a method to pack $K_4$ into $K_n$: i)First we will classify all nodes in the graph into mulitiple groups according degrees of the node and attach all those nodes with the correspondent group id (in the ascedent order of degrees) ii)We could classify all those $K_4$ in the graph according to the group id of nodes in the $K_4$ and sorting those groups by lexicographic ordering of group id's in the $K_4$ iii)We could classify all nodes in the graph by group id of $K_4$ using the node and sorting those gropus by lexicographic ordering of group id of $K_4$ too iv) repeat step ii) and iii) until the number of groups remains unchanged. v) If the multiple nodes are assign with same group id, select a group with mulitple nodes (heristic method to select the one with minimal group id) for each node in the group, attach a special feature to the node(let the node uses a unique group id) and repeat step iii)-V) until each node is assigned with a unique group id. For each choice of the node, we could result in a representation of the packing $K_4$ into $K_n$ configuration. (Each node is represented by the group id) Select the one which is minimal in lexicographic ordering as the unique representation of the method to packing $K_4$ into $K_n$. The C++ program could search for the problem up to15 trees. But it will fail when there're more trees because a temp file generated by the program will be more than 2G. In 93# the program is improved so that it will use multiple files to keep those temp data. and in 96# the program is successfully generate all method to pack at least 16 $K_4$ into $K_16$. (And it will generate temp files up to 10G) and the C++ equation simplifer(the 2nd program) proved that none of them could form a correspondent solution to plant no less than 16 rows with 16 trees. So we have proved that we could plant at most 15 rows with 16 trees. According to the result in wolfram, we could result in that the 13 to 16 items of A006065 are 9,10,12,15. A bug in source found at Ot 23, 2008. After the bug fix, the new code verified that the items 13 to 16 remains unchanged. In 167#, a result of planting 15 trees in 12 rows provided. In 171#and 172#, 2 result of planting 16 trees in 15 rows provided. And the coordinates of another two solutions are provided. And all solutions for planting 16 treesin 15 rows are here if result under projection is treated as same result.

mathe · 发表于 2008-8-21 10:00:53

The C++ source code to enumerate all non-equivalent method to packing at least u(n)-"MAX_FREEDOM" $K_4$ into $K_n$ Changing macro NUM_NODES to other value to search result for differrent 'n'. Make sure you have enough free hard disk space to run the code with larger n or MAX_FREEDOM. I have only tested the code with NODES_PER_EDGE equal to 4.

#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <dirent.h>
#include <algorithm>
#include <vector>
#include <set>
using namespace std;
template <int TNODES,int NODES_PER_EDGE=4>
class node_edge_set{
public:
static const int max_edges_per_node = (TNODES-1)/(NODES_PER_EDGE-1);
static const int max_total_edges = (((TNODES-1)/(NODES_PER_EDGE-1))*TNODES)/NODES_PER_EDGE;
struct final_edge_set{
static int num_edges;
char edges[max_total_edges][NODES_PER_EDGE];
static void set_curr_num_edges(int n){num_edges=n;}
bool operator<(const final_edge_set& s)const{
int i,j;
for(i=0;i<num_edges;i++)for(j=0;j<NODES_PER_EDGE;j++){
if(edges[i][j]<s.edges[i][j])return true;
if(edges[i][j]>s.edges[i][j])return false;
}
return false;
}
bool operator==(const final_edge_set& s)const{
int i,j;
for(i=0;i<num_edges;i++)for(j=0;j<NODES_PER_EDGE;j++){
if(edges[i][j]!=s.edges[i][j])return false;
}
return true;
}
void show()const{
int i,j;
for(i=0;i<num_edges;i++){
for(j=0;j<NODES_PER_EDGE;j++){
printf("%c",'A'+edges[i][j]);
}
printf(" ");
}
printf("\n");
}
void show(char *outbuf)const{
int i,j;
for(i=0;i<num_edges;i++){
for(j=0;j<NODES_PER_EDGE;j++){
outbuf[i*NODES_PER_EDGE+j]='A'+edges[i][j];
}
}
outbuf[i*NODES_PER_EDGE]='\0';
}
void show(int step_id)const{
char tname[100];
char obuf[1000];
char tbuf[1000];
int i,j;
int checksum=0;
for(i=0;i<num_edges;i++)for(j=0;j<NODES_PER_EDGE;j++){
checksum+=edges[i][j]*(i+1);
}
sprintf(tname,"data/t%d/%d",step_id,checksum);
FILE *pOutput=fopen(tname,"a+");
if(pOutput==NULL){
fprintf(stderr,"Cannot write to file %s\n",tname);
exit(-1);
}
for(i=0;i<num_edges;i++){
for(j=0;j<NODES_PER_EDGE;j++){
tbuf[i*NODES_PER_EDGE+j]='A'+edges[i][j];
}
}
tbuf[(num_edges)*NODES_PER_EDGE]='\n';
tbuf[(num_edges)*NODES_PER_EDGE+1]='\0';
if(fseek(pOutput,-(num_edges*NODES_PER_EDGE+2),SEEK_END)==0){
fread(obuf,num_edges*NODES_PER_EDGE+2,1,pOutput);
if(obuf[0]=='\n'&&strncmp(obuf+1,tbuf,(num_edges-2)*NODES_PER_EDGE)==0){
fclose(pOutput);
return;
}
}else if(fseek(pOutput,0,SEEK_SET)==0){
fread(obuf,num_edges*NODES_PER_EDGE+1,1,pOutput);
if(strncmp(obuf,tbuf,num_edges*NODES_PER_EDGE+1)==0){
fclose(pOutput);
return;
}
}
fseek(pOutput,0,SEEK_END);
fwrite(tbuf,num_edges*NODES_PER_EDGE+1,1,pOutput);
fclose(pOutput);
}
void show(FILE *pFile)const{
int i,j;
for(i=0;i<num_edges;i++){
for(j=0;j<NODES_PER_EDGE;j++){
fprintf(pFile,"%c",'A'+edges[i][j]);
}
}
fprintf(pFile,"\n");
}
};
struct node_group_type{
char groupid;
char group_counter;
char tedges;
char eids[max_edges_per_node];
bool operator<(const node_group_type& n)const{
int i;
if(tedges<n.tedges)return true;
if(tedges>n.tedges)return false;
for(i=0;i<tedges;i++){
if(eids[i]<n.eids[i])return true;
if(eids[i]>n.eids[i])return false;
}
return false;
}
bool operator==(const node_group_type& n)const{
int i;
if(tedges<n.tedges)return false;
if(tedges>n.tedges)return false;
for(i=0;i<tedges;i++){
if(eids[i]<n.eids[i])return false;
if(eids[i]>n.eids[i])return false;
}
return true;
}
void show()const{
int i;
printf("%d*{ ",group_counter);
if(eids[0]==-1){
printf("D=%d ",tedges);
}else{
for(i=0;i<tedges;i++){
printf("%d ",eids[i]);
}
}
printf("}");
}
};
struct edge_group_type{
char groupid;
char group_counter;
char nids[NODES_PER_EDGE];
bool operator<(const edge_group_type& e)const{
int i;
for(i=0;i<NODES_PER_EDGE;i++){
if(nids[i]<e.nids[i])return true;
if(nids[i]>e.nids[i])return false;
}
return false;
}
bool operator==(const edge_group_type& e)const{
int i;
for(i=0;i<NODES_PER_EDGE;i++){
if(nids[i]<e.nids[i])return false;
if(nids[i]>e.nids[i])return false;
}
return true;
}
void show()const{
int i;
printf("%d*{ ",group_counter);
for(i=0;i<NODES_PER_EDGE;i++)printf("%d ",nids[i]);
printf("}");
}
};
private:
int num_of_edges;///Number of total edges. The number of total nodes is TNODES
char edges_info[max_total_edges][NODES_PER_EDGE];///nodes in each edge.
char nodes_info[TNODES][max_edges_per_node];///edges in each node
char nodes_degree[TNODES];///number of edges in each node
class group_infos{
private:
int edge_groups_count;///number of edge groups
int node_groups_count;///number of node groups
node_group_type node_groups[TNODES];///the info of each group of node
edge_group_type edge_groups[max_total_edges];///the info of each group of edge
char node_group_id[TNODES];///The node group id of a node
char edge_group_id[max_total_edges];///The edge group id of an edge
node_edge_set *node_edges;
public:
void show_edges()const{
int i;
for(i=0;i<edge_groups_count;i++){
edge_groups[i].show();
}
printf("\n");
printf("id:{ ");
for(i=0;i<node_edges->num_of_edges;i++){
printf("%d ",edge_group_id[i]);
}
printf("}\n");
}
bool simple_case(int node_id)const{
int i,j;
if(node_edges->nodes_degree[node_id]==0)
return true;
if(node_edges->nodes_degree[node_id]>1)
return false;
int ngid=node_group_id[node_id];
int e=node_edges->nodes_info[node_id][0];
for(j=0;j<NODES_PER_EDGE;j++){
int n=node_edges->edges_info[e][j];
if(node_group_id[n]==ngid)continue;
if(node_groups[node_group_id[n]].group_counter>1)
break;
}
return j==NODES_PER_EDGE;
}
void show_nodes()const{
int i;
for(i=0;i<node_groups_count;i++){
node_groups[i].show();
}
printf("\n");
printf("id:{ ");
for(i=0;i<TNODES;i++){
printf("%d ",node_group_id[i]);
}
printf("}\n");
}
group_infos(node_edge_set *ne):node_edges(ne){}
group_infos(const group_infos& g){
int i;
node_edges=g.node_edges;
edge_groups_count=g.edge_groups_count;
node_groups_count=g.node_groups_count;
for(i=0;i<node_groups_count;i++)node_groups[i]=g.node_groups[i];
for(i=0;i<edge_groups_count;i++)edge_groups[i]=g.edge_groups[i];
for(i=0;i<TNODES;i++)node_group_id[i]=g.node_group_id[i];
for(i=0;i<node_edges->num_of_edges;i++)edge_group_id[i]=g.edge_group_id[i];
}
int get_edge_groups_count()const{return edge_groups_count;}
int get_node_groups_count()const{return node_groups_count;}
int get_node_group_id(int nid)const{
return node_group_id[nid];
}
int get_edge_group_id(int eid)const{
return edge_group_id[eid];
}
int get_edge_group_counter(int gid)const{
return edge_groups[gid].group_counter;
}
void assign_node_edge_id(char new_node_id[TNODES], int& next_node_id){
int i;
for(i=0;i<TNODES;i++){
if(new_node_id[i]==-1){
if(node_groups[node_group_id[i]].group_counter==1){
new_node_id[i]=next_node_id++;
}
}
}
}
int undetermined_node_group()const{
int i;
int smaller=TNODES+1;
int si=-1;
for(i=0;i<node_groups_count;i++){
if(node_groups[i].tedges>0&&node_groups[i].group_counter>1){
if(node_groups[i].group_counter<smaller){
smaller-node_groups[i].group_counter;
si=i;
}
}
}
return si;
}
private:
void sort_node_group()///sort node_groups
{
int i;
char tmp2[TNODES];
for(i=0;i<node_groups_count;i++)node_groups[i].groupid=i;
sort(&node_groups[0],&node_groups[node_groups_count]);
for(i=0;i<node_groups_count;i++)tmp2[node_groups[i].groupid]=i;
for(i=0;i<TNODES;i++)node_group_id[i]=tmp2[node_group_id[i]];
}
void sort_edge_group()///sort edge_groups
{
int i;
char tmp2[max_total_edges];
for(i=0;i<edge_groups_count;i++)edge_groups[i].groupid=i;
sort(&edge_groups[0],&edge_groups[edge_groups_count]);
for(i=0;i<edge_groups_count;i++)tmp2[edge_groups[i].groupid]=i;
for(i=0;i<node_edges->num_of_edges;i++)edge_group_id[i]=tmp2[edge_group_id[i]];
}
void set_edge_group()///reset edge group according to node group
{
int i,j;
edge_groups_count=0;
for(i=0;i<node_edges->num_of_edges;i++){
edge_group_type eg;
for(j=0;j<NODES_PER_EDGE;j++){
eg.nids[j]=node_group_id[node_edges->edges_info[i][j]];
}
sort(&eg.nids[0],&eg.nids[NODES_PER_EDGE]);
for(j=0;j<edge_groups_count;j++){
if(eg==edge_groups[j])
break;
}
if(j<edge_groups_count){///Matchs
edge_group_id[i]=j;
edge_groups[j].group_counter++;
}else{
edge_groups[edge_groups_count]=eg;
edge_groups[edge_groups_count].group_counter=1;
edge_group_id[i]=edge_groups_count++;
}
}
sort_edge_group();
}
void set_node_group()
{
int i,j;
char old_id[TNODES];
for(i=0;i<TNODES;i++)old_id[i]=node_group_id[i];
node_groups_count=0;
for(i=0;i<TNODES;i++){
node_group_type ng;
ng.tedges=node_edges->nodes_degree[i];
for(j=0;j<ng.tedges;j++){
ng.eids[j]=edge_group_id[node_edges->nodes_info[i][j]];
}
sort(&ng.eids[0],&ng.eids[ng.tedges]);
for(j=0;j<node_groups_count;j++){
if(ng==node_groups[j]&&
old_id[i]==node_groups[j].groupid)break;
}
if(j<node_groups_count){
node_group_id[i]=j;
node_groups[j].group_counter++;
}else{
node_groups[node_groups_count]=ng;
node_groups[node_groups_count].group_counter=1;
node_groups[node_groups_count].groupid = old_id[i];///save old id
node_group_id[i]=node_groups_count++;
}
}
sort_node_group();
}
public:
void promote_searched(int searched){
int old_node_id[TNODES];
int node_id_map[TNODES];
int i;
if(searched<=0)return;
int gcount=0;
for(i=0;i<TNODES;i++){
old_node_id[i]=node_group_id[i];
}
sort(&old_node_id[0],&old_node_id[searched]);
sort(&old_node_id[searched],&old_node_id[TNODES]);
for(i=1;i<TNODES;i++){
if(old_node_id[i]!=old_node_id[i-1]){
old_node_id[++gcount]=old_node_id[i];
}
}
gcount++;
for(i=0;i<gcount;i++){
node_id_map[old_node_id[i]]=i;
}
for(i=0;i<gcount;i++)old_node_id[i]=-1;
for(i=0;i<gcount;i++){
if(old_node_id[i]==-1){///Order not changed.
node_group_type backupn=node_groups[i];
int start=i;
if(node_id_map[start]!=start){
do{
int next=node_id_map[start];
node_group_type t=node_groups[next];
node_groups[next]=backupn;
old_node_id[start]=1;
backupn=t;
start=next;
}while(start!=i);
}
}
}
for(i=0;i<TNODES;i++){
node_group_id[i]=node_id_map[node_group_id[i]];
}
set_edge_group();
}
void init_group_info(int searched){
int i,j;
node_groups_count=0;
///first initialize node groups via degrees;
#ifdef OUTPUT
printf("start init group info\n");
#endif
for(i=0;i<searched;i++){
for(j=0;j<node_groups_count;j++){
if(node_groups[j].tedges==node_edges->nodes_degree[i])
break;
}
if(j<node_groups_count){
node_group_id[i]=j;
node_groups[j].group_counter++;
}else{
node_groups[node_groups_count].group_counter=1;
node_groups[node_groups_count].tedges=node_edges->nodes_degree[i];
node_groups[node_groups_count].eids[0]=-1;
node_group_id[i]=node_groups_count++;
}
}
int searched_groups=node_groups_count;
for(;i<TNODES;i++){
for(j=searched_groups;j<node_groups_count;j++){
if(node_groups[j].tedges==node_edges->nodes_degree[i])
break;
}
if(j<node_groups_count){
node_group_id[i]=j;
node_groups[j].group_counter++;
}else{
node_groups[node_groups_count].group_counter=1;
node_groups[node_groups_count].tedges=node_edges->nodes_degree[i];
node_groups[node_groups_count].eids[0]=-1;
node_group_id[i]=node_groups_count++;
}
}
sort_node_group();
#ifdef OUTPUT
printf("Init node info:\n");
show_nodes();
#endif
int old_node_group;
do{
old_node_group=node_groups_count;
set_edge_group();
#ifdef OUTPUT
printf("edge info:\n");
show_edges();
#endif
set_node_group();
#ifdef OUTPUT
printf("node info:\n");
show_nodes();
#endif
}while(node_groups_count!=old_node_group);
#ifdef OUTPUT
printf("end init\n");
#endif
}
bool unique_in_group(int nid){
int group=node_group_id[nid];
return (node_groups[group].group_counter==1);
}
void split_node(int nid){
int old_group=node_group_id[nid];
int i;
node_groups[old_group].group_counter--;
node_groups[node_groups_count].group_counter=1;
node_groups[node_groups_count].tedges=node_groups[old_group].tedges;
for(i=0;i<node_groups[old_group].tedges;i++){
node_groups[node_groups_count].eids[i]=node_groups[old_group].eids[i];
}
node_group_id[nid]=node_groups_count++;
int old_node_group;
do{
old_node_group=node_groups_count;
set_edge_group();
set_node_group();
}while(node_groups_count!=old_node_group);
#ifdef OUTPUT
printf("split node %d\n",nid);
printf("nodes:");
show_nodes();
set_edge_group();
printf("edges:");
show_edges();
#endif
}
};
void add_result(final_edge_set& final_edges, const group_infos& g)
{
int edges=g.get_edge_groups_count();
if(edges!=final_edges.num_edges){
int i;
for(i=0;i<final_edges.num_edges;i++){
printf("e[%d]=%d ",i,g.get_edge_group_id(i));
}
printf("\n");
for(i=0;i<edges;i++){
printf("c[%d]=%d ",i,g.get_edge_group_counter(i));
}
printf("\n");
for(i=0;i<TNODES;i++){
printf("n[%d]=%d ",i,g.get_node_group_id(i));
}
printf("\n");
g.show_edges();
throw exception();
}
int i;
final_edge_set e;
char rev[max_total_edges];
for(i=0;i<edges;i++){
rev[g.get_edge_group_id(i)]=i;
}
for(i=0;i<edges;i++){
int j;
for(j=0;j<NODES_PER_EDGE;j++){
e.edges[i][j]=g.get_node_group_id(edges_info[rev[i]][j]);
}
sort(&e.edges[i][0],&e.edges[i][NODES_PER_EDGE]);
}
#ifdef OUTPUT
printf("Find: ");
e.show();
#endif
if(e<final_edges)final_edges=e;
}
void find_all_edges(final_edge_set& final_edges,const group_infos& g, int searched)
{
int next_id = g.undetermined_node_group();
if(next_id<0){///All identify is found
group_infos ng(g);
if(searched>0){
ng.promote_searched(searched);
}
add_result(final_edges,ng);
}else{
int i;
for(i=0;i<TNODES;i++){
if(g.get_node_group_id(i)==next_id){
bool simple = g.simple_case(i);
group_infos ng(g);
ng.split_node(i);
find_all_edges(final_edges,ng,searched);
if(simple)
break;
}
}
}
}
public:
node_edge_set(int nedges, char info[][NODES_PER_EDGE]){
int i,j;
num_of_edges=nedges;
for(i=0;i<nedges;i++){
for(j=0;j<NODES_PER_EDGE;j++){
edges_info[i][j]=info[i][j];
}
}
for(i=0;i<TNODES;i++)nodes_degree[i]=0;
for(i=0;i<num_of_edges;i++){
for(j=0;j<NODES_PER_EDGE;j++){
char the_nid=info[i][j];
nodes_info[the_nid][nodes_degree[the_nid]++]=i;
}
}
}
void find_all_node_types( int searched, char group[TNODES]){
int i,j;
final_edge_set finals[TNODES];
for(i=0;i<TNODES;i++){
group[i]=-1;
for(j=0;j<num_of_edges;j++){
int k;
for(k=0;k<NODES_PER_EDGE;k++){
finals[i].edges[j][k]=TNODES;
}
}
}
group_infos g(this);
g.init_group_info(searched);
final_edge_set::set_curr_num_edges(num_of_edges);
for(i=searched;i<TNODES;i++){
if(g.unique_in_group(i))
continue;
group_infos ng(g);
ng.split_node(i);
find_all_edges(finals[i],ng,searched);
}
for(i=searched;i<TNODES;i++){
if(g.unique_in_group(i)){
group[i]=1;
continue;
}
for(j=searched;j<i;j++){
if(finals[j]==finals[i])break;
}
if(j==i){
group[i]=1;
}
}
}
final_edge_set normalize(int searched){
group_infos g(this);
g.init_group_info(searched);
int i;
final_edge_set::set_curr_num_edges(num_of_edges);
final_edge_set final_edges;
for(i=0;i<num_of_edges;i++){
int j;
for(j=0;j<NODES_PER_EDGE;j++)
final_edges.edges[i][j]=TNODES;//set it to maxiaml node
}
find_all_edges(final_edges, g,searched);
return final_edges;
}
};
template <int TNODES,int NODES_PER_EDGE>
int node_edge_set<TNODES,NODES_PER_EDGE>::final_edge_set::num_edges;
#if 0
#define NUM_EDGES 23
#define NUM_NODES 25
char init_value[NUM_EDGES][5]={
"BIFS","BOLG","BMPH","BQNJ","EBCD",
"AIET","EMOF","EQPS","EKNH","CKMA",
"CQOT","CNPG","CJHS","DQMI","DNOA",
"DPLF","DJGT","MRGS","NRFT","IORH",
"JPRA","AFGH","KQRL"
};
char init_v2[NUM_EDGES][4];
typedef node_edge_set<NUM_NODES> NESET;
int main()
{
int i,j;
for(i=0;i<NUM_EDGES;i++)for(j=0;j<4;j++){
init_v2[i][j]=init_value[i][j]-'A';
}
NESET nes(NUM_EDGES,init_v2);
NESET::final_edge_set e=nes.normalize(4);
printf("final result:\n");
e.show();
return 0;
}
#endif
#define NUM_NODES 16
#define NODES_PER_EDGE 4
#define MAX_FREEDOM 4
#define MAX_EDGES_PER_NODE ((NUM_NODES-1)/(NODES_PER_EDGE-1))
#define MAX_TOTAL_EDGES ((MAX_EDGES_PER_NODE*NUM_NODES)/NODES_PER_EDGE)
#define MAX_FREE_EDGES (MAX_FREEDOM*NODES_PER_EDGE)
char chessboard[NUM_NODES][NUM_NODES];
char used[NUM_NODES];
char line_buf[MAX_TOTAL_EDGES][NODES_PER_EDGE];
int line_count;
int cur_free_edge;
int cur_max_degree;
typedef node_edge_set<NUM_NODES,NODES_PER_EDGE> NESET;
void reorder_node(int orig_node, int new_node)
{
int i,j;
char tmp;
if(orig_node==new_node)return;
for(i=0;i<NUM_NODES;i++){
tmp=chessboard[orig_node][i];
chessboard[orig_node][i]=chessboard[new_node][i];
chessboard[new_node][i]=tmp;
}
for(i=0;i<NUM_NODES;i++){
tmp=chessboard[i][orig_node];
chessboard[i][orig_node]=chessboard[i][new_node];
chessboard[i][new_node]=tmp;
}
tmp=used[orig_node];
used[orig_node]=used[new_node];
used[new_node]=tmp;
for(i=0;i<line_count;i++){
for(j=0;j<NODES_PER_EDGE;j++){
if(line_buf[i][j]==orig_node){
line_buf[i][j]=new_node;
}else if(line_buf[i][j]==new_node){
line_buf[i][j]=orig_node;
}
}
}
}
int last_line_count;
NESET::final_edge_set last_e;
void output_cur(int step)
{
int i;
int changed=0;
for(i=step;i<NUM_NODES;i++){
if(used[i]>used[step-1])
return;
}
if(last_line_count!=line_count){
changed=1;
}
last_line_count=line_count;
{
NESET nes(line_count, line_buf);
NESET::final_edge_set e=nes.normalize(step);
if(!changed&&e==last_e)
return;
last_e=e;
e.show(step);
}
}
void init()
{
int i,j;
line_count=0;
for(i=0;i<NUM_NODES;i++)for(j=0;j<NUM_NODES;j++){
chessboard[i][j]=0;
}
for(i=0;i<NUM_NODES;i++)used[i]=0;
}
char test(int x,int y,int z, int w)
{
if(chessboard[x][y]||chessboard[x][z]||
chessboard[x][w]||chessboard[y][z]||
chessboard[y][w]||chessboard[z][w])
return 0;
return 1;
}
void set_to(int x,int y, int z,int w, int value)
{
chessboard[x][y]=chessboard[x][z]=chessboard[x][w]=
chessboard[y][x]=chessboard[y][z]=chessboard[y][w]=
chessboard[z][x]=chessboard[z][y]=chessboard[z][w]=
chessboard[w][x]=chessboard[w][y]=chessboard[w][z]=value;
}
void setv(int x,int y,int z,int w)
{
set_to(x,y,z,w,1);
used[x]++;
used[y]++;
used[z]++;
used[w]++;
line_buf[line_count][0]=x;
line_buf[line_count][1]=y;
line_buf[line_count][2]=z;
line_buf[line_count][3]=w;
line_count++;
}
void unsetv(int x,int y,int z,int w)
{
set_to(x,y,z,w,0);
used[x]--;
used[y]--;
used[z]--;
used[w]--;
line_count--;
}
void create_step_one()
{
char filename[100];
mkdir("data/s1",-1);
FILE *pFile=fopen("data/s1/0","w");
int i,j;
for(i=MAX_EDGES_PER_NODE;i>=1;i--){
if((MAX_EDGES_PER_NODE-i)*NUM_NODES>MAX_FREE_EDGES)
break;
for(j=0;j<i;j++){
fprintf(pFile,"A%c%c%c",j*3+'B',j*3+'C',j*3+'D');
}
fprintf(pFile,"\n");
}
fclose(pFile);
}
void search(int u)
{
int v,w,z;
if((MAX_EDGES_PER_NODE-used[u])*(NUM_NODES-u)<=cur_free_edge){
output_cur(u+1);
}
if(used[u]>=cur_max_degree)return;
for(v=u+1;v<NUM_NODES;v++){
if(used[v]>=cur_max_degree)continue;
if(chessboard[u][v])continue;
for(w=v+1;w<NUM_NODES;w++){
if(used[w]>=cur_max_degree)continue;
if(chessboard[u][w]||chessboard[v][w])continue;
for(z=w+1;z<NUM_NODES;z++){
if(used[z]>=cur_max_degree)continue;
if(chessboard[u][z]||chessboard[v][z]||chessboard[w][z])
continue;
setv(u,v,w,z);
search(u);
unsetv(u,v,w,z);
}
}
}
}
void process_one_line(int step, char *input)
{
int i;
init();
for(i=0;;i++){
char c=input[i*4];
if(c>='A'&&c<='Z'){
setv(input[i*4]-'A',input[i*4+1]-'A',input[i*4+2]-'A',input[i*4+3]-'A');
}else{
break;
}
}
cur_free_edge=0;
cur_max_degree=0;
for(i=0;i<step-1;i++){
cur_free_edge+=MAX_EDGES_PER_NODE - used[i];
if(cur_max_degree<used[i]){
cur_max_degree=used[i];
}
}
for(;i<NUM_NODES;i++){
if(used[i]>cur_max_degree)
return;
}
cur_free_edge = MAX_FREE_EDGES - cur_free_edge;
if(cur_free_edge<0){
fprintf(stderr,"Warning, invalid input\n");
return;
}
char group[NUM_NODES];
for(i=0;i<NUM_NODES;i++)group[i]=1;
/* if(step<(NUM_NODES-1)/3){
NESET nes(line_count, line_buf);
NESET::final_edge_set e=nes.normalize(step);
nes.find_all_node_types(step-1,group);
}*/
for(i=step-1;i<NUM_NODES;i++){
if(group[i]>0){
reorder_node(i,step-1);
search(step-1);
reorder_node(i,step-1);
}
}
}
void sort_step(int step)
{
char tname[100];
char cmdline[200];
sprintf(tname,"data/t%d",step);
DIR *dirs=opendir(tname);
if(dirs==NULL){
fprintf(stderr,"Failed to open dir %s\n",tname);
return;
}
struct dirent *flist;
while(flist=readdir(dirs)){
if(flist->d_name[0]<'0'||flist->d_name[0]>'9')
continue;
sprintf(cmdline,"sort -u data/t%d/%s -o data/s%d/%s",step,flist->d_name, step, flist->d_name);
system(cmdline);
}
closedir(dirs);
sprintf(cmdline,"rm -rf data/t%d",step);
system(cmdline);
}
#define LINE_BUF_SIZE (MAX_TOTAL_EDGES*(NODES_PER_EDGE+1)+20)
void create_step(int step)
{
char sname[100];
char line[LINE_BUF_SIZE];
FILE *pRead, *pWrite;
sprintf(sname,"data/t%d",step);
mkdir(sname,-1);
sprintf(sname,"data/s%d",step);
mkdir(sname,-1);
sprintf(sname,"data/s%d",step-1);
DIR *dirs=opendir(sname);
if(dirs==NULL){
fprintf(stderr,"Failed to open dir %s\n",sname);
return;
}
struct dirent *flist;
while(flist=readdir(dirs)){
if(flist->d_name[0]<'0'||flist->d_name[0]>'9')
continue;
sprintf(sname,"data/s%d/%s",step-1,flist->d_name);
pRead=fopen(sname,"r");
if(pRead==NULL){
fprintf(stderr,"Cannot read input file %s\n",sname);
exit(-1);
}
while(fgets(line,LINE_BUF_SIZE,pRead)){
process_one_line(step,line);
}
fclose(pRead);
}
closedir(dirs);
sort_step(step);
}
int main(char argc, char *argv[])
{
int step=atoi(argv[1]);
if(step<1){
fprintf(stderr,"step should be no less than 1\n");
return -1;
}
if(step==1){
mkdir("data",-1);
create_step_one();
}else{
create_step(step);
}
return 0;
}

复制代码

And a perl script is used to repeat the code by NUM_NODES times. Let's assuming that the binary file name is pm

#!/usr/bin/perl
\$NUM_NODES=16
for(\$i=1;\$i<=\$NUM_NODES;\$i++){
system("./pm \$i");
}
system("cat data/s\$NUM_NODES/* > data/output.all");

复制代码

BUGS found on Sep 10,2008, #define MAX_FREE_EDGES (MAX_FREEDOM*NODES_PER_EDGE) should be changed to #define MAX_FREE_EDGES (MAX_FREEDOM*NODES_PER_EDGE+MAX_EDGES_PER_NODE*NUM_NODES-MAX_TOTAL_EDGES*NODES_PER_EDGE) and it is lucky that the two macro generate same values for n between 12 and 16

mathe · 发表于 2008-8-21 10:02:17

The C++ source code to simplify/prove no solution/solve candidate equations for orchard planting problem. Changin the macro NUM_NODES to correspondent value too when trying to solve equations for different number of trees.

#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <gmp.h>
#include <vector>
using namespace std;
#define ASSERT(x)
typedef char cinfo[2];
class ssolver{
short _n;//number of variables
short _m;//number of equations
short _t;//original number of equations
int m;///number of columns
int lc;//left column count
bool _valid;
short **columnid;
short *left_column;
cinfo *column_info;
mpq_t tmp1,tmp2;
mpq_t **_matrix;
char *row_var;
char *row_evar;
short *ctmp;
short *ctmp2;
void init_memory(){
int i,j;
m=_n*(_n+1)/2;
columnid=new short *[_n];
columnid[0]=new short[_n*_n];
for(i=1;i<_n;i++)columnid[i]=columnid[i-1]+_n;
memset(columnid[0],0,sizeof(short)*m);
left_column=new short[m];
column_info=new cinfo[m];
_matrix = new mpq_t*[_m];///_m equations
row_var =new char[_m];
row_evar = new char[_m];
ctmp =new short[m];
ctmp2=new short[m];
_matrix[0]=new mpq_t[m*_m];
for(i=0;i<_m;i++){
if(i>0)
_matrix[i]=_matrix[i-1]+m;
for(j=0;j<m;j++){
mpq_init(_matrix[i][j]);
}
}
for(i=0;i<_m;i++){row_var[i]=-1;row_evar[i]=-1;}
mpq_init(tmp1);mpq_init(tmp2);
set_column_index();
_t=_m;
}
void set_column_index(){
int i,j,s=0;
for(i=0;i<_n;i++){
for(j=0;j<=i;j++){
column_info[s][0]=j;
column_info[s][1]=i;
columnid[i][j]=columnid[j][i]=s;
s++;
}
}
}
void init_left_columns(){
int i,j;
for(i=0,lc=0;i<m;i++){
for(j=0;j<_m;j++){
if(mpq_sgn(_matrix[j][i])!=0)
break;
}
if(j<_m){
left_column[lc++]=i;
}
}
_valid=true;
}
void remove_empty_rows(){
int i,j;
for(i=0;i<_m;i++){
for(j=0;j<lc;j++){
int cid=left_column[j];
if(mpq_sgn(_matrix[i][cid])!=0){
break;
}
}
if(j==lc){
swap_row(i,_m-1);
i--;
_m--;
}
}
}
int try_reduce_var(int varid);
int reduce_iterator();
int reduce_inorder();
void simplify_result(){
int i,j;
for(i=_m;i<_t;i++){
if(row_evar[i]==-1)continue;///empty line
int mvar=row_var[i];
int evar=row_evar[i];
int c=0;
for(j=0;j<_n;j++){
if(j==evar)continue;
int cid=columnid[mvar][j];
if(mpq_sgn(_matrix[i][cid])==0)continue;
ctmp2[c]=j;
ctmp[c++]=cid;
}
for(j=i+1;j<_t;j++){
int mmvar=row_var[j];
if(mmvar==-1)continue;
int ncid=columnid[mmvar][evar];
int k;
if(mpq_sgn(_matrix[j][ncid])==0)continue;
mpq_set(tmp1,_matrix[j][ncid]);
mpq_set_ui(_matrix[j][ncid],0,1);///clear it to zero
for(k=0;k<c;k++){
int ecid=columnid[mmvar][ctmp2[k]];
mpq_mul(tmp2,tmp1,_matrix[i][ctmp[k]]);
mpq_sub(_matrix[j][ecid],_matrix[j][ecid],tmp2);
}
}
}
if(row_evar[_t-1]!=-1){
for(i=_m;i<_t-1;i++){
if(row_evar[i]==-1){
swap_row(i, _t-1);///Make sure the last line is empty
break;
}
}
}
}
int eliminate_row(int row, int cid){///return v to be eliminated, return -1 if invalid row find
int i;
char mvar = row_var[row];
int evar;
int ecid;
if(column_info[cid][0]==mvar){
evar=column_info[cid][1];
}else{
ASSERT(column_info[cid][1]==mvar);
evar=column_info[cid][0];
}
row_evar[row]=evar;
mpq_set(tmp1,_matrix[row][cid]);
mpq_set_ui(_matrix[row][cid],1,1);
ecid=cid;
int c=0;
for(i=0;i<lc;i++){
cid=left_column[i];
if(cid==ecid)continue;
if(mpq_sgn(_matrix[row][cid])!=0){
ctmp[c++]=cid;
if(column_info[cid][0]==mvar){
ctmp2[c-1]=column_info[cid][1];
}else{
ASSERT(column_info[cid][1]==mvar);
ctmp2[c-1]=column_info[cid][0];
}
mpq_div(_matrix[row][cid],_matrix[row][cid],tmp1);///normalize the line
}
}
if(c==0){
_valid=false;
#ifdef OUTPUT
printf("Invalid row %d\n",row);
output();
#endif
return -1;///no solution
}
if(row!=_m-1){
swap_row(row,_m-1);
}
row=--_m;
for(i=0;i<lc;i++){
int cid=left_column[i];
int ovar=-1;
if(column_info[cid][0]==evar&&column_info[cid][1]==evar){
int u,v,j;
for(j=0;j<_m;j++){
mpq_set(tmp1, _matrix[j][cid]);
mpq_set_ui(_matrix[j][cid],0,1);//clear to zero
if(mpq_sgn(tmp1)!=0){
for(u=0;u<c;u++)for(v=u;v<c;v++){
int ncid=columnid[ctmp2[u]][ctmp2[v]];
mpq_mul(tmp2,_matrix[row][ctmp[u]],_matrix[row][ctmp[v]]);
if(u!=v){
mpq_add(tmp2,tmp2,tmp2);
}
mpq_mul(tmp2,tmp2,tmp1);
mpq_add(_matrix[j][ncid],_matrix[j][ncid],tmp2);
}
}
}
continue;
}else if(column_info[cid][0]==evar){
ovar=column_info[cid][1];
}else if(column_info[cid][1]==evar){
ovar=column_info[cid][0];
}
if(ovar>=0){
int u,j;
for(j=0;j<_m;j++){
mpq_set(tmp1,_matrix[j][cid]);
mpq_set_ui(_matrix[j][cid],0,1);
if(mpq_sgn(tmp1)!=0){
for(u=0;u<c;u++){
int ncid=columnid[ctmp2[u]][ovar];
mpq_mul(tmp2, _matrix[row][ctmp[u]],tmp1);
mpq_sub(_matrix[j][ncid],_matrix[j][ncid],tmp2);
}
}
}
}
}
init_left_columns();
}
void swap_row(int x,int y){
if(x==y)return;
int i;
for(i=0;i<lc;i++){
int cid=left_column[i];
mpq_set(tmp1, _matrix[x][cid]);
mpq_set(_matrix[x][cid],_matrix[y][cid]);
mpq_set(_matrix[y][cid],tmp1);
}
char c=row_var[x];
row_var[x]=row_var[y];
row_var[y]=c;
c=row_evar[x];
row_evar[x]=row_evar[y];
row_evar[y]=c;
}
public:
ssolver(short n,short m){_n=n,_m=m;init_memory();}
~ssolver(){
int i,j;
mpq_clear(tmp1);mpq_clear(tmp2);
for(i=0;i<_t;i++){
for(j=0;j<m;j++){
mpq_clear(_matrix[i][j]);
}
}
delete []_matrix[0];
delete []_matrix;
delete []columnid[0];
delete []columnid;
delete []left_column;
delete []column_info;
delete []row_var;
delete []row_evar;
delete []ctmp;
delete []ctmp2;
}
void output()const{
int i,j;
if(_m>0){
for(i=0;i<_n;i++){
ctmp[i]=0;
}
for(i=0;i<lc;i++){
int cid=left_column[i];
ctmp[column_info[cid][0]]=1;
ctmp[column_info[cid][1]]=1;
}
printf("solve([");
for(i=0;i<_m;i++){
int s;
if(i!=0)printf(",");
for(j=0;j<lc;j++){
int cid=left_column[j];
if((s=mpq_sgn(_matrix[i][cid]))!=0){
if(s>0){
printf("+");
}
mpq_out_str(stdout,10,_matrix[i][cid]);
if(column_info[cid][0]>0){
if(column_info[cid][0]&1){//X
printf("*%c_X",column_info[cid][0]/2+'A');
}else{
printf("*%c_Y",(column_info[cid][0]-2)/2+'A');
}
}
if(column_info[cid][1]>0){
if(column_info[cid][1]&1){//X
printf("*%c_X",column_info[cid][1]/2+'A');
}else{
printf("*%c_Y",(column_info[cid][1]-2)/2+'A');
}
}
}
}
}
printf("],[");
int start=1;
for(i=1;i<_n;i++){
if(ctmp[i]){
if(start)start=0;
else printf(",");
if(i&1){//X
printf("%c_X",i/2+'A');
}else{
printf("%c_Y",(i-2)/2+'A');
}
}
}
printf("])\n");
}
for(i=_m;i<_t;i++){
int varid=row_var[i];
if(varid==-1)continue;
for(j=0;j<_n;j++){
int cid=columnid[varid][j];
int s;
if((s=mpq_sgn(_matrix[i][cid]))!=0){
if(s>0)printf("+");
mpq_out_str(stdout,10,_matrix[i][cid]);
if(j>0){
if(j&1){//X
printf("*%c_X",j/2+'A');
}else{
printf("*%c_Y",(j-2)/2+'A');
}
}
}
}
printf("\n");
}
}
bool isValid()const{return _valid;}
void insert_item(int row, int v1, int v2, mpq_t value){
int cid=columnid[v1][v2];
mpq_add(_matrix[row][cid],_matrix[row][cid],value);
}
void simplify();
};
#define MAKE_ITEM(item_1, item_2) ((((item_1)+1)<<16)|((item_2)+1))
#define GET_ITEM_1(x) (((x)>>16)-1)
#define GET_ITEM_2(x) (((x)&0xFFFF)-1)
#define MAKE_SAFE_ITEM(item_1, item_2) ((item_1)>(item_2)?MAKE_ITEM(item_1,item_2):MAKE_ITEM(item_2,item_1))
#define NUM_NODES 16
#define NODES_PER_EDGE 4
#define MAX_EDGES_PER_NODE ((NUM_NODES-1)/(NODES_PER_EDGE-1))
#define MAX_TOTAL_EDGES ((MAX_EDGES_PER_NODE*NUM_NODES)/NODES_PER_EDGE)
char chessboard[NUM_NODES][NUM_NODES];
char used[NUM_NODES];
char line_buf[MAX_TOTAL_EDGES][NODES_PER_EDGE];
int line_count;
int cur_free_edge;
int cur_max_degree;
int curr_comb_id;
int curr_line_num;
int last_comb_id;
int last_line_num;
void output_cur()
{
int i,j;
for(i=0;i<line_count;i++)for(j=0;j<NODES_PER_EDGE;j++){
printf("%c",'A'+line_buf[i][j]);
}
printf("\n");
}
void init()
{
int i,j;
line_count=0;
for(i=0;i<NUM_NODES;i++)for(j=0;j<NUM_NODES;j++){
chessboard[i][j]=-1;
}
for(i=0;i<NUM_NODES;i++)used[i]=0;
curr_comb_id=0;
curr_line_num=0;
}
void init_mem(){
}
char test(int x,int y,int z, int w)
{
if(chessboard[x][y]>=0||chessboard[x][z]>=0||
chessboard[x][w]>=0||chessboard[y][z]>=0||
chessboard[y][w]>=0||chessboard[z][w]>=0)
return 0;
return 1;
}
void set_to(int x,int y, int z,int w, int value)
{
chessboard[x][y]=chessboard[x][z]=chessboard[x][w]=
chessboard[y][x]=chessboard[y][z]=chessboard[y][w]=
chessboard[z][x]=chessboard[z][y]=chessboard[z][w]=
chessboard[w][x]=chessboard[w][y]=chessboard[w][z]=value;
}
void setv(int x,int y,int z,int w)
{
set_to(x,y,z,w,line_count);
used[x]++;
used[y]++;
used[z]++;
used[w]++;
line_buf[line_count][0]=x;
line_buf[line_count][1]=y;
line_buf[line_count][2]=z;
line_buf[line_count][3]=w;
line_count++;
}
void unsetv(int x,int y,int z,int w)
{
set_to(x,y,z,w,-1);
used[x]--;
used[y]--;
used[z]--;
used[w]--;
line_count--;
}
enum value_type{ SET_CONST,UNSET,SET_INDEX,SET_INFTY};
void sort_by_feature(char nodes[NODES_PER_EDGE], int fx[NUM_NODES],int fy[NUM_NODES]){
int i,j;
for(i=0;i<NODES_PER_EDGE;i++){
for(j=i+1;j<NODES_PER_EDGE;j++){
if(fx[nodes[i]]>fx[nodes[j]]||
fx[nodes[i]]==fx[nodes[j]]&&fy[nodes[i]]>fy[nodes[j]]){
char t=nodes[i];
nodes[i]=nodes[j];
nodes[j]=t;
}
}
}
}
void output_node_x(char a, int fx[],int x[])
{
if(fx[a]==UNSET){
printf("%c_x",a+'A');
}else if(fx[a]==SET_CONST){
printf("%d",x[a]);
}else{
fprintf(stderr,"Invalid \n");
exit(-1);
}
}
void output_node_y(char a, int fy[],int y[])
{
if(fy[a]==UNSET){
printf("%c_y",a+'A');
}else if(fy[a]==SET_CONST){
printf("%d",y[a]);
}else{
fprintf(stderr,"Invaid\n");
exit(-1);
}
}
int ssolver::try_reduce_var(int varid)
{
int cur_row=0;
int i,j,k;
#ifdef OUTPUT
if(varid==0){
printf("try reduce constant\n");
}else if(varid&1){
printf("try reduce %c_X\n",'A'+varid/2);
}else{
printf("try reduce %c_Y\n",'A'+(varid-2)/2);
}
#endif
for(i=0;i<lc;i++){
int cid=left_column[i];
if(column_info[cid][0]==varid||column_info[cid][1]==varid)continue;
for(j=cur_row;j<_m;j++){
if(mpq_sgn(_matrix[j][cid])!=0)
break;
}
if(j==_m)continue;
if(j!=cur_row){
swap_row(j,cur_row);
}
mpq_set(tmp1, _matrix[cur_row][cid]);
for(j=0;j<lc;j++){///normalize the line
int cid2=left_column[j];
if(mpq_sgn(_matrix[cur_row][cid2])!=0){
mpq_div(_matrix[cur_row][cid2],_matrix[cur_row][cid2],tmp1);
}
}
for(k=0;k<_m;k++){///searching all rows except this
if(k==cur_row)continue;
if(mpq_sgn(_matrix[k][cid])!=0){
mpq_set(tmp1,_matrix[k][cid]);
for(j=0;j<lc;j++){
int cid2=left_column[j];
if(mpq_sgn(_matrix[cur_row][cid2])!=0){
mpq_mul(tmp2,tmp1,_matrix[cur_row][cid2]);
mpq_sub(_matrix[k][cid2],_matrix[k][cid2],tmp2);
}
}
}
}
cur_row++;
if(cur_row==_m)break;
}
if(cur_row==_m)return 0;
for(i=cur_row;i<_m;i++){
row_var[i]=varid;
}
int fr=0;
int end_row=cur_row;
for(i=_m-1;i>=end_row;i--){
int j;
int have_const=0,have_varid=0;
for(j=0;j<lc;j++){
int cid=left_column[j];
if(varid!=0&&(column_info[cid][0]==0||column_info[cid][1]==0)){
if(mpq_sgn(_matrix[i][cid])!=0){
have_const=1;
}
}else if(column_info[cid][0]==varid&&column_info[cid][1]==varid){
if(mpq_sgn(_matrix[i][cid])!=0){
have_varid=1;
}
}else if(mpq_sgn(_matrix[i][cid])!=0)break;
}
if(j==lc){
if(!have_const&&!have_varid){///empty row
if(i!=_m-1){
swap_row(i,_m-1);
}
row_var[_m-1]=-1;
--_m;
}else if(!have_const||!have_varid){
return -1;
}else{
if(i!=end_row){
swap_row(i,end_row);
}
end_row++;
}
continue;
}
int r=eliminate_row(i,left_column[j]);
if(r==-1)return -1;
fr=1;
}
if(cur_row<end_row){///Need to eliminate varid
int cid_const=columnid[0][varid];
int cid_var=columnid[varid][varid];
mpq_div(tmp1, _matrix[cur_row][cid_const], _matrix[cur_row][cid_var]);
for(i=cur_row+1;i<end_row;i++){
mpq_div(tmp2,_matrix[i][cid_const],_matrix[i][cid_var]);
if(mpq_cmp(tmp1,tmp2)!=0)
return -1;
}
int r=eliminate_row(end_row-1,cid_var);
if(r==-1)return -1;
fr=1;
if(end_row>cur_row+1){
remove_empty_rows();
}
}
#ifdef OUTPUT
printf("After reduce var\n");
output();
#endif
return fr;
}
int ssolver::reduce_inorder()
{
int i;
int max_c=-1,max_index=-1;
for(i=0;i<_n;i++){
ctmp[i]=0;
}
for(i=0;i<lc;i++){
int cid=left_column[i];
if(column_info[cid][0]==column_info[cid][1]){
ctmp[column_info[cid][0]]++;
}else{
ctmp[column_info[cid][0]]++;
ctmp[column_info[cid][1]]++;
}
}
for(i=0;i<_n;i++){
if(ctmp[i]>0){
int r=try_reduce_var(i);
if(r!=0)return r;
}
}
return 0;
}
int ssolver::reduce_iterator()
{
int i;
int max_c=-1,max_index=-1;
for(i=0;i<_n;i++){
ctmp[i]=0;
}
for(i=0;i<lc;i++){
int cid=left_column[i];
if(column_info[cid][0]==column_info[cid][1]){
ctmp[column_info[cid][0]]++;
}else{
ctmp[column_info[cid][0]]++;
ctmp[column_info[cid][1]]++;
}
}
for(i=0;i<_n;i++){
if(ctmp[i]>max_c){
max_c=ctmp[i];
max_index=i;
}
}
return try_reduce_var(max_index);
}
ssolver *solver;
void add_item_x(int x_n, int fx[], int x[], bool sub)
{
if(fx[x_n]==UNSET){
mpq_t q;mpq_init(q);mpq_set_ui(q,1,1);
if(sub)mpq_neg(q,q);
solver->insert_item(curr_line_num, 2*x_n+1,0,q);
mpq_clear(q);
}else if(fx[x_n]==SET_CONST){
mpq_t q;mpq_init(q);mpq_set_ui(q,x[x_n],1);
if(sub)mpq_neg(q,q);
solver->insert_item(curr_line_num, 0,0, q);
mpq_clear(q);
}else{
fprintf(stderr,"invaid\n");
}
}
void add_item_y(int y_n, int fy[], int y[], bool sub)
{
if(fy[y_n]==UNSET){
mpq_t q;mpq_init(q);mpq_set_ui(q,1,1);
if(sub)mpq_neg(q,q);
solver->insert_item(curr_line_num, 2*y_n+2, 0, q);
mpq_clear(q);
}else if(fy[y_n]==SET_CONST){
mpq_t q;mpq_init(q);mpq_set_ui(q,y[y_n],1);
if(sub)mpq_neg(q,q);
solver->insert_item(curr_line_num, 0, 0, q);
mpq_clear(q);
}else{
fprintf(stderr,"invaid\n");
}
}
void add_an_item(int x_n, int y_n, int fx[], int x[], int fy[], int y[] , bool sub)
{
if(fx[x_n]==UNSET&&fy[y_n]==UNSET){
mpq_t q;
mpq_init(q);mpq_set_ui(q,1,1);
if(sub)mpq_neg(q,q);
solver->insert_item(curr_line_num, 2*x_n+1, 2*y_n+2,q);
mpq_clear(q);
}else if(fx[x_n]==UNSET&&fy[y_n]==SET_CONST){
mpq_t q;mpq_init(q);mpq_set_ui(q,y[y_n],1);
if(sub)mpq_neg(q,q);
solver->insert_item(curr_line_num, 2*x_n+1,0,q);
mpq_clear(q);
}else if(fx[x_n]==SET_CONST&&fy[y_n]==UNSET){
mpq_t q;mpq_init(q);mpq_set_ui(q,x[x_n],1);
if(sub)mpq_neg(q,q);
solver->insert_item(curr_line_num, 2*y_n+2,0,q);
mpq_clear(q);
}else if(fx[x_n]==SET_CONST&&fy[y_n]==SET_CONST){
mpq_t q;mpq_init(q);mpq_set_ui(q,x[x_n]*y[y_n],1);
if(sub)mpq_neg(q,q);
solver->insert_item(curr_line_num, 0,0, q);
mpq_clear(q);
}else{
fprintf(stderr,"Invalid \n");
exit(-1);
}
}
#define ADD_AN_ITEM(xn,yn) add_an_item(xn,yn,fx,x,fy,y,false)
#define SUB_AN_ITEM(xn,yn) add_an_item(xn,yn,fx,x,fy,y,true)
#define ADD_X(xn) add_item_x(xn,fx,x,false)
#define SUB_X(xn) add_item_x(xn,fx,x,true)
#define ADD_Y(yn) add_item_y(yn,fy,y,false)
#define SUB_Y(yn) add_item_y(yn,fy,y,true)
void output_line(char a,char b,char c, int fx[],int fy[],int x[],int y[]){
if(fx[c]!=SET_INFTY){///Limit case
ADD_AN_ITEM(b,a);
ADD_AN_ITEM(c,b);
ADD_AN_ITEM(a,c);
SUB_AN_ITEM(b,c);
SUB_AN_ITEM(c,a);
SUB_AN_ITEM(a,b);
}else{
if(fy[c]==SET_INFTY){///parallel to x or y coord
if(x[c]==1&&y[c]==0){///y=const
ADD_Y(a);
SUB_Y(b);
}else if(x[c]==0&&y[c]==1){
ADD_X(a);
SUB_X(b);
}else{
fprintf(stderr,"invalid data");
exit(-1);
}
}else{///The ratio
ADD_Y(b);
SUB_Y(a);
SUB_AN_ITEM(b,c);
ADD_AN_ITEM(a,c);
}
}
curr_line_num++;
}
void ssolver::simplify()
{
init_left_columns();
#ifdef OUTPUT
printf("Before remove empty rows\n");
output();
#endif
remove_empty_rows();
#ifdef OUTPUT
printf("After init\n");
output();
#endif
do{
int r=reduce_iterator();
if(r==-1){
_valid=false;
return;
}
if(r==1)continue;
r=reduce_inorder();
if(r==-1){
_valid=false;
return;
}
if(r==0)break;
}while(1);
if(_valid){
simplify_result();
}
}
void process_one_line(char *input)
{
int i,j,k;
init();
for(i=0;;i++){
char c=input[i*4];
if(c>='A'&&c<='Z'){
setv(input[i*4]-'A',input[i*4+1]-'A',input[i*4+2]-'A',input[i*4+3]-'A');
}else{
break;
}
}
int max_degree=0;
char e_to_n[MAX_TOTAL_EDGES][MAX_TOTAL_EDGES];
char n_to_e[NUM_NODES][NUM_NODES];
for(i=0;i<MAX_TOTAL_EDGES;i++)for(j=0;j<MAX_TOTAL_EDGES;j++){
e_to_n[i][j]=-1;
}
for(i=0;i<NUM_NODES;i++)for(j=0;j<NUM_NODES;j++){
n_to_e[i][j]=-1;
}
for(i=0;i<line_count;i++){
for(j=0;j<NODES_PER_EDGE;j++){
for(k=j+1;k<NODES_PER_EDGE;k++){
n_to_e[line_buf[i][j]][line_buf[i][k]]=
n_to_e[line_buf[i][k]][line_buf[i][j]]=i;
}
}
}
for(i=0;i<NUM_NODES;i++){
for(j=0;j<NUM_NODES;j++){
if(n_to_e[i][j]==-1)continue;
for(k=j+1;k<NUM_NODES;k++){
if(n_to_e[i][k]==-1)continue;
if(n_to_e[i][j]==n_to_e[i][k])continue;
e_to_n[n_to_e[i][j]][n_to_e[i][k]]=
e_to_n[n_to_e[i][k]][n_to_e[i][j]]=i;
}
}
}
int max_sum=0;
int best_i,best_j,best_k;
int a,b,c;
for(i=0;i<line_count;i++){
for(j=i+1;j<line_count;j++){
if(e_to_n[i][j]==-1)continue;
for(k=j+1;k<line_count;k++){
if(e_to_n[i][k]==-1||e_to_n[j][k]==-1)continue;
if(e_to_n[i][j]==e_to_n[i][k])continue;
int a=e_to_n[j][k];
int b=e_to_n[i][k];
int c=e_to_n[i][j];
///Try edge i
int u;
int lsum=used[a];
for(u=0;u<NODES_PER_EDGE;u++){
lsum+=used[line_buf[i][u]];
}
if(lsum>max_sum){
max_sum=lsum;
best_i=i;best_j=j;best_k=k;
}
lsum=used[b];
for(u=0;u<NODES_PER_EDGE;u++){
lsum+=used[line_buf[j][u]];
}
if(lsum>max_sum){
max_sum=lsum;
best_i=j;best_j=k;best_k=i;
}
lsum=used[c];
for(u=0;u<NODES_PER_EDGE;u++){
lsum+=used[line_buf[k][u]];
}
if(lsum>max_sum){
max_sum=lsum;
best_i=k;best_j=i;best_k=j;
}
}
}
}
if(max_sum==0){
printf("%s\nerror: No triangle found\n",input);
return;
}
solver =new ssolver(2*NUM_NODES+1, 2*line_count);
a=e_to_n[best_j][best_k];///(0,0,1)
b=e_to_n[best_i][best_k];///(1,0,0)
c=e_to_n[best_i][best_j];///(0,1,0)
int x[NUM_NODES],y[NUM_NODES];
int fx[NUM_NODES],fy[NUM_NODES];
for(i=0;i<NUM_NODES;i++){
fx[i]=fy[i]=UNSET;///value unset
}
fx[a]=fy[a]=SET_CONST;x[a]=y[a]=0;
fx[b]=fy[b]=SET_INFTY;x[b]=1;y[b]=0;
fx[c]=fy[c]=SET_INFTY;x[c]=0;y[c]=1;
for(i=0;i<NODES_PER_EDGE;i++){///in best_k, y=0
if(line_buf[best_k][i]!=a&&line_buf[best_k][i]!=b){
int n=line_buf[best_k][i];
fx[n]=fy[n]=SET_CONST;
x[n]=1;y[n]=0;
break;
}
}
for(i++;i<NODES_PER_EDGE;i++){
if(line_buf[best_k][i]!=a&&line_buf[best_k][i]!=b){
int n=line_buf[best_k][i];
fy[n]=SET_CONST;
y[n]=0;
break;
}
}
for(i=0;i<NODES_PER_EDGE;i++){///in best_j, x=0
if(line_buf[best_j][i]!=a&&line_buf[best_j][i]!=c){
int n=line_buf[best_j][i];
fx[n]=fy[n]=SET_CONST;
x[n]=0;y[n]=1;
break;
}
}
for(i++;i<NODES_PER_EDGE;i++){
if(line_buf[best_j][i]!=a&&line_buf[best_j][i]!=c){
int n=line_buf[best_j][i];
fx[n]=SET_CONST;
x[n]=0;
break;
}
}
for(i=0;i<NODES_PER_EDGE;i++){///in best_i, INFTY LINE
if(line_buf[best_i][i]!=b&&line_buf[best_i][i]!=c){
int n=line_buf[best_i][i];
fx[n]=SET_INFTY;
fy[n]=UNSET;
x[n]=1;
}
}
for(i=0;i<line_count;i++){
if(i==best_i||i==best_j||i==best_k)
continue;
int s,t,z,w;
sort_by_feature(line_buf[i],fx,fy);
s=line_buf[i][0],t=line_buf[i][1],z=line_buf[i][2],w=line_buf[i][3];
if(fx[t]==SET_INFTY||fx[z]==SET_INFTY||fx[s]==SET_INFTY){
printf("%s\nInvalid lines data\n",input);
delete solver;
return;
}
if(fx[w]==SET_INFTY){
output_line(s,t,w,fx,fy,x,y);
output_line(s,z,w,fx,fy,x,y);
}else{
output_line(s,t,z,fx,fy,x,y);
output_line(s,t,w,fx,fy,x,y);
}
}
last_comb_id = curr_comb_id;
last_line_num = curr_line_num;
solver->simplify();
if(solver->isValid()){
printf("%s\n",input);
solver->output();
for(i=0;i<NUM_NODES;i++){
if(fx[i]==SET_INFTY&&fy[i]==SET_INFTY){
printf("%c=(%d,%d,0) ",i+'A',x[i],y[i]);
}else if(fx[i]==SET_INFTY){
printf("%c=(1,%c_y,0) ",i+'A',i+'A');
}else{
if(fx[i]==SET_CONST){
printf("%c_x=%d ",i+'A',x[i]);
}
if(fy[i]==SET_CONST){
printf("%c_y=%d ",i+'A',y[i]);
}
}
}
printf("\n");
}
delete solver;
}
#define LINE_BUF_SIZE (MAX_TOTAL_EDGES*(NODES_PER_EDGE+1)+20)
void process_file(char *filename)
{
char line[LINE_BUF_SIZE];
FILE *pRead;
pRead=fopen(filename,"r");
if(pRead==NULL){
fprintf(stderr,"Cannot read input file %s\n",filename);
exit(-1);
}
while(fgets(line,LINE_BUF_SIZE,pRead)){
process_one_line(line);
}
fclose(pRead);
}
int main(char argc, char *argv[])
{
init_mem();
process_file(argv[1]);
return 0;
}

复制代码

0→∞ · 发表于 2008-8-21 10:37:59

mathe可以发表论文了

0→∞ · 发表于 2008-8-21 10:48:50

为什么不允许游客访问呢，禁止他们回复和发帖还不行么

mathe · 发表于 2008-8-21 10:52:02

原帖由 0→∞ 于 2008-8-21 10:48 发表为什么不允许游客访问呢，禁止他们回复和发帖还不行么

游客，如果您要查看本帖隐藏内容请回复

发表论文我兴趣不大，比较花时间，网络上发表结果就可以了。

mathe · 发表于 2008-8-21 11:16:39

不过现在都只是证明了一些人家找出来的结果已经是最优结果，没有找出更好的结果。如果能够对n=17找出16条直线或者n=20找出24条直线的，那就更加好了。另外zgg和没一一问题有没有试一试用我那个解方程的程序过滤一下你们关于问题二找出来的一些候选解？

0→∞ · 发表于 2008-8-21 12:00:38

关于问题2，能否算出在20点24线时共有多少种不同类的解？

无心人 · 发表于 2008-8-21 12:49:58

连存在都不能证明不可能给出20点4点共线的24直线解的数量的

mathe · 发表于 2008-8-21 13:29:20

他说的是问题2的数量。这个关于17点不少于17条边情况的求解在我的计算机上产生几十G数据以后好像让我的Linux机器崩溃了，现在我不能远程登录了

账号		自动登录	找回密码
密码			欢迎注册

[求助] 果树问题讨论:这两个问题等价么？