MySQL实现算法:矩阵快速幂

前几日在地铁上和好友讨论起MySQL数据库，由于MySQL数据库不支持复合类型变量，就想到能否用MySQL去实现一些算法呢。最近在学习MySql存储程序，昨天闲来无事，就动手写写吧，借此巩固一下存储过程的用法。

这次要用MySQL实现的算法是矩阵快速幂，并借此解决斐波那契数列问题，首先我们要简单普及一下基本知识。一，菲波那切数列。即我们熟知的兔子繁殖问题，菲波那切数列的数学表达式:F(n) = F(n-1)+F(n-2),(n>2,F(1)=F(2)=1);第n项为前两项之和。二，快速幂。顾名思义幂运算的快速实现，对于普通的a的b次幂运算，需要进行b次乘法运算，而快速幂只用log2(b)级别的次数就可以做到。至于算法原理这里就不再赘述。三，矩阵乘法运算。线性代数里面的知识。这里只使用公式表达一下。A(x,y) = A(x,y)+B(x,z)*C(z,y);四，矩阵快速幂。其实就是对于矩阵A的b次幂，采用快速幂的原理实现而已。

好了，下面我们就开始用MySQL实现矩阵快速幂。但是迎面而来的问题是，MySQL不支持复合类型，而我们需要用二位数组来实现矩阵的表示，并进行乘法运算。想到这里，你恐怕感觉无法实现了，但是我们有一个变通的方案可用，我们可以用数据表去充当一维的数组，使用自增的id用来标识第几个元素。既然有了一维数组，那我们完全可以用一维数组去充当二维数组啊，只需对下标进行适当的转换就行了。下标id是用1开始递增的，所以转换规则可以这样表示 A(x,y)=A((x-1)*len+y),(len为数组维度的长度),将二维下标转换成一维下标。数据类型的问题已经解决，那我们就首先实现矩阵相乘。首先我们要做的就是创建三张表(matrix,matrix_a,matrix_b)用来存储数据，数据表结构如下:

DROP TABLE IF EXISTS matrix;
CREATE TABLE matrix(
  id INT(11) NOT NULL AUTO_INCREMENT,
  val INT(11) DEFAULT NULL,
  PRIMARY KEY (id)
);

由于MySQL的函数不支持数据表的操作，我们只能借用存储过程来实现。把两个表名以字符串的形式传递进去，但是在存储过程中变量值又不能作为查询语句中的表名或者字段名，这是一个问题。还好MySQL提供的有一个替代方案，那就是动态查询，只需要编写好查询字符串，然后动态编译查询就行了。此时的存储过程代码如下：

DROP PROCEDURE IF EXISTS p_multiply;
DELIMITER $$;
CREATE PROCEDURE p_multiply(IN matrix_a VARCHAR(16),IN matrix_b VARCHAR(16),IN lengths INT)
BEGIN
	SET @sql_a = CONCAT('SELECT val INTO @temp_aa FROM ',matrix_a,' WHERE id=(?-1)*?+?');
	SET @sql_b = CONCAT('SELECT val INTO @temp_bb FROM ',matrix_b,' WHERE id=(?-1)*?+?');
	SET @len = lengths;
END;

貌似MySQL动态查询功能不支持DECLARE声明的变量，只能使用用户自定义变量。接下来就是一个三层的循序进行矩阵乘法运算。在运算中，我们只需使用以上说过的转换规则，从matrix_a和matrix_b表示的表中去除所要的数据，然后相乘之后添加的matrix表相对性的位置上。以下就是完整的代码:

DROP PROCEDURE IF EXISTS p_multiply;
DELIMITER $$;
CREATE PROCEDURE p_multiply(IN matrix_a VARCHAR(16),IN matrix_b VARCHAR(16),IN lengths INT)
BEGIN
	-- MySQL动态查询貌似不支持DECLARE声明的变量，只能使用用户自定义变量，即: SET @var = value;
	-- DECLARE sql_str_a,sql_str_b VARCHAR(128);
	SET @temp_aa = 0,@temp_bb = 0;
	SET @xi=1,@yi=1,@zi=1;
	
	SET @sql_a = CONCAT('SELECT val INTO @temp_aa FROM ',matrix_a,' WHERE id=(?-1)*?+?');
	SET @sql_b = CONCAT('SELECT val INTO @temp_bb FROM ',matrix_b,' WHERE id=(?-1)*?+?');
	SET @len = lengths;
	
	-- 初始化一个新矩阵，用来存储运算结果
	TRUNCATE TABLE matrix;
	WHILE @xi<=lengths*lengths DO
		BEGIN
			INSERT INTO matrix(val) VALUES (0);
			SET @xi = @xi+1;
		END;
	END WHILE;
	
	-- 矩阵相乘运算，A(xy) = A(xy)+B(xz)*C(zy);
	SET @xi=1;
	WHILE @xi<=lengths DO
		SET @yi=1;
		WHILE @yi<=lengths DO
			SET @zi=1;
			WHILE @zi<=lengths DO
				BEGIN
					PREPARE val_a FROM @sql_a;
					PREPARE val_b FROM @sql_b;
					EXECUTE val_a USING @xi,@len,@zi; DEALLOCATE PREPARE val_a;
					EXECUTE val_b USING @zi,@len,@yi; DEALLOCATE PREPARE val_b;
					UPDATE matrix SET val=val+@temp_aa*@temp_bb WHERE id=(@xi-1)*@len+@yi;
					SET @zi=@zi+1;
				END;
			END WHILE;
			SET @yi=@yi+1;
		END WHILE;
		SET @xi=@xi+1;
	END WHILE;
