前4章以及其他补充已经整理如下：

本章我们将学习矩阵和通用函数（universal functions，即ufuncs）的相关内容。矩阵作为一种重要的数学概念，在NumPy中也有专门的表示方法。通用函数可以逐个处理数组中的元素，也可以直接处理标量。通用函数的输入是一组标量，输出也是一组标量，它们通常可以对应于基本数学运算，如加、减、乘、除等。我们还将介绍三角函数、位运算函数和比较函数。

本章涵盖以下内容：

矩阵创建；
矩阵运算；
基本通用函数；
三角函数；
位运算函数；
比较函数。

这一系列，我都是用Jupyter notebook进行学习。环境启动如下：
1
2
workon env1
jupyter notebook

5.1 矩阵

在NumPy中，矩阵是ndarray的子类，可以由专用的字符串格式来创建。与数学概念中的矩阵一样，NumPy中的矩阵也是二维的。如你所料，矩阵的乘法运算和NumPy中的普通乘法运算不同。幂运算当然也不一样。我们可以使用mat、matrix以及bmat函数来创建矩阵。

5.2 动手实践：创建矩阵

mat函数创建矩阵时，若输入已为matrix或ndarray对象，则不会为它们创建副本。因此，调用mat函数和调用matrix(data, copy=False)等价。我们还将展示矩阵转置和矩阵求逆的方法。

补充：在 NumPy 中，mat、matrix 和 bmat 是用于创建矩阵对象的函数。它们都属于 numpy.matlib 模块（或者是 numpy 本身的一部分，依赖于 NumPy 版本），但有些细节和功能差异。以下是它们的详细说明，包括函数签名、参数、返回值以及使用示例。

mat 函数：

mat 函数用于创建矩阵对象。矩阵是二维数组（ndarray）的一个特定类型，提供了一些额外的矩阵操作。

函数签名：

numpy.mat(data, dtype=None)

参数：

data: 用于创建矩阵的输入数据。可以是嵌套列表（列表的列表），二维数组，或其他类似对象。数据会被转化为一个矩阵。
dtype: 可选，指定输出矩阵的数据类型。如果没有指定，dtype 将从输入数据中推断。

输出：

返回一个矩阵对象（numpy.matrix）。

例如：

# mat创建矩阵
m1 = np.mat([[1, 2], [3, 4]])
m1

matrix 函数

matrix 函数与 mat 函数非常相似，实际上它们是等效的，都是创建矩阵对象的方式。在较早的版本中，matrix 是用于创建矩阵的函数，但随着 NumPy 的更新，mat 和 matrix 变得几乎没有差别。大多数情况下建议使用 np.array() 而不是 matrix，因为后者对矩阵进行的一些特定操作已不再推荐。

函数签名：

numpy.matrix(data, dtype=None)

参数：

data: 输入数据，可以是嵌套列表、二维数组等。
dtype: 可选，指定输出矩阵的数据类型。

输出：

返回一个矩阵对象（numpy.matrix）。

例如：

m2 = np.matrix([[1, 2], [3, 4]])
m2

进行矩阵运算：

print("m1*m1", m1*m1)
print("m2*m2", m2*m2)

bmat 函数

bmat 函数用于创建一个由多个矩阵块（block）组成的大矩阵。它是用来将多个矩阵按指定的位置（行和列）拼接成一个大的矩阵。

函数签名：

numpy.bmat(b)

参数：

b: 输入可以是一个列表、元组或字符串，表示多个矩阵块的排列。每个矩阵块也可以是矩阵、数组或其他支持矩阵操作的对象。矩阵将按照给定的顺序拼接在一起。
- 如果是字符串，它表示一个矩阵块的简单描述（如“AB; CD”表示两个矩阵块的组合）。
- 如果是列表或元组，它应包含可以拼接在一起的矩阵。

输出：返回一个拼接好的矩阵对象。

例如：

# 创建单独的矩阵块
A = np.matrix([[1, 2], [3, 4]])
B = np.matrix([[5, 6], [7, 8]])
C = np.matrix([[9, 10], [11, 12]])

# 使用 bmat 拼接矩阵
block_matrix = np.bmat([[A, B], [C, A]])
print("Block Matrix:")
print(block_matrix)

总结对比

函数	功能	参数	返回值	示例
`mat`	创建一个矩阵对象（与 `matrix` 类似）	`data`（数据） `dtype`（数据类型）	`numpy.matrix`	`np.mat([[1, 2], [3, 4]])`
`matrix`	创建一个矩阵对象（与 `mat` 类似）	`data`（数据） `dtype`（数据类型）	`numpy.matrix`	`np.matrix([[1, 2], [3, 4]])`
`bmat`	拼接多个矩阵块构成一个大矩阵	`b`（矩阵块的列表、元组或字符串）	拼接后的 `numpy.matrix`	`np.bmat([[A, B], [C, A]])`

mat 和 matrix：功能相同，用于创建矩阵对象，推荐使用 mat，因为 matrix 在新版的 NumPy 中已经不再推荐使用。

