If you are loading FIESTA from Mathematica, you should either use

SetDirectory[<path to FIESTA>]; Get["FIESTA3.m"];

or

FIESTAPath=<path to FIESTA>; Get[FIESTAPath<>"FIESTA3.m"];.

Do not load FIESTA by simply specifying a full path to it, it won’t work properly.
Now you can call the following commands:

SDEvaluate[{U,F,l},indices,order],

where U and F are the functions from formula 2, l is the number of loops,
indices is the set of indices and order is the required order of ε-expansion.
To avoid manual construction of U and F one can use a build-in function
UF and launch the evaluation as

SDEvaluate[UF[loop_momenta,propagators,subst],indices,order],

where subst is a set of substitutions for external momenta, masses and
other values (to remind: the code performs numerical integration so the functions U and F should not depend on anything external).

*Example:* SDEvaluate[UF[{k},{-k2,-(k+p1)2,-(k+p1+p2)2,-(k+p1+p2+p4)2},
{p1^2 → 0, p2^2 → 0, p4^2 → 0, p1 p2 →-S/2,p2 p4 →-T/2,p1 p4 →(S+T)/2,
S→3,T→ 1}], {1,1,1,1},0]

performs the evaluation of the massless on-shell box diagram.

In the following commands we will only provide the first version of the
syntax (with {U,F,l}). However, in all places this triple can be replaced with the UF generator.
Now if you are aiming to expand a Feynman integral by a small parameter
tt, then you should run

SDExpand[{U,F,l},indices,order,tt,order in tt].

*Example:* SDExpand[UF[{k, l}, {-k2, -(k + p1)2, -(k + p1 + p2)2, -l2,
-(l - k)2, -(l + p1 + p2)2, -(l + p1 + p2 + p4)2}, {p12 -> 0, p22 -> 0, p42 -> 0, p1*p2 -> s/2, p2*p4 -> t/2, p1*p4 -> -(s + t)/2, s -> -1, t -> -tt}], {1, 1, 1, 1, 1, 1, 1} , 0, tt, 0]

The new algorithm presented in this paper can also be used to expand this integral, however one has to provide a regularization variable and shift indices.

*Example:* RegVar=la; SDExpandAsy[UF[{k, l}, {-k2, -(k + p1)2, -(k + p1 + p2)2, -l2,
-(l - k)2, -(l + p1 + p2)2, -(l + p1 + p2 + p4)2}, {p12 -> 0, p22 -> 0, p42 -> 0, p1*p2 -> s/2, p2*p4 -> t/2, p1*p4 -> -(s + t)/2, s -> -1, t -> -tt}], {1+la, 1, 1, 1-la, 1, 1, 1} , 0, tt, 0]

If you use Speer sectors strategy, then you should use SDEvaluateG instead of SDEvaluate. The syntax is SDEvaluateG[graph_info,{U,F,l},indices,order] or SDExpandG[graph_info,{U,F,l},indices,order,tt,tt_order]. The graph information should be of the form {lines, external_vertices}, where lines is a list of pairs of vertices connected by this line. The vertices should be numbered from 1 without skipping numbers. It is also very important to have the order of lines coincide with the order of propagators in the U F input. For example, for the tetrahedral input is the following.

