MATALB大数的加减乘除幂运算
 
一、基本原理
二、MATLAB简单实现
1. 加（+）
2. 减（-）
3. 乘（*）
4. 除（/）
5. 幂（^）
  
  
三、测试
四、关于优化
五、其他

 
by HPC_ZY
 
最近需要做大数计算，且要求结果精度与原始数据一致（不能使用科学计数等近似解）。所以总结分享一些“超大正整数的基本运算方法”，如加、减、乘、除、幂等
 
一、基本原理
 
uint32（unsigned int）型可表示的最大数为$(2^{32}-1)$，uint64（unsigned long long）型可表示的最大数为$(2^{64}-1)$，超过这样的数就没法直接计算了。
 但是我们知道加减乘除的本质还是按位计算，所以我们可以将大数拆分为多个小位数。
 根据这个原理最合适的就是使用字符串的形式，这样就将上限从最大数据类型变成了计算机最大内存。
 下面我们就来研究字符串格式的加减乘除，首先分享MATLAB中的实现方法，然后改写成C。
 
注：本文仅研究正整数，暂不涉及负数与小数
 
二、MATLAB简单实现
 
1. 加（+）
 
任意大数相加，主要原理：
 1）按位相加
 2）满十进一
 
%% 大数加（输入为字符串格式）
function out = OLNadd(a,b)

