# Andriy Zhugayevych azh@ukr.net # BasicTools package # created 30.07.2005, modified - see version below # https://zhugayevych.me/maple/BasicTools/ # ################################################################################ #root: BasicTools #hfl: BasicTools #toc: _Overview #BasicTools,BasicTools/Setup BasicTools:=module() option package; global tex_pm, tex_cdot, tex_degree, tex_ll, tex_gg, tex_lnot, tex_mu, `index/symmetric12`, JP, `simplify/esu`,`simplify/emu`; export ModuleLoad, Setup, Round,mod1,Reduce2P,PeriodicMean,PeriodicStandardDeviation,angle2pair, ReduceFloat,ReduceFloat2,ReduceFloat3,ReduceFloatGrid, #Numerical conversions simplifyradical,coeffs2,idExch,idSym,idASym, #Expression tools SearchPos,Substitute,IdentifyPairs,Classify2,BlockSplit, #List tools Grid,Grid2,GetPath,convert2range,convert2repeat, #Grids and sequences PrintTree,Tree2Seq,Tree2Seq2,Tree2List, #Trees Dim2,idL,iidL,idI,iidI,idS,iidS,idA,iidA,idSS,iidSS,nextpointS, #Indexing isleq,Sort,Rort,SortM,SortIdx,MaxVal,MinVal,MaxIdx,MinIdx,SortEE,SortSign, #Sorting Squeeze, TrimLeft, TrimRight, Trim, StartNewLine, SearchKeySeq, ReadBody, ReadByMask, Split2, FormattedOperation, FormatFloat, AlignPoints, FormatTable, FormatTime2, PrintColumns, PrintMatrix, PrintUpperMatrix, PrintVector3d, ascii2unicode, # String tools CyclesAndTrees, defistree, CycleBasis2, # Combinatorics and algorithms Circumsphere, # Geometry fullparfrac2, fZero, fDblZero, fZeros, FZERO, NLPSolve2, HessianByGradient, PolynomialInterpolation1der, BezierControlPoints, BezierNodes, ExtrapolateMPE, # Calculus devpsolve, # Differential equations CheckEigenvector, NormalizeVector, NormalizeMatrix, RealJordanBlockMatrix, MatrixFunction2, SparseMM, SparseMHM, SparseMV, SparseMO, SparseMN, SparseMap, SuperimposeM, LocalizeVectors4, # Linear algebra RotationM, ReflectionM, RotationParam, RotationParamEulerX, RotationParamEulerY, RotationMEulerX, RotationMEulerY, EulerY2Rot, Rot2EulerY, RotationMp2q, RotationM2zx, TransSpherCoo, TransSpherCooI, TransTensor, ExtendedTransM, # Linear transformations Phim, NPhim2, Pnma, NPnma2, Ylm, Yxyz, YxyzM, YxyzR, YxyzL, Wignerd, WignerD, Wigner3j, Wigner6j, MatElemYlm, MatElemY, # Spherical harmonics LerchPhi2, polylogR, polylogR_init, # Other special functions ExpandPath, SimplifyPath, FoldersTree, LoadTextFile, UnloadTextFile, SearchFilePos, ReadLines, ReadValue, ReadPaginatedMatrix, WriteLines, ReplaceLines, ReadRecord, WriteRecord, ReadSet, ReadUTF8, WriteUTF8, ReadBIN, WriteBIN, PrintFile, convert2unix, convert2dos, OpenFolder, # File tools GetBox,GetBox2,Spline2,ErrorBar,Legend2,TextOnLine,Rainbow,WritePlotData,plotps,CorrectEPS,ImageCoreBox, # Graphics RootMeanSquare, LinearFit2, LinearFitTrig, MultiDimFit, GaussianFit, FindExtremum, ReorderXYdata, Intersection, # Data manipulation SqueezeHTML, SqueezeTable, ReadHref, ResolveURI, GetXMLChild, ReadXML, WriteXML, # HTML/XML tools GetConstantSU, ConvertSU, # Physics tools BuildHelp, ReadModules, AssignModule, ProcessSetupArgs, Check4Updates, TabulatedFunction, Timing, InitTiming; # System tools local version, MapleVersion, CyrCharMap, complexTrigBasis, conventionY, normalizedY, psplotoptions, phi,theta,psi,alpha,beta,gamma, polylogR_T, polylogR_n, polylogR_N, TextFiles,TextFilesAttr,TextFilesList,nTextFiles,sTextFiles, binformatidentifiers1, cc,hb,NA,kB,q2vu, s_TabulatedFunction, InitialTime,PrintTime,nDone,nLeft,nPerc,Percentages; ModuleLoad:=proc() local i,n,V; version:=20260215; # Check version MapleVersion:=sscanf(convert(kernelopts('version'),string),"Maple %f")[1]; if not(type(MapleVersion,numeric)) then error("Cannot determine Maple version: %1",kernelopts('version')) end; if (MapleVersion<10) then WARNING("Version 10 or later is needed for full operability. Detected version is %1",kernelopts('version')) end; # Special symbols tex_pm,tex_cdot,tex_degree,tex_ll,tex_gg,tex_lnot,tex_mu:=op(map(StringTools[Char],[177,183,176,171,187,172,181])); # CyrCharMap CyrCharMap:=table(): for i from 0 to 255 do CyrCharMap[i]:=cat("&#",i,";") end: for i from 192 to 255 do CyrCharMap[i]:=cat("&#",i+848,";") end: CyrCharMap[179]:="і": CyrCharMap[175]:="Ї": #³² CyrCharMap[191]:="ї": CyrCharMap[178]:="І": #¿¯ CyrCharMap[186]:="є": CyrCharMap[170]:="Є": #ºª CyrCharMap[184]:="ё": CyrCharMap[168]:="Ё": #¸¨ CyrCharMap[180]:="ґ": CyrCharMap[165]:="Ґ": #´¥ # Default settings Setup('InitTextFiles', 'complexTrigBasis'=true, 'conventionY'="standard", 'normalizedY'=true, 'nTextFiles'=1000, 'sTextFiles'=5*10^8, 'psplotoptions'="portrait,noborder,width=400,height=400"); # Indexing functions `index/symmetric12`:=proc(id::list(posint),A::Array,v::list) local i; if (nops(id)<2) then error("symmetric12 indexing requires at least 2 indices") end; i:=max(id[1],id[2]),min(id[1],id[2]),op(3..,id); if (nargs=2) then A[i] else A[i]:=op(v) end end; # Legendre and Jacobi polynomials JP:=orthopoly[P]; # LerchPhi patch polylogR_init(20,100); # Binary format identifiers binformatidentifiers1:=[1760568055,-148705176]; # System of units try Units[AddBaseUnit]('esu',context='SU',symbol='esu',dimension='electric_charge') catch: WARNING("Cannot add 'esu' unit because of the following error:"); WARNING(lasterror) end; try Units[AddBaseUnit]('emu',context='SU',symbol='emu',dimension='magnetic_flux') catch: WARNING("Cannot add 'emu' unit because of the following error:"); WARNING(lasterror) end; try Units[AddSystem]('SU',check=true,'cm'['SI'],'g','s','esu','emu','A','K','mol','cd') catch: WARNING("Cannot add 'SU' system because of the following error:"); WARNING(lasterror) end; Units[UseSystem]('SU',strict=true); cc:=ScientificConstants[GetValue](ScientificConstants[Constant]('c',system='SU')); hb:=ScientificConstants[GetValue](ScientificConstants[Constant]('hbar',system='SU')); NA:=ScientificConstants[GetValue](ScientificConstants[Constant]('N[A]',system='SU')); kB:=ScientificConstants[GetValue](ScientificConstants[Constant]('k',system='SU')); `simplify/esu`:=proc(expr) local u,v,e; v,u:=q2vu(simplify(expr)); e:=subs('A'=cc/10*'esu'/'s',parse(convert(op(u),string))); simplify(Units[Unit](v*e)) end; `simplify/emu`:=proc(expr) local u,v,e; v,u:=q2vu(simplify(expr)); e:=subs('A'=1/10*'emu'/'s',parse(convert(op(u),string))); simplify(Units[Unit](v*e)) end; NULL end: Setup:=proc({InitTextFiles::boolean:=false,printout::boolean:=false}) local sq; sq:=ProcessSetupArgs([_rest],['complexTrigBasis','conventionY','normalizedY','polylogR_n','polylogR_N','nTextFiles','sTextFiles','psplotoptions'],['version']); if printout then printf("complexTrigBasis=%a, conventionY=%a, normalizedY=%a\nnTextFiles=%d, sTextFiles=%.0fMB\npsplotoptions=%a\n", complexTrigBasis,conventionY,normalizedY,nTextFiles,sTextFiles*1e-6,psplotoptions) end; if InitTextFiles then TextFiles:=table(); TextFilesAttr:=table(); TextFilesList:=[] end; sq end: ################################################################################ #hfl: NumConvert #toc: Numerical conversions #BasicTools/Round,mod1,Reduce2P,PeriodicMean,PeriodicStandardDeviation,angle2pair,ReduceFloat,ReduceFloat2,ReduceFloat3,ReduceFloatGrid Round:=proc(v,n::integer,$) if type(v,realcons) then evalf(round(v*10^n)*10^(-n)) elif type(v,indexable) then map(Round,v,n) elif type(v,complex) then simplify(Round(Re(v),n)+Round(Im(v),n)*I) else error "Unrecognized format",v end end: mod1:=proc(e::integer,m::posint,$) modp(e-1,m)+1 end: Reduce2P:=proc(x,P:=1,x0:=0,$) local e,i; e:=indets(x); x-P*floor(`if`(e={},(x-x0)/P,eval((x-x0)/P,{seq(e[i]=exp(-2.07-0.5*i),i=1..nops(e))}))) end: PeriodicMean:=proc(V::indexable,i1orx1::{list,set}:=[1],P:=1,x0:=0,{maxiter::posint:=3,digits::integer:=Digits-2,printout::boolean:=false},$) local n,maxdev,x1,x2,iter,v; n:=add(1,v=V); maxdev:=evalf(n*10^(-digits)); x1:=`if`(type(i1orx1,list),V[op(i1orx1)],op(i1orx1)); for iter from 1 to maxiter do x2:=Reduce2P(add(Reduce2P(v,P,x1-P/2),v=V)/n,P,x0); if (abs(x1-x2)parse(TrimRight(sprintf("%a",evalf(f)),["0"]),'statement'): ReduceFloat2:=proc(f::atomic,$) local s; s:=TrimRight(sprintf("%a",evalf(f)),["0"]); if (s[-1]=".") then s:=s[..-2] end; parse(s,'statement') end: ReduceFloat3:=proc(f::atomic,base::posint,digits::posint:=Digits-2,$) local v; v:=round(f*base); `if`(abs(f-v/base)<10^(-digits),v/base,f) end: ReduceFloatGrid:=proc(G0::list(list),id::{posint,list(posint)},digits::posint,digits2::posint:=digits-2,$) local lso,no,kmins,lsn,n,G,i,o,v; lso:=`if`(type(id,posint),[$1..id],id); no:=nops(lso); kmins:=[seq(MinVal(remove(`<`,[seq(abs(v[o]),v=G0),1],10^(-digits))),o=lso)]; lsn:=map(v->1/v,map(convert,kmins,rational,digits2)); for n in lsn do if not(type(n,integer)) then error("Noninteger result: %1",lsn) end end; G:=map(v->subsop(seq(lso[o]=round(v[lso[o]]*lsn[o])/lsn[o],o=1..no),v),G0); for i from 1 to nops(G) do for o in lso do if (abs(G[i][o]-G0[i][o])>10^(-digits)) then error("Deviation is too large for newG[%1][%2]=%3<>%4",i,o,G[i][o],G0[i][o]) end end end; G,lsn end: ################################################################################ #hfl: Express #toc: Expression tools #BasicTools/simplifyradical,coeffs2,idExch,idSym,idASym simplifyradical:=proc(e) local i; if has(e,I) then subs(i=I,factor(simplifyradical(simplify(Re(e)))+i*simplifyradical(simplify(Im(e))))) else signum(e)*sqrt(expand(rationalize((radnormal(e)^2)))) end end: coeffs2:=proc(e,vars::list,$) local e2,C,M,x,i; if (vars=[]) then return [[e]] end; e2:=collect(e,vars,'distributed'); if not(type(e2,polynom(anything,vars))) then error("Not a polynom of %1: %2",vars,e) end; C:=[coeffs(e2,vars,M)]; M:=[M]; [seq([seq(degree(M[i],x),x=vars),C[i]],i=1..nops(C))] end: idExch:=proc(e,A::name,id1,id2,$) local s,pre,L,i,j,k,v,x; s:=convert(e,string); pre:=cat(convert(A,string),"["); L:=length(pre); i:=1; do k:=searchtext(pre,s,i..-1); if (k=0) then break end; i:=i+k+L-2; j:=searchtext("]",s,i..-1)+i-1; if (j=i-1) then error "There is no right bracket after left bracket at position %1",i end; v:=parse(s[i..j]); v:=convert(subs(id1=x,id2=id1,x=id2,v),string); s:=cat(s[1..i-1],v,s[j+1..-1]) end; parse(s) end: idSym:=proc(e,A::name,i,j) local e2,k; if (nargs=4) then (e+idExch(e,A,i,j))/2 elif (nargs=5) then k:=args[5]; e2:=idExch(e,A,i,j): (e+idExch(e,A,i,k)+idExch(e,A,j,k)+e2+idExch(e2,A,i,k)+idExch(e2,A,j,k))/6 else error "No more than three indices are supported" end end: idASym:=proc(e,A::name,i,j,$) (e-idExch(e,A,i,j))/2 end: ################################################################################ #hfl: List #toc: List tools #BasicTools/SearchPos,Substitute,IdentifyPairs,Classify2,BlockSplit SearchPos:=proc(s::{string,list},xorf,$) local f,y,i; if (type(s,string) and not(type(xorf,procedure))) then return SearchText(xorf,s) end; if type(xorf,procedure) then f:=xorf else f:=y->y=xorf end; for i from 1 to `if`(type(s,string),length(s),nops(s)) do if f(s[i]) then return i end end; 0 end: Substitute:=proc(s::{string,list},xorf) local f,y,i; if (type(s,string) and not(type(xorf,procedure))) then return StringTools[Substitute](s,xorf,`if`(_nrest>0,_rest,"")) end; if type(xorf,procedure) then f:=xorf else f:=y->y=xorf end; for i from 1 to `if`(type(s,string),length(s),nops(s)) do if f(s[i]) then break end end; if (type(s,list) and i<=nops(s)) then subsop(i=_rest,s) elif (type(s,string) and i<=length(s)) then cat(s[..i-1],_rest,s[i+1..]) else s end end: Classify2:=proc(ls,ForId::{procedure,list},{index::boolean:=false,nolist::boolean:=false},$) local F,tb,nolist2,lsi,F2,i,id,v; if type(ForId,list) then F:=v->v[op(ForId)] else F:=v->ForId(v) end; tb:=table(); if index then if type(ls,list) then nolist2:=true; lsi:=[$1..nops(ls)] elif type(ls,table) then nolist2:=nolist; lsi:=[indices(ls,`if`(nolist2,':-nolist',NULL))] elif type(ls,rtable) then nolist2:=evalb(rtable_num_dims(ls)=1); lsi:=`if`(nolist2,map(v->lhs(v),rtable_elems(ls)),map(v->[lhs(v)],rtable_elems(ls))) else error("Indexing is supported only for list,table,rtable") end; if nolist2 then F2:=i->F(ls[i]) else F2:=i->F(ls[op(i)]) end; for i in lsi do id:=F2(i); if assigned(tb[id]) then tb[id]:=[op(tb[id]),i] else tb[id]:=[i] end end; op(tb) elif type(ls,list) then for v in ls do id:=F(v); if assigned(tb[id]) then tb[id]:=[op(tb[id]),v] else tb[id]:=[v] end end; op(tb) else ListTools[Classify](F,ls) end end: IdentifyPairs:=proc(ls1::indexable, ls2::indexable, idf::procedure, output::string:="l", {nolist::boolean:=false}, $) local id1,id2,idv,i1,i2,n,V,i,k1,k2,R,c; id1,id2:=`if`(type(ls1,list),{seq([i],i=1..nops(ls1))},{indices(ls1)}),`if`(type(ls2,list),{seq([i],i=1..nops(ls2))},{indices(ls2)}); idv:=table(); for i1 in id1 do for i2 in id2 do idv[i1,i2]:=idf(ls1[op(i1)],ls2[op(i2)]) end end; n:=min(nops(id1),nops(id2)); V:=Vector(n); for i from 1 to n do k1,k2:=op(MinIdx(idv)); # V[i]:=[`if`(nolist,op(k1),k1),`if`(nolist,op(k2),k2),idv[k1,k2]]; V[i]:=[k1,k2,idv[k1,k2]]; for i2 in id2 do idv[k1,i2]:=evaln(idv[k1,i2]) end; id1:=`minus`(id1,{k1}); for i1 in id1 do idv[i1,k2]:=evaln(idv[i1,k2]) end; id2:=`minus`(id2,{k2}) end; R:=[op(id1),op(id2)]; seq(`if`(c="l", [seq(`if`(nolist,[op(V[i][1]),op(V[i][2]),V[i][3]],V[i]),i=1..n)], `if`(c="v", [seq([ls1[op(V[i][1])],ls2[op(V[i][2])],V[i][3]],i=1..n)], `if`(c="r", `if`(nolist,map(op,R),R), `if`(c="m", `if`(n=0,0,V[n][3]), NULL)))),c=output) end: BlockSplit:=proc(L::list,B::list(nonnegint),$) local N,n,A,i,R,v; N,n:=nops(L),nops(B); A:=Vector(B); if (A[1]>=N) then [L,[]$(n-1)] else for i from 2 to n do A[i]:=A[i-1]+A[i]; if (A[i]>N) then A[i]:=N; n:=i; break end end; R:=[[1,min(N,A[1])],seq([A[i-1]+1,A[i]],i=2..n)]; [seq(L[v[1]..v[2]],v=R)] end end: ################################################################################ #hfl: Grid #toc: Grids and sequences #BasicTools/Grid,Grid2,GetPath,convert2range,convert2repeat Grid:=proc(rng::complexcons..complexcons,n::integer,f::procedure:=(x->x),avoid::list:=[],$) local x1,x2,x,m,k,t,t1,t2,ls,lsex; x1,x2:=op(rng); m:=`if`(n<2,1,n-1); lsex:=avoid; if (f(x)=x) then ls:=[seq(x1+(x2-x1)*k/m,k=0..m)] else if (type(x1,float) or type(x2,float)) then lsex:=evalf(lsex); t1,t2:=max(fsolve(f(t)=x1,t)),max(fsolve(f(t)=x2,t)); else t1,t2:=max(solve(f(t)=x1,t)),max(solve(f(t)=x2,t)) end; ls:=[seq(f(t1+(t2-t1)*k/m),k=0..m)] end; remove(member,ls,lsex) end: Grid2:=proc(lsp::listlist,dx::numeric,c::numeric:=1.1,$) local e,n,ls,v,i; if (nops(lsp)=1) then [map(v->round(v/dx)*dx,lsp[1])] elif (nops(lsp)=2) then n:=round(c*sqrt(add(v^2,v=lsp[2]-lsp[1]))/dx); if (n=0) then Grid2([(lsp[1]+lsp[2])/2],dx,c) else e:=(lsp[2]-lsp[1])/n; ListTools[MakeUnique]([seq(map(v->round(v/dx)*dx,lsp[1]+e*i),i=0..n)]) end else ls:=[seq(Grid2([lsp[i-1],lsp[i]],dx,c),i=2..nops(lsp))]; [seq(op(..-2,v),v=ls[..-2]),op(ls[-1])] end end: GetPath:=proc(X::{list,Vector},Y::{list,Vector},N::list(nonnegint),$) local d,o,XY,n,i; d:=nops(N); if not(foldl(`and`,seq(type(X[o]*N[o],integer),o=1..d))) then error("Non-grid X: %1, but N=%2",X,N) end; if not(foldl(`and`,seq(type(Y[o]*N[o],integer),o=1..d))) then error("Non-grid Y: %1, but N=%2",Y,N) end; XY:=Y-X; n:=igcd(seq(XY[o]*N[o],o=1..d)); [seq(X+i/n*XY,i=0..n)] end: convert2range:=proc(ls0::list(integer),minnum::posint:=3,{asstring::boolean:=false},$) local ls,tb,i1,i,v; ls:=ls0; if (nops(ls)>1) then tb:=table(); i1:=1; for i from 2 to nops(ls) do if (ls[i]-ls[i-1]>1) then tb[i1]:=`if`(i=i1+1,ls[i1],[ls[i1],ls[i-1]]); i1:=i elif (ls[i]-ls[i-1]<1) then error("Nonmonotonic at i=%1 list %2",i,ls) end end; tb[i1]:=`if`(i=i1+1,ls[i1],[ls[i1],ls[i-1]]); ls:=map(v->`if`(type(v,list) and v[2]-v[1]+1`if`(type(v,list),sprintf("$%d..%d",op(v)),sprintf("%d",v)),ls))) else ls end end: convert2repeat:=proc(ls::list,minnum::posint:=3,format::string:="",$) local tb,i,L,k,v; if (ls=[]) then L:=[] else tb:=table([1=1]); for i from 2 to nops(ls) do if (ls[i]<>ls[i-1]) then tb[i]:=1 end end; L:=sort([indices(tb,'nolist'),nops(ls)+1]); L:=[seq([ls[L[k]],L[k+1]-L[k]],k=1..nops(L)-1)]; L:=map(v->`if`(v[2]`if`(type(v,list),cat(sprintf(format,v[1]),"$",v[2]),sprintf(format,v)),L))) end end: ################################################################################ #hfl: Tree #toc: Trees #BasicTools/PrintTree,Tree2Seq,Tree2Seq2,Tree2List PrintTree:=proc( tree::list, level::nonnegint:=0, str::procedure:=(x->x), { filedesc::{string,integer}:="", prefix::string:="." },$) local s,tree1; s:=cat(prefix$level,str(tree[1]),"\n"); if (filedesc="") then printf(s) else fprintf(filedesc,s) end; for tree1 in tree[2..-1] do PrintTree(tree1,level+1,str,':-filedesc'=filedesc,':-prefix'=prefix) end; NULL end: Tree2Seq:=proc( tree::list, prefix::string:="", str::procedure:=(x->x), {separator::string:="/"},$) local tree1; cat(prefix,str(tree[1])),seq(Tree2Seq(tree1,cat(prefix,str(tree[1]),separator),str,':-separator'=separator),tree1=tree[2..]) end: Tree2Seq2:=proc( tree::list, index::list(posint):=[], size::posint:=1, body::procedure:=(x->x),$) local i; [body(`if`(size=1,tree[1],tree[..size])),index,nops(tree)-size],seq(Tree2Seq2(tree[size+i],[op(index),i],size,body),i=1..nops(tree)-size) end: Tree2List:=proc( tree::list, size::posint:=1, body::procedure:=(x->x),$) local ls,id,i; ls:=[Tree2Seq2(tree,size,body)]; id:=table([seq(op(ls[i][2])=i,i=1..nops(ls))]); map(v->[v[1],`if`(v[2]=[],0,id[op(..-2,v[2])]),[seq(id[op(v[2]),i],i=1..v[3])]],ls) end: ################################################################################ #hfl: idx #toc: Indexing #BasicTools/Dim2,idL,iidL,idI,iidI,idS,iidS,idA,iidA,idSS,iidSS,nextpointS Dim2:=(A::{rtable,list})->`if`(type(A,list),nops(A),`if`(type(A,Array),op(map2(op,2,[op(2,A)])),op(1,A))): idL:=proc(x::list(integer),L::{posint,list(posint)},$) local d,LL,n,i; d:=nops(x); LL:=`if`(type(L,posint),[L$d],L); n:=x[d] mod LL[d]; for i from d-1 by -1 to 1 do n:=n*LL[i]+(x[i] mod LL[i]) end; n+1 end: iidL:=proc(n::posint,L::list(posint),$) local d,x,p,i; d:=nops(L); x:=array(1..d); p:=n-1; for i from 1 to d do x[i]:=irem(p,L[i],'p') end; convert(x,list) end: idI:=m->2*abs(m)+`if`(m>0,0,1): iidI:=n->`if`(type(n,odd),-1,1)*iquo(n,2): idS:=proc(x::list(posint),$) local i,j; i:=max(x[1],x[2]); j:=min(x[1],x[2]); i*(i-1)/2+j end: iidS:=proc(n::posint,$) local i; i:=round(sqrt(2*n)); [i,n-i*(i-1)/2] end: idA:=proc(x::list(posint),$) local i,j; if (x[1]=x[2]) then 0 else i:=max(x[1],x[2]); j:=min(x[1],x[2]); (i-1)*(i-2)/2+j end end: iidA:=proc(n::posint,$) local i; i:=round(sqrt(2*n))+1; [i,n-(i-1)*(i-2)/2] end: idSS:=proc(x::list(posint),$) local n,m; n:=idS(x[1..2]); m:=idS(x[3..4]); [max(n,m),min(m,n)] end: iidSS:=nm->[op(iidS(nm[1])),op(iidS(nm[2]))]: nextpointS:=proc(x::list(nonnegint),$) local d,i; d:=nops(x); for i from d by -1 to 2 while (x[i-1]=x[i]) do end; [op(..i-1,x),x[i]+1,0$(d-i)] end: ################################################################################ #hfl: Sort #toc: Sorting #BasicTools/isleq,Sort,Rort,SortM,SortIdx,MaxVal,MinVal,MaxIdx,MinIdx,SortEE,SortSign isleq:=proc(u,v) local n,i; if type(u,list) then n:=min(nops(u),nops(v)); for i from 1 to n while (u[i]=v[i]) do end; if (i>n) then true else isleq(u[i],v[i]) end else evalb(uu),id::list:=[],$)::list; local F2,V; if type(F,procedure(anything,anything)) then F2:=`if`(id=[],F,(u,v)->F(u[op(id)],v[op(id)])); sort(L,F2) else V:=`if`(id=[],map(F,L),map(v->F(v[op(id)]),L)); L[sort(V,'output=permutation')] end end: Rort:=proc(L::list,F::procedure:=(u->u),id::list:=[],$)::list; local F2,V; if type(F,procedure(anything,anything)) then F2:=`if`(id=[],(u,v)->F(v,u),(u,v)->F(v[op(id)],u[op(id)])); sort(L,F2) else V:=`if`(id=[],map(F,L),map(v->F(v[op(id)]),L)); L[sort(V,`>`,'output=permutation')] end end: #Sort:=proc(ls::list,F::procedure:=(u->u),id::list:=[],$)::list; #local F2; #if (id=[]) then F2:=`if`(type(F,procedure(anything,anything)),F,(u,v)->evalb(F(u)F(u[op(id)],v[op(id)]),(u,v)->evalb(F(u[op(id)])u),$)::list; #local F; #if type(ForId,list) then F:=(u,v)->evalb(u[op(ForId)]evalb(ForId(u)sort(L,proc(u,v) local i; for i in id while (u[i]=v[i]) do end; isleq(u[i],v[i]) end): SortIdx:=proc(L::indexable,F::procedure:=(u->u),id::list:=[],{nolist::boolean:=false},$)::list; local V,J,i,F2,P; if type(L,list) then V:=L; J:=[[i]$i=1..nops(V)] else V:=convert(L,table); J:=[indices(V)]; V:=[seq(V[op(i)],i=J)] end; if type(F,procedure(anything,anything)) then F2:=`if`(id=[],F,(u,v)->F(u[op(id)],v[op(id)])); P:=sort(V,F2,'output=permutation') else V:=`if`(id=[],map(F,V),map(v->F(v[op(id)]),V)); P:=sort(V,'output=permutation') end; J:=J[P]; `if`(nolist,map(op,J),J) end: #SortIdx:=proc(ls::indexable, ForId:=(u->u), {nolist::boolean:=false}, $)::list; #local tb,i,id,F,ls2,v; #if type(ls,list) then # tb:=ls; id:=[[i]$i=1..nops(tb)] #else # tb:=convert(ls,table); id:=[indices(tb)] end; #if type(ForId,list) then F:=(u,v)->evalb(u[1][op(ForId)]ForId(u[1],v[1]) #elif type(ForId,procedure) then F:=(u,v)->evalb(ForId(u[1])u),$)::list; local V,x,y; V:=`if`(type(ls,list),ls,[entries(convert(ls,table),'nolist')]); if (V=[]) then return [] end; x:=V[1]; if type(ForId,list) then for y in V do if (x[op(ForId)]u),$)::list; local V,x,y; V:=`if`(type(ls,list),ls,[entries(convert(ls,table),'nolist')]); if (V=[]) then return [] end; x:=V[1]; if type(ForId,list) then for y in V do if (y[op(ForId)]u),$)::list; local tb,id,i,j,x; if type(ls,list) then tb:=ls; id:=[[i]$i=1..nops(tb)] else tb:=convert(ls,table); id:=[indices(tb)] end; if (id=[]) then return [] end; i:=id[1]; x:=tb[op(i)]; if type(ForId,list) then for j in id do if (x[op(ForId)]u),$)::list; local tb,id,i,j,x; if type(ls,list) then tb:=ls; id:=[[i]$i=1..nops(tb)] else tb:=convert(ls,table); id:=[indices(tb)] end; if (id=[]) then return [] end; i:=id[1]; x:=tb[op(i)]; if type(ForId,list) then for j in id do if (tb[op(j)][op(ForId)]>, F::procedure:=(u->u),$) local P,v; P:=SortIdx([seq(v,v=ev)],F,'nolist'); ev[P],`if`(evc=<<0>>,NULL,evc[..,P]) end: SortSign:=proc(ls::{list,Vector},$) local V,v,i,j,s; V:=convert(ls,Vector); s:=0; for i from 2 to op(1,V) do v:=V[i]; for j from i-1 by -1 to 1 while (V[j]>v) do V[j+1]:=V[j] end; s:=s+i-1-j; V[j+1]:=v end; `if`(type(ls,list),convert(V,list),V), `if`(type(s,odd),-1,1) end: ################################################################################ #cat: String tools #hfl: specsymbols #specsymbols #hfl: Trim Squeeze:=proc(s::string,key::string:=" ", multiplicity::posint:=2, $) local dkey,s1,s2,i,j; dkey:=cat(key$multiplicity); s1:=s: s2:="": i:=searchtext(dkey,s1); while (i>0) do s2:=cat(s2,s1[1..i]); j:=i: while (s1[j+1]=key) do j:=j+1 end: s1:=s1[j+1..-1]; i:=searchtext(dkey,s1) end: cat(s2,s1); end: #hfl: Trim TrimLeft:=proc(s::string, ls::list:=[" "]) local i; for i from 1 to length(s) do if not(s[i] in ls) then break end end; s[i..-1] end: #hfl: Trim TrimRight:=proc(s::string, ls::list:=[" "]) local i; for i from length(s) by -1 to 1 do if not(s[i] in ls) then break end end; s[1..i] end: #hfl: Trim Trim:=proc(s::string, ls::list:=[" "]) TrimRight(TrimLeft(s,ls),ls) end: #hfl: StartNewLine StartNewLine:=proc(s::string,ls::list:=[],$) local s1,s2; s1:=s; for s2 in ls do s1:=StringTools[SubstituteAll](s1,s2,cat("\n",s2)) end; s1 end: #hfl: SearchKeySeq SearchKeySeq:=proc(s::string, keyseq::list(string), { start::posint:=1, skip::nonnegint:=0, sep::realcons:=infinity, fulloutput::boolean:=false, noerror::boolean:=false },$) local i,n,j,lsp,k; i:=start; for n from 0 to skip do j:=searchtext(keyseq[1],s,i..-1); if (j=0) then if noerror then return `if`(fulloutput,[],0) else printf(s); error("No key%1 \"%2\" is found after %3/%4 skips",1,keyseq[1],n,skip) end end; i:=i+j-1+length(keyseq[1]) end; lsp:=array(1..nops(keyseq)): lsp[1]:=i; for k from 2 to nops(keyseq) do j:=searchtext(keyseq[k],s,i..-1); if (j=0) then if noerror then return `if`(fulloutput,[],0) else printf(s); error("No key%1 \"%2\" is found",k,keyseq[k]) end end; if (j>sep+1) then if noerror then return `if`(fulloutput,[],0) else printf(s); error("key%1 \"%2\" is separated by %3>%4 positions",k,keyseq[k],j-1,sep) end end; i:=i+j-1+length(keyseq[k]); lsp[k]:=i end; `if`(fulloutput,convert(lsp,list),i) end: #hfl: ReadBody ReadBody:=proc( s::string, key1::string, key2::string:="", trim::list:=[], { noerror::boolean:=false, restoftheline::boolean:=false, format::string:="", wkeys::boolean:=false, pos1::name:=undefined, pos2::name:=undefined},$) local i,j,s2,v; i:=SearchText(key1,s); if (i=0) then if noerror then return "" else error( "No key1 (\"%1\") in %2",key1,s) end else if not(type(pos1,undefined)) then assign(pos1,i-1) end; i:=i+length(key1) end; j:=`if`(key2="",length(s),SearchText(key2,s,i..-1)+i-2); if (j=i-2 and not(key2="")) then if noerror then if restoftheline then j:=length(s) else return "" end else error("No key2 (\"%1\") in %2",key2,s) end end; if not(type(pos2,undefined)) then assign(pos2,j+length(key2)+1) end; s2:=`if`(trim=[],s[i..j],Trim(s[i..j],trim)); if (format="") then `if`(wkeys,cat(key1,s2,key2),s2) else v:=sscanf(s2,format); `if`(v=0 or v=[],s2,op(v)) end end: #hfl: ReadBody ReadByMask:=proc(mask::string,s::string,star::character:="*",$) local i,i1,s1; if (mask[1]="*") then if (mask="*") then s else i:=SearchText(star,mask[2..]); if (i=1) then error("Consecutive stars are not allowed in mask: %1",mask) end; s1:=mask[2..`if`(i=0,-1,i)]; i1:=SearchText(s1,s); if (i1=0) then error("No suffix from mask (%1) in search string (%2)",mask,s) end; s[..i1-1],`if`(i=0,NULL,ReadByMask(mask[i+1..],s[i1+length(s1)..],star)) end else i:=SearchText(star,mask); if (i=0) then if (mask=s) then NULL else error("No match: %1<>%2",mask,s) end else if (mask[..i-1]=s[..i-1]) then ReadByMask(mask[i..],s[i..],star) else error("No prefix from mask (%1) in search string (%2)",mask,s) end end end end: #hfl: Split2 Split2:=proc(s::string,c::string:=" ",{output::string:="s"},$) local i1,n,lsi,i,lss,k,o; if (length(c)<>1) then error("Unrecognized split character: \"%1\"",c) end; i1:=0; n:=length(s); lsi:=table([1=1,n+1=n+1]); do for i1 from i1+1 to n while (s[i1]=c) do end; if (i1>n) then break end; i:=SearchText(c,s,i1..n); if (i>0) then i1:=i1+i-1; lsi[i1]:=i1 else break end end; lsi:=convert(lsi,list); lss:=[seq(s[lsi[k-1]..lsi[k]-1],k=2..nops(lsi))]; seq(`if`(o="s",lss,`if`(o="i",lsi,NULL)),o=output) end: #hfl: FormattedOperation FormattedOperation:=proc(s::string,f::procedure,pos::integer,iorn::integer,$) local i2,i1,s1,ls,S1,n,V1,V2,S2,k,i,j,d,s2,v; i2:=length(s); if (iorn>0) then i2:=min(i2,iorn) end; if (pos>0) then i1,i2,s1:=pos,i2,s[pos..i2] else i1,i2,s1:=1,i2,cat(" "$(-pos),s[..i2]) end; S1:=Split2(s1); if (nops(S1)>1 and Trim(S1[-1])="") then S1:=[op(..-3,S1),cat(S1[-2],S1[-1])] end; if (iorn>0) then n:=nops(S1) elif (iorn=0) then return cat(s1,s[i2+1..]) else n:=-iorn; if (n>nops(S1)) then error("Expected %1 records, detected %2 in %3",n,nops(S1),s) end; S1:=S1[..n]; i2:=i1+add(length(v),v=S1)-1+min(0,pos) end; V1:=Vector(map(op@sscanf,S1,"%f")); V2:=f(V1); S2:=Vector(n); for k from 1 to n do i:=SearchText(".",S1[k]); for j from i+1 to length(S1[k]) while StringTools[IsDigit](S1[k][j]) do end; d:=j-i-1; s2:=sprintf("%.*f%s",d,V2[k],S1[k][j..]); j:=SearchText(".",s2); if (j>i) then WARNING("String %1 is expanded by %2 symbols",S1[k],j-i) end; S2[k]:=cat(" "$max(0,i-j),s2) end; cat(s[..i1-1],seq(v,v=S2),s[i2+1..]) end: #hfl: FormatFloat FormatFloat:=proc( v::realcons, d1::{nonnegint,realcons}, d2::posint:=`if`(type(d1,nonnegint),max(1,d1),9), { maxorder::nonnegint:=3, asfloat::boolean:=false, width::nonnegint:=0, alignsign::boolean:=false, showpoint::boolean:=false, showzero::boolean:=false, hidezero::boolean:=false, showsign::boolean:=false, hidenull::boolean:=false },$) local a,rmin,ppos,r,s; if type(d1,nonnegint) then if (not(asfloat) and type(v,integer)) then r:=sprintf("%d",abs(v)); if showpoint then r:=cat(r,"."); ppos:=0 else ppos:=-1 end elif (not(asfloat) and type(v,rational) and type(v*10^(d1+maxorder),integer)) then for ppos from d1 to d1+maxorder do if type(v*10^ppos,integer) then break end end; r:=FormatFloat(evalf(abs(v)),ppos,d2,':-maxorder'=0,':-showzero'=showzero,':-hidezero'=true) elif (not(asfloat) and type(v,float) and -op(2,v)<=d1+maxorder) then r:=sprintf("%.*f",-op(2,v),abs(v)); if (not(showzero) and r[..2]="0.") then r:=r[2..] end; ppos:=max(0,-op(2,v)) else a:=abs(evalf(v)); rmin:=10^(d2-1); for ppos from d1 to d1+maxorder do r:=round(a*10^ppos); if (r>=rmin) then break end end; if (ppos>d1+maxorder) then ppos:=d1+maxorder end; r:=sprintf("%0*d",ppos,r); if (length(r)>ppos) then if (ppos=0 and d1=0) then if showpoint then r:=cat(r,".") else if (width>0) then r:=cat(r," ") end end else r:=cat(r[..-ppos-1],".",r[-ppos..]) end else r:=cat(`if`(showzero,"0.","."),r) end end else for ppos from 0 to infinity while (round(d1*10^ppos)<=d2) do end; ppos:=ppos-1; r:=sprintf("%.*f(%.0f)",ppos,abs(v),d1*10^ppos); if (not(showzero) and r[..2]="0.") then r:=r[2..] end end; if (evalf(v)<0) then r:=cat("-",r) elif showsign then r:=cat("+",r) end; r:=cat(r," "$min(d1+maxorder-ppos,width-length(r))); s:=cat(""," "$(width-length(r))); s:=`if`(alignsign and (r[1]="+" or r[1]="-"),cat(r[1],s,r[2..]),cat(s,r)); if (hidenull and parse(s)=0) then s:=`if`(width=0,"",cat(" "$width)) end; `if`(hidezero and ppos>0,TrimRight(s,["0"]),s) end: #hfl: FormatFloat AlignPoints:=proc(ls::list(string),$) local n,V,i,s,p,l1,l2,f,v; n:=nops(ls); V:=Vector(n): for i from 1 to n do s:=ls[i]; p:=searchtext(".",s); V[i]:=`if`(p=0,[s," ",""],[s[..p-1],".",s[p+1..]]) end; l1,l2:=max(seq(length(V[i][1]),i=1..n)),max(seq(length(V[i][3]),i=1..n)); f:=cat("%",l1,"s%s%-",l2,"s"); [seq(sprintf(f,op(v)),v=V)] end: #hfl: FormatTable FormatTable:=proc( tb::table, format::string:=" %a", { prefix::string:="", align::boolean:=false, printout::boolean:=false },$) local ls,N,n,S,p,V,c,dp,i,islist,s,v; ls:=[indices(tb)]; N:=nops(ls); n:=min(seq(nops(v),v=ls)); ls:=SortM(ls,[$1..n]); S:=Vector(N,i->convert(ls[i],string)[2..-2]); if align then for i from 1 to N do S[i]:=cat(S[i],",") end; p:=1; do V:=[seq(`if`(p>length(s),0,SearchText(",",s,p..-1)),s=S)]; dp:=max(op(V)); if (dp=0) then break end; for i from 1 to N do if (V[i]"%a"); ls:=[seq(sprintf(cat("( %s )=",format),S[i],`if`(islist,op(tb[op(ls[i])]),tb[op(ls[i])])),i=1..N)]; if printout then for s in ls do printf("%s\n",s) end end; ls end: #hfl: FormatTime2 FormatTime2:=proc(t::numeric,{nodays::boolean:=false},$) local s,m,h,d; s:=irem(round(t),60,'m'); if (m=0) then return sprintf("%d",s) end; m:=irem(m,60,'h'); if (h=0) then return sprintf("%d:%02d",m,s) end; if (nodays or h<24) then return sprintf("%d:%02d:%02d",h,m,s) end; h:=irem(h,24,'d'); sprintf("%d:%02d:%02d:%02d",d,h,m,s) end: #hfl: PrintColumns PrintColumns:=proc(ls::list,n0::integer:=0,sep::string:=" ",{before::string:=""},$) local m,n,M,L,i,j; if (ls=[]) then return NULL end; if type(ls,list(string)) then if (n0>0) then n:=n0 else n:=round(sqrt(nops(ls))); if (n0<0) then n:=min(-n0,n) end end; m:=ceil(nops(ls)/n); n:=ceil(nops(ls)/m); M:=Matrix(m,n,(i,j)->`if`(i+(j-1)*m>nops(ls),"",ls[i+(j-1)*m])) elif type(ls,list(list(string))) then n:=nops(ls); m:=max(map(nops,ls)); M:=Matrix(m,n,(i,j)->`if`(i>nops(ls[j]),"",ls[j][i])) else error("Unrecognized ls: %1",ls) end; L:=Vector(n,j->max(seq(length(M[i,j]),i=1..m))); for i from 1 to m do printf(before); for j from 1 to n-1 do printf("%-*s%s",L[j],M[i,j],sep) end; printf("%-*s\n",L[n],M[i,n]) end; end: #hfl: PrintMatrix PrintMatrix:=proc( M::Matrix, fmt::{integer,string}:=-interface(displayprecision), corner::string:="", { maxsize::posint:=99, head::{boolean,[indexable,indexable]}:=false, left::boolean:=false, right::boolean:=false, Ichar::character:="i", subszero::string:="donotsubs", donotprint::boolean:=false },$) local pri,n,m,d,M1,M2,n1,d1,zero,print0,l,W,C1,R1,w1,v,i,j; pri:=not(donotprint); n,m:=Dim2(M); if (max(n,m)>maxsize) then pri:=false; WARNING("Will not print because max size of %1 is exceeded: [%2,%3]",maxsize,n,m) end; if type(fmt,negint) then d:=-fmt; M1:=map(ReduceFloat2,Matrix(map(Round,M,d),datatype=anything)); n1,d1:=interface(rtablesize,displayprecision); if pri then interface(rtablesize=max(n,m),displayprecision=d); print(M1); interface(rtablesize=n1,displayprecision=d1) end else if type(fmt,string) then d:=3; M1:=map2(sprintf,fmt,M) else d:=fmt; zero:=cat("."," "$d); print0:=s->`if`(member(TrimRight(Trim(s),["0"]),["0.","-0."]),zero,s); if (MatrixOptions(M,datatype)=complex[8]) then M1,M2:=seq(map(v->print0(sprintf("%.*f",d,v)),f(M)),f=[Re,Im]); M1:=Matrix(n,m,(i,j)->`if`(M2[i,j]=zero,cat(M1[i,j]," "),`if`(M1[i,j]=zero,cat(M2[i,j],Ichar),cat(M1[i,j],`if`(M2[i,j][1]="-","","+"),M2[i,j],Ichar)))) else M1:=map(v->print0(sprintf("%.*f",d,v)),M) end end; if (subszero<>"donotsubs") then M1:=map(v->`if`(`subset`(convert(v,set),{" ","+","-",".","0"}),subszero,v),M1) end; W:=Vector(m,j->max(map(length,M1[..,j])),datatype=integer); if (head=false) then M1:=Matrix(n,m,(i,j)->StringTools[PadLeft](M1[i,j],W[j])) else if type(head,list) then C1,R1:=op(head) else C1:=Vector(n,i->sprintf("%3d",i)); R1:=Vector(m,j->sprintf( "%d",j)) end; C1:=Vector(n+1,i->`if`(i=1,corner,C1[i-1])); w1:=max(map(length,C1)); W:=Vector(m,j->max(W[j],length(R1[j])),datatype=integer); M1:=Matrix(n+1,m+1,(i,j)->`if`(i=1, `if`(j=1,StringTools[PadRight](C1[1],w1),StringTools[PadLeft](`if`(right,R1[j-1],StringTools[PadRight](StringTools[Center](R1[j-1],3),d)),W[j-1])), `if`(j=1,`if`(left,StringTools[PadRight],StringTools[PadLeft])(C1[i],w1),StringTools[PadLeft](M1[i-1,j-1],W[j-1])))) end; if pri then printf("%s\n",M1) end end; M1 end: #hfl: PrintMatrix PrintUpperMatrix:=proc(M::Matrix,fmt::string,width::posint:=16,height::posint:=99,{step::posint:=1},$) local n,m,i1,i; n,m:=op(1,M); for i1 from 1 to step do for i from i1 by step to min(height,n,m) do printf(cat(fmt,"\n"),Vector(min(width,m-i+1),j->M[i,i+j-1])) end end end: #hfl: PrintVector3d PrintVector3d:=proc(V::indexable,$) local U,i1,i2,i3,k,v; U:=[seq(abs(v),v=V)]; i1,i2,i3:=op(SortIdx(U,'nolist')); k:=op(MaxIdx([U[i3],0.70710678*(U[i3]+U[i2]),0.57735027*(U[i3]+U[i2]+U[i1])])); if (k=1) then `if`(V[i3]>0,"XYZ","xyz")[i3] elif (k=2) then sort(cat(`if`(V[i3]>0,"XYZ","xyz")[i3],`if`(V[i2]>0,"XYZ","xyz")[i2])) else cat(`if`(V[1]>0,"X","x"),`if`(V[2]>0,"Y","y"),`if`(V[3]>0,"Z","z")) end end: #hfl: ascii2unicode ascii2unicode:=proc(s::string,CharMap::table:=CyrCharMap,imax::posint:=127,$) local ls,i; ls:=convert(s,'bytes'); for i in ls do if (i>imax) then return cat(seq(CharMap[i],i=ls)) end end; s end: ################################################################################ #cat: Combinatorics and algorithms #hfl: CyclesAndTrees CyclesAndTrees:=proc( G::Graph, istree::procedure:=defistree, extracycles::list(set):=[], {output::{1,2}:=1},$) local nn,cycles,nc,nc1,trees,nt,cbonds,tbonds,p,q,i,j,tb,b; nn:=GraphTheory[Neighbors](G); cycles:=map(convert,GraphTheory[CycleBasis](G),set); cycles:=[op(cycles),op(extracycles)]; nc:=nops(cycles); while (nc<>nc1) do cycles:=map(`union`@op,[ListTools[Categorize]((u,v)->`intersect`(u,v)<>{},cycles)]); nc1,nc:=nc,nops(cycles) end; cycles:=map(convert,cycles,list); trees:=GraphTheory[ConnectedComponents](GraphTheory[DeleteVertex](G,map(op,cycles))); nt:=nops(trees); cycles:=Vector(nc,p->cycles[p]); cbonds,tbonds:=Vector(nc,p->[]),Vector(nt,p->[]): for p from 2 to nc do for i in cycles[p] do for j in nn[i] do if not(member(j,cycles[p])) then for q from 1 to p-1 do if member(j,cycles[q]) then cbonds[p]:=[op(cbonds[p]),[i,j,q]]; cbonds[q]:=[op(cbonds[q]),[j,i,p]]; break end end end end end end; for p from 1 to nt do for i in trees[p] do for j in nn[i] do if not(member(j,trees[p] )) then for q from 1 to nc do if member(j,cycles[q]) then tbonds[p]:=[op(tbonds[p]),[i,j,q]]; break end end end end end end; tb:=table(): for p from 1 to nt do if (nops(tbonds[p])=1 and not(istree(trees[p],op(1..2,tbonds[p][1]),cycles[tbonds[p][1][3]]))) then cycles[tbonds[p][1][3]]:=[op(cycles[tbonds[p][1][3]]),op(trees[p])] else tb[p]:=[trees[p],op(tbonds[p])] end end; trees:=convert(tb,list); nt:=nops(trees); for p from 1 to nt do for b in trees[p][2..] do cbonds[b[3]]:=[op(cbonds[b[3]]),[b[2],b[1],nc+p]] end end; [seq([sort(cycles[p]),op(cbonds[p])],p=1..nc),op(trees)],`if`(output=1,NULL,nc) end: defistree:=((tree,i,j,cycle)->evalb(nops(tree)>1)): #hfl: CycleBasis2 CycleBasis2:=proc(G::Graph,{maxiter::posint:=9},$) local ff,C,B,N,T,iter,p,q,v; ff:=f->select(v->nops(select(member,C[v],f))>2,f); C:=table(zip(`=`,GraphTheory[Vertices](G),GraphTheory[Neighbors](G))); B:=GraphTheory[CycleBasis](G); N:=nops(B); B:=Vector(N,p->{op(B[p])}); T:=map(ff,B); for iter from 0 to maxiter while (max(map(nops,T))>0) do for p from 1 to N do if (T[p]<>{}) then for q from 1 to N do if (q<>p and `subset`(B[q],B[p])) then B[p]:=`minus`(B[p],`minus`(B[q],T[p])); T[p]:=ff(B[p]) end end end end end; B:=[seq([op(v)],v=B)]; if (iter>maxiter) then error("No convergence, B=%1",B) end; B end: ################################################################################ #cat: Geometry #hfl: Circumsphere Circumsphere:=proc() local n,p0,P,i,p,v,o; n:=_nrest-1; if (n<2) then error("At least three vectors must be provided") end; p0:=_rest[1]; P:=Matrix(n); for i from 1 to n do P[i,..]:=_rest[1+i]-p0 end; try p:=LinearAlgebra[LinearSolve](P,Vector(n,i->add(P[i,o]^2,o=1..n)/2)); catch: return undefined,infinity end; p+p0,sqrt(add(v^2,v=p)) end: ################################################################################ #cat: Calculus #hfl: fullparfrac2 fullparfrac2:=proc(e,z::name,$) local e1,e2,a,subex,indet,eq,sol,z1; e1:=convert(e,fullparfrac,z); e2:=0: if not(type(e1,`+`)) then e1:=[e1] end; for a in e1 do if (type(a,`*`) or type(a,`^`)) then e2:=e2+a else subex:=op(1,a); indet:=op(1,op(2,a)); eq:=op(2,op(2,a)); sol:=allvalues(eq); for z1 in sol do e2:=e2+subs(indet=z1,subex) end end end: e2 end: #hfl: fZero fZero:=proc(f::procedure, r::range(realcons), iniguess::realcons:=infinity, {abserror::nonnegative:=0, relerror::nonnegative:=0, maxiter::posint:=100, output::name:=undefined},$) local x1,x2,x0,iter,flag; x1,x2:=op(r); if (x1=-infinity) then x1:=-10^Digits else x1:=evalf(x1) end; if (x2= infinity) then x2:= 10^Digits else x2:=evalf(x2) end; x0:=evalf(iniguess); if not(x0>=x1 and x0<=x2) then x0:=(x1+x2)/2 end; flag,x0,x2,iter:=FZERO(f,x1,x2,x0,relerror,abserror,maxiter); if not(type(output,undefined)) then assign(output,[flag,x0,x2,iter]) end; if (flag=1 or flag=2) then x0 elif (flag=3) then x0 elif (flag=4) then error "No change in sign was found. After %1 iterations x0=%2, x2=%3, dx=%4, f(x0)=%5", iter,x0,x2,x2-x0,evalf(f(x0)) elif (flag=5) then error "Too many iterations: %1. At this moment x0=%2, x2=%3, dx=%4, f(x0)=%5", iter,x0,x2,x2-x0,evalf(f(x0)) else error "Unrecognized flag: %1",flag end end: #hfl: fZero fDblZero:=proc(f::procedure, r::range(realcons), {abserror::nonnegative:=0, relerror::nonnegative:=0, maxiter::posint:=100, noerror::boolean:=false, rootnumber::posint:=0},$) local x1,x2,f1,f2,s,dx,x0,f0,y1,y2,g1,g2; x1,x2:=min(op(r)),max(op(r)); if (x1=-infinity) then x1:=-10^Digits else x1:=evalf(x1) end; if (x2= infinity) then x2:= 10^Digits else x2:=evalf(x2) end; f1:=evalf(f(x1)); f2:=evalf(f(x2)); s:=signum(f1); if (s*f2<0) then if (rootnumber=1 or rootnumber=2) then fZero(f,x1..x2,':-abserror'=abserror,':-relerror'=relerror,':-maxiter'=maxiter) elif noerror then return 'undefined' else error "f(x1) and f(x2) have different signs" end end; dx:=`if`(abserror>0,abserror,(x2-x1)*max(relerror,10^(2-Digits))); x0:=0.5*(x1+x2); f0:=evalf(f(x0)); while (f0=0.) do x0:=max(0.5*(x0+x1),x0-dx); f0:=evalf(f(x0)) end; if (s*f0>0) then do if (x2-x1x), {abserror::nonnegative:=0, relerror::nonnegative:=0, maxiter::posint:=100},$) local x1,x2,x0,iter,i,lsx,lsy,tb,flag; x1,x2:=op(r); if (x1=-infinity) then x1:=-10^Digits else x1:=evalf(x1) end; if (x2= infinity) then x2:= 10^Digits else x2:=evalf(x2) end; lsx:=evalf(Grid(x1..x2,n,GridF)); lsy:=map(evalf@f,lsx); tb:=table(): for i from 2 to n do if lsy[i-1]*lsy[i]<0 then flag,x0,x2,iter:=FZERO(f,lsx[i-1],lsx[i],(lsx[i-1]+lsx[i])/2,relerror,abserror,maxiter); if (flag=1 or flag=2 or flag=3) then tb[i]:=x0 else error "Cannot find root on interval %1",lsx[i-1]..lsx[i] end end end; convert(tb,list) end: #hfl: fZero FZERO:=proc(F,B0,C0,R0,RE,AE,MAXITER) # based on FORTRAN subroutine of SLATEC package # original authors: L. F. Shampine and H. A. Watts local A,ACBS,ACMB,B,C,CMB,FA,FB,FC,FX,FZ,P,Q,R,T,TOL,Z,IC,ITER; B:=evalf(B0); C:=evalf(C0); R:=evalf(R0); TOL:=`if`( AE>0, AE, abs(B-C)*max(RE,10^(2-Digits)) ); Z:=`if`( R<=min(B,C) or R>=max(B,C), C, R ); IC:=0; FZ:=evalf(F(Z)); FC:=FZ; FB:=evalf(F(B)); if (signum(FZ)<>signum(FB)) then C:=Z elif (Z<>C) then FC:=evalf(F(C)); if (signum(FZ)<>signum(FC)) then B:=Z; FB:=FZ end end; A:=C; FA:=FC; ACBS:=abs(B-C); FX:=max(abs(FB),abs(FC)); ITER:=0; ##### MAIN LOOP do # Perform interchange if needed if (abs(FC)FX) then return 3,B,C,ITER else return 1,B,C,ITER end elif (FB=0) then return 2,B,C,ITER elif (ITER>=MAXITER) then return 5,B,C,ITER end; # Calculate new iterate implicitly as B+P/Q, where we arrange P>=0. The implicit form is used to prevent overflow P:=(B-A)*FB; Q:=FA-FB; if (P<0) then P:=-P; Q:=-Q end; # Update A and check for satisfactory reduction in the size of the bracketing interval. If not, perform bisection A:=B; FA:=FB; IC:=IC+1; if (IC>=4 and 8*ACMB>=ACBS) then # Use bisection (C+B)/2 B:=B+CMB else if (IC>=4) then IC:=0; ACBS:=ACMB end; # Test for too small a change if (P>abs(Q)*TOL) then # Root ought to be between B and (C+B)/2 if (P>=CMB*Q) then # Use bisection (C+B)/2 B:=B+CMB else # Use secant rule B:=B+P/Q end else # Increment by TOLerance. B:=B+signum(CMB)*TOL end end; # Have completed computation for new iterate B FB:=evalf(F(B)); ITER:=ITER+1; # Decide whether next step is extrapolation or interpolation if (signum(FB)=signum(FC)) then C:=A; FC:=FA end end do end: #hfl: NLPSolve2 NLPSolve2:=proc( n::posint, f::procedure, g::{0,procedure}:=0, dx::positive:=10^(-iquo(Digits,2)), {onesidegradient::boolean:=false} ) local grad; if type(g,procedure) then grad:=g elif onesidegradient then grad:=proc(x) local f0; f0:=evalf(f(x)); Vector(n,i->evalf((f(Vector(n,j->x[j]+`if`(j=i,dx,0)))-f0)/dx),datatype=float) end else grad:=x->Vector(n,i->evalf((f(Vector(n,j->x[j]+`if`(j=i,dx,0)))-f(Vector(n,j->x[j]-`if`(j=i,dx,0))))/(2*dx)),datatype=float) end; Optimization[NLPSolve]( n, f, _rest, objectivegradient=proc(V,G) local u,i; u:=grad(V); for i from 1 to n do G[i]:=u[i] end end) end: #hfl: HessianByGradient HessianByGradient:=proc( n::posint, g::procedure, x::indexable, dx::positive:=10^(-iquo(Digits,3)), {printout::boolean:=false},$) local x0,H,j,dxj,g1,g2,i; x0:=Vector(n,i->evalf(x[i]),datatype=float); H:=Matrix(n,datatype=float); for j from 1 to n do dxj:=Vector(n,datatype=float); dxj[j]:=dx; g1:=evalf(g(x0+dxj)); g2:=evalf(g(x0-dxj)); for i from 1 to n do H[i,j]:=(g1[i]-g2[i])/(2*dx) end end; if printout then printf("Frobenius asymmetry is %.2g, maximum asymmetry is %.2g\n", sqrt(add(add((H[i,j]-H[j,i])^2,j=i+1..n),i=1..n-1)), max(seq(seq(abs(H[i,j]-H[j,i]),j=i+1..n),i=1..n-1)) ) end; Matrix(n,(i,j)->(H[i,j]+H[j,i])/2,datatype=float,shape=symmetric); end: #hfl: PolynomInterp1der PolynomialInterpolation1der:=proc(xydata::list,x,k::nonnegint,ypk) local n,L,X,Y,i,j; n:=nops(xydata); L:=[seq(i,i=1..n)]; X,Y:=map2(op,1,xydata),map2(op,2,xydata); add(Y[i]*(x-X[k])/(X[i]-X[k])*mul((x-X[j])/(X[i]-X[j]),j=subsop(i=NULL,L)),i=subsop(k=NULL,L)) +(Y[k]+(ypk-Y[k]*add(1/(X[k]-X[j]),j=subsop(k=NULL,L)))*(x-X[k]))*mul((x-X[j])/(X[k]-X[j]),j=subsop(k=NULL,L)) end: #hfl: BezierControlPoints BezierControlPoints:=proc( x1::realcons,y1::realcons,p1::realcons,q1::realcons,k1::realcons, x4::realcons,y4::realcons,p4::realcons,q4::realcons,k4::realcons) local v1,v4,a,b1,b4,e,L; a:=2*(p1*q4-p4*q1); b1:=6*(q1*(x4-x1)-p1*(y4-y1)); b4:=6*(q4*(x4-x1)-p4*(y4-y1)); L:=evalf(sqrt((x4-x1)^2+(y4-y1)^2)); e:=fsolve({k1*v1^2+a*v4=-b1,k4*v4^2+a*v1=b4},{v1=L,v4=L},{v1=0..infinity,v4=0..infinity}); if (nops([e])>1) then WARNING("Multiple solution for velocities: %1. They are set to L=%2. The input was: %3", [e],L,[x1,y1,p1,q1,k1,x4,y4,p4,q4,k4]); e:={v1=L,v4=L} elif (type(e,function)) then WARNING("Cannot solve equations for velocities: %1. They are set to L=%2. The input was: %3", {k1*v1^2+a*v4=-b1,k4*v4^2+a*v1=b4},L,[x1,y1,p1,q1,k1,x4,y4,p4,q4,k4]); e:={v1=L,v4=L} end; evalf(subs(op(e),[[x1,y1],[x1+v1*p1/3,y1+v1*q1/3],[x4-v4*p4/3,y4-v4*q4/3],[x4,y4]])) end: #hfl: BezierNodes BezierNodes:=proc(Points::list, B::list:=[0,0]) local Nodes,n,i,i1,i2,calc,ypv,kv; # Subroutine calc:=proc(k::integer:=0) local x,dx,y,yp,c2,cosa,e; x:=Points[i][1]; e:=`if`(k=0, CurveFitting[PolynomialInterpolation](map(v->[v[1]-x,v[2]],Points[max(i-2,1)..min(i+2,n)]),dx), PolynomialInterpolation1der(map(v->[v[1]-x,v[2]],Points[max(i-2,1)..min(i+2,n)]),dx,k,ypv)); y,yp,c2:=seq(coeff(e,dx,i),i=0..2); cosa:=1/sqrt(1+yp^2); Nodes[i]:=evalf([x,y,cosa,yp*cosa,2*c2*cosa^3]) end; # Begin program n:=nops(Points); Nodes:=Array(1..n); for i from 1 to n do calc() end; i1,i2:=1,n; # Left boundary if (B[1] in [0,1,2]) then i1:=1+B[1] elif (B[1]=3) then ypv:=args[3][1]; for i from 1 to 2 do calc(1) end elif (B[1]=4) then ypv,kv:=args[3][1],args[4][1]; Nodes[1]:=evalf([op(Points[1]),1/sqrt(1+ypv^2),ypv/sqrt(1+ypv^2),kv]); i:=2; calc(1) else error "Unrecognized left boundary condition:",B[1] end; # Right boundary if (B[2] in [0,1,2]) then i2:=n-B[2] elif (B[2]=3) then ypv:=args[3][2]; for i from n-1 to n do calc(n) end elif (B[2]=4) then ypv,kv:=args[3][2],args[4][2]; Nodes[n]:=evalf([op(Points[n]),1/sqrt(1+ypv^2),ypv/sqrt(1+ypv^2),kv]); i:=n-1; calc(n) else error "Unrecognized right boundary condition:",B[2] end; [seq(Nodes[i],i=i1..i2)] end: #hfl: ExtrapolateMPE ExtrapolateMPE:=proc(X::Matrix,$) local n,m,dX,A,i,b; n,m:=op(1,X); m:=m-1; dX:=Matrix(n,m,(i,j)->X[i,j+1]-X[i,j],datatype=float): A:=LinearAlgebra[Transpose](dX).dX; for i from 1 to m do A[m,i]:=1 end; b:=Vector(m,datatype=float): b[m]:=1: X[..,..m].LinearAlgebra[LinearSolve](A,b) end: ################################################################################ #cat: Differential equations #hfl: devpsolve devpsolve:=proc( x1::realcons, x2::realcons, BC, U::procedure(anything), N::posint:=100, { numev::posint:=iquo(N+1,2), evonly::boolean:=false },$) local X,dx,dt,H,i,ev,evc0,evc,p,c,S; X:=Vector(N,i->x1+(x2-x1)*(i-1)/(N-1)); dx:=evalf((x2-x1)/(N-1)); if type(BC,complexcons) then dt:=`if`(BC=1 or BC=-1,float,complex); H:=Matrix(N-1,shape=hermitian); for i from 1 to N-1 do H[i,i]:=evalf(2+U(X[i])*dx^2) end; for i from 1 to N-2 do H[i,i+1]:=-1. end; H[N-1,1]:=-evalf(BC); if evonly then return (1/dx^2)*(LinearAlgebra[Eigenvalues](H))[1..numev] end; ev,evc0:=LinearAlgebra[Eigenvectors](H); evc:=Matrix(N,numev,datatype=dt); evc[1..-2,1..-1]:=evc0[1..-1,1..numev]; evc[N,1..-1]:=evalf(BC).evc0[1,1..numev] elif type(BC,[anything,anything]) then H:=Matrix(N-2,shape=symmetric,datatype=float); for i from 1 to N-2 do H[i,i]:=evalf(2+U(X[i+1])*dx^2) end; for i from 1 to N-3 do H[i,i+1]:=-1. end; for p from 1 to 2 do if (BC[p]=1) then c[p]:=0 elif (BC[p]=2) then c[p]:=1; elif (type(BC[p],list) and type(BC[p][1],positive)) then c[p]:=2/(1+BC[1][1]*dx)-1 else error "Unrecognized %2 boundary condition: %1",BC[p],`if`(p=1,"left","right") end end; H[ 1, 1]:=H[ 1, 1]-c[1]; H[-1,-1]:=H[-1,-1]-c[2]; if evonly then return (1/dx^2)*(LinearAlgebra[Eigenvalues](H))[1..numev] end; ev,evc0:=LinearAlgebra[Eigenvectors](H); evc:=Matrix(N,numev,datatype=float); evc[2..-2,1..-1]:=evc0[1..-1,1..numev]; evc[1,1..-1]:=c[1].evc0[ 1,1..numev]; evc[N,1..-1]:=c[2].evc0[-1,1..numev] else error "Unrecognized boundary condition: %1",BC end; S:=LinearAlgebra[DiagonalMatrix]([seq(`if`(evc[2,i]>=0,1,-1),i=1..numev)]); (1/dx^2)*ev[1..numev],(1/sqrt(dx))*evc.S,evalf(X) end: ################################################################################ #cat: Linear algebra #hfl: CheckEigenvector CheckEigenvector:=proc(M::Matrix,V::Vector,maxdev::numeric:=0,simplifyf::procedure:=simplify,{trueorfalse::boolean:=false},$) local v,R; v:=simplifyf(LinearAlgebra[HermitianTranspose](V).M.V/simplifyf(add(abs(v)^2,v=V))); R:=map(simplifyf,simplify((M-v).V)); if `if`(maxdev=0, verify(R,LinearAlgebra[ZeroVector](op(1,V)),Vector), evalb(LinearAlgebra[Norm](R,2)simplify(u/v),V); break end end end; try u:=ilcm(seq(`if`(type(v,`+`),denom(op(1,v)),denom(v)),v=V)) catch: u:=1 end; if (type(u,realcons) and u<>1) then LinearAlgebra[VectorScalarMultiply](V,u,'inplace') end; try u:=igcd(seq(subs(sbs1,v),v=V)) catch: u:=1 end; if (type(u,realcons) and u<>1) then LinearAlgebra[VectorScalarMultiply](V,1/u,'inplace') end elif (_nrest>0) then LinearAlgebra[Normalize](V,_rest,'inplace') end; if positive then w:=evalf(10^digits); LinearAlgebra[VectorScalarMultiply](V,signum(V[op(MaxIdx(evalf(V),v->round(abs(v)*w)))]),'inplace') end; V end: #hfl: NormalizeMatrix NormalizeMatrix:=proc(M::Matrix,{positive::boolean:=false,int::boolean:=false,rad::boolean:=false,digits::posint:=Digits-2}) local j; for j from 1 to Dim2(M)[2] do M[..,j]:=NormalizeVector(M[..,j],_passed[2..]) end; M end: #hfl: RealJordanBlockMatrix RealJordanBlockMatrix:=proc(lsJB::list,$) local ls,JB,ev,m,b1,b2,i; ls:=[]: for JB in lsJB do ev,m:=op(JB); if Im(ev)=0 then ls:=[op(ls),LinearAlgebra[BandMatrix]([[ev$m],[1$(m-1)]],0)] else b1:=[ Im(ev),seq(op([0, Im(ev)]),i=1..m-1)]; b2:=[-Im(ev),seq(op([0,-Im(ev)]),i=1..m-1)]; ls:=[op(ls),LinearAlgebra[BandMatrix]( [b1,[Re(ev)$(2*m)],b2,[1$(2*m-2)]],1)] end end; LinearAlgebra[DiagonalMatrix](ls) end: #hfl: MatrixFunction2 MatrixFunction2:=proc(M::Matrix,f::procedure,simplifyf::procedure:=(x->x),$) local ev,evc; ev,evc:=LinearAlgebra[Eigenvectors](M); map(simplifyf,evc.LinearAlgebra[DiagonalMatrix](map(f,ev)).evc^(-1)) end: #hfl: SparseMM SparseMM:=proc(M1::Matrix,M2::Matrix,$) local T1,K,T2,js,M,i,j,k,v; T1:=table(op(2,M1)); K:=map(v->map2(op,2,v),Classify2([indices(T1)],[1])); T2:=table(op(2,M2)); js:=map2(op,2,[indices(T2)]); M:=Matrix(Dim2(M1)[1],Dim2(M2)[2]); for i in indices(K,'nolist') do for j in js do M[i,j]:=add(`if`(assigned('T2[k,j]'),T1[i,k]*T2[k,j],0),k=K[i]) end end; M end: #hfl: SparseMM SparseMHM:=proc(M::Matrix,H::{Matrix,undefined}:=undefined,{conjugate::boolean:=false},$) local ff,T,J,js,S,Sshape,pi,i,j,k,l,v; if conjugate then ff:=:-conjugate else ff:=v->v end; T:=table(op(2,M)); J:=map(v->map2(op,1,v),Classify2([indices(T)],[2])); js:=sort([indices(J,'nolist')]); S:=Matrix(Dim2(M)[2]); Sshape:=`if`(type(H,Matrix) and MatrixOptions(H,shape)<>[symmetric] and MatrixOptions(H,shape)<>[hermitian],NULL,`if`(conjugate,hermitian,symmetric)); if (Sshape=NULL) then for i in js do for j in js do S[i,j]:=add(add(ff(T[k,i])*H[k,l]*T[l,j],k=J[i]),l=J[j]) end end else if type(H,Matrix) then for pi from 1 to nops(js) do i:=js[pi]; for j in js[pi..] do S[i,j]:=add(add(ff(T[k,i])*H[k,l]*T[l,j],k=J[i]),l=J[j]) end end else for pi from 1 to nops(js) do i:=js[pi]; for j in js[pi..] do S[i,j]:=add(`if`(assigned('T[k,j]'),ff(T[k,i])*T[k,j],0),k=J[i]) end end end; for pi from 1 to nops(js) do i:=js[pi]; for j in js[..pi-1] do S[i,j]:=ff(S[j,i]) end end end; S end: #hfl: SparseMM SparseMV:=proc(M::Matrix,V::Vector,$) local T,J,U,i,j,v; T:=table(op(2,M)); J:=map(v->map2(op,2,v),Classify2([indices(T)],[1])); U:=Vector(Dim2(M)[1]); for i in indices(J,'nolist') do U[i]:=add(T[i,j]*V[j],j=J[i]) end; U end: #hfl: SparseMM SparseMO:=proc(MM::Matrix,norm::{"","2","int"},sf::procedure:=(x->x),{sortbysize::boolean:=false},$) local ls,tb,j0,J,P,L,U,m,p,U1,U2,J2,i,j,v; ls:=LinearAlgebra[StronglyConnectedBlocks](MM,'outputform'='rowlist'); if sortbysize then ls:=Sort(ls,nops) end; tb:=table(); j0:=0; for J in ls do if (nops(J)=1) then U:=MM[J[1],J[1]]; if (U<>0) then tb[j0]:=[(J[1],j0+1)=`if`(norm="2",1/sqrt(U),1)]; j0:=j0+1 end else P,L,U:=LinearAlgebra[LUDecomposition](Matrix(MM[J,J])); m:=max(seq(lhs(v)[1],v=op(2,U))); #Rank p:=map2(op,1,Sort([seq([lhs(v)],v=op(2,P))],[2]))[..m]; #permutation U1:=U[..m,p]; if (norm="2") then for i from 1 to m do LinearAlgebra[RowOperation](U1,i,1/sqrt(U1[i,i]),'inplace') end end; U2:=map(sf,U1^(-1)); if (norm="int") then U2:=NormalizeMatrix(U2,'int') end; J2:=J[p]; # M2=M[..,J2].U2, HermitianTranspose(U2).MM[J2,J2].U2=`if`(norm="2",1,diagonalmatrix) tb[j0]:=[seq(seq(`if`(U2[i,j]=0,NULL,(J2[i],j0+j)=U2[i,j]),j=i..m),i=1..m)]; j0:=j0+m end end; Matrix(Dim2(MM)[1],j0,map(op,{entries(tb,'nolist')})) end: #hfl: SparseMM SparseMN:=proc(M::Matrix) local T,J,MN,K,V,i,j,k,v; T:=table(op(2,M)); J:=map(v->map2(op,1,v),Classify2([indices(T)],[2])); MN:=Matrix(Dim2(M)); for j in indices(J,'nolist') do K:=J[j]; V:=NormalizeVector(Vector([seq(M[i,j],i=K)]),_rest); for k from 1 to nops(K) do MN[K[k],j]:=V[k] end end; MN end: #hfl: SparseMM SparseMap:=proc(f,A) local v,u; u:=_rest; op(0,A)(op(1,A),map(v->lhs(v)=f(rhs(v),u),op(2,A)),op(3,A)) end: #hfl: SuperimposeM SuperimposeM:=proc( M1::Matrix, M2::Matrix, out::string:="m", {printout::boolean:=false},$) local n,m,U1,S1,Vt1,R,f,d,v,ijs,e,c; n,m:=Dim2(M1); if (Dim2(M1)<>Dim2(M2)) then error("Different matrix sizes: %1<>%2",Dim2(M1)<>Dim2(M2)) end; U1,S1,Vt1:=LinearAlgebra[SingularValues](M1.LinearAlgebra[Transpose](M2),'output'=['U','S','Vt']); R:=U1.Vt1; f:=R.M2-M1; d:=sqrt(add(v^2,v=f)/m); # Frobenius_norm/sqrt(m) ijs:=SortIdx(convert(f,array),v->-abs(v)); ijs:=ijs[..min(5,nops(ijs))]; e:=f[op(ijs[1])]; if printout then printf("RMSD = %.3g\n",d); printf("Deviations:%{c,}s\n",Vector(map(ij->sprintf(" [%d,%d]=%.2g",op(ij),f[op(ij)]),ijs))); printf("Rotation: det=%d, angle=%.3g, axis=<%{c,}.3g>\n",RotationParam(R,'degrees')) end; seq(`if`(c="m",R.M2,`if`(c="d",d,`if`(c="e",e,`if`(c="f",f,`if`(c="r",R,NULL))))),c=out) end: #hfl: LocalizeVectors4 LocalizeVectors4:=proc(M0::Matrix,tol::numeric,{method::{"sequential","parallel"}:="sequential",printout::{boolean,posint}:=false},$) local prn,freq,a0,d,n,m,M,S4,S3,A,i,j,k,l,a,c,s,pair, o,p; if type(printout,posint) then prn,freq:=true,printout else prn,freq:=printout,10 end; a0:=1; d:=1+max(0,floor(-log10(tol))); n,m:=op(1,M0); if (m=1) then return M0 end; M:=copy(M0); if prn then printf("%-*s localization\n",d+2,"angle") end; if (method="sequential") then S4:=Matrix(m,(i,j)->evalhf(add(M[p,i]^2*M[p,j]^2,p=1..n)),shape=symmetric,datatype=float); S3:=Matrix(m,(i,j)->evalhf(add(M[p,i]*M[p,j]*(M[p,i]^2-M[p,j]^2),p=1..n)),shape=antisymmetric,datatype=float); A:=Matrix(m,(i,j)->`if`(iabs(v))); a:=A[i,j]; if (prn and abs(a)|<-s,c>>; for k in [i,j] do for l from 1 to m do S4[k,l]:=evalhf(add(M[p,k]^2*M[p,l]^2,p=1..n)); S3[k,l]:=evalhf(add(M[p,k]*M[p,l]*(M[p,k]^2-M[p,l]^2),p=1..n)) end end; for k in [i,j] do for l from 1 to m do A[op(sort([k,l]))]:=evalhf(signum(l-k)*arctan(4*S3[k,l],(S4[k,k]+S4[l,l]-6*S4[k,l]))/4) end end end else do S4:=Matrix(m,(i,j)->add(M[p,i]^2*M[p,j]^2,p=1..n),shape=symmetric,datatype=float); A:=table([]): for i from 1 to m do for j from i+1 to m do A[i,j]:=arctan(4*add(M[p,i]*M[p,j]*(M[p,i]^2-M[p,j]^2),p=1..n),S4[i,i]+S4[j,j]-6*S4[i,j])/4 end end; a:=abs(MaxVal(A,v->abs(v))); if (prn and abs(a)abs(v))); c,s:=cos(A[i,j]),sin(A[i,j]); M[..,[i,j]]:=M[..,[i,j]].<|<-s,c>>; for k in [i,j] do for l from 1 to m do unassign('A[op(sort([k,l]))]') end end end end end; M end: ################################################################################ #cat: Linear transformations #hfl: RotationM RotationM:=proc(alpha)::Matrix; local n,L,R,theta,phi; if (nargs=1) then R:=<,> else if (nargs=2 and type(args[2],indexable)) then n:=LinearAlgebra[Normalize](convert(args[2],Vector),2,conjugate=false) elif (nargs=3) then if (type(args[2],indexable) and args[3]='normalized') then n:=args[2] else theta,phi:=args[2],args[3]; n:= end else error "Unrecognized arguments: %1",[args] end; L:=Matrix(3,{(3,2)=n[1],(1,3)=n[2],(2,1)=n[3]},shape=antisymmetric); R:=1+sin(alpha)*L+(1-cos(alpha))*L^2 end; `if`(hastype([args],'float'),evalf(R),R) end: #hfl: ReflectionM ReflectionM:=proc()::Matrix; local phi,R; if (nargs=1 and not(type(args[1],indexable))) then phi:=args[1]; R:=<<-cos(2*phi)|-sin(2*phi)>,<-sin(2*phi)|cos(2*phi)>>; `if`(type(phi,'float'),evalf(R),R) else -RotationM(Pi,args) end end: #hfl: RotationParam RotationParam:=proc( R::Matrix(3,3), P::Vector(3):=<0,0,1>, T::name:='undefined', { method::{"Eigenvectors","NullSpace"}:=`if`(hastype(R,'float'),"Eigenvectors","NullSpace"), degrees::boolean:=false, digits::posint:=2 },$) local det,ev,evc,i,n,alpha,e,e1,e2,T1,v; det:=simplify(LinearAlgebra[Determinant](R)); if (evalb(abs(abs(det)-1)>10^(digits-Digits))=true) then error("Impossible determinant: %1",det) else det:=signum(det) end; if (method="Eigenvectors") then ev,evc:=LinearAlgebra[Eigenvectors](det*R); i:=op(MinIdx(ev,v->abs(v-1))); n:=map(Re,evc[..,i]); alpha:=argument(simplify(ev[(i mod 3)+1])) else e:=LinearAlgebra[NullSpace](R-det); if (nops(e)<>1) then if (add(abs(v),v=R-det)=0) then return det,0,P else error "Impossible nullspace",e end end; n:=op(e); e:=simplify((det*LinearAlgebra[Trace](R)-1)/2); if (evalb(abs(e)>1+1e-4)=true) then error "Impossible cosine",e end; alpha:=arccos(e) end; e1:=LinearAlgebra[MatrixNorm](simplify(R-det*RotationM( alpha,n)),'Frobenius'); e2:=LinearAlgebra[MatrixNorm](simplify(R-det*RotationM(-alpha,n)),'Frobenius'); if (evalf(e1)>evalf(e2)) then alpha:=-alpha end; if (min(evalf(e1),evalf(e2))>10^(digits-Digits)) then error("Failed test: %1",min(evalf(e1),evalf(e2))) end; if (LinearAlgebra[DotProduct](n,P)<0) then n,alpha:=-n,-alpha end; if not(T='undefined') then T1:=LinearAlgebra[IsSimilar](RotationM(alpha,0,0),R,output='C'); if hastype(R,'float') then T1:=evalf(T1) end; assign(T,T1) end; if degrees then alpha:=180/Pi*alpha; if hastype(R,'float') then alpha:=evalf(alpha) end end; det,alpha,n end: #hfl: RotationParam RotationParamEulerX:=proc( R::Matrix(3,3), det::{1,-1}:=signum(LinearAlgebra[Determinant](R)), $) `if`(has(R,float),evalf,simplify)([arctan(det*R[1,3],-det*R[2,3]),arccos(det*R[3,3]),arctan(det*R[3,1],det*R[3,2])]) end: #hfl: RotationParam RotationParamEulerY:=proc( R::Matrix(3,3), det::{1,-1}:=signum(LinearAlgebra[Determinant](R)), $) `if`(has(R,float),evalf,simplify)([arctan(det*R[2,3],det*R[1,3]),arccos(det*R[3,3]),arctan(det*R[3,2],-det*R[3,1])]) end: unassign('phi,theta,psi,alpha,beta,gamma'); #hfl: RotationMEulerX RotationMEulerX:=unapply(RotationM(phi,<0,0,1>).RotationM(theta,<1,0,0>).RotationM(psi,<0,0,1>),phi,theta,psi): #hfl: RotationMEulerX RotationMEulerY:=unapply(RotationM(alpha,<0,0,1>).RotationM(beta,<0,1,0>).RotationM(gamma,<0,0,1>),alpha,beta,gamma): #hfl: EulerY2Rot EulerY2Rot:=proc(alpha,beta,gamma,$) simplify(2*arccos(cos(beta/2)*cos((alpha+gamma)/2))), `if`(beta=0,(1-signum(sin((alpha+gamma)/2)))/2*Pi,simplify(arccot(cot(beta/2)*sin((alpha+gamma)/2)))), simplify(Reduce2P((alpha-gamma+Pi)/2,2*Pi)) end: #hfl: EulerY2Rot Rot2EulerY:=proc(alpha,e) local theta,phi,s; if (_npassed=2) then theta:=simplify(arccos(e[3]/sqrt(e[1]^2+e[2]^2+e[3]^2))); phi:=simplify(arctan(e[2],e[1])) elif (_npassed=3) then theta:=e; phi:=_passed[3] else error("Unrecognized input: expected (alpha,e) or (apha,theta,phi), but received %1",_passed) end; s:=simplify(arctan(cos(theta)*sin(alpha/2),cos(alpha/2))); simplify(Reduce2P(s+phi-Pi/2,2*Pi)), simplify(2*arcsin(sin(theta)*sin(alpha/2))), simplify(Reduce2P(s-phi+Pi/2,2*Pi)) end: #hfl: RotationMp2q RotationMp2q:=proc(p::Vector,q::Vector,{normalized::boolean:=false},$) local npnq,L; npnq:=`if`(normalized,1,LinearAlgebra[Norm](p,2)*LinearAlgebra[Norm](q,2)); L:=(q.LinearAlgebra[Transpose](p)-p.LinearAlgebra[Transpose](q))/npnq; 1+L+L^2/(1+LinearAlgebra[Transpose](p).q/npnq) end: #hfl: RotationM2zx RotationM2zx:=proc(p::indexable,q::indexable,axes::{12,13,23,21,32,31}:=31,det::{0,1,-1}:=1,$) local i,j,k,e; i,j:=seq(parse(v),v=convert(axes,string)); k:=6-i-j; e:=Vector(3); e[i]:=LinearAlgebra[Normalize](convert(p,Vector),2); e[k]:=LinearAlgebra[Normalize](LinearAlgebra[CrossProduct](e[i],convert(q,Vector)),2); e[j]:=LinearAlgebra[CrossProduct](e[k],e[i]); if (det=1 and member(axes,[13,32,21]) or det=-1 and member(axes,[12,23,31])) then e[k]:=-e[k] end; Matrix(3,(i,j)->e[i][j]); end: #hfl: TransSpherCoo TransSpherCoo:=(t,p,a,b,g)->( arccos(cos(t)*cos(b)+sin(t)*sin(b)*cos(p-a)), arctan(sin(t)*sin(p-a),sin(t)*cos(b)*cos(p-a)-cos(t)*sin(b))-g ): #hfl: TransSpherCoo TransSpherCooI:=(t,p,a,b,g)->( arccos(cos(t)*cos(b)-sin(t)*sin(b)*cos(p+g)), arctan(sin(t)*sin(p+g),sin(t)*cos(b)*cos(p+g)+cos(t)*sin(b))+a ): #hfl: TransTensor TransTensor:=proc(A::indexable,T::Matrix) local dims,B,V,S,k,i,j; B:=Array(A); dims:=[ArrayDims(B)]; for k from 1 to nops(dims) do V:=Vector(dims[k],i->B[op(subsop(k=i,dims))]); for i in [$dims[k]] do S:=0; for j in [$dims[k]] do S:=S+T[i,j]*V[j] end; B[op(subsop(k=i,dims))]:=S end end; B end: #hfl: ExtendedTransM ExtendedTransM:=proc(T::Matrix) local ET,i,j; ET:=Matrix(6); for i from 1 to 3 do for j from 1 to 3 do ET[i,j]:=T[i,j]^2; ET[i,j+3]:=2*T[i,(j mod 3)+1]*T[i,(j+1 mod 3)+1]; ET[i+3,j]:=T[(i mod 3)+1,j]*T[(i+1 mod 3)+1,j]; ET[i+3,j+3]:=T[(i mod 3)+1,(j mod 3)+1]*T[(i+1 mod 3)+1,(j+1 mod 3)+1] +T[(i mod 3)+1,(j+1 mod 3)+1]*T[(i+1 mod 3)+1,(j mod 3)+1] end end; ET end: ################################################################################ #cat: #hfl: Ylm #toc: Spherical harmonics #BasicTools/Phim,NPhim2,Pnma,NPnma2,Ylm,Yxyz,YxyzM,YxyzR,YxyzL,Wignerd,WignerD,Wigner3j,Wigner6j,MatElemYlm,MatElemY Phim:=proc(m::integer,phi,{complex::boolean:=complexTrigBasis},$) if complex then exp(I*m*phi) else if (m=0) then 1 elif (m>0) then cos(m*phi) else sin(abs(m)*phi) end end end: NPhim2:=proc(m::integer,{complex::boolean:=complexTrigBasis},$) `if`(complex or m=0,2*Pi,Pi) end: Pnma:=proc(n::nonnegint,m::nonnegint,arg) local z,v,theta,a; if (m>n) then error "m must be not greater than n",m,n end; if (nargs=3) then theta:=args[3]; `if`(m=0, orthopoly['P'](n,cos(theta)), simplify(subs(z=cos(theta),diff(orthopoly['P'](n,z),z$m))*v)/v*sin(theta)^m) else a,theta:=args[3],args[4]; if (evalb(a<=-1/2)=true or evalb(evalf(a)<=-0.5)=true) then error "a must be greater than -1/2",a end; `if`(m=0, orthopoly['G'](n,a+1/2,cos(theta)), simplify(subs(z=cos(theta),diff(orthopoly['G'](n,a+1/2,z),z$m))*v)/v*sin(theta)^m) end end: NPnma2:=(n,m,a)->2*Pi/4^a/GAMMA(a+1/2)^2*GAMMA(2*a+1+n+m)/(2*a+1+2*n)/(n-m)!: Ylm:=proc(l0,m::integer,t0,phi, convention::string:=conventionY, {normalized::boolean:=normalizedY},$) local l,t,e,k,A; if type(l0,list) then l:=ListTools[Reverse](l0) else l:=[l0] end; if type(t0,list) then t:=ListTools[Reverse](t0) else t:=[t0] end; if (nops(l)<>nops(t)) then error "Inconsistent l and theta",l0,t0 end; for k from 1 to nops(l) do if not(type(l[k],nonnegint)) then error "l must be nonnegative integer",l0 end end; for k from 2 to nops(l) do if (l[k]0 and l[1]0,(-1)^m,1)*Phim(m,phi,'complex') elif (convention="complex") then A:=Phim(m,phi,'complex') elif (convention="real") then A:=Phim(m,phi,'complex'=false) else error "Convention can be one of [standard,complex,real] but received %1",convention end; e:=mul(factor(Pnma(l[k],`if`(k=1,abs(m),l[k-1]),(k-1)/2,t[k])),k=1..nops(l))*A; `if`(normalized, e*simplify(1/sqrt(mul(NPnma2(l[k],`if`(k=1,abs(m),l[k-1]),(k-1)/2),k=1..nops(l))*NPhim2(m,'complex'=evalb(convention<>"real")))), e) end: Yxyz:=proc(l::nonnegint,m::{integer,string},x,y,z,$) local X,Y,Z,q,u,C,d,c,n,r,v, m2,i; if type(m,integer) then if (abs(m)>l) then error("%1=|m|>l=%2",abs(m),l) end; u:=`if`(m=0,orthopoly['P'](l,Z),diff(orthopoly['P'](l,Z),Z$abs(m))); C:=[coeffs(u,Z)]; d:=max(map(denom,C)); c:=igcd(op(d*C))/d; n:=l-abs(m); u:=algsubs(r^2=X^2+Y^2+Z^2,expand(r^n*subs(Z=Z/r,u)/c/`if`(type(n,odd),Z,1))); v:=expand(`if`(m<0,((X+I*Y)^abs(m)-(X-I*Y)^abs(m))/2/I,((X+I*Y)^m+(X-I*Y)^m)/2)/`if`(type(m,odd),`if`(m<0,Y,X),1)); `if`(type(m,odd),`if`(m<0,y,x),1)*`if`(type(n,odd),z,1)*subs(X=x,Y=y,Z=z,v*u)*c/sqrt(NPnma2(l,abs(m),0))*`if`(m=0,1/sqrt(2),1)/sqrt(Pi) elif (m="trig") then [seq(Yxyz(l,m2,x,y,z),m2=[0,seq(op([i,-i]),i=1..l)])] elif (m="cubic") then [seq(op(Yxyz(l,m2,x,y,z)),m2=YxyzR(l))] elif (l=0) then if (m="A1") then [ 1/sqrt(4*Pi) ] else error("For l=0 m can be one of [A1], but received %1",m) end elif (l=1) then if (m="F1") then [ x*sqrt(3/4/Pi), y*sqrt(3/4/Pi), z*sqrt(3/4/Pi) ] else error("For l=1 m can be one of [F1], but received %1",m) end elif (l=2) then if (m="E" ) then [ (x^2-y^2)*sqrt(15/16/Pi), (2*z^2-x^2-y^2)*sqrt(5/16/Pi) ] elif (m="F2") then [ y*z*sqrt(15/4/Pi), z*x*sqrt(15/4/Pi), x*y*sqrt(15/4/Pi) ] else error("For l=2 m can be one of [E,F2], but received %1",m) end elif (l=3) then if (m="A2") then [ x*y*z*sqrt(105/4/Pi) ] elif (m="F1") then [ x*(2*x^2-3*y^2-3*z^2)*sqrt(7/16/Pi), y*(2*y^2-3*z^2-3*x^2)*sqrt(7/16/Pi), z*(2*z^2-3*x^2-3*y^2)*sqrt(7/16/Pi) ] elif (m="F2") then [ x*(y^2-z^2)*sqrt(105/16/Pi), y*(z^2-x^2)*sqrt(105/16/Pi), z*(x^2-y^2)*sqrt(105/16/Pi) ] else error("For l=3 m can be one of [A2,F1,F2], but received %1",m) end elif (l=4) then if (m="A1") then [ (x^4+y^4+z^4-3*x^2*y^2-3*y^2*z^2-3*z^2*x^2)*sqrt(21/16/Pi) ] elif (m="E" ) then [ (x^2-y^2)*(6*z^2-x^2-y^2)*sqrt(45/64/Pi), (x^4+y^4-2*z^4-12*x^2*y^2+6*x^2*z^2+6*y^2*z^2)*sqrt(15/64/Pi) ] elif (m="F1") then [ y*z*(y^2-z^2)*sqrt(315/16/Pi), z*x*(z^2-x^2)*sqrt(315/16/Pi), x*y*(x^2-y^2)*sqrt(315/16/Pi) ] elif (m="F2") then [ y*z*(6*x^2-y^2-z^2)*sqrt(45/16/Pi), z*x*(6*y^2-z^2-x^2)*sqrt(45/16/Pi), x*y*(6*z^2-x^2-y^2)*sqrt(45/16/Pi) ] else error("For l=4 m can be one of [A1,E,F1,F2], but received %1",m) end else if member(m,["A1","A2","E","F1","F2"]) then error("Max supported l is 4, but received %1",l) else error("m can be one of [trig,cubic,A1,A2,E,F1,F2], but received %1",m) end end end: YxyzM:=proc(l::nonnegint,m::{integer,string},convention::string:=conventionY,$) option remember; local e,t,p,b,M; e:=Yxyz(l,m,sin(t)*cos(p),sin(t)*sin(p),cos(t)); b:=type(e,list); if not(b) then e:=[e] end; M:=Matrix(nops(e),2*l+1,(i,j)->combine(simplify(int(int(Ylm(l,iidI(j),t,`if`(convention="real",p,-p),convention,'normalized')*e[i],p=0..2*Pi)*sin(t),t=0..Pi)),sqrt)); `if`(b,M,M[1,..]) end: YxyzR:=proc(l::nonnegint,$) if (l=0) then ["A1"] elif (l=1) then ["F1"] elif (l=2) then ["E","F2"] elif (l=3) then ["A2","F1","F2"] elif (l=4) then ["A1","E","F1","F2"] else error("Max supported l is 4, but received %1",l) end end: YxyzL:=proc(l::nonnegint,$) if (l=0) then ["S"] elif (l=1) then ["X","Y","Z"] elif (l=2) then ["X2","Z2","YZ","ZX","XY"] elif (l=3) then ["XYZ","X3","Y3","Z3","XY2","YZ2","ZX2"] elif (l=4) then ["S4","X2Z2","Z4","YZY2","ZXZ2","XYX2","YZX2","ZXY2","XYZ2"] else error("Max supported l is 4, but received %1",l) end end: Wignerd:=proc(l,m1,m2,c,s) local k; sqrt((l+m1)!*(l+m2)!*(l-m1)!*(l-m2)!)*add((-1)^(l-k-m2)/k!/(k+m1+m2)!/(l-k-m1)!/(l-k-m2)!*c^(2*k+m1+m2)*s^(2*l-2*k-m1-m2),k=max(0,-m1-m2)..l-max(m1,m2)) end: WignerD:=(l,m1,m2,alpha,beta,gamma)->exp(I*m1*alpha+I*m2*gamma)*Wignerd(l,m1,m2,cos(beta/2),sin(beta/2)): Wigner3j:=proc(j1,j2,j3,m1,m2,m3) # Adapted from the code of Dennis Isbister, dji@adfa.edu.au local s,ss,f,edmonds; if type([args],list(rational)) then if abs(m1)>j1 or abs(m2)>j2 or abs(m3)>j3 then RETURN (0) fi; if m1+m2+m3 <> 0 then RETURN (0) fi; if j3 < abs(j1-j2) or j3 > j1+j2 then RETURN (0) fi; ss:=0: for s from 0 to j3+m3 do if j1-m1-s < 0 or j3+m3-s < 0 or j2-j3+m1+s < 0 then f:=0 else f:=(j1+m1+s)!*(j2+j3-m1-s)!/(s!*(j1-m1-s)!*(j2-j3+m1+s)!*(j3+m3-s)!): fi; ss:=ss+(-1)**s*f; od; f:=(j1+j2-j3)!*(j1-m1)!*(j2-m2)!*(j3-m3)!*(j3+m3)! /((j1+j2+j3+1)!*(j1-j2+j3)!*(-j1+j2+j3)!*(j1+m1)!*(j2+m2)!): edmonds:=(-1)**(-2*j1-m1-j2-m3)*ss*sqrt(f); simplify(edmonds) else 'Wigner3j'(args) fi; end: Wigner6j:= proc(j1,j2,j3,l1,l2,l3) # Adapted from the code of Dennis Isbister, dji@adfa.edu.au local delta,s,smin,smax,ss,pref,f; delta:= proc(j1,j2,j3) local p; if abs(j1 - j2) <= j3 and j3 <= j1 + j2 then p:=(j1+j2-j3)!*(j1-j2+j3)!*(-j1+j2+j3)!/(j1+j2+j3+1)!; else RETURN(0) fi; RETURN(sqrt(p)) end: if type([args],list(rational)) then pref:=delta(j1,j2,j3)*delta(j1,l2,l3) *delta(l1,j2,l3)*delta(l1,l2,j3): if (pref = 0) then RETURN(0) fi; smin:= max(j1+j2+j3,j1+l2+l3,l1+j2+l3,l1+l2+j3); smax:= min(j1+j2+l1+l2,j2+j3+l2+l3,j3+j1+l3+l1); ss:=0: for s from smin to smax do if s-j1-j2-j3<0 or s-j1-l2-l3<0 or s-l1-j2-l3<0 or s-l1-l2-j3<0 or j1+j2+l1+l2-s<0 or j2+j3+l2+l3-s<0 or j3+j1+l3+l1-s<0 then f:= infinity else f:=(s-j1-j2-j3)!*(s-j1-l2-l3)!*(s-l1-j2-l3)! *(s-l1-l2-j3)!*(j1+j2+l1+l2-s)! *(j2+j3+l2+l3-s)!*(j3+j1+l3+l1-s)!: fi; ss:=ss+(-1)**s*(s+1)!/f: od; pref:=delta(j1,j2,j3)*delta(j1,l2,l3) *delta(l1,j2,l3)*delta(l1,l2,j3): simplify(pref*ss) else 'Wigner6j'(args) fi; end: MatElemYlm:=(l1,m1,l,m,l2,m2)->simplify((-1)^m1*sqrt((2*l1+1)*(2*l+1)*(2*l2+1)/4/Pi)*Wigner3j(l1,l,l2,-m1,m,m2)*Wigner3j(l1,l,l2,0,0,0)): MatElemY:=proc(l1,m1,f,l2,m2,convention::string:=conventionY,$) simplify(int(int( Ylm(l1,m1,theta,`if`(convention="real",phi,-phi),convention,'normalized') *f(theta,phi) *Ylm(l2,m2,theta,phi,convention,'normalized'),phi=0..2*Pi)*sin(theta),theta=0..Pi)) end: ################################################################################ #cat: #hfl: LerchPhi #toc: Lerch transcendent #BasicTools/LerchPhi2,polylogR,polylogR_init LerchPhi2:=proc(z,a,v,vmax::numeric:=10,$) local n; `if`(type(v,numeric) and v>vmax and z<>0 and z<>1, (-ln(z))^(a-1)/z^v*GAMMA(1-a,-v*ln(z))+add(pochhammer(a,n)/n!/v^n*polylogR(z,n),n=0..polylogR_n)/v^a, LerchPhi(z,a,v)) end: polylogR:=proc(z,n) local v,k; if (z=0) then `if`(n=0,1,0) elif (n>polylogR_n) then undefined elif (n=0) then `if`(z=1,1/2,1/ln(z)+1/(1-z)) elif (n<=10 and z<0.5 or n<=20 and z<0.1) then n!/ln(z)^(n+1)+(-1)^n*polylog(-n,z) else v:=0; for k from polylogR_N by -1 to 0 do v:=polylogR_T[n,k]+v*(1-z) end; v end end: polylogR_init:=proc(n::posint,N::posint,$) local x,e,k,n2; e:=convert(series(1/ln(1-x)+1/x,x,N+n+3),polynom); polylogR_T:=Array(0..N,0..N+n,datatype=rational); for k from 0 to N+n do polylogR_T[0,k]:=coeff(e,x,k) end; for n2 from 1 to n do for k from 0 to N+n-n2 do polylogR_T[n2,k]:=(k+1)*polylogR_T[n2-1,k+1]-k*polylogR_T[n2-1,k] end end; polylogR_n,polylogR_N:=n,N; NULL end: ################################################################################ #cat: File tools #hfl: ExpandPath ExpandPath:=proc( path::string, output::string:="apnx", { replacedirsepto::string:="/" },$) local i,p,nx,n,x,c,a,u,v; i:=max(StringTools[FirstFromRight]("/",path),StringTools[FirstFromRight]("\\",path)); p:=path[1..i]; nx:=path[i+1..-1]; i:=StringTools[FirstFromRight](".",nx); if (i=0) then n,x:=nx,"" else n,x:=nx[1..i-1],nx[i..-1] end; c:=cat(currentdir(),kernelopts(dirsep)); if (replacedirsepto<>"") then p:=StringTools[SubstituteAll](StringTools[SubstituteAll](p,"/",replacedirsepto),"\\",replacedirsepto); c:=StringTools[SubstituteAll](StringTools[SubstituteAll](c,"/",replacedirsepto),"\\",replacedirsepto) end; a:=`if`(p="/" or StringTools[Has](p,":"),"",c); seq(cat(seq(piecewise(u="a",a,u="c",c,u="p",p,u="n",n,u="x",x,""),u=v)),v=StringTools[Split](output,",")) end: #hfl: ExpandPath SimplifyPath:=proc( path::string, { replacedirsepto::string:="/" },$) local ls,exitdo,j,i,n,v; ls:=StringTools[Split](path,"/\\"); ls:=[seq(`if`(v="",NULL,v),v=ls[..-2]),ls[-1]]; do exitdo:=true; for j from nops(ls) by -1 to 1 do if (ls[j]="..") then for i from j by -1 to 1 while (ls[i]="..") do end; if (i=0) then break end; n:=min(j-i,i); ls:=[op(1..i-n,ls),op(i+1+n..-1,ls)]; exitdo:=false; break end end; if exitdo then break end end; StringTools[Join](ls,`if`(replacedirsepto="",kernelopts('dirsep'),replacedirsepto)) end: #hfl: FoldersTree FoldersTree:=proc(folder::string,depth::{infinity,nonnegint}:=infinity,$) local path; path:=cat(folder,kernelopts('dirsep')); [FileTools[Filename](folder),`if`(depth=0,NULL,op(map(FoldersTree,select(proc(v) local b; try b:=FileTools[IsDirectory](v) catch: WARNING("IsDirectory is failed on %1",v); b:=false end; b end, map2(cat,path,FileTools[ListDirectory](folder))),depth-1)))] end: #hfl: LoadTextFile LoadTextFile:=proc(f::string,keyfun::procedure:=(s->true),nlines::nonnegint:=infinity,{reload::boolean:=false},$) local fd,tb,s,i,size,v; if FileTools[IsOpen](f) then error("Cannot load open file: %1",f) end; if not(FileTools[Exists](f)) then error("File does not exists: %1",f) end; if assigned('TextFiles[f]') then if (reload or FileTools[ModificationTime](f)<>TextFilesAttr[f][1]) then UnloadTextFile(f) else return end end; fd:=fopen(f,READ,TEXT); tb:=table(); s:=readline(fd); for i from 1 to infinity while not(s=0 or keyfun(s)) do tb[i]:=s; s:=readline(fd) end; for i from i to i-1+nlines while not(s=0) do tb[i]:=s; s:=readline(fd) end; fclose(fd); size:=add(length(s),s=entries(tb,'nolist')); if (size>sTextFiles) then error("File size =%1 and is larger than sTextFiles=%2",size,sTextFiles) end; TextFiles[f]:=convert(tb,list); TextFilesAttr[f]:=[ FileTools[ModificationTime](f), size ]; # [modtime,size] TextFilesList:=[op(TextFilesList),f]; # Cleanup while (nops(TextFilesList)>nTextFiles) do UnloadTextFile(TextFilesList[1]) end; while (add(TextFilesAttr[v][2],v=TextFilesList)>sTextFiles) do for i from nops(TextFilesList)-1 by -1 to 2 while (TextFilesAttr[TextFilesList[i]][2]<=add(TextFilesAttr[v][2],v=TextFilesList[..i-1])) do end; UnloadTextFile(TextFilesList[i]) end; NULL end: #hfl: LoadTextFile UnloadTextFile:=proc(f::string,$) if assigned('TextFiles[f]') then unassign('TextFiles[f],TextFilesAttr[f]'); TextFilesList:=Substitute(TextFilesList,f) end; NULL end: #hfl: SearchFilePos SearchFilePos:=proc( f::string, key::string, keypos::range:=1..-1, number::integer:=1, { skiplines::integer:=0, linematch::boolean:=`if`(key="",true,false), reload::boolean:=false, noerror::boolean:=false },$) local j0,L,n,i1,tb,k,i,lsj,j; j0:=op(1,keypos)-1; if (j0<0) then error("No negative lower bound in keypos: %1",keypos) end; LoadTextFile(f,':-reload'=reload); L:=TextFiles[f]; n:=nops(L); i1:=`if`(skiplines<0,max(1,n+skiplines+1),skiplines+1); tb:=table(); k:=0; if linematch then if (number>0) then for i from i1 to n do if (L[i]=key) then k:=k+1; if (k=number) then return [i,1] end end end elif (number<0) then for i from n-i1+1 by -1 to 1 do if (L[i]=key) then k:=k-1; if (k=number) then return [i,1] end end end else return [seq(`if`(L[i]=key,[i,1],NULL),i=i1..n)] end else if (number>0) then for i from i1 to n do lsj:=[StringTools[SearchAll](key,L[i][keypos])]; if (lsj<>[]) then if (nops(lsj)[]) then if (nops(lsj)evalb(s="")), { skiplines::integer:=0, linematch::boolean:=`if`(key="",true,false), shift::integer:=`if`(key=undefined,0,1), nlines::{nonnegint,infinity}:=`if`(key=undefined,10^9,infinity), line::name:=undefined, reload::boolean:=false },$) local i0,p,i2,L,i; if (key=undefined) then LoadTextFile(f,':-reload'=reload); i0:=1 else p:=SearchFilePos(f,key,keypos,`if`(number=0,1,number),':-skiplines'=skiplines,':-linematch'=linematch,':-reload'=reload,'noerror'); if (p=[]) then return [] else i0:=p[1] end end; i2:=min(nops(TextFiles[f]),i0+shift+nlines-1); L:=TextFiles[f][max(1,i0+shift)..i2]; if (nlines=infinity) then for i from 1 to nops(L) do if endkeyfun(L[i]) then i2:=i2+i-nops(L); L:=L[..i-1]; break end end end; if not(line=undefined) then assign(line,i2) end; L end: #hfl: ReadValue ReadValue:=proc( f::string, key::{string,undefined}:=undefined, keypos::range:=1..-1, number::integer:=1, { skiplines::integer:=0, linematch::boolean:=`if`(key="",true,false), shift::integer:=0, startofline::boolean:=false, pos::range:=1..-1, format::string:="", line::name:=undefined, reload::boolean:=false, raiseerror::boolean:=false },$) local p,n,i,j,s,v; if (key=undefined) then LoadTextFile(f,':-reload'=reload); p:=[1,1] else p:=SearchFilePos(f,key,keypos,`if`(number=0,1,number),':-skiplines'=skiplines,':-linematch'=linematch,':-reload'=reload,'noerror'=not(raiseerror)); if (p=[]) then return end end; n:=nops(TextFiles[f]); i:=p[1]+shift; if (i<1 or i>n) then return end; if not(line=undefined) then assign(line,i) end; j:=`if`(startofline or shift<>0,1,p[2]+`if`(key=undefined,0,length(key))); s:=TextFiles[f][i][j..][pos]; if (format="") then s else v:=sscanf(s,format); if (v=0 or v=[]) then if raiseerror then error("Cannot read format=%1 in %2",format,s) else s end else op(v) end end end: #hfl: ReadPaginatedMatrix ReadPaginatedMatrix:=proc( f::string, key::{string,undefined}:=undefined, keypos::range:=1..-1, number::integer:=1, endkeyfun::procedure:=(s->evalb(s="")), { skiplines::integer:=0, linematch::boolean:=`if`(key="",true,false), shift::integer:=`if`(key=undefined,0,1), pos::range:=1..-1, format::string:="%f", sep::string:=" ", sym::boolean:=false, skip::nonnegint:=0, reload::boolean:=false },$) local p,L,k,n,V,m,j,i,M,ls; if (key=undefined) then LoadTextFile(f,':-reload'=reload); p:=[1,1] else p:=SearchFilePos(f,key,keypos,`if`(number=0,1,number),':-skiplines'=skiplines,':-linematch'=linematch,':-reload'=reload) end; L:=TextFiles[f][max(1,p[1]+shift)..]; for k from 1 to nops(L) while not(endkeyfun(L[k])) do end; n:=k-1; V:=Vector(n,i->L[i]); m:=nops(remove(`=`,StringTools[Split](V[n],sep),"")); if (m=0) then error("No entries in %1",V[n]) end; for j from m+1 by m to n do for k from k to nops(L) while endkeyfun(L[k]) do end; k:=k-1+skip; for i from `if`(sym,j,1) to n do k:=k+1; V[i]:=cat(V[i],sep,L[k][pos]) end end; M:=Matrix(n,`if`(sym,shape=symmetric,NULL),`if`(format="%f",datatype=float,`if`(format="%d",datatype=integer,NULL))); for i from 1 to n do ls:=remove(`=`,StringTools[Split](V[i],sep),""); if (nops(ls)=n or sym and nops(ls)=i) then for j from 1 to `if`(sym,i,n) do M[i,j]:=op(sscanf(ls[j],format)) end else error("Wrong number of entries in %1",V[i]) end end; M end: #hfl: WriteLines WriteLines:=proc(f::{string,integer},ls::list,{overwrite::boolean:=false,append::boolean:=false,printout::{boolean,nonnegint,[nonnegint,nonnegint]}:=false},$) local fd,s; if type(f,integer) then fd:=f elif (f="") then fd:='default' else if FileTools[Exists](f) then if (overwrite or append) then fd:=fopen(f,`if`(append,'APPEND','WRITE'),'TEXT') else error("File exists: %1",f) end else fd:=fopen(f,'WRITE','TEXT') end end; for s in ls do if type(s,string) then writeline(fd,s) end end; if (type(f,string) and f<>"") then fclose(fd); if not(printout=false) then PrintFile(f,`if`(printout=true,NULL,`if`(type(printout,posint),printout,op(printout)))) end end; NULL end: #hfl: ReplaceLines ReplaceLines:=proc(f::string,change::list,{nobackup::boolean:=false,forcebackup::boolean:=false},$) local ls; ls:=ReadLines(f); ls:=map(proc(s) local v; for v in change do if SearchText(v[1],s)>0 then return `if`(nops(v)=1,NULL,`if`(nops(v)=2,v[2],StringTools[Substitute](s,v[2],v[3]))) end end; s end,ls); FileTools[Rename](f,cat(f,".bak"),'force'=forcebackup); WriteLines(f,ls,'overwrite'); if nobackup then FileTools[Remove](cat(f,".bak")) end; NULL end: #hfl: ReadRecord ReadRecord:=proc( filename::string, records::string:="", linenum::nonnegint:=0, opt4ReadLines::list:=[], { comment::string:="#", assignment::string:="=", nlines::posint:=1, input::{"file","string","auto"}:="auto", nowarning::boolean:=false },$) local ls,tb,s,l,s1,i,out,dv,t,v; if (input="file" or (input="auto" and FileTools[Exists](filename))) then if (linenum=0) then ls:=remove(s->s[1]=comment,map(Trim,ReadLines(filename,op(opt4ReadLines)))) else ls:=ReadLines(filename,'shift'=linenum-1,':-nlines'=nlines); if (ls<>[]) then ls:=select(s->SearchText(assignment,s)>0,map(Trim,map(op@StringTools[Split],ls),[",",".",";",":","?","!"])) end end elif (input="string" or (input="auto" and SearchText(assignment,filename)>0)) then ls:=StringTools[Split](filename,"\n"); ls:=select(s->SearchText(assignment,s)>0,map(Trim,map(op@StringTools[Split],ls),[",",".",";",":","?","!"])) else if (input="auto" and not(nowarning)) then WARNING("File %1 does not exist and filename does not contain assignment symbol %2",filename,assignment) end; ls:=[] end; tb:=table(); s:=""; for l from nops(ls) by -1 to 1 do s1:=ls[l]; i:=SearchText(comment,s1); if (i>0) then s1:=TrimRight(s1[..i-1]) end; i:=SearchText(assignment,s1); if (i>0) then if (s1[i+1]=assignment) then s:=cat(s1[..i],s1[i+2..],s) else tb[TrimRight(s1[..i-1])]:=cat(TrimLeft(s1[i+1..]),s); s:="" end else s:=cat(s1,s) end end; if (s<>"") then error("Unprocessed string in %1, s=%2",filename,s) end; if (records="") then return op(tb) end; out:=[]; for s in StringTools[Split](records,",") do i:=SearchText(":=",s); if (i>0) then s,dv:=s[..i-1],s[i+2..] else s,dv:=s,"" end; i:=SearchText("::",s); if (i>0) then t:=parse(s[i+2..]); if not(type(t,type)) then error("Unrecognized type: %1",t) end; s:=s[..i-1] else s,t:=s,anything end; if assigned('tb[s]') then if (t=list(string)) then if not(tb[s][1]="[" and tb[s][-1]="]") then error("Unrecognized list(string) record: %1",tb[s]) else v:=map(Trim,StringTools[Split](tb[s][2..-2],",")) end else v:=`if`(t=string,tb[s],eval(parse(tb[s]))) end; if not(type(v,t) or t=table and type(v,list)) then error("Unrecognized value (must be of type %1): %2",t,v) end else if (dv="") then v:=`if`(t=string,"",`if`(t=list or t=list(string) or t=table,[],`if`(subtype(t,extended_numeric),0,undefined))) else v:=`if`(t=string,dv,eval(parse(dv))) end end; # out:=out,`if`(t=table,table(v),v) end; out:=[op(out),[t,v]] end; #`if`(nops([out])=1 and type(out,table),op(out),out) seq(`if`(v[1]=table,table(v[2]),v[2]),v=out) end: #hfl: ReadRecord WriteRecord:=proc(filename::string,records::list({string,[string,anything],[string,anything,string]}),{overwrite::boolean:=false,append::boolean:=false},$) local printdatatype,formatrecord; printdatatype:=v->`if`(rtable_options(v,datatype)=anything,"",sprintf(",datatype=%a",rtable_options(v,datatype))); formatrecord:=proc(v,$) local s; if type(v,string) then v elif type(v,[string,anything]) then `if`(type(v[2],list(string)), cat(v[1],"=[",StringTools[Join](v[2],","),"]"), sprintf("%s=%a",op(v))) else if type(v[2],atomic) then sprintf(cat("%s=",v[3]),v[1],v[2]) elif type(v[2],list(list)) then sprintf("%s=[%{c,}s]",v[1],Vector(nops(v[2]),i->sprintf(cat("[%{c,}",v[3][2..],"]"),Vector(nops(v[2][i]),j->v[2][i][j])))) elif type(v[2],list) then sprintf(cat("%s=[%{c,}",v[3][2..],"]"),v[1],Vector(nops(v[2]),i->v[2][i])) elif type(v[2],Vector) then sprintf(cat("%s=Vector(%d%s,[%{c,}",v[3][2..],"])"),v[1],op(1,v[2]),printdatatype(v[2]),v[2]) elif type(v[2],Matrix) then s:=StringTools[SubstituteAll](sprintf(cat("%{c,}",v[3][2..]),v[2]),"\n",",\n"); sprintf("%s=Matrix(%d,%d%s,[\n%s])",v[1],op(1,v[2]),printdatatype(v[2]),s) elif type(v[2],table) then s:=StringTools[Join](map(StringTools[Substitute],FormatTable(v[2],v[3],'align'),"=","=="),",\n"); sprintf("%s=[\n%s,\nNULL]",v[1],s) else error("Unsupported type in %1",v) end end end; WriteLines(filename,map(formatrecord,records),':-overwrite'=overwrite,':-append'=append) end: #hfl: ReadSet ReadSet:=proc(s::string,f::string,n::integer:=99,$) local ls,ls2,s2; ls:=StringTools[Split](s,","); if (n<0) then [] else if (nops(ls)=1) then ls:=ReadRecord(f,cat(s,"::list(string)"),'input'="file"); if (ls=[]) then [s] else ReadSet(StringTools[Join](ls,","),f,n-1) end else ListTools[MakeUnique]([seq(op(ReadSet(s2,f,n)),s2=ls)]) end end end: #hfl: ReadUTF8 ReadUTF8:=proc(fn::string,{unicodes::boolean:=false},$) local ls,tb,i,k; ls:=FileTools[Binary][Read](fn,integer[1],'byteorder'='native'); fclose(fn); ls:=map(v->`if`(v>0,v,v+256),ls); tb:=table(); for i from 1 to nops(ls) do if (ls[i]<128) then tb[i]:=ls[i] elif (ls[i]<192) then error("Unrecognized byte #%1: %2",i,ls[i]) elif (ls[i]<224) then tb[i]:=64*(ls[i]-192)+ls[i+1]-128; i:=i+1 elif (ls[i]<240) then tb[i]:=64^2*(ls[i]-224)+add(64^(2-k)*(ls[i+k]-128),k=1..2); i:=i+2 else tb[i]:=64^3*(ls[i]-240)+add(64^(3-k)*(ls[i+k]-128),k=1..3); i:=i+3 end end; ls:=convert(tb,list); if unicodes then ls else cat(seq(`if`(k<256,StringTools[Char](k),sprintf("&#%04d;",k)),k=ls)) end end: #hfl: ReadUTF8 WriteUTF8:=proc(fn::string,s::{string,list(nonnegint)},{overwrite::boolean:=false},$) local ls,tb,i,q,k; if (not(overwrite) and FileTools[Exists](fn)) then error("File exists: #1",fn) end; if type(s,list) then ls:=s else tb:=table(); for i from 1 to length(s) do if (s[i]="&" and s[i+1]="#" and s[i+6]=";") then tb[i]:=parse(s[i+2..i+5]); if type(tb[i],nonnegint) then i:=i+6 else tb[i]:=38 end else tb[i]:=StringTools[Ord](s[i]) end end; ls:=convert(tb,list) end; tb:=table(); for i from 1 to nops(ls) do if (ls[i]<128 ) then tb[4*i]:=ls[i] elif (ls[i]<2048 ) then tb[4*i]:=128+irem(ls[i],64,'q'); tb[4*i-1]:=192+q elif (ls[i]<65536) then tb[4*i]:=128+irem(ls[i],64,'q'); tb[4*i-1]:=128+irem(q,64,'q'); tb[4*i-2]:=224+q else tb[4*i]:=128+irem(ls[i],64,'q'); for k from 1 to 2 do tb[4*i-k]:=128+irem(q,64,'q') end; tb[4*i-3]:=240+q end end; ls:=Vector(map(v->`if`(v<128,v,v-256),convert(tb,list)),datatype=integer[1]); FileTools[Binary][Write](fn,ls,'byteorder'='native'); fclose(fn) end: #hfl: ReadBIN ReadBIN:=proc( filename::string, { code::[nonnegint,nonnegint,nonnegint,nonnegint]:=[0,0,0,0], forcecode::boolean:=false, dimensions::list:=[], datapos::integer:=-1, nodescription::boolean:=false, raw::{boolean,nonnegint}:=false, RMO::boolean:=false, defelementsize::boolean:=true, optrtable::list:=[], buffersize::posint:=10000, printout::boolean:=false },$) local fd,cod,pos,id,storedcode,i,dt,bo,dims,ls,j,dim2,dts,size,Vl,tb,V,p,k1,N,R,i1,i2,i3,i4, v,l,k; fd:=fopen(filename,READ,BINARY); try # code and pos if nodescription then id:=undefined; if (datapos<0) then error("Datapos<0, %1",pos) end; cod,pos:=code,datapos; if (dimensions=[]) then if (cod[3]=0) then cod[3]:=1 elif not(member(cod[3],[1,11,12])) then error("Dimensions are required for multidimensional arrays, code[2]=%1",code[2]) end end else id:=op(FileTools[Binary][Read](fd,integer[4],1,byteorder='native')); if not(member(id,binformatidentifiers1)) then error("Wrong file identifier") end; pos:=`if`(datapos<0,32,datapos); storedcode:=readbytes(fd,4); cod:=storedcode; for i from 1 to 4 do if (code[i]>0 and code[i]<>storedcode[i]) then cod[i]:=code[i] else cod[i]:=storedcode[i] end end; if (cod<>storedcode) then if forcecode then WARNING("Format code is changed from %1 to %2",storedcode,cod) else error("Stored format code %1<>%2 provided code",storedcode,cod) end end end; if not(member(cod[1],[1,2]) and member(cod[2],[1,2,4,8]) and member(cod[4],[1,2,3,4])) then error("Unrecognized format code, %1",cod) end; # datatype, byteorder, and dimensions dt:=[integer,float][cod[1]][cod[2]]; bo:=piecewise(cod[4]=1,'native',cod[4]=2,'network',cod[4]=3,'big',cod[4]=4,'little',NULL); dims:=`if`(dimensions=[], `if`(nodescription, [(FileTools[Size](fd)-pos)/cod[2]], FileTools[Binary][Read](fd,integer[4], `if`(member(cod[3],[1,11,12,13,14,15,3,31,32]),1, `if`(member(cod[3],[2,21,22,53]),2, `if`(cod[3]=54,3, cod[3]-40))),'byteorder'=bo)), dimensions); # special bytes if (cod[3]=13) then ls:=readbytes(fd,20); for j from 1 to nops(ls) while (ls[j]>0) do end; dim2:=j-1; dts:=[seq([integer,float][1+iquo(v,128)][irem(v,128)],v=ls[..dim2])] end; # read data size:=`if`(member(cod[3],[3,31,32,53,54]),dims[1]*(dims[1]+1)/2*mul(v,v=dims[2..]),mul(v,v=dims)); if printout then printf("id=%d, code=[%d,%d,%d,%d], pos=%d, size=%d\n",id,op(cod),pos,size) end; if type(raw,nonnegint) then size:=min(size,raw) end; filepos(fd,pos); if (cod[3]=13) then Vl:=FileTools[Binary][Read](fd,dt,dims[1],'byteorder'=bo); size:=add(l,l=Vl); tb:=table(); for j from 1 to dim2 do V:=Vector(size,datatype=dts[j]); for i from 0 by buffersize to size-1 do ls:=FileTools[Binary][Read](fd,dts[j],min(buffersize,size-i),'byteorder'=bo); V[i+1..i+nops(ls)]:=Vector(ls,datatype=dts[j]) end; tb[j]:=copy(V) end else V:=Vector(size,datatype=dt); for i from 0 by buffersize to size-1 do ls:=FileTools[Binary][Read](fd,dt,min(buffersize,size-i),'byteorder'=bo); V[i+1..i+nops(ls)]:=Vector(ls,datatype=dt) end end; catch: error finally fclose(fd) end try; # compose rtable if defelementsize then dt:=op(0,dt) end; if not(raw=false) then [id,cod,dims,`if`(size>buffersize,V,convert(V,list))] elif (cod[3]=1 ) then Vector (dims[1],V,datatype=dt,op(optrtable)) elif (cod[3]=11) then Vector[column](dims[1],V,datatype=dt,op(optrtable)) elif (cod[3]=12) then Vector[row] (dims[1],V,datatype=dt,op(optrtable)) elif (cod[3]=13) then if (dim2=1) then V:=Vector(dims[1]); k1:=0; for i from 1 to dims[1] do V[i]:=[seq(tb[1][k],k=k1+1..k1+Vl[i])]; k1:=k1+Vl[i] end; V else V:=Vector(dims[1]); k1:=0; for i from 1 to dims[1] do V[i]:=[seq([seq(tb[j][k],j=1..dim2)],k=k1+1..k1+Vl[i])]; k1:=k1+Vl[i] end; V end elif (cod[3]=2 or cod[3]=21) then Matrix(op(dims),(i,j)->V[j+dims[2]*(i-1)],datatype=dt,op(optrtable)) elif (cod[3]=22) then Matrix(op(dims),(i,j)->V[i+dims[1]*(j-1)],datatype=dt,op(optrtable)) elif (cod[3]=3 or cod[3]=31) then Matrix(dims[1],(i,j)->V[i*(i-1)/2+j],datatype=dt,shape=symmetric,storage=triangular[lower],op(optrtable)) elif (cod[3]=32) then Matrix(dims[1],(i,j)->V[j*(j-1)/2+i],datatype=dt,shape=symmetric,storage=triangular[lower],op(optrtable)) elif (cod[3]=41) then Array(1..dims[1],(i1)->V[i1],datatype=dt,op(optrtable)) elif (cod[3]=42) then if RMO then Array(1..dims[1],1..dims[2],(i1,i2)->V[i2+dims[2]*(i1-1)],datatype=dt,op(optrtable)) else Array(1..dims[1],1..dims[2],(i1,i2)->V[i1+dims[1]*(i2-1)],datatype=dt,op(optrtable)) end elif (cod[3]=43) then if RMO then Array(1..dims[1],1..dims[2],1..dims[3],(i1,i2,i3)->V[i3+dims[3]*(i2-1+dims[2]*(i1-1))],datatype=dt,op(optrtable)) else Array(1..dims[1],1..dims[2],1..dims[3],(i1,i2,i3)->V[i1+dims[1]*(i2-1+dims[2]*(i3-1))],datatype=dt,op(optrtable)) end elif (cod[3]=44) then if RMO then Array(1..dims[1],1..dims[2],1..dims[3],1..dims[4],(i1,i2,i3,i4)->V[i4+dims[4]*(i3-1+dims[3]*(i2-1+dims[2]*(i1-1)))],datatype=dt,op(optrtable)) else Array(1..dims[1],1..dims[2],1..dims[3],1..dims[4],(i1,i2,i3,i4)->V[i1+dims[1]*(i2-1+dims[2]*(i3-1+dims[3]*(i4-1)))],datatype=dt,op(optrtable)) end elif (cod[3]=53) then N:=dims[1]*(dims[1]+1)/2; R:=Array(symmetric12,1..dims[1],1..dims[1],1..dims[2],datatype=dt,op(optrtable)); for i1 from 1 to dims[1] do for i2 from 1 to i1 do for i3 from 1 to dims[2] do R[i1,i2,i3]:=V[i1*(i1-1)/2+i2+N*(i3-1)] end end end; R elif (cod[3]=54) then N:=dims[1]*(dims[1]+1)/2; R:=Array(symmetric12,1..dims[1],1..dims[1],1..dims[2],1..dims[3],datatype=dt,op(optrtable)); for i1 from 1 to dims[1] do for i2 from 1 to i1 do for i3 from 1 to dims[2] do for i4 from 1 to dims[3] do R[i1,i2,i3,i4]:=V[i1*(i1-1)/2+i2+N*(i3-1+dims[3]*(i4-1))] end end end end; R else error("Unrecognized shape, %1",cod[3]) end end: #hfl: ReadBIN WriteBIN:=proc( filename::string, R::{list,Vector,Matrix,Array}, opt::list:=[], { code::[nonnegint,nonnegint,nonnegint,nonnegint]:=[0,0,0,1], dimensions::list:=[], nodescription::boolean:=false, overwrite::boolean:=false },$) local cod,dt,v,dims,dim2,dts,bo,fd,pos,j, u,i; # determine datatype cod:=code; if (cod[1]*cod[2]>0) then dt:=[integer,float][cod[1]][cod[2]] else dt:=rtable_options(R,'datatype'); v:=op(0,dt); if (v='integer') then cod[1]:=1; cod[2]:=op(1,dt) elif (v='float') then cod[1]:=2; cod[2]:=op(1,dt) elif type(R[1],list) then cod[1]:=1; cod[2]:=4; dt:=integer[4] else error("Unsupported element type %1",v) end end; # determine shape and dimensions if (cod[3]>0 and dimensions<>[]) then dims:=dimensions else if type(R,Vector) then dims:=[op(1,R)]; cod[3]:=`if`(type(R[1],list),13,1) elif type(R,Matrix) then v:=MatrixOptions(R,'shape'); if (v=[]) then cod[3],dims:=2,[op(1,R)] elif (v=['symmetric']) then cod[3],dims:=3,[[op(1,R)][1]] else error("Unsupported shape %1",v) end elif type(R,Array) then dims:=[rtable_dims(R)]; if (nops(dims)>4) then error("Too large number of dimensions: %1",nops(dims)) end; if not(foldl(`and`,seq(evalb(op(1,v)=1),v=dims))) then error("Each dimension of Array must be enumeratted from 1, but received %1",dims) end; cod[3]:=40+nops(dims); dims:=map2(op,2,dims) end end; # process additional arguments if (cod[3]=13) then dim2:=`if`(type(R[1][1],list),nops(R[1][1]),1); if (dim2>20) then error("Too large secondary dimension (max 20), %1",dim2) end; if (nops(opt)=1) then dts:=[opt[1]$dim2] elif (nops(opt)=dim2) then dts:=opt else error("Unrecognized datatypes: %1",opt) end; dts:=map(v->`if`(v='integer',integer[4],`if`(v='float',float[8],v)),dts) end; # prepare writing bo:=['native','network','big','little'][cod[4]]; if (FileTools[Exists](filename) and not(overwrite)) then error("File exists, %1",filename) end; fd:=fopen(filename,WRITE,BINARY); # write description if not(nodescription) then FileTools[Binary][Write](fd,integer[4],[1760568055],'byteorder'=bo); writebytes(fd,cod); FileTools[Binary][Write](fd,integer[4],dims,'byteorder'=bo); pos:=4+4+4*nops(dims); if (cod[3]=13) then writebytes(fd,map(v->`if`(op(0,v)='float',128,0)+op(1,v),dts)); pos:=pos+dim2 end; writebytes(fd,[0$(32-pos)]) end; # write data if type(R,list) then FileTools[Binary][Write](fd,dt,R,'byteorder'=bo) elif (cod[3]=1) then FileTools[Binary][Write](fd, `if`(cod[1]*cod[2]>0, Vector(R,'datatype'=dt), R), 'byteorder'=bo) elif (cod[3]=13) then FileTools[Binary][Write](fd,Vector(dims[1],i->nops(R[i]),'datatype'=dt),'byteorder'=bo); if (dim2=1) then FileTools[Binary][Write](fd,dts[1],[seq(op(v),v=R)],'byteorder'=bo) else for j from 1 to dim2 do FileTools[Binary][Write](fd,dts[j],[seq(seq(u[j],u=v),v=R)],'byteorder'=bo) end end elif (cod[3]=2) then FileTools[Binary][Write](fd, `if`(cod[1]*cod[2]>0, Matrix(LinearAlgebra[Transpose](R),'datatype'=dt), LinearAlgebra[Transpose](R)), 'byteorder'=bo) elif (cod[3]=3) then FileTools[Binary][Write](fd, Vector([seq(seq(R[i,j],i=1..j),j=1..dims[1])],'datatype'=dt), 'byteorder'=bo) else FileTools[Binary][Write](fd,R,'byteorder'=bo) end; fclose(fd); FileTools[Size](filename) end: #hfl: PrintFile PrintFile:=proc(filename::{string,list(string)},nlines::nonnegint:=0,nendlines::nonnegint:=0,$) local ls,s,i; if (nlines=0 and nendlines=0) then if type(filename,string) then printf("%s\n",FileTools[Text][ReadFile](filename)) else seq(printf("%s\n",s),s=filename) end elif (nendlines=0) then ls:=`if`(type(filename,string),ReadLines(filename,':-nlines'='nlines'),filename[..min(nlines,nops(filename))]); for s in ls do printf("%s\n",s) end; if (nops(ls)=nlines) then printf("...\n") end else ls:=`if`(type(filename,string),StringTools[Split](FileTools[Text][ReadFile](filename),"\n"),filename); if (nlines+nendlinesf2 and not(overwrite) and FileTools[Exists](f2)) then error("File exists, %1",f2) end; FileTools[Binary][Open](f1); ls:=FileTools[Binary][Read](f1,integer[1]); FileTools[Binary][Close](f1); tb:=table(): for i from 1 to nops(ls) do if (ls[i]=13 and if2 and not(overwrite) and FileTools[Exists](f2)) then error("File exists, %1",f2) end; if force then WriteLines(f2,ReadLines(f1),':-overwrite'=true); true else fd:=fopen(f1,READ,TEXT); readline(fd); i:=filepos(fd); s:=readline(fd); try filepos(fd,i); s2:=readline(fd) catch: end; fclose(fd); if (s<>s2) then WriteLines(f2,ReadLines(f1),':-overwrite'=true) end; evalb(s<>s2) end end: #hfl: OpenFolder OpenFolder:=proc(filename::string,{noerror::boolean:=false},$) local s; s:=SimplifyPath(ExpandPath(filename),'replacedirsepto'=kernelopts('dirsep')); if not(FileTools[Exists](s)) then if noerror then ssystem(cat("explorer /select,",s)) else error("File or folder does not exist: %1",s) end else if FileTools[IsDirectory](s) then ssystem(cat("explorer /root,",s)) else ssystem(cat("explorer /select,",s)) end end; NULL end: ################################################################################ #cat: Data manipulation #hfl: RootMeanSquare RootMeanSquare:=proc(V::indexable,$) local n,a,v; n:=add(1,v=V); a:=add(v,v=V)/n; sqrt(add((v-a)^2,v=V)/n) end: #hfl: LinearFit2 LinearFit2:=proc( XorXY::{list,Vector,Matrix}, Yor::{list,Vector}:=[], o::posint:=1, ex::list(nonnegint):=[], output::string:="", {printout::boolean:=false, format::string:="%.3g"}, $) local X,Y,lsp,n,tb,x,p,d,std,Y0,v,r2,c,cstd,out,lo,e; if (Yor=[]) then if type(XorXY,list) then X,Y:=map2(op,1,XorXY),map2(op,2,XorXY) else X,Y:=XorXY[1,..],XorXY[2,..] end else X,Y:=XorXY,Yor end; lsp:=remove(member,[$0..o],ex); n:=nops(lsp); tb:=table(Statistics[LinearFit]([seq(x^p,p=lsp)],X,Y,x,':-output'='solutionmodule'):-Results()); d:=tb["degreesoffreedom"]; std:=`if`(assigned('tb["residualstandarddeviation"]'),tb["residualstandarddeviation"],undefined); Y0:=Statistics[Mean](Y); r2:=`if`(assigned('tb["residuals"]'),1-add(v^2,v=tb["residuals"])/add((v-Y0)^2,v=Y),undefined); c:=convert(tb["parametervalues"],list); cstd:=`if`(d=0,[0$n],convert(tb["standarderrors"],list)); if printout then printf("n=%d, o=%d, d=%d, std=%s, r2=%s\n",n,o,d,Trim(sprintf(format,std)),Trim(sprintf(format,r2))); printf("c =[%{c,}s]\ncstd=[%{c,}s]\n",Vector(map2(sprintf,format,c)),Vector(map2(sprintf,format,cstd))) end; out:=StringTools[Split](output,","); lo:=[]: if (out=[""]) then op(tb) else for e in out do if (e="f" ) then lo:=[op(lo),unapply(tb["leastsquaresfunction"],x)] elif (e="std" ) then lo:=[op(lo),std] elif (e="r2" ) then lo:=[op(lo),r2] elif (e="c" ) then lo:=[op(lo),c] elif (e="cstd") then lo:=[op(lo),cstd] elif (e="cint") then lo:=[op(lo),convert(tb["confidenceintervals"],list)] elif (e="res" ) then lo:=[op(lo),`if`(d=0,map(v->0,Y),tb["residuals"])] elif (e="cov" ) then lo:=[op(lo),`if`(d=0,Matrix(n),tb["variancecovariancematrix"])] end end; op(lo) end end: #hfl: LinearFitTrig LinearFitTrig:=proc( m1::{nonnegint,list(nonnegint)}, m2::{nonnegint,list(nonnegint)}, P::numeric, XY::list, c1::numeric:=1, XY1::list:=[], c2::numeric:=1, XY2::list:=[], output::string:="c,f,std", { weights::{list,procedure,numeric}:=((x,y,p)->1),maxcondnum::numeric:=1e-4, digits::[posint,nonnegint]:=[6,3], printout::boolean:=false },$) local cc,lsm1,lsm2,phi,bas,bas1,bas2,nb,n,n1,n2,m,np,W,B,b,i,XXi,XX,CN,Y,C,R,RMS,RMS0,RMS1,RMS2,f,out,lsm, V,v,s; cc:=evalf(2*Pi/P); lsm1:=`if`(type(m1,list),sort(m1),[$0..m1]); lsm2:=`if`(type(m2,list),sort(m2),[$1..m2]); bas :=[seq( cos(m*phi),m=lsm1),seq( sin(m*phi),m=lsm2)]; bas1:=[seq(-m *sin(m*phi),m=lsm1),seq( m *cos(m*phi),m=lsm2)]; bas2:=[seq(-m^2*cos(m*phi),m=lsm1),seq(-m^2*sin(m*phi),m=lsm2)]; nb:=nops(bas); n,n1,n2:=nops(XY),nops(XY1),nops(XY2); m:=n+n1; np:=m+n2; if type(weights,list) then if (nops(weights)<>np) then error("nops(weights)=%1<>np=%2",nops(weights),np) end; W:=weights else if type(weights,numeric) then if (np<>n) then error("Boltzmann weights are not supported for derivatives") else W:=[seq(evalf(exp(-weights*v[2])),v=XY)] end else W:=[seq(evalf(weights(op(v),0)),v=XY),seq(evalf(weights(op(v),1)),v=XY1),seq(evalf(weights(op(v),2)),v=XY2)] end; W:=W/max(W) end; W:=LinearAlgebra[DiagonalMatrix](W,datatype=float); B:=Matrix(np,nb,datatype=float); for b from 1 to nb do for i from 1 to n do B[ i,b]:=evalf(eval( bas [b],phi=cc*XY [i][1])) end; for i from 1 to n1 do B[n+i,b]:=evalf(eval(c1*bas1[b],phi=cc*XY1[i][1])) end; for i from 1 to n2 do B[m+i,b]:=evalf(eval(c2*bas2[b],phi=cc*XY2[i][1])) end end; XXi:=Matrix(LinearAlgebra[Transpose](B).W.B,shape=symmetric); XX:=Matrix(LinearAlgebra[MatrixInverse](XXi),shape=symmetric,datatype=float); CN:=LinearAlgebra[Norm](XX.XXi-LinearAlgebra[IdentityMatrix](nb),'Frobenius'); if (CN>maxcondnum) then error("Ill conditioned basis, condnum=%1",CN) end; Y:=Vector(np,datatype=float); for i from 1 to n do Y[ i]:=XY [i][2] end; for i from 1 to n1 do Y[n+i]:=XY1[i][2]*c1/cc end; for i from 1 to n2 do Y[m+i]:=XY2[i][2]*c2/cc^2 end; C:=XX.LinearAlgebra[Transpose](B).W.Y; R:=Y-B.C; RMS:=evalf(LinearAlgebra[Norm](R,2)/sqrt(np)); RMS0:=`if`(n =0 ,undefined,evalf(LinearAlgebra[Norm](R[ ..n],2)/sqrt(n))); RMS1:=`if`(n1=0 or c1=0,undefined,evalf(LinearAlgebra[Norm](R[n+1..m],2)/sqrt(n1)/c1)); RMS2:=`if`(n2=0 or c2=0,undefined,evalf(LinearAlgebra[Norm](R[-n2.. ],2)/sqrt(n2)/c2)); if printout then lsm:=convert(`union`({op(lsm1)},{op(lsm2)}),list); printf("%*d %*d\n",digits[1]-digits[2],0,digits[1],Vector(max(lsm),m->m)); V:=Vector(1+max(lsm1),m->""); for m from 1 to nops(lsm1) do V[1+lsm1[m]]:=sprintf("%.*f",digits[2],C[m]) end; printf("%*s\n",digits[1],V); if (lsm2<>[]) then V:=Vector(1+max(lsm2),m->""); for m from 1 to nops(lsm2) do V[1+lsm2[m]]:=sprintf("%.*f",digits[2],C[nops(lsm1)+m]) end; printf("%*s\n",digits[1],V) end; if (n=np) then printf("RMS=%.3g\n",RMS) else printf("RMS=%.3g (0=%.3g, 1=%.3g, 2=%.3g)\n",RMS,RMS0,RMS1,RMS2) end end; f:=unapply(subs(phi=cc*phi,add(C[b]*bas[b],b=1..nb)),phi); out:=StringTools[Split](output,","); seq(`if`(s="c",C,`if`(s="f",op(f),`if`(s="std",RMS,NULL))),s=out) end: #hfl: MultiDimFit MultiDimFit:=proc(X0::{listlist,Vector,Matrix}, Y::Vector, basis::string, order0::list(nonnegint), rng::list({numeric,numeric..numeric}):=[0$nops(order0)], output::string:="f", { cmin::numeric:=10^(2-Digits), printout::boolean:=false},$) local order,d,bf,l,m,x,nx,d2,X,a,b,xi1,xi2,x1,x2,nb,iidB,j,p,B,C,cmin2,CM, v,e; order:=order0; d:=nops(order); if (length(basis)<>d) then error("Inconsistent sizes of basis=%1 and order=%2",basis,order) end; for l from 1 to d do if (basis[l]="t" and type(order[l],odd)) then order[l]:=order[l]+1 end end; bf:=Array(1..d,0..max(order)); for l from 1 to d do for m from 0 to order[l] do bf[l,m]:=unapply( `if`(basis[l]="t",`if`(type(m,even),cos(m/2*x),sin((m+1)/2*x)), `if`(basis[l]="c",orthopoly['T'](m,x), undefined)),x) end end; if type(X0,Matrix) then nx,d2:=op(1,X0); X:=X0 else nx,d2:=Dim2(X0),Dim2(X0[1]); X:=Matrix(nx,d2,(i,l)->X0[i][l],datatype=float) end; if (d2<>d) then error("Inconsistent sizes of X and order") end; if (nops(rng)<>d) then error("Inconsistent sizes of rng and order") end; a,b:=Vector(d,datatype=float),Vector(d,datatype=float); for l from 1 to d do if (basis[l]="t") then xi1,xi2:=0,2*Pi elif (basis[l]="c") then xi1,xi2:=-1,1 else error("Unrecognized basis: %1",basis[l]) end; if (rng[l]=0) then x1,x2:=min(X[..,l]),max(X[..,l]) elif type(rng[l],numeric) then x1,x2:=0,rng[l] else x1,x2:=op(rng[l]) end; a[l],b[l]:=evalf((xi2-xi1)/(x2-x1)),evalf((x2*xi1-x1*xi2)/(x2-x1)) end; nb:=mul(v+1,v=order); iidB:=Matrix(nb,d,datatype=integer); for j from 1 to nb do p:=j-1; for l from 1 to d do iidB[j,l]:=irem(p,order[l]+1,'p') end end; B:=Matrix(nx,nb,(i,j)->evalf(mul(bf[l,iidB[j,l]](a[l]*X[i,l]+b[l]),l=1..d)),datatype=float); C:=LinearAlgebra[LeastSquares](B,Y); cmin2:=cmin*LinearAlgebra[Norm](C,infinity); for j from 1 to nb do if (abs(C[j])M0[i],datatype=float)); d,n:=Dim2(M); maxlen:=1+floor(log10(n)); maxr2:=`if`(type(maxr20,Vector),maxr20,Vector(n,k->maxr20,datatype=float)); effr2:=Mean(maxr2); cor:=evalf(d/2*(GAMMA(d/2)-GAMMA(d/2,effr2/2))/(GAMMA(d/2+1)-GAMMA(d/2+1,effr2/2))); cor0:=`if`(effr20>0,evalf(d/2*(GAMMA(d/2)-GAMMA(d/2,effr20/2))/(GAMMA(d/2+1)-GAMMA(d/2+1,effr20/2))),1); K2:=K0; for iter from 1 to niter do K:=K2; nk:=nops(K); V:=Vector(d,o->Mean(M[o,K]),datatype=float); dM:=Matrix(d,n,(o,k)->M[o,k]-V[o],datatype=float); dM2:=Matrix(d,(o1,o2)->dM[o1,K].dM[o2,K]/(nk-1)*`if`(iter=1,cor0,cor),shape=symmetric,datatype=float); dM2i:=dM2^(-1); r2:=Vector(n,k->Transpose(dM[..,k]).dM2i.dM[..,k],datatype=float); # add(v,v=r2[K])/(nk-1)=d if printout then printf("nk=%*d, mean=[%+{c,}.3f], sigma=[%{c,}.3f]\n",maxlen,nk,V,map(sqrt,Eigenvalues(dM2))) end; K2:=select(k->r2[k][op(v),v[3]-v[2]],[seq([v,evalf(nops(K2)*GAMMA(d/2,v/2)/(GAMMA(d/2)-GAMMA(d/2,effr2/2))),add(`if`(r2[k]evalb(v[1]>x1 and v[1]0..0 and printout) then printf("After xrng(%a): np=%d, x-range=%g..%g\n",xrng,nops(xy1),min(seq(v[1],v=xy1)),max(seq(v[1],v=xy1))) end; if (cythr<>0.5) then y1,y2:=min(seq(v[2],v=xy1)),max(seq(v[2],v=xy1)); ythr:=y1+cythr*(y2-y1); if (cythr>0 and cythr<0.5) then xy2:=select(v->evalb(v[2]0.5) then xy2:=select(v->evalb(v[2]>ythr),xy1) else error("Expected 00) then xy2:=`if`(cythr<0.5,Sort,Rort)(xy2,[2]); if (nops(xy2)<3) then error("3 points are required for pre-interpolation, %2 points are selected by ythr=%3",3,nops(xy2),ythr) end; if (nops(xy2)evalb(v[1]>x1 and v[1]1) then WARNING("Multiple extrema %1",sol); x0:=sol[op(MinIdx([seq(abs(v-x0),v=sol)]))] else x0:=op(sol) end; if not(type(x0,numeric)) then error("No extremum") end; ypp:=add(n*(n-1)*c[n+1]*x0^(n-2),n=2..o); [x0, sqrt(evalf[2*Digits](add(add(n*m*cov[n+1,m+1]*x0^(n+m-2),n=1..o),m=1..o)))/abs(ypp), add(c[n+1]*x0^n,n=0..o), sqrt(evalf[2*Digits](add(add(cov[n+1,m+1]*x0^(n+m),n=0..o),m=0..o))), ypp, unapply(add(c[n+1]*x^n,n=0..o),x), fitrng] end: #hfl: ReorderXYdata ReorderXYdata:=proc( xydata, n::posint:=Dim2(xydata), m::posint:=Dim2(xydata[1][2]), { order::posint:=2, xindexing::procedure:=((A,i)->A[i][1]), yindexing::procedure:=((A,i,j)->A[i][2][j]) },$) local X,Y,y,i,j,id; X:=Vector(n,i->xindexing(xydata,i)); Y:=Matrix(n,m,(i,j)->yindexing(xydata,i,j)); y:=Vector(m); for i from 1 to n do Y[i,..]:=sort(Y[i,..]) end; for i from order+2 to n do id:=SortIdx(Vector(m,j->evalf(CurveFitting[PolynomialInterpolation](X[i-order-1..i-1],Y[i-order-1..i-1,j],evalf(X[i])))),'nolist'); for j from 1 to m do y[id[j]]:=Y[i,j] end; Y[i,..]:=y end; [seq([X[i],Y[i,..]],i=1..n)] end: #hfl: Intersection Intersection:=proc(xydata, n::posint:=Dim2(xydata), { order::nonnegint:=1, xindexing::procedure:=((A,i)->A[i][1]), yindexing::procedure:=((A,i,j)->A[i][2][j]) },$) local X,Y,i1,i2,s,i,x,e1,e2; X:=Vector(n,i->xindexing(xydata,i)); Y:=Matrix(n,2,(i,j)->yindexing(xydata,i,j)); if (order=0) then i1,i2:=1,n else s:=signum(Y[1][1]-Y[1][2]); for i from 2 to n while s*(Y[i][1]-Y[i][2])>0 do end: if (i>n) then error "No intersection is found" end; i1,i2:=i-iquo(order+2,2),i+iquo(order,2); if (type(order,even) and i1=0 and i2<=n ) then i1:=1 elif (type(order,even) and i1>=1 and i2=n+1) then i2:=n elif (i1<1 or i2>n) then error "For order=%1 provide at least %2 interpolation points",order,2*iquo(order+2,2) elif (type(order,even)) then if (abs(Y[i1][1]-Y[i1][2])0) then dx:=(x2-x1)*rdx; x1,x2:=x1-dx,x2+dx end; if (rdy<>0) then dy:=(y2-y1)*rdy; y1,y2:=y1-dy,y2+dy end; x1,x2,y1,y2 end: GetBox2:=proc(L::list,$) local R,A; R:=map(plottools[getdata],L,'rangesonly'); if (R=[]) then error("No range data") end; A:=Array(1..2,1..2,1..nops(R),(i,j,k)->op(i,R[k][j])); min(A[1,1,..]),max(A[2,1,..]),min(A[1,2,..]),max(A[2,2,..]) end: Spline2:=proc(XY,cod::{0,1,2,12}:=0,{degree::posint:=3}) local G,e,x; if (cod>0) then G:=`if`(cod=12,Matrix(2,4*degree+1,{(1,1)=1,(2,2*degree+2)=1}),Matrix(1,4*degree+1,{(1,`if`(cod=1,1,2*degree+2))=1})) end; e:=CurveFitting[Spline](XY,x,':-degree'=degree,`if`(cod>0,'endpoints'=G,NULL),_rest); unapply(simplify(piecewise(xXY[-1][1],'undefined',e)),x) end: ErrorBar:=proc(x0,r::range,w::realcons:=0,opt::list:=[],$) local x,y1,y2,dx,e; x,y1,y2:=evalf(x0),op(evalf(r)); e:=plottools[line]([x,y1],[x,y2],op(opt)): if (evalf(w)>0) then dx:=evalf(w/2); CURVES([[x-dx,y1],[x+dx,y1]],[[x-dx,y2],[x+dx,y2]],op(e)) else e end end: Legend2:=proc( L::list(list), x::numeric, y::numeric, a::string, h::numeric, h2w::numeric:=0.55, symb::string:="line", {symbw::numeric:=`if`(symb="line",2,1), sep::numeric:=0, imx::numeric:=3/4, imy::numeric:=1/4, omx::numeric:=1/2, omy::numeric:=1/4, opts::list:=[thickness=2], optt::list:=[], optb::list:=[color="Gray",thickness=0], nobox::boolean:=false},$) local n,ax,ay,c,w,maxw,x0,y0,x1,y1,x2,y2,Opts,Optt,k; n:=nops(L); ax,ay:=0,0; for c in a do if (c="L") then ax:=-1 elif (c="R") then ax:=1 elif (c="A") then ay:=1 elif (c="B") then ay:=-1 else error("Unrecognized alignment %1",a) end end; w:=h*h2w; maxw:=symbw+sep+max(map(v->length(v[1]),L));# max width x0:=piecewise(ax=+1,x+w*(imx+omx ),ax=-1,x-w*(imx+omx+maxw ),x-w*maxw/2);# determine position of the top left corner of the first symbol y0:=piecewise(ay=-1,y-h*(imy+omy+1/2),ay=+1,y+h*(imy+omy+n-1/2),y+h*(n-1)/2); x1,x2,y1,y2:=x0-w*imx,x0+w*(maxw+imx),y0+h*(1/2+imy),y0-h*(n-1/2+imy); Opts:=`if`(nops(opts)=n and type(opts,list(list)),opts,[opts$n]); Optt:=`if`(nops(optt)=n and type(optt,list(list)),optt,[optt$n]); op(`if`(symb="line", [seq(plot([[x0,y0-h*(k-1)],[x0+w*symbw,y0-h*(k-1)]],color=L[k][2], op(Opts[k])),k=1..n)], [seq(plot([[x0+w*symbw/2,y0-h*(k-1)]] ,color=L[k][2],style="point",symbol=symb,op(Opts[k])),k=1..n)])), seq(textplot([x0+w*(symbw+sep),y0-h*(k-1),L[k][1]],color=L[k][2],align={"right"},op(Optt[k])),k=1..n), `if`(nobox,NULL,plot([[x1,y1],[x2,y1],[x2,y2],[x1,y2],[x1,y1]],op(optb))) end: TextOnLine:=proc(x1::realcons,y1::realcons,x2::realcons,y2::realcons,s::string, h::numeric,h2v0::numeric:=0.55,g0::numeric:=0.03,rxy::list(numeric):=[0.5], {textover::boolean:=false}) local h2v,g,w,rx,ry,x0,y0,Dx,Dy,sx,sy,dx,dy; h2v,g:=max(h2v0,g0),min(h2v0,g0); w:=evalf(length(s)*h*h2v); if (nops(rxy)=1) then rx,ry:=rxy[1],rxy[1] elif (nops(rxy)=2) then rx,ry:=rxy[1],rxy[2] else error("Unrecognized rxy: %1",rxy) end; x0,y0:=x1*(1-rx)+x2*rx,y1*(1-ry)+y2*ry; Dx,Dy:=evalf(x2-x1),evalf(y2-y1); if textover then [plot([[x1+g*Dx,y1+g*Dy],[x2-g*Dx,y2-g*Dy]],_rest), plottools[rectangle]([x0-w/2,y0-h/2],[x0+w/2,y0+h/2],'color'="White",'style'="polygon"), plots[textplot]([x0,y0,s],_rest)] else sx,sy:=`if`(Dx<0,-1,1),`if`(Dy<0,-1,1); if (w*Dy^211) then ls:=[seq(op(ls),i=1..ceil(n/11))][..n] end end; if reverse then ls:=ListTools[Reverse](ls) end; [op(before),op(ls),op(after)] end end: WritePlotData:=proc( filename::string, ls::list(list(list)), Xformat::string, Yformat::string, com::string:="", { Xlabel::string:="Xvalue", legend::list(string):=[], emptycell::string:="", sep::character:="\t", overwrite::boolean:=false },$) local X,Y,m,n,i,j,v,fmt,u; n:=nops(ls); X:=Sort(convert({seq(seq(sprintf(Xformat,u[1]),u=v),v=ls)},list),v->parse(v)); m:=nops(X); Y:=Matrix(m,n,(i,j)->emptycell); for j from 1 to n do for v in ls[j] do i:=SearchPos(X,sprintf(Xformat,v[1])); if (i>0) then Y[i,j]:=sprintf(Yformat,v[2]) else error("Cannot find index of x-value for data point %1 in %2",v,ls) end end end; if (nops(legend)>0 and nops(legend)<>n) then error("Inconsistent legend %1 for data %2",legend,ls) end; fmt:=cat("%s",cat(sep,"%s")$n); WriteLines(filename,[ `if`(com="",NULL,cat("# ",com)), `if`(legend=[],NULL,sprintf("%s%s%{c*}s",Xlabel,sep,sep,Vector(legend))), seq(sprintf(fmt,X[i],seq(Y[i,j],j=1..n)),i=1..m)], ':-overwrite'=overwrite) end: plotps:=proc( filename0::string, pl, opt::string:=psplotoptions, cbb::[integer,integer,integer,integer]:=[0,0,0,0], { waittime::numeric:=0.1, overwrite::boolean:=false, eps2pdf::string:="eps2pdf.bat" }) local p,filename,e,cd,ans; p,filename,e:=ExpandPath(filename0,"p,pn,x"); filename:=`if`(e="" or e=".pdf",cat(filename,".eps"),filename0); if FileTools[Exists](filename) then if overwrite then fremove(filename) else error("File exists: %1",filename) end end; plotsetup('ps','plotoutput'=filename,'plotoptions'=opt); print(plots[display](pl,_rest)); plotsetup('default'); if (cbb<>[0,0,0,0]) then Threads[Sleep](waittime); CorrectEPS(filename,cbb) end; if (e=".pdf") then if FileTools[Exists](filename0) then if overwrite then fremove(filename0) else error("File exists: %1",filename0) end end; if (p<>"") then cd:=currentdir(); currentdir(p) end; ans:=ssystem(cat(eps2pdf," ",ExpandPath(filename,"nx"))); if (p<>"") then currentdir(cd) end; if (ans[1]=0) then fremove(filename) else error("Converter's error: %1",ans) end end; NULL end: CorrectEPS:=proc(filename::string,cbb::list:=[0,0,0,0],$) local ls,fd,s,bb; ls:=ReadLines(filename); fd:=fopen(filename,WRITE,TEXT): for s in ls do if (s[1..8]="%%Pages:") then elif (s[1..14]="%%BoundingBox:") then bb:=sscanf(s[15..-1],"%d %d %d %d"); fprintf(fd,"%%%%BoundingBox: %d %d %d %d\n",op(bb+cbb)) else writeline(fd,s) end end: fclose(fd); NULL end: ImageCoreBox:=proc(img::Array,threshold::numeric:=1,$) local nx,ny,v,x1,y1,x2,y2; ny,nx,v:=Dim2(img); for x1 from 1 to nx while (min(img[..,x1,..])>=threshold) do end; for y1 from 1 to ny while (min(img[y1,..,..])>=threshold) do end; for x2 from nx by -1 to x1 while (min(img[..,x2,..])>=threshold) do end; for y2 from ny by -1 to y1 while (min(img[y2,..,..])>=threshold) do end; [x1,y1,x2,y2] end: ################################################################################ #cat: HTML and XML tools #hfl: SqueezeHTML SqueezeHTML:=proc(s::string,ls::list:=[],$) local s1,i,j,ls1,v; s1:=StringTools[Trim](Squeeze(s)); i:=1; while (i<=length(s1)) do j:=searchtext("<",s1,i..-1); if (j=0) then break else i:=i+j end; if (s1[i]=" ") then s1:=cat(s1[1..i-1],s1[i+1..-1]) end end; i:=1; while (i<=length(s1)) do j:=searchtext(">",s1,i..-1); if (j=0) then break else i:=i+j end; if (i>2 and s1[i-2]=" ") then s1:=cat(s1[1..i-3],s1[i-1..-1]) end end; ls1:=[seq([length(v),cat(v,">")],v=ls)]; i:=1; while (i<=length(s1)) do j:=searchtext("<",s1,i..-1); if (j=0) then break else i:=i+j end; if (i>2 and s1[i-2]=" " and foldl(`or`,false,seq(s1[i..i+v[1]]=v[2],v=ls1))) then s1:=cat(s1[1..i-3],s1[i-1..-1]) end end; ls1:=[seq([length(v),cat("<",v)],v=ls)]; i:=1; while (i<=length(s1)) do j:=searchtext(">",s1,i..-1); if (j=0) then break else i:=i+j end; if (s1[i]=" " and foldl(`or`,false,seq(s1[max(1,i-2-v[1])..i-2]=v[2],v=ls1))) then s1:=cat(s1[1..i-1],s1[i+1..-1]) end end; s1 end: #hfl: SqueezeHTML SqueezeTable:=proc(s0::string,$) local s,i,s1,pos,s2; s:=SqueezeHTML(s0,["tr","/tr","td","/td","th","/th"]); i:=SearchKeySeq(s,[""])-4; s1:=Trim(s[1..i-1],[" ","\n"]); do pos:=SearchKeySeq(s,["",""],'start'=i,'fulloutput','noerror'); if (pos=0 or pos=[]) then break end; s2:=Trim(s[i..pos[1]-5],[" ","\n"]); s1:=cat( s1, `if`(s2="","",cat("\n",s2)), "\n", StringTools[SubstituteAll](Trim(s[pos[1]..pos[2]-6],[" ","\n"]),"\n"," "), ""); i:=pos[2] end; s1:=cat(s1,"\n",Trim(s[i..-1],[" ","\n"])); SqueezeHTML(s1,["tr","/tr","td","/td","th","/th"]) end: #hfl: ReadHref ReadHref:=proc(s0::string) local i,s; s:=SqueezeHTML(StringTools[SubstituteAll](s0,"\n",""),["a","/a"]); i:=searchtext("\">",s,10..-5); if (s[1..9]<>""" or i=0) then error "Incorrect href structure: %1",s else [s[i+11..-5],s[10..i+8]] end end: #hfl: ReadHref ResolveURI:=proc(URI::string,thisfile::string,$)::string; local s,n,ls; if (URI[1..8]="file:///") then URI[9..-1] elif (URI[1..7]="http://") then URI else s:=`if`(URI=".","",URI); n:=0; while (s[1..3]="../") do s:=s[4..-1]; n:=n+1 end; ls:=StringTools[Split](thisfile,"/"); if (nops(ls)true),{all::boolean:=false},$) local tb,i,v; if not(XMLTools[IsElement](xml) or XMLTools[IsTree](xml)) then error("First argument must be XML element or tree but received %1",xml) end; tb:=table(): for i from 1 to XMLTools[ContentModelCount](xml) do v:=XMLTools[GetChild](xml,i); if (XMLTools[IsElement](v) and XMLTools[ElementTypeName](v)=nam and f(v)) then if (XMLTools[ContentModelCount](v)=1 and XMLTools[IsText](XMLTools[GetChild](v,1))) then v:=XMLTools[TextText](XMLTools[GetChild](v,1)) end; if not(all) then return v end; tb[i]:=v end end; `if`(all,convert(tb,list),NULL) end: #hfl: ReadXML ReadXML:=proc(filename::string,{reload::boolean:=false},$) local ff,ls,i,XMLdeclaration,DOCTYPEdeclaration,xml; ff:=proc(x) local t,u,v; if type(x,function) then t:=op(0,x); if (t=_XML_Element) then u:=[XMLTools[ElementName](x),op(XMLTools[AttributeTable](x))]; v:=map(ff,op(3,x)); if (v=[]) then [u] elif (nops(v)=1 and type(v[1],string)) then [[op(u),v[1]]] else [u,op(v)] end elif (t=_XML_Text) then op(1,x) else XMLTools[PrintToString](x) end else x end end; ls:=ReadLines(filename,':-reload'=reload); i:=1; if (TrimLeft(ls[i][..5])="", dtd::string:="", overwrite::boolean:=false},$) local ff,XMLdeclaration,DOCTYPEdeclaration; ff:=proc(x) local v,att,u; v:=x[1]; att:=cat(seq(sprintf(" %s='%s'",op(u)),u=v[2])); if (nops(x)=1) then if (nops(v)=3) then sprintf("<%s%s>%s\n",v[1],att,v[-1],v[1]) else sprintf("<%s%s />\n",v[1],att) end else cat(sprintf("<%s%s>\n",v[1],att),seq(ff(u),u=x[2..]),sprintf("\n",v[1])) end end; WriteLines(filename,[`if`(xmld="",NULL,xmld),`if`(dtd="",NULL,dtd),ff(xml)],':-overwrite'=overwrite) end: ################################################################################ #cat: Physics tools #hfl: SU #SU q2vu:=proc(expr) local q,v,u,e; q:=simplify(expr); e:=selectfun(q,Units[Unit]); if nops(e)=0 then u:=1; v:=q elif nops(e)=1 then u:=op(e); v:=simplify(q/u) else error "Unrecognized format in %1",q end; v,u end: #hfl: GetConstantSU GetConstantSU:=proc(const::name) local opts,e,v,u; opts:=[_passed][2..-1]; e:=ScientificConstants[Constant](const,'system'='SU'); e:=ScientificConstants[GetValue](e)*ScientificConstants[GetUnit](e); if ('esu' in opts) then e:=simplify(e,'esu') elif ('emu' in opts) then e:=simplify(e,'emu') else e:=simplify(e) end; v,u:=q2vu(e); if ('value' in opts) then v elif ('unit' in opts) then u else e end end: #hfl: ConvertSU ConvertSU:=proc(expr,sys::name,$) local e,v,u; if (sys='SI') then v,u:=q2vu(expr); e:=subs('esu'=10/cc*'A'*'s','emu'=10*'A'*'s',parse(convert(op(u),string))) elif (sys='CGS') then v,u:=q2vu(expr); e:=parse(convert(op(u),string)); if has(e,'A') then error "A unit is present in %1",expr end; e:=subs( 'esu'=sqrt('g'*'cm'^3/'s'^2), 'emu'=sqrt('g'*'cm'^3/'s'^2), 'K'=kB*'g'*'cm'^2/'s'^2, 'mol'=NA, 'cd'=10^7/683*'g'*'cm'^2/'s'^3, e) elif (sys='cm') then v,u:=q2vu(ConvertSU(expr,'CGS')); e:=subs('s'=cc*'cm','g'=cc/hb/'cm',parse(convert(op(u),string))) else error "Unrecognized system of units: %1",sys end; simplify(Units[Unit](v*e)) end: ################################################################################ #cat: #hfl: System #toc: System tools #BasicTools/BuildHelp,ReadModules,AssignModule,ProcessSetupArgs,Check4Updates,TabulatedFunction,ReduceFloat,ReduceFloat2,ReduceFloat3,ReduceFloatGrid,Timing,InitTiming BuildHelp:=proc( source::string, hdbfile::string:="", { rootkey::string:="#root:", catkey::string:="#cat:", tockey::string:="#toc:", hflkey::string:="#hfl:", ext::string:=".mw", test::boolean:=false, printout::boolean:=test},$) local lib,folder,s,root,category,n1,n2,n3,n4,S,N,l,helpfile,tocentry,A,topic,bro,v,i; if (hdbfile="") then lib:=NULL else lib:=cat(hdbfile,`if`(MapleVersion<18,".hdb",".help")); if FileTools[Exists](lib) then WARNING("Help-file exists and can not be deleted during Maple session. Build over. For sure, remove help-file.") else HelpTools[Database][Create](lib) end end; folder,s:=ExpandPath(source,"p,n"); if printout then printf("folder=%s\n",folder) end; root:=cat(s,"/"); category:=""; n1,n2,n3,n4:=length(rootkey),length(catkey),length(tockey),length(hflkey); S:=ReadLines(source); N:=nops(S); for l from 1 to N do s:=S[l]; if (s[..n1]=rootkey) then root:=cat(Trim(s[n1+1..]),"/"); if printout then printf("%s\n",s) end elif (s[..n2]=catkey) then v:=Trim(s[n2+1..]); category:=`if`(v="",v,cat(v,"/")); if printout then printf("%s\n",s) end elif (s[..n4]=hflkey) then helpfile:=cat(folder,Trim(s[n4+1..]),ext); if printout then printf("%s\n",s) end; if not(FileTools[Exists](helpfile)) then error("Helpfile \"%1\" does not exist",Trim(s[n4+1..])) end; if (l0) then error "Module in module or no module end" end; s1:=s[i..j]: for k from i-2 by -1 to 1 while (s[k]<>"\n") do end; tb[i]:=cat(prefix,Trim(s[k+1..i-2],[" ",":","="])); assign(convert(tb[i],name),s1) end; convert(tb,list) end: AssignModule:=proc(moduledef::string) local s,m; if (SearchText(" ",moduledef)>0) then s:=moduledef else s:=eval(convert(cat("s_",moduledef),name)); if not(type(s,string)) then ReadModules(moduledef); s:=eval(s) end end; m:=parse(s,'statement'); if member('ModuleInit',m) then m:-ModuleInit(_rest) end; m end: ProcessSetupArgs:=proc(lsargs::list,variables::list(name),constants::list(name):=[],suffix::string:="",$) local ls1,ls2,tb1,i,tb2,sq,e,s,v; ls1:=map(convert,variables,string); ls2:=map(convert,constants,string); tb1:=table(): for i from 1 to nops(ls1) do tb1[ls1[i]]:=variables[i] end; tb2:=copy(tb1): for i from 1 to nops(ls2) do tb2[ls2[i]]:=constants[i] end; sq:=NULL; for e in lsargs do if type(e,equation) then s:=cat(convert(lhs(e),string),suffix); if assigned(tb1[s]) then assign(tb1[s],rhs(e)) else error("Unrecognized variable, %1 (use printout for the list of correct keywords)",s) end else s:=cat(convert(e,string),suffix); if assigned(tb2[s]) then v:=eval(tb2[s]); sq:=sq,`if`(type(v,table),op(v),v) else error("Unrecognized variable or constant, %1 (use printout for the list of correct keywords)",s) end end end; sq end: Check4Updates:=proc( root::string:="http://zhugayevych.me/maple/", packages::list(string):=["BasicTools","SSH","FiniteGroups","LatticeTools","MolMod","AstroTools"], {openlinks::boolean:=false},$) local n,p,v1,v2,s,s2; n:=max(map(length,packages)); for p in packages do v1,v2:=0,0; s:=sprintf("%-*s ",n,p); try v1:=eval(parse(cat(p,":-Setup('version')"))); s:=sprintf("%s%d",s,v1) catch: s:=cat(s,"not installed") end; try s2:=HTTP:-Get(cat(root,p,"/",p,".ini"))[2]; v2:=ReadBody(s2,"version:=",";",'format'="%d"); if (v1=v2) then s:=sprintf("%s=%d",s,v2) elif (v1%d - you have newer version",s,v2) end catch: s:=cat(s," - no online information") end; printf("%s\n",s) end; end: TabulatedFunction:=proc() AssignModule(s_TabulatedFunction,_rest) end: s_TabulatedFunction:="module()\n export ModuleInit,ModuleApply,XV,Remove,Read,Save;\n local f,rho,defdr,x,y;\n \n ModuleInit:=proc(ForXVorFilename,dr::positive,rhof::procedure:=proc(x,y) local v; `if`(type(x,atomic),abs(x-y),evalf(sqrt(add(v^2,v=x-y)))) end,{Fdim::posint:=1},$)\n if type(ForXVorFilename,list) then\n XV:=ForXVorFilename;\n f:=x->'undefined'\n elif type(ForXVorFilename,string) then\n XV:=[];\n Read(filename,Fdim);\n f:=x->'undefined'\n else\n XV:=[];\n f:=ForXVorFilename end;\n defdr:=dr;\n rho:=rhof;\n NULL\n end:\n \n ModuleApply:=proc(x,dr::nonnegative:=defdr,$)\n local ls,i,xv;\n if (XV=[]) then\n xv:=[evalf(x),evalf(f(x))]; XV:=[xv]\n else\n ls:=[seq(rho(x,xv[1]),xv=XV)];\n i:=op(MinIdx(ls));\n if (ls[i]>dr) then xv:=[evalf(x),evalf(f(x))]; XV:=[op(XV),xv] else xv:=XV[i] end end;\n xv\n end:\n \n Remove:=proc(x,$)\n local i,xv;\n if (XV<>[]) then\n i:=op(MinIdx([seq(rho(x,xv[1]),xv=XV)])):\n XV:=subsop(i=NULL,XV) end;\n NULL\n end:\n \n Read:=proc(filename::string,Fdim::posint:=1,convert::procedure:=(V->Vector[column](V)),$)\n local M,n,m,Xdim,i;\n M:=op(fscanf(filename,""%{;h}fm""));\n n,m:=Dim2(M);\n Xdim:=m-Fdim;\n XV:=[op(XV),seq([`if`(Xdim=1,M[i,1],convert(M[i,1..Xdim])),`if`(Fdim=1,M[i,-1],convert(M[i,-Fdim..-1]))],i=1..n)];\n ls\n end:\n \n Save:=proc(filename::string,{digits::posint:=min(10,Digits-2),overwrite::boolean:=false},$)\n local fd,xv;\n if (XV=[]) then\n WARNING(""Nothing to save"")\n else\n if (filename<>"""") then\n if (FileTools[Exists](filename) and not(overwrite)) then error ""File exists, %1"",filename end;\n fd:=fopen(filename,WRITE,TEXT)\n else\n fd:=fopen('terminal',WRITE) end;\n for xv in XV do fprintf(fd,""%.*f %.*f\n"",\n digits, `if`(type(xv[1],float),xv[1],convert(xv[1],Vector)),\n digits, `if`(type(xv[2],float),xv[2],convert(xv[2],Vector))) end;\n if (filename<>"""") then fclose(fd) end end;\n nops(XV)\n end:\n \n end": Timing:=proc(nleft::integer:=nLeft,{nmin::posint:=3,timestep::positive:=10},$) local t,tleft,perc; nDone:=nDone+1; nLeft:=nleft-1; t:=time[real](); if (t>PrintTime and nDone>=nmin or nDone>=nPerc) then if (nLeft<0) then tleft:=infinity; PrintTime:=t+timestep else tleft:=(t-InitialTime)*nLeft/nDone; PrintTime:=t+max(tleft/2,timestep) end; if (nDone>=nPerc) then if (nLeft<0) then nPerc:=2000000000 else Percentages:=Percentages[2..]; nPerc:=round(Percentages[1]*(nDone+nLeft)/100) end end; if (tleft>timestep) then perc:=`if`(nLeft>0,round(nDone/(nDone+nLeft)*100),0); if (tleft<100) then printf("%d%% done, %ds left to finish at %ds\n",perc,round(tleft),round(t+tleft)) elif (tleft<3000) then printf("%d%% done, %dm left to finish at %s\n",perc,round(tleft/60),StringTools[FormatTime]("%R",round(tleft))) elif (tleft<80000) then printf("%d%% done, %dh left to finish at %s\n",perc,round(tleft/3600),StringTools[FormatTime]("%R",round(tleft))) elif (tleft