bmat：用于按指定的块结构拼接多个矩阵块，特别适合处理更复杂的矩阵拼接问题。

(1) 在创建矩阵的专用字符串中，矩阵的行与行之间用分号隔开，行内的元素之间用空格隔开。使用如下的字符串调用mat函数创建阵：

A = np.mat('1 2 3;4 5 6;7 8 9')
A

(2) 用T属性获取转置矩阵：

A.T

(3) 用I属性获取逆矩阵：

A.I

但是出现报错：

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
Cell In[11], line 1
----> 1 A.I

File /Library/Python/3.9/site-packages/numpy/matrixlib/defmatrix.py:836, in matrix.I(self)
    834 else:
    835     from numpy.linalg import pinv as func
--> 836 return asmatrix(func(self))

File /Library/Python/3.9/site-packages/numpy/linalg/linalg.py:561, in inv(a)
    559 signature = 'D->D' if isComplexType(t) else 'd->d'
    560 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 561 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    562 return wrap(ainv.astype(result_t, copy=False))

File /Library/Python/3.9/site-packages/numpy/linalg/linalg.py:112, in _raise_linalgerror_singular(err, flag)
    111 def _raise_linalgerror_singular(err, flag):
--> 112     raise LinAlgError("Singular matrix")

LinAlgError: Singular matrix

这个错误提示表明你尝试对一个“奇异矩阵”（singular matrix）进行求逆。奇异矩阵指的是行列式为零的矩阵，它没有逆矩阵，无法进行求逆运算。简单来说，当一个矩阵的行或列是线性相关的时，它是奇异的，无法进行求逆。

换一个非奇异矩阵求逆：

(4) 除了使用字符串创建矩阵以外，我们还可以使用NumPy数组进行创建：

m = np.mat(np.arange(9).reshape(3, 3))
m

刚才做了些什么 : 我们使用mat函数创建了矩阵，用T属性获取了转置矩阵，用I属性获取了逆矩阵。

5.3 从已有矩阵创建新矩阵

有些时候，我们希望利用一些已有的较小的矩阵来创建一个新的大矩阵。这可以用bmat函数来实现。这里的b表示“分块”，bmat即分块矩阵（block matrix）。

5.4 动手实践：从已有矩阵创建新矩阵

我们将利用两个较小的矩阵创建一个新的矩阵，步骤如下。

(1) 首先，创建一个2×2的单位矩阵：

n = np.eye(2)
n

创建另一个与A同型的矩阵，并乘以2：

n1 = 2*n
n1

(2) 使用字符串创建复合矩阵，该字符串的格式与mat函数中一致，只是在这里你可以用矩阵变量名代替数字：

np.bmat('n n1;n n1')

刚才做了些什么 : 我们使用bmat函数，从两个小矩阵创建了一个分块复合矩阵。我们用矩阵变量名替代了数字，并将字符串传给bmat函数。

突击测验：使用字符串定义矩阵

问题1 在使用mat和bmat函数创建矩阵时，需要输入字符串来定义矩阵。在字符串中，以下哪一个英文标点符号是矩阵的行分隔符？ (1) 分号 “;” (2) 冒号 “:” (3) 逗号 “,” (4) 空格 “ ”

答案：（1），解释：在使用 mat 和 bmat 函数创建矩阵时，如果使用字符串来定义矩阵，矩阵的行分隔符是 分号 (;)。

5.5 通用函数

通用函数的输入是一组标量，输出也是一组标量，它们通常可以对应于基本数学运算，如加、减、乘、除等。

5.6 动手实践：创建通用函数

我们可以使用NumPy中的frompyfunc函数，通过一个Python函数来创建通用函数，步骤如下。

(1) 定义一个回答宇宙、生命及万物的终极问题的Python函数（问题和答案来源于《银河系漫游指南》，如果你没看过，可以忽略）：

def ultimate_answer(a):

到这里为止还没有什么特别的。我们将这个函数命名为ultimate_answer，并为之定义了一个参数a。

(2) 使用zeros_like函数创建一个和a形状相同，并且元素全部为0的数组result：

results = np.zeros_like(a)

(3) 现在，我们将刚刚生成的数组中的所有元素设置为“终极答案”其值为42，并返回这个结果。完整的函数代码如下所示。flat属性为我们提供了一个扁平迭代器，可以逐个设置数组元素的值：

def ultimate_answer(a):
    results = np.zeros_like(a)
    results.flat = 42
    return results

(4) 使用frompyfunc创建通用函数。指定输入参数的个数为1，随后的1为输出参数的个数：

ufunc = np.frompyfunc(ultimate_answer, 1, 1)
print(ufunc(np.arange(4)))

我们可以对二维数组进行完全一样的操作，代码如下：

print("The answer", ufunc(np.arange(4).reshape(2, 2)))

刚才做了些什么 : 我们定义了一个Python函数。其中，我们使用zeros_like函数根据输入参数的形状初始化一个全为0的数组，然后利用ndarray对象的flat属性将所有的数组元素设置为“终极答案”其值为42。

完整代码：