% 转数字矩阵(为了逻辑方便，我们将数字反序）
mata = a(end:-1:1)-'0';
matb = b(end:-1:1)-'0';

% 计算位数
la = length(a);
lb = length(b);

% 初始化
if la>=lb
    digitStart = lb;
    matout = [mata,0];
    mattmp = matb;
else
    digitStart = la;
    matout = [matb,0];
    mattmp = mata;
end

% 循环求解
for k = digitStart:-1:1
    % 位求和，满十进一
    tmp = matout(k)+mattmp(k);
    if tmp<10
        matout(k) = tmp;
    else
        matout(k) = tmp-10;
        matout(k+1) = matout(k+1)+1;
    end
end

% 转字符
out = char(matout(end:-1:1)+'0');
if out(1) == '0'
    out(1) = [];
end

end
 
 
2. 减（-）
 
任意大数相减，主要原理：
 1）比较大小，以大减小
 2）符号由大小关系决定
 3）按位相减
 4）不够向前借位
 
%% 大数减（输入为字符串格式）
function out = OLNsub(a,b)

% 转数字矩阵(为了逻辑方便，我们将数字反序）
mata = a(end:-1:1)-'0';
matb = b(end:-1:1)-'0';

% 计算位数
la = length(a);
lb = length(b);

% 初始化
isEqual = 0;
if la>lb % 默认用大的减小的（正负号由大小关系决定）
    digitEnd = lb;
    matout = [mata,0];
    mattmp = matb;
    flag = 1; % 结果为正
elseif la<lb
    digitEnd = la;
    matout = [matb,0];
    mattmp = mata;
    flag = 0;   % 结果为负
else
    % 对比第一个不相等的位数，谁大谁小
    idx = find(a~=b,1);
    if ~isempty(idx)
        if a(idx)>b(idx)
            digitEnd = lb;
            matout = [mata,0];
            mattmp = matb;
            flag = 1; % 结果为正
        else
            digitEnd = la;
            matout = [matb,0];
            mattmp = mata;
            flag = 0;   % 结果为负
        end
    else % 如果为空，说明两个数完全一致
        isEqual = 1;
    end
end

if ~isEqual 
    % 循环求解
    for k = 1:digitEnd
        % 位求差，不足借一
        if matout(k)>=mattmp(k)
            matout(k) = matout(k)-mattmp(k);
        else
            matout(k) = 10+matout(k)-mattmp(k);
            borrow = 1;
            while 1 % 不够减时往前借，直到借到为止。
                if matout(k+borrow)>=1 % 够借ok
                    matout(k+borrow) = matout(k+borrow)-1;
                    break
                else % 不够，再往前
                    matout(k+borrow) = 10+matout(k+borrow)-1;
                    borrow = borrow+1;
                end
            end
        end
    end
    % 转字符
    out = char(matout(end:-1:1)+'0');
    while out(1)=='0'      
            out(1) = [];
    end
    if ~flag
        out = ['-',out];
    end
else
    out = '0';
end

end
 
 
3. 乘（*）
 
任意大数相乘，主要原理：
 1）按位求积
 2）结果累加
 
%% 大数乘（输入为字符串格式）
function out = OLNmult(a,b)

% 转数字矩阵(为了逻辑方便，我们将数字反序）
mata = a(end:-1:1)-'0';
matb = b(end:-1:1)-'0';

% 计算位数
la = length(a);
lb = length(b);

% 初始化
matout = zeros(1,la+lb);

% 循环求解
for i = 1:la
    if mata(i)==0
        continue
    end
    for j = 1:lb
        if matb(j)==0
            continue
        end       
        % 位求积，分高低位
        valtmp = mata(i)*matb(j);
        valup = floor(valtmp/10);
        vallow = valtmp-valup*10;
        % 低位求和
        idx = i+j-1;
        val = matout(idx)+vallow;
        if val<10
            matout(idx) = val;
        else
            matout(idx) = val-10;
            matout(idx+1) = matout(idx+1)+1;
        end
        % 求高位
        idx = i+j;
        val = matout(idx)+valup;
        if val<10
            matout(idx) = val;
        else
            matout(idx) = val-10;
            matout(idx+1) = matout(idx+1)+1;
        end
    end
end

% 转字符
out = char(matout(end:-1:1)+'0');
while out(1)=='0'
    out(1) = [];
end

if ~flag
    out = ['-',out];
end

end
 
 
4. 除（/）
 
被除数任意大，除数$<99999999$。
 （原因是在除法中，被除数符合分配律，除数不符合。对于超过类型允许的除数，我不知道怎么处理。所以采取的思路是将被除数拆分为类型允许的大小，即$<2^{32}-1=4294967295$（十位数），拆分的被除数最大只能为$999999999$（九位数），进而得到除数最多为$99999999$（八位数），主要原理：
 1）按位求商，不足移位（即用最高位/除数，不够除则用最高两位，还不够用最高三位，以此类推）
 2）高位余数，加至低位
 
%% 大数除（输入为字符串格式）
function [out,remainder] = OLNdiv(a,b)

% 转数字矩阵（为了逻辑方便，我们将数字反序）
mata = a(end:-1:1)-'0';
matb = b(end:-1:1)-'0';

% 计算位数
la = length(a);
lb = length(b);

% 转数字（按四位分组）
numa = mata(end:-1:1);
numb = 0;
for k = 1:lb
    numb = numb+matb(k)*10^(k-1);
end

% 分组除保存结果
res = zeros(1,la); %
remainder = 0; % 每次计算的余数
for i = 1:la
    divisor = numa(i)+remainder*10; % 被除数 = 上一组的余数*10000+这一组
    remainder = rem(divisor,numb);
    res(i) = (divisor-remainder)/numb; 
end

% 转字符
remainder = num2str(remainder);
out = char(res+'0');
while out(1)=='0'
    out(1) = [];
end
if ~flag
    out = ['-',out];
end

end
 
 
5. 幂（^）
 
底数指数任意大。主要原理：
 1）循环大数乘
 
%% 大数幂（输入为字符串格式）
function out = OLNpow(a,b)

% 转数字矩阵（为了逻辑方便，我们将数字反序）
matb = b(end:-1:1)-'0';

% 计算位数
lb = length(matb);

% 初始化
n = -1;
for k = 1:lb
n = n+matb(k)*10^(k-1);
end
% 循环乘（直接利用大数乘法，偷下懒）
tmp = a;
for k = 1:n 
   tmp = OLNmult(tmp,a);
end
out = tmp;

end
 
 
三、测试
 
准确性测试
 首先测试各个算法结果是否正确，我们用int64能支持的数字进行，如下
 
%% 加
ca = '68924355';
cb = '411533';

cc = OLNadd(ca,cb); % 字符计算
nc = str2double(ca)+str2double(cb); % 数字计算

disp('加：')
disp([ca,' + ', cb,' = ',cc])
disp(['标准答案 = ',num2str(nc)])
disp(' ')


%% 减
% 基本测试
ca = '985355';
cb = '411533';

cc = OLNsub(ca,cb); % 字符计算
nc = str2double(ca)-str2double(cb); % 数字计算

disp('减-基本测试：');
disp([ca,' - ', cb,' = ',cc])
disp(['标准答案 = ',num2str(nc)])
disp(' ')

% 小减大测试
ca = '85355';
cb = '411533';

cc = OLNsub(ca,cb); % 字符计算
nc = str2double(ca)-str2double(cb); % 数字计算

disp('减-小减大测试：')
disp([ca,' - ', cb,' = ',cc])
disp(['标准答案 = ',num2str(nc)])
disp(' ')

% 借位测试
ca = '100000';
cb = '1';

cc = OLNsub(ca,cb); % 字符计算
nc = str2double(ca)-str2double(cb); % 数字计算

disp('减-借位测试：')
disp([ca,' - ', cb,' = ',cc])
disp(['标准答案 = ',num2str(nc)])
disp(' ')


%% 乘
ca = '996655';
cb = '872233';

cc = OLNmult(ca,cb); % 字符计算
nc = str2double(ca)*str2double(cb); % 数字计算

disp('乘：')
disp([ca,' * ', cb,' = ',cc])
disp(['标准答案 = ',num2str(nc)])
disp(' ')


%% 除
% 整除测试
ca = '869315380615';
cb = '872233';

[cc,cr] = OLNdiv(ca,cb); % 字符计算
nr = rem(str2double(ca),str2double(cb)); % 余数计算
nc = (str2double(ca)-nr)/str2double(cb); % 数字计算

disp('除-整除测试：')
disp([ca,' / ', cb,' = ', cc,' 余 ',cr])
disp(['标准答案 = ',num2str(nc),' 余 ',num2str(nr)])
disp(' ')

% 非整除测试
ca = '88772211996655';
cb = '872233';

[cc,cr] = OLNdiv(ca,cb); % 字符计算
nr = rem(str2double(ca),str2double(cb)); % 余数计算
nc = (str2double(ca)-nr)/str2double(cb); % 数字计算

disp('除-非整除测试：')
disp([ca,' / ', cb,' = ', cc,' 余 ',cr])
disp(['标准答案 = ',num2str(nc),' 余 ',num2str(nr)])
disp(' ')

%% 幂
ca = '16';
cb = '11';

cc = OLNpow(ca,cb); % 字符计算
nc = str2double(ca)^str2double(cb); % 数字计算

disp('幂：')
disp([ca,' ^ ', cb,' = ',cc])
disp(['标准答案 = ',num2str(nc)])
disp(' ')

 
结果如下，完全一致，接下来可以放心测试大数
 
 
 
大数测试
 使用超大数字进行测试
 
%% 加
ca = '22368936111988924355';
cb = '118897773322411533';

cc = OLNadd(ca,cb); % 字符计算

disp('加：')
disp([ca,' + ', cb]),disp([' = ',cc])
disp(' ')


%% 减
ca = '22368936111988924355';
cb = '118897773322411533';

cc = OLNsub(ca,cb); % 字符计算

disp('减：');
disp([ca,' - ', cb]),disp([' = ',cc])
disp(' ')


%% 乘
ca = '789123996655';
cb = '432987872233';

cc = OLNmult(ca,cb); % 字符计算

disp('乘：')
disp([ca,' * ', cb]),disp([' = ',cc])
disp(' ')


%% 除
ca = '12345678988772211996655';
cb = '32872233';

[cc,cr] = OLNdiv(ca,cb); % 字符计算

disp('除')
disp([ca,' / ', cb]),disp([' = ',cc,' 余 ',cr])
disp(' ')


%% 幂
ca = '541';
cb = '25';

cc = OLNpow(ca,cb); % 字符计算

disp('幂：')
disp([ca,' ^ ', cb]),disp([' = ',cc])
disp(' ')
 
结果如下
 
 
 
四、关于优化
 
为了理解方便所以函数写法有累赘的地方，自己使用时可合理优化。
有大佬提出，加减乘除可以在二进制计算中完成。
关于幂运算中，指数太大计算时间会爆表，可以改写为快速幂方法。
我没有加入判断语句，即使输入的不是数字还是会计算，可以加上。
可以对输入加上符号，这样就能进行负数计算。
关于小数计算，加减乘很简单，只要判断小数点位置即可。除法则需要在此基础上加上计算精度（小数后保留几位）。关于幂，底数为小数类似乘法，但指数为小数就要配合开方。
 
五、其他
 
1、该方法的C实现可看另一篇文章《C/C++大数的加减乘除幂运算》。撰写中……
 2、函数与测试代码文中已全部公开，如果嫌懒得复制且已开通会员，可以在这里下载demo

☞ ░ 前往老猿Python博客 https://blog.csdn.net/LaoYuanPython ░
 
一、引言
 
在写该文之前，老猿就图像的一些运算已经单独边学边发了，在写这些文的过程中，发现这些运算函数共同点很多，例如大部分参数一样、部分处理方法一样等，另外还有些函数可以实现相同或近似的效果，因此在前面那些文章的基础上，将其综合为一个整体来发布，更方便大家阅读。
 
OpenCV中图像存储为矩阵，因此图像的运算其实就是矩阵的运算。图像的运算主要包括图像基础算术运算、图像加权运算（又称为图像融合）、按位运算等类别。这些运算可以直接通过numpy矩阵进行，也可以通过opencv的专用方法进行，但opencv的矩阵运算和numpy矩阵运算还是有些不同。例如在加法处理上，OpenCV加法是饱和运算（超过255即按255），而Numpy加法是模运算（超过255按256取模的结果作为结果）。对加法来说，颜色值越大OpenCV 的结果会更好，因此推荐使用opencv的进行算术运算。本文的内容全部基于OpenCV的方法进行介绍。
 
二、OpenCV图像运算语法
 
2.1、图像运算函数列表
 
OpenCV图像运算包括如下函数：
 
加法运算：add(src1, src2, dst=None, mask=None, dtype=None)
减法运算：subtract(src1, src2, dst=None, mask=None, dtype=None)
乘法运算：multiply(src1, src2, dst=None, scale=None, dtype=None)
除法运算：divide(src1, src2, dst=None, scale=None, dtype=None)
幂运算：pow(src, power, dst=None)
开方运算：sqrt(src, dst=None)
自然常数e为底的指数函数：exp(src, dst=None)
对数运算：log(src, dst=None)
融合权重加法：addWeighted(src1, alpha, src2, beta, gamma, dst=None, dtype=None)
位与运算：bitwise_and(src1, src2, dst=None, mask=None)
位或运算：bitwise_or(src1, src2, dst=None, mask=None)
位异或运算：bitwise_xor(src1, src2, dst=None, mask=None)
位非运算：bitwise_not(src, dst=None, mask=None)
 
2.2、图像运算常用参数说明
 
 src：输入图像矩阵
 src1、src2：两副大小和通道数相等的输入图像或一副输入图像和一个标量（关于标量请参考《OpenCV-Python中的标量Scalar是什么》）
 dst：目标图像输出，要求与输入图像大小相同，如果传值则可以直接以实参作为目标图像存储变量，否则可以用函数返回值作为目标图像存储变量，在实参非None传入的情况下，返回值与该实参值相同
 scale：缩放因子，图像运算时，先执行src1*scale，再以该乘积进行后续运算
 mask：图像掩膜，可选参数，为8位单通道的灰度图像，用于指定要更改的输出图像数组的元素，即输出图像像素只有mask对应位置元素不为0的部分才输出，否则该位置像素的所有通道分量都设置为0。更多关于掩膜的内容请参考《OpenCV学习02-矩阵的掩膜操作》或《Opencv图像处理之详解掩膜mask》
 dtype：可选参数，输出图像数组的深度，即图像单个像素值的位数（如RGB用三个字节表示，则为24位）。在帮助文档中介绍，参数src1和src2可以有不同的图像深度（即图像像素位数，如8位、16位、24位和32位），如可以将16位图像和一个8位图像相加将输出结果保存在32位输出数组中。关于这个参数老猿研究可很长时间，因为认为如果图像深度不一样，意味着图像的通道数也不一样，通道数不一样，意味着数组大小不一样，这样的两个数组无法进行运算。后来机缘契合下自认为理解了，这里说的图像深度不是通道数*8，而是单个通道值表示的位数，即单通道不一定是8位的，也可以是16位、24位或32位的。看如下代码就明白了：
 
def main():
    img1 = cv2.imread(r'F:\pic\shape1.png').astype(np.float32)
    img2 = cv2.imread(r'F:\pic\shape2.png')
    img = cv2.add(img1,img2,dtype=24)
 
上述代码将以两种不同表示方法读入两副图像，第一幅图像是以float32来表示图像单通道值，第二幅图像是缺省值uint8来表示图像单通道值，二者机器位数不同，但相加之后转为了24位图像，即单通道为8字节图像。
 
注：
 
上面介绍的dtype取值说明总体是对的，但dtype=24并不表示24位，具体含义请参考《OpenCV-Python图像像素位深表示法以及应用：https://blog.csdn.net/LaoYuanPython/article/details/109574736》。
 
2.3、部分图像运算函数详解
 
针对部分重要的图像运算，老猿在前面已经单独进行了介绍，包括：
 
OpenCV-Python图像的加法运算cv2.add函数详解
OpenCV-Python图像的减法运算cv2.subtract函数详解以及和矩阵减法的差异对比
OpenCV-Python图像乘法运算cv2.multiply函数详解及像素值溢出归一化处理
OpenCV-Python图像除法运算cv2.divide函数及图像相除处理
OpenCV-Python图像融合cv2.addWeighted权重加法函数详解
OpenCV-Python图像位与运算bitwise_and函数详解
 
2.4、其他图像运算函数简介
 
 
幂运算：pow(src, power, dst=None)，对图像的每个通道值计算power参数对应的幂作为结果图像的通道值，如果power为整数，则直接计算幂值，如果power为浮点数，则取通道值的绝对值参与计算，即：
 
 
开方运算：sqrt(src, dst=None)，对图像的每个通道值开方作为结果图像的通道值，即：dst(I)=sqrt（src1(I)）
 
 
自然常数e为底的指数函数：exp(src, dst=None)，以e为底对图像的每个通道值作为幂值计算结果图像的通道值，即：
 
 
对数运算：log(src, dst=None)，计算图像的每个通道值的自然对数作为结果图像的通道值，即：dst(I)=ln(src(I))
 
 
位或运算：bitwise_or(src1, src2, dst=None, mask=None)，计算两副图像每个通道值或一图像通道值与一个标量的按位或的结果作为结果图像的通道值
 
 
位异或运算：bitwise_xor(src1, src2, dst=None, mask=None)，计算两副图像每个通道值或一图像通道值与一个标量的按位异或的结果作为结果图像的通道值
 
 
位非运算：bitwise_not(src, dst=None, mask=None)，将src图像的每个通道值按位取反作为结果图像的值。
 
 
三、相关函数的作用分析
 
图像的减法、除法以及异或都可以用于分析图像的差异点，但减法和异或更准确 ；
加法、权重加法以及图像或运算都可以用于合并图像，但由于运算方法与差异，效果会有不同，针对不同图像的融合选择哪种方法与图像数据及应用的目标相关；
图像的自乘、幂运算都可以用于调整图像的对比度；
非运算及异或运算可以将图像的一部分图像的视觉效果提高；
图像乘法、位与都可以用于提取图像的感兴趣部分或者屏蔽某些部分，在这方面的功能与图像掩码的效果相同，只是实现方式不同；
通过阈值处理、图像求反、与处理等可以提取图像中的结构特征
 
四、按位运算的简单案例
 
前面部分图像运算函数详解中举例介绍了图像加减乘除权重加及位与的相关功能，在此补充一个简单的按位运算的案例。代码如下：
 
import numpy as np
import cv2
def main():
    img1 = cv2.imread(r'F:\pic\shape1.png').astype(np.float32)
    img2 = cv2.imread(r'F:\pic\shape2.png')

    resultImgAnd = cv2.bitwise_and(img1, img2)
    resultImgOr = cv2.bitwise_or(img1, img2)
    resultImgXor = cv2.bitwise_xor(img1, img2)
    resultImgNot = cv2.bitwise_not(img1)
    resultImgXorScalar = cv2.bitwise_xor(img1, (255,255,255,255))

    cv2.imshow('img1',img1)
    cv2.imshow('img2', img2)
    cv2.imshow('resultImgAnd', resultImgAnd)
    cv2.imshow('resultImgOr', resultImgOr)
    cv2.imshow('resultImgXor', resultImgXor)
    cv2.imshow('resultImgNot', resultImgNot)
    cv2.imshow('resultImgXorScalar', resultImgXorScalar)
    cv2.waitKey(0)

main()
 
运行显示的图片截屏：
 
 
五、小结
 
本文详细介绍了OpenCV-Python图像的加减乘除幂开方对数及位运算相关的函数及语法，并总结了相关函数的作用。OpenCV中图像存储为矩阵，因此图像的运算其实就是矩阵的运算。图像的运算主要包括图像基础算术运算、图像加权运算（又称为图像融合）、按位运算等类别。这些运算可以直接通过numpy矩阵进行，也可以通过opencv的专用方法进行，但opencv的矩阵运算是饱和运算，其运算效果比纯粹的矩阵运算效果更好。
 
如果觉得本文可以给您带来帮助，请大家帮忙点个赞、加个收藏，谢谢！
 
更多OpenCV-Python介绍请参考专栏《OpenCV-Python图形图像处理 》
 专栏网址：https://blog.csdn.net/laoyuanpython/category_9979286.html
 
关于老猿的付费专栏
 
老猿的付费专栏《使用PyQt开发图形界面Python应用 》(https://blog.csdn.net/laoyuanpython/category_9607725.html)专门介绍基于Python的PyQt图形界面开发基础教程，付费专栏《moviepy音视频开发专栏》 （https://blog.csdn.net/laoyuanpython/category_10232926.html)详细介绍moviepy音视频剪辑合成处理的类相关方法及使用相关方法进行相关剪辑合成场景的处理，两个专栏都适合有一定Python基础但无相关知识的小白读者学习。
 
付费专栏文章目录：《moviepy音视频开发专栏文章目录》(https://blog.csdn.net/LaoYuanPython/article/details/107574583）、《使用PyQt开发图形界面Python应用专栏目录 》（https://blog.csdn.net/LaoYuanPython/article/details/107580932）。
 
对于缺乏Python基础的同仁，可以通过老猿的免费专栏《专栏：Python基础教程目录》（https://blog.csdn.net/laoyuanpython/category_9831699.html）从零开始学习Python。
 
如果有兴趣也愿意支持老猿的读者，欢迎购买付费专栏。
 
 
跟老猿学Python、学OpenCV!
 
☞ ░ 前往老猿Python博文目录 https://blog.csdn.net/LaoYuanPython ░

Python/Numpy之加减乘除开方平方幂值三角函数等算数运算
 
 
 
import numpy as np
 a = np.array([1,2,3])
 b = np.array([2,3,4])
 
 
加法：
 add1 = a + b
 add2 = np.add(a,b)
减法：
 sub1 = a - b
 sub2 = np.subtract(a,b)
乘法：
 mul1 = a * b
 mul2 = np.multiply(a,b)
除法：
 div1 = a / b
 div2 = np.divide(a,b)
幂指运算开方、平方、指数、对数：
 sqrtx = np.sqrt(a) # 开方
 powerx = np.power(a,2) # 平方
 expx = np.exp(a)
 logx = np.log(a)
三角函数：
 sinx = np.sin(a)
 cosx = np.cos(a)
 
>>> import numpy as np
>>> a = np.array([1,2,3])
>>> b = np.array([2,3,4])
# 减法运算
>>> sub1 = a - b
>>> sub1
array([-1, -1, -1])
>>> sub2 = np.subtract(a,b)
>>> sub2
array([-1, -1, -1])
# 加法运算
>>> add1 = a + b
>>> add1
array([3, 5, 7])
>>> add2 = np.add(a,b)
>>> add2
array([3, 5, 7])
# 除法运算
>>> div1 = a / b
>>> div1
array([0.5       , 0.66666667, 0.75      ])
>>> div2 = np.divide(a,b)
>>> div2
array([0.5       , 0.66666667, 0.75      ])
# 乘法运算
>>> mul1 = a * b
>>> mul1
array([ 2,  6, 12])
>>> mul2 = np.multiply(a,b)
>>> mul2
array([ 2,  6, 12])
# 三角函数
>>> sinx = np.sin(a)
>>> sinx
array([0.84147098, 0.90929743, 0.14112001])
>>> cosx = np.cos(a)
>>> cosx
array([ 0.54030231, -0.41614684, -0.9899925 ])
>>> tanx = np.tan(a)
>>> tanx
array([ 1.55740772, -2.18503986, -0.14254654])
# 幂指数运算
>>> expx = np.exp(a)
>>> expx
array([ 2.71828183,  7.3890561 , 20.08553692])
>>> logx = np.log(a)
>>> logx
array([0.        , 0.69314718, 1.09861229])
# 开方运算
>>> sqrtx = np.sqrt(a)
>>> sqrtx
array([1.        , 1.41421356, 1.73205081])
# 平方运算
>>> powerx = np.power(a,2)
>>> powerx
array([1, 4, 9], dtype=int32)
 
 
 
三角函数采用的是弧度制，如np.sin(np.pi/2)=1
 更多Python数据科学相关见专题Python数据科学技能索引
 

运算公式
 
（这将是百度上关于复数运算较为全面的一篇）
 
加法
 
 
 (a+bi)+(c+di)=(a+c)+(b+d)i
 
 
减法
 
 
 (a+bi)+(c+di)=(a-c)+(b-d)i
 
 
乘法
 
 
 (a+bi)(c+di)=(ac-bd)+(bc+ad)i
 
 
除法
 
 
 (a+bi)/(c+di)=(ac+bd)/(c ^ 2   + d ^ 2) +((bc-ad)/(c ^ 2   + d ^ 2)) i
 
 
exp
 
 
log
 

 对于其它对数计算可使用换低公式：
 
 
sin & cos
 

 
 
幂
 
 
结果检验
 
 
代码
 
（注意参数可能为空指针，由于使用浮点数作为基础数据类型，对于复数为零的判断要设置一个阈值，这里使用了10倍的Double最小值）
 
package complex;
import java.util.Objects;
public class Complex {
    private double re;   //real
    private double im;   //imaginary
    static final Complex I=new Complex(0,1);
    static final Complex NI=new Complex(0,-1);
    public Complex(double real, double imag) {
        re = real;
        im = imag;
    }
    public static Complex add(Complex a,Complex b) {
    	Objects.requireNonNull(a);
    	Objects.requireNonNull(b);
        return new Complex(a.re + b.re,a.im + b.im);
    }
    public static Complex sub(Complex a,Complex b) {
    	Objects.requireNonNull(a);
    	Objects.requireNonNull(b);
        return new Complex(a.re - b.re,a.im - b.im);
    }
    public static Complex mul(Complex a,Complex b) {
    	Objects.requireNonNull(b);
    	Objects.requireNonNull(b);
        double real = a.re * b.re - a.im * b.im;
        double imag = a.re * b.im + a.im * b.re;
        return new Complex(real,imag);
    }
    public static Complex div(Complex a,Complex b) {
    	Objects.requireNonNull(b);
    	Objects.requireNonNull(b);
        if(b.isZero()) {
        	throw new ArithmeticException("/ by zero");
        }
        double den=b.re * b.re + b.im * b.im;
        double real = a.re * b.re + a.im * b.im;
        double imag = a.im * b.re - a.re * b.im;
        return new Complex(real/den, imag/den);
    }
    public static Complex log(Complex com) {
    	double real = com.re * com.re + com.im * com.im;
    	double imag = Math.atan(com.im * com.re );
        return new Complex(Math.log(real)/2,imag);
    }
    public static  Complex exp(Complex com) {
    	double real = Math.cos(com.im);
    	double imag = Math.sin(com.im);
    	double expx= Math.exp(com.re);
        return new Complex(expx*real,expx*imag); 
    }
    public static  Complex sin(Complex com) {
    	final Complex cf=new Complex(0,2);//coefficient
    	Complex e1=exp(mul(I,com));
    	Complex e2=exp(mul(NI,com));
        return div(sub(e1,e2),cf); 
    }
    public static  Complex cos(Complex com) {
    	final Complex cf=new Complex(0,2);//coefficient
    	Complex e1=exp(mul(I,com));
    	Complex e2=exp(mul(NI,com));
        return div(add(e1,e2),cf); 
    }
    public static  Complex pow(Complex a,Complex b) {;
        return Complex.exp( Complex.mul(b, Complex.log(a))); 
    }
    public boolean isZero() {
    	if(Math.abs(re)<10*Double.MIN_VALUE
    			&&Math.abs(im)<10*Double.MIN_VALUE) {
    		 return true;
    	}
        return false;
    }
    public Complex conjugate() {
        return new Complex(re,-im);
    }
    public double length() {
        return Math.sqrt(re * re + im * im);
    }
    public String  toString() {
    	if(im<0)
    		return "( "+re+im+"i )";
        return "( "+re+"+"+im+"i )";
    }
}
 
测试数据
 
    public static void main(String[] args) {
    	Complex a=new Complex(1,2);
    	Complex b=new Complex(0,1);
    	System.out.println("add "+a+" + " +b +" = "+add(a,b));
    	System.out.println("sub "+a+" - " +b +" = "+sub(a,b));
    	System.out.println("mul "+a+" * " +b +" = "+mul(a,b));
    	System.out.println("div "+a+" / " +b +" = "+div(a,b));
    	Complex lg=new Complex(1,1);
    	System.out.println("log "+lg+" = "+log(lg));
    	System.out.println("value :"+0.5*Math.log(2)+" "+Math.PI/4+"i");
    	
    	Complex ex=new Complex(2,1);
    	System.out.println("exp "+ex+" = "+exp(ex));
    	System.out.println("value :"+(Math.exp(2)*Math.cos(1))+" "+(Math.exp(2)*Math.sin(1))+"i");
    	
    	Complex sinc=new Complex(0,1);
    	System.out.println("sin "+sinc+" = "+sin(sinc));
    	System.out.println("value :"+0+" "+(Math.exp(1)-Math.exp(-1))/2+"i");
    	
    	Complex cosc=new Complex(0,1);
    	System.out.println("cos "+sinc+" = "+cos(cosc));
    	System.out.println("value :"+0+" "+(Math.exp(1)+Math.exp(-1))/2+"i");
    	
    	Complex p1=new Complex(1,1);
    	Complex p2=new Complex(0,1);
    	System.out.println("pow  "+p1+","+p2 +" = "+pow(p1,p2));
    	System.out.println("value :"+0.428829006+" "+0.154871752+"i");
    }
 
add ( 1.0+2.0i ) + ( 0.0+1.0i ) = ( 1.0+3.0i )
 sub ( 1.0+2.0i ) - ( 0.0+1.0i ) = ( 1.0+1.0i )
 mul ( 1.0+2.0i ) * ( 0.0+1.0i ) = ( -2.0+1.0i )
 div ( 1.0+2.0i ) / ( 0.0+1.0i ) = ( 2.0-1.0i )
 log ( 1.0+1.0i ) = ( 0.34657359027997264+0.7853981633974483i )
 value :0.34657359027997264 0.7853981633974483i
 exp ( 2.0+1.0i ) = ( 3.992324048441272+6.217676312367968i )
 value :3.992324048441272 6.217676312367968i
 sin ( 0.0+1.0i ) = ( 0.0+1.1752011936438014i )
 value :0 1.1752011936438014i
 cos ( 0.0+1.0i ) = ( 0.0-1.5430806348152437i )
 value :0 1.5430806348152437i
 pow ( 1.0+1.0i ),( 0.0+1.0i ) = ( 0.4288290062943678+0.15487175246424675i )
 value :0.428829006 0.154871752i