END;

紧接着就到到快速幂的核心部分，快速幂有好几种写法，这里使用我最喜欢的一种，无需递归无需数组。由于是矩阵快速幂，因此里面的运算需要改为基于矩阵的运算。需要一个单位矩阵开始，由于没有编程语言里面的数据类型的灵活性。一个简单的表达式matrix_b = matrix_b*matrix_a;我们就需要分成好几部实现，这里 matrix_b*matrix_a的运算已经有上面得到存储过程给出。但是赋值语句是必须的，每次矩阵相乘都会重置matrix表，必须将表中存储的运算结果复制出来，防止丢失。使用者两条语句来进行进行数据复制，同时保证id的正确性。TRUNCATE TABLE matrix_b;INSERT INTO matrix_b(val) SELECT val FROM matrix ORDER BY id;

DROP PROCEDURE IF EXISTS p_matrix_pow;
DELIMITER $$;
CREATE PROCEDURE p_matrix_pow(IN num INT,IN lengths INT)
BEGIN
	-- 初始化matrix_b为单位矩阵
	DECLARE i INT DEFAULT 1;
	TRUNCATE TABLE matrix_b;
	WHILE i<=lengths*lengths DO
		BEGIN
			IF (i-1)%lengths = (i-1) DIV lengths THEN
				INSERT INTO matrix_b(val) VALUES (1);
			ELSE
				INSERT INTO matrix_b(val) VALUES (0);
			END IF;
			SET i = i+1;
		END;
	END WHILE;
	
	-- 快速幂运算
	WHILE num>0 DO
		BEGIN
			IF num%2!=0 THEN
				BEGIN
					-- matrix_b = matrix_b*matrix_a;
					CALL p_multiply('matrix_b','matrix_a',lengths);
					TRUNCATE TABLE matrix_b;
					INSERT INTO matrix_b(val) SELECT val FROM matrix ORDER BY id;
				END;
			END IF;
			SET num = num DIV 2;
			-- matrix_a = matrix_a*matrix_a;
			CALL p_multiply('matrix_a','matrix_a',lengths);
			TRUNCATE TABLE matrix_a;
			INSERT INTO matrix_a(val) SELECT val FROM matrix ORDER BY id;
		END;
	END WHILE;
	-- return matrix_b;
END;

当快速幂进行结束，最终结果就在matrix_b表中。

说到这里，细心的可能会有疑问，矩阵快速幂和斐波那契数列有什么关系啊。斐波那契数列的数学表达式又不是幂运算，连乘法运算都没有。很明显，直接参照那个递推公式是不能使用矩阵快速幂的。这里需要用到一个不可思议的转换，此转换能适用很多递推表达式。我们再来看看这个表达式F(n) = F(n-1)+F(n-2)。此时将等号量表的位置替换一下，F(n-2)+F(n-1) = F(n);如果F(n)参与下次递推的话，就是这样的F(n-1)+F(n) = F(n+1);转换成矩阵乘法就如下所示:

  | 1  1 |   |F(n)  |   |F(n+1)|
  |      | * |      | = |      |
  | 1  0 |   |F(n-1)|   |F(n)  |

F(n+2)就可以这样表示:

 | 1  1 |   | 1  1 |   |F(n)  |   | 1  1 |   |F(n+1)|   |F(n+2)|
 |      | * |      | * |      | = |      | * |      | = |      |
 | 1  0 |   | 1  0 |   |F(n-1)|   | 1  0 |   |F(n)  |   |F(n+1)|

这样我们就可以依照矩阵乘法运算的结合律，先使用快速幂求出前面的常数矩阵的连乘部分，然后再和最后的一个2*1的矩阵相乘，求出最终的结果。

下面是执行求解的MySQL命令:

TRUNCATE TABLE matrix_a;
INSERT matrix_a(val) VALUES (1),(1),(1),(0);
CALL p_matrix_pow(42,2);  -- F(n) = F(n-1)+F(n-2),(n>2,F(1)=F(2)=1)
SELECT SUM(val) FROM matrix_b WHERE id=1 OR id=2;

快速幂的运算结果最终保存在matrix_b表中，同时此处省去了最后一步 矩阵相乘的运算，直接相加求的结果。 p_matrix_pow(42,2);实际求解的是F(44);F(1)和F(2)并适用这里的快速幂算法。

以上就是用MySQL实现矩阵快速幂的全过程。再次强调一下，MySQL并没有C/C++,java等编程语言在算法处理上的的灵活性，这里只是尝试实现一下，纯属娱乐。