def ultimate_answer(a):
    results = np.zeros_like(a)
    results.flat = 42
    return results
ufunc = np.frompyfunc(ultimate_answer, 1, 1)
# print("The answer", ufunc(np.arange(4)) )
print("The answer", ufunc(np.arange(4).reshape(2, 2)))

补充：frompyfunc函数的使用方法

frompyfunc 是 NumPy 中的一个函数，它允许你将普通的 Python 函数转换为可以在 NumPy 数组上逐元素操作的“ufunc”（即用户自定义函数）。ufunc 是一种高效的函数，支持在整个数组上进行元素级运算。

函数签名：

numpy.frompyfunc(func, nin, nout)

参数：

func：一个普通的 Python 函数，它将被转换为 ufunc。此函数可以接受多个输入，并返回一个输出或多个输出。
nin：输入参数的数量（即该函数需要接受的输入数量）。nin 是一个整数，表示函数的输入参数个数。
nout：输出参数的数量（即该函数返回的输出数量）。nout 也是一个整数，表示该函数应该返回的输出数量。

返回值：

返回一个 ufunc 对象，这个对象可以像 NumPy 中的其他 ufunc 一样，被用来逐元素地应用到数组。

例如：

def add(a, b):
    return a + b

func_add = np.frompyfunc(add, 2, 1)
n1 = np.arange(8).reshape(2, 4)
n2 = np.array([[1, 3, 5, 7], [2, 4, 6, 8]])

print(n1)
print(n2)
func_add(n1, n2)

5.7 通用函数的方法

函数竟然也可以拥有方法？如前所述，其实通用函数并非真正的函数，而是能够表示函数的对象。通用函数有四个方法，不过这些方法只对输入两个参数、输出一个参数的ufunc对象有效，例如add函数。其他不符合条件的ufunc对象调用这些方法时将抛出ValueError异常。因此只能在二元通用函数上调用这些方法。以下将逐一介绍这4个方法：

reduce
accumulate
reduceat
outer

5.8 动手实践：在 add 上调用通用函数的方法

我们将在add函数上分别调用4个方法。

(1) 沿着指定的轴，在连续的数组元素之间递归调用通用函数，即可得到输入数组的规约（reduce）计算结果。对于add函数，其对数组的reduce计算结果等价于对数组元素求和。调用reduce方法：

A = np.arange(9)
print("A", A)
print(np.add.reduce(A))

(2) accumulate方法同样可以递归作用于输入数组。但是与reduce方法不同的是，它将存储运算的中间结果并返回。因此在add函数上调用accumulate方法，等价于直接调用cumsum函数。在add函数上调用accumulate方法：

np.add.accumulate(A)

(3) reduceat方法解释起来有点复杂，我们先运行一次，再一步一步来看它的算法。reduceat方法需要输入一个数组以及一个索引值列表作为参数。

np.add.reduceat(A,[0, 5, 2, 7])

第一步用到索引值列表中的0和5，实际上就是对数组中索引值在0到5之间的元素进行reduce操作。

np.add.reduce(A[0:5])

第二步用到索引值5和2。由于2比5小，所以直接返回索引值为5的元素：

A[5]

第三步用到索引值2和7。这一步是对索引值在2到7之间的数组元素进行reduce操作：

np.add.reduce(A[2:7])

第四步用到索引值7。这一步是对索引值从7开始直到数组末端的元素进行reduce操作：

np.add.reduce(A[7:])

(4) outer方法返回一个数组，它的秩（rank）等于两个输入数组的秩的和。它会作用于两个输入数组之间存在的所有元素对。在add函数上调用outer方法：

np.add.outer(np.arange(3),A)

刚才做了些什么 : 我们在通用函数add上调用了四个方法：reduce、accumulate、reduceat以及outer。

总结补充：在 NumPy 中，np.add 是一个常用的 ufunc（通用函数），它用于执行逐元素的加法运算。除了基本的逐元素加法，np.add 还提供了四种额外的操作方法：reduce、accumulate、reduceat 和 outer，它们可以用来进行更复杂的操作。

reduce 方法将给定的函数（这里是加法）逐步应用到数组的元素上，最终返回一个单一的结果。它可以看作是 归约操作，即通过重复应用操作（如加法）来将一个多元素数组“压缩”成一个标量。

accumulate 方法对数组执行逐步的累加操作，返回一个新的数组，其中每个元素是从原始数组开始的累加结果。它可以看作是 前缀和（prefix sum），即每个位置的元素是从数组开始到该位置的所有元素的加和。

reduceat 方法允许你在指定的“切片”位置上应用 reduce 操作。它能在给定的索引位置对数组进行局部归约操作，返回的结果是一个数组，它包含了每个“切片”的归约结果。

outer 方法用于计算两个数组的外积。具体来说，它会对两个数组中的每对元素应用 np.add 函数，生成一个矩阵，其中第 i 行第 j 列的元素是 arr1[i] + arr2[j]。

5.9 算术运算

在NumPy中，基本算术运算符+、-和*隐式关联着通用函数add、subtract和multiply。

