魔术球问题

首先诚挚感谢BYvoid大神博客里提供的好题和数据

此为线性规划与网络流24题的第四题。

«问题描述：
假设有n根柱子，现要按下述规则在这n根柱子中依次放入编号为1，2，3，¼的球。
（1）每次只能在某根柱子的最上面放球。
（2）在同一根柱子中，任何2个相邻球的编号之和为完全平方数。
试设计一个算法，计算出在n根柱子上最多能放多少个球。例如，在4 根柱子上最多可
放11 个球。
«编程任务：
对于给定的n，计算在n根柱子上最多能放多少个球。
«数据输入：
由文件input.txt提供输入数据。文件第1 行有1个正整数n，表示柱子数。
«结果输出:
程序运行结束时，将n 根柱子上最多能放的球数以及相应的放置方案输出到文件
output.txt中。文件的第一行是球数。接下来的n行，每行是一根柱子上的球的编号。
输入文件示例输出文件示例
input.txt output.txt
4 　　　　11
　　　　　1 8
　　　　 2 7 9
　　　　　3 6 10
　　　　　4 5 11

建模分析：

实际上是最小路径覆盖问题，存在有向边a->b的条件是 (a<b) a+b为完全平方数.

我用匈牙利和dinic两种办法做了下。dinic：

  1 #include<cstdio>
  2 #include<cstring>
  3 #include<cstdlib>
  4 #include<queue>
  5 #include<iostream>
  6 #define rep(i,a,b) for(int i=a;i<+b;++i)
  7 using namespace std;
  8 const int MAXN=5010;
  9 const int maxn=1000;
 10 const int INF=~0U>>1;
 11 int n,f,tot,S,T;
 12 int h[MAXN*2];
 13 int square[MAXN],pointer[MAXN*2],que[MAXN],path[MAXN];
 14 bool visit[MAXN];
 15 void square_number()
 16 {
 17     square[0]=maxn;
 18     rep(i,1,maxn)
 19     {
 20         square[i]=i*i;
 21     }
 22 }
 23 struct edge{
 24     int to,next;
 25     int cap,f;
 26     int op;
 27 }edge[MAXN*100];
 28 void addedge(int a,int b)
 29 {
 30     edge[tot].to=b;
 31     edge[tot].next=pointer[a];
 32     edge[tot].cap=1;
 33     edge[tot].op=tot^1;
 34     edge[tot].f=0;
 35     pointer[a]=tot;
 36     tot++;
 37     edge[tot].to=a;                 //建立反向边
 38     edge[tot].next=pointer[b];
 39     edge[tot].cap=0;
 40     edge[tot].op=tot^1;
 41     edge[tot].f=0;
 42     pointer[b]=tot;
 43     tot++;
 44 }
 45 void newpoint(int k)
 46 {
 47     rep(i,1,square[0])
 48     {
 49         if(square[i]>=2*k) break;
 50         if(square[i]>k)
 51         {
 52             addedge(square[i]-k,MAXN+k);
 53         }
 54     }
 55     addedge(S,k);
 56     addedge(MAXN+k,T);
 57 }
 58 void Input()
 59 {
 60     scanf("%d",&n);
 61     memset(pointer,-1,sizeof(pointer));
 62     S=0;T=2*MAXN-1;
 63     square_number();
 64 }
 65 
 66 bool dinic_restruct()       //bfs建层次图
 67 {
 68     memset(h,-1,sizeof(h));
 69     queue <int>q;
 70     h[S]=0;
 71     q.push(S);
 72     while(!q.empty())
 73     {
 74         int u=q.front();q.pop();
 75         for(int j=pointer[u];j!=-1;j=edge[j].next)
 76         {
 77             int v=edge[j].to;
 78             if(h[v]==-1&&edge[j].cap-edge[j].f)
 79             {
 80                 h[v]=h[u]+1;
 81                 q.push(v);
 82                 if(v==T) return true;
 83             }
 84         }
 85     }
 86     return false;
 87 
 88 }
 89 void dinic_augment()            //dfs找增广路
 90 {
 91     int Stop=1;
 92     que[Stop]=S;
 93     while(Stop)
 94     {
 95         int u=que[Stop];
 96         if(u!=T)
 97         {
 98             int j;
 99             for(j=pointer[u];j!=-1;j=edge[j].next)
100             {
101                 int v=edge[j].to;
102                 if(h[v]==h[u]+1&&edge[j].cap-edge[j].f>0)
103                 {
104                     que[++Stop]=v;
105                     break;
106                 }
107             }
108             if(j!=-1)
109             {
110                 path[Stop]=j;       //记录边的序号
111             }
112             else            //回退
113             {
114                 Stop--;
115                 h[u]=-1;
116             }
117         }
118         else
119         {
120             int flow=INF;
121             for(int j=Stop;j>=2;--j)
122             {
123                 int &tmp=path[j];
124                 if(edge[tmp].cap-edge[tmp].f<flow) {
125                         flow=edge[tmp].cap-edge[tmp].f;
126                 }
127             }
128            // printf("flow=%d f=%d
",flow,f);
129             f+=flow;
130             for(int j=Stop;j>=2;--j)
131             {
132                 int &tmp=path[j];
133                 edge[tmp].f+=flow;
134                 edge[tmp^1].f-=flow;
135                 if(edge[tmp].cap-edge[tmp].f==0) Stop=j-1;            }
136         }
137     }
138 }
139 void Dinic(int i)
140 {
141     newpoint(i);
142     while(dinic_restruct())
143     {
144         dinic_augment();
145     }
146 }
147 void Pout(int ans)
148 {
149     printf("%d
",ans);
150     memset(path,0,sizeof(path));
151     rep(i,1,ans)
152     {
153         for(int j=pointer[i];j!=-1;j=edge[j].next)
154         {
155             if(edge[j].f==1&&edge[j].to!=S)
156             {
157                     path[i]=edge[j].to;
158                     break;
159             }
160         }
161     }
162     rep(i,1,ans)
163     {
164         int t=i;
165         if(!visit[i])
166         {
167             while(path[t])
168             {
169                 printf("%d ",t);
170                 t=path[t]-MAXN;
171                 visit[t]=1;
172             }
173             printf("%d
",t);
174         }
175     }
176 }
177 void work()
178 {
179     f=0;
180     rep(i,1,maxn*maxn)
181     {
182         Dinic(i);
183         if(i-f>n)
184         {
185             Pout(i-1);
186             break;
187         }
188     }
189 }
190 int main()
191 {
192     freopen("ball10.in","r",stdin);
193     Input();
194     work();
195     return 0;
196 }

匈牙利：

  1 #include<cstdio>
  2 #include<cstring>
  3 #define rep(i,a,b) for(int i=a;i<=b;++i)
  4 using namespace std;
  5 const int MAXN=2010;
  6 const int maxn=1000;
  7 int graph[MAXN][MAXN];          //邻接表建图
  8 int match[MAXN];
  9 int visit[MAXN];
 10 int path[MAXN];
 11 int square[maxn+10];
 12 int n,q,S,T,f,z;
 13 void square_number()
 14 {
 15     square[0]=maxn;
 16     rep(i,1,maxn)
 17     {
 18         square[i]=i*i;
 19     }
 20 }
 21 void newpoint(int k)
 22 {
 23     rep(i,1,square[0])
 24     {
 25         if(square[i]>=2*k) break;
 26         if(square[i]>k)
 27         {
 28             int t=square[i]-k;
 29             graph[t][++graph[t][0]]=k;
 30         }
 31     }
 32 
 33 
 34 }
 35 void init()
 36 {
 37     f=0;
 38     S=0;T=2*MAXN-1;
 39     square_number();
 40     scanf("%d",&n);
 41 }
 42 void pout()
 43 {
 44     printf("%d
",q-1);
 45     rep(i,1,q-1)
 46     {
 47         int  t=i;
 48         if(!visit[i])
 49         {
 50 
 51             while(path[t]) t=path[t];
 52             while(match[t])
 53             {
 54                 visit[t]=1;
 55                 printf("%d ",t);
 56                 t=match[t];
 57             }
 58             printf("%d
",t);
 59         }
 60     }
 61 
 62 }
 63 bool crosspath(int k)
 64 {
 65     rep(i,1,graph[k][0])
 66     {
 67         int &tmp=graph[k][i];
 68         if(!visit[tmp])
 69         {
 70             visit[tmp]=1;
 71             if(match[tmp]==0||match[tmp]!=k&&crosspath(match[tmp]))
 72             {
 73                 match[tmp]=k;
 74               //  printf("k=%d tmp=%d f=%d q=%d z=%d
",k,tmp,f,q,z);
 75                 path[k]=tmp;
 76                 return true;
 77             }
 78         }
 79     }
 80     return false;
 81 }
 82 void hungary()
 83 {
 84     for(q=1;q<=maxn*maxn;++q)
 85     {
 86         newpoint(q);
 87         for(z=1;z<=q-1;++z)
 88         {
 89             if(!path[z]&&crosspath(z))      //如果这个点还是单身就试试能不能从它出发找到一条交替路
 90             {
 91                     f++;
 92                    // rep(j,1,q-1)
 93                    // printf("q=%d match[%d]=%d f=%d
",q,j+MAXN,match[j],f);
 94             }
 95             memset(visit,0,sizeof(visit));
 96         }
 97         if(q-f>n)           //注意：虽然输出的时候跑过一遍q+1的流程,但match数组没有更新过，因为q+1轮中f必然没有+1，从而crosspath一直在返回0
 98         {
 99             pout();
100             break;
101         }
102     }
103 
104 }
105 
106 int main()
107 {
108     freopen("ball0.in","r",stdin);
109     init();
110     hungary();
111     return 0;
1