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
.
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
))
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
},
19 * pages
= {285 -- 293},
24 val print
: string -> unit
= print
25 type 'a array1
= 'a Array
.array
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
))))
35 fun fold f l b
= List.foldl f b l
36 fun revfold f l b
= List.foldr f b l
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
45 val smlnj_mod
= op mod
46 val smlnj_div
= op div
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
60 fun arrayoflist (hd
::tl
) =
61 let val a
= array1((length tl
) + 1,hd
)
63 |
al(hd
::tl
,i
) = (update1(a
,i
,hd
); al(tl
,i
+1))
67 |
arrayoflist ([]) = raise ArrayofList
70 structure Util
= struct
71 datatype relation
= Less | Equal | Greater
73 exception NotImplemented
of string
74 exception Impossible
of string (* flag
"impossible" condition
*)
75 exception Illegal
of string (* flag function use violating precondition
*)
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
82 (* arr
[i
] := obj
:: arr
[i
]; extend non
-empty arr
if necessary
*)
83 fun insert (obj
,i
,arr
) = let
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
));
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
)
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
.
109 fun stripSort compare
= fn a
=> let
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
,[])
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
)
145 (* unused code unless P
.display is used
149 fun coerceInt n
= F (n smlnj_mod p
)
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
.
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
169 fun divide (n
,m
) = multiply (n
, reciprocal m
)
172 (* unused code unless power is used
179 if k
<=3 then case k
of
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)))
189 fun isZero (F n
) = n
=0
190 (* unused codeunless P
.display is used
191 fun equal (F n
,F m
) = n
=m
193 fun display (F n
) = if n
<=p smlnj_div
2 then Int.toString n
194 else "-" ^
Int.toString (p
-n
)
198 structure M
= struct (* MONO
*)
202 (* val op << = Bits
.lshift
and op >> = Bits
.rshift
and op andb
= Bits
.andb
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
211 datatype mono
= M
of int list
213 fun show (M x
) = (print
"<"; app (fn i
=> (print (Int.toString i
); print
",")) x
; print
">")
215 exception DoesntDivide
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
)
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
229 (* x^k
> y^l
if x
>k or x
=y
and k
>l
*)
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
239 fun display (M (l
: int list
)) : string =
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
245 (dv v
)::(String.explode (Int.toString p
)) @ acc
247 in String.implode(fold d l
[]) end
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
)
259 (String.concat
["Mono.multiply overflow: ",
260 display (M(u
::us
)),", ",
261 display (M(v
::vs
))]))
263 else if u
>v
then u
:: mul(us
,v
::vs
)
264 else (* u
<v
*) v
:: mul(u
::us
,vs
)
266 in fn (M m
,M m
') => M (mul (m
,m
')) end
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
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
) =
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
288 case tryDivide(m
,m
') of SOME q
=> q | NONE
=> raise DoesntDivide
290 end end (* local, structure M
*)
292 structure MI
= struct (* MONO_IDEAL
*)
295 * index first by increasing order
of vars
296 * children listed
in increasing degree order
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
*)
303 fun rev ([],l
) = l |
rev (x
::xs
,l
) = rev(xs
,x
::l
)
305 fun tl (_
::l
) = l | tl
[] = raise (Util
.Impossible
"MONO_IDEAL.tl")
306 fun hd (x
::_
) = x | hd
[] = raise (Util
.Impossible
"MONO_IDEAL.hd")
308 val emptyTrie
= MT(NONE
,[])
309 fun mkEmpty () = MI(ref (0,emptyTrie
))
311 (* unused code unless searchDeg is used
312 fun maxDeg (MI(x
)) = #
1(!x
)
316 (* unused code unless decode is used
323 fun encode (var
,pwr
) = lshift(var
,16)+pwr
325 fun decode vp
= (rshift(vp
,16),andb(vp
,65535))
327 fun grabVar vp
= andb(vp
,~
65536)
328 fun grabPwr vp
= andb(vp
,65535)
329 fun smallerVar (vp
,vp
') = vp
< andb(vp
',~
65536)
332 fun search (MI(x
),M
.M m
') = let
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
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
))
350 in s(rev(m
',[]),[],mt
)
351 handle Found (* (m
,a
) => SOME(m
,a
) *) => !result
354 (* assume m is a new generator
, i
.e
. not a multiple
of an existing one
*)
355 fun insert (MI (mi
),m
,a
) = let
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
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
)
372 in mi
:= (Int.max(d
,M
.deg m
),i (rev(map
encode(M
.explode m
),[]),mt
)) end
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
383 fun sort i
= if i
>=n
then mi
else let
384 fun redundant (m
,_
) = case search(mi
,m
) of NONE
=> false
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
)
389 in filter(sub1(buckets
,i
),[]);
390 update1(buckets
,i
,[]);
395 fun fold
g (MI(x
)) init
= let
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
)
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
[]
412 end (* structure MI
*)
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
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]
423 fun resetCounts() = app(fn i
=> app (fn j
=> update1(sub1(counts
,i
),j
,0)) indices
) indices
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
433 map (fn i
=> map (fn j
=> sub1(sub1(counts
,i
),j
)) indices
) indices
438 datatype poly
= P
of (F
.field
*M
.mono
) list (* descending mono order
*)
440 fun show (P x
) = (print
"[ ";
442 (print
"("; F
.show f
; print
","; M
.show m
; print
") ")) x
;
446 (* unused code unless power is used
447 val one
= P
[(F
.one
,M
.one
)]
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
)]
454 fun coerce (a
,m
) = P
[(a
,m
)]
456 fun cons (am
,P p
) = P (am
::p
)
459 fun neg p
= (map (fn (a
,m
) => (F
.negate a
,m
)) p
)
460 fun plus ([],p2
) = p2
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
)
469 fun minus ([],p2
) = neg p2
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
)
478 fun termMult (a
,m
,p
) =
479 (map (fn (a
',m
') => (F
.multiply(a
,a
'),M
.multiply(m
,m
'))) p
)
482 fun negate (P p
) = P (neg p
)
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
)))
487 (* unused code unless power is used
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
))
497 fun singleReduce (P y
,a
,m
,P x
) =
498 (pair(length y
,length x
); P(minus(y
,termMult(a
,m
,x
))))
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
))
506 (* unused code unless power is used
516 fun scalarMult (a
,P p
) = P (map (fn (b
,m
) => (F
.multiply(a
,b
),m
)) p
)
521 if k
<=3 then case k
of
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)))
531 fun isZero (P
[]) = true |
isZero (P (_
::_
)) = false
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
541 in fn (P p
,P q
) => eq (p
,q
) end
544 (* these should only be called
if there is a leading term
, i
.e
. poly
<>0 *)
546 fun leadTerm (P(am
::_
)) = am
547 |
leadTerm (P
[]) = Util
.illegal
"POLY.leadTerm"
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"
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
562 (* only used
if r is used
563 fun display (P
[]) = F
.display F
.zero
564 |
display (P p
) = let
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
576 structure HP
= struct
577 datatype hpoly
= HP
of P
.poly array1
579 fun log(n
,l
) = if n
<8 then l
else log((n
>> 2),1+l
)
580 in fn n
=> log(n
,0) end
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
589 fn i
=> if i
<n
then sub1(ps
,i
)
590 else if i
=l
then p
else P
.zero
))
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
))
598 fun leadAndRest (HP ps
) = let
600 fun lar (m
,indices
,i
) = if i
>=n
then lar
'(m
,indices
) else let
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)
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
))
615 in if F
.isZero a
then lar(M
.one
,[],0) else SOME(a
,m
,HP ps
)
617 in lar(M
.one
,[],0) end
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
628 fun reset () = (maybePairs
:=0; primePairs
:=0; usedPairs
:=0; newGens
:=0)
630 fun inc r
= r
:= !r
+ 1
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
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
642 fun mkMonic f
= P
.scalarMult(F
.reciprocal(P
.leadCoeff f
),f
)
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
)
654 fun addPairs (h
,mi
,pairs
) = let
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
)
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
[]))
666 fun ins (~
1,pairs
) = pairs
667 |
ins (i
,pairs
) = case sub1(buckets
,i
) of
669 | gs
=> ins(i
-1,Util
.insert(arrayoflist(h
::gs
),i
,pairs
))
670 in ins(length1 buckets
- 1,pairs
) end
673 fun pr l
= print (String.concat (l@
["\n"]))
674 val fs
= revfold (fn (f
,fs
) => Util
.insert(f
,P
.deg f
,fs
))
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
))
684 val ((_
,m
),rh
) = P
.leadAndRest h
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
))
695 MI
.insert (mi
,P
.leadMono h
,rrh
);
696 app (fn f
=> f
:=autoReduce(!f
)) (!newDegGens
);
697 newDegGens
:= rrh
:: !newDegGens
699 val tasksleft
= ref
0
700 fun feedback () = let
703 if (n
&& 15)=0 then print (Int.toString n
) else ();
705 TextIO.flushOut
TextIO.stdOut
;
715 else let val h
= mkMonic h
716 val _
= (print
"#"; TextIO.flushOut
TextIO.stdOut
)
717 in pairs
:= addPairs(h
,mi
,!pairs
);
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
))
729 try (P
.spair(b
,M
.divide(k
,m
),f
,a
,M
.divide(k
,n
),g
));
732 in tryPair (length1 fgs
-1) end
734 fun numPairs ([],n
) = n
735 |
numPairs (p
::ps
,n
) = numPairs(ps
,n
-1+length1 p
)
737 fun gb d
= if d
>=length1(!pairs
) then mi
else
738 (* note
: i nullify entries to reclaim space
*)
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 ()
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
)]
767 var
::= a |
... | z | A |
... | Z
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
774 datatype char
= Dig
of int | Var
of int | Sign
of int
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
)))
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
)
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
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
803 fun parsePoly s
= poly (P
.zero
,map
char(String.explode s
))
806 fun readIdeal stream
= let
807 fun readLine () = let
808 val s
= input_line stream
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
815 | _
=> map
char (String.explode(substring(s
,i
,n
-i
)))
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 #"
823 val _
= 1=num() orelse Util
.illegal
"stream doesn't start w/ `1'"
826 val _
= length i
= n
orelse Util
.illegal
"wrong # poly's"
831 fun read filename
= let
832 val stream
= open_in filename
833 val i
= readIdeal stream
834 val _
= close_in stream
839 end (* structure G
*)
843 val _
= G
.maxDeg
:=1000000
845 fun grab mi
= MI
.fold (fn ((m
,g
),l
) => P
.cons((F
.one
,m
),!g
)::l
) mi
[]
849 val p
= G
.parsePoly s
850 in print (P
.display p
); print
"\n";
851 print (P
.display(G
.reduce(p
,mi
))); print
"\n"
855 (* unused code unless printCounts is used
856 fun p6 i
= let val s
= Int.toString (i
:int)
858 in print(substring(" ",0,6-n
)); print s
end
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
871 fun maxCount () = revfold (fn (l
,m
) => revfold
Int.max l m
) (getCounts()) 0
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
)
880 (* Unused code unless
sort (analyze
) is used
882 val a
= 16807.0 and m
= 2147483647.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
)
892 (* Unused code unless analyze is used
896 val a
= arrayoflist a
897 val b
= tabulate(length1 a
,fn i
=> i
)
898 val sub
= sub1
and update
= update1
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
914 (* Unused code unless analyze is used
915 fun sum f l
= revfold
op + (map f l
) 0
918 (* Unused code included
in the benchmark
920 val aa
= revfold terms gb
[]
921 val bb
= map M
.explode 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
)
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
) [])))
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
}
938 val g
= G
.grobner fs
handle (Util
.Illegal s
) => (print s
; raise Div
)
940 fun info f
= app print
941 [M
.display(P
.leadMono f
),
942 " + ", Int.toString(P
.numTerms f
- 1), " terms\n"]
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
)
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"]
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"]
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"]
980 * val u4
= map G
.parsePoly
["abcd-e4","a+b+c+d","ab+bc+cd+da","abc+bcd+cda+dab"]
986 * val _
= (print
"Enter fs, u7, u6, u5, or u4: ";
987 * TextIO.flushOut
TextIO.stdOut
)
988 * val s
= TextIO.inputN(TextIO.stdIn
,2)
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"))
995 * gb data
handle e
=> report e
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"]
1012 else (gb u6
; loop (n
- 1))