Défi N°4 : Combinatoire et matrices binaires

[Jean-Marc.Bourguet]

J'aurais vraissemblablement pas le temps de faire quelque chose de sérieux dans un proche avenir, donc voici mon code.

L'idée que j'aimerais explorer, c'est le fait que le DAG est symétrique si on impose que la feuille soit au niveau n. Ca devrait permettre d'éliminer a priori de ne pas calculer les feuilles qui sont aux niveaux inférieurs.

Code C++ :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
#include <exception>
#include <iostream>
#include <map>
#include <stddef.h>
#include <vector>
 
#ifndef ACCUM
typedef unsigned long long AccumType;
#else
typedef ACCUM AccumType;
#endif
 
long eval(char* s)
{
  char* last;
  if (s == NULL) {
    throw "Missing argument.";
  }
  long result = strtol(s, &last, 0);
  if (*last != '\0') {
    std::cerr << "'" << s << "' is not a valid number.\n";
    throw EXIT_FAILURE;
  }
  return result;
}
 
class Array
{
public:
  Array(int n) : value_(n) {}
  // Array(Array const& o);
  // Array& operator=(Array const& o);
  // ~Array();
  int& operator[](int n)        { return value_[n]; }
  int operator[](int n) const   { return value_[n]; }
  int size() const              { return value_.size(); }
private:
  std::vector<int> value_;
};
 
bool operator<(Array const& l, Array const& r)
{
  if (l.size() < r.size()) {
    return true;
  } else if (r.size() < l.size()) {
    return false;
  } else {
    for (int i = 0; i < l.size(); ++i) {
      if (l[i] < r[i]) {
        return true;
      } else if (r[i] < l[i]) {
        return false;
      }
    }
    return false;
  }
}
 
typedef std::map<Array, AccumType> Map;
 
AccumType comb(int n, int k)
{
  AccumType result = 1;
  for (int i = 0; i < k; ++i) {
    result *= n-i;
    result /= i+1;
  }
  return result;
}
 
void computeNext(Map const& current, Map& next, int k)
{
  Map().swap(next);
  Array state(k+1);
  for (Map::const_iterator i = current.begin(), e = current.end();
       i != e;
       ++i)
  {
    state = i->first;
    Array diff(k+1);
    diff[0] = 0;
 
    int pos = 1;
    int remainder = k;
    int available = 0;
    for (int j = 1; j <= k; ++j) {
      available += state[j];
    }
    int last = k+1;
    AccumType count = i->second;
 
    while (pos > 0) {
      int cur = std::min(std::min(remainder, state[pos]), last-1);
 
      diff[pos] = cur;
      available -= state[pos];
      count *= comb(state[pos], cur);
      state[pos] -= cur;
      state[pos-1] += cur;
      remainder -= cur;
 
      pos++;
      last = k+1;
 
      if (pos > k && remainder == 0) {
        next[state] += count;
        assert(remainder == 0);
      }
      if ((pos > k) || (remainder > available)) {
        while ((pos > 0) && ((remainder >= available) || (last == 0))) {
          pos--;
          last = diff[pos];
          remainder += last;
          state[pos] += last;
          state[pos-1] -= last;
          count /= comb(state[pos], last);
          available += state[pos];
        }
      }
    }
  }
}
 
int main(int argc, char* argv[])
{
  int status = EXIT_SUCCESS;
  try {
    int n = eval(argv[1]);
    int k = eval(argv[2]);
 
    Array state(k+1);
    Map current;
    Map next;
 
    state[k] = n;
    next[state] = 1;
 
    for (int i = 0; i < n; ++i) {
      current.swap(next);
      computeNext(current, next, k);
    }
    assert(next.size() == 1);
    std::cout << k << ", " << n << " -> "
              << next.begin()->second << " solutions\n";
  } catch (std::exception const& e) {
    std::cerr << e.what() << std::endl;
    status = EXIT_FAILURE;
  } catch (char const* m) {
    std::cerr << m << std::endl;
    status = EXIT_FAILURE;
  } catch (int s) {
    status = s;
  } catch (...) {
    std::cerr << "Unknown exception." << std::endl;
    status = EXIT_FAILURE;
  }
  return status;
}

[pseudocode]

La version récursive java ressemble ENORMEMENT, même dans le nom des variables

Code java :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
 
/**
 *  Défi 4 : Combinatoire et matrices binaires
 *
 * @author xavier philippeau
 */
public class Defi4 {
 
	int N=0, K=0;
	int[] state = null;
 
	// C(n,p) calculation (using a cache)
	long[][] CNP = null;
	long Cnp(int n,int p) {
		if (CNP[n][p]==0) { // not in cache -> compute
			long c = 1, fact=1;
			for(int i=0;i<p;i++) {c*=(n-i); fact*=i+1;}
			CNP[n][p] = c/fact;
		}
		return CNP[n][p];
	}
 
	long explore(int remainingtoken, long currentRowPossibilities, int starttoken) {
		// done all rows
		if (remainingtoken==0 && state[0]==N) {
			return currentRowPossibilities;
		}
 
		// done 1 row -> compute next rows, then return total
		if (remainingtoken==0 && state[0]<N) {
			long nextRowsPossibilities = explore(K,1,1);
			return currentRowPossibilities*nextRowsPossibilities;
		}		
 
		// place remaining token for the current row
		long newPossibilities = 0;
		for(int i=starttoken;i<=K;i++) {
			for(int p=Math.min(remainingtoken,state[i]);p>=1;p--) {
				long cnp = Cnp(state[i],p);
				state[i]-=p;state[i-1]+=p;
				newPossibilities += explore(remainingtoken-p,currentRowPossibilities*cnp,i+1);
				state[i-1]-=p;state[i]+=p;
			}
		}
		return newPossibilities;
	}
 
