| 1 | (* Modified by sweeks@sweeks.com 2001-10-03 to go in the MLton benchmark suite. |
| 2 | * Hardwired in the u6 list of polynomials and added a loop. |
| 3 | *) |
| 4 | (* tyan.sml |
| 5 | * A Grobner Basis calculation for polynomials over F17 |
| 6 | * Adapted from the TIL benchmark suite by Allyn Dimock: |
| 7 | * update to SML '97, Standard Basis Library, comment out unreachable code |
| 8 | * Original code from Thomas Yan, who has given his permission for this |
| 9 | * code be used as a benchmarking code for SML compilers |
| 10 | * (e-mail message Tue, 10 Apr 2001 13:07:44 -0400 (EDT)) |
| 11 | * |
| 12 | * The data structure for the intermediate results is described in |
| 13 | * @article{Yan:1998:GBSP, |
| 14 | * author = {Yan, Thomas}, |
| 15 | * title = {The Geobucked Data Structure For Polynomials}, |
| 16 | * journal = {Journal Of Symbolic Computation}, |
| 17 | * volume = 23, |
| 18 | * number = 3, |
| 19 | * pages = {285 -- 293}, |
| 20 | * year = 1998, |
| 21 | * } |
| 22 | *) |
| 23 | |
| 24 | val print : string -> unit = print |
| 25 | type 'a array1 = 'a Array.array |
| 26 | val sub1 = Array.sub |
| 27 | val update1 = Array.update |
| 28 | val array1 = Array.array |
| 29 | val length1 = Array.length |
| 30 | val op && = fn (i1, i2) => (Word.toInt (Word.andb (Word.fromInt (i1), Word.fromInt (i2)))) |
| 31 | val op || = fn (i1, i2) => (Word.toInt (Word.orb (Word.fromInt (i1), Word.fromInt (i2)))) |
| 32 | val op << = fn (i1, i2) => (Word.toInt (Word.<< (Word.fromInt (i1), Word.fromInt (i2)))) |
| 33 | val op >> = fn (i1, i2) => (Word.toInt (Word.>> (Word.fromInt (i1), Word.fromInt (i2)))) |
| 34 | infix && || << >> |
| 35 | fun fold f l b = List.foldl f b l |
| 36 | fun revfold f l b = List.foldr f b l |
| 37 | |
| 38 | val input_line = TextIO.inputLine |
| 39 | val end_of_stream = TextIO.endOfStream |
| 40 | val open_in = TextIO.openIn |
| 41 | val close_in = TextIO.closeIn |
| 42 | |
| 43 | nonfix smlnj_mod |
| 44 | nonfix smlnj_div |
| 45 | val smlnj_mod = op mod |
| 46 | val smlnj_div = op div |
| 47 | infix 7 smlnj_mod |
| 48 | infix 7 smlnj_div |
| 49 | |
| 50 | exception Tabulate |
| 51 | fun tabulate (i,f) = |
| 52 | if i <= 0 then raise Tabulate else |
| 53 | let val a = array1(i,f 0) |
| 54 | fun tabify j = if j < i then (update1(a,j,f j); tabify (j+1)) else a |
| 55 | in |
| 56 | tabify 1 |
| 57 | end |
| 58 | |
| 59 | exception ArrayofList |
| 60 | fun arrayoflist (hd::tl) = |
| 61 | let val a = array1((length tl) + 1,hd) |
| 62 | fun al([],_) = a |
| 63 | | al(hd::tl,i) = (update1(a,i,hd); al(tl,i+1)) |
| 64 | in |
| 65 | al(tl,1) |
| 66 | end |
| 67 | | arrayoflist ([]) = raise ArrayofList |
| 68 | |
| 69 | |
| 70 | structure Util = struct |
| 71 | datatype relation = Less | Equal | Greater |
| 72 | |
| 73 | exception NotImplemented of string |
| 74 | exception Impossible of string (* flag "impossible" condition *) |
| 75 | exception Illegal of string (* flag function use violating precondition *) |
| 76 | |
| 77 | fun error exn msg = raise (exn msg) |
| 78 | fun notImplemented msg = error NotImplemented msg |
| 79 | fun impossible msg = error Impossible msg |
| 80 | fun illegal msg = error Illegal msg |
| 81 | |
| 82 | (* arr[i] := obj :: arr[i]; extend non-empty arr if necessary *) |
| 83 | fun insert (obj,i,arr) = let |
| 84 | val len = length1 arr |
| 85 | val res = if i<len then (update1(arr,i,obj::sub1(arr,i)); arr) |
| 86 | else let val arr' = array1(Int.max(i+1,len+len),[]) |
| 87 | fun copy ~1 = (update1(arr',i,[obj]); arr') |
| 88 | | copy j = (update1(arr',j,sub1(arr,j)); |
| 89 | copy(j-1)) |
| 90 | in copy(len-1) end |
| 91 | in res |
| 92 | end |
| 93 | |
| 94 | (* |
| 95 | fun arrayoflists [] = arrayoflist [] |
| 96 | | arrayoflists ([]::ls) = arrayoflists ls |
| 97 | | arrayoflists [l] = arrayoflist l |
| 98 | | arrayoflists (ls as (obj0::_)::_) = let |
| 99 | val a = array1(revfold (fn (l,n) => length l + n) ls 0,obj0) |
| 100 | fun ins (i,[]) = i | ins (i,x::l) = (update1(a,i,x); ins(i+1,l)) |
| 101 | fun insert (i,[]) = a | insert (i,l::ll) = insert(ins(i,l),ll) |
| 102 | in insert(0,ls) end |
| 103 | *) |
| 104 | |
| 105 | (* given compare and array a, return list of contents of a sorted in |
| 106 | * ascending order, with duplicates stripped out; which copy of a duplicate |
| 107 | * remains is random. NOTE that a is modified. |
| 108 | *) |
| 109 | fun stripSort compare = fn a => let |
| 110 | infix sub |
| 111 | |
| 112 | |
| 113 | val op sub = sub1 and update = update1 |
| 114 | fun swap (i,j) = let val ai = a sub i |
| 115 | in update(a,i,a sub j); update(a,j,ai) end |
| 116 | (* sort all a[k], 0<=i<=k<j<=length a *) |
| 117 | fun s (i,j,acc) = if i=j then acc else let |
| 118 | val pivot = a sub ((i+j) smlnj_div 2) |
| 119 | fun partition (lo,k,hi) = if k=hi then (lo,hi) else |
| 120 | case compare (pivot,a sub k) of |
| 121 | Less => (swap (lo,k); partition (lo+1,k+1,hi)) |
| 122 | | Equal => partition (lo,k+1,hi) |
| 123 | | Greater => (swap (k,hi-1); partition (lo,k,hi-1)) |
| 124 | val (lo,hi) = partition (i,i,j) |
| 125 | in s(i,lo,pivot::s(hi,j,acc)) end |
| 126 | val res = s(0,length1 a,[]) |
| 127 | |
| 128 | in |
| 129 | res |
| 130 | end |
| 131 | end |
| 132 | |
| 133 | structure F = struct |
| 134 | val p = 17 |
| 135 | |
| 136 | datatype field = F of int (* for (F n), always 0<=n<p *) |
| 137 | (* exception Div = Integer.Div *) |
| 138 | (* unused code unless P.show, commented out in earlier version, is used |
| 139 | fun show (F x) = print (Int.toString x) |
| 140 | *) |
| 141 | (* unused code |
| 142 | val char = p |
| 143 | *) |
| 144 | |
| 145 | (* unused code unless P.display is used |
| 146 | val zero = F 0 |
| 147 | *) |
| 148 | val one = F 1 |
| 149 | fun coerceInt n = F (n smlnj_mod p) |
| 150 | |
| 151 | fun add (F n,F m) = let val k = n+m in if k>=p then F(k-p) else F k end |
| 152 | fun subtract (F n,F m) = if n>=m then F(n-m) else F(n-m+p) |
| 153 | fun negate (F 0) = F 0 | negate (F n) = F(p-n) |
| 154 | fun multiply (F n,F m) = F ((n*m) smlnj_mod p) |
| 155 | fun reciprocal (F 0) = raise Div |
| 156 | | reciprocal (F n) = let |
| 157 | (* consider euclid gcd alg on (a,b) starting with a=p, b=n. |
| 158 | * if maintain a = a1 n + a2 p, b = b1 n + b2 p, a>b, |
| 159 | * then when 1 = a = a1 n + a2 p, have a1 = inverse of n mod p |
| 160 | * note that it is not necessary to keep a2, b2 around. |
| 161 | *) |
| 162 | fun gcd ((a,a1),(b,b1)) = |
| 163 | if b=1 then (* by continued fraction expansion, 0<|b1|<p *) |
| 164 | if b1<0 then F(p+b1) else F b1 |
| 165 | else let val q = a smlnj_div b |
| 166 | in gcd((b,b1),(a-q*b,a1-q*b1)) end |
| 167 | in gcd ((p,0),(n,1)) end |
| 168 | (* unused code |
| 169 | fun divide (n,m) = multiply (n, reciprocal m) |
| 170 | *) |
| 171 | |
| 172 | (* unused code unless power is used |
| 173 | val andb = op && |
| 174 | val rshift = op >> |
| 175 | *) |
| 176 | |
| 177 | (* unused code |
| 178 | fun power(n,k) = |
| 179 | if k<=3 then case k of |
| 180 | 0 => one |
| 181 | | 1 => n |
| 182 | | 2 => multiply(n,n) |
| 183 | | 3 => multiply(n,multiply(n,n)) |
| 184 | | _ => reciprocal (power (n,~k)) (* know k<0 *) |
| 185 | else if andb(k,1)=0 then power(multiply(n,n),rshift(k,1)) |
| 186 | else multiply(n,power(multiply(n,n),rshift(k,1))) |
| 187 | *) |
| 188 | |
| 189 | fun isZero (F n) = n=0 |
| 190 | (* unused codeunless P.display is used |
| 191 | fun equal (F n,F m) = n=m |
| 192 | |
| 193 | fun display (F n) = if n<=p smlnj_div 2 then Int.toString n |
| 194 | else "-" ^ Int.toString (p-n) |
| 195 | *) |
| 196 | end |
| 197 | |
| 198 | structure M = struct (* MONO *) |
| 199 | local |
| 200 | val andb = op && |
| 201 | infix sub << >> andb |
| 202 | (* val op << = Bits.lshift and op >> = Bits.rshift and op andb = Bits.andb |
| 203 | *) |
| 204 | in |
| 205 | |
| 206 | (* encode (var,pwr) as a long word: hi word is var, lo word is pwr |
| 207 | masks 0xffff for pwr, mask ~0x10000 for var, rshift 16 for var |
| 208 | note that encoded pairs u, v have same var if u>=v, u andb ~0x10000<v |
| 209 | *) |
| 210 | |
| 211 | datatype mono = M of int list |
| 212 | (* |
| 213 | fun show (M x) = (print "<"; app (fn i => (print (Int.toString i); print ",")) x; print">") |
| 214 | *) |
| 215 | exception DoesntDivide |
| 216 | |
| 217 | (* unused code |
| 218 | val numVars = 32 |
| 219 | *) |
| 220 | |
| 221 | val one = M [] |
| 222 | fun x_i v = M [(v<<16)+1] |
| 223 | fun explode (M l) = map (fn v => (v>>16,v andb 65535)) l |
| 224 | fun implode l = M (map (fn (v,p) => (v<<16)+p) l) |
| 225 | |
| 226 | val deg = let fun d([],n) = n | d(u::ul,n) = d(ul,(u andb 65535) + n) |
| 227 | in fn (M l) => d(l,0) end |
| 228 | |
| 229 | (* x^k > y^l if x>k or x=y and k>l *) |
| 230 | val compare = let |
| 231 | fun cmp ([],[]) = Util.Equal |
| 232 | | cmp (_::_,[]) = Util.Greater |
| 233 | | cmp ([],_::_) = Util.Less |
| 234 | | cmp ((u::us), (v::vs)) = if u=v then cmp (us,vs) |
| 235 | else if u<v then Util.Less |
| 236 | else (* u>v *) Util.Greater |
| 237 | in fn (M m,M m') => cmp(m,m') end |
| 238 | |
| 239 | fun display (M (l : int list)) : string = |
| 240 | let |
| 241 | fun dv v = if v<26 then chr (v+ord #"a") else chr (v-26+ord #"A") |
| 242 | fun d (vv,acc) = let val v = vv>>16 and p = vv andb 65535 |
| 243 | in if p=1 then dv v::acc |
| 244 | else |
| 245 | (dv v)::(String.explode (Int.toString p)) @ acc |
| 246 | end |
| 247 | in String.implode(fold d l []) end |
| 248 | |
| 249 | val multiply = let |
| 250 | fun mul ([],m) = m |
| 251 | | mul (m,[]) = m |
| 252 | | mul (u::us, v::vs) = let |
| 253 | val uu = u andb ~65536 |
| 254 | in if uu = (v andb ~65536) then let |
| 255 | val w = u + (v andb 65535) |
| 256 | in if uu = (w andb ~65536) then w::mul(us,vs) |
| 257 | else |
| 258 | (Util.illegal |
| 259 | (String.concat ["Mono.multiply overflow: ", |
| 260 | display (M(u::us)),", ", |
| 261 | display (M(v::vs))])) |
| 262 | end |
| 263 | else if u>v then u :: mul(us,v::vs) |
| 264 | else (* u<v *) v :: mul(u::us,vs) |
| 265 | end |
| 266 | in fn (M m,M m') => M (mul (m,m')) end |
| 267 | |
| 268 | val lcm = let |
| 269 | fun lcm ([],m) = m |
| 270 | | lcm (m,[]) = m |
| 271 | | lcm (u::us, v::vs) = |
| 272 | if u>=v then if (u andb ~65536)<v then u::lcm(us,vs) |
| 273 | else u::lcm(us,v::vs) |
| 274 | else if (v andb ~65536)<u then v::lcm(us,vs) |
| 275 | else v::lcm(u::us,vs) |
| 276 | in fn (M m,M m') => M (lcm (m,m')) end |
| 277 | val tryDivide = let |
| 278 | fun rev([],l) = l | rev(x::xs,l)=rev(xs,x::l) |
| 279 | fun d (m,[],q) = SOME(M(rev(q,m))) |
| 280 | | d ([],_::_,_) = NONE |
| 281 | | d (u::us,v::vs,q) = |
| 282 | if u<v then NONE |
| 283 | else if (u andb ~65536) = (v andb ~65536) then |
| 284 | if u=v then d(us,vs,q) else d(us,vs,u-(v andb 65535)::q) |
| 285 | else d(us,v::vs,u::q) |
| 286 | in fn (M m,M m') => d (m,m',[]) end |
| 287 | fun divide (m,m') = |
| 288 | case tryDivide(m,m') of SOME q => q | NONE => raise DoesntDivide |
| 289 | |
| 290 | end end (* local, structure M *) |
| 291 | |
| 292 | structure MI = struct (* MONO_IDEAL *) |
| 293 | |
| 294 | (* trie: |
| 295 | * index first by increasing order of vars |
| 296 | * children listed in increasing degree order |
| 297 | *) |
| 298 | datatype 'a mono_trie = MT of 'a option * (int * 'a mono_trie) list |
| 299 | (* tag, encoded (var,pwr) and children *) |
| 300 | datatype 'a mono_ideal = MI of (int * 'a mono_trie) ref |
| 301 | (* int maxDegree = least degree > all elements *) |
| 302 | |
| 303 | fun rev ([],l) = l | rev (x::xs,l) = rev(xs,x::l) |
| 304 | (* unused code |
| 305 | fun tl (_::l) = l | tl [] = raise (Util.Impossible "MONO_IDEAL.tl") |
| 306 | fun hd (x::_) = x | hd [] = raise (Util.Impossible "MONO_IDEAL.hd") |
| 307 | *) |
| 308 | val emptyTrie = MT(NONE,[]) |
| 309 | fun mkEmpty () = MI(ref (0,emptyTrie)) |
| 310 | |
| 311 | (* unused code unless searchDeg is used |
| 312 | fun maxDeg (MI(x)) = #1(!x) |
| 313 | *) |
| 314 | |
| 315 | val lshift = op << |
| 316 | (* unused code unless decode is used |
| 317 | val rshift = op >> |
| 318 | *) |
| 319 | val andb = op && |
| 320 | (* unused code |
| 321 | val orb = op || |
| 322 | *) |
| 323 | fun encode (var,pwr) = lshift(var,16)+pwr |
| 324 | (* unused code |
| 325 | fun decode vp = (rshift(vp,16),andb(vp,65535)) |
| 326 | *) |
| 327 | fun grabVar vp = andb(vp,~65536) |
| 328 | fun grabPwr vp = andb(vp,65535) |
| 329 | fun smallerVar (vp,vp') = vp < andb(vp',~65536) |
| 330 | |
| 331 | exception Found |
| 332 | fun search (MI(x),M.M m') = let |
| 333 | val (d,mt) = !x |
| 334 | val result = ref NONE |
| 335 | (* exception Found of M.mono * '_a *) |
| 336 | (* s works on remaining input mono, current output mono, tag, trie *) |
| 337 | fun s (_,m,MT(SOME a,_)) = |
| 338 | raise(result := SOME (M.M m,a); Found) |
| 339 | | s (m',m,MT(NONE,trie)) = s'(m',m,trie) |
| 340 | and s'([],_,_) = NONE |
| 341 | | s'(_,_,[]) = NONE |
| 342 | | s'(vp'::m',m,trie as (vp,child)::children) = |
| 343 | if smallerVar(vp',vp) then s'(m',m,trie) |
| 344 | else if grabPwr vp = 0 then (s(vp'::m',m,child); |
| 345 | s'(vp'::m',m,children)) |
| 346 | else if smallerVar(vp,vp') then NONE |
| 347 | else if vp<=vp' then (s(m',vp::m,child); |
| 348 | s'(vp'::m',m,children)) |
| 349 | else NONE |
| 350 | in s(rev(m',[]),[],mt) |
| 351 | handle Found (* (m,a) => SOME(m,a) *) => !result |
| 352 | end |
| 353 | |
| 354 | (* assume m is a new generator, i.e. not a multiple of an existing one *) |
| 355 | fun insert (MI (mi),m,a) = let |
| 356 | val (d,mt) = !mi |
| 357 | fun i ([],MT (SOME _,_)) = Util.illegal "MONO_IDEAL.insert duplicate" |
| 358 | | i ([],MT (NONE,children)) = MT(SOME a,children) |
| 359 | | i (vp::m,MT(a',[])) = MT(a',[(vp,i(m,emptyTrie))]) |
| 360 | | i (vp::m,mt as MT(a',trie as (vp',_)::_)) = let |
| 361 | fun j [] = [(vp,i(m,emptyTrie))] |
| 362 | | j ((vp',child)::children) = |
| 363 | if vp<vp' then (vp,i(m,emptyTrie))::(vp',child)::children |
| 364 | else if vp=vp' then (vp',i(m,child))::children |
| 365 | else (vp',child) :: j children |
| 366 | in |
| 367 | if smallerVar(vp,vp') then |
| 368 | MT(a',[(grabVar vp,MT(NONE,trie)),(vp,i(m,emptyTrie))]) |
| 369 | else if smallerVar(vp',vp) then i(grabVar vp'::vp::m,mt) |
| 370 | else MT(a',j trie) |
| 371 | end |
| 372 | in mi := (Int.max(d,M.deg m),i (rev(map encode(M.explode m),[]),mt)) end |
| 373 | |
| 374 | fun mkIdeal [] = mkEmpty() |
| 375 | | mkIdeal (orig_ms : (M.mono * '_a) list)= let |
| 376 | fun ins ((m,a),arr) = Util.insert((m,a),M.deg m,arr) |
| 377 | val msa = arrayoflist orig_ms |
| 378 | val ms : (M.mono * '_a) list = |
| 379 | Util.stripSort (fn ((m,_),(m',_)) => M.compare (m,m')) msa |
| 380 | val buckets = revfold ins ms (array1(0,[])) |
| 381 | val n = length1 buckets |
| 382 | val mi = mkEmpty() |
| 383 | fun sort i = if i>=n then mi else let |
| 384 | fun redundant (m,_) = case search(mi,m) of NONE => false |
| 385 | | SOME _ => true |
| 386 | fun filter ([],l) = app (fn (m,a) => insert(mi,m,a)) l |
| 387 | | filter (x::xx,l) = if redundant x then filter(xx,l) |
| 388 | else filter(xx,x::l) |
| 389 | in filter(sub1(buckets,i),[]); |
| 390 | update1(buckets,i,[]); |
| 391 | sort(i+1) |
| 392 | end |
| 393 | in sort 0 end |
| 394 | |
| 395 | fun fold g (MI(x)) init = let |
| 396 | val (_,mt) = !x |
| 397 | fun f(acc,m,MT(NONE,children)) = f'(acc,m,children) |
| 398 | | f(acc,m,MT(SOME a,children)) = |
| 399 | f'(g((M.M m,a),acc),m,children) |
| 400 | and f'(acc,m,[]) = acc |
| 401 | | f'(acc,m,(vp,child)::children) = |
| 402 | if grabPwr vp=0 then f'(f(acc,m,child),m,children) |
| 403 | else f'(f(acc,vp::m,child),m,children) |
| 404 | in f(init,[],mt) end |
| 405 | |
| 406 | (* unused code |
| 407 | fun searchDeg (mi,d) = |
| 408 | if d>maxDeg mi then [] |
| 409 | else fold (fn ((m,a),l) => if M.deg m=d then (m,a)::l else l) mi [] |
| 410 | *) |
| 411 | |
| 412 | end (* structure MI *) |
| 413 | |
| 414 | |
| 415 | val log = let fun log(n,l) = if n<=1 then l else log((n >> 1),1+l) |
| 416 | in fn n => log(n,0) end |
| 417 | val maxLeft = ref 0 |
| 418 | val maxRight = ref 0 |
| 419 | val counts = tabulate(20,fn _ => array1(20,0)) |
| 420 | val indices = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] |
| 421 | |
| 422 | (* unused code |
| 423 | fun resetCounts() = app(fn i => app (fn j => update1(sub1(counts,i),j,0)) indices) indices |
| 424 | *) |
| 425 | |
| 426 | fun pair(l,r) = let |
| 427 | val l = log l and r = log r |
| 428 | val _ = maxLeft := Int.max(!maxLeft,l) and _ = maxRight := Int.max(!maxRight,r) |
| 429 | val a = sub1(counts,l) |
| 430 | in update1(a,r,sub1(a,r)+1) end |
| 431 | (* unused code unless printCounts is used |
| 432 | fun getCounts () = |
| 433 | map (fn i => map (fn j => sub1(sub1(counts,i),j)) indices) indices |
| 434 | *) |
| 435 | |
| 436 | structure P = struct |
| 437 | |
| 438 | datatype poly = P of (F.field*M.mono) list (* descending mono order *) |
| 439 | (* |
| 440 | fun show (P x) = (print "[ "; |
| 441 | app (fn (f, m) => |
| 442 | (print "("; F.show f; print ","; M.show m; print ") ")) x; |
| 443 | print " ]") |
| 444 | *) |
| 445 | val zero = P [] |
| 446 | (* unused code unless power is used |
| 447 | val one = P [(F.one,M.one)] |
| 448 | *) |
| 449 | (* unused code |
| 450 | fun coerceInt n = P [(F.coerceInt n,M.one)] |
| 451 | fun coerceField a = P [(a,M.one)] |
| 452 | fun coerceMono m = P [(F.one,m)] |
| 453 | *) |
| 454 | fun coerce (a,m) = P [(a,m)] |
| 455 | fun implode p = P p |
| 456 | fun cons (am,P p) = P (am::p) |
| 457 | |
| 458 | local |
| 459 | fun neg p = (map (fn (a,m) => (F.negate a,m)) p) |
| 460 | fun plus ([],p2) = p2 |
| 461 | | plus (p1,[]) = p1 |
| 462 | | plus ((a,m)::ms,(b,n)::ns) = case M.compare(m,n) of |
| 463 | Util.Less => (b,n) :: plus ((a,m)::ms,ns) |
| 464 | | Util.Greater => (a,m) :: plus (ms,(b,n)::ns) |
| 465 | | Util.Equal => let val c = F.add(a,b) |
| 466 | in if F.isZero c then plus(ms,ns) |
| 467 | else (c,m)::plus(ms,ns) |
| 468 | end |
| 469 | fun minus ([],p2) = neg p2 |
| 470 | | minus (p1,[]) = p1 |
| 471 | | minus ((a,m)::ms,(b,n)::ns) = case M.compare(m,n) of |
| 472 | Util.Less => (F.negate b,n) :: minus ((a,m)::ms,ns) |
| 473 | | Util.Greater => (a,m) :: minus (ms,(b,n)::ns) |
| 474 | | Util.Equal => let val c = F.subtract(a,b) |
| 475 | in if F.isZero c then minus(ms,ns) |
| 476 | else (c,m)::minus(ms,ns) |
| 477 | end |
| 478 | fun termMult (a,m,p) = |
| 479 | (map (fn (a',m') => (F.multiply(a,a'),M.multiply(m,m'))) p) |
| 480 | in |
| 481 | (* unused code |
| 482 | fun negate (P p) = P (neg p) |
| 483 | *) |
| 484 | fun add (P p1,P p2) = (pair(length p1,length p2); P (plus(p1,p2))) |
| 485 | fun subtract (P p1,P p2) = (pair(length p1,length p2); P (minus(p1,p2))) |
| 486 | |
| 487 | (* unused code unless power is used |
| 488 | val multiply = let |
| 489 | fun times (p1,p2) = |
| 490 | revfold (fn ((a,m),tot) => plus (termMult(a,m,p2),tot)) p1 [] |
| 491 | in fn (P p1,P p2) => if length p1 > length p2 then P(times (p2,p1)) |
| 492 | else P(times (p1,p2)) |
| 493 | end |
| 494 | *) |
| 495 | |
| 496 | (* unused code |
| 497 | fun singleReduce (P y,a,m,P x) = |
| 498 | (pair(length y,length x); P(minus(y,termMult(a,m,x)))) |
| 499 | *) |
| 500 | |
| 501 | fun spair (a,m,P f,b,n,P g) = |
| 502 | (pair(length f,length g); P(minus(termMult(a,m,f),termMult(b,n,g)))) |
| 503 | val termMult = fn (a,m,P f) => P(termMult(a,m,f)) |
| 504 | end |
| 505 | |
| 506 | (* unused code unless power is used |
| 507 | val rshift = op >> |
| 508 | val lshift = op << |
| 509 | *) |
| 510 | |
| 511 | (* unused code |
| 512 | val andb = op && |
| 513 | val orb = op || |
| 514 | *) |
| 515 | |
| 516 | fun scalarMult (a,P p) = P (map (fn (b,m) => (F.multiply(a,b),m)) p) |
| 517 | |
| 518 | |
| 519 | (* unused code |
| 520 | fun power(p,k) = |
| 521 | if k<=3 then case k of |
| 522 | 0 => one |
| 523 | | 1 => p |
| 524 | | 2 => multiply(p,p) |
| 525 | | 3 => multiply(p,multiply(p,p)) |
| 526 | | _ => Util.illegal "POLY.power with k<0" |
| 527 | else if andb(k,1)=0 then power(multiply(p,p),rshift(k,1)) |
| 528 | else multiply(p,power(multiply(p,p),rshift(k,1))) |
| 529 | *) |
| 530 | |
| 531 | fun isZero (P []) = true | isZero (P (_::_)) = false |
| 532 | |
| 533 | (* unused code |
| 534 | val equal = let |
| 535 | fun eq ([],[]) = true |
| 536 | | eq (_::_,[]) = false |
| 537 | | eq ([],_::_) = false |
| 538 | | eq ((a,m)::p,(b,n)::q) = |
| 539 | F.equal(a,b) andalso M.compare(m,n)=Util.Equal |
| 540 | andalso eq (p,q) |
| 541 | in fn (P p,P q) => eq (p,q) end |
| 542 | *) |
| 543 | |
| 544 | (* these should only be called if there is a leading term, i.e. poly<>0 *) |
| 545 | (* unused code |
| 546 | fun leadTerm (P(am::_)) = am |
| 547 | | leadTerm (P []) = Util.illegal "POLY.leadTerm" |
| 548 | *) |
| 549 | fun leadMono (P((_,m)::_)) = m |
| 550 | | leadMono (P []) = Util.illegal "POLY.leadMono" |
| 551 | fun leadCoeff (P((a,_)::_)) = a |
| 552 | | leadCoeff (P []) = Util.illegal "POLY.leadCoeff" |
| 553 | fun rest (P (_::p)) = P p |
| 554 | | rest (P []) = Util.illegal "POLY.rest" |
| 555 | fun leadAndRest (P (lead::rest)) = (lead,P rest) |
| 556 | | leadAndRest (P []) = Util.illegal "POLY.leadAndRest" |
| 557 | |
| 558 | fun deg (P []) = Util.illegal "POLY.deg on zero poly" |
| 559 | | deg (P ((_,m)::_)) = M.deg m (* homogeneous poly *) |
| 560 | fun numTerms (P p) = length p |
| 561 | |
| 562 | (* only used if r is used |
| 563 | fun display (P []) = F.display F.zero |
| 564 | | display (P p) = let |
| 565 | fun dsp (a,m) = let |
| 566 | val s = |
| 567 | if M.deg m = 0 then F.display a |
| 568 | else if F.equal(F.one,F.negate a) then "-" ^ M.display m |
| 569 | else if F.equal(F.one,a) then M.display m |
| 570 | else F.display a ^ M.display m |
| 571 | in if substring(s,0,1)="-" then s else "+" ^ s end |
| 572 | in String.concat(map dsp p) end |
| 573 | *) |
| 574 | end |
| 575 | |
| 576 | structure HP = struct |
| 577 | datatype hpoly = HP of P.poly array1 |
| 578 | val log = let |
| 579 | fun log(n,l) = if n<8 then l else log((n >> 2),1+l) |
| 580 | in fn n => log(n,0) end |
| 581 | fun mkHPoly p = let |
| 582 | val l = log(P.numTerms p) |
| 583 | in HP(tabulate(l+1,fn i => if i=l then p else P.zero)) end |
| 584 | fun add(p,HP ps) = let |
| 585 | val l = log(P.numTerms p) |
| 586 | in if l>=length1 ps then let |
| 587 | val n = length1 ps |
| 588 | in HP(tabulate(n+n, |
| 589 | fn i => if i<n then sub1(ps,i) |
| 590 | else if i=l then p else P.zero)) |
| 591 | end |
| 592 | else let |
| 593 | val p = P.add(p,sub1(ps,l)) |
| 594 | in if l=log(P.numTerms p) then (update1(ps,l,p); HP ps) |
| 595 | else (update1(ps,l,P.zero); add (p,HP ps)) |
| 596 | end |
| 597 | end |
| 598 | fun leadAndRest (HP ps) = let |
| 599 | val n = length1 ps |
| 600 | fun lar (m,indices,i) = if i>=n then lar'(m,indices) else let |
| 601 | val p = sub1(ps,i) |
| 602 | in if P.isZero p then lar(m,indices,i+1) |
| 603 | else if null indices then lar(P.leadMono p,[i],i+1) |
| 604 | else case M.compare(m,P.leadMono p) of |
| 605 | Util.Less => lar(P.leadMono p,[i],i+1) |
| 606 | | Util.Equal => lar(m,i::indices,i+1) |
| 607 | | Util.Greater => lar(m,indices,i+1) |
| 608 | end |
| 609 | and lar' (_,[]) = NONE |
| 610 | | lar' (m,i::is) = let |
| 611 | fun extract i = case P.leadAndRest(sub1(ps,i)) of |
| 612 | ((a,_),rest) => (update1(ps,i,rest); a) |
| 613 | val a = revfold (fn (j,b) => F.add(extract j,b)) |
| 614 | is (extract i) |
| 615 | in if F.isZero a then lar(M.one,[],0) else SOME(a,m,HP ps) |
| 616 | end |
| 617 | in lar(M.one,[],0) end |
| 618 | end |
| 619 | |
| 620 | structure G = struct |
| 621 | val autoReduce = ref true |
| 622 | val maxDeg = ref 10000 |
| 623 | val maybePairs = ref 0 |
| 624 | val primePairs = ref 0 |
| 625 | val usedPairs = ref 0 |
| 626 | val newGens = ref 0 |
| 627 | |
| 628 | fun reset () = (maybePairs:=0; primePairs:=0; usedPairs:=0; newGens:=0) |
| 629 | |
| 630 | fun inc r = r := !r + 1 |
| 631 | |
| 632 | fun reduce (f,mi) = if P.isZero f then f else let |
| 633 | (* use accumulator and reverse at end? *) |
| 634 | fun r hp = case HP.leadAndRest hp of |
| 635 | NONE => [] |
| 636 | | (SOME(a,m,hp)) => case MI.search(mi,m) of |
| 637 | NONE => (a,m)::(r hp) |
| 638 | | SOME (m',p) => r (HP.add(P.termMult(F.negate a,M.divide(m,m'),!p),hp)) |
| 639 | in P.implode(r (HP.mkHPoly f)) end |
| 640 | |
| 641 | (* assume f<>0 *) |
| 642 | fun mkMonic f = P.scalarMult(F.reciprocal(P.leadCoeff f),f) |
| 643 | |
| 644 | (* given monic h, a monomial ideal mi of m's tagged with g's representing |
| 645 | * an ideal (g1,...,gn): a poly g is represented as (lead mono m,rest of g). |
| 646 | * update pairs to include new s-pairs induced by h on g's: |
| 647 | * 1) compute minimal gi1...gik so that <gij:h's> generate <gi:h's>, i.e. |
| 648 | * compute monomial ideal for gi:h's tagged with gi |
| 649 | * 2) toss out gij's whose lead mono is rel. prime to h's lead mono (why?) |
| 650 | * 3) put (h,gij) pairs into degree buckets: for h,gij with lead mono's m,m' |
| 651 | * deg(h,gij) = deg lcm(m,m') = deg (lcm/m) + deg m = deg (m':m) + deg m |
| 652 | * 4) store list of pairs (h,g1),...,(h,gn) as vector (h,g1,...,gn) |
| 653 | *) |
| 654 | fun addPairs (h,mi,pairs) = let |
| 655 | val m = P.leadMono h |
| 656 | val d = M.deg m |
| 657 | fun tag ((m' : M.mono,g' : P.poly ref),quots) = (inc maybePairs; |
| 658 | (M.divide(M.lcm(m,m'),m),(m',!g'))::quots) |
| 659 | fun insert ((mm,(m',g')),arr) = (* recall mm = m':m *) |
| 660 | if M.compare(m',mm)=Util.Equal then (* rel. prime *) |
| 661 | (inc primePairs; arr) |
| 662 | else (inc usedPairs; |
| 663 | Util.insert(P.cons((F.one,m'),g'),M.deg mm+d,arr)) |
| 664 | val buckets = MI.fold insert (MI.mkIdeal (MI.fold tag mi [])) |
| 665 | (array1(0,[])) |
| 666 | fun ins (~1,pairs) = pairs |
| 667 | | ins (i,pairs) = case sub1(buckets,i) of |
| 668 | [] => ins(i-1,pairs) |
| 669 | | gs => ins(i-1,Util.insert(arrayoflist(h::gs),i,pairs)) |
| 670 | in ins(length1 buckets - 1,pairs) end |
| 671 | |
| 672 | fun grobner fs = let |
| 673 | fun pr l = print (String.concat (l@["\n"])) |
| 674 | val fs = revfold (fn (f,fs) => Util.insert(f,P.deg f,fs)) |
| 675 | fs (array1(0,[])) |
| 676 | (* pairs at least as long as fs, so done when done w/ all pairs *) |
| 677 | val pairs = ref(array1(length1 fs,[])) |
| 678 | val mi = MI.mkEmpty() |
| 679 | val newDegGens = ref [] |
| 680 | val addGen = (* add and maybe auto-reduce new monic generator h *) |
| 681 | if not(!autoReduce) then |
| 682 | fn h => MI.insert (mi,P.leadMono h,ref (P.rest h)) |
| 683 | else fn h => let |
| 684 | val ((_,m),rh) = P.leadAndRest h |
| 685 | fun autoReduce f = |
| 686 | if P.isZero f then f |
| 687 | else let val ((a,m'),rf) = P.leadAndRest f |
| 688 | in case M.compare(m,m') of |
| 689 | Util.Less => P.cons((a,m'),autoReduce rf) |
| 690 | | Util.Equal => P.subtract(rf,P.scalarMult(a,rh)) |
| 691 | | Util.Greater => f |
| 692 | end |
| 693 | val rrh = ref rh |
| 694 | in |
| 695 | MI.insert (mi,P.leadMono h,rrh); |
| 696 | app (fn f => f:=autoReduce(!f)) (!newDegGens); |
| 697 | newDegGens := rrh :: !newDegGens |
| 698 | end |
| 699 | val tasksleft = ref 0 |
| 700 | fun feedback () = let |
| 701 | val n = !tasksleft |
| 702 | in |
| 703 | if (n && 15)=0 then print (Int.toString n) else (); |
| 704 | print "."; |
| 705 | TextIO.flushOut TextIO.stdOut; |
| 706 | tasksleft := n-1 |
| 707 | end |
| 708 | |
| 709 | fun try h = |
| 710 | let |
| 711 | val _ = feedback () |
| 712 | val h = reduce(h,mi) |
| 713 | in if P.isZero h |
| 714 | then () |
| 715 | else let val h = mkMonic h |
| 716 | val _ = (print "#"; TextIO.flushOut TextIO.stdOut) |
| 717 | in pairs := addPairs(h,mi,!pairs); |
| 718 | addGen h; |
| 719 | inc newGens |
| 720 | end |
| 721 | end |
| 722 | |
| 723 | fun tryPairs fgs = let |
| 724 | val ((a,m),f) = P.leadAndRest (sub1(fgs,0)) |
| 725 | fun tryPair i = if i=0 then () else let |
| 726 | val ((b,n),g) = P.leadAndRest (sub1(fgs,i)) |
| 727 | val k = M.lcm(m,n) |
| 728 | in |
| 729 | try (P.spair(b,M.divide(k,m),f,a,M.divide(k,n),g)); |
| 730 | tryPair (i-1) |
| 731 | end |
| 732 | in tryPair (length1 fgs -1) end |
| 733 | |
| 734 | fun numPairs ([],n) = n |
| 735 | | numPairs (p::ps,n) = numPairs(ps,n-1+length1 p) |
| 736 | |
| 737 | fun gb d = if d>=length1(!pairs) then mi else |
| 738 | (* note: i nullify entries to reclaim space *) |
| 739 | ( |
| 740 | pr ["DEGREE ",Int.toString d," with ", |
| 741 | Int.toString(numPairs(sub1(!pairs,d),0))," pairs ", |
| 742 | if d>=length1 fs then "0" else Int.toString(length(sub1(fs,d))), |
| 743 | " generators to do"]; |
| 744 | tasksleft := numPairs(sub1(!pairs,d),0); |
| 745 | if d>=length1 fs then () |
| 746 | else tasksleft := !tasksleft + length (sub1(fs,d)); |
| 747 | if d>(!maxDeg) then () |
| 748 | else ( |
| 749 | reset(); |
| 750 | newDegGens := []; |
| 751 | app tryPairs (sub1(!pairs,d)); |
| 752 | update1(!pairs,d,[]); |
| 753 | if d>=length1 fs then () |
| 754 | else (app try (sub1(fs,d)); update1(fs,d,[])); |
| 755 | pr ["maybe ",Int.toString(!maybePairs)," prime ", |
| 756 | Int.toString (!primePairs), |
| 757 | " using ",Int.toString (!usedPairs), |
| 758 | "; found ",Int.toString (!newGens)] |
| 759 | ); |
| 760 | gb(d+1) |
| 761 | ) |
| 762 | in gb 0 end |
| 763 | |
| 764 | local |
| 765 | (* grammar: |
| 766 | dig ::= 0 | ... | 9 |
| 767 | var ::= a | ... | z | A | ... | Z |
| 768 | sign ::= + | - |
| 769 | nat ::= dig | nat dig |
| 770 | mono ::= | var mono | var num mono |
| 771 | term ::= nat mono | mono |
| 772 | poly ::= term | sign term | poly sign term |
| 773 | *) |
| 774 | datatype char = Dig of int | Var of int | Sign of int |
| 775 | fun char ch = |
| 776 | let val och = ord ch in |
| 777 | if ord #"0"<=och andalso och<=ord #"9" then Dig (och - ord #"0") |
| 778 | else if ord #"a"<=och andalso och<=ord #"z" then Var (och - ord #"a") |
| 779 | else if ord #"A"<=och andalso och<=ord #"Z" then Var (och - ord #"A" + 26) |
| 780 | else if och = ord #"+" then Sign 1 |
| 781 | else if och = ord #"-" then Sign ~1 |
| 782 | else Util.illegal ("bad ch in poly: " ^ (Char.toString(ch))) |
| 783 | end |
| 784 | |
| 785 | fun nat (n,Dig d::l) = nat(n*10+d,l) | nat (n,l) = (n,l) |
| 786 | fun mono (m,Var v::Dig d::l) = |
| 787 | let val (n,l) = nat(d,l) |
| 788 | in mono(M.multiply(M.implode[(v,n)],m),l) end |
| 789 | | mono (m,Var v::l) = mono(M.multiply(M.x_i v,m),l) |
| 790 | | mono (m,l) = (m,l) |
| 791 | |
| 792 | fun term l = let |
| 793 | val (n,l) = case l of (Dig d::l) => nat(d,l) | _ => (1,l) |
| 794 | val (m,l) = mono(M.one,l) |
| 795 | in ((F.coerceInt n,m),l) end |
| 796 | fun poly (p,[]) = p |
| 797 | | poly (p,l) = let |
| 798 | val (s,l) = case l of Sign s::l => (F.coerceInt s,l) | _ => (F.one,l) |
| 799 | val ((a,m),l) = term l |
| 800 | in poly(P.add(P.coerce(F.multiply(s,a),m),p),l) end |
| 801 | |
| 802 | in |
| 803 | fun parsePoly s = poly (P.zero,map char(String.explode s)) |
| 804 | |
| 805 | (* unused code |
| 806 | fun readIdeal stream = let |
| 807 | fun readLine () = let |
| 808 | val s = input_line stream |
| 809 | val n = size s |
| 810 | val n = if n>0 andalso substring(s,n-1,1)="\n" then n-1 else n |
| 811 | fun r i = if i>=n then [] |
| 812 | else case substring(s,i,1) of |
| 813 | ";" => r(i+1) |
| 814 | | " " => r(i+1) |
| 815 | | _ => map char (String.explode(substring(s,i,n-i))) |
| 816 | in r 0 end |
| 817 | fun r () = if end_of_stream stream then [] |
| 818 | else poly(P.zero,readLine()) :: r() |
| 819 | fun num() = if end_of_stream stream then Util.illegal "missing #" |
| 820 | else case nat(0,readLine()) of |
| 821 | (_,_::_) => Util.illegal "junk after #" |
| 822 | | (n,_) => n |
| 823 | val _ = 1=num() orelse Util.illegal "stream doesn't start w/ `1'" |
| 824 | val n = num() |
| 825 | val i = r() |
| 826 | val _ = length i = n orelse Util.illegal "wrong # poly's" |
| 827 | in i end |
| 828 | *) |
| 829 | |
| 830 | (* unused code |
| 831 | fun read filename = let |
| 832 | val stream = open_in filename |
| 833 | val i = readIdeal stream |
| 834 | val _ = close_in stream |
| 835 | in i end |
| 836 | *) |
| 837 | end (* local *) |
| 838 | |
| 839 | end (* structure G *) |
| 840 | |
| 841 | |
| 842 | |
| 843 | val _ = G.maxDeg:=1000000 |
| 844 | |
| 845 | fun grab mi = MI.fold (fn ((m,g),l) => P.cons((F.one,m),!g)::l) mi [] |
| 846 | |
| 847 | (* unused code |
| 848 | fun r mi s = let |
| 849 | val p = G.parsePoly s |
| 850 | in print (P.display p); print "\n"; |
| 851 | print (P.display(G.reduce(p,mi))); print "\n" |
| 852 | end |
| 853 | *) |
| 854 | |
| 855 | (* unused code unless printCounts is used |
| 856 | fun p6 i= let val s= Int.toString (i:int) |
| 857 | val n= size s |
| 858 | in print(substring(" ",0,6-n)); print s end |
| 859 | *) |
| 860 | |
| 861 | (* unused code |
| 862 | fun hex n = let |
| 863 | fun h n = if n=0 then "" |
| 864 | else h(n smlnj_div 16) ^ substring("0123456789ABCDEF",n smlnj_mod 16,1) |
| 865 | in if n=0 then "0" else h n end |
| 866 | fun printCounts () = map (fn l => (map p6 l; print "\n")) (getCounts()) |
| 867 | fun totalCount () = revfold (fn (l,c) => revfold op + l c) (getCounts()) 0 |
| 868 | *) |
| 869 | |
| 870 | (* unused code |
| 871 | fun maxCount () = revfold (fn (l,m) => revfold Int.max l m) (getCounts()) 0 |
| 872 | *) |
| 873 | |
| 874 | (* unused code unless analyze is used |
| 875 | fun terms (p,tt) = if P.isZero p then tt else terms(P.rest p,P.leadMono p::tt) |
| 876 | fun tails ([],tt) = tt |
| 877 | | tails (t as _::t',tt) = tails (t',t::tt) |
| 878 | *) |
| 879 | |
| 880 | (* Unused code unless sort (analyze) is used |
| 881 | local |
| 882 | val a = 16807.0 and m = 2147483647.0 |
| 883 | in |
| 884 | val seed = ref 1.0 |
| 885 | fun random n = let val t = a*(!seed) |
| 886 | in seed := t - m * real(floor(t/m)); |
| 887 | floor(real n * !seed/m) |
| 888 | end |
| 889 | end |
| 890 | *) |
| 891 | |
| 892 | (* Unused code unless analyze is used |
| 893 | fun sort [] = [] |
| 894 | | sort a = |
| 895 | let |
| 896 | val a = arrayoflist a |
| 897 | val b = tabulate(length1 a,fn i => i) |
| 898 | val sub = sub1 and update = update1 |
| 899 | infix sub |
| 900 | fun swap (i,j) = let val ai = a sub i in update(a,i,a sub j); update(a,j,ai) end |
| 901 | (* sort all k, 0<=i<=k<j<=length a *) |
| 902 | fun s (i,j,acc) = if i=j then acc else let |
| 903 | val pivot = a sub (b sub (i+random(j-i))) |
| 904 | fun partition (dup,lo,k,hi) = if k=hi then (dup,lo,hi) else |
| 905 | (case M.compare (pivot, a sub (b sub k)) of |
| 906 | Util.Less => (swap (lo,k); partition (dup,lo+1,k+1,hi)) |
| 907 | | Util.Equal => partition (dup+1,lo,k+1,hi) |
| 908 | | Util.Greater => (swap (k,hi-1); partition (dup,lo,k,hi-1))) |
| 909 | val (dup,lo,hi) = partition (0,i,i,j) |
| 910 | in s(i,lo,(dup,pivot)::s(hi,j,acc)) end |
| 911 | in s(0,length1 a,[]) end |
| 912 | *) |
| 913 | |
| 914 | (* Unused code unless analyze is used |
| 915 | fun sum f l = revfold op + (map f l) 0 |
| 916 | *) |
| 917 | |
| 918 | (* Unused code included in the benchmark |
| 919 | fun analyze gb = let |
| 920 | val aa = revfold terms gb [] |
| 921 | val bb = map M.explode aa |
| 922 | val aa = sort aa |
| 923 | fun len m = length (M.explode m) |
| 924 | fun prt (s:string) (i:int) = (print s; print(Int.toString i); print "\n"; i) |
| 925 | val m= sum #1 aa |
| 926 | val u= length aa |
| 927 | val cm =sum (fn (d,l) => d*len l) aa |
| 928 | val cu =sum (len o #2) aa |
| 929 | val for=length(sort(map M.implode (revfold tails bb []))) |
| 930 | val bak=length(sort(map (M.implode o rev) (revfold tails (map rev bb) []))) |
| 931 | in |
| 932 | {m=prt "m = " m, u=prt "u = " u, cm =prt "cm = " cm, cu =prt "cu = " cu, for=prt "for= " for, bak=prt "bak= " bak} |
| 933 | end |
| 934 | *) |
| 935 | |
| 936 | fun gb fs = |
| 937 | let |
| 938 | val g = G.grobner fs handle (Util.Illegal s) => (print s; raise Div) |
| 939 | val fs = grab g |
| 940 | fun info f = app print |
| 941 | [M.display(P.leadMono f), |
| 942 | " + ", Int.toString(P.numTerms f - 1), " terms\n"] |
| 943 | in app info fs end |
| 944 | |
| 945 | |
| 946 | fun report (e as Tabulate) = (print "exn: Tabulate\n"; raise e) |
| 947 | | report (e as ArrayofList) = (print "exn: ArrayofList\n"; raise e) |
| 948 | | report (e as (Util.NotImplemented s)) = |
| 949 | (print ("exn: NotImplemented " ^ s ^ "\n"); raise e) |
| 950 | | report (e as (Util.Impossible s)) = |
| 951 | (print ("exn: Impossible " ^ s ^ "\n"); raise e) |
| 952 | | report (e as (Util.Illegal s)) = |
| 953 | (print ("exn: Illegal " ^ s ^ "\n"); raise e) |
| 954 | | report (e as (M.DoesntDivide)) = (print ("exn: DoesntDivide\n"); raise e) |
| 955 | | report (e as (MI.Found)) = (print ("exn: Found\n"); raise e) |
| 956 | |
| 957 | |
| 958 | (* rather long running test case *) |
| 959 | (* val fs = map G.parsePoly |
| 960 | * ["-El-Dh-Cd+Bo+xn+tm","-Fo+Ep-Ek-Dg-Cc+Ao+wn+sm","-Fn-Ej+Dp-Df-Cb+zo+vn+rm", |
| 961 | * "-Fm-Ei-De+Cp-Ca+yo+un+qm","Fl-Bp+Bk-Al-zh-yd+xj+ti","El-Bo-zg-yc+wj+si", |
| 962 | * "Dl-Bn-Aj+zk-zf-yb+vj+ri","Cl-Bm-Ai-ze+yk-ya+uj+qi", |
| 963 | * "Fh+Bg-xp+xf-wl-vh-ud+te","Eh+Ag-xo-wk+wf-vg-uc+se","Dh+zg-xn-wj-ub+re", |
| 964 | * "Ch+yg-xm-wi-ve+uf-ua+qe","Fd+Bc+xb-tp+ta-sl-rh-qd", |
| 965 | * "Ed+Ac+wb-to-sk+sa-rg-qc","Dd+zc+vb-tn-sj-rf+ra-qb","Cd+yc+ub-tm-si-re"] |
| 966 | *) |
| 967 | |
| 968 | (* rather long running test case *) |
| 969 | (* val u7 = map G.parsePoly |
| 970 | * ["abcdefg-h7","a+b+c+d+e+f+g","ab+bc+cd+de+ef+fg+ga", |
| 971 | * "abc+bcd+cde+def+efg+fga+gab","abcd+bcde+cdef+defg+efga+fgab+gabc", |
| 972 | * "abcde+bcdef+cdefg+defga+efgab+fgabc+gabcd", |
| 973 | * "abcdef+bcdefg+cdefga+defgab+efgabc+fgabcd+gabcde"] |
| 974 | *) |
| 975 | |
| 976 | |
| 977 | (* val u5 = map G.parsePoly ["abcde-f5","a+b+c+d+e","ab+bc+cd+de+ea", |
| 978 | * "abc+bcd+cde+dea+eab","abcd+bcde+cdea+deab+eabc"] |
| 979 | * |
| 980 | * val u4 = map G.parsePoly ["abcd-e4","a+b+c+d","ab+bc+cd+da","abc+bcd+cda+dab"] |
| 981 | * |
| 982 | *) |
| 983 | |
| 984 | (* fun runit () = |
| 985 | * let |
| 986 | * val _ = (print "Enter fs, u7, u6, u5, or u4: "; |
| 987 | * TextIO.flushOut TextIO.stdOut) |
| 988 | * val s = TextIO.inputN(TextIO.stdIn,2) |
| 989 | * val data = |
| 990 | * if (s = "fs") then fs else if (s = "u7") then u7 |
| 991 | * else if (s = "u6") then u6 else if (s = "u5") then u5 |
| 992 | * else if (s = "u4") then u4 else |
| 993 | * (print "no such data\n"; raise (Util.Impossible "no such data")) |
| 994 | * in |
| 995 | * gb data handle e => report e |
| 996 | * end |
| 997 | *) |
| 998 | |
| 999 | structure Main = |
| 1000 | struct |
| 1001 | fun doit n = |
| 1002 | let |
| 1003 | val u6 = |
| 1004 | map G.parsePoly |
| 1005 | ["abcdef-g6","a+b+c+d+e+f","ab+bc+cd+de+ef+fa", |
| 1006 | "abc+bcd+cde+def+efa+fab", |
| 1007 | "abcd+bcde+cdef+defa+efab+fabc", |
| 1008 | "abcde+bcdef+cdefa+defab+efabc+fabcd"] |
| 1009 | fun loop n = |
| 1010 | if n = 0 |
| 1011 | then () |
| 1012 | else (gb u6; loop (n - 1)) |
| 1013 | in |
| 1014 | loop n |
| 1015 | end |
| 1016 | end |