C++开源算法库OpenSAL(Open Standardized Algorithm Library)

http://www.oschina.net/code/OpenSAL

包括了《算法导论》中的几乎所有数据结构和算法（标准库中已有的、不通用的或太简单的除外）。包含算法导论中数据结构：一般堆、二项堆、斐波那契堆、红黑树、通用散列（采用全域散列和完全散列技术）、不相交集合；包含算法导论中的算法：15个常用图论算法、20多个常用代数方面的算法、若干其他算法。包含自己发明的一个编程技术（任意维数组）、一个数据结构（多维对称数组）、一个算法（快速方幂和算法）；该算法库采用安全的智能指针技术，并且尽量使用了泛型编程。图论算法（兼容有向图，无向图）包括：广度和深度优先遍历、确定图是否存在回路、拓扑排序、强连通分支、欧拉环（欧拉路径）、最小生成树（Kruskal、Prim）、单源最短路径（3种）、每对顶点间最短路径（2种）、最大流（2种）等等。代数算法：方幂和、霍纳法则计算多项式和、矩阵乘法（2种)、方阵的LUP分解、解线性方程组（2种）、矩阵求逆（2种）、求伪逆矩阵（2种）、解正态方程组（2种）、最小二乘估计（2种）、快速傅里叶变换、快速傅里叶逆变换、多维快速傅里叶变换、多维快速傅里叶逆变换、快速向量求卷积（单变量多项式乘积）、快速张量求卷积（多变量多项式乘积）等等。及最长公共子序列、简单求大质数、随机实数、键值分离排序等其他算法。欢迎算法爱好者批评指正。

][C/C++]部分代码   

 

001/*//Copyright Statements:

002    This source code file is completed as part of the Open SAL(Open Standardized Algorithm Library) by Ming Li during the summer vocation of 2013 at College of Information Science and Engineering, Central South University. All rights reserved. No part of these files may be reproduced or transmitted in any form or by any means without permission from the author. While, you are permitted to download this Library and its instruction manual freely at http://www.oschina.net/code/snippet_125

0039068_24335# and use it by yourself. If you find some bugs or problems of this Library,you are welcomed to email me([e-mail]). If you want to know any details or updates about this Library or would like to discuss issues about algorithms with enthusiasts,you can apply to join the QQ group(222576709).

004*///

005 

006//此文件（Algorithms_Lmtc.cpp）对本算法库中大部分算法进行测试。

007//通过浏览此文件可以大致了解该算法库中的各个数据结构和算法的使用方式。

008 

009#include "stdafx.h"

010//C++类库

011#include <iostream>

012#include <sstream>

013#include <string>

014#include <vector>

015#include <list>

016#include <queue>

017#include <set>

018#include <bitset>

019#include <exception>

020#include <memory>  //含有auto_ptr类，用于安全的管理动态申请的对象（不能管理动态数组，不能用于容器）

021#include <algorithm> //标准库算法

022#include <numeric> //算术型算法

023#include <functional> //函数对象

024 

025//C函数库

026#include <cctype>

027#include <cstddef>

028#include <cstring>

029#include <ctime>

030 

031//自定义头文件

032#include "smartPtr.h"

033#include "myException.h"

034#include "array.h"

035#include "symmetryArray.h"

036#include "coordinateMapping.h"

037#include "algebra.h"

038#include "heap.h"

039#include "redBlackTree.h"

040#include "hash.h"

041#include "numberTheory.h"

042#include "myMath.h"

043#include "sequence.h"

044#include "binomialHeap.h"

045#include "fibonacciHeap.h"

046#include "nonIntersectSet.h"

047#include "graph.h"

048#include "number.h"

049#include "fastFourierTransform.h"

050 

051//作为函数指针参数的若干函数定义

052bool canExpandTo(const int &item,conststd::vector<int> &st){returntrue;}

053void expandTo(const int &item,std::vector<int> &st){st.push_back(item);}

054using std::string;

055bool g(const int &a,constint &b){returna<b;}

056void f(double x){std::cout<<x<<std::endl;}

057unsigned longkeyToNumber(constunsigned long&t){return t;}

058unsigned longstringToNumber(conststd::string &str){

059    unsignedlong t=0;

060    for(std::string::const_iterator c=str.begin();c!=str.end();c++)

061    {  

062        doubletem=256.0*t+(unsigned long)(*c);

063        unsignedlong tem1=(unsignedlong)(tem/lmtc::MAX_PRIME);

064        tem-=(1.0*tem1*lmtc::MAX_PRIME);

065        t=(unsignedlong)tem;

066    }

067    returnt;

068}

069double polynomial(unsigned int i,double x){

070    returnstd::pow(x,(int)i);

071}

072 

073int _tmain(int argc, _TCHAR* argv[])

