← Back to branch summary

~ubuntu-branches/debian/experimental/sks/experimental

« back to all changes in this revision

Viewing changes to linearAlg.ml

  • Committer: Package Import Robot
  • Author(s): Daniel Kahn Gillmor
  • Date: 2013-06-27 16:39:02 UTC
  • mfrom: (1.1.4)
  • Revision ID: package-import@ubuntu.com-20130627163902-qqic4va2187boeji
Tags: 1.1.4-1
http://bugs.debian.org/690135
http://bugs.debian.org/709322
* New Upstream Release (Closes: #690135)
* added myself to Uploaders.
* convert to dh 9
* Standards-Version: bump to 3.9.4 (no changes needed)
* debian/rules: clean up
* refresh and clean up debian/patches
* switch packaging vcs to git
* avoid trying to upgrade DB_CONFIG (Closes: #709322)

Show diffs side-by-side

added added

removed removed

Lines of Context:
linearAlg.ml
1
 
(************************************************************************)
2
 
(* This file is part of SKS.  SKS is free software; you can
3
 
   redistribute it and/or modify it under the terms of the GNU General
4
 
   Public License as published by the Free Software Foundation; either
5
 
   version 2 of the License, or (at your option) any later version.
6
 
 
7
 
   This program is distributed in the hope that it will be useful, but
8
 
   WITHOUT ANY WARRANTY; without even the implied warranty of
9
 
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
10
 
   General Public License for more details.
11
 
 
12
 
   You should have received a copy of the GNU General Public License
13
 
   along with this program; if not, write to the Free Software
14
 
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
15
 
   USA *)
 
1
(***********************************************************************)
 
2
(* linearAlg.ml                                                        *)
 
3
(*                                                                     *)
 
4
(* Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, *)
 
5
(*               2011, 2012  Yaron Minsky and Contributors             *)
 
6
(*                                                                     *)
 
7
(* This file is part of SKS.  SKS is free software; you can            *)
 
8
(* redistribute it and/or modify it under the terms of the GNU General *)
 
9
(* Public License as published by the Free Software Foundation; either *)
 
10
(* version 2 of the License, or (at your option) any later version.    *)
 
11
(*                                                                     *)
 
12
(* This program is distributed in the hope that it will be useful, but *)
 
13
(* WITHOUT ANY WARRANTY; without even the implied warranty of          *)
 
14
(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU   *)
 
15
(* General Public License for more details.                            *)
 
16
(*                                                                     *)
 
17
(* You should have received a copy of the GNU General Public License   *)
 
18
(* along with this program; if not, write to the Free Software         *)
 
19
(* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 *)
 
20
(* USA or see <http://www.gnu.org/licenses/>.                          *)
16
21
(***********************************************************************)
17
22
 
18
23
open StdLabels
31
36
    riter ~f (low + 1) high
32
37
  )
33
38
 
34
 
let rec rfind ~f low high = 
 
39
let rec rfind ~f low high =
35
40
  if low >= high then raise Not_found
36
41
  else if f(low) then low
37
42
  else rfind ~f (low + 1) high
41
46
(*********************************************************************)
42
47
(*********************************************************************)
43
48
(*********************************************************************)
44
 
    
45
 
module MatrixSlow = 
 
49
 
 
50
module MatrixSlow =
46
51
struct
47
52
 
48
53
  type t = { columns: int;
49
 
             rows: int;
50
 
             array: ZZp.zz array;
51
 
           }
 
54
             rows: int;
 
55
             array: ZZp.zz array;
 
56
           }
52
57
 
53
58
  let columns m = m.columns
54
59
  let rows m = m.rows
55
60
  let dims t = (t.columns,t.rows)
56
 
                 
 
61
 
57
62
  let copy m = { m with array = Array.copy m.array; }
58
63
 
59
 
  let make ~columns ~rows init = 
 
64
  let make ~columns ~rows init =
60
65
    let array = Array.create (columns * rows) init in
61
66
    { columns = columns;
62
67
      rows = rows;
63
68
      array = array;
64
69
    }
65
70
 
66
 
  let init ~columns ~rows ~f = 
 
71
  let init ~columns ~rows ~f =
67
72
    { columns = columns;
68
73
      rows = rows;
69
 
      array = 
70
 
        Array.init (columns * rows)
71
 
          ~f:(fun i -> 
72
 
                let (i,j) = i mod columns, i / columns in
73
 
                f i j)
 
74
      array =
 
75
        Array.init (columns * rows)
 
76
          ~f:(fun i ->
 
77
                let (i,j) = i mod columns, i / columns in
 
78
                f i j)
74
79
    }
75
80
 
76
 
  let get m i j = 
 
81
  let get m i j =
77
82
    m.array.(i + j * m.columns)
78
83
 
79
 
  let set m i j v = 
 
84
  let set m i j v =
80
85
    m.array.(i + j * m.columns) <- v
81
86
 
82
 
  let scmult_ip m sc = 
 
87
  let scmult_ip m sc =
83
88
    for i = 0 to Array.length m.array - 1 do
84
89
      m.array.(i) <- ZZp.mult m.array.(i) sc
85
90
    done
86
91
 
87
 
  let scmult m v = 
 
92
  let scmult m v =
88
93
    { m with
89
 
        array = Array.map ~f:(fun x -> ZZp.mult v x) m.array 
 
94
        array = Array.map ~f:(fun x -> ZZp.mult v x) m.array
90
95
    }
91
96
 
92
97
  let scmult_row m j sc =
95
100
      m.array.(start + i) <- ZZp.mult m.array.(start + i) sc
96
101
    done
97
102
 
98
 
  let swap_rows m j1 j2 = 
99
 
    let start1 = j1 * m.columns 
 
103
  let swap_rows m j1 j2 =
 
104
    let start1 = j1 * m.columns
100
105
    and start2 = j2 * m.columns in
101
 
    riter 0 m.columns 
102
 
      ~f:(fun i -> 
103
 
            let tmp = m.array.(start1 + i) in
104
 
            m.array.(start1 + i) <- m.array.(start2 + i);
105
 
            m.array.(start2 + i) <- tmp)
 
106
    riter 0 m.columns
 
107
      ~f:(fun i ->
 
108
            let tmp = m.array.(start1 + i) in
 
109
            m.array.(start1 + i) <- m.array.(start2 + i);
 
110
            m.array.(start2 + i) <- tmp)
106
111
 
107
 
  let add_ip m1 m2 = 
 
112
  let add_ip m1 m2 =
108
113
    if m1.columns <> m2.columns || m1.rows <> m2.rows then
109
114
      raise LayoutMismatch;
110
115
    for i = 0 to Array.length m1.array - 1 do
111
116
      m1.array.(i) <- ZZp.add m1.array.(i) m2.array.(i)
112
117
    done
113
118
 
114
 
  let add m1 m2 = 
 
119
  let add m1 m2 =
115
120
    if m1.columns <> m2.columns || m1.rows <> m2.rows then
116
121
      raise LayoutMismatch;
117
122
    { m1 with
118
 
        array = Array.init (m1.columns * m1.rows)
119
 
                  ~f:(fun i -> ZZp.add m1.array.(i) m2.array.(i))
 
123
        array = Array.init (m1.columns * m1.rows)
 
124
                  ~f:(fun i -> ZZp.add m1.array.(i) m2.array.(i))
120
125
    }
121
126
 
122
 
  let rec idot_rec m1 m2 ~i ~pos1 ~pos2 sum = 
 
127
  let rec idot_rec m1 m2 ~i ~pos1 ~pos2 sum =
123
128
    if i >= m1.columns then sum
124
 
    else 
 
129
    else
125
130
      idot_rec m1 m2 ~i:(i+1) ~pos1:(pos1 + 1) ~pos2:(pos2 + m2.columns)
126
 
        (ZZp.add sum (ZZp.mult m1.array.(pos1) m2.array.(pos2)))
 
131
        (ZZp.add sum (ZZp.mult m1.array.(pos1) m2.array.(pos2)))
127
132
 
128
 
  let idot m1 m2 i j = 
 
133
  let idot m1 m2 i j =
129
134
    idot_rec m1 m2 ~i:0 ~pos1:(m1.columns * i) ~pos2:j ZZp.zero
130
135
 
131
 
  let mult m1 m2  = 
 
136
  let mult m1 m2  =
132
137
    if m1.columns <> m2.rows then
133
138
      raise LayoutMismatch;
134
139
    init ~columns:m2.columns ~rows:m1.rows
135
140
      ~f:(fun i j -> idot m1 m2 i j)
136
141
 
137
142
 
138
 
  let transpose m = 
 
143
  let transpose m =
139
144
    init ~columns:m.rows ~rows:m.columns ~f:(fun i j -> get m j i)
140
145
 
141
146
 
142
 
  let rowadd m ~src ~dst ~scmult = 
 
147
  let rowadd m ~src ~dst ~scmult =
143
148
    for i = 0 to m.columns - 1 do
144
149
      let newval = ZZp.add (ZZp.mult (get m i src) scmult) (get m i dst) in
145
150
      set m i dst newval
146
151
    done
147
152
 
148
 
  let rowsub m ~src ~dst ~scmult = 
 
153
  let rowsub m ~src ~dst ~scmult =
149
154
    if scmult <>: ZZp.one then
150
155
      for i = 0 to m.columns - 1 do
151
 
        let sval = get m i src in
152
 
        if sval <>: ZZp.zero then
153
 
          let newval = ZZp.sub (get m i dst) (ZZp.mult_fast sval scmult) in
154
 
          set m i dst newval
 
156
        let sval = get m i src in
 
157
        if sval <>: ZZp.zero then
 
158
          let newval = ZZp.sub (get m i dst) (ZZp.mult_fast sval scmult) in
 
159
          set m i dst newval
155
160
      done
156
161
    else
157
162
      for i = 0 to m.columns - 1 do
158
 
        let sval = get m i src in
159
 
        if sval <>: ZZp.zero then
160
 
          let newval = ZZp.sub (get m i dst) sval in
161
 
          set m i dst newval
 
163
        let sval = get m i src in
 
164
        if sval <>: ZZp.zero then
 
165
          let newval = ZZp.sub (get m i dst) sval in
 
166
          set m i dst newval
162
167
      done
163
168
 
164
169
  let print m =
165
170
    for j = 0 to m.rows - 1 do
166
171
      print_string "| ";
167
172
      for i = 0 to m.columns - 1 do
168
 
        ZZp.print (get m i j);
169
 
        print_string " "
 
173
        ZZp.print (get m i j);
 
174
        print_string " "
170
175
      done;
171
176
      print_string " |\n"
172
177
    done
178
183
(*********************************************************************************)
179
184
 
180
185
(* Does everything in-place, using the in-place numerix operators *)
181
 
module Matrix = 
 
186
module Matrix =
182
187
struct
183
188
 
184
189
  type t = { columns: int;
185
 
             rows: int;
186
 
             array: ZZp.zzref array;
187
 
           }
 
190
             rows: int;
 
191
             array: ZZp.zzref array;
 
192
           }