也就是说，当你对NumPy数组使用这些算术运算符时，对应的通用函数将自动被调用。除法包含的过程则较为复杂，在数组的除法运算中涉及三个通用函数divide、true_divide和floor_division，以及两个对应的运算符/和//。

补充：

np.divide 是执行标准除法的函数，它对两个数组的每一对元素进行除法运算，返回的结果通常是浮动类型的结果，即使是整数类型的输入数组也会得到浮动结果。

函数签名：

np.divide(x1, x2)

x1 和 x2：输入数组或标量，表示要进行除法的操作数。

返回值：返回一个与 x1 和 x2 形状相同的数组，其中每个元素是 x1[i] / x2[i] 的结果。

例子：

x1 = np.array([8, 6, 4])
x2 = np.array([2, 3, 2])
np.divide(x1, x2)

np.true_divide

np.true_divide 也是执行除法操作，但是它始终返回 浮动类型 结果。无论输入的数组元素类型是整数还是浮动类型，np.true_divide 都会强制转换输出为浮动类型。

函数签名：

np.true_divide(x1, x2)

返回值：返回的结果总是浮动类型，甚至当两个整数相除时。

x1 = np.array([8, 6, 4])
x2 = np.array([2, 3, 2])
np.true_divide(x1, x2)

np.floor_divide

np.floor_divide 执行的是 向下取整的除法。与普通除法不同，它会返回 整数类型的结果，并且结果是通过对标准除法的结果向下取整得到的。

函数签名：

np.floor_divide(x1, x2)

返回值：返回的结果是整数类型，且是向下取整的商。

举例：

arr1 = np.array([7, 10, 5])
arr2 = np.array([2, 5, 2])
np.floor_divide(arr1, arr2)

/ 是除法运算符，功能与 np.divide 函数相同。它执行标准的除法运算，返回浮动类型的结果。

// 是整数除法运算符，功能与 np.floor_divide 函数相同。它执行 向下取整的除法，结果是整数类型，且结果是通过将标准除法结果向下取整得到的。

5.10 动手实践：数组的除法运算

让我们在实践中了解数组的除法运算。

(1) divide函数在整数和浮点数除法中均只保留整数部分：

注：原书写错了，np.divide 是执行标准除法的函数，它对两个数组的每一对元素进行除法运算，返回的结果通常是浮动类型的结果，即使是整数类型的输入数组也会得到浮动结果。

arr1 = np.array([2, 6, 5])
arr2 = np.array([1, 2, 3])
np.divide(arr1, arr2)
np.divide(arr2, arr1)

(2) true_divide函数与数学中的除法定义更为接近，即返回除法的浮点数结果而不作截断：

arr1 = np.array([2, 6, 5])
arr2 = np.array([1, 2, 3])
print(np.true_divide(arr1, arr2))
print(np.true_divide(arr2, arr1))

(3) floor_divide函数总是返回整数结果，相当于先调用divide函数再调用floor函数。floor函数将对浮点数进行向下取整并返回整数：

arr1 = np.array([2, 6, 5])
arr2 = np.array([1, 2, 3])
print(np.floor_divide(arr1, arr2))
print(np.floor_divide(arr2, arr1))

(4) 默认情况下，使用/运算符相当于调用divide函数：

arr1 = np.array([2, 6, 5])
arr2 = np.array([1, 2, 3])
print(arr1/arr2)
print(arr2/arr1)

注，原书可能应该基于Python2，但是在Python3中，无论是否引用 __future__ 模块，除法行为都是一致的.

from __future__ import division
arr1 = np.array([2, 6, 5])
arr2 = np.array([1, 2, 3])
print(arr1/arr2)
print(arr2/arr1)

(5) 运算符//对应于floor_divide函数。例如下面的代码：

