Menu
▾
▴
[MLton-commit] r4986
|
From: Wesley T. <we...@ml...> - 2006-12-18 18:59:20
|
import didn't work, try adding and commiting
----------------------------------------------------------------------
A mltonlib/trunk/ca/terpstra/math/README
A mltonlib/trunk/ca/terpstra/math/algebra.fun
A mltonlib/trunk/ca/terpstra/math/algebra.sig
A mltonlib/trunk/ca/terpstra/math/c.fun
A mltonlib/trunk/ca/terpstra/math/c.sml
A mltonlib/trunk/ca/terpstra/math/factor.fun
A mltonlib/trunk/ca/terpstra/math/factor.sml
A mltonlib/trunk/ca/terpstra/math/galois.fun
A mltonlib/trunk/ca/terpstra/math/galois.sml
A mltonlib/trunk/ca/terpstra/math/gcd.fun
A mltonlib/trunk/ca/terpstra/math/groups.fun
A mltonlib/trunk/ca/terpstra/math/groups.sig
A mltonlib/trunk/ca/terpstra/math/log.fun
A mltonlib/trunk/ca/terpstra/math/math.mlb
A mltonlib/trunk/ca/terpstra/math/mersenne.fun
A mltonlib/trunk/ca/terpstra/math/mersenne.sml
A mltonlib/trunk/ca/terpstra/math/ops.sig
A mltonlib/trunk/ca/terpstra/math/ops.sml
A mltonlib/trunk/ca/terpstra/math/order.sml
A mltonlib/trunk/ca/terpstra/math/permutation.sml
A mltonlib/trunk/ca/terpstra/math/polynomial.fun
A mltonlib/trunk/ca/terpstra/math/q.fun
A mltonlib/trunk/ca/terpstra/math/q.sml
A mltonlib/trunk/ca/terpstra/math/r.fun
A mltonlib/trunk/ca/terpstra/math/r.sml
A mltonlib/trunk/ca/terpstra/math/rings.fun
A mltonlib/trunk/ca/terpstra/math/rings.sig
A mltonlib/trunk/ca/terpstra/math/test/
A mltonlib/trunk/ca/terpstra/math/test/test.mlb
A mltonlib/trunk/ca/terpstra/math/test/test.sml
A mltonlib/trunk/ca/terpstra/math/test/test2.sml
A mltonlib/trunk/ca/terpstra/math/test/test3.mlb
A mltonlib/trunk/ca/terpstra/math/test/test3.sml
A mltonlib/trunk/ca/terpstra/math/test/test4.mlb
A mltonlib/trunk/ca/terpstra/math/test/test4.sml
A mltonlib/trunk/ca/terpstra/math/z.fun
A mltonlib/trunk/ca/terpstra/math/z.sml
----------------------------------------------------------------------
Added: mltonlib/trunk/ca/terpstra/math/README
===================================================================
--- mltonlib/trunk/ca/terpstra/math/README 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/README 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,44 @@
+This is an incomplete library for number theory.
+
+What there is:
+ small-size galois fields
+ mersenne fields
+ polynomials with Karatsuba (no FFT)
+ permutations
+ generic methods (gcd, discrete logarithm, ...)
+ a reasonably fast integer factorization algorithm
+
+Basically, it's a framework in which I dumped any algorithms or mathematical
+structures I needed as I needed them.
+
+Using it is pretty nice, since you can bind any structure to an operation.
+Generally, one combines functors to create the desired mathematical object.
+
+Some objects already exist (Z = integers, Galois8 = GF(8))
+Some need to be created, eg: ComplexOfField(FieldOfReal(LargeReal))
+
+Structures contain their mathematical operations and substructures.
+The structure Z includes substructures Z.Addition (an abelian group) and
+Z.Multiplication (an abelian monoid). These include relevant operations.
+
+Once a structure has been created, it can be bound to an operation.
+For example:
+ structure Binding = AbelianGroupAddPercent(Z.Addition)
+ open Binding
+will bind =%, +%, -%, ~%, and ++% operations for manipulating Z elements.
+
+To bind the entire structure of Z, one can use:
+ structure Binding = EuclideanDomainDollar(Z)
+ open Binding
+This binds operations: =$, !=$, *$, **$, /$, %$, //$, +$, -$, ~$
+
+There are also some generic algorithms, for example:
+ structure G = GCD(Z)
+This gives you a GCD method that works over integers.
+
+Another example:
+ structure P = PolynomialOverField(ComplexOfField(FieldOfReal(LargeReal)))
+ structure B = Polynomial(P)
+ open B
+Here, the operations with % on the end operate on complex values, while $
+operates on polynomials over the complex numbers.
Added: mltonlib/trunk/ca/terpstra/math/algebra.fun
===================================================================
--- mltonlib/trunk/ca/terpstra/math/algebra.fun 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/algebra.fun 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,69 @@
+(****************************************************************** Algebra ops *)
+
+functor ScalarMultiply(S : SCALAR_MULTIPLY) =
+ struct
+ val (op *&) = S.MUL
+ end
+
+functor Module(M : MODULE) =
+ struct
+ local structure S = RingPercent(M.Base) in open S end
+ local structure S = AbelianGroupAddDollar(M.Addition) in open S end
+ local structure S = ScalarMultiply(M.ScalarMultiplication) in open S end
+ end
+
+functor UnitaryModule(U : UNITARY_MODULE) =
+ struct
+ local structure S = UnitaryRingPercent(U.Base) in open S end
+ local structure S = AbelianGroupAddDollar(U.Addition) in open S end
+ local structure S = ScalarMultiply(U.ScalarMultiplication) in open S end
+ end
+
+functor VectorSpace(V : VECTOR_SPACE) =
+ struct
+ local structure S = FieldPercent(V.Base) in open S end
+ local structure S = AbelianGroupAddDollar(V.Addition) in open S end
+ local structure S = ScalarMultiply(V.ScalarMultiplication) in open S end
+ end
+
+functor Algebra(A : ALGEBRA) =
+ struct
+ local structure S = CommutativeRingPercent(A.Base) in open S end
+ local structure S = NonAssociativeRingDollar(A) in open S end
+ local structure S = ScalarMultiply(A.ScalarMultiplication) in open S end
+ end
+
+functor FieldAlgebra(A : FIELD_ALGEBRA) =
+ struct
+ local structure S = FieldPercent(A.Base) in open S end
+ local structure S = NonAssociativeRingDollar(A) in open S end
+ local structure S = ScalarMultiply(A.ScalarMultiplication) in open S end
+ end
+
+functor AssociativeAlgebra(A : ASSOCIATIVE_ALGEBRA) =
+ struct
+ local structure S = CommutativeRingPercent(A.Base) in open S end
+ local structure S = RingDollar(A) in open S end
+ local structure S = ScalarMultiply(A.ScalarMultiplication) in open S end
+ end
+
+functor AssociativeFieldAlgebra(A : ASSOCIATIVE_FIELD_ALGEBRA) =
+ struct
+ local structure S = FieldPercent(A.Base) in open S end
+ local structure S = RingDollar(A) in open S end
+ local structure S = ScalarMultiply(A.ScalarMultiplication) in open S end
+ end
+
+functor UnitaryAssociativeAlgebra(A : UNITARY_ASSOCIATIVE_ALGEBRA) =
+ struct
+ local structure S = CommutativeRingPercent(A.Base) in open S end
+ local structure S = UnitaryRingDollar(A) in open S end
+ local structure S = ScalarMultiply(A.ScalarMultiplication) in open S end
+ end
+
+functor UnitaryAssociativeFieldAlgebra(A : UNITARY_ASSOCIATIVE_FIELD_ALGEBRA) =
+ struct
+ local structure S = FieldPercent(A.Base) in open S end
+ local structure S = UnitaryRingDollar(A) in open S end
+ local structure S = ScalarMultiply(A.ScalarMultiplication) in open S end
+ end
Added: mltonlib/trunk/ca/terpstra/math/algebra.sig
===================================================================
--- mltonlib/trunk/ca/terpstra/math/algebra.sig 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/algebra.sig 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,103 @@
+signature SCALAR_MULTIPLY =
+ sig
+ type e
+ type t
+
+ (* if one exists: one*v = v*one = v *)
+ val MUL: e * t -> t
+
+ (* a*(b*v) = (a*.b)*v *)
+ val associative : unit
+
+ (* c*(v+w) = c*v + v*w
+ * (c+.d)*v = c*v + d*v
+ *)
+ val distributive : unit
+ end
+
+(* All MODULEs here are left modules *)
+signature MODULE =
+ sig
+ type e
+ type t
+ structure Base : RING where type t = e
+ structure Addition : ABELIAN_GROUP where type t = t
+ structure ScalarMultiplication : SCALAR_MULTIPLY where type e = e and type t = t
+ end
+
+signature UNITARY_MODULE =
+ sig
+ type e
+ type t
+ structure Base : UNITARY_RING where type t = e
+ structure Addition : ABELIAN_GROUP where type t = t
+ structure ScalarMultiplication : SCALAR_MULTIPLY where type e = e and type t = t
+ end
+
+(* same as a UNITARY_MODULE, but over a field *)
+signature VECTOR_SPACE =
+ sig
+ type e
+ type t
+ structure Base : FIELD where type t = e
+ structure Addition : ABELIAN_GROUP where type t = t
+ structure ScalarMultiplication : SCALAR_MULTIPLY where type e = e and type t = t
+ end
+
+signature ALGEBRA =
+ sig
+ include NON_ASSOCIATIVE_RING
+ type e
+ structure Base : COMMUTATIVE_RING where type t = e
+ structure ScalarMultiplication : SCALAR_MULTIPLY where type e = e and type t = t
+
+ (* (ax)y = a(xy)
+ * a(xy) = x(ay)
+ *)
+ val bilinear : unit
+ end
+
+signature FIELD_ALGEBRA =
+ sig
+ include NON_ASSOCIATIVE_RING
+ type e
+ structure Base : FIELD where type t = e
+ structure ScalarMultiplication : SCALAR_MULTIPLY where type e = e and type t = t
+ val bilinear : unit
+ end
+
+signature ASSOCIATIVE_ALGEBRA =
+ sig
+ include RING
+ type e
+ structure Base : COMMUTATIVE_RING where type t = e
+ structure ScalarMultiplication : SCALAR_MULTIPLY where type e = e and type t = t
+ val bilinear : unit
+ end
+
+signature ASSOCIATIVE_FIELD_ALGEBRA =
+ sig
+ include RING
+ type e
+ structure Base : FIELD where type t = e
+ structure ScalarMultiplication : SCALAR_MULTIPLY where type e = e and type t = t
+ val bilinear : unit
+ end
+
+signature UNITARY_ASSOCIATIVE_ALGEBRA =
+ sig
+ include UNITARY_RING
+ type e
+ structure Base : COMMUTATIVE_RING where type t = e
+ structure ScalarMultiplication : SCALAR_MULTIPLY where type e = e and type t = t
+ val bilinear : unit
+ end
+
+signature UNITARY_ASSOCIATIVE_FIELD_ALGEBRA =
+ sig
+ include UNITARY_RING
+ type e
+ structure Base : FIELD where type t = e
+ structure ScalarMultiplication : SCALAR_MULTIPLY where type e = e and type t = t
+ val bilinear : unit
+ end
Added: mltonlib/trunk/ca/terpstra/math/c.fun
===================================================================
--- mltonlib/trunk/ca/terpstra/math/c.fun 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/c.fun 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,47 @@
+(* This is still a field only if (x^2 + 1) is irreducible *)
+functor ComplexOfField(F : FIELD) : FIELD =
+ struct
+ local
+ structure B = FieldPercent(F)
+ structure O = EuclideanDomainDollar(Order)
+ open B
+ open O
+ in
+ type t = F.t * F.t
+ val characteristic = F.characteristic
+
+ structure Addition =
+ struct
+ type t = F.t * F.t
+ val order = F.Addition.order *$ F.Addition.order
+
+ val associative = ()
+ val commutative = ()
+ val one = (#%0, #%0)
+
+ val EQ = fn ((ar, ai), (br, bi)) => (ar =% br) andalso (ai =% bi)
+ val MUL = fn ((ar, ai), (br, bi)) => (ar +% br, ai +% bi)
+ val DIV = fn ((ar, ai), (br, bi)) => (ar -% br, ai -% bi)
+ val INV = fn (a, b) => (~%a, ~%b)
+ end
+
+ structure Multiplication =
+ struct
+ type t = F.t * F.t
+ val order = F.Addition.order *$ F.Addition.order -$ #$1
+
+ val associative = ()
+ val commutative = ()
+ val one = (#%1, #%0)
+
+ val EQ = fn ((ar, ai), (br, bi)) => (ar =% br) andalso (ai =% bi)
+ val INV = fn (r, i) =>
+ let val e = !%(r*%r +% i*%i) in (r*%e, ~%i*%e) end
+ val MUL = fn ((ar, ai), (br, bi)) => (ar*%br -% ai*%bi, ar*%bi +% ai*%br)
+ val DIV = fn (a, c) => MUL (a, INV c)
+ end
+
+ val distributive = ()
+ val no_zero_divisors = ()
+ end
+ end
Added: mltonlib/trunk/ca/terpstra/math/c.sml
===================================================================
--- mltonlib/trunk/ca/terpstra/math/c.sml 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/c.sml 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,4 @@
+structure C = ComplexOfField(R)
+
+structure C32 = ComplexOfField(R32)
+structure C64 = ComplexOfField(R64)
Added: mltonlib/trunk/ca/terpstra/math/factor.fun
===================================================================
--- mltonlib/trunk/ca/terpstra/math/factor.fun 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/factor.fun 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,22 @@
+functor Factor(E : EUCLIDEAN_DOMAIN) =
+ struct
+ local
+ structure B = EuclideanDomainPercent(E)
+ open B
+
+ fun push (z, l) =
+ case l of
+ nil => (z, #%1) :: nil
+ | (y, e) :: l =>
+ if y =% z
+ then (y, e +% #%1) :: l
+ else (z, #%1) :: (y, e) :: l
+
+ fun niaveFactor (z, l) a =
+ if z <% a*%a then push (z, l) else
+ if z %% a =% #%0 then niaveFactor (z /% a, push (a, l)) a else
+ niaveFactor (z, l) (a +% #%1)
+ in
+ fun factor z = niaveFactor (z, nil) (#%2)
+ end
+ end
Added: mltonlib/trunk/ca/terpstra/math/factor.sml
===================================================================
--- mltonlib/trunk/ca/terpstra/math/factor.sml 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/factor.sml 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,54 @@
+structure Factor =
+ struct
+ local
+ open LargeInt
+ in
+ fun isPrime n =
+ let
+ fun chop (c, e) = if c mod 2 = 0 then chop (c div 2, e+1) else (c, e)
+ val (c, e) = chop (n-1, 0)
+
+ fun exp (x, 0) = 1
+ | exp (x, 1) = x
+ | exp (x, e) =
+ let val y = exp (x*x mod n, e div 2)
+ in if e mod 2 = 0 then y else y * x mod n end
+
+ fun checkPow (w, e) =
+ e > 0 andalso (w+1 = n orelse checkPow (w*w mod n, e-1))
+ fun millerRabin w =
+ let val v = exp (w, c)
+ in v = 1 orelse checkPow (v, e) end
+ in
+ List.foldl (fn (w, a) => a andalso millerRabin (1 + w mod (n-1)))
+ true [62151, 7444, 40814, 49239, 71708759, 4481, 665652934 ]
+ end
+
+ fun factor n =
+ let
+ fun sort nil = nil
+ | sort (p :: r) =
+ let val (s, b) = List.partition (fn x => x < p) r
+ in sort s @ p :: sort b end
+
+ fun gcd (a, b) = if a = 0 then b else gcd (b mod a, a)
+
+ fun findafactor n =
+ let
+ val start = Word.toLargeInt (MLton.Random.rand ()) mod n
+ fun f x = (x*x + 2) mod n
+ fun loop x y =
+ let val g = gcd (n, n+x-y)
+ in if g = 1 then loop (f x) (f (f y)) else g end
+ in loop start (f start) end
+
+ fun suck n l =
+ if n = 1 then l else
+ if isPrime n then n :: l else
+ let val x = findafactor n
+ in suck x (suck (n div x) l) end
+ in
+ sort (suck n [])
+ end
+ end
+ end
Added: mltonlib/trunk/ca/terpstra/math/galois.fun
===================================================================
--- mltonlib/trunk/ca/terpstra/math/galois.fun 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/galois.fun 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,93 @@
+signature GaloisParam =
+ sig
+ structure W : WORD
+ val base : W.word
+ end
+
+functor GaloisFromTable(P : GaloisParam) : FIELD =
+ struct
+ local
+ open P.W
+ open Int
+
+ val fromInt = P.W.fromInt
+ val toInt = P.W.toInt
+
+ val zero = fromInt 0
+ val one = fromInt 1
+
+ val msize = toInt (notb zero)
+ val size = msize + 1
+ val highbit = << (one, Word.-(Word.fromInt wordSize, 0w1))
+
+ fun nastymul x y =
+ if y = zero then zero else
+ if y = one then x else
+ let val h = nastymul x (>> (y, 0w1))
+ val x = if andb (y, one) = one then x else zero
+ val b = if andb (h, highbit) = highbit then P.base else zero
+ in xorb( << (h, 0w1), xorb(x, b)) end
+
+ val gen = orb (one, << (one, 0w1)) (* x^1 + x^0 <- a generator *)
+
+ val log = Array.array (size, zero)
+ fun set _ l ~1 = l
+ | set a l i = (
+ Array.update (log, toInt a, fromInt i)
+ ; set (nastymul a gen) (a :: l) (i - 1))
+
+ val exp = Vector.fromList (set gen nil (msize + msize - 1))
+ val log = Array.vector log
+(*
+ val _ = Vector.app (fn x => print (P.W.toString x ^ " ")) exp
+ val _ = print "\n"
+ val _ = Vector.app (fn x => print (P.W.toString x ^ " ")) log
+ val _ = print "\n"
+*)
+ val exp = fn x => Vector.sub (exp, x)
+ val log = fn x => toInt (Vector.sub (log, toInt x))
+ in
+ type t = word
+ val characteristic = LargeInt.fromInt 2
+
+ structure Addition =
+ struct
+ type t = word
+ val order = FINITE (Int.toLarge msize)
+
+ val associative = ()
+ val commutative = ()
+ val one = fromInt 0
+
+ val EQ = (op =)
+ val MUL = fn (x, y) => xorb (x, y)
+ val DIV = MUL
+ val INV = fn x => x
+ end
+
+ structure Multiplication =
+ struct
+ type t = word
+ val order = FINITE (Int.toLarge size)
+
+ val associative = ()
+ val commutative = ()
+ val one = fromInt 1
+
+ val EQ = (op =)
+ val MUL = fn (x, y) =>
+ if x = zero orelse y = zero then zero else
+ let val (lx, ly) = (log x, log y)
+ in exp (lx + ly) end
+ val DIV = fn (x, y) =>
+ let val (lx, ly) = (log x, log y)
+ in exp (lx + msize - ly) end
+ val INV = fn x =>
+ let val lx = log x
+ in exp (msize - lx) end
+ end
+
+ val distributive = ()
+ val no_zero_divisors = ()
+ end
+ end
Added: mltonlib/trunk/ca/terpstra/math/galois.sml
===================================================================
--- mltonlib/trunk/ca/terpstra/math/galois.sml 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/galois.sml 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,21 @@
+local
+ structure P4 =
+ struct
+ structure W = Word4
+ val base : W.word = 0wx3
+ end
+ structure P5 =
+ struct
+ structure W = Word5
+ val base : W.word = 0wx9
+ end
+ structure P8 =
+ struct
+ structure W = Word8
+ val base : W.word = 0wx1B
+ end
+in
+ structure Galois4 = GaloisFromTable(P4)
+ structure Galois5 = GaloisFromTable(P5)
+ structure Galois8 = GaloisFromTable(P8)
+end
Added: mltonlib/trunk/ca/terpstra/math/gcd.fun
===================================================================
--- mltonlib/trunk/ca/terpstra/math/gcd.fun 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/gcd.fun 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,49 @@
+exception NoSolution
+
+functor GCD(E : EUCLIDEAN_DOMAIN) =
+ struct
+ local
+ structure B = EuclideanDomainPercent(E)
+ open B
+ in
+ (* (g, i, j, x) = egcd(a, b)
+ * gcd(a,b) = g = a*i + b*j
+ * x = -1 -> i <= 0 and j >= 0
+ * x = +1 -> j <= 0 and i >= 0
+ *)
+ fun egcd (a, b) =
+ if (b =% #%0) then (a, #%1, #%0, #%1) else
+ let
+ val (q, r) = a //% b
+ val (g, i, j, x) = egcd (b, r)
+ in
+ (g, j, i -% q*%j, ~%x)
+ end
+
+ fun gcd (a, b) = case egcd (a, b) of (g, _, _, _) => g
+ fun lcm (a, b) = a*%b /% gcd(a, b)
+
+ fun gcdinv (n, m) =
+ let
+ val (g, ni, mi, x) = egcd(n, m)
+ in
+ if x =% #%1
+ then (g, ni, n+%mi)
+ else (g, m+%ni, mi)
+ end
+
+ fun inv (a, n) = case gcdinv (n, a) of (_, _, ai) => ai
+
+ fun crt nil = (#%0, #%1)
+ | crt ((a, n) :: x) =
+ let
+ val (b, m) = crt x
+ val (g, ni, mi) = gcdinv (n, m)
+ val l = n*%m/%g
+ in
+ if a %% g =% b %% g
+ then ((a*%m*%mi +% b*%n*%ni)/%g %% l, l)
+ else raise NoSolution
+ end
+ end
+ end
Added: mltonlib/trunk/ca/terpstra/math/groups.fun
===================================================================
--- mltonlib/trunk/ca/terpstra/math/groups.fun 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/groups.fun 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,149 @@
+(****************************************************************** Group ops *)
+
+exception ImpossibleConstant
+
+functor SetPercent(S : BINARY_OPERATION) =
+ struct
+ val (op =% ) = S.EQ
+ val (op !=%) = not o S.EQ
+ end
+
+functor SetDollar(S : BINARY_OPERATION) =
+ struct
+ val (op =$ ) = S.EQ
+ val (op !=$) = not o S.EQ
+ end
+
+functor BinaryOperationMulPercent(B : BINARY_OPERATION) =
+ struct
+ local structure S = SetPercent(B) in open S end
+ val (op *% ) = B.MUL
+ local open LargeInt in
+ fun _ **% 0 = raise ImpossibleConstant
+ | x **% 1 = x
+ | x **% e =
+ let val h = x **% (e div 2)
+ in if (e mod 2) = 0 then h *% h else h *% h *% x end
+ end
+ end
+functor BinaryOperationAddPercent(B : BINARY_OPERATION) =
+ struct
+ local structure S = SetPercent(B) in open S end
+ val (op +% ) = B.MUL
+ local open LargeInt in
+ fun _ ++% 0 = raise ImpossibleConstant
+ | x ++% 1 = x
+ | x ++% e =
+ let val h = x ++% (e div 2)
+ in if (e mod 2) = 0 then h +% h else h +% h +% x end
+ end
+ end
+
+functor BinaryOperationMulDollar(B : BINARY_OPERATION) =
+ struct
+ local structure S = SetDollar(B) in open S end
+ val (op *$ ) = B.MUL
+ local open LargeInt in
+ fun _ **$ 0 = raise ImpossibleConstant
+ | x **$ 1 = x
+ | x **$ e =
+ let val h = x **$ (e div 2)
+ in if (e mod 2) = 0 then h *$ h else h *$ h *$ x end
+ end
+ end
+functor BinaryOperationAddDollar(B : BINARY_OPERATION) =
+ struct
+ local structure S = SetDollar(B) in open S end
+ val (op +$ ) = B.MUL
+ local open LargeInt in
+ fun _ ++$ 0 = raise ImpossibleConstant
+ | x ++$ 1 = x
+ | x ++$ e =
+ let val h = x ++$ (e div 2)
+ in if (e mod 2) = 0 then h +$ h else h +$ h +$ x end
+ end
+ end
+
+functor SemiGroupMulPercent(S : SEMIGROUP) = BinaryOperationMulPercent(S)
+functor SemiGroupAddPercent(S : SEMIGROUP) = BinaryOperationAddPercent(S)
+functor SemiGroupMulDollar (S : SEMIGROUP) = BinaryOperationMulDollar (S)
+functor SemiGroupAddDollar (S : SEMIGROUP) = BinaryOperationAddDollar (S)
+
+functor MonoidMulPercent(M : MONOID) =
+ struct
+ local structure S = SemiGroupMulPercent(M) in open S end
+ val (op **%) = fn (_, 0) => M.one | (x, e) => x **% e
+ val #% = fn 1 => M.one | _ => raise ImpossibleConstant
+ end
+functor MonoidAddPercent(M : MONOID) =
+ struct
+ local structure S = SemiGroupAddPercent(M) in open S end
+ val (op ++%) = fn (_, 0) => M.one | (x, e) => x ++% e
+ val #% = fn 0 => M.one | _ => raise ImpossibleConstant
+ end
+
+functor MonoidMulDollar(M : MONOID) =
+ struct
+ local structure S = SemiGroupMulDollar(M) in open S end
+ val (op **$) = fn (_, 0) => M.one | (x, e) => x **$ e
+ val #$ = fn 1 => M.one | _ => raise ImpossibleConstant
+ end
+functor MonoidAddDollar(M : MONOID) =
+ struct
+ local structure S = SemiGroupAddDollar(M) in open S end
+ val (op ++$) = fn (_, 0) => M.one | (x, e) => x ++$ e
+ val #$ = fn 0 => M.one | _ => raise ImpossibleConstant
+ end
+
+functor AbelianMonoidMulPercent(A : ABELIAN_MONOID) = MonoidMulPercent(A)
+functor AbelianMonoidAddPercent(A : ABELIAN_MONOID) = MonoidAddPercent(A)
+functor AbelianMonoidMulDollar (A : ABELIAN_MONOID) = MonoidMulDollar (A)
+functor AbelianMonoidAddDollar (A : ABELIAN_MONOID) = MonoidAddDollar (A)
+
+functor GroupMulPercent(G : GROUP) =
+ struct
+ local structure S = MonoidMulPercent(G) in open S end
+ val (op !% ) = G.INV
+ val (op **%) = fn (x, e) => if e < 0 then !%x **% ~e else x **% e
+ end
+functor GroupAddPercent(G : GROUP) =
+ struct
+ local structure S = MonoidAddPercent(G) in open S end
+ val (op ~% ) = G.INV
+ val (op ++%) = fn (x, e) => if e < 0 then ~%x ++% ~e else x ++% e
+ end
+
+functor GroupMulDollar(G : GROUP) =
+ struct
+ local structure S = MonoidMulDollar(G) in open S end
+ val (op !$ ) = G.INV
+ val (op **$) = fn (x, e) => if e < 0 then !$x **$ ~e else x **$ e
+ end
+functor GroupAddDollar(G : GROUP) =
+ struct
+ local structure S = MonoidAddDollar(G) in open S end
+ val (op ~$ ) = G.INV
+ val (op ++$) = fn (x, e) => if e < 0 then ~$x ++$ ~e else x ++$ e
+ end
+
+functor AbelianGroupMulPercent(A : ABELIAN_GROUP) =
+ struct
+ local structure S = GroupMulPercent(A) in open S end
+ val (op /%) = A.DIV
+ end
+functor AbelianGroupAddPercent(A : ABELIAN_GROUP) =
+ struct
+ local structure S = GroupAddPercent(A) in open S end
+ val (op -%) = A.DIV
+ end
+
+functor AbelianGroupMulDollar(A : ABELIAN_GROUP) =
+ struct
+ local structure S = GroupMulDollar(A) in open S end
+ val (op /$) = A.DIV
+ end
+functor AbelianGroupAddDollar(A : ABELIAN_GROUP) =
+ struct
+ local structure S = GroupAddDollar(A) in open S end
+ val (op -$) = A.DIV
+ end
Added: mltonlib/trunk/ca/terpstra/math/groups.sig
===================================================================
--- mltonlib/trunk/ca/terpstra/math/groups.sig 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/groups.sig 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,65 @@
+signature SEMIGROUP =
+ sig
+ include BINARY_OPERATION
+ (* a*(b*c) = (a*b)*c *)
+ val associative : unit
+ end
+
+signature MONOID =
+ sig
+ include SEMIGROUP
+
+ (* one*x = x*one = x *)
+ val one : t
+ end
+
+signature ABELIAN_MONOID =
+ sig
+ include MONOID
+
+ (* a*b = b*a *)
+ val commutative : unit
+ end
+
+signature GROUP =
+ sig
+ include MONOID
+ (* y * (y^-1) = (y^-1) * y = one *)
+ val INV : t -> t
+ end
+
+signature PERMUTATION =
+ sig
+ include GROUP
+ (* include ENDOFUNCTION *)
+ type v
+ val EVAL: t -> v -> v
+ end
+signature AUTOFUNCTION = PERMUTATION
+
+signature ABELIAN_GROUP =
+ sig
+ include GROUP
+
+ (* a*b = b*a *)
+ val commutative : unit
+
+ (* x/y = x*(y^-1)) -- could be defined in GROUP, but is confusing *)
+ val DIV : (t * t) -> t
+ end
+
+signature CYCLIC_GROUP =
+ sig
+ include ABELIAN_GROUP
+ (* x = g^i *)
+ val generator : t
+ end
+
+(****************************************************************** Unity *)
+
+signature FINITE_SUBGROUPS =
+ sig
+ type t
+ (* !!! Somehow figure out how to grab unity roots *)
+ end
+
Added: mltonlib/trunk/ca/terpstra/math/log.fun
===================================================================
--- mltonlib/trunk/ca/terpstra/math/log.fun 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/log.fun 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,28 @@
+(*
+functor DiscreteLogarithm(G : GROUP) =
+ struct
+ local
+ structure B = GroupMulPercent(G)
+ structure D = GCD(Z)
+ open B D Factor LargeInt
+
+ exception NotDiscrete
+ val order =
+ case G.order of
+ FINITE x => x
+ | _ => raise NotDiscrete
+
+ datatype Factorization =
+ PRIME of LargeInt.int
+ | POWER of Order * LargeInt.int
+ | COPRIME of Order * Order
+ withtype Order = LargeInt.int * Factorization
+
+ val factors = factor order
+
+ fun log (n, f)
+ in
+ fun log = crtLog
+ end
+ end
+*)
Added: mltonlib/trunk/ca/terpstra/math/math.mlb
===================================================================
--- mltonlib/trunk/ca/terpstra/math/math.mlb 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/math.mlb 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,44 @@
+ann
+ "deadCode true"
+ "sequenceNonUnit warn"
+in
+ local
+ $(MLTON_ROOT)/basis/mlton.mlb
+ $(MLTON_ROOT)/basis/basis.mlb
+ in
+ (* base types *)
+ ops.sig
+ ops.sml
+ groups.sig
+ groups.fun
+ rings.sig
+ rings.fun
+ algebra.sig
+ algebra.fun
+
+ polynomial.fun
+
+ (* generic algorithms need for construction *)
+ gcd.fun
+
+ (* concrete constructions *)
+ order.sml
+ z.fun
+ z.sml
+ r.fun
+ r.sml
+ c.fun
+ c.sml
+ q.fun
+ q.sml
+ mersenne.fun
+ mersenne.sml
+ galois.fun
+ galois.sml
+ permutation.sml
+
+ (* other algorithms *)
+ factor.sml
+ log.fun
+ end
+end
Added: mltonlib/trunk/ca/terpstra/math/mersenne.fun
===================================================================
--- mltonlib/trunk/ca/terpstra/math/mersenne.fun 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/mersenne.fun 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,77 @@
+signature MERSENNE_BASE =
+ sig
+ structure Z : WORD (* fits the numbers with 1 bit to spare *)
+ structure ZZ : WORD (* fits the product *)
+ val bits : word
+ end
+
+exception MersenneOverflow
+
+functor Mersenne(M : MERSENNE_BASE) : FIELD =
+ struct
+ local
+ structure Z = M.Z
+ structure ZZ = M.ZZ
+ structure O = EuclideanDomainDollar(Order)
+ structure G = GCD(EuclideanDomainOfWord(Z))
+ open Z
+ open O
+ open G
+
+ val mbits = M.bits
+ val mmbits = Word.+ (mbits, mbits)
+ val rbits = Word.fromInt wordSize
+ val rrbits = Word.fromInt ZZ.wordSize
+ val rbitsm1 = Word.- (rbits, 0w1)
+
+ val _ = if Word.>= ( mbits, rbits) then raise MersenneOverflow else ()
+ val _ = if Word.>= (mmbits, rrbits) then raise MersenneOverflow else ()
+
+ val mmask = << (fromInt 1, mbits) - fromInt 1
+ in
+ type t = word
+ val characteristic = toLargeInt mmask
+
+ structure Addition =
+ struct
+ type t = word
+ val order = FINITE characteristic
+
+ val associative = ()
+ val commutative = ()
+ val one = fromInt 0
+
+ val EQ = (op =)
+ val MUL = fn (x, y) =>
+ let val z = x + y in andb (z, mmask) + >> (z, mbits) end
+ val DIV = fn (x, y) =>
+ let val z = x - y in andb (z, mmask) - >> (z, rbitsm1) end
+ val INV = fn x => mmask - x
+ end
+
+ structure Multiplication =
+ struct
+ type t = word
+ val order = Addition.order -$ #$1
+
+ val associative = ()
+ val commutative = ()
+ val one = fromInt 1
+
+ val EQ = (op =)
+ val MUL = fn (x, y) =>
+ let
+ val ZZ = ZZ.fromLarge o toLarge
+ val Z = fromLarge o ZZ.toLarge
+ val z = ZZ.* (ZZ x, ZZ y)
+ in
+ Addition.MUL (Z (ZZ.>> (z, mbits)), andb (Z z, mmask))
+ end
+ val INV = fn x => inv (x, mmask)
+ val DIV = fn (x, y) => MUL (x, INV y)
+ end
+
+ val distributive = ()
+ val no_zero_divisors = ()
+ end
+ end
Added: mltonlib/trunk/ca/terpstra/math/mersenne.sml
===================================================================
--- mltonlib/trunk/ca/terpstra/math/mersenne.sml 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/mersenne.sml 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,27 @@
+(* all Mersenne primes up to 64 bits: 2, 3, 5, 7, 13, 17, 19, 31, 61 *)
+(* In order to use 61 we need a Word128 on 64bit machines *)
+
+local
+ structure P7 =
+ struct
+ structure Z = Word8
+ structure ZZ = Word16
+ val bits = 0w7
+ end
+ structure P13 =
+ struct
+ structure Z = Word16
+ structure ZZ = Word32
+ val bits = 0w13
+ end
+ structure P31 =
+ struct
+ structure Z = Word32
+ structure ZZ = Word64
+ val bits = 0w31
+ end
+in
+ structure Mersenne7 = Mersenne(P7)
+ structure Mersenne13 = Mersenne(P13)
+ structure Mersenne31 = Mersenne(P31)
+end
Added: mltonlib/trunk/ca/terpstra/math/ops.sig
===================================================================
--- mltonlib/trunk/ca/terpstra/math/ops.sig 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/ops.sig 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,23 @@
+(****************************************************************** Groups *)
+
+datatype order = FINITE of LargeInt.int | COUNTABLE | UNCOUNTABLE
+
+signature SET =
+ sig
+ type t
+ val order: order
+ val EQ: (t * t) -> bool
+ end
+
+signature BINARY_OPERATION =
+ sig
+ include SET
+ val MUL: (t * t) -> t
+ end
+
+signature ENDOFUNCTION =
+ sig
+ include BINARY_OPERATION
+ type v
+ val EVAL: t -> v -> v
+ end
Added: mltonlib/trunk/ca/terpstra/math/ops.sml
===================================================================
--- mltonlib/trunk/ca/terpstra/math/ops.sml 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/ops.sml 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,16 @@
+(****************************************************************** Operators *)
+
+infix 8 ++% **% (* additive and multiplicative exponentiation *)
+infix 7 *% /% %% //% (* //% = (/%, %%) = (div, mod) *)
+infix 6 +% -%
+infix 4 =% !=% <% (* <>% would imply <% which does not always exist *)
+nonfix ~% !% #%
+
+(* the operators in an algebra are always these (above for scalars) *)
+infix 8 ++$ **$
+infix 7 *$ /$ %$ //$
+infix 6 +$ -$
+infix 4 =$ !=$ <$
+nonfix ~$ !$ #$
+
+infix 7 *& (* for scalar operations with a vector *)
Added: mltonlib/trunk/ca/terpstra/math/order.sml
===================================================================
--- mltonlib/trunk/ca/terpstra/math/order.sml 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/order.sml 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,59 @@
+structure Order : EUCLIDEAN_DOMAIN =
+ struct
+ local
+ exception Undefined
+ open LargeInt
+ in
+ type t = order
+ val characteristic = LargeInt.fromInt 0
+
+ structure Addition =
+ struct
+ type t = order
+ val order = UNCOUNTABLE
+
+ val associative = ()
+ val commutative = ()
+ val one = FINITE 0
+
+ val EQ = (op =)
+ fun MUL (FINITE x, FINITE y) = FINITE (x+y)
+ | MUL (UNCOUNTABLE, _) = UNCOUNTABLE
+ | MUL (_, UNCOUNTABLE) = UNCOUNTABLE
+ | MUL _ = COUNTABLE
+ fun INV (FINITE x) = FINITE (~x)
+ | INV _ = raise Undefined
+ fun DIV (x, y) = MUL (x, INV y)
+ end
+
+ structure Multiplication =
+ struct
+ type t = order
+ val order = UNCOUNTABLE
+
+ val associative = ()
+ val commutative = ()
+ val one = FINITE 1
+
+ val EQ = (op =)
+ fun MUL (FINITE x, FINITE y) = FINITE (x*y)
+ | MUL (_, FINITE 0) = raise Undefined
+ | MUL (FINITE 0, _) = raise Undefined
+ | MUL (UNCOUNTABLE, _) = UNCOUNTABLE
+ | MUL (_, UNCOUNTABLE) = UNCOUNTABLE
+ | MUL _ = COUNTABLE
+ end
+
+ val distributive = ()
+ val no_zero_divisors = ()
+
+ fun QUO (FINITE a, FINITE b) = (FINITE (a div b), FINITE (a mod b))
+ | QUO _ = raise Undefined
+
+ fun LT (UNCOUNTABLE, _) = false
+ | LT (_, UNCOUNTABLE) = true
+ | LT (COUNTABLE, _) = false
+ | LT (_, COUNTABLE) = true
+ | LT (FINITE x, FINITE y) = x < y
+ end
+ end
Added: mltonlib/trunk/ca/terpstra/math/permutation.sml
===================================================================
--- mltonlib/trunk/ca/terpstra/math/permutation.sml 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/permutation.sml 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,210 @@
+signature FINITE_PERMUTATION =
+ sig
+ include PERMUTATION
+
+ (* a bijection exists between naturals and permutations *)
+ val toInt: t -> int
+ val fromInt: int -> t
+ val toLargeInt: t -> LargeInt.int
+ val fromLargeInt: LargeInt.int -> t
+
+ (* an alternate representation exists as cycles *)
+ val fromCycle: int list -> t
+ val fromCycles: int list list -> t
+ val toCycles: t -> int list list
+ end
+
+structure FinitePermutation : FINITE_PERMUTATION =
+ struct
+ type v = int
+ type t = int Vector.vector
+ val order = COUNTABLE
+
+ val associative = ()
+ val one = Vector.tabulate (0, fn _ => 0)
+
+ fun EVAL f x =
+ if x >= Vector.length f then x else
+ Vector.sub (f, x)
+
+ fun EQ (x, y) =
+ let
+ val lx = Vector.length x
+ val ly = Vector.length y
+ in
+ if lx < ly
+ then Vector.foldli (fn (i, v, a) => a andalso v = EVAL x i) true y
+ else Vector.foldli (fn (i, v, a) => a andalso v = EVAL y i) true x
+ end
+
+ fun MUL (x, y) =
+ let
+ val lx = Vector.length x
+ val ly = Vector.length y
+ val l = if lx < ly then ly else lx
+ in
+ Vector.tabulate (l, fn i => EVAL x (EVAL y i))
+ end
+
+ fun INV x =
+ let
+ val y = Array.array (Vector.length x, 0)
+ in
+ Vector.appi (fn (i, v) => Array.update (y, v, i)) x;
+ Array.vector y
+ end
+
+ fun toInt p =
+ let
+ fun helper p =
+ let
+ val l = Vector.length p
+ in
+ if l < 2 then (1, 0) else
+ let
+ val i = l - 1
+ val s = VectorSlice.slice (p, 0, SOME i)
+ val x = Vector.sub (p, i)
+ val v = VectorSlice.map (fn y => if y > x then y-1 else y) s
+ val (fact, out) = helper v
+ in
+ (fact * l, out + (i - x) * fact)
+ end
+ end
+ in
+ #2 (helper p)
+ end
+
+ fun toLargeInt p =
+ let
+ fun helper p =
+ let
+ val l = Vector.length p
+ in
+ if l < 2 then (1, 0) else
+ let
+ val i = l - 1
+ val s = VectorSlice.slice (p, 0, SOME i)
+ val x = Vector.sub (p, i)
+ val v = VectorSlice.map (fn y => if y > x then y-1 else y) s
+ val (fact, out) = helper v
+ val y = i - x
+ open LargeInt
+ in
+ (fact * fromInt l, out + fromInt y * fact)
+ end
+ end
+ in
+ #2 (helper p)
+ end
+
+ fun fromInt x =
+ if x = 0 then Vector.tabulate (0, fn _ => 0) else
+ let
+ fun grow (l, f) =
+ if f > x then (l, f) else
+ grow (l+1, f*(l+1))
+ val (l, f) = grow (0, 1)
+ fun helper (1, _, _) = Vector.tabulate (1, fn _ => 0)
+ | helper (l, f, x) =
+ let
+ val (l1, f) = (l - 1, f div l)
+ val (q, r) = (x div f, x mod f)
+ val p = helper (l1, f, r)
+ val z = (l1 - q)
+ fun build i =
+ if i = l1 then z else
+ let val y = Vector.sub (p, i) in
+ if y >= z then y+1 else y end
+ in
+ Vector.tabulate (l, build)
+ end
+ in
+ helper (l, f, x)
+ end
+
+ fun fromLargeInt x =
+ if x = 0 then Vector.tabulate (0, fn _ => 0) else
+ let
+ fun grow (l, f) =
+ if f > x then (l, f) else
+ grow (l + 1, LargeInt.* (f, LargeInt.fromInt (l + 1)))
+ val (l, f) = grow (0, 1)
+ fun helper (1, _, _) = Vector.tabulate (1, fn _ => 0)
+ | helper (l, f, x) =
+ let
+ val (l1, f) = (l - 1, LargeInt.div (f, LargeInt.fromInt l))
+ val (q, r) = IntInf.quotRem (x, f)
+ val p = helper (l1, f, r)
+ val z = (l1 - LargeInt.toInt q)
+ fun build i =
+ if i = l1 then z else
+ let val y = Vector.sub (p, i) in
+ if y >= z then y+1 else y end
+ in
+ Vector.tabulate (l, build)
+ end
+ in
+ helper (l, f, x)
+ end
+
+ fun fromCycle c =
+ let
+ val m = List.foldl (fn (x, a) => if x < a then a else x) 0 c
+ val a = Array.tabulate (m+1, fn i => i)
+ fun helper [] = ()
+ | helper (x :: []) = Array.update (a, x, List.hd c)
+ | helper (x :: y :: r) = (
+ Array.update (a, x, y);
+ helper (y :: r))
+ in
+ helper c;
+ Array.vector a
+ end
+
+ fun fromCycles l =
+ List.foldl (fn (x, a) => MUL (fromCycle x, a)) one l
+
+ fun toCycles p =
+ let
+ val a = Array.tabulate (Vector.length p, fn i => false)
+ fun trace i =
+ if Array.sub (a, i) then [] else
+ (Array.update (a, i, true);
+ i :: trace (Vector.sub (p, i)))
+ fun scan (i, l) =
+ if i = l then [] else
+ case trace i of
+ [] => scan (i+1, l)
+ | _ :: [] => scan (i+1, l)
+ | x => x :: scan (i+1, l)
+ in
+ scan (0, Vector.length p)
+ end
+ end
+
+(*
+fun test (x, e) =
+ if x = e then () else
+ let
+ val p = FinitePermutation.fromLargeInt (LargeInt.fromInt x)
+ val y = FinitePermutation.toLargeInt p
+ in
+ print (LargeInt.toString y ^ ":");
+ Vector.app (fn y => print (" " ^ Int.toString y)) p;
+ print "\n";
+ test (x+1, e)
+ end
+
+val () = test (0, 720)
+
+val p = FinitePermutation.fromCycles [[5, 4], [3], [], [2, 6, 3]]
+val () = Vector.app (fn y => print (" " ^ Int.toString y)) p
+val () = print "\n"
+
+val p = FinitePermutation.fromCycles [[5, 4], [3], [], [2, 6, 3]]
+val cs = FinitePermutation.toCycles p
+val pc = List.app (fn y => print (" " ^ Int.toString y))
+val pcs = List.app (fn y => (pc y; print "\n"))
+val () = pcs cs
+*)
Added: mltonlib/trunk/ca/terpstra/math/polynomial.fun
===================================================================
--- mltonlib/trunk/ca/terpstra/math/polynomial.fun 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/polynomial.fun 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,311 @@
+signature POLYNOMIAL =
+ sig
+ include EUCLIDEAN_DOMAIN
+ (* include UNITARY_ASSOCIATIVE_FIELD_ALGEBRA: *)
+ type e
+ structure Base : FIELD where type t = e
+ structure ScalarMultiplication : SCALAR_MULTIPLY where type e = e and type t = t
+ val bilinear : unit
+
+ val fromList: e list -> t
+ val eval: t -> e -> e
+ val primitive: t -> bool
+ end
+
+functor PolynomialOverField(F : FIELD) : POLYNOMIAL =
+ struct
+ local
+ structure B = FieldPercent(F)
+ structure V = Vector
+ open B
+ open V
+ open VectorSlice
+
+ val zero = #%0
+ val one = #%1
+
+ fun constant c = full (tabulate (1, fn _ => c))
+ fun trim p n = subslice (p, 0, if n < length p then SOME n else NONE)
+ fun smul (x, y) = map (fn a => x *% a) y
+
+ (* reverse the polnomial's coeffs; p(1/x)*x^deg(p) *)
+ fun rev p =
+ let
+ val b = length p - 1
+ in
+ full (tabulate (length p, fn i => sub (p, b - i)))
+ end
+
+ fun normalize p =
+ let
+ val (v, s, l) = base p
+ fun lastzero i =
+ if i = s then 0 else
+ let
+ val i1 = i - 1
+ in
+ if V.sub (v, i1) =% zero then lastzero i1 else i-s
+ end
+ in
+ subslice (p, 0, SOME (lastzero (s+l)))
+ end
+
+ fun padd (x, y) =
+ let
+ val (v1, v2) =
+ if length x < length y then (x, y) else (y, x)
+ val len1 = length v1 (* shorter length *)
+ val len2 = length v2
+ fun add i =
+ if i < len1
+ then sub (v1, i) +% sub (v2, i)
+ else sub (v2, i)
+ in
+ full (tabulate (len2, add))
+ end
+
+ fun psub (x, y) =
+ let
+ val tab =
+ if length x < length y
+ then (length y,
+ let
+ val len = length x
+ in
+ fn i =>
+ if i < len
+ then sub (x, i) -% sub (y, i)
+ else ~% (sub (y, i))
+ end)
+ else (length x,
+ let
+ val len = length y
+ in
+ fn i =>
+ if i < len
+ then sub (x, i) -% sub (y, i)
+ else sub (x, i)
+ end)
+ in
+ full (tabulate tab)
+ end
+
+ (* long multiplication *)
+ fun long (x, y) = (* length x < length y *)
+ if length x = 1 then smul (sub (x, 0), y) else
+ let
+ val l = smul (sub (x, 0), y)
+ val hx = subslice (x, 1, NONE)
+ val h = long (hx, y)
+ fun join i =
+ if i = 0 then V.sub (l, 0) else
+ if i >= V.length l then V.sub (h, i-1) else
+ V.sub (l, i) +% V.sub (h, i-1)
+ in
+ tabulate (length x + length y - 1, join)
+ end
+
+ (* karatsuba's trick *)
+ fun kara (x, y) = (* length x < length y *)
+ let val (lx, ly) = (length x, length y) (* example: 5, 8 *) in
+ if lx < 5 then long (x, y) else
+ (* (a*x+b)*(c*x+d) = a*c*x + (a*d+b*c)*x + b*d
+ * (a+b)*(c+d) = (a*c + (a*d + b*c) + b*d)
+ *)
+ let
+ (* lx2 <= rx <= ry *)
+ val lx2 = lx div 2 (* example: 2 *)
+ val (rx, ry) = (lx - lx2, ly - lx2) (* example: 3, 6 *)
+ val lxy = lx + ly - 1 (* example: 12 *)
+ val lxh = lx2 + lx2 (* example: 4 *)
+ val lxf = lxh - 1 (* example: 3 *)
+ val rxy = rx + ry - 1 (* example: 8 *)
+ val lx2r = lx2 + rxy (* example: 10 *)
+
+ val a = subslice (x, lx2, NONE) (* length = rx (3) *)
+ val b = subslice (x, 0, SOME lx2) (* length = lx2 (2) *)
+ val c = subslice (y, lx2, NONE) (* length = ry (6) *)
+ val d = subslice (y, 0, SOME lx2) (* length = lx2 (2) *)
+
+ val ab = padd (a, b) (* length = rx (3) *)
+ val cd = padd (c, d) (* length = ry (6) *)
+ val abcd = kara (ab, cd) (* length = rxy (8) *)
+ val ac = kara (a, c) (* length = rxy (8) *)
+ val bd = kara (b, d) (* length = lxf (3) *)
+
+ val (adbc, _, _) =
+ base (psub(
+ full abcd, padd(full ac, full bd)))
+
+ (* v-- lxh
+ * ---ac---
+ * ^--adbc--
+ * | ^bd-
+ * | |^ 0
+ * | |lx2
+ * lx2r lxf
+ *)
+ (*!!! all these tests are grossly inefficient *)
+ fun join i =
+ case Int.compare (i, lxf) of
+ LESS => if i < lx2
+ then V.sub (bd, i)
+ else V.sub (adbc, i-lx2) +% V.sub (bd, i)
+ | EQUAL => V.sub (adbc, i-lx2)
+ | GREATER =>
+ if i < lx2r
+ then V.sub (ac, i-lxh) +% V.sub (adbc, i-lx2)
+ else V.sub (ac, i-lxh)
+ in
+ V.tabulate (lxy, join)
+ end end
+
+ (* !!! No FFT yet *)
+ fun pmul (x, y) =
+ if length x < length y
+ then if length x = 0 then x else full (kara (x, y))
+ else if length y = 0 then y else full (kara (y, x))
+
+ (* !!! write an inner product method *)
+
+ (* !!! don't be stupid; actually use the shortcut *)
+ fun pmuls (x, y, l) =
+ let
+ val xs = trim x l
+ val ys = trim y l
+ val z = pmul (xs, ys)
+ in
+ trim z l
+ end
+
+ (* p(x)*q(x) = 1 + r(x)*x^2^n
+ * => 2-p(x)*q(x) = 1 - r(x)*x^2^n
+ * => (2-p(x)*q(x))*p(x)*q(x) = 1 - r(x)^2*x^2^(n+1)
+ * => q' = (2-pq)q
+ *)
+ fun invmodn p n =
+ let
+ val init = constant (!% (sub (p, 0)))
+ val two = constant (#%2)
+ fun grow q =
+ if length q >= n then q else
+ let
+ val nl = length q + length q
+ val ps = trim p nl
+ val pq = pmul (q, ps) (* !!! we don't need high terms... *)
+ val neg = psub (two, pq)
+ (* !!! use inner product; we know low terms and don't need high *)
+ val nq = pmul (neg, q)
+ in
+ if length p = 1 then init else grow (trim nq nl)
+ end
+ in
+ trim (grow init) n
+ end
+ in
+ type e = F.t
+ type t = e slice
+ val characteristic = F.characteristic
+
+ fun eval p x = foldr (fn (a, out) => out *% x +% a) zero p
+ val fromList = normalize o full o fromList
+ fun primitive p = sub (p, length p - 1) =% one
+
+ structure Base = F
+
+ structure Addition =
+ struct
+ type t = t
+ val order = case Base.Addition.order of
+ UNCOUNTABLE => UNCOUNTABLE
+ | _ => COUNTABLE
+
+ val associative = ()
+ val commutative = ()
+ val one = full (tabulate (0, fn _ => #%0))
+
+ val EQ = fn (x, y) =>
+ let
+ val (v1, s1, l1) = base x
+ val (v2, s2, l2) = base y
+ val e1 = l1 + s1
+ fun compare (i1, i2) =
+ if i1 = e1 then true else
+ if V.sub (v1, i1) =% V.sub (v2, i2)
+ then compare (i1 + 1, i2 + 1)
+ else false
+ in
+ if l1 <> l2 then false else
+ compare (s1, s2)
+ end
+
+ val MUL = normalize o padd
+ val DIV = normalize o psub
+ val INV = fn x => full (map (fn a => ~%a) x)
+ end
+
+ structure Multiplication =
+ struct
+ type t = t
+ val order = Addition.order
+
+ val associative = ()
+ val commutative = ()
+ val one = full (tabulate (1, fn _ => #%1))
+
+ val EQ = Addition.EQ
+ val MUL = normalize o pmul
+ end
+
+ structure ScalarMultiplication =
+ struct
+ type e = e
+ type t = t
+
+ val associative = ()
+ val distributive = ()
+
+ val MUL = normalize o full o smul
+ end
+
+ val QUO = fn (z, y) =>
+ let
+ (* q*y + r = z, where deg(r) < deg(y) and deg(q)+deg(y)=deg(z)
+ * => rev(q*y + r) = rev(z)
+ * => rev(q)rev(y) + rev(r)*x^(deg(z)-deg(r)) = rev(z)
+ * => rev(q)rev(y) = rev(z) mod x^(1+deg(z)-deg(y))
+ * => rev(q) = rev(z)/rev(y) mod x^(1+deg(q))
+ *)
+ val degz = length z - 1
+ val degy = length y - 1
+ in
+ if degz < degy then (Addition.one, z) else
+ let
+ val degq = degz - degy
+ val lenq = degq + 1
+ val ry = rev y
+ val iry = invmodn ry lenq
+ val rz = rev z
+ val rq = pmuls (rz, iry, lenq)
+ val q = rev rq
+ val qy = pmuls (q, y, degy)
+ val r = psub (trim z degy, qy)
+ in
+ (q, normalize r)
+ end
+ end
+
+ fun LT (x, y) = length x < length y
+
+ val distributive = ()
+ val no_zero_divisors = ()
+ val bilinear = ()
+ end
+ end
+
+functor Polynomial(P : POLYNOMIAL) =
+ struct
+ local structure S = FieldPercent(P.Base) in open S end
+ local structure S = EuclideanDomainDollar(P) in open S end
+ local structure S = ScalarMultiply(P.ScalarMultiplication) in open S end
+ end
Added: mltonlib/trunk/ca/terpstra/math/q.fun
===================================================================
--- mltonlib/trunk/ca/terpstra/math/q.fun 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/q.fun 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,97 @@
+functor QuotientOfIntegralDomain(I : INTEGRAL_DOMAIN) : FIELD =
+ struct
+ local
+ structure B = IntegralDomainPercent(I)
+ structure O = EuclideanDomainDollar(Order)
+ open B
+ open O
+ in
+ type t = I.t * I.t
+ val characteristic = I.characteristic
+
+ structure Addition =
+ struct
+ type t = I.t * I.t
+ val order = COUNTABLE (* !!! no idea how to compute ... *)
+
+ val associative = ()
+ val commutative = ()
+ val one = (#%0, #%1)
+
+ val MUL = fn ((an, ad), (bn, bd)) =>
+ (an *% bd +% bn *% ad, ad *% bd)
+ val DIV = fn ((an, ad), (bn, bd)) =>
+ (an *% bd -% bn *% ad, ad *% bd)
+ val INV = fn (an, ad) => (~%an, ad)
+
+ val EQ = fn ((an, ad), (bn, bd)) =>
+ an *% bd =% bn *% ad
+ end
+
+ structure Multiplication =
+ struct
+ type t = I.t * I.t
+ val order = COUNTABLE (* !!! no idea *)
+
+ val associative = ()
+ val commutative = ()
+ val one = (#%1, #%1)
+
+ val EQ = Addition.EQ
+ val MUL = fn ((an, ad), (bn, bd)) =>
+ (an *% bn, ad *% bd)
+ val INV = fn (an, ad) => (ad, an)
+ val DIV = fn (a, b) => MUL (a, INV b)
+ end
+
+ val distributive = ()
+ val no_zero_divisors = ()
+ end
+ end
+
+functor QuotientOfEuclideanDomain(E : EUCLIDEAN_DOMAIN) : FIELD =
+ struct
+ local
+ structure B = EuclideanDomainPercent(E)
+ structure G = GCD(E)
+ structure Q = QuotientOfIntegralDomain(E)
+ open B
+ open G
+ in
+ type t = E.t * E.t
+ val characteristic = E.characteristic
+
+ fun simplify (n, d) =
+ if d =% #%0 then Q.Addition.one else
+ let
+ val g = gcd (n, d);
+ in
+ (n /% g, d /% g)
+ end
+
+ structure Addition =
+ struct
+ open Q.Addition
+
+ (* easier to compute *)
+ val EQ = fn ((an, ad), (bn, bd)) =>
+ ((an =% bn) andalso (ad =% bd)) orelse
+ ((an =% ~%bn) andalso (ad =% ~%bd))
+
+ val MUL = simplify o MUL
+ val DIV = simplify o DIV
+ end
+
+ structure Multiplication =
+ struct
+ open Q.Multiplication
+
+ val EQ = Addition.EQ
+ val MUL = simplify o MUL
+ val DIV = simplify o DIV
+ end
+
+ val distributive = ()
+ val no_zero_divisors = ()
+ end
+ end
Added: mltonlib/trunk/ca/terpstra/math/q.sml
===================================================================
--- mltonlib/trunk/ca/terpstra/math/q.sml 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/q.sml 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,7 @@
+structure Q = QuotientOfEuclideanDomain(Z)
+
+(* these are vulnerable to overflow *)
+structure Q64 = QuotientOfEuclideanDomain(Z64)
+structure Q32 = QuotientOfEuclideanDomain(Z32)
+structure Q16 = QuotientOfEuclideanDomain(Z16)
+structure Q8 = QuotientOfEuclideanDomain(Z8)
Added: mltonlib/trunk/ca/terpstra/math/r.fun
===================================================================
--- mltonlib/trunk/ca/terpstra/math/r.fun 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/r.fun 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,38 @@
+functor FieldOfReal(Real : REAL) : FIELD =
+ struct
+ type t = Real.real
+ val characteristic = LargeInt.fromInt 0
+
+ structure Addition =
+ struct
+ type t = Real.real
+ val order = UNCOUNTABLE
+
+ val associative = ()
+ val commutative = ()
+ val one = Real.fromInt 0
+
+ val EQ = Real.==
+ val MUL = Real.+
+ val DIV = Real.-
+ val INV = Real.~
+ end
+
+ structure Multiplication =
+ struct
+ type t = Real.real
+ val order = UNCOUNTABLE
+
+ val associative = ()
+ val commutative = ()
+ val one = Real.fromInt 1
+
+ val EQ = Real.==
+ val MUL = Real.*
+ val DIV = Real./
+ val INV = fn x => DIV (one, x)
+ end
+
+ val distributive = ()
+ val no_zero_divisors = ()
+ end
Added: mltonlib/trunk/ca/terpstra/math/r.sml
===================================================================
--- mltonlib/trunk/ca/terpstra/math/r.sml 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/r.sml 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,4 @@
+structure R = FieldOfReal(LargeReal)
+
+structure R32 = FieldOfReal(Real32)
+structure R64 = FieldOfReal(Real64)
Added: mltonlib/trunk/ca/terpstra/math/rings.fun
===================================================================
--- mltonlib/trunk/ca/terpstra/math/rings.fun 2006-12-19 02:58:00 UTC (rev 4985)
+++ mltonlib/trunk/ca/terpstra/math/rings.fun 2006-12-19 02:59:11 UTC (rev 4986)
@@ -0,0 +1,86 @@
+(****************************************************************** Ring ops *)
+
+functor NonAssociativeRingPercent(R : NON_ASSOCIATIVE_RING) =
+ struct
+ local structure S = BinaryOperationMulPercent(R.Multiplication) in open S end
+ local structure S = AbelianGroupAddPercent(R.Addition) in open S end
+ end
+functor NonAssociativeRingDollar(R : NON_ASSOCIATIVE_RING) =
+ struct
+ local structure S = BinaryOperationMulDollar(R.Multiplication) in open S end
+ local structure S = AbelianGroupAddDollar(R.Addition) in open S end
+ end
+
+functor RingPercent(R : RING) =
+ struct
+ local structure S = SemiGroupMulPercent(R.Multiplication) in open S end
+ local structure S = AbelianGroupAddPercent(R.Addition) in open S end
+ end
+functor RingDollar(R : RING) =
+ struct
+ local structure S = SemiGroupMulDollar(R.Multiplication) in open S end
+ local structure S = AbelianGroupAddDollar(R.Addition) in open S end
+ end
+
+functor UnitaryRingPercent(U : UNITARY_RING) =
+ struct
+ local structure S = MonoidMulPercent(U.Multiplication) in open S end
+ local structure S = AbelianGroupAddPercent(U.Addition) in open S end
+ val #% = fn e => U.Multiplication.one ++% e
+ end
+functor UnitaryRingDollar(U : UNITARY_RING) =
+ struct
+ local structure S = MonoidMulDollar(U.Multiplication) in open S end
+ local structure S = AbelianGroupAddDollar(U.Addition) in open S end
+ val #$ = fn e => U.Multiplication.one ++$ e
+ end
+
+functor CommutativeRingPercent(C : COMMUTATIVE_RING) =
+ struct
+ local structure S = AbelianMonoidMulPercent(C.Multiplication) in open S end
+ local structure S = AbelianGroupAddPercent(C.Addition) in open S end
+ val #% = fn e => C.Multiplication.one ++% e
+ end
+functor CommutativeRingDollar(C : COMMUTATIVE_RING) =
+ struct
+ local structure S = AbelianMonoidMulDollar(C.Multiplication) in open S end
+ local structure S = AbelianGroupAddDollar(C.Addition) in open S end
+ val #$ = fn e => C.Multiplication.one ++$ e
+ end
+
+functor IntegralDomainPercent(I : INTEGRAL_DOMAIN) = CommutativeRingPercent(I)
+functor IntegralDomainDollar (I : INTEGRAL_DOMAIN) = CommutativeRingDollar (I)
+
+functor EuclideanDomainPercent(E : EUCLIDEAN_DOMAIN) =
+ struct
+ local structure S = AbelianMonoidMulPercent(E.Multiplication) in open S end
+ local structure S = AbelianGroupAddPercent(E.Addition) in open S end
+ val #% = fn e => E.Multiplication.one ++% e
+ val (op /%) = #1 o E.QUO
+ val (op %%) = #2 o E.QUO
+ val ...
[truncated message content] |