188
193
 
189
194
  let columns m = m.columns
190
195
  let rows m = m.rows
191
196
  let dims t = (t.columns,t.rows)
192
 
                 
 
197
 
193
198
  let copy m = { m with array = Array.copy m.array; }
194
199
 
195
 
  let init ~columns ~rows ~f = 
 
200
  let init ~columns ~rows ~f =
196
201
    { columns = columns;
197
202
      rows = rows;
198
 
      array = 
199
 
        Array.init (columns * rows)
200
 
          ~f:(fun i -> 
201
 
                let (i,j) = i mod columns, i / columns in
202
 
                ZZp.make_ref (f i j))
 
203
      array =
 
204
        Array.init (columns * rows)
 
205
          ~f:(fun i ->
 
206
                let (i,j) = i mod columns, i / columns in
 
207
                ZZp.make_ref (f i j))
203
208
    }
204
209
 
205
 
  let make ~columns ~rows x = 
 
210
  let make ~columns ~rows x =
206
211
    init ~columns ~rows ~f:(fun i j -> x)
207
212
 
208
 
  let lget m i j = 
 
213
  let lget m i j =
209
214
    ZZp.look (m.array.(i + j * m.columns))
210
215
 
211
 
  let rget m i j = 
 
216
  let rget m i j =
212
217
    m.array.(i + j * m.columns)
213
218
 
214
219
  let get m i j = ZZp.copy_out m.array.(i + j * m.columns)
215
220
 
216
 
  let set m i j v = 
 
221
  let set m i j v =
217
222
    ZZp.copy_in m.array.(i + j * m.columns) v
218
223
 
219
224
  let scmult_row ?(scol=0) m j sc =
223
228
      ZZp.mult_in v (ZZp.look v) sc
224
229
    done
225
230
 
226
 
  let swap_rows m j1 j2 = 
227
 
    let start1 = j1 * m.columns 
 
231
  let swap_rows m j1 j2 =
 
232
    let start1 = j1 * m.columns
228
233
    and start2 = j2 * m.columns in
229
 
    riter 0 m.columns 
230
 
      ~f:(fun i -> 
231
 
            let tmp = ZZp.copy_out m.array.(start1 + i) in
232
 
            ZZp.copy_in m.array.(start1 + i) (ZZp.look m.array.(start2 + i));
233
 
            ZZp.copy_in m.array.(start2 + i) tmp)
 
234
    riter 0 m.columns
 
235
      ~f:(fun i ->
 
236
            let tmp = ZZp.copy_out m.array.(start1 + i) in
 
237
            ZZp.copy_in m.array.(start1 + i) (ZZp.look m.array.(start2 + i));
 
238
            ZZp.copy_in m.array.(start2 + i) tmp)
234
239
 
235
 
  let transpose m = 
 
240
  let transpose m =
236
241
    init ~columns:m.rows ~rows:m.columns ~f:(fun i j -> lget m j i)
237
242
 
238
 
  let rowsub ?(scol=0) m ~src ~dst ~scmult = 
 
243
  let rowsub ?(scol=0) m ~src ~dst ~scmult =
239
244
    if scmult <>: ZZp.one then
240
245
      for i = scol to m.columns - 1 do
241
 
        let sval = rget m i src in
242
 
        if ZZp.look sval <>: ZZp.zero then
243
 
          let v = rget m i dst in
244
 
          ZZp.sub_in v (ZZp.look v) (ZZp.mult_fast (ZZp.look sval) scmult)
 
246
        let sval = rget m i src in
 
247
        if ZZp.look sval <>: ZZp.zero then
 
248
          let v = rget m i dst in
 
249
          ZZp.sub_in v (ZZp.look v) (ZZp.mult_fast (ZZp.look sval) scmult)
245
250
      done
246
251
    else
247
252
      for i = scol to m.columns - 1 do
248
 
        let sval = rget m i src in
249
 
        if ZZp.look sval <>: ZZp.zero then
250
 
          let v = rget m i dst in
251
 
          ZZp.sub_in v (ZZp.look v) (ZZp.look sval)
 
253
        let sval = rget m i src in
 
254
        if ZZp.look sval <>: ZZp.zero then
 
255
          let v = rget m i dst in
 
256
          ZZp.sub_in v (ZZp.look v) (ZZp.look sval)
252
257
      done
253
258
 
254
259
  let print m =
255
260
    for j = 0 to m.rows - 1 do
256
261
      print_string "| ";
257
262
      for i = 0 to m.columns - 1 do
258
 
        ZZp.print (lget m i j);