	void compute(int n, int k) {
		this.N = n;
		this.K = k;
		this.state = new int[K+1];
		this.CNP = new long[N+1][N+1];
 
		// initialize state array
		for(int i=0;i<K;i++) state[i]=0;
		state[K]=N;
 
		// explore possibilities
		long result = explore(K,1,1);
 
		System.out.println("result="+result);
	}
 
	public static void main(String[] args) {
		new Defi4().compute(9, 4);
	}
}

[Jean-Marc.Bourguet]

Envoyé par pseudocode

La version récursive java ressemble ENORMEMENT, même dans le nom des variables

Si j'ai bien suivi, tu parcours tous les chemins du DAG, j'evite de le faire en utilisant de la programmation dynamique. Tu aurais le meme resultat en memorisant les resultats des appels a explore, sans oublie que cette fonction a un parametre cache.

[pseudocode]

Envoyé par Jean-Marc.Bourguet

Si j'ai bien suivi, tu parcours tous les chemins du DAG, j'evite de le faire en utilisant de la programmation dynamique. Tu aurais le meme resultat en memorisant les resultats des appels a explore, sans oublie que cette fonction a un parametre cache.

Dans la version "originale" j'ai effectivement mis un cache d'appel à explore(), ce qui donne le temps de calcul de 35ms que j'ai annoncé précédemment.

La version que j'ai postée ici (donc sans le cache) donne 250ms.

[Jean-Marc.Bourguet]

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
#include <exception>
#include <iostream>
#include <map>
#include <stddef.h>
#include <vector>

#ifndef ACCUM
typedef unsigned long long AccumType;
#else
typedef ACCUM AccumType;
#endif

long eval(char* s)
{
  char* last;
  if (s == NULL) {
    throw "Missing argument.";
  }
  long result = strtol(s, &last, 0);
  if (*last != '\0') {
    std::cerr << "'" << s << "' is not a valid number.\n";
    throw EXIT_FAILURE;
  }
  return result;
}

class Array
{
public:
  Array(int n) : value_(n) {}
  // Array(Array const& o);
  // Array& operator=(Array const& o);
  // ~Array();
  int& operator[](int n)        { return value_[n]; }
  int operator[](int n) const   { return value_[n]; }
  int size() const              { return value_.size(); }
private:
  std::vector<int> value_;
};

bool operator<(Array const& l, Array const& r)
{
  if (l.size() < r.size()) {
    return true;
  } else if (r.size() < l.size()) {
    return false;
  } else {
    for (int i = 0; i < l.size(); ++i) {
      if (l[i] < r[i]) {
        return true;
      } else if (r[i] < l[i]) {
        return false;
      }
    }
    return false;
  }
}

typedef std::map<Array, std::pair<AccumType, AccumType> > Map;

AccumType comb(int n, int k)
{
  AccumType result = 1;
  for (int i = 0; i < k; ++i) {
    result *= n-i;
    result /= i+1;
  }
  return result;
}

void computeNext(Map const& current, Map& next, int k)
{
  Map().swap(next);
  Array state(k+1);
  for (Map::const_iterator i = current.begin(), e = current.end();
       i != e;
       ++i)
  {
    state = i->first;
    Array diff(k+1);
    diff[0] = 0;

    int pos = 1;
    int remainder = k;
    int available = 0;
    for (int j = 1; j <= k; ++j) {
      available += state[j];
    }
    int last = k+1;
    AccumType countUp = i->second.first;
    AccumType countDown = i->second.second;
    
    while (pos > 0) {
      int cur = std::min(std::min(remainder, state[pos]), last-1);
      
      diff[pos] = cur;
      available -= state[pos];
      countUp *= comb(state[pos], cur);
      state[pos] -= cur;
      state[pos-1] += cur;
      countDown *= comb(state[pos-1], cur);
      remainder -= cur;
      
      pos++;
      last = k+1;
      
      if (pos > k && remainder == 0) {
        next[state].first += countUp;
        next[state].second += countDown;
        assert(remainder == 0);
      }
      if ((pos > k) || (remainder > available)) {
        while ((pos > 0) && ((remainder >= available) || (last == 0))) {
          pos--;
          last = diff[pos];
          remainder += last;
          countDown /= comb(state[pos-1], last);
          state[pos] += last;
          state[pos-1] -= last;
          countUp /= comb(state[pos], last);
          available += state[pos];
        }
      }
    }
  }
}

AccumType consolidate(Map const& current, Map const& next, int k)
{
  AccumType result = 0;
  Array state(k+1);
  for (Map::const_iterator i = current.begin(), e = current.end();
       i != e;
       ++i)
  {
    for (int j = 0; j <= k; ++j) {
      state[j] = i->first[k-j];
    }
    Map::const_iterator n = next.find(state);
    assert(n != next.end());
    result += i->second.first * n->second.second;
  }
  return result;
}

int main(int argc, char* argv[])
{
  int status = EXIT_SUCCESS;
  try {
    int n = eval(argv[1]);
    int k = eval(argv[2]);

    Array state(k+1);
    Map current;
    Map next;

    state[k] = n;
    next[state].first = 1;
    next[state].second = 1;
    for (int i = 0; i < (n+1)/2; ++i) {
      current.swap(next);
      computeNext(current, next, k);
    }
    if (n % 2 == 0) {
      std::cout << k << ", " << n << " -> "
                << consolidate(next, next, k)
                << " solutions\n";
    } else {
      std::cout << k << ", " << n << " -> "
                << consolidate(current, next, k)
                << " solutions\n";
    }
  } catch (std::exception const& e) {
    std::cerr << e.what() << std::endl;
    status = EXIT_FAILURE;
  } catch (char const* m) {
    std::cerr << m << std::endl;
    status = EXIT_FAILURE;
  } catch (int s) {
    status = s;
  } catch (...) {
    std::cerr << "Unknown exception." << std::endl;
    status = EXIT_FAILURE;
  }
  return status;
}

