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 154 : matpermanent(GEN M)
27 : {
28 : pari_sp av;
29 154 : long n = lg(M)-1, i, x, upper;
30 : GEN p, in;
31 154 : if (typ(M) != t_MAT) pari_err_TYPE("matpermanent", M);
32 154 : if (!n) return gen_1;
33 147 : if (n != nbrows(M)) pari_err_DIM("matpermanent");
34 : #ifdef LONG_IS_64BIT /* because of vals(long x) => x <= LONG_MAX */
35 120 : 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 133 : if (n == 1) return gcopy(gcoeff(M,1,1));
40 :
41 112 : av = avma;
42 112 : if (RgM_is_QM(M))
43 : {
44 : GEN cM;
45 98 : M = Q_primitive_part(M, &cM);
46 98 : p = ZM_permanent(M);
47 98 : if (cM) p = gerepileupto(av, gmul(p, gpowgs(cM,n)));
48 98 : return p;
49 : }
50 :
51 14 : p = gen_0;
52 14 : in = zerovec(n);
53 14 : upper = 1L << n;
54 112 : for (x = 1; x < upper; x++)
55 : {
56 98 : long gray = x ^ (x>>1), k = vals(x);
57 98 : GEN col = gel(M,k+1);
58 98 : if (gray & (1L<<k))
59 224 : { for (i=1; i <= n; ++i) gel(in, i) = gadd(gel(in, i), gel(col, i)); }
60 : else
61 168 : { for (i=1; i <= n; ++i) gel(in, i) = gsub(gel(in, i), gel(col, i)); }
62 98 : if (hammingl(gray)&1)
63 56 : p = gsub(p, RgV_prod(in));
64 : else
65 42 : p = gadd(p, RgV_prod(in));
66 98 : if (gc_needed(av, 1)) gerepileall(av, 2, &in, &p);
67 : }
68 14 : if (n&1) p = gneg(p);
69 14 : return gerepileupto(av, p);
70 : }
71 :
72 : /* ||M||_oo = max_i \sum_j | M[i,j] | */
73 : static GEN
74 98 : ZM_normoo(GEN M)
75 : {
76 98 : long i, j, m, l = lg(M);
77 98 : GEN N = NULL;
78 98 : if (l == 1) return gen_0;
79 98 : m = lgcols(M);
80 938 : for (i = 1; i < m; i++)
81 : {
82 840 : GEN s = mpabs_shallow(gcoeff(M,i,1));
83 11214 : for (j = 2; j < l; j++) s = addii(s, mpabs_shallow(gcoeff(M,i,j)));
84 840 : if (!N || abscmpii(N, s) < 0) N = s;
85 : }
86 98 : return N;
87 : }
88 :
89 : /* Assumes M square; dimensions <= 31x31 (32-bit) or 63x63 (64-bit). */
90 : GEN
91 98 : ZM_permanent(GEN M)
92 : {
93 98 : pari_sp av = avma;
94 : GEN p, in;
95 98 : long i, x, upper, n = lg(M)-1;
96 98 : if (!is_bigint(ZM_normoo(M)))
97 91 : return gerepileuptoint(av, zm_permanent(ZM_to_zm(M)));
98 7 : p = gen_0;
99 7 : in = zerocol(n);
100 7 : upper = (1L<<n);
101 7340032 : for (x = 1; x < upper; x++)
102 : {
103 7340025 : long gray = x ^ (x>>1), k = vals(x);
104 7340025 : GEN c, col = gel(M, k+1);
105 7340025 : if (gray & (1L<<k))
106 77070336 : { for (i=1; i <= n; ++i) gel(in, i) = addii(gel(in,i), gel(col,i)); }
107 : else
108 77070189 : { for (i=1; i <= n; ++i) gel(in, i) = subii(gel(in,i), gel(col,i)); }
109 7340025 : c = ZV_prod(in);
110 7340025 : if (hammingl(gray)&1) togglesign_safe(&c);
111 7340025 : p = addii(p, c);
112 7340025 : if (gc_needed(av, 1)) gerepileall(av, 2, &in, &p);
113 : }
114 7 : if (n&1) togglesign_safe(&p);
115 7 : return gerepilecopy(av, p);
116 : }
117 :
118 : /* Assumes M square; dimensions <= 31x31 (32-bit) or 63x63 (64-bit). */
119 : GEN
120 91 : zm_permanent(GEN M)
121 : {
122 91 : pari_sp av = avma;
123 91 : long n = lg(M)-1;
124 91 : ulong x, upper = (1UL<<n);
125 91 : GEN p = gen_0, in = const_vecsmall(n, 0);
126 91 : pari_sp av2 = avma;
127 14694484 : for (x = 1; x < upper; x++)
128 : {
129 14694393 : ulong gray = x ^ (x>>1);
130 14694393 : long i, k = vals(x);
131 14694393 : GEN c, col = gel(M, k+1);
132 14694393 : if (gray & (1UL<<k))
133 154212562 : { for (i = 1; i <= n; ++i) in[i] += col[i]; }
134 : else
135 154211771 : { for (i = 1; i <= n; ++i) in[i] -= col[i]; }
136 14694393 : c = vecsmall_prod(in);
137 14694393 : if (hammingl(gray)&1) togglesign_safe(&c);
138 14694393 : p = addii(p, c);
139 14694393 : if (gc_needed(av2, 1)) p = gerepileuptoint(av2, p);
140 : }
141 91 : if (n&1) togglesign_safe(&p);
142 91 : return gerepileuptoint(av, p);
143 : }
|