259
 
        print_string " "
 
263
        ZZp.print (lget m i j);
 
264
        print_string " "
260
265
      done;
261
266
      print_string " |\n"
262
267
    done
271
276
(****** Gauss-Jordan Reduction *****************)
272
277
 
273
278
let process_row m j =
274
 
  try 
275
 
    let v = 
 
279
  try
 
280
    let v =
276
281
      let v = Matrix.rget m j j in
277
282
      if ZZp.look v <>: ZZp.zero then v
278
283
      else
279
 
        let jswap = 
280
 
          try
281
 
            rfind (j + 1) (Matrix.rows m) 
282
 
              ~f:(fun jswap -> Matrix.lget m j jswap <>: ZZp.zero)
283
 
          with Not_found -> raise Exit
284
 
        in
285
 
        Matrix.swap_rows m j jswap;
286
 
        Matrix.rget m j j 
 
284
        let jswap =
 
285
          try
 
286
            rfind (j + 1) (Matrix.rows m)
 
287
              ~f:(fun jswap -> Matrix.lget m j jswap <>: ZZp.zero)
 
288
          with Not_found -> raise Exit
 
289
        in
 
290
        Matrix.swap_rows m j jswap;
 
291
        Matrix.rget m j j
287
292
    in
288
293
    if ZZp.look v <>: ZZp.one then Matrix.scmult_row m j (ZZp.inv (ZZp.look v));
289
294
    for j2 = 0 to Matrix.rows m - 1 do
290
 
      if j2 <> j 
 
295
      if j2 <> j
291
296
      then Matrix.rowsub m ~src:j ~dst:j2 ~scmult:(Matrix.get m j j2)
292
297
    done
293
298
  with
304
309
(****** Gaussian Reduction *****************)
305
310
 
306
311
let process_row_forward m j =
307
 
  try 
308
 
    let v = 
 
312
  try
 
313
    let v =
309
314
      let v = Matrix.rget m j j in
310
315
      if ZZp.look v <>: ZZp.zero then v
311
316
      else
312
 
        let jswap = 
313
 
          try
314
 
            rfind (j + 1) (Matrix.rows m) 
315
 
              ~f:(fun jswap -> Matrix.lget m j jswap <>: ZZp.zero)
316
 
          with Not_found -> raise Exit
317
 
        in
318
 
        Matrix.swap_rows m j jswap;
319
 
        Matrix.rget m j j 
 
317
        let jswap =
 
318
          try
 
319
            rfind (j + 1) (Matrix.rows m)
 
320
              ~f:(fun jswap -> Matrix.lget m j jswap <>: ZZp.zero)
 
321
          with Not_found -> raise Exit
 
322
        in
 
323
        Matrix.swap_rows m j jswap;
 
324
        Matrix.rget m j j
320
325
    in
321
326
    if ZZp.look v <>: ZZp.one then Matrix.scmult_row ~scol:j m j (ZZp.inv (ZZp.look v));
322
327
    for j2 = j + 1 to Matrix.rows m - 1 do
325
330
  with
326
331
      Exit -> ()
327
332
 
328
 
let backsubstitute m j = 
329
 
  if Matrix.lget m j j =: ZZp.one 
 
333
let backsubstitute m j =
 
334
  if Matrix.lget m j j =: ZZp.one
330
335
  then (
331
336
    let last = Matrix.rows m - 1 in
332
337
    for j2 = j - 1 downto 0 do
333
338
      Matrix.rowsub ~scol:last m ~src:j ~dst:j2 ~scmult:(Matrix.get m j j2);
334
 
      Matrix.set m j j2 ZZp.zero  
335
 
    done 
 
339
      Matrix.set m j j2 ZZp.zero
 
340
    done
336
341
  )
337
342
 
338
343
let greduce m =
344
349
  for j = Matrix.rows m - 1 downto 1 do
345
350
    backsubstitute m j;
346
351
  done
347
 
  
 
352
 
348
353
 
349
354
let reduce = greduce