I’ve been been on a kick learning about some cool math topics. In particular, I busted out my copy of Concrete Mathematics (awesome book) and was reading up on the number theory section. Bezout’s identity is that if you have ma+nb=d for some m,n and d divides a and b, then d is the greatest divisor of a and b. Bezout’s identity is a certificate theorem like the dual solution to a linear program. It doesn’t matter how you found m and n or d, Euclid’s algorithm or brute search, but once you’ve found that certificate, you know you have the biggest divisor. Pretty cool. Suppose you have some other common divisor d’. That means xd’=a and yd’=b. Subsitute this into the certificate. m(xd’)+n(yd’)=(mx+ny)d’ =d. This means d’ is a divisor of d. But the divisor relation implies the inequality relation (a divisor is always less than or equal to the thing it divides). Therefore d'<=d.

But another thing I’ve been checking out is algebraic geometry and Groebner bases. Groebner bases are a method for solving multivariate polynomial equations. The method is a relative of euclidean division and also in a sense Gaussian elimination. Groebner basis are a special recombination of a set of polynomial equations such that polynomial division works uniquely on them (in many variables polynomial division doesn’t work like you’d expect from the one variable case anymore). Somewhat surprisingly other good properties come along for the ride. More on this in posts to come

Some interesting stuff in the Haskell arena:

Gröbner bases in Haskell: Part I

https://konn.github.io/computational-algebra/

I love Doug McIlroy’s powser. It is my favorite functional pearl. https://www.cs.dartmouth.edu/~doug/powser.html

One interesting aspect not much explored is that you can compose multiple layers of [] like [[[Double]]] to get multivariate power series. You can then also tuple up m of these series to make power series from R^n -> R^m. An interesting example of a category and a nonlinear relative of the matrix category going from V^n -> V^m. I’m having a hell of a time getting composition to work automatically though. I’m a little stumped.

I made a version of division that uses the Integral typeclass rather than the fractional typeclass with an eye towards these applications. It was kind of annoying but I think it is right now.

I thought that in order to make things more Groebner basis like, I should make polynomials that start with the highest power first. It also makes sense for a division algorithm, but now I’m not so sure. I should also have probably used Vector and not List.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 |
bezout x 0 = (1, 0, x) bezout x y = (n, m - n*z , d) where r = x `mod` y z = x `div` y (m, n, d) = bezout y r -- my + nr = d -- my + n(x - zy) = d -- (m-nz)y + nx = d -- newtype Poly a = P [a] instance Num a => Num [a] where (+) xs [] = xs (+) [] ys = ys (+) xss@(x : xs) yss@(y : ys) | nx > ny = x : (xs + yss) | nx < ny = y : (xss + ys) | otherwise = x + y : (xs + ys) where nx = length xs ny = length ys (*) _ [] = [] (*) [] _ = [] (*) (x:xs) ys = ((* x) <$> ys ++ (replicate nx 0)) + (xs * ys) where nx = length xs -- (pure x * ys) * xpow nx + (xs * ys) where nx = length xs -- negate = fmap negate abs = error "no don't do this" signum = error "or this" fromInteger 0 = [] fromInteger x = (pure . fromInteger) x xpow n = 1 : replicate n 0 {- class DivMod a where divmod :: a -> a -> (a, a) instance Integral a => DivMod a where divmod = quotRem -} -- instance Integral a => Integral [a] where quotRem [] _ = ([], []) quotRem xss@(x : xs) yss@(y : ys) | nx < ny = ([], xss) | otherwise = (d' + d'', r' + r'') where -- | nx > ny = (d' + d'',r' + r'') -- | otherwise = ( pure d, r : xs - (pure d) * ys) where -- (pure d , xss - (pure d) * yss) -- ( pure d, r : xs - (pure d) * ys) -- | x == 0 = quotRem xs yss nx = length xs ny = length ys (d , r) = x `quotRem` y d' = d : (replicate (nx - ny) 0) r' = r : (replicate (nx - ny) 0) -- (d'' ,r'') = quotRem (xss - d' * yss) yss (d'' , r'') = quotRem (xs - (pure d) * ys) yss toInteger = error "don't do this" instance (Num a, Ord a) => Real [a] where toRational = error "screw you haskell" instance (Enum a) => Enum [a] where toEnum = error " crew you" -- pure . toEnum fromEnum = error "screw yo " -- fromEnum . head {- -- actually this is a different Ord instance from the main library instance Ord a => Ord [a] where compare xs yy = case compare nx ny of EQ -> compare (head xs) (head ys) z -> z where nx = length xs ny = length ys -} -- if we use powser with Vector, -- we can just grab the tail for the remainer. --it makes sense to grab the end -- we can just reuse the powser infracsture -- with a tail. -- tail view in agda for that purpose also. example :: [[Double]] example = [1, 1] var :: Num a => [a] var = [1,0] x :: Num a => [a] x = var y :: Num a => [[a]] y = pure x z :: Num a => [[[a]]] z = pure y x' (x, _, _) = x y' (_, y, _) = y f :: [[Double]] f = 1 + 2 * x + 3 * x * y + 3 * y --(a , (a ,( a ,()))) -- '[a , b, c, d] f1 :: (Num a) => (a, a, a) -> (a, a) f1 (x, y, z) = (x, y + z) f2 :: (Num a) => (a, a) -> a f2 (x, y) = x*x g = f2 . f1 example2 = g (x, y, z) -- only a single variable -- we need something else use :: Num a => [a] -> a -> a use [] z = 0 use (x : xs) z = x * (z ^ (length xs)) + (use xs z) |