Code shell :

1
2
3
4
5
6
7
[~/src]$ g++ -O3 -DACCUM="long double" defi4.cpp -o defi4
[~/src]$ time ./defi4  18 9
9, 18 -> 2.18826e+72 solutions

real    0m23.928s
user    0m23.861s
sys     0m0.008s

Je n'ai plus d'idées algorithmique. Il reste à soigner l'implémentation.

[Saint_Louis]

J'ai découvert par hasard votre forum. Je réponds en Caml Light.
Pour éviter le débordement des entiers. La réponse est donnée sous forme d'un produit.
Le dernier calcul pouvant être effectué est n=8, k=4 (environ 5 ks) sinon la réponse devient négative pour raison de débordement d'entier.
J'ai effectué une permutation de la première colonne et aussi de la première ligne. Ce qui explique la multiplication non effectuée dans la réponse.

Ma méthode consiste en deux fonctions récursives qui s'appellent mutuellement. On parcourt ainsi un graphe virtuel (programmation dynamique). Voici le code qui a l'avantage d'être court:

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
let compteur=ref 0;;

let notest x y (i,j) k p=x > k ||  y > k || y< k+i-p||x < k+j-p ;;

let prec (i,j) p=if i<> 0 then (i-1,j) else (p,j-1);;
let suiv (i,j) p=if i<> p then (i+1,j) else (0,j+1);;

let rec binomial n=fun
|0 -> 1
|p -> (binomial (n-1) (p-1))*n/p;;

let rec avance (h::q) (i,j) c l k p=
    match h with 
	|[] -> if q <> [] then (let (u,v)=prec (i,j) p in
                     l.(u) <-l.(u)- hd (hd q);c.(v) <-c.(v)- hd (hd q));
              raise Exit
             |(a::b) ->let e=c.(j) and f=l.(i) in let x1= e+a and y1=f+a in
              if (notest y1 x1 (i,j) k p) then avance (b::q) (i,j) c l k p
              else  begin
                     let (u,v)=suiv (i,j) p in
                     c.(j) <- x1;l.(i) <- y1;
                     try if (v=0 && u< k) || (u=0 && v<k) 
                         then resout ([1]::h::q) (u,v) c l k p
                         else resout ([0;1]::h::q) (u,v) c l k p
                    with Exit -> c.(j) <- e;l.(i) <- f;avance (b::q) (i,j) c l k p 
            end
and resout (h::q) (i,j) c l k p=
    match h with
        |[] ->raise Exit
        |_ when (i,j)<> (p,p) -> avance (h::q) (i,j) c l k p
        |_ ->if not(notest (l.(p)+hd h) (c.(p)+hd h) (p,p) k p) 
                then incr compteur;
            resout ((tl h)::q) (i,j) c l k p;;      

let n=6 and k=3 in
    let  c=make_vect n 0 and l=make_vect n 0 in
    try if k=0 || k=n then print_int 1 
        else resout [[1]] (0,0) c l k (n-1)
    with Exit -> let s=binomial (n-1) (k-1) in
        print_int (!compteur * s * s*n/k);;

[Sve@r]

Envoyé par Saint_Louis

J'ai découvert par hasard votre forum. Je réponds en Caml Light.

Moi aussi, j'ai découvert les défis hier et je viens de coder une soluce en Python

Envoyé par Saint_Louis

Ma méthode consiste en deux fonctions récursives qui s'appellent mutuellement. On parcourt ainsi un graphe virtuel.Voici le code qui a l'avantage d'être court:

J'ai d'abord pensé à la récursivité puis je suis parti sur un bête incrément.
En effet, imaginons une demande sur le couple (5, 3), on peut alors partir de l'hypothèse que chaque ligne n'aura que 3 jetons (comme l'a dit Jean-Marc). Mais si on part de la ligne "0 0 1 1 1", on s'aperçoit qu'il suffit d'incrémenter la ligne de 1 jusqu'à trouver la ligne suivante qui convient. Ainsi on passe par dessus "0 1 0 0 0", "0 1 0 0 1", "0 1 0 1 0" et on arrive jusqu'à "0 1 0 1 1".

Ensuite si on regarde la matrice complète, alors elle démarre de cette façon
0 0 1 1 1
0 0 1 1 1
0 0 1 1 1
0 0 1 1 1
0 0 1 1 1

Ensuite il suffit d'incrémenter la matrice, c.a.d. incrémenter la 5° ligne. Dès que la 5° ligne a atteint son max, soit "1 1 1 0 0", alors on la repasse à "0 0 1 1 1" mais on incrémente la 4° ligne qui passe alors à 0 1 0 1 1". Et ainsi de suite. Chaque ligne qui ne convient pas s'incrémente automatiquement jusqu'à convenir. Et une fois qu'on arrive sur une position finalisée (où toutes les lignes conviennent), il suffit de regarder si la somme de chaque colonne convient. Et on termine dès que les 5 lignes sont à "1 1 1 0 0"

Voici le code Python

Code Python :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
#!/usr/bin/env python
# coding: Latin-1 -*-
 
class cElem:
    # Objet pour gérer un élément
 
    # Constructeur
    def __init__(
            self,                       # Référence à l'objet (obligatoire)
            val=0):                     # Valeur
 
        self.set(val)
    # __init__()
 
    # Accesseur set
    def set(
            self,                       # Référence à l'objet (obligatoire)
            val):                       # Valeur
        self.val=val
    # set()
 
    # Accesseur get
    def get(
            self):                      # Référence à l'objet (obligatoire)
        return self.val
    # get()
 
    # L'affichage de l'élément
    def __str__(
            self):                      # Référence à l'objet (obligatoire)
 
        return "%d" % self.val
    # __str__()
# class cElem
 
# Objet pour gérer une ligne ou une colonne
class cLigCol:
    # Objet pour gérer une ligne/colonne
 
    # Constructeur
    def __init__(
            self,                       # Référence à l'objet (obligatoire)
            size,                       # Taille de l'ensemble
            but):                       # But à atteindre
 
        self.tab=[]
        self.size=size
        self.but=but
    # __init__()
 
    # Ajout d'un élément
    def append(
            self,                       # Référence à l'objet (obligatoire)
            elem):                      # Element à ajouter
 
        self.tab.append(elem)
    # append()
 
    # Initialisation ensemble
    def init(
            self):                      # Référence à l'objet (obligatoire)
        # Un ensemble initialisé commence par "k - 1" 0
        for i in range(self.size - self.but + 1):
            self.tab[i].set(0)
 
        # Un ensemble initialisé se termine par "k" 1
        for i in range(self.size - self.but, self.size):
            self.tab[i].set(1)
    # init()
 
    # Somme
    def somme(
            self):                      # Référence à l'objet (obligatoire)
        res=0
        for elem in self.tab:
            res+=elem.get()
        return res
    #somme()
 
    # Ensemble ok
    def isOk(
            self):                      # Référence à l'objet (obligatoire)
        return self.somme() == self.but
    # isOk()
 
    # Ensemble à son maximum
    def isMax(
            self):                      # Référence à l'objet (obligatoire)
 
        # Un ensemble est au max si il commence par "k" 1
        for i in range(self.but):
            if self.tab[i].get() != 1: return False
        return True
    # isMax()
 
    # Incrément ensemble
    def add(
            self):                      # Référence à l'objet (obligatoire)
 
        # Boucle infinie (sortie si ok)
        while True:
            # Incrémenter le dernier chiffre et gérer la retenue
            flag_inc=True
            i=self.size - 1
 
            # Tant qu'il faut incrémenter
            while flag_inc:
                # Si le chiffre est déjà à 1
                if self.tab[i].get() == 1:
                    # On le met à 0
                    self.tab[i].set(0)
 
                    # L'indice se déplace au chiffre précédent
                    i=i - 1
                else:
                    # On le met à 1
                    self.tab[i].set(1)
 
                    # On n'a pas de retenue
                    flag_inc=False
 
                    # on replace l'indice sur le dernier chiffre
                    i=self.size - 1
 
            # Si l'ensemble est ok on sort
            if self.isOk(): break
    # add()
 
    # Affichage de l'élément
    def __str__(
            self):                      # Référence à l'objet (obligatoire)
 
        str=""
        for elem in self.tab:
            str+="%s " % elem
        str+="(%d)" % self.somme()
        return str
    # __str__()
# class cLigCol
 
class cMatrice:
    # Objet pour gérer la matrice
 
    # Constructeur
    def __init__(
            self,                       # Référence à l'objet (obligatoire)
            size,                       # Taille de la matrice
            but,                        # Nombre à atteindre
            val=0):                     # Valeur d'initialisation
 
        # Mémorisation des paramètres
        self.size=size
        self.but=but
        self.surface=size * size
 
        # Tableau pour stocker les chiffres
        self.tabElem=[]
        for i in range(size * size):
            self.tabElem.append(cElem(val))
 
        # Tableau pour gérer les lignes
        self.tabLig=[]
        for i in range(size):
            self.tabLig.append(cLigCol(size, but))
            for j in range(size):
                self.tabLig[i].append(self.tabElem[i * size + j])
 
        # Tableau pour gérer les colonnes
        self.tabCol=[]
        for i in range(size):
            self.tabCol.append(cLigCol(size, but))
            for j in range(size):
                self.tabCol[i].append(self.tabElem[i + j * size])
 
        # Initialisation des lignes (si il ya  au-moins un but)
        if but != 0:
            for lig in self.tabLig:
                lig.init()
    # __init__()
 
    # Matrice ok
    def isOk(
            self):                      # Référence à l'objet (obligatoire)
 
        # On sait que les lignes seront ok donc on ne regarde que les colonnes
        for col in self.tabCol:
            # Si la colonne n'est pas ok
            if not col.isOk(): return False
 
        # La matrice est ok
        return True
    # isOk()
 
    # Matrice à son maximum
    def isMax(
            self):                      # Référence à l'objet (obligatoire)
 
        # La matrice est au max si toutes ses lignes sont au max
        for lig in self.tabLig:
            if not lig.isMax(): return False
        return True
    # isMax()
 
    # Incrémente la matrice
    def add(
            self):                      # Référence à l'objet (obligatoire)
 
        # Incrémenter la dernière ligne et gérer la retenue
        flag_inc=True
        i=self.size - 1
 
        # Tant qu'il faut incrémenter
        while flag_inc:
            # Si la ligne est déjà au max
            if self.tabLig[i].isMax() == 1:
                # On la réinitialise
                self.tabLig[i].init()
 
                # L'indice se déplace sur la ligne précédente
                i=i - 1
            else:
                # On l'incrémente
                self.tabLig[i].add()
 
                # On n'a pas de retenue
                flag_inc=False
 
                # on replace l'indice sur la dernière ligne
                i=self.size - 1
    # add()
 
    # Affichage de la matrice
    def __str__(
            self,                       # Référence à l'objet (obligatoire)
            verbose=False):             # Mode verbeux
 
        str=""
 
        # Matrice
        str+="\n"
        for i in range(self.size):
            for j in range(self.size):
                str+="%1s " % self.tabElem[i * self.size + j]
            str+="\n"
 
        if verbose:
            str+="\n"
            # Lignes:
            for i in range(self.size):
                str+="Ligne %d: %s\n" % (i + 1, self.tabLig[i])
            str+="\n"
 
            # Colonnes:
            for i in range(self.size):
                str+="Colonne %d: %s\n" % (i + 1, self.tabCol[i])
            str+="\n"
 
        return str
    # __str__()
 
    # Recherche des matrices correspondantes
    def find(
            self,                       # Référence à l'objet (obligatoire)
            verbose=False):             # Mode verbose
 
        res=0
        # Tant qu'on n'est pas en fin de matrice
        while not self.isMax():
            # On incrémente la matrice
            self.add()
            # Si la matrice est ok, on compte
            if self.isOk():
                res+=1
                if verbose: print "Matrice: %s" % self
        return res
    # find()
# class cMatrice
 
# Pour tester le module
import sys
if __name__ == "__main__":
    if len(sys.argv) > 2:
        # Récupération des paramètres
        size=int(sys.argv[1])
        but=int(sys.argv[2])
 
        # Les deux petits tests du crétin
        if but > size:
            print "Résultat: Aucune matrice de taille (%d/%d) ne peut donner %d\n" % (size, size, but)
            sys.exit(0)
        if but == 0 or but == size:
            if but == size:
                n=1
            else:
                n=0
            print "Résultat: Une seule matrice de taille (%d/%d) donne %d:%s" % (size, size, but, cMatrice(size, but, n))
            sys.exit(0)
 
        # Création de la matrice
        mat=cMatrice(size, but)
 
        # Recherche
        print mat.find()

[EDIT]Bon, mauvais algo. Je trouve les bonnes valeurs pour (4, 2), (5, 2) et (6, 2) mais la recherche (6, 2) a pris 5 bonnes minutes. Je vais maintenant examiner la solution de Jean-Marc...

Envoyé par Fractal LLG

(n, k) = (6, 2) -> 67 956

67950...

[EDIT 2]Super mauvais algo. J'ai lancé (9, 4) il y a plus d'une heure et toujours pas de réponse...

[Jean-Marc.Bourguet]

Envoyé par Sve@r

[EDIT 2]Super mauvais algo. J'ai lancé (9, 4) il y a plus d'une heure et toujours pas de réponse...

La réponse est 315 031 400 802 720, donc un algo qui génère les solutions et ne se contente pas de les compter en tenant compte de symétries diverses (ce que fait le mien) devrait en générer près de 10 milliards par seconde pour y arriver en une heure...

[Sve@r]

Envoyé par Jean-Marc.Bourguet

La réponse est 315 031 400 802 720, donc un algo qui génère les solutions et ne se contente pas de les compter en tenant compte de symétries diverses (ce que fait le mien) devrait en générer près de 10 milliards par seconde pour y arriver en une heure...

Ou 1 milliard par seconde pour y arriver en 10h (mon pgm n'a toujours pas trouvé actuellement). Mais comme le problème parle d'une énumération de matrices je suis parti sur une génération individuelle (avec même une option d'affichage pour contrôle). Ca dépend évidemment de la façon qu'on a de définir le verbe "énumérer"...

[Saint_Louis]

J'ai repris les idées déjà exposées par les brillants cerveaux de ce forum et j'en ai tiré un code Caml Light très rapide.
Ainsi, il calcule C(9,4) en environ 7s.
L'idée est de bien suivre les permutations de colonnes possibles. Si on tient compte de la permutation des lignes sur la première colonne, on doit pouvoir encore gagner un facteur binomial(n-1,k-1) mais je ne l'ai pas fait.
J'ai ouvert le module des grands nombres ce qui permet de faire sauter la limitation des entiers. Le programme tel qu'il est présenté est correct jusqu'à n=16 sinon il faut changer cette valeur dans le tableau des coefficients du binome. J'ai essayé de documenter mon programme mais ce n'est pas simple.
Tout est récursif sauf la construction du tableau de Pascal

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
#open "num";;(*Pour ouvrir le module des grands nombres*)

type paquet_col={largeur:int;hauteur:int};;

let binome n= let u=make_matrix (n+1) (n+1) (Int 1) in
    for i=1 to n do
        for j=1 to (i-1) do
        u.(i).(j) <- u.(i-1).(j) +/ u.(i-1).(j-1)
        done
    done;
    u;;

let cnp=binome 16;;    
    
(*Un paquet colonne est un groupe de colonnes de même hauteur, la hauteur désigne le nombre de 1 encore à placer.
La liste des paquets colonnes correspond aux paquets colonne de hauteurs toutes distinctes.
Les hauteurs vont décroitre au fur et à mesure de l'avancement de l'algorithme.
Quand une hauteur devient nulle, c'est qu'on a rempli la colonne avec les  k valeurs 1
Au depart, on a un seul paquet {largeur=n;hauteur=k}.
Aprés une opération, on aura deux paquets {largeur=k;hauteur=k-1} et {largeur=n-k;hauteur=k}
Plus loin on regroupera les colonnes car [{largeur=8;hauteur=1};{largeur=2;hauteur=1}]
est en fait la même chose que [{largeur=10;hauteur=1}]
*)
let cons a li=a::li;;
(*la fonction solution repartit la valeur k sur la liste des paquets colonnes.
On obtient ici toutes les possibilités sous forme de liste de liste d'entiers.
on vérifie en même temps que les hauteurs restent positives pour 
toute solution proposée.*)

let rec solutionbis liste_sol (pa::queue as liste_paquet) t =
        match queue with
        |[] when pa.hauteur=0 -> if t=0 then [[0]] else raise Exit
        |[] ->if t<= pa.largeur then (
            if liste_sol=[] then [[t]] else map (cons t) liste_sol
            ) else raise Exit
        |_  when pa.hauteur=0 -> map (cons 0) (solutionbis liste_sol queue t)
        |_ ->let h= if t<= pa.largeur then t else pa.largeur in
            solutionter h liste_sol liste_paquet t
            
    and
        solutionter s liste_sol liste_paquet t=match s with
        | -1 -> liste_sol
        |_ -> try let u=map (cons s) (solutionbis liste_sol (tl liste_paquet) (t-s)) in 
            u@solutionter (s-1) liste_sol liste_paquet t 
            with Exit -> solutionter (s-1) liste_sol liste_paquet t;;
    
    let solution liste_paquet t=solutionbis [] liste_paquet t;;
    
    (*Exemple: solution [ {largeur=3;hauteur=1};{largeur=1;hauteur=0};{largeur=2;hauteur=1}] 4;;
    renvoie [[3; 0; 1]; [2; 0; 2]], 
    en deuxième position, obligatoirement 0 car la hauteur du groupe était 0 donc il était rempli,
    si pas de solution renvoie Exit.
    solution [ {largeur=3;hauteur=4};{largeur=1;hauteur=3};{largeur=2;hauteur=2}] 4;;
    renvoie [[3; 1; 0]; [3; 0; 1]; [2; 1; 1]; [2; 0; 2]; [1; 1; 2]]*)

(*On regroupe maintenant les paquets colonnes contigus de même hauteur (les hauteurs croissent
avec l'indice de paquets colonne).
On utilise la fonction miroir pour garder cet ordre croissant*)
let rec miroiraux a=fun
	| [] -> a
	|(t::q) -> miroiraux (t::a) q;;

let rec regroupe=fun
    (u::v) (a::q) -> if u.hauteur=a.hauteur then regroupe ({largeur=u.largeur+a.largeur;hauteur=u.hauteur}:: v) q 
	else regroupe (a::u:: v) q
    |[] (a::q)-> regroupe [a] q
    |w [] -> miroiraux [] w;;

(* Exemple: regroupe [] [{largeur=3;hauteur=1};{largeur=1;hauteur=1};{largeur=2;hauteur=2}];;
 renvoie [{largeur = 4; hauteur = 1}; {largeur = 2; hauteur = 2}]*)


type etat={coeff:num;li_pa:paquet_col list};;

(*La fonction avance part d'un découpage possible, et d'une liste de paquet_colonnes
et renvoie une nouvelle liste de paquets_colonnes avec un coefficient multiplicateur
Cette fonction avance depose en fait une ligne de 1 suivant le découpage paramétre
Il y a donc k hauteurs qui diminuent de 1, on précise plus loin le résultat obtenu*)

let rec avance sol1 liste_paquet=match sol1 with
    |[]->{coeff=Int 1;li_pa=[]}
    |(a::q) -> let {largeur=s;hauteur=t} as v=hd liste_paquet 
        and res= avance q (tl liste_paquet) in
        if a=s then {coeff=res.coeff;li_pa={largeur=a;hauteur=t -1}::res.li_pa}
        else if a=0 then {coeff=res.coeff;li_pa=v::res.li_pa}
            else {coeff=(cnp.(s).(a)) */ res.coeff;li_pa=
                {largeur=a;hauteur=t -1}::{largeur=s -a;hauteur=t}::res.li_pa};; 


(* Exemple: avance [2;1;0] [{largeur=3;hauteur=2};{largeur=1;hauteur=3};{largeur=2;hauteur=4}];;  
renvoie - : etat =
 {coeff = 3;
  li_pa =[{largeur = 2; hauteur = 1}; {largeur = 1; hauteur = 2};
    {largeur = 1; hauteur = 2}; {largeur = 2; hauteur = 4}]} *)

(*La fonction sommehaut permettra de vérifier qu'il y a au moins k colonnes encore à remplir*)
let rec sommehaut=function
|[] ->0
|{largeur=u;hauteur=v}::q ->if v>0 then u+sommehaut q else sommehaut q;;

(* La fonction valide termine la fonction avance, elle renvoie Exit si ce n'est pas possible
sinon elle renvoie le booléen true si toutes les colonnes sont remplies et false sinon*)
let valide sol1 liste_paquet k=let {coeff=u;li_pa=t} as res=avance sol1 liste_paquet in
    let tot=sommehaut t in if tot=0 then true,res
    else if tot >=k then false,{coeff=u;li_pa=(regroupe [] t)} else raise Exit;;

(*Exemple: valide [3;5] [{largeur=4;hauteur=2};{largeur=6;hauteur=3}] 8;;
renvoie:  false, {coeff = 24;  li_pa =
   [{largeur = 3; hauteur = 1}; {largeur = 6; hauteur = 2};
    {largeur = 1; hauteur = 3}]}*)
    
let rec totpartiel liste_paquet k=function
|[] -> Int 0
|(s::queue) -> try let (boo,{coeff=coeffi;li_pa=liste_paquet_res})=valide s liste_paquet k in
if boo then coeffi else (coeffi */ total liste_paquet_res k) +/ totpartiel liste_paquet k queue 
	 with Exit -> totpartiel liste_paquet k queue
and total liste_paquet k=let u=solution liste_paquet k in
	totpartiel liste_paquet k u;;

let compter n k=total [{largeur=n;hauteur=k}] k;;

compter 9 4;;

(*compter 16 2:
num = 36574751938491748341360000
#- : temps = 12.375
compter 9 4
num = 315031400802720
#- : temps = 7.422*)

[Saint_Louis]

Pour tenir compte de la permutation possible des lignes sur la première colonne (mettre tous les 1 en bas à gauche), j'ai légèrement amélioré les deux dernières fonctions de mon précédent algorithme. J'édite ici les modifications (en Caml-Light) de deux ou trois dernières fonctions.

Code :

1
2
3
4
5
6
7
8
9
10
11
let rec totpartiel liste_paquet t q k=function
|[] -> Int 0
|(s::queue) -> try let (boo,{coeff=coeffi;li_pa=liste_paquet_res})=valide s liste_paquet t in
if boo then coeffi else (coeffi */ total liste_paquet_res (q-1) k) +/ totpartiel liste_paquet t q k queue 
	 with Exit -> totpartiel liste_paquet t q k queue
and total liste_paquet q k=if q>0 then (let u=solution liste_paquet (k-1) in
	totpartiel liste_paquet (k-1) q k u) else (let u=solution liste_paquet k in
	totpartiel liste_paquet k q k u);;

let compter n k=(cnp.(n).(k)) */ (cnp.(n-1).(k-1)) */ (total [{largeur=k-1;hauteur=k-1};{largeur=n-k;hauteur=k}] (k-1) k);;

On gagne un facteur 3 par rapport à l'ancienne mouture.
Voici quelques réponses pour k=2.
n=11:num = 158815387962000
#- : temps = 0.047
n=12:num = 21959547410077200
#- : temps = 0.11
n=13:num = 3574340599104475200
#- : temps = 0.313
n=14:num = 676508133623135814000
#- : temps = 0.828

Pour celui qui s'intéresse à la valeur k=2, cette suite a déjà été étudiée: On trouvera une relation de récurrence, un lien avec les séries entières, le comportement asymptotique et la valeur exacte à l'adresse suivante:
http://www.research.att.com/~njas/se...lish&go=Search. On y trouvera aussi de nombreuses références bibliographiques

[Sve@r]

Envoyé par Saint_Louis

Pour tenir compte de la permutation possible des lignes sur la première colonne (mettre tous les 1 en bas à gauche), j'ai légèrement amélioré les deux dernières fonctions de mon précédent algorithme. J'édite ici les modifications (en Caml-Light) de deux ou trois dernières fonctions.

Code :

On gagne un facteur 3 par rapport à l'ancienne mouture.
Voici quelques réponses pour k=2.
n=11:num = 158815387962000
#- : temps = 0.047
n=12:num = 21959547410077200
#- : temps = 0.11
n=13:num = 3574340599104475200
#- : temps = 0.313
n=14:num = 676508133623135814000
#- : temps = 0.828

Pour celui qui s'intéresse à la valeur k=2, cette suite a déjà été étudiée: On trouvera une relation de récurrence, un lien avec les séries entières, le comportement asymptotique et la valeur exacte à l'adresse suivante:
http://www.research.att.com/~njas/se...lish&go=Search. On y trouvera aussi de nombreuses références bibliographiques

Tiens... me semble reconnaître un certain saintlouis qui écume le forum maths...

[egery]

Bonjour,

voici une version qui ne vaudra que pour sa relative concision et quelques singularités (par exemple, simplification des coefficients binomiaux, de sorte qu'il n'en reste aucun à calculer !).

On repart du principe introduit par Jean-Marc :
la matrice se voit divisée en blocs verticaux. A chaque nouvelle ligne, on injecte les k "1" (les "jetons") dans les différents blocs. Avec par exemple k=9, on décompose les solutions ainsi :
- Insérer 9 jetons dans le 1er bloc, et 0 dans les autres.
- Insérer 8 jetons dans le 1er bloc, et 1 dans les autres.
- ...
- Insérer 0 jeton dans le 1er bloc, et 9 dans les autres.
Ceci se prête donc bien à une récursion où les blocs seront typiquement organisés en liste.

En outre, l'itération 9, 8, 7... n'est pas implémentée en tant que boucle. La fonction d'insertion se charge de renvoyer le nouveau itérateur.
Cela permet quelques raccourcis :
- si le bloc ne fait que 3 de largeur, impossible d'insérer 9 jetons ! On en insère 3, et on renvoie 2 comme itérateur suivant.
- si le bloc contient déjà trop de "1", on invalide le cas, et on fixe l'itérateur à 0.

Bref, voici ce que cela donne :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
(*                                    *)
(* Defi developpez.com : denombrement *)
(*                                    *)
(* denombrement-defi4.ml  (OCaml)     *)
(*                                    *)
(* (c) and left Egery 28/06/2008      *)
(*                                    *)
(* Entree : n k                       *)
(* Sortie : nombre de matrices n*n    *)
(*          ayant k elements a 1 dans *)
(*          chaque ligne et chaque    *)
(*          colonne (les autres a 0)  *)
(*                                    *)

(* Module precision arbritraire *)
#load "nums.cma";;
open Num;;
let big = num_of_int ;; (* Alias fonction ! *)

(* Pour chaque partie (sous-matrice), on recense              *
 * sa largeur et le nombre de 1 qu'elle contient en hauteur.  *)
type partie = {taille: int ; cumul: int ; } 

(* Fonction factorielle *)
let fact n = 
    let rec loop n acc =
        assert (n>0) ;
        match n with
        | 0 | 1 -> acc
        | _     -> loop (n-1) acc */ big n
    in loop n (big 1) ;;

(* Donne le nombre de permutations pour une partition donnee *)
let evalue partition =
    let rec loop partition i acc =
        match partition with
        | []         ->  acc */ fact i
        | hd::tl     ->  loop tl (i + hd.taille) (acc // fact hd.taille) 
    in loop partition 0 (big 1) ;;

(* Fontion de partitionnement, selon nombre jetons. Renvoie :  *)
(* sous-partition, iterateur, jetons non places, cas ok, flag iteration *)
let insert n k partie  jetons  restants  ligne =
    let nb1      = min jetons partie.taille in assert (nb1 > 0) ;
    let nb0      = partie.taille - nb1      in
    match nb0, partie.cumul with
    (* Cas 1 : Tous les 1 sont deja places *)
    (* Cas 2 : Tous les 0 sont deja places, nb jetons suffisant *)
    (* Cas 3 ; Tous les 0 sont deja places, nb jetons insuffisant *)
    (* Cas 4 : Aucun 0 injecte *)
    (* Cas 5 : Des 1 et des 0 injectes : subdivision en 2 matrices bloc *)
    | _, c when c = k           ->    ([],      0, restants,    false, true)
    | _, c when (ligne-c = n-k) && (jetons >= partie.taille)
        -> ( [{partie with cumul = c+1}],   nb1-1, restants-nb1, true, true)
    | _, c when (ligne-c = n-k) ->    ([],      0, 0,          false, false)
    | 0, c -> ( [{partie with cumul=c+1}],  nb1-1, restants-nb1, true, true)
    | _, c -> ([{taille=nb1 ; cumul=c+1} ;
                {taille=nb0 ; cumul=c  }],  nb1-1, restants-nb1, true, true)

(* Notre fonction principale, qui amorce le decompte *)
let denombre_matrices n k =

   (* On construit les matrices convenables,              *
    * Et on decompte pour chaque le nombre de combinaison *)
   let rec loop prefixe suffixe jetons_hd jetons_restant ligne acc =
        match suffixe, jetons_hd, jetons_restant, ligne with  
        (* Cas 1 : matrice construite, on decompte *)
        (* Cas 2 : Jetons tous places, on passe a ligne suivante *)
        (* Cas 3 : Jetons restent a inserer, mais fin de ligne atteinte *)
        (* Cas 4 : Jetons places dans partie suivante *)
        (* Cas 5 : Sous partitionne selon jetons inseres, et itere cas suivant *)
        | _ ,     _,  _,  l when l = n -> acc +/ evalue ( prefixe @ suffixe )  
        | _ ,     _,  0,  _ -> loop [] (prefixe @ suffixe) k  k  (ligne+1) acc
        | [],     _,  _,  _ -> acc   
        | hd::tl, 0,  jr, _ -> loop (prefixe @ [hd])   tl  jr jr  ligne    acc
        | hd::tl, jh, jr, _
        ->  let (hd', jh', jr', ok, itere) = insert n k hd jh jr ligne in
        let acc2 = if ok then loop (prefixe@hd') tl jr' jr' ligne acc  else acc
            in  if itere then loop prefixe  suffixe jh' jr  ligne acc2 else acc2
                
   in loop [] [{taille = n ; cumul = 0}] k k 0 (big 0) ;;

let test n k = string_of_num (denombre_matrices n k) ;;

Le résultat est lent :

- un des tests d'arrêt n'est pas effectué aussi tôt qu'il le pourrait (je laisse en guise d'exercice la détection de ce soucis !). Une autre version de ce programme corrige ce point, sans résoudre le problème de lenteur.
- il s'agit de mon premier programme en Caml. Je mets peut-être en défaut la machine OCaml (pas de récursion terminale pour le parcours de l'arbre des solutions, par exemple ?).

De toute façon, l'idéal serait de trouver les éventuelles relations de récurrence pour k>2.

Yves

[KindPlayer]

J'ai fait une version scheme a base de liste, mais le probleme c'est que je ne sais pas trop a quels endroit je pourrais eviter certaines enumérations. Il faut analyser plus finement le pb

[bredelet]

Bien, voila qui a l'air interessant

Quelqu'un a parle d'une solution 'fractale'

Je vais commencer par quelques observations:

- On peut inverser les 0 et les 1, ce qui fait que defi4(n, k) = defi4(n, n-k)
- Ce faisant il suffit de s'interesser aux grilles qui ont au moins autant de 0 que de 1
- Une solution triviale de (n, k) se compose d'une solution de (m, k) et une solution de (p, k) avec m+p = n et m>=k et p>=k. On remplit le reste de zeros. Par exemple (5, 2):

Code :

1
2
3
4
5
6
[ 1  1  0  0  0 ]          [ a  a  0  0  0 ]
[ 1  1  0  0  0 ]          [ a  a  0  0  0 ]
[ 0  0  1  1  0 ]   soit   [ 0  0  b  b  b ]
[ 0  0  1  0  1 ]          [ 0  0  b  b  b ]
[ 0  0  0  1  1 ]          [ 0  0  b  b  b ]

La on voit que [a] est une solution de (2, 2) et [b] est une solution de (3, 2)

- On peut echanger un 1 de [a] avec un 0 dans la zone des zeros, pour peu que l'on fasse de meme avec [b]. Par exemple:

Code :

1
2
3
4
5
6
[ 0  1  1  0  0 ]          [ 0  a  1  0  0 ]
[ 1  1  0  0  0 ]          [ a  a  0  0  0 ]
[ 1  0  0  1  0 ]   soit   [ 1  0  0  b  b ]
[ 0  0  1  0  1 ]          [ 0  0  b  b  b ]
[ 0  0  0  1  1 ]          [ 0  0  b  b  b ]

Je *pense* qu'on peut arriver a une definition recursive de defi4(n, k) sur ces bases, ou au moins profiter de ces proprietes pour reduire le nombre de calculs. Je n'y suis pas encore.

[KindPlayer]

Ouai il doit y avoir des trucs de combinatoire comme ça qui évitent de tout énumérer. L'ennui avec les "trucs" c'est qu'on a vite fait d'oublier des solutions en route même si on gagne beaucoup en vitesse. Le seul truc pour lequel j'étais vraiment sur, c'est qu'il y a C_n^k facons de remplir la premiere ligne par exemple. Mais la combinatoire c'est moyennement mon fort, même si c'est un domaine que je trouve parmi les plus interessants.