arr1 = np.array([2, 6, 5])
arr2 = np.array([1, 2, 3])
print(arr1//arr2)
print(arr2//arr1)

刚才做了些什么 : 我们学习了NumPy中三种不同的除法函数。其中，divide函数在整数和浮点数除法中均只保留整数部分(其实Python3不是这样的，是保留小数的)，true_divide函数不作截断返回浮点数结果，而floor_divide函数同样返回整数结果并等价于先调用divide函数再调用floor函数。

勇敢出发：尝试`future.division`

动手实验，验证引入__future__.division的效果。

上面已有说明。不再重复

5.11 模运算

计算模数或者余数，可以使用NumPy中的mod、remainder和fmod函数。当然，也可以使用%运算符。这些函数的主要差异在于处理负数的方式。fmod函数在这方面异于其他函数。

补充：

numpy.mod()

mod 函数用于计算元素除法的余数，结果与 Python 的 % 操作符类似。numpy.mod() 和 Python 的 % 操作符一样，遵循如下规则：

余数的符号与除数相同。

函数签名：

numpy.mod(x1, x2)

x1 和 x2 是要进行操作的数组或数值。
返回值是 x1 除以 x2 的余数，符号与除数 (x2) 相同。

例子：

# 基本使用
print(np.mod(10, 3))  # 输出 1
print(np.mod(-10, 3))  # 输出 2
print(np.mod(10, -3))  # 输出 -2
print(np.mod(-10, -3))  # 输出 -1

注：如何理解上面的结果：

参考公式：remainder(x1,x2)=x1−x2×floor(x1/x2)

例如：-10 / 3 = -3.333...，向下取整为 -4，然后 -10 - 3 * (-4) = -10 + 12 = 2。

numpy.remainder()

numpy的取余运算和取模运算是一样的

【题目】：下列代码中np.remainder(m,2)输出的结果是？

import numpy as np
m = np.array([4, 5, 6])
print("【显示】m =", m)
print("【显示】np.remainder(m,2)：", np.remainder(m,2))
print("【显示】np.mod(m,2)：", np.mod(m,2))
print("【显示】m%2：", m%2)
print("【显示】np.mod.__name__：", np.mod.__name__)

A选项：[0 1 0]
B选项：[1 1 0]
C选项：[0 1 1]
D选项：[0 0 0]

正确答案是：A

例子：

print(np.remainder(10, 3))  # 输出1
print(np.remainder(-10, 3)) # 输出2
print(np.remainder(10, -3)) # 输出-2 
print(np.remainder(-10, -3)) # 输出-1

numpy.fmod():和上面两个不一样。所得余数的正负由被除数决定，与除数的正负无关：

# 基本使用
print(np.fmod(10, 3))  # 输出 1.0
print(np.fmod(-10, 3))  # 输出 -1.0
print(np.fmod(10, -3))  # 输出 1.0
print(np.fmod(-10, -3))  # 输出 -1.0

Python 中的 % 运算符也是计算余数的运算符，它和 numpy.mod() 基本一致。两者在大多数情况下是可以互换的，区别在于 numpy.mod() 支持广播（broadcasting），而 % 只能用于两个标量值或者对应形状的数组。

arr1 = np.array([10, -10, 10, -10])
arr2 = np.array([3, 3, -3, -3])

print(arr1 % arr2)  # 使用 % 运算符
print(np.mod(arr1, arr2))  # 使用 np.mod()

5.12 动手实践：模运算

我们将逐一调用前面提到的函数。

(1) remainder函数逐个返回两个数组中元素相除后的余数。如果第二个数字为0，则直接返回0：

n = np.arange(-4, 4)
print(n)
print(np.remainder(n, 2))

(2) mod函数与remainder函数的功能完全一致：

n = np.arange(-4, 4)
print(n)
print(np.mod(n, 2))

(3) %操作符仅仅是remainder函数的简写：

n = np.arange(-4, 4)
print(n)
print(n % 2)

(4) fmod函数处理负数的方式与remainder、mod和%不同。所得余数的正负由被除数决定，与除数的正负无关：

n = np.arange(-4, 4)
print(n)
print(np.fmod(n, 2))

刚才做了些什么 : 我们学习了NumPy中的mod、remainder和fmod等模运算函数。

5.13 斐波那契数列

斐波那契（Fibonacci）数列是基于递推关系生成的。直接用NumPy代码来解释递推关系是比较麻烦的，不过我们可以用矩阵的形式或者黄金分割公式来解释它。因此，我们将介绍matrix 和rint函数。使用matrix函数创建矩阵，rint函数对浮点数取整，但结果仍为浮点数类型。

5.14 动手实践：计算斐波那契数列

斐波那契数列的递推关系可以用矩阵来表示。斐波那契数列的计算等价于矩阵的连乘。

(1) 创建斐波那契矩阵：

m = np.matrix([[1, 1], [1, 0]])
m

(2) 计算斐波那契数列中的第8个数，即矩阵的幂为8减去1。计算出的斐波那契数位于矩阵的对角线上：

m_8 = m**7
m_8
(m**7)[0, 0]

(3) 利用黄金分割公式或通常所说的比奈公式（Binet’ s Formula），加上取整函数，就可以直接计算斐波那契数。计算前8个斐波那契数：

n = np.arange(1, 9)
sqrt5 = np.sqrt(5)
phi = (1 + sqrt5)/2
fibonacci = np.rint((phi**n - (-1/phi)**n)/sqrt5)
fibonacci

刚才做了些什么 : 我们分别用两种方法计算了斐波那契数列。在这个过程中，我们学习使用matrix函数创建矩阵，以及使用rint函数对浮点数取整但不改变浮点数类型。

补充：rint函数的使用

numpy.rint() 是 NumPy 中用于对数组或数字中的每个元素进行四舍五入操作的函数。它的作用是将每个元素四舍五入到最接近的整数，返回的是一个与输入形状相同的数组。

函数签名：

numpy.rint(x)

x: 输入数组或数值。可以是一个标量值、列表、NumPy 数组等。

返回值：

返回一个与输入数组或标量形状相同的数组，其中每个元素被四舍五入到最近的整数。

arr = np.array([3.6, 2.4, 5.5, -4.5, -2.3])
print(np.rint(arr))  # 输出 [ 4.  2.  6. -4. -2.]

勇敢出发：分析计算耗时

你可能很想知道究竟哪种方法计算更快，那就分析一下它们的耗时吧。使用frompyfunc 创建一个计算斐波那契数列的通用函数，并进行计时。

回顾：使用 frompyfunc 创建通用函数： np.frompyfunc 是一个用于将普通 Python 函数转换为通用函数的方法，这样就能在数组的每个元素上应用该函数，并且支持 NumPy 数组广播。我们将用它来定义一个计算斐波那契数列的函数。

GPT写的答案：

import numpy as np
import time

# 方法1：矩阵法
def fib_matrix(n):
    m = np.matrix([[1, 1], [1, 0]])
    return (m**(n-1))[0, 0]

# 方法2：Binet公式
def fib_binet(n):
    sqrt5 = np.sqrt(5)
    phi = (1 + sqrt5) / 2
    return np.rint((phi**n - (-1/phi)**n) / sqrt5)

# 使用 np.frompyfunc 创建通用函数来计算斐波那契数列
def fib_generic(n):
    sqrt5 = np.sqrt(5)
    phi = (1 + sqrt5) / 2
    return np.rint((phi**n - (-1/phi)**n) / sqrt5)

# 将普通函数转换为通用函数
fib_frompyfunc = np.frompyfunc(fib_generic, 1, 1)

# 计时对比
n_vals = np.arange(1, 1000001)  # 测试的斐波那契数列项（1 到 10000）

# 测量方法1（矩阵法）耗时
start_time = time.time()
matrix_results = np.array([fib_matrix(n) for n in n_vals])
method1_time = time.time() - start_time
print(f"方法1（矩阵法）耗时: {method1_time:.6f}秒")

# 测量方法2（Binet公式法）耗时
start_time = time.time()
binet_results = np.array([fib_binet(n) for n in n_vals])
method2_time = time.time() - start_time
print(f"方法2（Binet公式法）耗时: {method2_time:.6f}秒")

# 测量方法3（通用函数法）耗时
start_time = time.time()
generic_results = fib_frompyfunc(n_vals)
method3_time = time.time() - start_time
print(f"方法3（通用函数法）耗时: {method3_time:.6f}秒")

5.15 利萨茹曲线

在NumPy中，所有的标准三角函数如sin、cos、tan等均有对应的通用函数。利萨茹曲线（Lissajous curve）是一种很有趣的使用三角函数的方式。我至今仍记得在物理实验室的示波器上显示出利萨茹曲线时的情景。利萨茹曲线由以下参数方程定义：

x = A sin(at + n/2)  
y = B sin(bt) 

5.16 动手实践：绘制利萨茹曲线

利萨茹曲线的参数包括A、B、a和b。为简单起见，我们令A和B均为1。

(1) 使用linspace函数初始化变量t，即从-pi到pi上均匀分布的201个点：

import sys

# a = float(sys.argv[1])
# b = float(sys.argv[2])
a = 1
b = 1
t = np.linspace(-np.pi, np.pi, 201)

(2) 使用sin函数和NumPy常量pi计算变量x：

x = np.sin(a * t + np.pi / 2)

(3) 使用sin函数计算变量y：

y = np.sin(b * t)

(4) 我们将在第9章中详细讲解Matplotlib的用法。绘制的曲线如下所示。

plot(x ,y)
show()

刚才做了些什么 : 我们根据参数方程的定义，以参数A=B=1、a=9和b=8绘制了利萨茹曲线。我们还使用了sin和linspace函数，以及NumPy常量pi。

完整代码如下：

import sys
from matplotlib.pyplot import plot
from matplotlib.pyplot import show

# a = float(sys.argv[1])
# b = float(sys.argv[2])
a = 1
b = 1
t = np.linspace(-np.pi, np.pi, 201)
x = np.sin(a * t + np.pi / 2)
y = np.sin(b * t)
plot(x ,y)
show()

5.17 方波

方波也是一种可以在示波器上显示的波形。方波可以近似表示为多个正弦波的叠加。事实上，任意一个方波信号都可以用无穷傅里叶级数来表示。

傅里叶级数（Fourier series）是以正弦函数和余弦函数为基函数的无穷级数，以著名的法国数学家Jean-Baptiste Fourier命名。

方波可以表示为如下的傅里叶级数。

\[\sum_{k = 1}^{\infty} \frac{4 sin(2k-1)t}{(2k-1) \pi}\]

5.18 动手实践：绘制方波

与前面的教程中一样，我们仍将以相同的方式初始化t和k。我们需要累加很多项级数，且级数越多结果越精确，这里取k=99以保证足够的精度。绘制方波的步骤如下。

(1) 我们从初始化t和k开始，并将函数值初始化为0：

import sys
t = np.linspace(-np.pi, np.pi, 201)
# k = np.arange(1, float(sys.argv[1]))
k = np.arange(1, 99)
k = 2 * k -1
f = np.zeros_like(t)

(2) 接下来，直接使用sin和sum函数进行计算：

for i in range(len(t)): 
    f[i] = np.sum(np.sin(k * t[i])/k)  
f = (4 / np.pi) * f 

(3) 绘制波形的代码和前面的教程中几乎一模一样：

plot(t, f)
show

刚才做了些什么 : 我们使用sin函数生成了一个方波，或者起码是非常接近于方波的波形。函数的输入值是用linspace产生的，而一组k值是用arange函数生成的。

完整代码如下：

import sys
from matplotlib.pyplot import plot
from matplotlib.pyplot import show

t = np.linspace(-np.pi, np.pi, 201)
# k = np.arange(1, float(sys.argv[1]))
k = np.arange(1, 99)
k = 2 * k -1
f = np.zeros_like(t)

for i in range(len(t)): 
    f[i] = np.sum(np.sin(k * t[i])/k)  
f = (4 / np.pi) * f 

plot(t, f)
show

勇敢出发：摆脱循环语句

你可能已经注意到，在代码中有一个循环语句。使用NumPy函数摆脱循环，并确保你的代码性能因此而得到提升。

方法1:

import numpy as np
from matplotlib.pyplot import plot, show

# 生成 t 值
t = np.linspace(-np.pi, np.pi, 201)

# 定义 k
k = np.arange(1, 99)
k = 2 * k - 1  # 只保留奇数

# 使用 NumPy 向量化操作代替循环
# 对于每个 t[i]，我们不需要逐个处理，而是直接通过广播机制计算整个数组
f = np.sum(np.sin(k[:, None] * t) / k[:, None], axis=0)  # k[:, None] 用于扩展 k 维度

# 归一化
f = (4 / np.pi) * f

# 绘制图形
plot(t, f)
show()

方法2:vectorize函数

import numpy as np
from matplotlib.pyplot import plot, show

# 生成 t 值
t = np.linspace(-np.pi, np.pi, 201)

# 定义 k
k = np.arange(1, 99)
k = 2 * k - 1  # 只保留奇数

# 使用 np.vectorize 创建一个函数来替代循环
def compute_f(t_value):
    return np.sum(np.sin(k * t_value) / k)

# 使用 np.vectorize 将上面的函数向量化
vectorized_compute_f = np.vectorize(compute_f)

# 计算 f 数组
f = vectorized_compute_f(t)

# 归一化
f = (4 / np.pi) * f

# 绘制图形
plot(t, f)
show()

5.19 锯齿波和三角波

在示波器上，锯齿波和三角波也是常见的波形。和方波类似，我们也可以将它们表示成无穷傅里叶级数。对锯齿波取绝对值即可得到三角波。锯齿波的无穷级数表达式如下：

\[\sum_{k = 1}^{\infty} \frac{-2 sin (2 \pi kt)}{k \pi}\]

5.20 动手实践：绘制锯齿波和三角波

与前面的教程中一样，我们仍将以相同的方式初始化t和k。同样，取k=99以保证足够的精度。绘制锯齿波和三角波的步骤如下。

(1) 将函数值初始化为0：

import sys
from matplotlib.pyplot import plot
from matplotlib.pyplot import show

t = np.linspace(-np.pi, np.pi, 201)
# k = np.arange(1, float(sys.argv[1]))
k = np.arange(1, 99)
f = np.zeros_like(t)

(2) 同样，直接使用sin和sum函数进行计算：

for i in range(len(t)): 
    f[i] = np.sum(np.sin(2 * np.pi * k * t[i])/k)  
f = (-2 / np.pi) * f 

(3) 同时绘制锯齿波和三角波并不难，因为三角波函数的取值恰好是锯齿波函数值的绝对值。使用如下代码绘制波形：

plot(t, f, lw=1.0)  
plot(t, np.abs(f), lw=2.0)  
show() 

在下图中，较粗的曲线为三角波。

刚才做了些什么 : 我们使用sin函数绘制了锯齿波。函数的输入值是用linspace产生的，而一组k值是用arange函数生成的。三角波则是对锯齿波取绝对值得到的。

完整代码如下：

import sys
from matplotlib.pyplot import plot
from matplotlib.pyplot import show

t = np.linspace(-np.pi, np.pi, 201)
# k = np.arange(1, float(sys.argv[1]))
k = np.arange(1, 99)
f = np.zeros_like(t)
for i in range(len(t)): 
    f[i] = np.sum(np.sin(2 * np.pi * k * t[i])/k)  
f = (-2 / np.pi) * f 

plot(t, f, lw=1.0)  
plot(t, np.abs(f), lw=2.0)  
show() 

勇敢出发：摆脱循环语句

你是否愿意接受挑战，摆脱代码中的循环语句？使用NumPy函数应该可以完成这个任务，并且代码性能可以翻倍。

与5.18节的勇敢出发类似，不再赘述。

5.21 位操作函数和比较函数

位操作函数可以在整数或整数数组的位上进行操作，它们都是通用函数。^、&、|、<<、>>等位操作符在NumPy中也有对应的部分，<、>、==等比较运算符也是如此。有了这些操作符，你可以在代码中玩一些高级技巧以提升代码的性能。不过，它们会使代码变得难以理解，因此需谨慎使用。

5.22 动手实践：玩转二进制位

我们将学习三个小技巧——检查两个整数的符号是否一致，检查一个数是否为2的幂数，以及计算一个数被2的幂数整除后的余数。我们会分别展示两种方法，即使用位操作符和使用相应的NumPy函数。

(1) 第一个小技巧需要用XOR或者^操作符。XOR操作符又被称为不等运算符，因此当两个操作数的符号不一致时，XOR操作的结果为负数。在NumPy中，^操作符对应于bitwise_xor函数，<操作符对应于less函数。

x = np.arange(-9, 9)  
y = -x 
print("Sign different?", (x ^ y) < 0 )
print("Sign different?", np.less(np.bitwise_xor(x, y), 0))

不出所料，除了都等于0的情况，所有整数对的符号均相异。

(2) 在二进制数中，2的幂数表示为一个1后面跟一串0的形式，例如10、100、1000等。而比2的幂数小1的数表示为一串二进制的1，例如11、111、1111（即十进制里的3、7、15）等。如果我们在2的幂数以及比它小1的数之间执行位与操作AND，那么应该得到0。NumPy中，&操作符对应于bitwise_and函数，==操作符对应于equal函数。

print( "Power of 2?\n", x, "\n", (x & (x - 1)) == 0)
print( "Power of 2?\n", x, "\n", np.equal(np.bitwise_and(x, (x - 1)), 0) )

(3) 计算余数的技巧实际上只在模为2的幂数（如4、8、16等）时有效。二进制的位左移一位，则数值翻倍。在前一个小技巧中我们看到，将2的幂数减去1可以得到一串1组成的二进制数，如11、111、1111等。这为我们提供了掩码（mask），与这样的掩码做位与操作AND即可得到以2 的幂数作为模的余数。在NumPy中，<<操作符对应于left_shift函数。

print("Modulus 4\n", x, "\n", x & ((1 << 2) - 1) )
print("Modulus 4\n", x, "\n", np.bitwise_and(x, np.left_shift(1, 2) - 1) )

刚才做了些什么 : 我们学习了三个运用位操作的小技巧——检查两个整数的符号是否一致，检查一个数是否为2的幂数，以及计算一个数被2的幂数整除后的余数。我们看到了NumPy中对应于^、&、<<、<等操作符的通用函数。

5.23 本章小结

在本章中，我们学习了NumPy中的矩阵和通用函数，包括如何创建矩阵以及通用函数的工作方式。我们还简单介绍了算术运算函数、三角函数、位操作函数和比较函数等通用函数。

下一章中，我们将开始学习NumPy模块的相关内容。

numpy(5)矩阵和通用函数

第 5 章 NumPy矩阵和通用函数

5.1 矩阵

5.2 动手实践：创建矩阵

5.3 从已有矩阵创建新矩阵

5.4 动手实践：从已有矩阵创建新矩阵

突击测验：使用字符串定义矩阵

5.5 通用函数

5.6 动手实践：创建通用函数

5.7 通用函数的方法

5.8 动手实践：在 add 上调用通用函数的方法

5.9 算术运算

5.10 动手实践：数组的除法运算

勇敢出发：尝试`future.division`

5.11 模运算

5.12 动手实践：模运算

5.13 斐波那契数列

5.14 动手实践：计算斐波那契数列

勇敢出发：分析计算耗时

5.15 利萨茹曲线

5.16 动手实践：绘制利萨茹曲线

5.17 方波

5.18 动手实践：绘制方波

勇敢出发：摆脱循环语句

5.19 锯齿波和三角波

5.20 动手实践：绘制锯齿波和三角波

勇敢出发：摆脱循环语句

5.21 位操作函数和比较函数

5.22 动手实践：玩转二进制位

5.23 本章小结

CATALOG

FEATURED TAGS

FRIENDS

5.1 矩阵

5.2 动手实践：创建矩阵

5.3 从已有矩阵创建新矩阵

5.4 动手实践：从已有矩阵创建新矩阵

突击测验 ： 使用字符串定义矩阵

5.5 通用函数

5.6 动手实践：创建通用函数

5.7 通用函数的方法

5.8 动手实践：在 add 上调用通用函数的方法

5.9 算术运算

5.10 动手实践：数组的除法运算

勇敢出发 ： 尝试__future__.division

5.11 模运算

5.12 动手实践：模运算

5.13 斐波那契数列

5.14 动手实践：计算斐波那契数列

勇敢出发：分析计算耗时

5.15 利萨茹曲线

5.16 动手实践：绘制利萨茹曲线

5.17 方波

5.18 动手实践：绘制方波

勇敢出发：摆脱循环语句

5.19 锯齿波和三角波

5.20 动手实践：绘制锯齿波和三角波

勇敢出发：摆脱循环语句

5.21 位操作函数和比较函数

5.22 动手实践：玩转二进制位

5.23 本章小结

CATALOG

FEATURED TAGS

FRIENDS

突击测验：使用字符串定义矩阵

勇敢出发：尝试`future.division`