# A.I.R. # # main file containing all the routines and subroutines. # last updated 23.04.04 #------------------- global variables and read input routines MAXDIR:=400: #MAXDIR is the maximum number of files per directory that we allow #for reasons of O.S. functionality TOPOLOGIES:=convert( map(xx->op(0, xx), indets(ibp_equations, function)), list): #TOPOLOGIES is a nameholder for the integrals, like DX,PB,V etc if nops(TOPOLOGIES)<>1 then ERROR("Wrong format for the IBP equations"); fi: #check of consistency TOPOLOGIES:=[op(TOPOLOGIES), xb]: #xb is an additional name needed internally. prop_entries:=[op(map(nops, select(has, indets(ibp_equations, function), TOPOLOGIES[1])))]: if nops(prop_entries)<>1 then ERROR("Wrong format for the IBP equations"); else prop_entries:=prop_entries[1]: fi: #another consistency check on the syntax of the ibp's prop_indices:=indets(ibp_equations, function): prop_indices:=map( proc(xxx) RETURN(map( proc(yyy) RETURN(op(indets(yyy+1000, symbol))); end, xxx)): end, prop_indices): if nops(prop_indices)<>1 then ERROR("Wrong format for the IBP equations"); else prop_indices:=op(op(prop_indices)): fi: NPROPS:=nops([prop_indices]): ThetaF:=proc(x) if type(x, numeric) then if x>0 then RETURN(1); fi: RETURN(0); fi: RETURN('procname'(x)); end: #ThetaF is a declaration for the step function zero_clauses:=unapply(ZERO_TOPOLOGIES, prop_indices): geneq:=unapply(ibp_equations, prop_indices): assign(TOPOLOGIES[1], (proc() local test_zero: test_zero:=map(is, zero_clauses(args)); test_zero:=has(test_zero, true): if test_zero=true then RETURN(0); fi: RETURN('procname'(args)); end) ): #the ibp's as templates are read. kernelopts(printbytes=false): #forces maple9 not to give info about memory, cpu and time used with(student): #loaading the student package with(combinat): #loading the combinat package (used in SEEDGEN) #------------------ side routines find_top:=proc() # routine that finds the topology of the integral. #It returns the index of positive elements in its argument (the argument must be a sequence) local topol, x, i; if whattype(args)=exprseq then x:=[]; for i from 1 to nargs do if args[i]>0 then x:=[op(x), i]; fi; od: RETURN(x); fi; if type(args, {list, function}) then RETURN(procname(op(args))); fi; #print(args); ERROR(`Wrong type of aguments-EGW`); end: check_values:=map(proc(expr) #routine that converts floating numbers to fractions, used in checking numerical values of coefficients RETURN(map(x->lhs(x)=convert(evalf(rhs(x)), fraction), expr)); end, check_values): Mycollect:=proc(expr, function_symbol) local inds, subi, restterms, activeterms, coeff_inds; #this is an alternative for the collect() routine of maple. It's much faster and more efficient for our purposes inds:=select(has, indets(expr, function), function_symbol): restterms:=subs(map(xxx-> xxx=0, inds), expr): inds:=convert(inds, list): coeff_inds:=map(zzz->coeff(expr, zzz), inds): if nargs=3 then coeff_inds:=map(args[3], coeff_inds): fi: activeterms:=zip((x, y)->x*y, coeff_inds, inds): activeterms:=convert(activeterms, `+`): RETURN(activeterms+restterms); end: #-------------------- ROUTINES FOR THE PRIORITY #Finds the number of propagators of an integral nprops:=proc(expr) local val, i, elem; if not type(expr, specfunc(anything,{op(TOPOLOGIES)})) then ERROR(`Wrong kind of arguments`); fi; val:=0; for i from 1 to nops(expr) do elem:=op(i, expr); if elem>0 then val:=val+1; fi; od: RETURN(val); end: #Finds excess of positive powers of props of an integral sumpp:=proc(expr) local val, i, elem; if not type(expr, specfunc(anything,{op(TOPOLOGIES)})) then ERROR(`Wrong kind of arguments`); fi; val:=0; for i from 1 to nops(expr) do elem:=op(i, expr); if elem>0 then val:=val+(elem-1); fi; od: RETURN(val); end: #Find sum of powers of irred. numer. of an integral sumnp:=proc(expr) local val, i, elem; if not type(expr, specfunc(anything, {op(TOPOLOGIES)})) then ERROR(`Wrong kind of arguments`); fi; val:=0; for i from 1 to nops(expr) do elem:=op(i, expr); if not elem>0 then val:=val-elem; fi; od: RETURN(val); end: val_top:=proc(expr) local val, i, elem; # checks whether the name of the integral is correct if not type(expr, specfunc(anything, {op(TOPOLOGIES)})) then ERROR(`Wrong kind of arguments`); fi; val:=op(0, expr); for i from 1 to nops(TOPOLOGIES) do if val=TOPOLOGIES[i] then RETURN(nops(TOPOLOGIES)-i+1); fi; od: ERROR("Not known function"); end: #Finds integral with the biggest priority, through the laporta criteria priority:=proc(equat) local eqs, ints, aux,elem, topo, i; if type(equat, equation) then eqs:=equat; else eqs:=equat=0; fi; ints:=[op(indets(eqs, specfunc(anything, {op(TOPOLOGIES)})))]; #Apply zero criterion: More important topology... aux:= max(op(map(x->val_top(x), ints))); ints:=map(proc(x) if val_top(x)<>aux then RETURN(NULL); else RETURN(x); fi; end, ints); # Apply first criterion: Bigger number of propagators aux:= max(op(map(x->nprops(x), ints))); ints:=map(proc(x) if nprops(x)<>aux then RETURN(NULL); else RETURN(x); fi; end, ints); # Apply second criterion: Bigger number of extra-positive powers aux:= max(op(map(x->sumpp(x), ints))); ints:=map(proc(x) if sumpp(x)<>aux then RETURN(NULL); else RETURN(x); fi; end, ints); # Apply third criterion: Bigger number of irreducible numerators aux:= max(op(map(x->sumnp(x), ints))); ints:=map(proc(x) if sumnp(x)<>aux then RETURN(NULL); else RETURN(x); fi; end, ints); if ints = [] then RETURN([]): else RETURN(ints[1]); fi: end: #----------------GENERATES SEEDS SEEDGEN:=proc(filename, maxtop, prop_range, tensor_range, dot_range) local n, k, Pnum, Pprop, topol, j, seed, ALLSEED, SEED, alltops, i, l, newseed, NEWS, OLDNEWS, NEWSEED, fp, pp, pj, oldseed, a_min, min_N, N, min_a, a, min_b, b, mmm; min_N:=prop_range[1]: N:=prop_range[2]: min_a:=tensor_range[1]: a:=tensor_range[2]: min_b:=dot_range[1]: b:=dot_range[2]: fp:=fopen(filename, WRITE); fprintf(fp, "Seeds for the %a topology.\n", maxtop); fprintf(fp, "Tensor depth= %a-%a \n",min_a, a); fprintf(fp, "Power depth= %a-%a \n",min_b, b); SEED:=[seq(0, mmm=1..NPROPS)]: for n from min_N to N do alltops:=choose(maxtop, n); for i from 1 to nops(alltops) do topol:=alltops[i]; OLDNEWS:={}: for Pprop from 0 to b do NEWS:=composition(Pprop+n, n); NEWS:=NEWS minus OLDNEWS: # print(NEWS); for pp from 1 to nops(NEWS) do NEWSEED:=SEED: for pj from 1 to nops(topol) do NEWSEED:=subsop(topol[pj]=op(pj, op(pp, NEWS)), NEWSEED); od: # print(NEWSEED); oldseed:={}: for Pnum from 0 to abs(a) do newseed:=map(proc(x) local trelos; trelos:=-(x-[seq(1, mmm=1..NPROPS)]); for j from 1 to nops(topol) do trelos:=subsop(topol[j]=0, trelos): od: RETURN(trelos); end, composition(Pnum+NPROPS, NPROPS) ): newseed:=newseed minus oldseed: #print("newseed", newseed); for l from 1 to nops(newseed) do seed:=NEWSEED+op(l, newseed); if (TOPOLOGIES[1])(op(seed))<>0 and abs(sumnp(xb(op(seed))))- min_a >=0 and abs(sumpp(xb(op(seed))))- min_b >=0 then fprintf(fp, "\n %a", seed); fi: od: oldseed:=oldseed union newseed: od: od: OLDNEWS:=OLDNEWS union NEWS: od: od: od: #fprintf(fp, "NULL ];"); fclose(fp); # RETURN(ALLSEED); end: #----------------------> MAIN ROUTINE <---------- Reducer:=proc(SEED_FILE, CALC_SEED, RESULTS_DIR) #<-change local fp_seed, seed, counter, eqtns, eqs, xxxx, i, fake1, fake2; counter:=0: # Read in the seed fp_seed:=readline(SEED_FILE); while fp_seed<>0 do seed:=traperror(parse(fp_seed)); if type(seed, list) then counter:=counter+1: if (TOPOLOGIES[1])(op(seed))=0 then fp_seed:=readline(SEED_FILE); xxxx:=fopen(CALC_SEED, APPEND); fprintf(xxxx,"%a,%a\n", seed, counter); fflush(xxxx); fclose(xxxx); next; fi; #Create all equations eqtns:=geneq(op(seed)); for i from 1 to nops(eqtns) do # Pick the i equation eqs:=eqtns[i]: # Substitute all known stuff eqs:=sub_eqs(eqs, RESULTS_DIR); eqs:=freeze_k(eqs, MASTERS, RESULTS_DIR); #print(eqs); eqs:=freeze_q(eqs, check_values, RESULTS_DIR); # Equation solver eqs:=eqs_solver(eqs, RESULTS_DIR); if eqs=NEXTT then next; fi: #Store eqs store_eqs(eqs, RESULTS_DIR); #Change information about dependencies mod_incl(eqs, RESULTS_DIR); #Make the backsubstitution back_sub(eqs, RESULTS_DIR); od: xxxx:=fopen(CALC_SEED, APPEND); fprintf(xxxx,"%a %a\n", seed, counter); fflush(xxxx); fclose(xxxx); if VERBOSE=TRUE then print(seed," has been treated.");fi; fi: fp_seed:=readline(SEED_FILE); od: end: #---------------------- routine for back substitution back_sub:=proc(eqs, RESULTS_DIR) local i, lefteq, leftincl, eqslist, filei, old_eqs, new_eqs, INT_SUB, funcss; lefteq:=lhs(eqs); leftincl:=locate_int_inc(lefteq, RESULTS_DIR); if touchfile(leftincl)=true then eqslist:=parsefile(leftincl); for i in eqslist do #read old equation in, and make the substitution #print(i); filei:=locate_int(i, RESULTS_DIR); old_eqs:=parsefile(filei); fclose(filei); # new_eqs:=lhs(old_eqs)=subs(eqs, rhs(old_eqs)); new_eqs:=lhs(old_eqs)=sub_eqs(rhs(old_eqs), RESULTS_DIR); funcss:=map(x->op(0, x), indets(new_eqs, function)): new_eqs:=lhs(new_eqs)=freeze_k(Mycollect(rhs(new_eqs), funcss, factor),MASTERS, RESULTS_DIR); #print(new_eqs); new_eqs:=lhs(new_eqs)=freeze_q(Mycollect(rhs(new_eqs), funcss, factor),check_values, RESULTS_DIR); #store the new equation store_eqs(new_eqs, RESULTS_DIR); #update the include files update_inc(new_eqs, old_eqs, RESULTS_DIR); od: fi: end: #---------------------- side routine for the back-substitution process update_inc:=proc(neweq, oldeq, RESULTS_DIR) local i, target, oldinds, newinds, addtarget, rmtarget, fp, filedir,incdir, incfile, INT_INCL; #updates the include files of the integrals in its arguments if lhs(neweq)<>lhs(oldeq) then ERROR("Wrong kind of arguments"); fi: target:=lhs(neweq); oldinds:=indets(rhs(oldeq), specfunc(anything, {op(TOPOLOGIES)})); newinds:=indets(rhs(neweq), specfunc(anything, {op(TOPOLOGIES)})); addtarget:=newinds minus oldinds; rmtarget:=oldinds minus newinds; for i in addtarget do filedir:=locate_dir(i, RESULTS_DIR); incdir:=locate_dir_inc(i, RESULTS_DIR); incfile:=locate_int_inc(i, RESULTS_DIR); traperror(mkdir(filedir)); traperror(mkdir(incdir)); if touchfile(incfile)=true then INT_INCL:=parsefile(incfile); fclose(incfile); else INT_INCL:={}: fi: INT_INCL:=INT_INCL union {target}; fp:=fopen(incfile, WRITE); fprintf(fp,"%a\n", INT_INCL); fflush(fp); fclose(fp); od: for i in rmtarget do filedir:=locate_dir(i, RESULTS_DIR); incdir:=locate_dir_inc(i, RESULTS_DIR); incfile:=locate_int_inc(i, RESULTS_DIR); traperror(mkdir(filedir)); traperror(mkdir(incdir)); if touchfile(incfile)=true then INT_INCL:=parsefile(incfile); fclose(incfile); else INT_INCL:={}: fi: INT_INCL:=INT_INCL minus {target}; fp:=fopen(incfile, WRITE); fprintf(fp,"%a\n", INT_INCL); fflush(fp); fclose(fp); od: end: mod_incl:=proc(eqs, RESULTS_DIR) local lefth, righth, i, incfile, incdir, filedir, fp, INT_INCL; #does the actual read/write process in modifying the include files lefth:=lhs(eqs); righth:=indets(rhs(eqs), specfunc(anything, {op(TOPOLOGIES)})); for i in righth do filedir:=locate_dir(i, RESULTS_DIR); incdir:=locate_dir_inc(i, RESULTS_DIR); incfile:=locate_int_inc(i, RESULTS_DIR); traperror(mkdir(filedir)); traperror(mkdir(incdir)); if touchfile(incfile)=true then INT_INCL:=parsefile(incfile); fclose(incfile); else INT_INCL:={}: fi: INT_INCL:=INT_INCL union {lefth}; fp:=fopen(incfile, WRITE); fprintf(fp,"%a\n", INT_INCL); fflush(fp); fclose(fp); od: end: #----------------------- side routines for the MAIN ROUTINE store_eqs:=proc(eqs, RESULTS_DIR) local aux, filename,filedir, fp; #stores the equation in its correct file filedir:=locate_dir(lhs(eqs), RESULTS_DIR); traperror(mkdir(filedir)); filename:=locate_int(lhs(eqs), RESULTS_DIR); fp:=fopen(filename, WRITE); fprintf(fp,"%a = %a\n",lhs(eqs), rhs(eqs)); fflush(fp); fclose(fp); end: sub_eqs:=proc(eqs, RESULTS_DIR) local inds, i,j, symbs, outp, flag, funcs, filename, elem, subi, fp, newaux; # performs the FRONT substitution inds:=indets(eqs, function); subi:=[]: for i from 1 to nops(inds) do elem:=op(i, inds): filename:=locate_int(elem, RESULTS_DIR); if touchfile(filename)=true then subi:=[op(subi), parsefile(filename)]: fclose(filename); fi: od: outp:=subs(op(subi), eqs): inds:=map(x->op(0, x), indets(outp, function)): outp:=Mycollect(outp, inds, xx->simplify(xx, size)): RETURN(outp); end: touchfile:=proc(filename) local fp; #checks if a file exists without throwing an error if it doesn't fp:=traperror(readline(filename)); traperror(fclose(filename)); if fp= "file or directory does not exist" then RETURN(false); else RETURN(true); fi: end: put_dash:=proc() local i; #puts a dash in between its arguments, useful when dealing with the directories of integrals #Then one needs strings like DX_1_2_1_-2_1_0_-1 RETURN( cat(seq(cat("_", args[i]) ,i=1..nargs)) ); end: locate_int:=proc(int, RESULTS_DIR) local topol, aux, outp; #finds the adress of an integral aux:=traperror(find_top(int)); outp:=cat(RESULTS_DIR,"/", op(0, int), put_dash(op(aux)),"/", op(0, int) , put_dash(op(int)), ".map"); RETURN(outp); end: locate_dir:=proc(int, RESULTS_DIR) local topol, aux, outp; #finds the directory of an integral aux:=traperror(find_top(int)); outp:=cat(RESULTS_DIR, "/", op(0, int), put_dash(op(aux))); RETURN(outp); end: locate_int_inc:=proc(int, RESULTS_DIR) local topol, aux, outp; #finds the adress of the include file of an integral aux:=traperror(find_top(int)); outp:=cat(RESULTS_DIR,"/", op(0, int), put_dash(op(aux)),"/", "INCLUDE","/", op(0, int) , put_dash(op(int)), ".inc"); RETURN(outp); end: locate_dir_inc:=proc(int, RESULTS_DIR) local topol, aux, outp; #finds the directory of the include file of an integral aux:=traperror(find_top(int)); outp:=cat(RESULTS_DIR, "/", op(0, int), put_dash(op(aux)), "/INCLUDE"); RETURN(outp); end: find_top:=proc() local topol, x, i; #returns the topology of an integral, i.e. a list with the positions of propagators with positive exponents. if whattype(args)=exprseq then x:=[]; for i from 1 to nargs do if args[i]>0 then x:=[op(x), i]; fi; od: RETURN(x); fi; if type(args, {list, function}) then RETURN(procname(op(args))); fi; #print(args); ERROR(`Wrong type of aguments-EGW`); end: eqs_solver:=proc(eqs) local maxprior, outp, funcs, coeffi, neweq, fake1, fake2; #solves the ibp equation for the priority integral (which is returned by the routine priority(...)) global NEXTT; fake1:=0=0; fake2:=0; if eqs=fake1 or eqs=fake2 then RETURN(NEXTT); fi: maxprior:=priority(eqs); #maxprior is the priority integral if maxprior = [] then RETURN(NEXTT): fi: outp:=solve(eqs, maxprior); funcs:=map(x->op(0, x), indets(outp, function)); RETURN(maxprior=Mycollect(outp, funcs, factor)); end: parsefile:=proc(filename) local aaa, str0, str; #parsefile is parses a file and returns the output as a string str0:=" "; str:=str0; aaa:=readline(filename); while aaa<>0 do str:=cat(str, aaa); aaa:=readline(filename); od: RETURN(parse(str)); end: #----------------- routines for Freezing (or Masking) of master integrals and large coefficients freeze_k:=proc(eqs, MASTERS, RESULTS_DIR) local i, activ, stor, work, known, subi0, Ks, counterfile, counter, outp,exprfile, fp, funcs, nonKs; #freezes the master integrals with their coefficients if MASTERS=[] then RETURN(eqs); fi: funcs:=indets(eqs, function): #funcs is a set with K's and integrals Ks:=map(x->if op(0, x)=K then RETURN(x); fi, funcs): nonKs:= map(x->if op(0, x)<>K then RETURN(x); fi, funcs): known:={op(MASTERS), op(Ks)}: subi0:=map(x->x=0, known): work:=Mycollect(eqs, known): activ:=Mycollect(subs(subi0, work), nonKs, xxx->simplify(xxx, size)): stor:=Mycollect(work-activ, {op(nonKs), op(known)}, xxx->simplify(xxx, size)): if stor=0 then RETURN(activ); fi: counterfile:=cat(RESULTS_DIR, "/ICED/kcounter.map"): if touchfile(counterfile)=true then counter:=parsefile(counterfile); else counter:=0: traperror(mkdir(cat(RESULTS_DIR,"/ICED"))): fi: fp:=fopen(counterfile,WRITE): counter:=counter+1: fprintf(fp,"%a",counter): fclose(fp): outp:=activ+K(counter): exprfile:=cat(choose_iced(RESULTS_DIR, counter, ON), "/KEXPR/"): if touchfile(exprfile)=false then traperror(mkdir(exprfile)); fi: exprfile:=cat(choose_iced(RESULTS_DIR, counter), "/KEXPR/kexpr_", counter, ".map"): fp:=fopen(exprfile, WRITE); fprintf(fp, "%a = %a", K(counter), stor); fclose(fp); RETURN(outp); end: freeze_q:=proc(eqs, check_values, RESULTS_DIR) local funcs, knownq, inte, coeffi, final_eqs, allvalues, counter, counterfile, exprfile, numfile, fp, outp, i, coeff_values,gg; #freezes the large coefficients of non-master integrals after checking for possible zero coefficients funcs:=convert(indets(eqs, function), list): #funcs is the list of integrals DX if check_values=[] then RETURN(eqs); fi: #the routine is skipped if there are no checkvalues provided allvalues:=[]: for i from 1 to nops(funcs) do inte:=funcs[i]; #inte is the current integral coeffi:=coeff(eqs, inte): #coeffi is the coefficient of this integral knownq:=map(x->op(1, x), convert(indets(coeffi, indexed), list)): #knownq is a list with the indices of the known q in coeffi coeff_values:=[seq(subs(check_values[j],map(x->findq(x, j, RESULTS_DIR), knownq),coeffi),j=1..nops(check_values))]: #coeff_values is a list of the numerical values of the coeffi for every check value set #of the kinematical parameters. allvalues:=[op(allvalues), [inte, coeffi, simplify(coeff_values)]]: #allvalues is a list of [integral,coefficient,coefficientvalues] od: allvalues:=map(proc(x) local yy; yy:=sum(simplify(abs(x[3][kk])), kk=1..nops(check_values)): for gg from 1 to nops(check_values) do if (abs(x[3][gg])=0 and yy>0) then print("WARNING:possible error. A non-zero coefficient has an exceptional zero check_value. If you get DIVISION BY ZERO, please change the check_values");fi;od; if yy>0 then RETURN(x); fi; RETURN(NULL); end ,allvalues): #allvalues elements are mapped to themselves if all of their values are nonzero. Otherwise they are eliminated from allvalues counterfile:=cat(RESULTS_DIR, "/ICED/counter.map"): if touchfile(counterfile)=true then counter:=parsefile(counterfile); else counter:=0: traperror(mkdir(cat(RESULTS_DIR,"/ICED"))): fi: #print(allvalues); outp:=0: for i from 1 to nops(allvalues) do if length(allvalues[i][2]) < MAXLENGTH then #if the length of the coefficient is smaller than MAXLENGTH outp:=outp+allvalues[i][2]*allvalues[i][1]: #then it is not frozen else fp:=fopen(counterfile,WRITE): counter:=counter+1: fprintf(fp,"%a",counter): fclose(fp): outp:=outp+q[counter]*allvalues[i][1]: #else it is frozen exprfile:=cat(choose_iced(RESULTS_DIR, counter, ON), "/EXPR/"): if touchfile(exprfile)=false then traperror(mkdir(exprfile)); fi: exprfile:=cat(choose_iced(RESULTS_DIR, counter), "/EXPR/qexpr_", counter, ".map"): fp:=fopen(exprfile, WRITE); fprintf(fp, "%a = %a", q[counter], allvalues[i][2]); fclose(fp); exprfile:=cat(choose_iced(RESULTS_DIR, counter), "/NUM/"): if touchfile(exprfile)=false then traperror(mkdir(exprfile)); fi: exprfile:=cat(choose_iced(RESULTS_DIR, counter) ,"/NUM/qnum_", counter, ".map"): fp:=fopen(exprfile, WRITE); fprintf(fp, "%a",allvalues[i][3]); fclose(fp); fi: od: RETURN(outp); end: #------------------------- file r/w routines for freezing findq:=proc(num, j, RESULTS_DIR) local filename, vals; #finds the numerical values of a q filename:=cat(choose_iced(RESULTS_DIR, num), "/NUM/qnum_", num, ".map"): vals:=parsefile(filename): RETURN(q[num]=vals[j]); end: choose_iced:=proc(RESULTS_DIR, counter) local filename; #finds the directory where the next q will be saved filename:=cat(RESULTS_DIR, "/", convert(ICED, string), "/", convert(ICED, string), floor(counter/MAXDIR) ); if nargs=3 and args[3]=ON then mymkdir(filename): fi: RETURN(filename); end: mymkdir:=proc(dirname) local restname, currname, pos; #alternate version of mkdir for making a directory restname:=dirname: currname:="": while length(restname) >0 do pos:=searchtext("/", restname): if pos = 0 then pos:= length(restname): fi: currname:=cat(currname, substring(restname, 1..pos)); restname:=substring(restname, pos+1..length(restname)); if touchfile(currname)=false then traperror(mkdir(currname));fi: od: end: mysimplify:=proc() #alternative for simplify (simpler and more efficient for our purposes) if length(args[1]) > MAXSIMPLIFY then RETURN(args[1]); fi: RETURN(simplify(args)); end: #----------------------- routine for Intermediate freezing during melting #during the melting procedure q's and k's are melt in nested combination. Then through this routine #the large symbolic expressions are frozen again as f's. These f's will be finally melt in the last #stage of the program freeze_f:=proc(expr, RESULTS_DIR) local counterfile, counter, lll, filename, num, fp; lll:=length(expr): #lll is the length of the expression if lll < MAPLEMAXSUB then RETURN(expr);fi: #if the length is smaller than the GLOBAL MAPLEMAXSUB the expr is returned untouched counterfile:=cat(RESULTS_DIR, "/ICED/fcounter.map"): if touchfile(counterfile)=true then #the counterfile is opened or created and the counter is read or set to 0 counter:=parsefile(counterfile); else counter:=0: traperror(mkdir(cat(RESULTS_DIR,"/ICED"))): fi: num:=counter+1: #num is the new counter filename:=cat(choose_iced(RESULTS_DIR, num), "/FEXPR/"): #the f_file directory is created if touchfile(filename)=false then traperror(mkdir(filename)); fi: filename:=cat(choose_iced(RESULTS_DIR, num), "/FEXPR/fexpr_", num, ".map"): #the f_file is created and opened fp:=fopen(filename, WRITE); fprintf(fp, "%a = %a", f[num], expr); #the expression is written there fclose(fp); fp:=fopen(counterfile, WRITE); fprintf(fp, "%a", num); #the new counter is written fclose(fp); RETURN(f[num]); #the alias f[num] is returned end: #------------------ melting all f-coefficients melt_all_f:=proc(RESULTS_DIR) local i,j,counter,flist,curf,curcurf,fp,filename,curnum; counter:=parsefile(cat(RESULTS_DIR,"/ICED/fcounter.map")); if nargs()=2 then counter:=args[2];fi; for i from 1 to counter do filename:=cat(choose_iced(RESULTS_DIR, i), "/FEXPR/fexpr_", i, ".map"): curf:=rhs(parsefile(filename)); fp:=fopen(filename, WRITE); fprintf(fp, "%a = %a", f[i], simplify(curf,size)); #the expression is simplified and rewritten in its file fclose(fp); od; for i from 1 to counter do filename:=cat(choose_iced(RESULTS_DIR, i), "/FEXPR/fexpr_", i, ".map"): curf:=rhs(parsefile(filename)); flist:=convert(select(has,indets(curf,indexed),f),list); if nops(flist)>0 then for j from 1 to nops(flist) do curnum:=op(1,op(j,flist)); curcurf:=parsefile(cat(choose_iced(RESULTS_DIR, curnum), "/FEXPR/fexpr_", curnum, ".map")); curf:=subs(curcurf,curf); gc(); od; fi; if VERBOSE=TRUE then print(cat("melting" ,i,"/",counter));fi; curf:=normal(curf); flist:=convert(select(has,indets(curf,indexed),f),list); #- - > if nops(flist)>0 then print("error:not all f are melt");fi; fp:=fopen(filename, WRITE); fprintf(fp, "%a = %a", f[i], curf); #the expression is written there fclose(fp); gc(); od; end: #----------------------- routines for melting K-coefficients show_kcounter:=proc(RESULTS_DIR) #returns the number of freezed coefficients. local filename; filename:=cat(RESULTS_DIR, "/ICED/kcounter.map"); if touchfile(filename) = false then RETURN(0); fi; RETURN(parsefile(filename)); end: show_K:=proc(num, RESULTS_DIR) local filename,outp; #returns the K(num) checking both in melt and normal directories outp:=show_mK(num,RESULTS_DIR); if outp<>false then RETURN(outp);fi; filename:=cat(choose_iced(RESULTS_DIR, num), "/KEXPR/kexpr_", num, ".map"); if touchfile(filename)=false then RETURN(K(num)); fi: RETURN(rhs(parsefile(filename))); end: choose_meltK:=proc(RESULTS_DIR, counter) local filename; #returns the adress where the melt k will be saved filename:=cat(RESULTS_DIR, "/", convert(MELTK, string), "/",convert(MELTK, string),floor(counter/MAXDIR)); if nargs=3 and args[3]=ON then mymkdir(filename): fi: RETURN(filename); end: show_mK:=proc(num, RESULTS_DIR) local filename; #returns the melt K if it exists or false otherwise filename:=cat(choose_meltK(RESULTS_DIR, num), "/kmelt_", num, ".map"); if touchfile(filename)=false then RETURN(false); fi: RETURN(rhs(parsefile(filename))); end: melt_K:=proc(num, RESULTS_DIR) #melts K's the top down way local outp, inds, fp, filename, subi, indsq, funcs; outp:=show_mK(num, RESULTS_DIR): if outp <> false then RETURN(outp); fi: #checks for melt K. if so the melt K is returned outp:=show_K(num, RESULTS_DIR); #the K(i) is read inds:=map(x->K(x), sort(map(op, convert(select(has, indets(outp, function), K), list)))); #inds is a list with the K's contained in the specific K if inds=[] then #if there is no K indsq:=sort(map(op, convert(select(has, indets(outp, indexed), q), list))): #indsq is a list with q's in K indsq:=map(xxx->q[xxx]=melt_q(xxx, RESULTS_DIR), indsq): #all the q's are melt in a list outp:=subs(op(indsq), outp): #the q's are substituted in outp funcs:=map(xxx->op(0, xxx), indets(outp, function)): #funcs are the functions contained in outp,i.e. DX,PB or whatever is the name of the integrals # You must refreeze here # outp:=Mycollect(outp, funcs,xxx->freeze_f(normal(xxx), RESULTS_DIR)): #outp is collected with the integral coefficients refrozen by freeze_f filename:=cat(choose_meltK(RESULTS_DIR, num, ON), "/kmelt_", num, ".map"); fp:=fopen(filename, WRITE); fprintf(fp, "K(%a)=%a",num,outp); #the outp is written in MELT_K directory fclose(fp); RETURN(procname(num, RESULTS_DIR)); #the routine is rerun but now the K(i) will be found in the MELT_K dir fi: #if the inds list is not empty, i.e. if there is a nested K indsq:=sort(map(op, convert(select(has, indets(outp, indexed), q), list))): #indsq is a list if q's indsq:=map(xxx->q[xxx]=melt_q(xxx, RESULTS_DIR), indsq): #the q's are melt outp:=subs(op(indsq), outp): #and substituted in outp subi:=map(x->x=eval(melt_K(op(1, x), RESULTS_DIR)), inds): #every nested K is melt by recursively calling the current routine outp:=subs(subi, outp): #and it is substituted in outp funcs:=map(xxx->op(0, xxx), indets(outp, function)): outp:=Mycollect(outp, funcs,xxx->freeze_f(normal(xxx), RESULTS_DIR)): #outp is collected with the coefficient of integrals refrozen by freeze_f filename:=cat(choose_meltK(RESULTS_DIR, num, ON), "/kmelt_", num, ".map"); fp:=fopen(filename, WRITE); fprintf(fp, "K(%a)=%a",num,outp); #the K is written in its directory fclose(fp); RETURN(procname(num, RESULTS_DIR)); #the routine is recalled but now the K is found melt. end: #------------------------- routines for q-melting show_counter:=proc(RESULTS_DIR) local filename; #returns the q-counter from its file filename:=cat(RESULTS_DIR, "/ICED/counter.map"); if (touchfile(filename)=false) then RETURN(0); fi; RETURN(parsefile(filename)); end: show_q:=proc(num, RESULTS_DIR) #returns the saved coefficient local filename,outp; #returns q[num] checking in the melt_q and the normal directory outp:=show_m(num,RESULTS_DIR); if outp<>false then RETURN(outp);fi; filename:=cat(choose_iced(RESULTS_DIR, num), "/EXPR/qexpr_", num, ".map"); if touchfile(filename)=false then RETURN(q[num]); fi: RETURN(rhs(parsefile(filename))); end: choose_melt:=proc(RESULTS_DIR, counter) local filename; #returns the adress where the next melt q should be saved filename:=cat(RESULTS_DIR, "/", convert(MELT, string), "/",convert(MELT, string),floor(counter/MAXDIR)); if nargs=3 and args[3]=ON then mymkdir(filename): fi: RETURN(filename); end: show_m:=proc(num, RESULTS_DIR) local filename; #returns the melt q[num] or false if that q is not melt filename:=cat(choose_melt(RESULTS_DIR, num), "/qmelt_", num, ".map"); if touchfile(filename)=false then RETURN(false); fi: RETURN(rhs(parsefile(filename))); end: melt_q:=proc(num, RESULTS_DIR) local outp, inds, fp, filename, subi; #routine for q-melting. It uses f_freeze i.e. it automatically refreezes the large expressions outp:=show_m(num, RESULTS_DIR): if outp <> false then RETURN(outp); fi: outp:=show_q(num, RESULTS_DIR); inds:=map(x->q[x], sort(map(op, convert(select(has, indets(outp, indexed), q), list)))); if inds=[] then filename:=cat(choose_melt(RESULTS_DIR, num, ON), "/qmelt_", num, ".map"); fp:=fopen(filename, WRITE); fprintf(fp, "q[%a]=%a",num,freeze_f(normal(outp), RESULTS_DIR)); fclose(fp); RETURN(procname(num, RESULTS_DIR)); fi: subi:=map(x->x=eval(melt_q(op(1, x), RESULTS_DIR)), inds): outp:=normal(subs(subi, outp)): filename:=cat(choose_melt(RESULTS_DIR, num, ON), "/qmelt_", num, ".map"); fp:=fopen(filename, WRITE); fprintf(fp, "q[%a]=%a",num,freeze_f(outp, RESULTS_DIR)); fclose(fp); RETURN(procname(num, RESULTS_DIR)); end: #----------------- main melting routines for the user show_int:=proc(inte, RESULTS_DIR) local topol, dirloc, fileloc, nam, i, outp, inds, subi, fp; # returns the final form of the wanted integral with f-aliases for large coefficients #it is prefered that the routine tidy_list runs BEFORE show_int if not type(inte, function) then RETURN(inte); fi: # the given integral must be a function of integers #i.e. of the type DX(1,0,0,1,1,0,1,..) for i in inte do if not type(i, integer) then RETURN(inte); fi: od: topol:=traperror(find_top(inte)): #the topology of the integral is found nam:=op(0, inte): dirloc:=cat(RESULTS_DIR, "/RESULTS/", nam, put_dash(op(topol)) ); # the directory of the integral is found in the already melted integrals fileloc:=cat(dirloc, "/", nam, put_dash(op(inte)), ".map"); # the file where the integral is located is found in the already melted integrals if touchfile(fileloc)=true then #this is the case when the demanded integral has already been melt before #and archived in result_dir/RESULTS RETURN(rhs(parsefile(fileloc))); #in this case the integral in just parsede and returned. fi: dirloc:=cat(RESULTS_DIR, "/", nam, put_dash(op(topol)) ); # the directory of the integral is found fileloc:=cat(dirloc, "/", nam, put_dash(op(inte)), ".map"); # the file where the integral is located is found if touchfile(fileloc)=true then outp:=rhs(parsefile(fileloc)); #the integral is parsed from its file. inds:=select(has, indets(outp, indexed), q): #frozen coefficients are located subi:=map( x-> x=melt_q(op(1,x), RESULTS_DIR), inds): #they are melt outp:=subs(subi, outp): #they are substituted in the equation for the integral inds:=select(has, indets(outp, function), K): #master integrals are located subi:=map( x-> x=melt_K(op(1,x), RESULTS_DIR), inds): # they are melt outp:=subs(subi, outp): # they are substituted in the equation for the integral inds:=map(x->op(0, x), indets(outp, function)): outp:=Mycollect(outp, inds, xxx->mysimplify(subs(subi, xxx), size)): dirloc:=cat(RESULTS_DIR, "/RESULTS/", nam, put_dash(op(topol)) ); # the calculated integral is archived in RESULTS mymkdir(dirloc); fileloc:=cat(dirloc, "/", nam, put_dash(op(inte)), ".map"); fp:=fopen(fileloc, WRITE); fprintf(fp, "%a=%a",inte,outp); fclose(fp); RETURN(procname(inte, RESULTS_DIR)); # the routine is run again, but now the integral WILL be found in Results ! fi: RETURN(inte); #this is a safety return : if there is no file for it in Melt directory # or in RESULTS_DIR then the integral is not reduced. Hence it is returned intact. end: #-------------------------> main melting routine <-------------- tidy_list:=proc(SEED_FILE) local i, filename, aux, inte, fp_seed, seed, RESULTS_DIR; # it melts all the integrals that are produced as seeds for the ibp's . It uses refreezing through #the melt_K and melt_q routines above. if nargs=2 then RESULTS_DIR:=args[2]: else RESULTS_DIR:=("."): fi: #if there is a second argument, that's the RESULTS_DIR. otherwise RESULTS_DIR is the current dir fp_seed:=readline(SEED_FILE); while fp_seed<>0 do seed:=traperror(parse(fp_seed)); #the seed is read if type(seed, list) then if VERBOSE=TRUE then print("treating seed ",seed);fi; inte:=(TOPOLOGIES[1])(op(seed)); show_int(inte, RESULTS_DIR); #the seed integral is melt fi: fp_seed:=readline(SEED_FILE); #next seed is read od: end: show_full_int:=proc(int,RESULTS_DIR) local outp,i,fs; #returns the wanted integral fully reduced and without f-aliases outp:=show_int(int,RESULTS_DIR); fs:=convert(select(has,indets(outp,indexed),f),list); for i to nops(fs) do outp:=subs(op(i,fs)=show_f(op(op(i,fs)),RESULTS_DIR),outp); od; return(outp); end: #--------------------------- file handling routines for f-coefficients show_fcounter:=proc(RESULTS_DIR) #returns the number of freezed coefficients. local filename; filename:=cat(RESULTS_DIR, "/ICED/fcounter.map"); if touchfile(filename) = false then RETURN(0); fi; RETURN(parsefile(filename)); end: show_f:=proc(num, RESULTS_DIR) #returns the saved coefficient local filename,outp; filename:=cat(choose_iced(RESULTS_DIR, num), "/FEXPR/fexpr_", num, ".map"); if touchfile(filename)=false then RETURN(f[num]); fi: RETURN(rhs(parsefile(filename))); end: