Line data Source code
1 : /* Copyright (C) 2016 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /********************************************************************/
16 : /** **/
17 : /** MATRIX PERMANENT, via RYSER'S FORMULA **/
18 : /** (initial implementation: C. Greathouse) **/
19 : /** **/
20 : /********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : /* Ryser's formula */
25 : GEN
26 132 : matpermanent(GEN M)
27 : {
28 : pari_sp av;
29 132 : long n = lg(M)-1, i, x, upper;
30 : GEN p, in;
31 132 : if (typ(M) != t_MAT) pari_err_TYPE("matpermanent", M);
32 132 : if (!n) return gen_1;
33 126 : if (n != nbrows(M)) pari_err_DIM("matpermanent");
34 : #ifdef LONG_IS_64BIT /* because of vals(long x) => x <= LONG_MAX */
35 100 : if (n > 63) pari_err_IMPL("large matrix permanent");
36 : #else
37 20 : if (n > 31) pari_err_IMPL("large matrix permanent");
38 : #endif
39 114 : if (n == 1) return gcopy(gcoeff(M,1,1));
40 :
41 96 : av = avma;
42 96 : if (RgM_is_QM(M))
43 : {
44 : GEN cM;
45 84 : M = Q_primitive_part(M, &cM);
46 84 : p = ZM_permanent(M);
47 84 : if (cM) p = gc_upto(av, gmul(p, gpowgs(cM,n)));
48 84 : return p;
49 : }
50 :
51 12 : p = gen_0;
52 12 : in = zerovec(n);
53 12 : upper = 1L << n;
54 96 : for (x = 1; x < upper; x++)
55 : {
56 84 : long gray = x ^ (x>>1), k = vals(x);
57 84 : GEN col = gel(M,k+1);
58 84 : if (gray & (1L<<k))
59 192 : { for (i=1; i <= n; ++i) gel(in, i) = gadd(gel(in, i), gel(col, i)); }
60 : else
61 144 : { for (i=1; i <= n; ++i) gel(in, i) = gsub(gel(in, i), gel(col, i)); }
62 84 : if (hammingu(gray)&1)
63 48 : p = gsub(p, RgV_prod(in));
64 : else
65 36 : p = gadd(p, RgV_prod(in));
66 84 : if (gc_needed(av, 1)) (void)gc_all(av, 2, &in, &p);
67 : }
68 12 : if (n&1) p = gneg(p);
69 12 : return gc_upto(av, p);
70 : }
71 :
72 : /* ||M||_oo = max_i \sum_j | M[i,j] | */
73 : static GEN
74 84 : ZM_normoo(GEN M)
75 : {
76 84 : long i, j, m, l = lg(M);
77 84 : GEN N = NULL;
78 84 : if (l == 1) return gen_0;
79 84 : m = lgcols(M);
80 804 : for (i = 1; i < m; i++)
81 : {
82 720 : GEN s = mpabs_shallow(gcoeff(M,i,1));
83 9612 : for (j = 2; j < l; j++) s = addii(s, mpabs_shallow(gcoeff(M,i,j)));
84 720 : if (!N || abscmpii(N, s) < 0) N = s;
85 : }
86 84 : return N;
87 : }
88 :
89 : /* Assumes M square; dimensions <= 31x31 (32-bit) or 63x63 (64-bit). */
90 : GEN
91 84 : ZM_permanent(GEN M)
92 : {
93 84 : pari_sp av = avma;
94 : GEN p, in;
95 84 : long i, x, upper, n = lg(M)-1;
96 84 : if (!is_bigint(ZM_normoo(M)))
97 78 : return gc_INT(av, zm_permanent(ZM_to_zm(M)));
98 6 : p = gen_0;
99 6 : in = zerocol(n);
100 6 : upper = (1L<<n);
101 6291456 : for (x = 1; x < upper; x++)
102 : {
103 6291450 : long gray = x ^ (x>>1), k = vals(x);
104 6291450 : GEN c, col = gel(M, k+1);
105 6291450 : if (gray & (1L<<k))
106 66060288 : { for (i=1; i <= n; ++i) gel(in, i) = addii(gel(in,i), gel(col,i)); }
107 : else
108 66060162 : { for (i=1; i <= n; ++i) gel(in, i) = subii(gel(in,i), gel(col,i)); }
109 6291450 : c = ZV_prod(in);
110 6291450 : if (hammingu(gray)&1) togglesign_safe(&c);
111 6291450 : p = addii(p, c);
112 6291450 : if (gc_needed(av, 1)) (void)gc_all(av, 2, &in, &p);
113 : }
114 6 : if (n&1) togglesign_safe(&p);
115 6 : return gc_GEN(av, p);
116 : }
117 :
118 : /* Assumes M square; dimensions <= 31x31 (32-bit) or 63x63 (64-bit). */
119 : GEN
120 78 : zm_permanent(GEN M)
121 : {
122 78 : pari_sp av = avma;
123 78 : long n = lg(M)-1;
124 78 : ulong x, upper = (1UL<<n);
125 78 : GEN p = gen_0, in = const_vecsmall(n, 0);
126 78 : pari_sp av2 = avma;
127 12595272 : for (x = 1; x < upper; x++)
128 : {
129 12595194 : ulong gray = x ^ (x>>1);
130 12595194 : long i, k = vals(x);
131 12595194 : GEN c, col = gel(M, k+1);
132 12595194 : if (gray & (1UL<<k))
133 132182196 : { for (i = 1; i <= n; ++i) in[i] += col[i]; }
134 : else
135 132181518 : { for (i = 1; i <= n; ++i) in[i] -= col[i]; }
136 12595194 : c = vecsmall_prod(in);
137 12595194 : if (hammingu(gray)&1) togglesign_safe(&c);
138 12595194 : p = addii(p, c);
139 12595194 : if (gc_needed(av2, 1)) p = gc_INT(av2, p);
140 : }
141 78 : if (n&1) togglesign_safe(&p);
142 78 : return gc_INT(av, p);
143 : }
|
