File: poly.ml

package info (click to toggle)
sks 1.1.6-14
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 2,296 kB
  • sloc: ml: 15,228; ansic: 1,069; sh: 358; makefile: 347
file content (228 lines) | stat: -rw-r--r-- 6,338 bytes parent folder | download | duplicates (5)
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
(***********************************************************************)
(* poly.ml - Simple polynomial implementation                          *)
(*                                                                     *)
(* Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, *)
(*               2011, 2012, 2013  Yaron Minsky and Contributors       *)
(*                                                                     *)
(* This file is part of SKS.  SKS is free software; you can            *)
(* redistribute it and/or modify it under the terms of the GNU General *)
(* Public License as published by the Free Software Foundation; either *)
(* version 2 of the License, or (at your option) any later version.    *)
(*                                                                     *)
(* This program is distributed in the hope that it will be useful, but *)
(* WITHOUT ANY WARRANTY; without even the implied warranty of          *)
(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU   *)
(* General Public License for more details.                            *)
(*                                                                     *)
(* You should have received a copy of the GNU General Public License   *)
(* along with this program; if not, write to the Free Software         *)
(* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 *)
(* USA or see <http://www.gnu.org/licenses/>.                          *)
(***********************************************************************)

open StdLabels
open MoreLabels
module Unix = UnixLabels
open Printf
open Scanf
open ZZp.Infix
module Map = PMap.Map

let rec rfind ~f low high =
  if low >= high then raise Not_found
  else if f(low) then low
  else rfind ~f (low + 1) high

type t = { a : ZZp.zz array;
           (** coefficients, listed from lowest to highest degree *)
           degree : int; (** degree of polynomial *)
         }

let compute_degree a =
  let rec loop a i =
    if i <= 0 then 0
    else (
      if a.(i) =: ZZp.zero
      then loop a (i - 1)
      else i
    )
  in
  loop a (Array.length a - 1)

let init degree ~f =
  let a = Array.init (degree + 1) ~f:(fun i -> f i) in
  let degree = compute_degree a in
  { a = (if degree + 1 < Array.length a
         then Array.sub a ~pos:0 ~len:(degree + 1)
         else a);
    degree = degree;
  }

let make degree x =
  if x =: ZZp.zero then { a = [| ZZp.zero |]; degree = 0; }
  else
    { a = Array.init (degree + 1) ~f:(fun i -> x);
      degree = degree;
    }

let zero = make 0 ZZp.zero
let one = make 0 ZZp.one

(* Get and set coeffs *)
(*let getc x i =  x.a.(i)
  let setc x i v = x.a.(i) <- v
  let lgetc x i = x.a.(i)
  let rgetc x i = x.a.(i) *)
let degree x = x.degree
let length x = Array.length x.a

let copy x = { x with a = Array.copy x.a }

let to_string x =
  let buf = Buffer.create 0 in
  for i = degree x downto 1 do
    bprintf buf "%s z^%d + " (ZZp.to_string x.a.(i)) i;
  done;
  if degree x >= 0
  then bprintf buf "%s" (ZZp.to_string x.a.(0))
  else bprintf buf "0";
  Buffer.contents buf

let splitter = Str.regexp "[ \t]+\\+[ \t]+"

let parse_digit s =
  try sscanf s "%s z^%d" (fun digit degree -> (degree,ZZp.of_string digit))
  with End_of_file -> (0,ZZp.of_string s)

let map_keys map =
  Map.fold ~init:[] ~f:(fun ~key ~data keylist -> key::keylist) map


let of_string s =
  let digits = List.map ~f:parse_digit (Str.split splitter s) in
  let digitmap = Map.of_alist digits in
  let degree = MList.reduce ~f:max (map_keys digitmap) in
  init degree ~f:(fun deg ->
                    try Map.find deg digitmap
                    with Not_found -> ZZp.zero)



let print x =
  for i = degree x downto 1 do
    ZZp.print x.a.(i);
    printf " z^%d + " i;
  done;
  if degree x >= 0 then
    ZZp.print x.a.(0)
  else
    print_string "0"

exception NotEqual

let eq x y =
  try
    if x.degree <> y.degree then raise NotEqual;
    for i = 0 to x.degree do
      if x.a.(i) <>: y.a.(i)
      then raise NotEqual
    done;
    true
  with
      NotEqual -> false


let of_array array =
  if Array.length array = 0 then zero
  else
    let deg = compute_degree array in
    { a = Array.init (deg + 1) ~f:(fun i -> array.(i));
      degree = deg;
    }

let term deg c =
  init ~f:(fun i -> if i = deg then c else ZZp.zero) deg

let set_length length x =
  assert (length + 1 > degree x);
  { a = Array.init (length + 1)
            ~f:(fun i ->
                  if i <= x.degree
                  then x.a.(i)
                  else ZZp.zero);
    degree = x.degree
  }

let to_array x = Array.copy x.a
let is_monic x = x.a.(degree x) =: ZZp.one

let eval poly z =
  let zd = ref ZZp.one
  and sum = ref ZZp.zero in
  for deg = 0 to degree poly do
    sum := !sum +: poly.a.(deg) *: !zd;
    zd := !zd *: z
  done;
  !sum

let mult x y =
  let mdegree = degree x + degree y in
  let prod = { a = Array.make ( mdegree + 1 ) ZZp.zero;
               degree = mdegree ;
             }
  in
  for i = 0 to degree x  do
    for j = 0 to degree y do
      prod.a.(i + j) <- prod.a.(i + j) +: x.a.(i) *: y.a.(j)
    done
  done;
  prod

(** scalar multiplication *)
let scmult x c =
  { x with a = Array.map ~f:(fun z -> z *: c) x.a; }

let add x y =
  let deg = max x.degree y.degree in
  init deg
    ~f:(fun i ->
          (if i <= x.degree then x.a.(i) else ZZp.zero) +:
          (if i <= y.degree then y.a.(i) else ZZp.zero))

let neg x = { x with a = Array.map ~f:(fun c -> ZZp.neg c) x.a }

let sub x y = add x (neg y)

let rec divmod x y =
  if eq x zero then (zero,zero)
  else if degree y > degree x then (zero,x)
  else
    let degdiff = degree x - degree y in
    assert (degdiff >= 0);
    let c = x.a.(degree x) /: y.a.(degree y) in
    let m = term degdiff c in
    let new_x = sub x (mult m y) in
    assert (degree new_x < degree x || degree x = 0);
    let (q,r) = divmod new_x y in
    (add q m,r)

let modulo x y = let (q,r) = divmod x y in r
let div x y = let (q,r) = divmod x y in q

let const_coeff x = x.a.(0)
let nth_coeff x n = x.a.(n)
let const c = make 0 c


let rec gcd_rec x y =
  if eq y zero then x
  else
    let (q,r) = divmod x y in
    gcd_rec y r

let gcd x y =
  let result = gcd_rec x y in
  (* force the GCD to be monic *)
  mult result (const (ZZp.inv result.a.(degree result)))