*Example:* SDEvaluateG[{{{1, 2}, {2, 3}, {3, 1}, {4, 2}, {4, 3}, {4, 1}, UF[{k1, k2, k3}, {-k12 + 1, -k22 + 1, -k32 + 1, -(k1 - k2)2 + 1, -(k2 - k3)2 + 1, -(k3 - k4)2 + 1, }, {}], {1, 1, 1, 1, 1, 1, 1, 1}, -1]

There is one command that makes possible to apply sectors decomposition to itegrals different from Feynman integrals: SDEvaluateDirect[var,function,degrees,order,deltas_optional]. Here var stands for the integration variable used in functions (for example, x goes for x[1], x[2], ...), functions is a list of functions and degrees is the list of their exponents. order is the requested order is epsilon, and deltas goes for the list of delta functions attached to the integrand. By default is is empty. For example, {{1,3},{2,4}} goes for the product of Delta[x[1]+x[3]] and Delta[x[2]+x[4]].

*Example:* SDEvaluateDirect[x, {1, x[1] x[2] x[3] + x[1] x[2] x[4] + x[1] x[3] x[4] + x[1] x[2] x[5], x[1] x[3] + x[2] x[3] + x[1] x[4] + x[2] x[4] + x[3] x[4] + x[1] x[5] + x[2] x[5] + x[3] x[5]}, {1, -1 - 2 ep, -1 + 3 ep}, 0, {{1, 2, 3, 4, 5}}]

A similar syntax works for the expansion: SDExpandDirect[var,function,degrees,expand_var,deltas].

*Example:* MathematicaBinary=math; SDEvaluateDirect[x, {1, x[1] x[2] x[3] + x[1] x[2] x[4] + x[1] x[3] x[4] + x[1] x[2] x[5], x[1] x[3] + x[2] x[3] + x[1] x[4] + x[2] x[4] + t (x[3] x[4] + x[1] x[5] + x[2] x[5] + x[3] x[5])}, {1, -1 - 2 ep, -1 + 3 ep}, 0, t, 0, {{1, 2, 3, 4, 5}}]

There is one more way to use FIESTA. An analytic Feynman integral evaluation method suggested by Lee needs to know for which values of space-time dimension d the integral can have poles. FIESTA can locate those values with the use of the following command: SDAnalyze[{U,F,l},indices,dmin,dmax]. Here dmin and dmax are the ends of the interval where poles should be located. This syntax used only algebraic transformations, so the result is exact. However, the program might miss some pole cancellations, so some of the returned values might be “false alerts”.

*Example:* SDAnalyze[UF[{k},{-k2,-(k+p1)2,-(k+p1+p2)2,-(k+p1+p2+p4)2},
{p1^2 → 0, p2^2 → 0, p4^2 → 0, p1 p2 →-S/2,p2 p4 →-T/2,p1 p4 →(S+T)/2,
S→3,T→ 1}], {1,1,1,1},1,8]

Returned answer is {2,4} which means that the integrand has poles for d equal to 2 and 4.

The code has the following options:

**DataPath**: by default FIESTA stores databases in the temp subfolder. However you might wish to direct it elsewhere, especially if the folder with FIESTA is on a network disk. The database should be preferably stored on a fast local disk;**NumberOfSubkernels**: the number of subkernels that Mathematica launches. Might be set equal to the number of cores on your computer but should not exceed the number of licensed subkernels. This setting can speed up the integrand preparation;**NumberOfLinks**: the number of CIntegrate processes that will be launched. The name of this option is left for compatibility with old versions of FIESTA where each CIntegrate process was called by a separate MathLink connection. This setting corresponds to the parallelization during integration;**CubaCores**: the internal parallelization option of the integrators inside CIntegrate. By default the integrator uses one core, but it can be changed with this option. This setting also corresponds to the parallelization during integration. Normally NumberOfLinks is more efficient, but there might be situations when increasing CubaCores leads to better results;**CubaBatch**(0 by default): request the integrand to generate integration points in bunches. Both CubaBatch = 0 and CubaBatch = 1 result in requesting one point per cycle, however the CubaBatch = 1 uses indirect storing of data during evaluation optimized for using processor caches. Bigger values might increase performance a lot, however this increases the usage of RAM (and the value of CubaBatch should not approach the maxeval setting of the current integrator);**STRATEGY**: sector decomposition strategy. By default we use STRATEGY_S (our strategy), but there are also such options as STRATEGY_B (Bogner and Weinzierl), STRATEGY_X (Binoth and Heinrich), STRATEGY_SS (Speer sectors) and STRATEGY_KU, STRATEGY_KU0, STRATEGY_KU2 (Kaneko and Ueda). Among the last three strategies the last variant is the full implementation of the algorithm, the first two are faster but might result in more sectors;**QHullPath**: if you use strategies KU and similar or are using the new SDExpandAsy syntax, you need to provide a correct path to the qhull executable. By default it is set to qhull assuming that the package is installed on your system, however you might provide a specific path;**CIntegratePath**: by default the integration pool chooses itself the integration binary, however you might provide another path; 16**UsingC**: by default this option is set to True. This means that FIESTA uses the c++ integration. If set to False, it switches to Mathematica integration, however this is not recommended. With UsingC set to False the option ExactIntegrationOrder (if set) specifies the order in epsilon where FIESTA tries exact integration for some time. The default time is 10 seconds per sector and can be modified by the ExactIntegrationTimeout option;**UsingQLink**: by default this option is set to True. Switching it off will turn off database usage, however in FIESTA3 this is possible only together with UsingC=False;**OnlyPrepare**: by default this option is set to False. In this case the calculation is performed completely, otherwise a database is prepared for integration and a shell command to run CIntegrate without Mathematica is printed. This can be used if you are preparing a database on one computer and are integrating elsewhere, or if you are willing to try different integrators or precision requests;**SeparateTerms**: if True, the integrator receives integrable terms independently, not whole expressions for each sector; on one hand, this simplifies the integrands and the integrators return better precision, on the other hand the error grows after summing up the results, so we can’t recommend whether to use this option or not;**ComplexMode**(False by default): with this setting set to True FIESTA performs a contour deformation in order to avoid poles in physical regions. The deformation size depends on a parameter, that is either set manually by giving a value to the ComplexShift variable or is tuned automatically in the interval from zero to MaxComplexShift (1 by default). Increasing the option LambdaSplit (4 by default) might result in better tuning but slows the preparation; same is true for the search option LambdaIterations that is set by default to 1000; The contour transformation cannot deal with cases when F turns to zero for the reason that some variables are equal to 1. In order to trace those cases, set TestF1=True. In order to handle those singularities, set the BisectionVariables equal to the list of variables such that the integration cube is divided into two parts. By default, the separation point is equal to 1/2, however this can be changed either by setting the BisectionPoint value, or by providing a list of BisectionPoints; setting ConservativeComplexShift to True makes the Lambda estimate more safe;**CurrentIntegrator**: (vegasCuba by default): the integrator used at the final stage. The options allowed in the current setup are vegasCuba, suaveCuba, divonneCuba and cuhreCuba. It is also possible to add your own integrators by modifying the integrators.c source file, however this is far beyond the standart usage of FIESTA. Related parameter (CurrentIntegratorOptions) presents a list of options of the currently chosen integrator (for details see [31]). By default it has no value and the actual options are printed out when you start the evaluation (the defaults are stored within the c++ part). The most commonly used integrator option is maxeval, which can be set, for example, by CurrentIntegratorOptions = {{"maxeval","500000"}}. The default value is 50000. One more virtual integrator is justEvaluate. It simply evaluates the integrand at a given point, by default it is the point where all coordinates are equal to 1/2. Its options are x1, x2 and so on representing the integration coordinates;**SmallX, MPThreshold, PrecisionShift, MPPrecision, MPMin**: the options that allow to fine-tune the MPFR subsystem. For details see the paper on FIESTA2;**RegVar**: an option introduced for SDExpandAsy but usable also in other situations. Sets an extra regularization variable;**AnalyticIntegration**: an option used only for SDExpandAsy, True by default, tells FIESTA to try to take some integrations analytically after introducing regions;**FastASY**: (False by default) specifies the region search mode (used in SDExpandAsy). With FastASY set to False the polynomial U × F is analyzed, with True — the F polynomial. The FastASY variant might work significantly faster and will produce correct results almost all the time, but use it at your own risk;**PolesMultiplicity**: an option used only for SDAnalyze, False by default, changes the answer so that it returns not only values of d but also maximal pole multiplicities;**MathematicaBinary**: a path to the executable Mathematica kernel. If set, it is passed to the integration pool, so it can request Mathematica for values of functions it cannot evaluate itself (currently this feature is used only for PolyGamma);**BucketSize**(25 by default): an option tuning the database usage (for details see the documentation on KyotoCabinet); increasing this variable might result is faster database access, but increases the RAM usage;**MixSectors**(0 by default): lets FIESTA to sum up integrands in different sectors before integration;**RemoveDatabases**(True by default): specifies whether the database files should be removed after the integration;**d0**(4 by default): specifies the space-time dimension;**ReturnErrorWithBrackets**(False by default) changes the output — with True the error estimates are printed as pm[NUMBER] instead of pmNUMBER;**FixSectors**: (True by default) — for the reason of fixing sector numbers and easier debugging we perform the variable substitutions stage on the main kernel. If you set this option to false, this stage will be also made parallel, however it normally does not influence the total time much;**PrimarySectorCoefficients**: you might specify the list of coefficients before primary sectors. A zero means that this sector will be ignored. With this setting you can split your problem into parts and also take diagram symmetries into account.**NoDatabaseLock**: prevents FIESTA from locking the database. This may be away to avoid restrictions on some file systems but might result in corrupted databases.