diff --git a/README.md b/README.md
index 5212721..ddc6da4 100644
--- a/README.md
+++ b/README.md
@@ -3,10 +3,16 @@
 Automatical Generator of Conformal Bootstrap Equation
 
 For more information, see
-[autoboot: A generator of bootstrap equations with global symmetry](https://arxiv.org/abs/1903.10522).
+[autoboot: A generator of bootstrap equations with global symmetry](https://arxiv.org/abs/1903.10522),
+or [An Automated Generation of Bootstrap Equations for Numerical Study of Critical Phenomena](https://arxiv.org/abs/2006.04173).
 
 Some usages also can be checked by typing `?someSymbolName` (for example, `?getGroup`) in [Mathematica](http://reference.wolfram.com/language/tutorial/GettingInformationAboutWolframLanguageObjects.html).
 ``?somePackageName`*`` (for example, ``?ClebschGordan`*``) will give package-lebel information.
+These usages are also documented in [inv.md](doc/inv.md), [group.md](doc/group.md).
+
+- [Setup](#setup)
+- [Usage](#usage)
+- [Example](#example)
 
 ## Setup
 
diff --git a/doc/CustomGroup.md b/doc/CustomGroup.md
index caf06a2..210bf19 100644
--- a/doc/CustomGroup.md
+++ b/doc/CustomGroup.md
@@ -6,6 +6,14 @@ You also attach labels of irrep-objects for all (unique) irreps.
 [`getO[3]`](https://github.com/selpoG/autoboot/tree/master/grouplie.m) is
 a good example for creating custom group-objects.
 
+- [`isrep`](#isrep)
+- [`minrep`](#minrep)
+- [`dim`](#dim)
+- [`dual`](#dual)
+- [`id`](#id)
+- [`prod`](#prod)
+- [`gG` and `gA`](#gg-and-ga)
+
 ## `isrep`
 
 `g[isrep[r]]` must give `True` or `False` (`True` if and only if `r` is a valid irrep-object for `g`).
diff --git a/doc/IrrepLabels.md b/doc/IrrepLabels.md
index 6ee659e..ffd5883 100644
--- a/doc/IrrepLabels.md
+++ b/doc/IrrepLabels.md
@@ -1,5 +1,17 @@
 # Labels for Irreps
 
+- [`G=group[g,i]`](#ggroupgi)
+- [`G=so[2]`](#gso2)
+- [`G=o[2]`](#go2)
+- [`G=so[3]`](#gso3)
+- [`G=o[3]`](#go3)
+- [`G=su[2]`](#gsu2)
+- [`G=su[4]`](#gsu4)
+- [`G=dih[N]` (`N`: even)](#gdihn-n-even)
+- [`G=dih[N]` (`N`: odd)](#gdihn-n-odd)
+- [`G=dic[N]`](#gdicn)
+- [`G=pGroup[G1,G2]`](#gpgroupg1g2)
+
 All labels of irreps are defined as follows:
 
 ## `G=group[g,i]`
diff --git a/doc/group.md b/doc/group.md
new file mode 100644
index 0000000..2c9ed99
--- /dev/null
+++ b/doc/group.md
@@ -0,0 +1,212 @@
+# `group.m`
+
+Supported groups are direct products of finite groups and some Lie groups.
+Supported finite groups are those whose irreps was calculated by `GAP` (in `sgd/`),
+and any dihedral and quartenion groups.
+Supported Lie groups are `su[2]`, `su[4]`, `so[2]`, `o[2]`, `so[3]`, `o[3]`.
+We support only compact groups, so we can assume any finite dimensional irrep can be unitarized.
+
+This package imports `groupd.m` and `grouplie.m`.
+
+- [Get Groups](#get-groups)
+- [Group Data](#group-data)
+- [Irrep-Objects](#irrep-objects)
+- [`groupd.m`](#groupdm)
+- [`grouplie.m`](#groupliem)
+
+---
+
+## Get Groups
+
+### `getGroup`
+
+`getGroup[g,i]` loads data from `sgd/sg.g.i.m` and returns group-object `group[g,i]`.
+`g` is the order of the finite group, `i` is the number of the group assined by `GAP`.
+
+### `product`
+
+`product[g1,g2]` returns group-object `pGroup[g1,g2]`
+which represents direct product of two group-object `g1`, `g2`.
+
+### `group`
+
+`group[g,i]` is a group-object whose order is `g` and whose number assigned by `GAP` is `i`.
+Before using this value, you have to call `getGroup[g,i]` to get proper group-object.
+
+### `pGroup`
+
+`pGroup[g1,g2]` is a group-object which is a direct product of `g1`, `g2`.
+Before using this value, you have to call `product[g1,g2]` to get proper group-object.
+
+### `setGroup`
+
+`setGroup[G]` loads `inv.m` with global symmetry `G`.
+This action clears all values calculated by `inv.m` previously.
+
+### `available`
+
+`available[g,i]` gives whether `group[g,i]` are supported or not.
+
+### `setPrecision` (only in `ngroup.md`, not in `group.md`)
+
+After calling `setPrecision[prec]`, all calculation in this package will be done in precision `prec`
+and any number less than `1/10^(prec-10)` will be choped.
+
+It is assumed that `prec` is sufficiently bigger than `10` and `setPrecision` is called just once just after loading this package.
+
+---
+
+## Group Data
+
+A group-object `g` has attributes `ncg`, `ct`, `id`, `dim`, `prod`, `dual`, `isrep`, `gG`, `gA`, `minrep`.
+You can evaluate attributes in putting it in `g[...]`.
+For example, `g[dim[r]]` gives the dimension of irrep `r`.
+
+### `ncg`
+
+`ncg` is the number of conjugacy classes, which is also the number of inequivalent irreps.
+This is not defined for Lie groups.
+
+### `ct`
+
+`ct` is the character table. This is not defined for Lie groups.
+
+### `id`
+
+`id` is the trivial representation.
+
+### `dim`
+
+`dim[r]` is the dimension of irrep `r`.
+
+### `prod`
+
+`prod[r,s]` gives a list of all irreps arising
+in irreducible decomposition of direct product representation of `r` and `s`.
+`prod[r,s]` may not be duplicate-free.
+
+### `dual`
+
+`dual[r]` gives dual representation of irrep `r`.
+
+### `isrep`
+
+`isrep[r]` gives whether `r` is recognised as a irrep-object of the group-object or not.
+
+### `gG`
+
+`gG` is a list of all generator-objects of finite group part of the group-object.
+
+### `gA`
+
+`gA` is a list of all generator-objects of Lie algebra part of the group-object.
+
+### `minrep`
+
+`minrep[r,s]` gives `r` if `r < s` else `s`. `r` and `s` are irrep-objects.
+
+---
+
+## Irrep-Objects
+
+We need all irreps to be sorted in some linear order.
+All irrep-objects of `G=group[g,i]` are `rep[1]`, `rep[2]`, ..., `rep[n]` (`n=G[ncg]`).
+all irrep-objects of `pGroup[g1,g2]` are `rep[r1,r2]`,
+where `r1` is a irrep-object of `g1` and `r2` is a irrep-object of `g2`.
+`minrep` compares irreps in lexical order.
+
+### `rep`
+
+`rep[n]` is `n`-th irrep-object (`n` is assined by `GAP` and corresponds to the index of `ct`).
+This is recognised only by `group[g,i]`.
+
+`rep[r1,r2]` is natural irrep-object of `pGroup[g1,g2]`
+where `r1` is irrep-object of `g1`, `r2` is irrep-object of `g2`.
+This is recognised only by `pGroup[g1,g2]`.
+
+### `v`
+
+`v[n]` is spin-`n` irrep-object.
+This is recognised only by `dih[n]`, `dic[n]`, `su[2]`, `so[3]`, `o[2]` and `so[2]`.
+
+`v[n,s]` is spin-`n` irrep-object with sign `s`. This is recognised only by `o[3]`.
+
+### `i`
+
+`i[a]` is one-dimensional irrep-object with sign `a`.
+This is recognised only by `dih[n]` (`n`: odd) and `o[2]`.
+
+`i[a,b]` is one-dimensional irrep-object with sign `a,b`.
+This is recognised only by `dih[n]` (`n`:even), `dic[n]`.
+
+---
+
+## `groupd.m`
+
+If `n` is even, all irrep-objects of `dih[n]` are `i[1,1], i[1,-1], i[-1,1], i[-1,-1], v[1], ..., v[n/2-1]`.
+If `n` is odd, all irrep-objects of `dih[n]` are `i[1], i[-1] v[1], ..., v[(n-1)/2]`.
+All irrep-objects of `dic[n]` are `i[1,1], i[1,-1], i[-1,1], i[-1,-1], v[1], ..., v[n-1]`.
+
+### `getDihedral`
+
+`getDihedral[n]` returns group-object `dih[n]` which represents the dihedral group of order `2n`.
+
+### `getDicyclic`
+
+`getDicyclic[n]` returns group-object `dic[n]` which represents the dicyclic group of order `4n`.
+
+### `dih`
+
+`dih[n]` is a group-object which is the dihedral group of order `2n`.
+Before using this value, you have to call `getDihedral[n]` to get proper group-object.
+
+### `dic`
+
+`dic[n]` is a group-object which is the dicyclic group of order `4n`.
+Before using this value, you have to call `getDicyclic[n]` to get proper group-object.
+
+---
+
+## `grouplie.m`
+
+All irrep-objects of `G=su[2]` are `v[0], v[1/2], v[1], v[3/2], ...`.
+
+All irrep-objects of `G=su[4]` are `v[n, m, l]` (`n,m,l=0,1,2,...` and `n >= m >= l`).
+
+All irrep-objects of `G=o[3]` are `v[0,1], v[0,-1], v[1,1], v[1,-1], v[2,1], v[2,-1], v[3,1], v[3,-1], ...`.
+
+All irrep-objects of `G=so[3]` are `v[0], v[1], v[2], v[3], ...`.
+
+All irrep-objects of `G=o[2]` are `i[1], i[-1], v[1], v[2], v[3], ...`.
+
+All irrep-objects of `G=so[2]` are `v[x]` (`x \in \mathbb{R}`).
+
+### `getSU`
+
+`getSU[n]` returns group-object `su[n]` which represents the special unitary group of rank `n`.
+`n` must be `2,4`.
+
+### `getO`
+
+`getO[n]` returns group-object `o[n]` which represents the orthogonal group of rank `n`.
+`n` must be `2,3`.
+
+### `getSO`
+
+`getSO[n]` returns group-object `su[n]` which represents the special orthogonal group of rank `n`.
+`n` must be `2,3`.
+
+### `su`
+
+`su[n]` is a group-object which is the special unitary group of rank `n`.
+Before using this value, you have to call `getSU[n]` to get proper group-object.
+
+### `o`
+
+`o[n]` is a group-object which is the orthogonal group of rank `n`.
+Before using this value, you have to call `getO[n]` to get proper group-object.
+
+### `so`
+
+`so[n]` is a group-object which is the special orthogonal group of rank `n`.
+Before using this value, you have to call `getSO[n]` to get proper group-object.
diff --git a/doc/group.pdf b/doc/group.pdf
new file mode 100644
index 0000000..dcd26cf
Binary files /dev/null and b/doc/group.pdf differ
diff --git a/doc/inv.md b/doc/inv.md
new file mode 100644
index 0000000..2030ca9
--- /dev/null
+++ b/doc/inv.md
@@ -0,0 +1,281 @@
+# inv.m
+
+Do not import this package directly. `inv.m` memoise many values automatically and is imported or cleared by `group.m`.
+
+- [Symmetry](#symmetry)
+- [Clebsch-Gordan Coefficients](#clebsch-gordan-coefficients)
+- [Symmetries of CG Coefficients](#symmetries-of-cg-coefficients)
+- [Conformal Blocks](#conformal-blocks)
+- [OPE Coefficients](#ope-coefficients)
+- [Bootstrap Equations](#bootstrap-equations)
+- [Conversions of Bootstrap Equations](#conversions-of-bootstrap-equations)
+
+---
+
+## Symmetry
+
+### `symmetryGroup`
+
+`symmetryGroup` represents the global symmetry of CFT.
+
+### `clearCG`
+
+`clearCG[]` clears all values calculated by this package.
+
+---
+
+## Clebsch-Gordan Coefficients
+
+### irrep-object
+
+An irrep-object `r` is a symbol such that `symmetryGroup[isrep[r]]` is True.
+We use `r,s,t,r1,r2,...` to denote irrep-objects, and `id` is the irrep-object of the trivial irrep.
+`a,b,c,a1,a2,...` are indices of irreps. An index `a` of irrep `r` must sutisfy `1 <= a <= symmetryGroup[dim[r]]`.
+`n,m` represents multiplicity of Clebsch-Gordan coefficients.
+For more information, please refer to [`arXiv:1903.10522`](https://arxiv.org/abs/1903.10522).
+
+### `inv`
+
+`inv[r,s,t]` gives multiplicity of `id` in decomposition of `r \[TensorProduct] s \[TensorProduct] t`.
+
+`inv[{r,s},{t}]` gives multiplicity of `t` in decomposition of `r \[TensorProduct] s`.
+
+`inv[r1,r2,r3,r4]` gives an association
+whose key is a irrep `s` such that `inv[{r1,r2},{s}]>0` and `inv[r3,r4,s]>0`,
+whose value is a list `{inv[{r1,r2},{s}], inv[r3,r4,s]}`.
+
+### `invs`
+
+`invs[r1,r2,r3,r4]` gives a list of `{s,n,m}` such that `cor[r1,r2,r3,r4][s,n,m]` is defined.
+
+### `ope`
+
+`ope[r,s,t][n][a,b,c]` gives Clebsch-Gordan coefficient defined by
+`|r,a\[RightAngleBracket] \[TensorProduct] |s,b\[RightAngleBracket] = \[Sum] ope[r,s,t][n][a,b,c] |t,c,n\[RightAngleBracket]`.
+
+`ope[r][a,b]` gives `Sqrt[dim[r]] ope[r,dual[r],id][1][a,b,1]`.
+
+### `cor`
+
+`cor[r,s,t][n][a,b,c]` gives `\[Sum] ope[r,s,dual[t]][n][a,b,c2] ope[dual[t]][c2,c] / Sqrt[dim[t]]`.
+
+`cor[r1,r2,r3,r4][s,n,m][a1,a2,a3,a4]` gives`\[Sum] cor[r1,r2,dual[s]][n][a1,a2,b] ope[r3,r4,dual[s]][a3,a4,b].`
+
+---
+
+## Symmetries of CG Coefficients
+
+### `\[Sigma]`
+
+`\[Sigma][r,s,t][n]` gives the sign change in swap between `r` and `s` in `cor`,
+i.e. `cor[r,s,t][n][a,b,c] = \[Sigma][r,s,t][n] cor[s,r,t][n][b,a,c]`.
+
+`\[Sigma][r]` gives `\[Sigma][r,dual[r],id]`.
+
+`\[Sigma][op[x,r,p,q]]` gives `p`.
+
+### `\[Tau]`
+
+`\[Tau][r,s,t][n,m]` describes the behavior of `cor` under cyclic permutation of `r,s,t`,
+i.e. `cor[s,t,r][m][b,c,a] = \[Sum] cor[r,s,t][n][a,b,c] \[Tau][r,s,t][n,m]`.
+
+### `\[Omega]`
+
+`\[Omega][r,s,t][n,m]` describes the behavior of `cor` under complex conjugate,
+i.e. `\[Sum] ope[dual[r]][a2,a] ope[dual[s]][b2,b] ope[dual[t]][c2,c] cor[dual[r],dual[s],dual[t]][m][a2,b2,c2]\[Conjugate] = \[Sum] cor[r,s,t][n][a,b,c] \[Omega][r,s,t][n,m]`.
+
+### `six`
+
+`six[r1,r2,r3,r4][s,n,m,t,k,l]` describes the behavior of `cor` under swap between `r2` and `r4`,
+i.e. `cor[r1,r4,r3,r2][t,k,l][a1,a4,a3,a2] = \[Sum] cor[r1,r2,r3,r4][s,n,m][a1,a2,a3,a4] six[r1,r2,r3,r4][s,n,m,t,k,l]`.
+
+### `isReal`
+
+`isReal[r]` gives whether `r` is a real representation or not.
+
+### `isComplex`
+
+`isComplex[r]` gives whether `r` is a complex representation or not.
+
+### `isPseaudo`
+
+`isPseaudo[r]` gives whether `r` is a pseaudo real representation or not.
+
+---
+
+## Conformal Blocks
+
+### `op`
+
+`op[x,r,p,q]` represents a primary operator object O whose name is `x`, which belongs to irrep `r`,
+which sign (=`\[Sigma][O]`) is `p` and whose parity of spin is `q` (=`±1`).
+If O appears in summation, we use 'op' as its name.
+
+`op[x,r]` is a short-hand for `op[x,r,1,1]`.
+
+`op[x]` is a short-hand for `op[x,r,1,1]` if you registered `op[x,r,1,1]` as a fundamental scalar.
+
+### `dualOp`
+
+`dualOp[op[x,r,p,q]]` gives dual operator object of `op[x,r,p,q]`.
+
+### `format`
+
+`format[eq]` gives readable format of `eq` with little loss of information.
+We need much redundancy to calculate properly, so formatted value cannot be used for any argument of our function.
+`format[eq]` is assumed to be used only for human-readability of last output.
+
+### `sum`
+
+`sum[x,op[op,r,1,q]]` represents sum of `x` over all intermediate primary operator `O=op[op,r,1,q]`
+which belongs to irrep `r` and whose parity of spin is `q`.
+`sum[...]` is automatically expanded so as `x` to be `(F or H)*\[Beta]^2`.
+
+### `single`
+
+`single[x]` represents `x`. We need to wrap `x` for redundancy.
+`single[...]` is automatically expanded so as `x` to be `(Fp or Hp)*\[Beta]^2`, or `(Fp or Hp)`.
+
+### `F`
+
+`F[o1,o2,o3,o4,s,p]` represents generalized conformal block `F_{p,s}^{1,2,3,4}`, where `o1,...,o4` are primary scalars.
+
+`F[a,b,c,d]` represents normal conformal block of type-F.
+
+### `H`
+
+`H[a,b,c,d]` represents normal conformal block of type-H.
+
+### `Fp`
+
+`Fp[o1,o2,o3,o4,o]` represents generalized conformal block `F_{o}^{1,2,3,4}`, where `o1,...,o4` and intermediate `o` are primary scalars.
+
+`Fp[a,b,c,d,o]` represents normal conformal block of type-F with intermediate `o`.
+
+### `Hp`
+
+`Hp[a,b,c,d,o]` represents normal conformal block of type-H with intermediate `o`.
+
+---
+
+## OPE Coefficients
+
+### `\[Lambda]`
+
+`\[Lambda][o1,o2,o3][n]` gives `n`-th OPE coefficient of `o1\[Times]o2 \[Rule] OverBar[o3]`.
+
+### `\[Alpha]`
+
+`\[Alpha][o1,o2,o3][n]` gives `n`-th coefficient of three-point function `\[LeftAngleBracket]o1 o2 o3\[RightAngleBracket]`.
+
+### `\[Mu]`
+
+`\[Mu][o1,o2,o3][n]` gives n-th OPE coefficient of `o1\[Times]o2 \[Rule] OverBar[o3]`, where `o3` is registered as a fundamental scalar.
+
+### `\[Nu]`
+
+`\[Nu][o1,o2,o3][n]` gives `n`-th coefficient of three-point function `\[LeftAngleBracket]o1 o2 o3\[RightAngleBracket]`,
+where `o3` is registered as a fundamental scalar."
+
+### `\[Beta]`
+
+`\[Beta][o1,o2,o3][n]` gives minimal basis to describe `\[Lambda]`, `\[Alpha]`, etc.
+
+`\[Lambda]`, `\[Alpha]`, `\[Mu]` and `\[Nu]` are linear combination of `\[Beta]`.
+
+---
+
+## Bootstrap Equations
+
+### `eqn`
+
+`eqn[{a,b,...}]` represents bootstrap equation that claims `a,b,...` must equals to `0`,
+where `a,b,...` are real-linear combination of sum and single.
+
+`eqn[sec,{a,b,...}]` represents extracted part of `eqn[{...}]` which contains only `sum` or `single` related to `sec`.
+
+### `bootAll`
+
+`bootAll[ops]` generates bootstrap equation from all four-point functions of `ops`.
+
+`bootAll[]` generates bootstrap equation from all four-point functions of fundamental scalars.
+
+### `setOps`
+
+`setOps[ops]` registers `ops` and duals of `ops` as fundamental scalars.
+
+Fundamental scalars are used to seperate sum of conformal blocks over scalars to `single[...]` and `sum[...]`.
+
+### `one`
+
+`one` represents the unit operator. This is implicitly registered as a fundamental scalar.
+
+### `extract`
+
+`extract[x,op[op,r,1,p]]` extracts terms of the form `sum[...,op[op,r,1,p]]` from `x`.
+
+`extract[x,scalar]` extracts terms of the form `single[(Fp or Hp)*\[Beta]^2]` from `x`.
+
+`extract[x,unit]` extracts terms of the form `single[Fp or Hp]` from `x` (contribution of the unit operator).
+
+`extract` is a projection:
+`x == extract[x,unit] + extract[x,scalar] + \[Sum]_{r,p} extract[x,op[op,r,1,p]]` and `extract[x,sec] == extract[extract[x,sec]]`.
+
+### `unit`
+
+`unit` is a option for `extract` and means contribution of unit operator.
+
+### `scalar`
+
+`scalar` is a option for `extract` and means contribution of fundamental scalars but unit operator.
+
+### `sector`
+
+`sector[eq]` gives a list of all nontrivial option `sec` for extract applied to `eq`, i.e. `{sec | extract[eq,sec] != 0}`.
+
+---
+
+## Conversions of Bootstrap Equations
+
+### `makeG`
+
+`makeG[eqn[sec,{a,b,...}]]` gives an undirected graph
+whose vertices are OPE coefficients `\[Beta]` in extracted bootstrap equation `eqn[sec,{a,b,...}]`.
+
+### `makeMat`
+
+`makeMat[eqn[sec,{a,b,...}]]` gives a matrix-representation of extracted bootstrap equation `eqn[sec,{a,b,...}]`.
+
+### `makeSDP`
+
+`makeSDP[eqn[{a,b,...}]]` converts whole bootstrap equation `eqn[{a,b,...}]` into sdp-object.
+
+### `sdpobj`
+
+`sdpobj[secs,scalarnum,vals,mats]` is a sdp-object.
+`secs` is section data of bootstrap equation.
+`scalarnum` is the number of connected components in scalar sections.
+`vals` are real constants in bootstrap equation.
+`mats` are matrix-representation of bootstrap equation.
+
+### `toQboot`
+
+`toQboot[sdp]` converts sdp-object into `c++` code for `qboot`.
+
+### `toCboot`
+
+`toCboot[sdp]` converts sdp-object into `python` code for `cboot`.
+
+### `toTeX`
+
+`toTeX[eq]` gives `LaTeX` string of `eq` (you need call `Print[toTeX[eq]]` to paste to your `tex` file).
+
+`toTeX[eqn[{a,b,...}]]` gives `LaTeX` string of `eq` with `align` environment (you need call `Print[toTeX[eq]]` to paste to your `tex` file).
+
+### `repToTeX`
+
+`repToTeX[r]` is needed to transrate irrep-object `r` as `LaTeX` string. Please set appropriate value.
+
+### `opToTeX`
+
+`opToTeX[o]` is needed to transrate operator name `o` as `LaTeX` string. Please set appropriate value.
diff --git a/doc/inv.pdf b/doc/inv.pdf
new file mode 100644
index 0000000..450321e
Binary files /dev/null and b/doc/inv.pdf differ
diff --git a/group.m b/group.m
index 9ebbba4..b29b334 100644
--- a/group.m
+++ b/group.m
@@ -4,7 +4,7 @@
 
 (* Supported groups are direct products of finite groups and some lie groups. *)
 (* Supported finite groups are those whose irreps was calculated by GAP (in sgd folder), and any dihedral and quartenion groups. *)
-(* Supported lie groups are su[2], so[2], o[2], so[3], o[3]. We support only compact groups, so we can assume any finite dim. irrep can be unitarized. *)
+(* Supported lie groups are su[2], su[4], so[2], o[2], so[3], o[3]. We support only compact groups, so we can assume any finite dim. irrep can be unitarized. *)
 getGroup::usage = "getGroup[g,i] loads data from sg.g.i.m and returns group-object group[g,i]. g is the order of the finite group, i is the number of the group assined by GAP."
 product::usage = "product[g1,g2] returns group-object pGroup[g1,g2] which represents direct product of two group-object g1, g2."
 group::usage = "group[g,i] is a group-object whose order is g and whose number assigned by GAP is i. Before using this value, you have to call getGroup[g,i] to get proper group-object."
diff --git a/groupd.m b/groupd.m
index 2044a8a..ba8a1a7 100644
--- a/groupd.m
+++ b/groupd.m
@@ -4,8 +4,8 @@
 
 getDihedral::usage = "getDihedral[n] returns group-object dih[n] which represents the dihedral group of order 2n."
 getDicyclic::usage = "getDicyclic[n] returns group-object dic[n] which represents the dicyclic group of order 4n."
-dih::usage = "dih[n] is a group-object which is the dihedral group of order 2n. Before using this value, you have to call dih[n] to get proper group-object."
-dic::usage = "dic[n] is a group-object which is the dicyclic group of order 4n. Before using this value, you have to call dic[n] to get proper group-object."
+dih::usage = "dih[n] is a group-object which is the dihedral group of order 2n. Before using this value, you have to call getDihedral[n] to get proper group-object."
+dic::usage = "dic[n] is a group-object which is the dicyclic group of order 4n. Before using this value, you have to call getDicyclic[n] to get proper group-object."
 
 (* if n is even, all irrep-objects of dih[n] are i[1,1], i[1,-1], i[-1,1], i[-1,-1], v[1], ..., v[n/2-1]. *)
 (* if n is odd, all irrep-objects of dih[n] are i[1], i[-1] v[1], ..., v[(n-1)/2]. *)
diff --git a/inv.m b/inv.m
index aa0f4cf..0c80bfc 100644
--- a/inv.m
+++ b/inv.m
@@ -1,6 +1,7 @@
 (* Do not import this package directly. *)
 Needs["CommonFunctions`", "common.m"]
 Needs["ToPython`", "topy.m"]
+Needs["ToCpp`", "tocpp.m"]
 ClebschGordan`clearCG[];
 
 BeginPackage["ClebschGordan`"]
@@ -64,7 +65,7 @@
 
 (* fundamental scalars are used to seperate sum of conformal blocks over scalars to single[...] and sum[...]. *)
 setOps::usage = "setOps[ops] registers ops and duals of ops as fundamental scalars."
-one::usage = "one represents unit operator. This is implicitly registered as a fundamental scalar."
+one::usage = "one represents the unit operator. This is implicitly registered as a fundamental scalar."
 
 (* extract is a projection. x==extract[x,unit]+extract[x,scalar]+\[Sum]_{r,p}extract[x,op[op,r,1,p]] and extract[x,sec]==extract[extract[x,sec]]. *)
 extract::usage = "extract[x,op[op,r,1,p]] extracts terms of the form sum[...,op[op,r,1,p]] from x.
@@ -72,12 +73,13 @@
 extract[x,unit] extracts terms of the form single[Fp or Hp] from x (contribution of unit operator)."
 unit::usage = "unit is a option for extract and means contribution of unit operator."
 scalar::usage = "scalar is a option for extract and means contribution of fundamental scalars but unit operator."
-sector::usage = "sector[eq] gives a list of all nontrivial option sec for extract applied to eq, i.e. {sec | extract[eq,sec]!=0}. "
+sector::usage = "sector[eq] gives a list of all nontrivial option sec for extract applied to eq, i.e. {sec | extract[eq,sec]!=0}."
 
 makeG::usage = "makeG[eqn[sec,{a,b,...}]] gives an undirected graph whose vertices are OPE coefficients \[Beta] in extracted bootstrap equation eqn[sec,{a,b,...}]."
 makeMat::usage = "makeMat[eqn[sec,{a,b,...}]] gives a matrix-representation of extracted bootstrap equation eqn[sec,{a,b,...}]."
 makeSDP::usage = "makeSDP[eqn[{a,b,...}]] converts whole bootstrap equation eqn[{a,b,...}] into sdp-object."
 sdpobj::usage = "sdpobj[secs,scalarnum,vals,mats] is a sdp-object. secs is section data of bootstrap equation. scalarnum is the number of connected components in scalar sections. vals are real constants in bootstrap equation. mats are matrix-representation of bootstrap equation."
+toQboot::usage = "toQboot[sdp] converts sdp-object into c++ code for qboot."
 toCboot::usage = "toCboot[sdp] converts sdp-object into python code for cboot."
 toTeX::usage = "toTeX[eq] gives latex string of eq (you need call Print[toTeX[eq]] to paste to your tex file).
 toTeX[eqn[{a,b,...}]] gives latex string of eq with align environment (you need call Print[toTeX[eq]] to paste to your tex file)."
@@ -164,7 +166,7 @@
 	l = Length[sol = NullSpace @ eq[r, r, t]];
 	If[l == 0, Message[setOPE::imcmpt, r, r, t]; Return[]];
 	If[l != inv[{r, r}, {t}], Message[setOPE::diff, r, r, t, l, inv[{r, r}, {t}]]; Return[]];
-	even = Select[simp @ Orthogonalize[simp[sym[#, +1] & /@ sol], simp[Conjugate[#1].#2] &], AnyTrue[#, simp @ # != 0 &] &];
+	even = Select[simp @ Orthogonalize[simp[sym[#,  1] & /@ sol], simp[Conjugate[#1].#2] &], AnyTrue[#, simp @ # != 0 &] &];
 	odd  = Select[simp @ Orthogonalize[simp[sym[#, -1] & /@ sol], simp[Conjugate[#1].#2] &], AnyTrue[#, simp @ # != 0 &] &];
 	l = Length[even];
 	myAbortProtect[
@@ -337,10 +339,10 @@
 eqList[x_ == y_] := {x == y}
 eqList[True] = {}
 eqList[a_And] := List @@ a
-sum[a_?NumericQ b_, x_op] := a sum[b, x]
+sum[Times[a_?NumericQ, b_], x_op] := a sum[b, x]
 sum[a_Plus, x_op] := Plus @@ (sum[#, x] &) /@ List @@ a
 sum[0, x_op] := 0
-single[a_?NumericQ b_] := a single[b]
+single[Times[a_?NumericQ, b_]] := a single[b]
 single[a_Plus] := Plus @@ single /@ List @@ a
 single[0] = 0
 
@@ -365,7 +367,7 @@
 setOPE[o1 : op[_, _, -1, 1], o2 : op[_, _, 1 | -1, 1], o3 : op[_, t_, 1, 1 | -1]] /; ! isPseaudo[t] := setOPE[dualOp[o1], dualOp[o2], dualOp[o3]]
 setOPE[o1 : op[_, r_, 1, 1], o2 : op[_, _, -1, 1], o3 : op[_, t_, 1, 1 | -1]] /; ! isPseaudo[t] && ! isPseaudo[r] := setOPE[dualOp[o1], dualOp[o2], dualOp[o3]]
 setOPE[o1 : op[_, r_, 1 | -1, 1], o2 : op[_, s_, 1 | -1, 1], o3 : op[_, t_, 1 | -1, 1 | -1]] :=
-	Module[{re, im, a, b, c, d, e, v, cnt, d1 = dualOp[o1], d2 = dualOp[o2], d3 = dualOp[o3], n, N = inv[r, s, t], sol, add, eq, ev0, res},
+	Module[{re, im, a, b, c, d, cnt, n, N = inv[r, s, t], sol, add, eq, ev0, res},
 		add[x_] := Sow[Join[Array[D[x, re[#]] &, N], Array[D[x, im[#]] &, N]]];
 		a = Array[L[L[o1, o2, o3, #], re[#] + I im[#]] &, N];
 		b = flip12[a];
@@ -383,7 +385,7 @@
 		Scan[ev0, Keys[cnt]];
 		myAbortProtect @ Scan[ev, Keys[cnt]]]
 setOPES[o1 : op[_, r_, 1 | -1, 1], o2 : op[_, s_, 1 | -1, 1], o3 : op[_, t_, 1 | -1, 1]] :=
-	Module[{re, im, eq, a123, b123, e, v, cnt, x, N = inv[r, s, t], a213, a132, a231, a312, a321, b213, b132, b231, b312, b321, sol, add, n, ev0, res},
+	Module[{re, im, eq, a123, b123, cnt, x, N = inv[r, s, t], a213, a132, a231, a312, a321, b213, b132, b231, b312, b321, sol, add, n, ev0, res},
 		add[x_] := Sow[Join[Array[D[x, re[#]] &, N], Array[D[x, im[#]] &, N]]];
 		a123 = Array[L[L[o1, o2, o3, #], re[#] + I im[#]] &, N];
 		a312 = rotate[a123];
@@ -462,7 +464,7 @@
 sign[-1, x_] := SuperMinus[x]
 
 extract2[a_Plus, o_] := extract2[#, o] & /@ a
-extract2[x_?NumericQ y_, o_] := x extract2[y, o]
+extract2[Times[x_?NumericQ, y_], o_] := x extract2[y, o]
 extract2[sum[x_, o_op], o_] := sum[x, o]
 extract2[sum[_, _op], _] := 0
 extract2[x_single, unit | scalar] := x
@@ -486,7 +488,7 @@
 
 tmpsec = newSet[]
 tmpf[a_ + b_] := (tmpf[a]; tmpf[b];)
-tmpf[a_?NumericQ b_] := (tmpf[b];)
+tmpf[Times[a_?NumericQ, b_]] := (tmpf[b];)
 tmpf[sum[x_, o_op]] := add[tmpsec, o]
 sector[eqn[eq_List]] := Module[{pr},
 	clear[tmpsec];
@@ -589,7 +591,7 @@
 				TemplateApply["\\frac{`x``s`}{`y`}", <|"x" -> If[n == 1, "", " " <> numToTeX[n]], "y" -> numToTeX[d], "s" -> s|>]
 			, r != 1,
 				TemplateApply["`x` `s`", <|"x" -> numToTeX[r], "s" -> s|>]
-			, _,
+			, True,
 				s]]
 
 termToTeX[1, sum[x_, op[op, rp_, 1, l_]]] :=
@@ -625,11 +627,11 @@
 connectedComponents[eqn[_, z_]] :=
 	Module[{f, uf = newUF[]},
 		f[a_Plus] := Scan[f, a];
-		f[x_?NumericQ y_] := f[y];
-		f[sum[(_F | _H) a_ b_, _]] := unite[uf, a, b];
-		f[sum[(_F | _H) a_^2, _]] := add[uf, a];
-		f[single[(_Fp | _Hp) a_ b_]] := unite[uf, a, b];
-		f[single[(_Fp | _Hp) a_^2]] := add[uf, a];
+		f[Times[x_?NumericQ, y_]] := f[y];
+		f[sum[Times[_F | _H, a_, b_], _]] := unite[uf, a, b];
+		f[sum[Times[_F | _H, a_^2], _]] := add[uf, a];
+		f[single[Times[_Fp | _Hp, a_, b_]]] := unite[uf, a, b];
+		f[single[Times[_Fp | _Hp, a_^2]]] := add[uf, a];
 		Scan[f, z];
 		classify[uf]]
 
@@ -639,36 +641,104 @@
 		ind[_] := 0;
 		mat = Array[0&, {size, size}];
 		f[a_Plus] := Scan[f, a];
-		f[x_?NumericQ y_] := f[x, y];
+		f[Times[x_?NumericQ, y_]] := f[x, y];
 		f[x : _sum | _single] := f[1, x];
 		add[_, _, ___, 0, ___] := Null;
 		add[c_, x_, n_, m_] := (mat[[n, m]] += c x / 2; mat[[m, n]] += c x / 2;);
 		add[c_, x_, n_] := (mat[[n, n]] += c x;);
-		f[c_, sum[(x:_F | _H) a_ b_, _]] := add[c, x, ind[a], ind[b]];
-		f[c_, sum[(x:_F | _H) a_^2, _]] := add[c, x, ind[a]];
-		f[c_, single[(x:_Fp | _Hp) a_ b_]] := add[c, x, ind[a], ind[b]];
-		f[c_, single[(x:_Fp | _Hp) a_^2]] := add[c, x, ind[a]];
+		f[c_, sum[Times[x:_F | _H, a_, b_], _]] := add[c, x, ind[a], ind[b]];
+		f[c_, sum[Times[x:_F | _H, a_^2], _]] := add[c, x, ind[a]];
+		f[c_, single[Times[x:_Fp | _Hp, a_, b_]]] := add[c, x, ind[a], ind[b]];
+		f[c_, single[Times[x:_Fp | _Hp, a_^2]]] := add[c, x, ind[a]];
 		f[z];
 		mat
 	]
 
 makeMat[z:eqn[_, eq_List]] := Module[{con = connectedComponents[z], c},
-	If[Length[con] == 0, {{{# /. single[x_] :> x}} & /@ eq}, Table[makeMatHelper[#, c]& /@ eq, {c, con}]]]
+	If[Length[con] == 0, {{{1}, {{# /. single[x_] :> x}} & /@ eq}}, Table[{c, makeMatHelper[#, c]& /@ eq}, {c, con}]]]
 
 makeSDP[z:eqn[eq_List]] := Module[{mat, tmp, sec = sector[z], secs = <||>, scalnum, val = newSet[], f, s, i},
 		f[x_Plus] := f /@ x;
 		f[0] = 0;
-		f[x_?NumericQ y_] := If[x > 0, add[val, x]; block[1, x, y], If[x != -1, add[val, -x]]; block[-1, -x, y]];
+		f[Times[x_?NumericQ, y_]] := If[x > 0, add[val, x]; block[1, x, y], If[x != -1, add[val, -x]]; block[-1, -x, y]];
 		f[y_] := block[1, 1, y];
 		SetAttributes[f, Listable];
-		mat[unit, 1] = f[makeMat[extract[z, unit]][[1]]];
+		tmp = makeMat[extract[z, unit]][[1]];
+		mat[unit, 1] = {tmp[[1]], f[tmp[[2]]]};
 		tmp = makeMat[extract[z, scalar]];
 		scalnum = Length[tmp];
-		Do[mat[scalar, i] = f[tmp[[i]]], {i, scalnum}];
-		Do[tmp = makeMat[extract[z, s]]; secs[s] = Length[tmp]; Do[mat[s, i] = f[tmp[[i]]], {i, secs[s]}];, {s, Drop[sec, 2]}];
+		Do[mat[scalar, i] = {tmp[[i, 1]], f[tmp[[i, 2]]]}, {i, scalnum}];
+		Do[tmp = makeMat[extract[z, s]]; secs[s] = Length[tmp]; Do[mat[s, i] = {tmp[[i, 1]], f[tmp[[i, 2]]]}, {i, secs[s]}];, {s, Drop[sec, 2]}];
 		sdpobj[secs, scalnum, keys[val], mat]
 	]
 
+betaToString[1] = 1
+betaToString[op[x_, r_, p : 1 | -1, 1 | -1]] /; !isComplex[r] := TemplateApply["`sg``x`", <|"sg"->If[p > 0, "", "~"], "x"->ToString@x|>]
+betaToString[op[x_, r_?isComplex, 1, 1 | -1]] := TemplateApply["`sg``x`", <|"sg"->If[r === minrep[r, dual[r]], "", "~"], "x"->ToString@x|>]
+betaToString[\[Beta][a_, b_, c_][n_]] := TemplateApply["beta[`a`, `b`, `c`][`n`]", <|"a"->betaToString@a, "b"->betaToString@b, "c"->betaToString@c, "n"->ToString[n]|>]
+toCString[F[x_, y_, z_, w_]] := TemplateApply["ext(\"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
+toCString[H[x_, y_, z_, w_]] := TemplateApply["ext(\"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
+toCString[Fp[x_, y_, z_, w_, 0]] := TemplateApply["one, ext(\"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
+toCString[Hp[x_, y_, z_, w_, 0]] := TemplateApply["one, ext(\"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
+toCString[Fp[x_, y_, z_, w_, o_]] := TemplateApply["ops.at(\"`o`\"), ext(\"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w, "o"->o|>];
+toCString[Hp[x_, y_, z_, w_, o_]] := TemplateApply["ops.at(\"`o`\"), ext(\"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w, "o"->o|>];
+
+toQboot[sdpobj[secs_, scalarnum_, vals_, mat_]] :=
+	Module[{s, v, revval = reverseIndex[vals], numeq = Length @ mat[unit, 1][[2]], secsStr, valsStr, filename, extdelta, exts, extops, scalarsecs, secname, matsize, terms, syms, n, i, e, r, c, f, regsym, eqStr, convert},
+		terms = ConstantArray[{}, numeq];
+		syms = ConstantArray[None, numeq];
+		Do[matsize[scalar, i] = Length @ mat[scalar, i][[2, 1]], {i, scalarnum}];
+		Do[matsize[s, i] = Length @ mat[s, i][[2, 1]], {s, Keys[secs]}, {i, secs[s]}];
+		secname[unit] = secname[unit, 1] = "\"unit\"";
+		secname[scalar, id_] := If[scalarnum > 1, TemplateApply["\"(scalar, `n`)\"", <|"n" -> id - 1|>], "\"scalar\""];
+		secname[op[op, r_, 1, p_], id_] := TemplateApply["\"(`r`, `p``n`)\"", <|"r" -> r, "p" -> If[p > 0, "even", "odd"], "n" -> If[secs[op[op, r, 1, p]] > 1, ", " <> ToString[id - 1], ""]|>];
+		secsStr = StringRiffle[Flatten[Function[s, Array[Function[i, TemplateApply[secTemplate, <|"sec" -> secname[s, i], "p" -> (1 - s[[4]])/2, "sz" -> matsize[s, i], "opes" -> StringRiffle[betaToString /@ mat[s, i][[1]], ", "]|>]], secs[s]]] /@ Keys[secs]], "\n\t"];
+		valsStr = StringRiffle[Array[TemplateApply["val[`i`] = `v`;", <|"i" -> # - 1, "v" -> ToCpp`cppeval @ vals[[#]]|>] &, Length[vals]], "\n\t"];
+		exts = DeleteDuplicatesBy[Select[keys[allopsForBoot], #[[2]] === minrep[#[[2]], dual[#[[2]]]] && #[[3]] > 0 &], {#[[1]], #[[2]]} &];
+		filename = StringRiffle[TemplateApply["deltas.at(\"`k`\").str('#')", <|"k" -> #[[1]]|>] & /@ exts, " + \"-\" + "];
+		extdelta = StringRiffle[Array[TemplateApply["deltas[\"`k`\"] = parse(args[`i`]).value();", <|"i" -> #, "k" -> exts[[#, 1]]|>] &, Length[exts]], "\n\t"];
+		extops = StringRiffle[TemplateApply["ops.emplace(\"`k`\", Op(real(deltas.at(\"`k`\")), 0, c));  // `r`", <|"k" -> #[[1]], "r" -> #[[2]]|>] & /@ exts, "\n\t"];
+		scalarsecs = StringRiffle[Array[TemplateApply["// {`opes`}\n\tsecs.emplace_back(`sec`, `s`);", <|"sec" -> secname[scalar, #], "s" -> matsize[scalar, #], "opes" -> StringRiffle[betaToString /@ mat[scalar, #][[1]], ", "]|>] &, scalarnum], "\n\t"];
+		regsym[e_, F] := (syms[[e]] = "Odd";);
+		regsym[e_, H] := (syms[[e]] = "Even";);
+		regsym[e_, Fp] := (syms[[e]] = "Odd";);
+		regsym[e_, Hp] := (syms[[e]] = "Even";);
+		(* diagonal *)
+		convert[0, block[1, 1, bl_]] := toCString[bl];
+		convert[0, block[-1, 1, bl_]] := TemplateApply["real(-1), `bl`", <|"bl" -> toCString[bl]|>];
+		convert[0, block[1, v_, bl_]] := TemplateApply["val[`i`], `bl`", <|"i" -> revval[v], "bl" -> toCString[bl]|>];
+		convert[0, block[-1, v_, bl_]] := TemplateApply["-val[`i`], `bl`", <|"i" -> revval[v], "bl" -> toCString[bl]|>];
+		(* off-diagonal *)
+		convert[1, block[1, 1, bl_]] := TemplateApply["real(2), `bl`", <|"bl" -> toCString[bl]|>];
+		convert[1, block[-1, 1, bl_]] := TemplateApply["real(-2), `bl`", <|"bl" -> toCString[bl]|>];
+		convert[1, block[1, v_, bl_]] := TemplateApply["2 * val[`i`], `bl`", <|"i" -> revval[v], "bl" -> toCString[bl]|>];
+		convert[1, block[-1, v_, bl_]] := TemplateApply["-2 * val[`i`], `bl`", <|"i" -> revval[v], "bl" -> toCString[bl]|>];
+		convert[s_, r_, r_, x_] := TemplateApply["eq.add(`sec`, `r`, `c`, `block`);", <|"sec" -> s, "r" -> r - 1, "c" -> r - 1, "block" -> convert[0, x]|>];
+		convert[s_, r_, c_, x_] := TemplateApply["eq.add(`sec`, `r`, `c`, `block`);", <|"sec" -> s, "r" -> r - 1, "c" -> c - 1, "block" -> convert[1, x]|>];
+		f[s_, n_, e_, r_, c_, 0] := None;
+		f[s_, n_, e_, r_, c_, x_Plus] := Scan[f[s, n, e, r, c, #] &, List @@ x];
+		f[s_, n_, e_, r_, c_, x_] := (regsym[e, x[[3, 0]]]; If[r <= c, AppendTo[terms[[e]], convert[secname[s, n], r, c, x]]]);
+		Do[f[unit, 1, e, 1, 1, mat[unit, 1][[2, e, 1, 1]]], {e, numeq}];
+		Do[f[scalar, n, e, r, c, mat[scalar, n][[2, e, r, c]]], {n, scalarnum}, {e, numeq}, {r, matsize[scalar, n]}, {c, matsize[scalar, n]}];
+		Do[f[s, n, e, r, c, mat[s, n][[2, e, r, c]]], {s, Keys[secs]}, {n, secs[s]}, {e, numeq}, {r, matsize[s, n]}, {c, matsize[s, n]}];
+		eqStr = StringRiffle[Array[TemplateApply[eqTemplate, <|"sym" -> syms[[#]], "terms" -> StringRiffle[terms[[#]], "\n\t\t"]|>] &, numeq], "\n\t"];
+		ToCpp`createCpp[secsStr, scalarsecs, valsStr, Length[vals], eqStr, Length[exts] + 1, extops, extdelta, filename]
+	]
+
+eqTemplate = "{
+		Eq eq(boot, `sym`);
+		`terms`
+		boot.add_equation(eq);
+	}"
+
+secTemplate = "{
+		// {`opes`}
+		Sector s(`sec`, `sz`, ContinuousType);
+		for (const auto& spin : spins)
+			if (spin % 2 == `p`) s.add_op(spin);
+		secs.push_back(s);
+	}"
+
 toString[F[x_, y_, z_, w_]] := TemplateApply["get(F, \"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
 toString[H[x_, y_, z_, w_]] := TemplateApply["get(H, \"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
 toString[Fp[x_, y_, z_, w_, 0]] := TemplateApply["get(F, \"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
@@ -680,15 +750,16 @@
 secs = <|op[op, rep[1], 1, 1] -> 3, ...|>
 scalarnum = 2
 vals = {1/2, Sqrt[2], Pi, ...}
-mat[sec, num] = {{{block[-1,Sqrt[2],F[e,v,e,v]]+block[1,Sqrt[3]/2,F[e,e,v,v]], ...}, ...}, ...}
+mat[sec, num] = {{\[Beta][op[], op[], op[]][1], \[Beta][op[], op[], op[]][1]}, {{{block[-1,Sqrt[2],F[e,v,e,v]]+block[1,Sqrt[3]/2,F[e,e,v,v]], ...}, ...}, ...}}
 *)
 toCboot[sdpobj[secs_, scalarnum_, vals_, mat_]] :=
-	Module[{s, v, o, revval = reverseIndex[vals], blk, f, convert, make, tmp, secsStr, valsStr, rmats, smats, umats, filename},
+	Module[{s, v, revval = reverseIndex[vals], blk, f, convert, make, tmp, secsStr, valsStr, rmats, smats, umats, filename, opinfo},
 		secsStr = StringRiffle[TemplateApply["(\"`r`\", `p`): `n`", <|"r" -> #[[2]], "p" -> (1 - #[[4]])/2, "n" -> secs[#]|>] & /@ Keys[secs], ", "];
 		valsStr = tensorToString[vals, ToPython`pyeval];
 		rmats = newSet[];
 		smats = newSet[];
 		filename = StringRiffle[TemplateApply["{0[`k`]}", <|"k" -> #|>] & /@ DeleteDuplicates[#[[1]] & /@ keys[allopsForBoot]], "-"];
+		opinfo = StringRiffle[TemplateApply["# `o`: `r`", <|"o" -> #[[1]], "r" -> #[[2]]|>] & /@ DeleteDuplicates[{#[[1]], #[[2]]} & /@ keys[allopsForBoot]], "\n"];
 		(* convert elements of mats into string. *)
 		convert[0] = "0";
 		convert[block[1, 1, bl_]] := TemplateApply["bl[`i`]", <|"i" -> blk[bl]|>];
@@ -702,16 +773,17 @@
 		f[block[_, _, bl_]] := If[!KeyExistsQ[blk, bl], blk[bl] = Length[blk];];
 		make[m_, init_] := Module[{tmp = <||>},
 			Scan[(tmp[#] = init[#]) &, Keys[init]];
-			blk = <||>; f[m];
+			blk = <||>; f[m[[2]]];
 			tmp["bl"] = tensorToString[Keys[blk], toString];
-			tmp["mats"] = If[Length[m[[1]]] != 1, tensorToString[m, convert], tensorToString[#[[1, 1]] & /@ m , convert]];
+			tmp["opes"] = StringRiffle[betaToString /@ m[[1]], ", "];
+			tmp["mats"] = If[Length[m[[2, 1]]] != 1, tensorToString[m[[2]], convert], tensorToString[#[[1, 1]] & /@ m[[2]], convert]];
 			tmp];
 		Do[add[rmats, make[mat[s, n], <|"r" -> s[[2]], "p" -> (1 - s[[4]])/2, "n" -> n - 1|>]], {s, Keys[secs]}, {n, secs[s]}];
 		Do[add[smats, make[mat[scalar, n], <|"n" -> n - 1|>]], {n, scalarnum}];
 		rmats = StringRiffle[TemplateApply[rmatstemplate, #] & /@ keys[rmats], "\n"];
 		smats = StringRiffle[TemplateApply[smatstemplate, #] & /@ keys[smats], "\n"];
 		umats = TemplateApply[umatstemplate, make[mat[unit, 1], <||>]];
-		ToPython`createPython[secsStr, scalarnum, valsStr, rmats, smats, umats, filename]
+		ToPython`createPython[secsStr, scalarnum, valsStr, rmats, smats, umats, filename, opinfo]
 	]
 
 tensorToString[tensor_List, conv_] := "[" <> StringRiffle[tensorToString[#, conv] & /@ tensor, ", "] <> "]"
@@ -719,10 +791,12 @@
 
 rmatstemplate = "	if sector == \"`r`\" and spin % 2 == `p` and num == `n`:
 		bl = `bl`
+		# {`opes`}
 		return `mats`"
 
 smatstemplate = "	if num == `n`:
 		bl = `bl`
+		# {`opes`}
 		return `mats`"
 
 umatstemplate = "bl = `bl`
@@ -736,11 +810,11 @@
 		add[a_, a_] := add[a];
 		add[a_, b_] := (vs[a] = 1; vs[b] = 1; es[Sort[{a, b}]] = 1);
 		f[a_Plus] := Scan[f, List @@ a];
-		f[_?NumericQ y_] := f[y];
-		f[sum[(_F | _H) a_ b_, _]] := add[a, b];
-		f[single[(_Fp | _Hp) a_ b_]] := add[a, b];
-		f[sum[(_F | _H) a_^2, _]] := add[a, a];
-		f[single[(_Fp | _Hp) a_^2]] := add[a, a];
+		f[Times[_?NumericQ, y_]] := f[y];
+		f[sum[Times[_F | _H, a_, b_], _]] := add[a, b];
+		f[single[Times[_Fp | _Hp, a_, b_]]] := add[a, b];
+		f[sum[Times[_F | _H, a_^2], _]] := add[a, a];
+		f[single[Times[_Fp | _Hp, a_^2]]] := add[a, a];
 		Scan[f, eq];
 		Graph[Labeled[#, format[#]] & /@ Keys[vs], #[[1]] <-> #[[2]] & /@ Keys[es], PlotLabel -> ToString[label, StandardForm]]]
 
diff --git a/ngroup.m b/ngroup.m
index 12bbe63..c5904f8 100644
--- a/ngroup.m
+++ b/ngroup.m
@@ -31,7 +31,7 @@
 minrep::usage = "minrep[r,s] gives r if r < s else s. r and s are irrep-objects."
 setGroup::usage = "setGroup[G] loads inv.m with global symmetry G. This action clears all values calculated by inv.m previously."
 available::usage = "available[g,i] gives whether group[g,i] are supported or not."
-setPrecision::usage = "After calling setPrecision[prec], all calculation in this package will be done in precision prec and any number less than "<>ToString[Superscript[10,"prec-10"]]<>" will be choped. "<>
+setPrecision::usage = "After calling setPrecision[prec], all calculation in this package will be done in precision prec and any number less than 1/"<>ToString[Superscript[10,"prec-10"]]<>" will be choped. "<>
 "It is assumed that prec is sufficiently bigger than 10 and setPrecision is called just once just after loading this package."
 
 (* all irrep-objects of G=group[g,i] are rep[1], rep[2], ..., rep[n] (n=G[ncg]). *)
diff --git a/ngroupd.m b/ngroupd.m
index 0e09490..6eaab0c 100644
--- a/ngroupd.m
+++ b/ngroupd.m
@@ -4,8 +4,8 @@
 
 getDihedral::usage = "getDihedral[n] returns group-object dih[n] which represents the dihedral group of order 2n."
 getDicyclic::usage = "getDicyclic[n] returns group-object dic[n] which represents the dicyclic group of order 4n."
-dih::usage = "dih[n] is a group-object which is the dihedral group of order 2n. Before using this value, you have to call dih[n] to get proper group-object."
-dic::usage = "dic[n] is a group-object which is the dicyclic group of order 4n. Before using this value, you have to call dic[n] to get proper group-object."
+dih::usage = "dih[n] is a group-object which is the dihedral group of order 2n. Before using this value, you have to call getDihedral[n] to get proper group-object."
+dic::usage = "dic[n] is a group-object which is the dicyclic group of order 4n. Before using this value, you have to call getDicyclic[n] to get proper group-object."
 
 (* if n is even, all irrep-objects of dih[n] are i[1,1], i[1,-1], i[-1,1], i[-1,-1], v[1], ..., v[n/2-1]. *)
 (* if n is odd, all irrep-objects of dih[n] are i[1], i[-1] v[1], ..., v[(n-1)/2]. *)
diff --git a/ninv.m b/ninv.m
index 94de5c9..f5fb7bf 100644
--- a/ninv.m
+++ b/ninv.m
@@ -1,6 +1,7 @@
 (* Do not import this package directly. *)
 Needs["CommonFunctions`", "common.m"]
 Needs["ToPython`", "topy.m"]
+Needs["ToCpp`", "tocpp.m"]
 NClebschGordan`clearCG[];
 
 BeginPackage["NClebschGordan`"]
@@ -64,7 +65,7 @@
 
 (* fundamental scalars are used to seperate sum of conformal blocks over scalars to single[...] and sum[...]. *)
 setOps::usage = "setOps[ops] registers ops and duals of ops as fundamental scalars."
-one::usage = "one represents unit operator. This is implicitly registered as a fundamental scalar."
+one::usage = "one represents the unit operator. This is implicitly registered as a fundamental scalar."
 
 (* extract is a projection. x==extract[x,unit]+extract[x,scalar]+\[Sum]_{r,p}extract[x,op[op,r,1,p]] and extract[x,sec]==extract[extract[x,sec]]. *)
 extract::usage = "extract[x,op[op,r,1,p]] extracts terms of the form sum[...,op[op,r,1,p]] from x.
@@ -72,23 +73,25 @@
 extract[x,unit] extracts terms of the form single[Fp or Hp] from x (contribution of unit operator)."
 unit::usage = "unit is a option for extract and means contribution of unit operator."
 scalar::usage = "scalar is a option for extract and means contribution of fundamental scalars but unit operator."
-sector::usage = "sector[eq] gives a list of all nontrivial option sec for extract applied to eq, i.e. {sec | extract[eq,sec]!=0}. "
+sector::usage = "sector[eq] gives a list of all nontrivial option sec for extract applied to eq, i.e. {sec | extract[eq,sec]!=0}."
 
 makeG::usage = "makeG[eqn[sec,{a,b,...}]] gives an undirected graph whose vertices are OPE coefficients \[Beta] in extracted bootstrap equation eqn[sec,{a,b,...}]."
 makeMat::usage = "makeMat[eqn[sec,{a,b,...}]] gives a matrix-representation of extracted bootstrap equation eqn[sec,{a,b,...}]."
 makeSDP::usage = "makeSDP[eqn[{a,b,...}]] converts whole bootstrap equation eqn[{a,b,...}] into sdp-object."
 sdpobj::usage = "sdpobj[secs,scalarnum,vals,mats] is a sdp-object. secs is section data of bootstrap equation. scalarnum is the number of connected components in scalar sections. vals are real constants in bootstrap equation. mats are matrix-representation of bootstrap equation."
+toQboot::usage = "toQboot[sdp] converts sdp-object into c++ code for qboot."
 toCboot::usage = "toCboot[sdp] converts sdp-object into python code for cboot."
 toTeX::usage = "toTeX[eq] gives latex string of eq (you need call Print[toTeX[eq]] to paste to your tex file).
 toTeX[eqn[{a,b,...}]] gives latex string of eq with align environment (you need call Print[toTeX[eq]] to paste to your tex file)."
 repToTeX::usage = "repToTeX[r] is needed to transrate irrep-object r as latex string. Please set appropriate value."
-opToTeX::usage = "opToTeX[o] is needed to transrate operator object o as latex string. Please set appropriate value."
+opToTeX::usage = "opToTeX[o] is needed to transrate operator name o as latex string. Please set appropriate value."
 
 Begin["`Private`"]
 
 setOPE::usage = "setOPE[r,s,t] calculates all values of ope[r,s,t].
 setOPE[r,s] calls setOPE[r,s,t] for all t in prod[r,s]."
 setAllOPE::usage = "setAllOPE[reps] calls setOPE[r,s] and setOPE[t,dual[t],id] for all r,s in reps and t in prod[r,s]."
+opToTeX2::usage = "opToTeX2[op[...]] is needed to transrate operator object op[...] as latex string."
 precision::usage = "precision represents internal precision for calculation."
 
 allPublicSymbol = {clearCG, inv, ope, cor,
@@ -135,6 +138,7 @@
 x:invs[r1_, r2_, r3_, r4_] := myAbortProtect[x = Module[{t = inv[r1, r2, r3, r4], s, n, m},
 	MyReap @ Do[Sow[{s, n, m}], {s, Keys[t]}, {n, t[s][[1]]}, {m, t[s][[2]]}]]]
 
+(* eq[r,s,t] gives equations which Clebsch-Gordan coefficients {a,b/r,s|c/t} must satisfy. *)
 eq[r_, s_, t_] /; inv[{r, s}, {t}] > 0 :=
 Module[{a, b, c, x, y, z, X, g, e, d1 = dim[r], d2 = dim[s], d3 = dim[t], d},
 	d = d1 d2 d3;
@@ -164,7 +168,7 @@
 	l = Length[sol = NullSpace @ num @ eq[r, r, t]];
 	If[l == 0, Message[setOPE::imcmpt, r, r, t]; Return[]];
 	If[l != inv[{r, r}, {t}], Message[setOPE::diff, r, r, t, l, inv[{r, r}, {t}]]; Return[]];
-	even = Select[num @ Orthogonalize[num @ sym[#, +1] & /@ sol, Conjugate[#1].#2 &], AnyTrue[#, # != 0 &] &];
+	even = Select[num @ Orthogonalize[num @ sym[#,  1] & /@ sol, Conjugate[#1].#2 &], AnyTrue[#, # != 0 &] &];
 	odd  = Select[num @ Orthogonalize[num @ sym[#, -1] & /@ sol, Conjugate[#1].#2 &], AnyTrue[#, # != 0 &] &];
 	l = Length[even];
 	myAbortProtect[
@@ -220,9 +224,14 @@
 	1 <= n <= inv[{r1, r2}, {s}] && 1 <= m <= inv[r3, r4, s] := (setCor[r1, r2, r3, r4, s]; x)
 cor[_, _, _, _][_, _, _][_, _, _, _] := 0
 
+(* decompose invariant tensor f into linear combination of cor[r,s,t][n]. *)
 dec[r_, s_, t_][f_] := Module[{vec},
 	decPrep[r, s, t]; vec = Array[{f @@ bas[r, s, t][#]} &, inv[r, s, t]]; num @ Flatten[invmat[r, s, t].vec]]
 
+(* decPrep[r,s,t] prepares bas[r,s,t][n] and invmat for dec.
+dec could be done by just taking inner-product, but most of cor[r,s,t][n][a,b,c] are zero, we need decPrep for more efficiency.
+bas[r,s,t] are sufficient components to distinguish invariant tensors.
+mat is a matrix whose components are values of cor[r,s,t][n] evaluated at bas[r,s,t], and invmat is its inverse. *)
 x:decPrep[r_, s_, t_] /; inv[r, s, t] > 0 := x =
 Module[{n, f = cor[r, s, t], max = inv[r, s, t], mat, rev, new, old, sc, a, b, c, tmp, tmp2, res},
 	Do[If[f[1][a, b, c] != 0, res[1] = {a, b, c}; mat = {{num @ f[1][a, b, c]}}; rev = Inverse[mat]; Break[]]
@@ -312,8 +321,8 @@
 
 x:setOmega[r_, s_, t_] := x = Module[{n, m, l = inv[r, s, t], v, f, res},
 	Do[
-		f[a_, b_, c_] :=
-			num @ Sum[ope[dual[r]][a2, a] ope[dual[s]][b2, b] ope[dual[t]][c2, c] Conjugate[cor[dual[r], dual[s], dual[t]][m][a2, b2, c2]]
+		f[a_, b_, c_] := num @
+			Sum[ope[dual[r]][a2, a] ope[dual[s]][b2, b] ope[dual[t]][c2, c] Conjugate[cor[dual[r], dual[s], dual[t]][m][a2, b2, c2]]
 			, {a2, dim[r]}, {b2, dim[s]}, {c2, dim[t]}];
 		v = dec[r, s, t][f];
 		Do[res[n, m] = v[[n]], {n, l}]
@@ -358,10 +367,10 @@
 eqList[x_ == y_] := {x == y}
 eqList[True] = {}
 eqList[a_And] := List @@ a
-sum[a_?NumericQ b_, x_op] := a sum[b, x]
+sum[Times[a_?NumericQ, b_], x_op] := a sum[b, x]
 sum[a_Plus, x_op] := Plus @@ (sum[#, x] &) /@ List @@ a
 sum[0, x_op] := 0
-single[a_?NumericQ b_] := a single[b]
+single[Times[a_?NumericQ, b_]] := a single[b]
 single[a_Plus] := Plus @@ single /@ List @@ a
 single[0] = 0
 
@@ -386,7 +395,7 @@
 setOPE[o1 : op[_, _, -1, 1], o2 : op[_, _, 1 | -1, 1], o3 : op[_, t_, 1, 1 | -1]] /; ! isPseaudo[t] := setOPE[dualOp[o1], dualOp[o2], dualOp[o3]]
 setOPE[o1 : op[_, r_, 1, 1], o2 : op[_, _, -1, 1], o3 : op[_, t_, 1, 1 | -1]] /; ! isPseaudo[t] && ! isPseaudo[r] := setOPE[dualOp[o1], dualOp[o2], dualOp[o3]]
 setOPE[o1 : op[_, r_, 1 | -1, 1], o2 : op[_, s_, 1 | -1, 1], o3 : op[_, t_, 1 | -1, 1 | -1]] :=
-	Module[{re, im, a, b, c, d, e, v, cnt, d1 = dualOp[o1], d2 = dualOp[o2], d3 = dualOp[o3], n, N = inv[r, s, t], sol, add, eq, ev0, res},
+	Module[{re, im, a, b, c, d, cnt, n, N = inv[r, s, t], sol, add, eq, ev0, res},
 		add[x_] := Sow[Join[Array[D[x, re[#]] &, N], Array[D[x, im[#]] &, N]]];
 		a = Array[L[L[o1, o2, o3, #], re[#] + I im[#]] &, N];
 		b = flip12[a];
@@ -404,7 +413,7 @@
 		Scan[ev0, Keys[cnt]];
 		myAbortProtect @ Scan[ev, Keys[cnt]]]
 setOPES[o1 : op[_, r_, 1 | -1, 1], o2 : op[_, s_, 1 | -1, 1], o3 : op[_, t_, 1 | -1, 1]] :=
-	Module[{re, im, eq, a123, b123, e, v, cnt, x, N = inv[r, s, t], a213, a132, a231, a312, a321, b213, b132, b231, b312, b321, sol, add, n, ev0, res},
+	Module[{re, im, eq, a123, b123, cnt, x, N = inv[r, s, t], a213, a132, a231, a312, a321, b213, b132, b231, b312, b321, sol, add, n, ev0, res},
 		add[x_] := Sow[Join[Array[D[x, re[#]] &, N], Array[D[x, im[#]] &, N]]];
 		a123 = Array[L[L[o1, o2, o3, #], re[#] + I im[#]] &, N];
 		a312 = rotate[a123];
@@ -469,7 +478,7 @@
 reduceBoot[x_ == y_] := {Expand[x - y]}
 reduceBoot[eq_And] := Module[{add, bls = <||>, ks, scan, scan2, mat = <||>, row, n, b, min = Infinity, max = 0},
 	scan[x_, n_] := scan2[#, n]& /@ x;
-	scan2[x_?NumericQ y_, n_] := (add[y]; mat[n][bls[y]] = x);
+	scan2[Times[x_?NumericQ, y_], n_] := (add[y]; mat[n][bls[y]] = x);
 	add[x:(_sum|_single)] := If[!KeyExistsQ[bls, x], bls[x] = 1 + Length[bls]];
 	Do[mat[n] = <||>; scan[Expand[eq[[n, 1]] - eq[[n, 2]]], n], {n, Length[eq]}];
 	mat = Table[If[KeyExistsQ[mat[n], b], mat[n][b], 0], {n, Length[eq]}, {b, Length[bls]}];
@@ -497,7 +506,7 @@
 sign[-1, x_] := SuperMinus[x]
 
 extract2[a_Plus, o_] := extract2[#, o] & /@ a
-extract2[x_?NumericQ y_, o_] := x extract2[y, o]
+extract2[Times[x_?NumericQ, y_], o_] := x extract2[y, o]
 extract2[sum[x_, o_op], o_] := sum[x, o]
 extract2[sum[_, _op], _] := 0
 extract2[x_single, unit | scalar] := x
@@ -521,7 +530,7 @@
 
 tmpsec = newSet[]
 tmpf[a_ + b_] := (tmpf[a]; tmpf[b];)
-tmpf[a_?NumericQ b_] := (tmpf[b];)
+tmpf[Times[a_?NumericQ, b_]] := (tmpf[b];)
 tmpf[sum[x_, o_op]] := add[tmpsec, o]
 sector[eqn[eq_List]] := Module[{pr},
 	clear[tmpsec];
@@ -552,7 +561,7 @@
 Hp[a_, b_, c_, d_, o_] /; lexComp[{c, d, a, b}, {a, b, c, d}] < 0 := Hp[c, d, a, b, o]
 addOp[o : op[x_, r_, 1 | -1, 1]] := myAbortProtect[allops[o] = 1; add[allopsForBoot, o];
 	If[!KeyExistsQ[allreps, r], allreps[r] = <|o->1|>, allreps[r][o] = 1];
-	op[x] = op[x, r];
+	If[Length[op[x]] == 1, op[x] = op[x, r]];
 	If[!KeyExistsQ[ord, x], ord[x] = Length[ord]];]
 setOps[ops_List] := myAbortProtect[allops = <|one->1|>; allreps = <|id-><|one->1|>|>; ord = <||>; allopsForBoot = newSet[]; Scan[(addOp[#]; addOp[dualOp[#]]) &, ops]]
 
@@ -566,6 +575,9 @@
 
 opToTeX[op] := "\\op{O}"
 
+opToTeX2[op[x_, r_, p : 1 | -1, 1 | -1]] /; !isComplex[r] := TemplateApply[If[p > 0, "``", "\\overline{``}"], opToTeX[x]]
+opToTeX2[op[x_, r_?isComplex, 1, 1 | -1]] := TemplateApply[If[r === minrep[r, dual[r]], "``", "\\overline{``}"], opToTeX[x]]
+
 FtoTeX[F[a_, b_, c_, d_]] :=
 	TemplateApply[
 		"F^{`o1` `o2`,`o3` `o4`}_{-,\\op{O}}",
@@ -595,11 +607,11 @@
 \[Beta]toTeX[\[Beta][a_, b_, c_][n_]] :=
 	TemplateApply[
 		"\\lambda_{`o1` `o2` `o3`}^{(`n`)}",
-		<|"o1" -> opToTeX[a[[1]]], "o2" -> opToTeX[b[[1]]], "o3" -> opToTeX[c[[1]]], "n" -> numToTeX[n]|>]
+		<|"o1" -> opToTeX2[a], "o2" -> opToTeX2[b], "o3" -> opToTeX2[c], "n" -> numToTeX[n]|>]
 \[Beta]toTeX[\[Beta][a_, b_, c_][1]] :=
 	TemplateApply[
 		"\\lambda_{`o1` `o2` `o3`}",
-		<|"o1" -> opToTeX[a[[1]]], "o2" -> opToTeX[b[[1]]], "o3" -> opToTeX[c[[1]]]|>]
+		<|"o1" -> opToTeX2[a], "o2" -> opToTeX2[b], "o3" -> opToTeX2[c]|>]
 
 \[Beta]FtoTeX[(x : \[Beta][__][_]) (y : \[Beta][__][_]) (z_F | z_H | z_Fp | z_Hp)] :=
 	TemplateApply[
@@ -619,10 +631,10 @@
 		Which[
 			r[[0]] === Rational,
 				n = Numerator[r]; d = Denominator[r];
-				TemplateApply["\\frac{`x``s`}{`y`}", <|"x" -> If[n == 1, "", " "<>numToTeX[n]], "y" -> numToTeX[d], "s" -> s|>]
+				TemplateApply["\\frac{`x``s`}{`y`}", <|"x" -> If[n == 1, "", " " <> numToTeX[n]], "y" -> numToTeX[d], "s" -> s|>]
 			, r != 1,
 				TemplateApply["`x` `s`", <|"x" -> numToTeX[r], "s" -> s|>]
-			, _,
+			, True,
 				s]]
 
 termToTeX[one_, sum[x_, op[op, rp_, 1, l_]]] /; num[Abs[one - 1]] < epsilon :=
@@ -658,11 +670,11 @@
 connectedComponents[eqn[_, z_]] :=
 	Module[{f, uf = newUF[]},
 		f[a_Plus] := Scan[f, a];
-		f[x_?NumericQ y_] := f[y];
-		f[sum[(_F | _H) a_ b_, _]] := unite[uf, a, b];
-		f[sum[(_F | _H) a_^2, _]] := add[uf, a];
-		f[single[(_Fp | _Hp) a_ b_]] := unite[uf, a, b];
-		f[single[(_Fp | _Hp) a_^2]] := add[uf, a];
+		f[Times[x_?NumericQ, y_]] := f[y];
+		f[sum[Times[_F | _H, a_, b_], _]] := unite[uf, a, b];
+		f[sum[Times[_F | _H, a_^2], _]] := add[uf, a];
+		f[single[Times[_Fp | _Hp, a_, b_]]] := unite[uf, a, b];
+		f[single[Times[_Fp | _Hp, a_^2]]] := add[uf, a];
 		Scan[f, z];
 		classify[uf]]
 
@@ -672,36 +684,104 @@
 		ind[_] := 0;
 		mat = Array[0&, {size, size}];
 		f[a_Plus] := Scan[f, a];
-		f[x_?NumericQ y_] := f[x, y];
+		f[Times[x_?NumericQ, y_]] := f[x, y];
 		f[x : _sum | _single] := f[1, x];
 		add[_, _, ___, 0, ___] := Null;
 		add[c_, x_, n_, m_] := (mat[[n, m]] += c x / 2; mat[[m, n]] += c x / 2;);
 		add[c_, x_, n_] := (mat[[n, n]] += c x;);
-		f[c_, sum[(x:_F | _H) a_ b_, _]] := add[c, x, ind[a], ind[b]];
-		f[c_, sum[(x:_F | _H) a_^2, _]] := add[c, x, ind[a]];
-		f[c_, single[(x:_Fp | _Hp) a_ b_]] := add[c, x, ind[a], ind[b]];
-		f[c_, single[(x:_Fp | _Hp) a_^2]] := add[c, x, ind[a]];
+		f[c_, sum[Times[x:_F | _H, a_, b_], _]] := add[c, x, ind[a], ind[b]];
+		f[c_, sum[Times[x:_F | _H, a_^2], _]] := add[c, x, ind[a]];
+		f[c_, single[Times[x:_Fp | _Hp, a_, b_]]] := add[c, x, ind[a], ind[b]];
+		f[c_, single[Times[x:_Fp | _Hp, a_^2]]] := add[c, x, ind[a]];
 		f[z];
-		mat
+		N[mat, precision - 30]
 	]
 
 makeMat[z:eqn[_, eq_List]] := Module[{con = connectedComponents[z], c},
-	If[Length[con] == 0, {{{# /. single[x_] :> x}} & /@ eq}, Table[makeMatHelper[#, c]& /@ eq, {c, con}]]]
+	If[Length[con] == 0, {{{1}, {{# /. single[x_] :> x}} & /@ eq}}, Table[{c, makeMatHelper[#, c]& /@ eq}, {c, con}]]]
 
-makeSDP[z:eqn[eq_List]] := Module[{mat, tmp, sec = sector[z], secs = <||>, scalnum, val = newSet[], f, s, i},
+makeSDP[z:eqn[eq_List]] := Module[{mat, tmp, sec = sector[z], secs = <||>, scalnum, val = newSet[], f, s, i, add2},
+		add2[x_] := If[num @ Abs[x - 1] >= 10^21 epsilon, add[val, x]];
 		f[x_Plus] := f /@ x;
 		f[0] = 0;
-		f[x_?NumericQ y_] := If[x > 0, add[val, x]; block[1, x, y], add[val, -x]; block[-1, -x, y]];
+		f[Times[x_?NumericQ, y_]] := If[x > 0, add2[x]; block[1, x, y], add2[-x]; block[-1, -x, y]];
 		f[y_] := block[1, 1, y];
 		SetAttributes[f, Listable];
-		mat[unit, 1] = f[makeMat[extract[z, unit]][[1]]];
+		tmp = makeMat[extract[z, unit]][[1]];
+		mat[unit, 1] = {tmp[[1]], f[tmp[[2]]]};
 		tmp = makeMat[extract[z, scalar]];
 		scalnum = Length[tmp];
-		Do[mat[scalar, i] = f[tmp[[i]]], {i, scalnum}];
-		Do[tmp = makeMat[extract[z, s]]; secs[s] = Length[tmp]; Do[mat[s, i] = f[tmp[[i]]], {i, secs[s]}];, {s, Drop[sec, 2]}];
+		Do[mat[scalar, i] = {tmp[[i, 1]], f[tmp[[i, 2]]]}, {i, scalnum}];
+		Do[tmp = makeMat[extract[z, s]]; secs[s] = Length[tmp]; Do[mat[s, i] = {tmp[[i, 1]], f[tmp[[i, 2]]]}, {i, secs[s]}];, {s, Drop[sec, 2]}];
 		sdpobj[secs, scalnum, keys[val], mat]
 	]
 
+betaToString[op[x_, r_, p : 1 | -1, 1 | -1]] /; !isComplex[r] := TemplateApply["`sg``x`", <|"sg"->If[p > 0, "", "~"], "x"->ToString@x|>]
+betaToString[op[x_, r_?isComplex, 1, 1 | -1]] := TemplateApply["`sg``x`", <|"sg"->If[r === minrep[r, dual[r]], "", "~"], "x"->ToString@x|>]
+betaToString[\[Beta][a_, b_, c_][n_]] := TemplateApply["beta[`a`, `b`, `c`][`n`]", <|"a"->betaToString@a, "b"->betaToString@b, "c"->betaToString@c, "n"->ToString[n]|>]
+toCString[F[x_, y_, z_, w_]] := TemplateApply["ext(\"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
+toCString[H[x_, y_, z_, w_]] := TemplateApply["ext(\"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
+toCString[Fp[x_, y_, z_, w_, 0]] := TemplateApply["one, ext(\"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
+toCString[Hp[x_, y_, z_, w_, 0]] := TemplateApply["one, ext(\"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
+toCString[Fp[x_, y_, z_, w_, o_]] := TemplateApply["ops.at(\"`o`\"), ext(\"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w, "o"->o|>];
+toCString[Hp[x_, y_, z_, w_, o_]] := TemplateApply["ops.at(\"`o`\"), ext(\"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w, "o"->o|>];
+
+toQboot[sdpobj[secs_, scalarnum_, vals_, mat_]] :=
+	Module[{s, v, revval = reverseIndex[vals], numeq = Length @ mat[unit, 1][[2]], secsStr, valsStr, filename, extdelta, exts, extops, scalarsecs, secname, matsize, terms, syms, n, i, e, r, c, f, regsym, eqStr, convert},
+		terms = ConstantArray[{}, numeq];
+		syms = ConstantArray[None, numeq];
+		Do[matsize[scalar, i] = Length @ mat[scalar, i][[2, 1]], {i, scalarnum}];
+		Do[matsize[s, i] = Length @ mat[s, i][[2, 1]], {s, Keys[secs]}, {i, secs[s]}];
+		secname[unit] = secname[unit, 1] = "\"unit\"";
+		secname[scalar, id_] := If[scalarnum > 1, TemplateApply["\"(scalar, `n`)\"", <|"n" -> id - 1|>], "\"scalar\""];
+		secname[op[op, r_, 1, p_], id_] := TemplateApply["\"(`r`, `p``n`)\"", <|"r" -> r, "p" -> If[p > 0, "even", "odd"], "n" -> If[secs[op[op, r, 1, p]] > 1, ", " <> ToString[id - 1], ""]|>];
+		secsStr = StringRiffle[Flatten[Function[s, Array[Function[i, TemplateApply[secTemplate, <|"sec" -> secname[s, i], "p" -> (1 - s[[4]])/2, "sz" -> matsize[s, i], "opes" -> StringRiffle[betaToString /@ mat[s, i][[1]], ", "]|>]], secs[s]]] /@ Keys[secs]], "\n\t"];
+		valsStr = StringRiffle[Array[TemplateApply["val[`i`] = `v`;", <|"i" -> # - 1, "v" -> ToCpp`cppeval @ vals[[#]]|>] &, Length[vals]], "\n\t"];
+		exts = DeleteDuplicatesBy[Select[keys[allopsForBoot], #[[2]] === minrep[#[[2]], dual[#[[2]]]] && #[[3]] > 0 &], {#[[1]], #[[2]]} &];
+		filename = StringRiffle[TemplateApply["deltas.at(\"`k`\").str('#')", <|"k" -> #[[1]]|>] & /@ exts, " + \"-\" + "];
+		extdelta = StringRiffle[Array[TemplateApply["deltas[\"`k`\"] = parse(args[`i`]).value();", <|"i" -> #, "k" -> exts[[#, 1]]|>] &, Length[exts]], "\n\t"];
+		extops = StringRiffle[TemplateApply["ops.emplace(\"`k`\", Op(real(deltas.at(\"`k`\")), 0, c));  // `r`", <|"k" -> #[[1]], "r" -> #[[2]]|>] & /@ exts, "\n\t"];
+		scalarsecs = StringRiffle[Array[TemplateApply["// {`opes`}\n\tsecs.emplace_back(`sec`, `s`);", <|"sec" -> secname[scalar, #], "s" -> matsize[scalar, #], "opes" -> StringRiffle[betaToString /@ mat[scalar, #][[1]], ", "]|>] &, scalarnum], "\n\t"];
+		regsym[e_, F] := (syms[[e]] = "Odd";);
+		regsym[e_, H] := (syms[[e]] = "Even";);
+		regsym[e_, Fp] := (syms[[e]] = "Odd";);
+		regsym[e_, Hp] := (syms[[e]] = "Even";);
+		(* diagonal *)
+		convert[0, block[1, one_, bl_]] /; num[Abs[one - 1]] < 10^21 epsilon := toCString[bl];
+		convert[0, block[-1, one_, bl_]] /; num[Abs[one - 1]] < 10^21 epsilon := TemplateApply["real(-1), `bl`", <|"bl" -> toCString[bl]|>];
+		convert[0, block[1, v_, bl_]] := TemplateApply["val[`i`], `bl`", <|"i" -> revval[v], "bl" -> toCString[bl]|>];
+		convert[0, block[-1, v_, bl_]] := TemplateApply["-val[`i`], `bl`", <|"i" -> revval[v], "bl" -> toCString[bl]|>];
+		(* off-diagonal *)
+		convert[1, block[1, one_, bl_]] /; num[Abs[one - 1]] < 10^21 epsilon := TemplateApply["real(2), `bl`", <|"bl" -> toCString[bl]|>];
+		convert[1, block[-1, one_, bl_]] /; num[Abs[one - 1]] < 10^21 epsilon := TemplateApply["real(-2), `bl`", <|"bl" -> toCString[bl]|>];
+		convert[1, block[1, v_, bl_]] := TemplateApply["2 * val[`i`], `bl`", <|"i" -> revval[v], "bl" -> toCString[bl]|>];
+		convert[1, block[-1, v_, bl_]] := TemplateApply["-2 * val[`i`], `bl`", <|"i" -> revval[v], "bl" -> toCString[bl]|>];
+		convert[s_, r_, r_, x_] := TemplateApply["eq.add(`sec`, `r`, `c`, `block`);", <|"sec" -> s, "r" -> r - 1, "c" -> r - 1, "block" -> convert[0, x]|>];
+		convert[s_, r_, c_, x_] := TemplateApply["eq.add(`sec`, `r`, `c`, `block`);", <|"sec" -> s, "r" -> r - 1, "c" -> c - 1, "block" -> convert[1, x]|>];
+		f[s_, n_, e_, r_, c_, 0] := None;
+		f[s_, n_, e_, r_, c_, x_Plus] := Scan[f[s, n, e, r, c, #] &, List @@ x];
+		f[s_, n_, e_, r_, c_, x_] := (regsym[e, x[[3, 0]]]; If[r <= c, AppendTo[terms[[e]], convert[secname[s, n], r, c, x]]]);
+		Do[f[unit, 1, e, 1, 1, mat[unit, 1][[2, e, 1, 1]]], {e, numeq}];
+		Do[f[scalar, n, e, r, c, mat[scalar, n][[2, e, r, c]]], {n, scalarnum}, {e, numeq}, {r, matsize[scalar, n]}, {c, r, matsize[scalar, n]}];
+		Do[f[s, n, e, r, c, mat[s, n][[2, e, r, c]]], {s, Keys[secs]}, {n, secs[s]}, {e, numeq}, {r, matsize[s, n]}, {c, r, matsize[s, n]}];
+		eqStr = StringRiffle[Array[TemplateApply[eqTemplate, <|"sym" -> syms[[#]], "terms" -> StringRiffle[terms[[#]], "\n\t\t"]|>] &, numeq], "\n\t"];
+		ToCpp`createCpp[secsStr, scalarsecs, valsStr, Length[vals], eqStr, Length[exts] + 1, extops, extdelta, filename]
+	]
+
+eqTemplate = "{
+		Eq eq(boot, `sym`);
+		`terms`
+		boot.add_equation(eq);
+	}"
+
+secTemplate = "{
+		// {`opes`}
+		Sector s(`sec`, `sz`, ContinuousType);
+		for (const auto& spin : spins)
+			if (spin % 2 == `p`) s.add_op(spin);
+		secs.push_back(s);
+	}"
+
 toString[F[x_, y_, z_, w_]] := TemplateApply["get(F, \"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
 toString[H[x_, y_, z_, w_]] := TemplateApply["get(H, \"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
 toString[Fp[x_, y_, z_, w_, 0]] := TemplateApply["get(F, \"`x`\", \"`y`\", \"`z`\", \"`w`\")", <|"x"->x, "y"->y, "z"->z, "w"->w|>];
@@ -713,36 +793,40 @@
 secs = <|op[op, rep[1], 1, 1] -> 3, ...|>
 scalarnum = 2
 vals = {1/2, Sqrt[2], Pi, ...}
-mat[sec, num] = {{{block[-1,Sqrt[2],F[e,v,e,v]]+block[1,Sqrt[3]/2,F[e,e,v,v]], ...}, ...}, ...}
+mat[sec, num] = {{\[Beta][op[], op[], op[]][1], \[Beta][op[], op[], op[]][1]}, {{{block[-1,Sqrt[2],F[e,v,e,v]]+block[1,Sqrt[3]/2,F[e,e,v,v]], ...}, ...}, ...}}
 *)
 toCboot[sdpobj[secs_, scalarnum_, vals_, mat_]] :=
-	Module[{s, v, o, revval = reverseIndex[vals], blk, f, convert, make, tmp, secsStr, valsStr, rmats, smats, umats, filename},
+	Module[{s, v, revval = reverseIndex[vals], blk, f, convert, make, tmp, secsStr, valsStr, rmats, smats, umats, filename, opinfo},
 		secsStr = StringRiffle[TemplateApply["(\"`r`\", `p`): `n`", <|"r" -> #[[2]], "p" -> (1 - #[[4]])/2, "n" -> secs[#]|>] & /@ Keys[secs], ", "];
 		valsStr = tensorToString[vals, ToPython`pyeval];
 		rmats = newSet[];
 		smats = newSet[];
 		filename = StringRiffle[TemplateApply["{0[`k`]}", <|"k" -> #|>] & /@ DeleteDuplicates[#[[1]] & /@ keys[allopsForBoot]], "-"];
+		opinfo = StringRiffle[TemplateApply["# `o`: `r`", <|"o" -> #[[1]], "r" -> #[[2]]|>] & /@ DeleteDuplicates[{#[[1]], #[[2]]} & /@ keys[allopsForBoot]], "\n"];
+		(* convert elements of mats into string. *)
 		convert[0] = "0";
-		convert[block[1, 1, bl_]] := TemplateApply["bl[`i`]", <|"i" -> blk[bl]|>];
-		convert[block[-1, 1, bl_]] := TemplateApply["-bl[`i`]", <|"i" -> blk[bl]|>];
+		convert[block[1, one_, bl_]] /; num[Abs[one - 1]] < 10^21 epsilon := TemplateApply["bl[`i`]", <|"i" -> blk[bl]|>];
+		convert[block[-1, one_, bl_]] /; num[Abs[one - 1]] < 10^21 epsilon := TemplateApply["-bl[`i`]", <|"i" -> blk[bl]|>];
 		convert[block[1, v_, bl_]] := TemplateApply["val[`j`] * bl[`i`]", <|"i" -> blk[bl], "j" -> revval[v]|>];
 		convert[block[-1, v_, bl_]] := TemplateApply["-val[`j`] * bl[`i`]", <|"i" -> blk[bl], "j" -> revval[v]|>];
 		convert[x_Plus] := Module[{lx = List @@ x}, convert[lx[[1]]] <> StringJoin[Table[If[y[[1]] > 0, " + " <> convert[y], " - " <> convert[block[1, y[[2]], y[[3]]]]], {y, Rest[lx]}]]];
+		(* assign indices to block *)
 		f[x_List] := Scan[f, x];
 		f[x_Plus] := Scan[f, List @@ x];
 		f[block[_, _, bl_]] := If[!KeyExistsQ[blk, bl], blk[bl] = Length[blk];];
 		make[m_, init_] := Module[{tmp = <||>},
 			Scan[(tmp[#] = init[#]) &, Keys[init]];
-			blk = <||>; f[m];
+			blk = <||>; f[m[[2]]];
 			tmp["bl"] = tensorToString[Keys[blk], toString];
-			tmp["mats"] = If[Length[m[[1]]] != 1, tensorToString[m, convert], tensorToString[#[[1, 1]] & /@ m , convert]];
+			tmp["opes"] = StringRiffle[betaToString /@ m[[1]], ", "];
+			tmp["mats"] = If[Length[m[[2, 1]]] != 1, tensorToString[m[[2]], convert], tensorToString[#[[1, 1]] & /@ m[[2]], convert]];
 			tmp];
 		Do[add[rmats, make[mat[s, n], <|"r" -> s[[2]], "p" -> (1 - s[[4]])/2, "n" -> n - 1|>]], {s, Keys[secs]}, {n, secs[s]}];
 		Do[add[smats, make[mat[scalar, n], <|"n" -> n - 1|>]], {n, scalarnum}];
 		rmats = StringRiffle[TemplateApply[rmatstemplate, #] & /@ keys[rmats], "\n"];
 		smats = StringRiffle[TemplateApply[smatstemplate, #] & /@ keys[smats], "\n"];
 		umats = TemplateApply[umatstemplate, make[mat[unit, 1], <||>]];
-		ToPython`createPython[secsStr, scalarnum, valsStr, rmats, smats, umats, filename]
+		ToPython`createPython[secsStr, scalarnum, valsStr, rmats, smats, umats, filename, opinfo]
 	]
 
 tensorToString[tensor_List, conv_] := "[" <> StringRiffle[tensorToString[#, conv] & /@ tensor, ", "] <> "]"
@@ -750,10 +834,12 @@
 
 rmatstemplate = "	if sector == \"`r`\" and spin % 2 == `p` and num == `n`:
 		bl = `bl`
+		# {`opes`}
 		return `mats`"
 
 smatstemplate = "	if num == `n`:
 		bl = `bl`
+		# {`opes`}
 		return `mats`"
 
 umatstemplate = "bl = `bl`
@@ -767,11 +853,11 @@
 		add[a_, a_] := add[a];
 		add[a_, b_] := (vs[a] = 1; vs[b] = 1; es[Sort[{a, b}]] = 1);
 		f[a_Plus] := Scan[f, List @@ a];
-		f[_?NumericQ y_] := f[y];
-		f[sum[(_F | _H) a_ b_, _]] := add[a, b];
-		f[single[(_Fp | _Hp) a_ b_]] := add[a, b];
-		f[sum[(_F | _H) a_^2, _]] := add[a, a];
-		f[single[(_Fp | _Hp) a_^2]] := add[a, a];
+		f[Times[_?NumericQ, y_]] := f[y];
+		f[sum[Times[_F | _H, a_, b_], _]] := add[a, b];
+		f[single[Times[_Fp | _Hp, a_, b_]]] := add[a, b];
+		f[sum[Times[_F | _H, a_^2], _]] := add[a, a];
+		f[single[Times[_Fp | _Hp, a_^2]]] := add[a, a];
 		Scan[f, eq];
 		Graph[Labeled[#, format[#]] & /@ Keys[vs], #[[1]] <-> #[[2]] & /@ Keys[es], PlotLabel -> ToString[label, StandardForm]]]
 
diff --git a/tocpp.m b/tocpp.m
new file mode 100644
index 0000000..2ab8965
--- /dev/null
+++ b/tocpp.m
@@ -0,0 +1,158 @@
+BeginPackage["ToCpp`"]
+
+cppeval::usage = "cppeval[x]"
+createCpp::usage = "createCpp[secs,scalarsecs,vals,numval,eqs,extops,extdelta,filename]"
+
+Begin["`Private`"]
+
+cppeval[n_Integer] := TemplateApply["real(\"``\")", n]
+cppeval[a_Times] := StringRiffle[cppeval /@ List @@ a, {"(", " * ", ")"}]
+cppeval[a_Plus] := StringRiffle[cppeval /@ List @@ a, {"(", " + ", ")"}]
+cppeval[Power[a_, b_]] := TemplateApply["mp::pow(``, ``)", {cppeval[a], cppeval[b]}]
+cppeval[Surd[x_, n_]] /; x >= 0 := TemplateApply["mp::pow(``, ``)", {cppeval[x], cppeval[1/n]}]
+cppeval[Surd[x_, n_]] := cppeval[-Surd[-x, n]]
+cppeval[Sqrt[a_]] := TemplateApply["mp::sqrt(``)", cppeval[a]]
+cppeval[Rational[a_, b_]] := TemplateApply["(`` / ``)", {cppeval[a], cppeval[b]}]
+cppeval[Sin[a_]] := TemplateApply["mp::sin(``)", cppeval[a]]
+cppeval[Cos[a_]] := TemplateApply["mp::cos(``)", cppeval[a]]
+cppeval[Tan[a_]] := TemplateApply["mp::tan(``)", cppeval[a]]
+cppeval[Sec[a_]] := TemplateApply["mp::sec(``)", cppeval[a]]
+cppeval[Csc[a_]] := TemplateApply["mp::csc(``)", cppeval[a]]
+cppeval[Cot[a_]] := TemplateApply["mp::cot(``)", cppeval[a]]
+cppeval[Exp[a_]] := TemplateApply["mp::exp(``)", cppeval[a]]
+cppeval[Log[a_]] := TemplateApply["mp::log(``)", cppeval[a]]
+cppeval[Log[b_, a_]] := cppeval[Log[a] / Log[b]]
+cppeval[Pi] = "mp::const_pi()";
+cppeval[E] = cppeval[Exp[1]];
+cppeval[x_?ExactNumberQ] := TemplateApply["real(\"``\")", ToString @ FortranForm @ N[x, 400]]
+cppeval[x_?NumberQ] := Module[{r = Rationalize[x, 0]}, If[Abs[Denominator[r]] + Abs[Numerator[r]] < 10000, cppeval[r], TemplateApply["real(\"``\")", ToString @ FortranForm[x]]]]
+
+createCpp[secs_, scalarsecs_, vals_, numval_, eqs_, numext_, extops_, extdelta_, filename_] :=
+	TemplateApply[template, <|"secs"->secs, "scalarsecs"->scalarsecs, "vals"->vals, "numval"->numval, "eqs"->eqs, "numext"->numext, "extops"->extops, "extdelta"->extdelta, "filename"->filename|>]
+
+(* filename: "deltas[\"e\"].str(8) + \"-\" + ..." *)
+(* numext: 3 *)
+(* extdelta: "deltas[\"s\"] = real(args[1]);\n..." *)
+(* extops: "ops.emplace(\"s\", Op(deltas[\"s\"], 0, c));\n..." *)
+(* scalarsecs: "secs.emplace_back(\"(scalar, 0)\", 2);\n..." *)
+(* secs: "{\nSector s(\"even\", 2, ContinuousType);\nfor (const auto& spin: spins) if (spin % 2 == 0) s.add_op(spin);\nsecs.push_back(s);\n}\n..." *)
+(* vals: "val[0] = real(-1);\n..." *)
+(* eqs: "{\nEq eq(boot, Odd);\neq.add(\"scalar;0\", 0, 0, ops[\"s\"], ext(\"e\", \"s\", \"e\", \"s\"));\neq.add(\"odd+\", ext(\"e\", \"s\", \"e\", \"s\"));\neq.add(\"odd-\", ext(\"e\", \"s\", \"e\", \"s\"));\nboot.add_equation(eq);\n}..." *)
+
+template = "#include <array>
+#include <functional>
+#include <map>
+#include <memory>
+#include <set>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include \"qboot/qboot.hpp\"
+
+namespace fs = qboot::fs;
+namespace mp = qboot::mp;
+using mp::real, mp::rational, mp::parse;
+using qboot::Context, qboot::PolynomialProgram, qboot::BootstrapEquation, qboot::Sector;
+using qboot::algebra::Vector;
+using std::array, std::map, std::set, std::move, std::vector, std::string, std::unique_ptr;
+
+template <class T>
+using dict = map<string, T, std::less<>>;
+using Op = qboot::PrimaryOperator;
+using Eq = qboot::Equation;
+constexpr auto ContinuousType = qboot::SectorType::Continuous;
+constexpr auto Odd = qboot::algebra::FunctionSymmetry::Odd;
+constexpr auto Even = qboot::algebra::FunctionSymmetry::Even;
+
+static string name(const dict<rational>& deltas);
+static BootstrapEquation create(const Context& c, const dict<rational>& deltas, uint32_t numax,
+                                const set<uint32_t>& spins);
+
+string name(const dict<rational>& deltas)
+{
+	return string(\"sdp-\") + `filename`;
+}
+
+BootstrapEquation create(const Context& c, const dict<rational>& deltas, uint32_t numax, const set<uint32_t>& spins)
+{
+	dict<Op> ops;
+	`extops`
+	auto ext = [&ops](auto o1, auto o2, auto o3, auto o4) {
+		return array{ops.at(o1), ops.at(o2), ops.at(o3), ops.at(o4)};
+	};
+	Op one{c};
+	vector<Sector> secs;
+	// you can add discrete sectors
+	// for example, to add sector \"hoge\" whose size of the matrix is 4,
+	//   secs.emplace_back(\"hoge\", 4);
+	// if you know OPE coefficients of the sector, for example, {1.2, 0.7, -0.1, 0.3},
+	//   secs.emplace_back(\"hoge\", 4, Vector{real(\"1.2\"), real(\"0.7\"), real(\"-0.1\"), real(\"0.3\")});
+	secs.emplace_back(\"unit\", 1, Vector{real(1)});
+	`scalarsecs`
+	// you can customize spectrums
+	// example:
+	//   customize
+	//     Sector s(\"hoge\", 2, ContinuousType);
+	//     for (const auto& spin: spins)
+	//         if (spin % 2 == 0) s.add_op(spin);
+	//     secs.push_back(s);
+	//   to
+	//     Sector s(\"hoge\", 2, ContinuousType);
+	//     s.add_op(0, real(\"3\"));    // first scalar: 3 <= delta (irrelevance)
+	//     s.add_op(2);               // use unitarity bound for spin-2
+	//     s.add_op(4, real(\"8.1\"));  // first spin-4: 8.1 <= delta
+	//     // sector \"hoge\" consists of even-spin operators, and we customized spin-0, 2, 4
+	//     // use unitarity bound for other operators
+	//     for (const auto& spin: spins)
+	//         if (spin % 2 == 0 && spin >= 6) s.add_op(spin);
+	//     secs.push_back(s);
+	`secs`
+	// do not edit from here
+	Vector<real> val(`numval`);
+	`vals`
+	BootstrapEquation boot(c, secs, numax);
+	`eqs`
+	boot.finish();
+	return boot;
+}
+
+int main(int argc, char* argv[])
+{
+	// internal precision (in binary digits)
+	mp::global_prec = 1000;
+	mp::global_rnd = MPFR_RNDN;
+	// n_Max: the order of taylor expansion of gBlock (we recommend n_Max >= 0.4 * global_prec)
+	// lambda: controls the number of derivatives (z = x + sqrt(y), (der x) ^ m (der y) ^ n for m + 2 n <= lambda)
+	// dim: the dimension of the physical space
+	// numax: controls the number of poles picked in the gBlock (numax + min(spin, numax) / 2)
+	// parallel: the number of internal threads
+	constexpr uint32_t n_Max = 400, lambda = 14, numax = 6, parallel = 8;
+	const rational dim(\"3\");
+	// spins: spins for the continuous sectors
+	set<uint32_t> spins;
+	for (uint32_t s = 0; s < 27; ++s) spins.insert(s);
+	spins.merge(set{49u, 50u});
+	assert(argc == `numext`);
+	unique_ptr<char*[]> args(argv);
+	dict<rational> deltas;
+	// external scalars
+	`extdelta`
+	args.release();
+	Context c(n_Max, lambda, dim, parallel);
+	auto boot = create(c, deltas, numax, spins);
+	// to maximize (resp. minimize) OPE in \"hoge\" sector,
+	// replace the first argument from qboot::FindContradiction(...)
+	// to qboot::ExtremalOPE(true (resp. false), \"hoge\", \"unit\")
+	auto prob = boot.convert(qboot::FindContradiction(\"unit\"), parallel);
+	auto dir = fs::current_path() / name(deltas);
+	move(prob).create_input(parallel).write(dir, parallel);
+	return 0;
+}
+"
+
+Protect[cppeval, createCpp]
+
+End[ ]
+
+EndPackage[ ]
diff --git a/topy.m b/topy.m
index 097ce6f..e425e25 100644
--- a/topy.m
+++ b/topy.m
@@ -1,7 +1,7 @@
 BeginPackage["ToPython`"]
 
 pyeval::usage = "pyeval[x]"
-createPython::usage = "createPython[secs,scalarnum,vals,rmats,smats,umats,filename]"
+createPython::usage = "createPython[secs,scalarnum,vals,rmats,smats,umats,filename,opinfo]"
 
 Begin["`Private`"]
 
@@ -48,10 +48,10 @@
 pyeval[Log[b_, a_]] := TemplateApply["log(``, ``)", {pyeval[a], pyeval[b]}]
 pyeval[I] = "I";
 pyeval[x_?ExactNumberQ] := TemplateApply["context(\"``\")", ToString @ FortranForm @ N[x, 300]]
-pyeval[x_?NumberQ] := TemplateApply["context(\"``\")", ToString @ FortranForm[x]]
+pyeval[x_?NumberQ] := Module[{r = Rationalize[x, 0]}, If[Abs[Denominator[r]] + Abs[Numerator[r]] < 10000, pyeval[r], TemplateApply["context(\"``\")", ToString @ FortranForm[x]]]]
 
-createPython[secs_, scalarnum_, vals_, rmats_, smats_, umats_, filename_] :=
-	TemplateApply[template, <|"secs"->secs, "scalarnum"->scalarnum, "vals"->vals, "rmats"->rmats, "smats"->smats, "umats"->umats, "filename"->filename|>]
+createPython[secs_, scalarnum_, vals_, rmats_, smats_, umats_, filename_, opinfo_] :=
+	TemplateApply[template, <|"secs"->secs, "scalarnum"->scalarnum, "vals"->vals, "rmats"->rmats, "smats"->smats, "umats"->umats, "filename"->filename, "opinfo"->opinfo|>]
 
 template = "# -*- coding: utf-8 -*-
 from __future__ import print_function
@@ -66,6 +66,7 @@
 context = cb.context_for_scalar(epsilon=0.5, Lambda=15, Prec=800)
 spins = list(range(27)) + [49, 50]
 nu_max = 14
+`opinfo`
 # {(r, s): d, ...} means operators in irrep r and spin s must have Delta >= d.
 mygap = {}
 def gaps(deltas):