074{

075    srand((long)time(NULL));

076 

077    ///*对称数组测试

078    std::cout<<"对称数组测试："<<std::endl;

079    std::vector<unsignedint>symArrDimV(4);

080    symArrDimV[0]=7;symArrDimV[1]=8;symArrDimV[2]=6;symArrDimV[3]=9;

081    lmtc::SymmetryArray<int> symArr(4,symArrDimV);

082    symArrDimV[0]=1;symArrDimV[1]=2;symArrDimV[2]=3;symArrDimV[3]=4;

083    symArr(symArrDimV)=12;

084    symArrDimV[0]=4;symArrDimV[1]=2;symArrDimV[2]=3;symArrDimV[3]=1;

085    constlmtc::SymmetryArray<int> symArrConst=symArr;

086    std::cout<<symArrConst(1,2,4,3)<<"  "<<symArr(1,4,2,3)<<"  "<<symArrConst(1,3,2,4)<<"  "<<symArr(1,4,3,2)<<"  "<<symArrConst(symArrDimV)<<"  "<<symArr(symArrDimV)<<std::endl;

087    std::cout<<"7*8*6*9="<<7*8*6*9<<std::endl;

088    std::cout<<"actual memory needed ="<<symArr.size()<<std::endl;

089    std::cout<<std::endl;//*/

090 

091    ///*一般堆测试

092    std::cout<<"一般堆测试："<<std::endl;

093    intarrHeap[]={5,4,3,2,1,6,8,10,-1};

094    lmtc::Heap<int> hp(arrHeap,arrHeap+9);

095    hp.printHeap();

096    hp.increaseKey(3,5);

097    hp.insert(9);

098    hp.printHeap();

099    std::vector<int> st=hp.sort();

100    for(size_ti=0;i<st.size();i++)

101        std::cout<<st[i]<<"  ";

102    std::cout<<std::endl;

103    hp.printHeap();

104    std::cout<<std::endl;//*/

105 

106    ///*全域完全散列测试

107    std::cout<<"全域完全散列测试："<<std::endl;

108    lmtc::CompleteHash<unsignedlong> hash(50,lmtc::MAX_PRIME,keyToNumber);

109    hash.insert(20);

110    for(unsignedint i=0;i<1000;i++)

111        hash.insert((unsignedlong)(lmtc::averageRandom()*50000));

112    hash.completeHashOptimize();

113    unsignedlong *hashItem=hash.search(20);

114    hash.remove(*hashItem);

115    std::cout<<std::endl;//*/

116   

117 

118    ///*红黑树测试

119    std::cout<<"红黑树测试："<<std::endl;

120    lmtc::RedBlackTree<int> rbtree;

121    std::cout<<"随机插入"<<std::endl;

122    for(inti=1;i<20;i++)

123    {

124        inttemp=(int)(lmtc::averageRandom()*50000);//(int)(lmtc::averageRandom()*100000);

125        if(temp!=rbtree.getItemValue(rbtree.insert(temp)))

126            std::cout<<"error"<<std::endl;

127        if(temp!=rbtree.getItemValue(rbtree.search(temp)))

128            std::cout<<"error"<<std::endl;

129        std::cout<<temp<<std::endl;

130    }

131 

132    std::cout<<"size= "<<rbtree.size()<<std::endl;

133   

134    rbtree.insert(12);

135    lmtc::RedBlackTree<int>::ItemType &po=rbtree.deleteItem(rbtree.search(12));

136 

137    for(unsignedint i=0;i<rbtree.size();i++){

138        lmtc::RedBlackTree<int>::ItemType t=rbtree.searchItmOfOrder(i);

139        if(rbtree.getOrderOfItm(t)!=i)

140            std::cout<<"顺序统计错误"<<std::endl;

141    }

142 

143    std::cout<<"size= "<<rbtree.size()<<std::endl;

144 

145    std::cout<<"随机删除"<<std::endl;

146    for(inti=1;i<20;i++)

147    {

148        inttemp=(int)(lmtc::averageRandom()*50000);//(int)(lmtc::averageRandom()*100000);

149        lmtc::RedBlackTree<int>::ItemType del=rbtree.deleteItem(rbtree.search(temp));

150        if(!del.isEmpty())

151            std::cout<<rbtree.getItemValue(del)<<std::endl;

152    }

153 

154    std::cout<<"size= "<<rbtree.size()<<std::endl;

155 

156    std::cout<<"逆序输出："<<std::endl;

157    lmtc::RedBlackTree<int>::ItemType x=rbtree.maximum();

158    while(true){

159        if(x.isEmpty())

160            break;

161        std::cout<<rbtree.getItemValue(x)<<std::endl;

162        x=rbtree.predecessor(x);

163    }

164 

165    std::cout<<"中序遍历"<<std::endl;

166    rbtree.traver_inOrder(f);

167    std::cout<<"前序遍历"<<std::endl;

168    rbtree.traver_preOrder(f);

169    std::cout<<"遍历结束"<<std::endl;

170 

171    std::cout<<"开始验证"<<std::endl;

172    rbtree.asertTree(rbtree.getRoot());

173 

174    for(unsignedint i=0;i<rbtree.size();i++){

175        lmtc::RedBlackTree<int>::ItemType t=rbtree.searchItmOfOrder(i);

176        if(rbtree.getOrderOfItm(t)!=i)

177            std::cout<<"顺序统计错误"<<std::endl;

178    }

179    std::cout<<rbtree.size()<<std::endl;

180    rbtree.setEmpty();

181    std::cout<<"over"<<rbtree.size()<<std::endl;

182    std::cout<<std::endl;//*/

183 

184///*利用标准库顺序统计测试

185std::cout<<"标准库顺序统计测试："<<std::endl;

186int ntha1[]={2,5,8,10,12,7,99,3,54,46};

187int ntha2[]={3,5,6,13,9,19,8,7,54,28,99,46,3,5,2,5,8,10,12,7,99,3,54,46,2,5,8,10,12,7,99,3,54,46,2,5,8,10,12,7,99,3,54,46};

188std::vector<int> nthv1(ntha1,ntha1+10);

189std::vector<int> nthv2(ntha2,ntha2+44);

190std::vector<int> nthvec=lmtc::longestCommonSubsequence<int>(nthv1.begin(),nthv1.end(),nthv2.begin(),nthv2.end());

191for(unsignedint i=0;i<nthvec.size();i++)

192    std::cout<<nthvec[i]<<std::endl;

193std::random_shuffle(ntha2,ntha2+44);

194std::nth_element(ntha2,ntha2+20,ntha2+44);

195std::cout<<"顺序统计"<<std::endl;

196for(inti=0;i<44;i++)

197    std::cout<<ntha2[i]<<std::endl;

198std::cout<<std::endl;//*/

199 

200///*二项堆测试

201std::cout<<"二项堆测试："<<std::endl;

202lmtc::BinomialHeap<int> bihp;

203lmtc::BinomialHeap<int> bihp1;

204 

205for(inti=0;i<100;i++){

206    inttemp=(int)(lmtc::averageRandom()*5000000);

207    inttemp1=(int)(lmtc::averageRandom()*5000000);

208    bihp.increaseKey(bihp.insert(temp),(int)(lmtc::averageRandom()*5000000));

209    bihp1.increaseKey(bihp1.insert(temp1),(int)(lmtc::averageRandom()*5000000));

210    if(bihp.get_prio().isEmpty()||bihp1.get_prio().isEmpty())

211        std::cout<<"error"<<std::endl;

212    if(lmtc::averageRandom()<0.5){

213        intpio=bihp.getItemValue(bihp.get_prio());

214        intpio1=bihp.getItemValue(bihp.extract_prio());

215        if(pio!=pio1)

216            std::cout<<"error1"<<std::endl;

217    }

218    if(lmtc::averageRandom()<0.5){

219        intpio=bihp1.getItemValue(bihp1.get_prio());

220        intpio1=bihp1.getItemValue(bihp1.extract_prio());

221        if(pio!=pio1)

222            std::cout<<"error1"<<std::endl;

223    }

224 

225    if(lmtc::averageRandom()<0.1)

226        bihp.asertBinomialHeap();

227 

228    if(lmtc::averageRandom()<0.1)

229        bihp1.asertBinomialHeap();

230}

231 

232bihp.heapUnion(bihp1);

233std::cout<<"验证开始"<<std::endl;

234bihp.asertBinomialHeap();

235std::cout<<"验证结束"<<std::endl;

236 

237std::cout<<bihp.size()<<"end"<<bihp1.size()<<std::endl;

238 

239std::cout<<"释放资源1"<<std::endl;

240bihp.setEmpty();

241std::cout<<"释放资源2"<<std::endl;

242bihp.traver_inOrder(f);

243std::cout<<std::endl;//*/

244 

245///*斐波那契堆测试

246std::cout<<"斐波那契堆测试："<<std::endl;

247lmtc::FibonacciHeap<int> fbnqhp;

248lmtc::FibonacciHeap<int> fbnqhp1;

249for(inti=0;i<100;i++){

250    inttemp=(int)(lmtc::averageRandom()*5000000);

251    inttemp1=(int)(lmtc::averageRandom()*5000000);

252    fbnqhp.increaseKey(fbnqhp.insert(temp),(int)(lmtc::averageRandom()*5000000));

253    fbnqhp1.increaseKey(fbnqhp1.insert(temp1),(int)(lmtc::averageRandom()*5000000));

254    if(fbnqhp.get_prio().isEmpty()||fbnqhp1.get_prio().isEmpty())

255        std::cout<<"error"<<std::endl;

256    if(lmtc::averageRandom()<0.5){

257        intpio=fbnqhp.getItemValue(fbnqhp.get_prio());

258        intpio1=fbnqhp.getItemValue(fbnqhp.extract_prio());

259        if(pio!=pio1)

260            std::cout<<"error1"<<std::endl;

261    }

262    if(lmtc::averageRandom()<0.5){

263        intpio=fbnqhp1.getItemValue(fbnqhp1.get_prio());

264        intpio1=fbnqhp1.getItemValue(fbnqhp1.extract_prio());

265        if(pio!=pio1)

266            std::cout<<"error1"<<std::endl;

267    }

268 

269    if(lmtc::averageRandom()<0.1)

270        fbnqhp.asertFibonacciHeap();

271 

272    if(lmtc::averageRandom()<0.1)

273        fbnqhp1.asertFibonacciHeap();

274}

275fbnqhp.heapUnion(fbnqhp1);

276 

277std::cout<<"验证开始"<<std::endl;

278fbnqhp.asertFibonacciHeap();

279std::cout<<"验证结束"<<std::endl;

280 

281std::cout<<fbnqhp.size()<<"end"<<fbnqhp1.size()<<std::endl;

282 

283std::cout<<"释放资源1"<<std::endl;

284fbnqhp.setEmpty();

285std::cout<<"释放资源2"<<std::endl;

286std::cout<<std::endl;//*/

287 

288///*不相交集合测试

289std::cout<<"不相交集合测试："<<std::endl;

290std::vector<lmtc::NonIntersectSet<int>> setVec;

291for(inti=0;i<10;i++){

292    setVec.push_back(lmtc::NonIntersectSet<int>());

293    setVec[i].setValue(i);

294    if(setVec[i].getValue()!=i)

295        std::cout<<"error1"<<std::endl;

296}

297for(inti=0;i<10;i++){

298    unsignedint a=(unsignedint)(lmtc::averageRandom()*10);

299    unsignedint b=(unsignedint)(lmtc::averageRandom()*10);

300    if(setVec[a].find()==setVec[b].find())

301        std::cout<<a<<"===="<<b<<std::endl;

302    setVec[a].unionSet(setVec[b]);

303}

304for(inti=0;i<10;i++){

305    std::cout<<i<<"->"<<setVec[i].find()->getValue()<<std::endl;

306    setVec[i].find()->setValue(setVec[i].find()->getValue());

307}

308std::cout<<std::endl;//*/

309 

310int arr[][5]={{0,1,0,0,0},{0,0,0,0,0},{1,1,0,1,1},{1,0,0,0,0},{1,0,0,1,0}};

311lmtc::Array<int> ajacencyMatrix(2,5,5);

312for(inti=0;i<5;i++)

313    for(intj=0;j<5;j++)

314        ajacencyMatrix(i,j)=arr[i][j];

315 

316//邻接矩阵转邻接表

317std::vector<std::list<lmtc::Edge<int>>> ajacencyList;

318lmtc::Graph::ajacencyMatrixToList(ajacencyMatrix,ajacencyList);

319 

320//邻接表的转置

321std::vector<std::list<lmtc::Edge<int>>> transposedAjacencyList;

322lmtc::Graph::transposeAjacencyList(ajacencyList,transposedAjacencyList);

323 

324 

325///*拓扑排序测试

326std::cout<<"拓扑排序测试："<<std::endl;

327std::vector<unsigned int> order=lmtc::Graph::topologicalSort(transposedAjacencyList);

328std::cout<<"拓扑序列为："<<std::endl;

329for(std::vector<unsignedint>::iterator iter=order.begin();iter!=order.end();iter++)

330    std::cout<<*iter<<"  ";

331std::cout<<std::endl;//*/

332 

333///*强连通分支测试

334std::cout<<"强连通分支测试："<<std::endl;

335std::vector<std::vector<unsigned int>> strongConnectComponents;

336lmtc::Graph::computeStrngConctComps(ajacencyList,strongConnectComponents);

337for(unsignedint i=0;i<strongConnectComponents.size();i++){

338    std::cout<<"连通分支"<<i<<std::endl;

339    for(std::vector<unsignedint>::iterator iter=strongConnectComponents[i].begin();iter!=strongConnectComponents[i].end();iter++)

340        std::cout<<*iter<<" ; ";

341    std::cout<<std::endl;

342}

343std::cout<<std::endl;//*/

344 

345///*欧拉回路测试

346std::cout<<"欧拉回路测试："<<std::endl;

347int arr1[][5]={{1,1,1,0,0},{0,1,1,1,0},{0,0,0,1,1},{1,0,0,0,1},{1,1,0,0,0}};

348//int arr1[][5]={{2,0,1,1,0},{0,2,0,1,1},{1,0,0,0,1},{1,1,0,0,0},{0,1,1,0,0}};

349lmtc::Array<int> ajacencyMatrix1(2,5,5);

350for(inti=0;i<5;i++)

351    for(intj=0;j<5;j++)

352        ajacencyMatrix1(i,j)=arr1[i][j];

353std::vector<std::list<lmtc::Edge<int>>> ajacencyList1;

354lmtc::Graph::ajacencyMatrixToList(ajacencyMatrix1,ajacencyList1);

355 

356std::list<lmtc::Edge<int>> EulerCircuit;

357bool existEuler=lmtc::Graph::computeEulerCircuit(ajacencyList1,EulerCircuit,true);

358if(existEuler==false&&EulerCircuit.empty())

359    std::cout<<"Euler circuit or path not existed"<<std::endl;

360else{

361    std::cout<<"欧拉回路为："<<std::endl;

362    for(std::list<lmtc::Edge<int>>::iterator iter=EulerCircuit.begin();iter!=EulerCircuit.end();iter++)

363        std::cout<<iter->vSt<<"-"<<iter->vEd<<"   -->   ";

364    std::cout<<std::endl;

365}

366std::cout<<std::endl;//*/

367 

368///*最小生成树测试

369std::cout<<"最小生成树测试："<<std::endl;

370int arr2[][8]={{1,1,2,7,3,0,0,0},{1,1,6,2,0,0,0,0},{2,6,1,4,0,0,0,0},{7,2,4,1,5,0,0,0},{3,0,0,5,1,0,0,0},{0,0,0,0,0,1,2,3},{0,0,0,0,0,2,1,1},{0,0,0,0,0,3,1,1}};

371lmtc::Array<int> ajacencyMatrix2(2,8,8);

372for(inti=0;i<8;i++)

373    for(intj=0;j<8;j++)

374        ajacencyMatrix2(i,j)=arr2[i][j];

375std::vector<std::list<lmtc::Edge<int>>> ajacencyList2;

376lmtc::Graph::ajacencyMatrixToList(ajacencyMatrix2,ajacencyList2);

377std::vector<std::list<lmtc::Edge<int>>> mstAjacencyList;

378int weightKruskal=lmtc::Graph::mstKruskal(ajacencyList2,mstAjacencyList);

379std::cout<<"Kruskal最小生成树为：权重："<<weightKruskal<<std::endl;

380for(unsignedint i=0;i<mstAjacencyList.size();i++){

381    for(std::list<lmtc::Edge<int>>::const_iterator iter=mstAjacencyList[i].begin();iter!=mstAjacencyList[i].end();iter++){

382        std::cout<<iter->vSt<<"-"<<iter->vEd<<":"<<iter->data<<"   ";

383    }

384    std::cout<<std::endl;

385}

386std::vector<std::list<lmtc::Edge<int>>> mstAjacencyListPrim;

387int weightPrim=lmtc::Graph::mstPrim(ajacencyList2,mstAjacencyListPrim);

388std::cout<<"Prim最小生成树为：权重："<<weightPrim<<std::endl;

389for(unsignedint i=0;i<mstAjacencyListPrim.size();i++){

390    for(std::list<lmtc::Edge<int>>::const_iterator iter=mstAjacencyListPrim[i].begin();iter!=mstAjacencyListPrim[i].end();iter++){

391        std::cout<<iter->vSt<<"-"<<iter->vEd<<":"<<iter->data<<"   ";

392    }

393    std::cout<<std::endl;

394}

395std::cout<<std::endl;//*/

396 

397///*判断有向图或者无向图是否存在回路

398std::cout<<"判断有向图或者无向图是否存在回路测试："<<std::endl;

399int arr3[][4]={{1,1,1,1},{0,0,0,0},{0,0,0,0},{0,0,0,0}};

400lmtc::Array<int> ajacencyMatrix3(2,4,4);

401for(inti=0;i<4;i++)

402    for(intj=0;j<4;j++)

403        ajacencyMatrix3(i,j)=arr3[i][j];

404std::vector<std::list<lmtc::Edge<int>>> ajacencyList3;

405lmtc::Graph::ajacencyMatrixToList(ajacencyMatrix3,ajacencyList3);

406if(lmtc::Graph::hasLoop(ajacencyList3,true)==true)

407    std::cout<<"存在回路"<<std::endl;

408else

409    std::cout<<"不存在回路"<<std::endl;

410std::cout<<std::endl;//*/

411 

412///*计算最短路径

413std::cout<<"最短路径测试："<<std::endl;

414int arr4[][4]={{0,1,3,4},{0,3,1,2},{2,1,6,0},{6,3,0,1}};

415lmtc::Array<int> ajacencyMatrix4(2,4,4);

416for(inti=0;i<4;i++)

417    for(intj=0;j<4;j++)

418        ajacencyMatrix4(i,j)=arr4[i][j];

419std::vector<std::list<lmtc::Edge<int>>> ajacencyList4;

420lmtc::Graph::ajacencyMatrixToList(ajacencyMatrix4,ajacencyList4);

421unsigned s=0;

422std::vector<unsigned int> p;

423std::vector<int>d;

424//if(true==lmtc::Graph::shortestPathBellmanFord(ajacencyList4,s,p,d))shortestPathDijkstra

425//if(true==lmtc::Graph::shortestPathOfDag(ajacencyList4,s,p,d))

426if(true==lmtc::Graph::shortestPathDijkstra(ajacencyList4,s,p,d))

427    for(unsignedint i=0;i<p.size();i++)

428        std::cout<<"parent of "<<i<<" is "<<p[i]<<" , and d of the path is "<<d[i]<<std::endl;

429else

430    std::cout<<"存在负边"<<std::endl;

431std::cout<<std::endl;//*/

432 

433///*计算每对顶点间的最短路径

434std::cout<<"每对顶点间的最短路径测试："<<std::endl;

435int arr5[][4]={{1,1,3,4},{0,1,1,2},{-2,-1,1,0},{6,3,0,1}};

436lmtc::Array<int> ajacencyMatrix5(2,4,4);

437for(inti=0;i<4;i++)

438    for(intj=0;j<4;j++)

439        ajacencyMatrix5(i,j)=arr5[i][j];

440std::vector<std::list<lmtc::Edge<int>>> ajacencyList5;

441lmtc::Graph::ajacencyMatrixToList(ajacencyMatrix5,ajacencyList5);

442lmtc::Array<unsigned int> pMatrix;

443lmtc::Array<int>dMatrix;

444//lmtc::Graph::shortestPathAllFloydWarshall(ajacencyList5,pMatrix,dMatrix);

445if(lmtc::Graph::shortestPathAllJohnson(ajacencyList5,pMatrix,dMatrix))

446for(unsignedint k=0;k<pMatrix.getDimLen(0);k++){

447    std::cout<<"以k="<<k<<"为源点的最短路径树"<<std::endl;

448    for(unsignedint i=0;i<pMatrix.getDimLen(1);i++)

449        std::cout<<"parent of "<<i<<" is "<<pMatrix(k,i)<<" , and d of the path is "<<dMatrix(k,i)<<std::endl;

450}

451std::cout<<std::endl;//*/

452 

453///*计算最大流

454std::cout<<"最大流测试："<<std::endl;

455int arr6[][6]={{-6,16,13,0,0,0},{0,0,10,12,0,0},{-6,4,0,0,14,0},{0,0,9,0,0,20},{0,0,0,7,0,4},{-4,0,0,-4,0,-4}};

456lmtc::Array<int> ajacencyMatrix6(2,6,6);

457for(inti=0;i<6;i++)

458    for(intj=0;j<6;j++)

459        ajacencyMatrix6(i,j)=arr6[i][j];

460std::vector<std::list<lmtc::Edge<int>>> ajacencyList6;

461lmtc::Graph::ajacencyMatrixToList(ajacencyMatrix6,ajacencyList6);

462std::vector<std::list<lmtc::Edge<int>>> flow;

463//int maxF=lmtc::Graph::maximumFlowFordFulkerson_EdmondsKarp(ajacencyList6,0,5,flow);

464int maxF=lmtc::Graph::maximumFlowPushRelabelToFront(ajacencyList6,0,5,flow);

465std::cout<<"最大流="<<maxF<<std::endl;

466for(unsignedint i=0;i<flow.size();i++)

467    for(std::list<lmtc::Edge<int>>::const_iterator iter=flow[i].begin();iter!=flow[i].end();iter++)

468        std::cout<<iter->vSt<<"-"<<iter->vEd<<":"<<iter->data<<"   ";

469std::cout<<std::endl;

470std::cout<<std::endl;//*/

471 

472 

473/*矩阵运算*/

474 

475///*矩阵乘法测试

476std::cout<<"矩阵乘法测试："<<std::endl;

477int arr7[][3]={{0,1,2},{3,1,5},{4,2,6},{1,8,3}};

478lmtc::Array<int> Matrix1(2,4,3);

479for(inti=0;i<4;i++)

480    for(intj=0;j<3;j++)

481        Matrix1(i,j)=arr7[i][j];

482 

483int arr8[][5]={{4,8,2,1,3},{3,2,6,4,7},{2,5,6,3,5}};

484lmtc::Array<int> Matrix2(2,3,5);

485for(inti=0;i<3;i++)

486    for(intj=0;j<5;j++)

487        Matrix2(i,j)=arr8[i][j];

488 

489//lmtc::Array<int> matrixRs=lmtc::Algebra::matrixMultiplySimple(Matrix1,Matrix2);

490lmtc::Array<int> matrixRs=lmtc::Algebra::matrixMultiplyStrassen(Matrix1,Matrix2);

491std::cout<<"矩阵乘法结果"<<std::endl;

492for(unsignedint i=0;i<matrixRs.getDimLen(0);i++){

493    for(unsignedint j=0;j<matrixRs.getDimLen(1);j++)

494        std::cout<<matrixRs(i,j)<<"   ";

495    std::cout<<std::endl;

496}

497std::cout<<std::endl;//*/

498 

499///*LUP 分解测试

500std::cout<<"LUP 分解测试："<<std::endl;

501double arr9[][4]={{2,0,2,0.6},{3,3,4,-2},{5,5,4,2},{-1,-2,3.4,-1}};

502lmtc::Array<double> A(2,4,4);

503for(inti=0;i<4;i++)

504    for(intj=0;j<4;j++)

505        A(i,j)=arr9[i][j];

506lmtc::Array<double> L,U;

507std::vector<unsigned int> P;

508bool b=lmtc::Algebra::lupDecompose(A,L,U,P);

509std::cout<<"LUP分解结果"<<std::endl;

510if(b==true){

511    std::cout<<"L:"<<std::endl;

512    for(unsignedint i=0;i<L.getDimLen(0);i++){

513        for(unsignedint j=0;j<L.getDimLen(1);j++)

514            std::cout<<L(i,j)<<"   ";

515        std::cout<<std::endl;

516    }

517    std::cout<<"U:"<<std::endl;

518    for(unsignedint i=0;i<U.getDimLen(0);i++){

519        for(unsignedint j=0;j<U.getDimLen(1);j++)

520            std::cout<<U(i,j)<<"   ";

521        std::cout<<std::endl;

522    }

523    std::cout<<"P:"<<std::endl;

524    for(unsignedint i=0;i<P.size();i++)

525        std::cout<<P[i]<<"   ";

526    std::cout<<std::endl;

527}

528else

529    std::cout<<"LUP分解失败，奇异方阵！"<<std::endl;

530std::cout<<std::endl;//*/

531 

532///*解线性方程组

533std::cout<<"解线性方程组测试："<<std::endl;

534double arr10[][3]={{1,2,0},{3,4,4},{5,6,3}};

535std::vector<double> b1;

536b1.push_back(3);

537b1.push_back(7);

538b1.push_back(8);

539lmtc::Array<double> A1(2,3,3);

540for(inti=0;i<3;i++)

541    for(intj=0;j<3;j++)

542        A1(i,j)=arr10[i][j];

543std::vector<double> X1;

544if(lmtc::Algebra::solveLinearEquationsFast(A1,b1,X1))

545//if(lmtc::Algebra::solveLinearEquationsByLUP(A1,b1,X1))

546    for(unsignedint i=0;i<X1.size();i++)

547        std::cout<<X1[i]<<"   ";

548else

549    std::cout<<"无解";

550std::cout<<std::endl;//*/

551 

552///*求逆矩阵测试

553std::cout<<"矩阵求逆测试："<<std::endl;

554double arr11[][3]={{1,2,0},{2,4,5},{5,6,3}};

555lmtc::Array<double> A2(2,3,3);

556for(inti=0;i<3;i++)

557    for(intj=0;j<3;j++)

558        A2(i,j)=arr11[i][j];

559lmtc::Array<double> _A;

560std::cout<<"矩阵求逆"<<std::endl;

561if(lmtc::Algebra::inverseMatrixFast(A2,_A)){

562//if(lmtc::Algebra::inverseMatrixByLUP(A2,_A)){

563    for(unsignedint i=0;i<_A.getDimLen(0);i++){

564        for(unsignedint j=0;j<_A.getDimLen(1);j++)

565            std::cout<<_A(i,j)<<"   ";

566        std::cout<<std::endl;

567    }

568    lmtc::Array<double> matrixRs=lmtc::Algebra::matrixMultiplyStrassen(A2,_A);

569    std::cout<<"与逆矩阵相乘的单位矩阵"<<std::endl;

570    for(unsignedint i=0;i<matrixRs.getDimLen(0);i++){

571        for(unsignedint j=0;j<matrixRs.getDimLen(1);j++)

572            std::cout<<matrixRs(i,j)<<"   ";

573        std::cout<<std::endl;

574    }

575}

576else

577    std::cout<<"不可逆"<<std::endl;

578std::cout<<std::endl;//*/

579 

580///*伪逆矩阵测试

581std::cout<<"伪逆矩阵测试："<<std::endl;

582double arr12[][3]={{1,-1,1},{1,1,1},{1,2,4},{1,3,9},{1,5,25}};

583lmtc::Array<double> A3(2,5,3);

584for(inti=0;i<5;i++)

585    for(intj=0;j<3;j++)

586        A3(i,j)=arr12[i][j];

587lmtc::Array<double> A3T=lmtc::Algebra::transposeMatrix(A3);

588lmtc::Array<double> _Af;

589lmtc::Array<double> _Ab;

590std::cout<<"求前向伪逆矩阵"<<std::endl;

591if(lmtc::Algebra::pseudoInverseMatrixForward(A3,_Af)){

592    for(unsignedint i=0;i<_Af.getDimLen(0);i++){

593        for(unsignedint j=0;j<_Af.getDimLen(1);j++)

594            std::cout<<_Af(i,j)<<"   ";

595        std::cout<<std::endl;

596    }

597    lmtc::Array<double> matrixRs=lmtc::Algebra::matrixMultiplyStrassen(_Af,A3);

598    std::cout<<"与前向伪逆矩阵相乘的单位矩阵"<<std::endl;

599    for(unsignedint i=0;i<matrixRs.getDimLen(0);i++){

600        for(unsignedint j=0;j<matrixRs.getDimLen(1);j++)

601            std::cout<<matrixRs(i,j)<<"   ";

602        std::cout<<std::endl;

603    }

604}

605else

606    std::cout<<"不可前向伪逆"<<std::endl;

607 

608std::cout<<"求后向伪逆矩阵"<<std::endl;

609if(lmtc::Algebra::pseudoInverseMatrixBackward(A3T,_Ab)){

610    for(unsignedint i=0;i<_Ab.getDimLen(0);i++){

611        for(unsignedint j=0;j<_Ab.getDimLen(1);j++)

612            std::cout<<_Ab(i,j)<<"   ";

613        std::cout<<std::endl;

614    }

615    lmtc::Array<double> matrixRs=lmtc::Algebra::matrixMultiplyStrassen(A3T,_Ab);

616    std::cout<<"与后向伪逆矩阵相乘的单位矩阵"<<std::endl;

617    for(unsignedint i=0;i<matrixRs.getDimLen(0);i++){

618        for(unsignedint j=0;j<matrixRs.getDimLen(1);j++)

619            std::cout<<matrixRs(i,j)<<"   ";

620        std::cout<<std::endl;

621    }

622}

623else

624    std::cout<<"不可后向伪逆"<<std::endl;

625std::cout<<std::endl;//*/

626 

627///*最小二乘估计与解正态方程组测试

628std::cout<<"解正态方程组测试："<<std::endl;

629double arr13[][3]={{1,-1,1},{1,1,1},{1,2,4},{1,3,9},{1,5,25}};

630lmtc::Array<double> A4(2,5,3);

631for(inti=0;i<5;i++)

632    for(intj=0;j<3;j++)

633        A4(i,j)=arr13[i][j];

634double yArr[]={2,1,1,0,3};

635double xArr[]={-1,1,2,3,5};

636std::vector<double> X(xArr,xArr+5);

637std::vector<double> Y(yArr,yArr+5);

638std::vector<double> C;

639std::vector<double> C1;

640if(lmtc::Algebra::solveNormalityEquationsByLUP(A4,Y,C)){

641//if(lmtc::Algebra::solveNormalityEquationsFast(A4,Y,C)){

642    std::cout<<"正态方程组的解为：";

643    for(unsignedint i=0;i<C.size();i++)

644        std::cout<<C[i]<<"   ";

645    std::cout<<std::endl;

646}

647else

648    std::cout<<"正态方程组不可解"<<std::endl;

649std::cout<<"最小二乘估计测试："<<std::endl;

650//if(lmtc::Algebra::leastSquaresEstimationByLUP(polynomial,3,X,Y,C1)){

651if(lmtc::Algebra::leastSquaresEstimationFast(polynomial,3,X,Y,C1)){

652    std::cout<<"最小二乘估计的解为：";

653    for(unsignedint i=0;i<C1.size();i++)

654        std::cout<<C1[i]<<"   ";

655    std::cout<<std::endl;

656}

657else

658    std::cout<<"最小二乘不可解"<<std::endl;

659std::cout<<std::endl;//*/

660 

661///*方幂和测试

662std::cout<<"方幂和测试："<<std::endl;

663const unsigned int K=6;

664const unsigned int N=100;

665lmtc::Array<double> S=lmtc::Algebra::powerSumFormula(K);

666std::cout<<"方幂和公式如下"<<std::endl;

667for(unsignedint i=0;i<=K;i++){

668    std::cout<<"幂为"<<i<<"时的公式：";

669    for(unsignedint j=0;j<i+2;j++)

670        std::cout<<"c"<<j<<" = "<<S(i,j)<<", ";

671    std::cout<<std::endl;

672}

673std::cout<<"when K="<<K<<" ,N="<<N<<"  时，0^K+1^K+2^K+...+N^K = "<<lmtc::Algebra::powerSum(N,K)<<std::endl;

674std::cout<<std::endl;//*/

675 

676///*复数运算测试

677std::cout<<"复数运算测试："<<std::endl;

678lmtc::ComplexNumber c1(0.34,0.69);

679lmtc::ComplexNumber c2(0.53,0.49);

680c1+=c2;

681std::cout<<(std::string)c1<<std::endl;

682c1-=c2;

683std::cout<<(std::string)c1<<std::endl;

684c1*=c2;

685std::cout<<(std::string)c1<<std::endl;

686c1/=c2;

687std::cout<<(std::string)c1<<std::endl;

688std::cout<<std::endl;//*/

689 

690///*快速傅里叶变换,及逆变换，数值稳定性非常好，当cN=1000时精度为0.00000000001，当cN=10000时为0.0000000001，当cN=100000时为0.000000001，当cN=1000000时为0.00000001，

691std::cout<<"快速傅里叶变换测试："<<std::endl;

692std::vector<lmtc::ComplexNumber> a;

693unsigned intcN=100;

694for(unsigned i=0;i<cN;i++)

695    a.push_back(lmtc::ComplexNumber(lmtc::averageRandom(0,10),lmtc::averageRandom(0,10)));

696std::vector<lmtc::ComplexNumber> y=lmtc::FastFourierTransform::fft(a);

697std::vector<lmtc::ComplexNumber> a1=lmtc::FastFourierTransform::_fft(y);

698std::cout<<"系数向量a：";

699for(unsigned i=0;i<a.size();i++)

700    std::cout<<(std::string)a[i]<<"   ";

701std::cout<<std::endl;

702std::cout<<"点值向量y：";

703for(unsigned i=0;i<y.size();i++)

704    std::cout<<(std::string)y[i]<<"   ";

705std::cout<<std::endl;

706std::cout<<"系数向量a1：";

707for(unsigned i=0;i<a1.size();i++)

708    std::cout<<(std::string)a1[i]<<"   ";

709std::cout<<std::endl;

710 

711for(unsignedint i=0;i<a1.size();i++){

712    if(i<a.size()&&!a[i].equal(a1[i],0.00001))

713        std::cout<<"error:"<<(std::string)a[i]<<"==?"<<(std::string)a1[i]<<"   "<<std::endl;

714    if(i>=a.size()&&!a1[i].equal(0,0.01))

715        std::cout<<"error1"<<std::endl;

716}

717std::cout<<std::endl;//*/

718 

719///*求卷积

720std::cout<<"卷积测试："<<std::endl;

721double Aarr[]={1,2,3,23};

722double Barr[]={5,4,34};

723std::vector<lmtc::ComplexNumber> Avec(Aarr,Aarr+4);

724std::vector<lmtc::ComplexNumber> Bvec(Barr,Barr+3);

725std::vector<lmtc::ComplexNumber>Cvec=lmtc::Algebra::convolution(Avec,Bvec);

726std::cout<<"卷积为："<<std::endl;

727for(unsignedint i=0;i<Cvec.size();i++)

728    std::cout<<(std::string)Cvec[i]<<"   ";

729std::cout<<std::endl;//*/

730 

731///*多维傅里叶变换测试

732std::cout<<"高维傅里叶变换测试："<<std::endl;

733lmtc::Array<lmtc::ComplexNumber> aMultiFFT(2,3,4);

734for(unsignedint i=0;i<aMultiFFT.getDimLen(0);i++)

735    for(unsignedint j=0;j<aMultiFFT.getDimLen(1);j++)

736        aMultiFFT(i,j)=lmtc::ComplexNumber(lmtc::averageRandom(0,10),lmtc::averageRandom(0,10));

737lmtc::Array<lmtc::ComplexNumber> yMultiFFT=lmtc::FastFourierTransform::fft(aMultiFFT);

738lmtc::Array<lmtc::ComplexNumber> a1MultiFFT=lmtc::FastFourierTransform::_fft(yMultiFFT);

739std::cout<<"高维傅里叶变换"<<std::endl;

740std::cout<<"aMultiFFT:"<<std::endl;

741for(unsignedint i=0;i<aMultiFFT.getDimLen(0);i++){

742    for(unsignedint j=0;j<aMultiFFT.getDimLen(1);j++)

743        std::cout<<(std::string)aMultiFFT(i,j)<<", ";

744    std::cout<<std::endl;

745}

746std::cout<<"yMultiFFT:"<<std::endl;

747for(unsignedint i=0;i<yMultiFFT.getDimLen(0);i++){

748    for(unsignedint j=0;j<yMultiFFT.getDimLen(1);j++)

749        std::cout<<(std::string)yMultiFFT(i,j)<<", ";

750    std::cout<<std::endl;

751}

752std::cout<<"a1MultiFFT:"<<std::endl;

753for(unsignedint i=0;i<a1MultiFFT.getDimLen(0);i++){

754    for(unsignedint j=0;j<a1MultiFFT.getDimLen(1);j++)

755        std::cout<<(std::string)a1MultiFFT(i,j)<<", ";

756    std::cout<<std::endl;

757}

758for(unsignedint i=0;i<a1MultiFFT.getDimLen(0);i++)

759    for(unsignedint j=0;j<a1MultiFFT.getDimLen(1);j++){

760        if(i<aMultiFFT.getDimLen(0)&&j<aMultiFFT.getDimLen(1)){

761            if(!aMultiFFT(i,j).equal(a1MultiFFT(i,j),0.00001))

762                std::cout<<"error:"<<(std::string)a1MultiFFT(i,j)<<"==?"<<(std::string)aMultiFFT(i,j)<<"   "<<std::endl;

763        }else{

764            if(!a1MultiFFT(i,j).equal(0,0.00001))

765                std::cout<<"error1:"<<(std::string)a1MultiFFT(i,j)<<"==?"<<0<<"   "<<std::endl;

766        }

767    }

768std::cout<<std::endl;//*/

769 

770///*多维卷积测试

771std::cout<<"多维卷积测试："<<std::endl;

772lmtc::Array<lmtc::ComplexNumber> mutiAvec(2,3,2);

773lmtc::Array<lmtc::ComplexNumber> mutiBvec(3,2,3,3);

774mutiAvec(1,1)=1;

775mutiAvec(2,1)=2;

776mutiAvec(0,1)=4;

777mutiBvec(1,1,1)=3;

778mutiBvec(1,2,0)=6;

779mutiBvec(0,1,2)=5;

780lmtc::Array<lmtc::ComplexNumber>mutiCvec=lmtc::Algebra::convolution(mutiAvec,mutiBvec);

781std::cout<<"多维卷积为："<<std::endl;

782std::vector<unsigned int> dimC(mutiCvec.getDimNum(),0);

783while(dimC[mutiCvec.getDimNum()-1]<mutiCvec.getDimLen(mutiCvec.getDimNum()-1)){//将a中元素复制到正合幂数组a1

784    if(!mutiCvec(dimC).equal(0,0.001)){

785        std::cout<<"mutiCvec(";

786        for(unsignedint i=0;i<dimC.size();i++)

787            std::cout<<dimC[i]<<"  ";

788        std::cout<<")   =   "<<(std::string)mutiCvec(dimC)<<std::endl;

789    }

790    for(unsignedint i=0;i<mutiCvec.getDimNum();i++){//计算多维数组遍历的下一个坐标

791        if(dimC[i]!=(mutiCvec.getDimLen(i)-1)){

792            dimC[i]++;

793            break;

794        }elseif(i!=(mutiCvec.getDimNum()-1))

795            dimC[i]=0;

796        else

797            dimC[i]++;

798    }

799}

800std::cout<<std::endl;//*/

801 

802 

803std::cout<<"all test over!!!"<<std::endl;

804    return0;//表示成功执行

805}