commit aa0f3ae2988d372b190b9bde2e75a6d17e744e93 Author: Ian Jauslin Date: Sun Jun 14 00:52:45 2015 +0000 Initial commit diff --git a/INSTALL b/INSTALL new file mode 100644 index 0000000..9594767 --- /dev/null +++ b/INSTALL @@ -0,0 +1,41 @@ +####################################################################### +## ## +## Compiling and installing meankondo. ## +## ## +####################################################################### + +* meankondo should work on any POSIX compliant system, such as GNU/Linux or OSX. + +* Compiling: + Run + make + + The default paths can be modified by passing the appropriate arguments to + make, as specified in the following table + compiler : CC=/usr/bin/gcc + linker : LD=/usr/bin/gcc + archive : AR=/usr/bin/ar + include : INCLUDE= + lib paths : LIB= + optimize : OPT=-O3 + + For example, + make CC=/usr/local/bin/gcc INCLUDE=/usr/local/include LIBDIR=-L/usr/local/lib + + In addition, libkondo can be linked dynamically or statically, which is + controlled by the STATIC option of the makefile. + If STATIC=0 then link dynamically + If STATIC=2 then link statically + If STATIC=1 (default) then link libkondo statically but other libraries dynamically. + + +* Installing: + Run + make install + + The default install prefix (/usr) can be changed by changing the PREFIX + variable. + + For example + make install PREFIX=/usr/local + diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..d645695 --- /dev/null +++ b/LICENSE @@ -0,0 +1,202 @@ + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..4931175 --- /dev/null +++ b/Makefile @@ -0,0 +1,161 @@ +## Copyright 2015 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# whether to link dynamically +# if static=0 then link dynamically +# if static=2 then link statically +# if static=1 then link libkondo statically but other libraries dynamically +STATIC=1 + +VERSION=1.2 + +# products of the compilation +PROJECT_BINS= meankondo numkondo meantools kondo_preprocess meantools-convert +PROJECT_SO=libkondo.so.$(VERSION) +PROJECT_LIBS=libkondo.a +PROJECT_MANS=meankondo.1 numkondo.1 meantools.1 meantools-convert.1 kondo_preprocess.1 +PROJECT_SCRIPTS=config_functions.sh +PROJECT_DOCS=meankondo-doc.html + +# debug and optimization flags +#DB= -ggdb +OPT= -O3 + +# warning flags +WARNINGS= -Wall -Wextra -Wno-strict-overflow -std=c99 -pedantic + +# installation dirs +PREFIX=/usr +BINDIR=$(PREFIX)/bin +LIBDIR=$(PREFIX)/lib +MANDIR=$(PREFIX)/share/man/man1 +SCRIPTDIR=$(PREFIX)/share/meankondo/scripts +DOCDIR=$(PREFIX)/share/doc/meankondo-doc + +# compiler +CC=/usr/bin/gcc +LD=$(CC) +AR=/usr/bin/ar + +# directories +INCLUDE = +LIB = + +# flags +override LDFLAGS +=$(LIB) +override CFLAGS +=$(INCLUDE)$(DB) $(OPT) $(WARNINGS) + +# build directories +BUILDDIR=./build +SRCDIR=./src +OBJDIR=./objs + +# objects +LIBKONDO_OBJS = $(addprefix $(OBJDIR)/,array.o cli_parser.o coefficient.o fields.o grouped_polynomial.o idtable.o istring.o number.o parse_file.o polynomial.o rational_float.o rational_int.o rcc.o tools.o) +MEANKONDO_OBJS = $(addprefix $(OBJDIR)/,meankondo.o mean.o) +NUMKONDO_OBJS = $(addprefix $(OBJDIR)/,numkondo.o flow.o) +MEANTOOLS_OBJS = $(addprefix $(OBJDIR)/,meantools.o meantools_exp.o meantools_deriv.o meantools_eval.o) +KONDO_PP_OBJS = $(addprefix $(OBJDIR)/,kondo_preprocess.o kondo.o) + + +# flags which depend on whether to link statically or dynamically +# lib flag for libkondo +LIBKONDO_FLAG=-lkondo +# additional library required for static linking +XTRA_LIBS= + +ifeq ($(STATIC),1) + # compile libkondo.a + PREREQ=static + # libkondo is linked against libm + XTRA_LIBS=-lm + # link binaries using the static library + LIBKONDO_FLAG=-l:libkondo.a + # install static lib + INSTALLLIB=install-static +else ifeq ($(STATIC),2) + # compile libkondo.a + PREREQ=static + # libkondo is linked against libm + XTRA_LIBS=-lm + # link binaries statically + override LDFLAGS += -static + INSTALLLIB=install-static +else + # compile libkondo.so.$(VERSION) + PREREQ=shared + # required flag for subsequent dynamic linking + override CFLAGS += -fPIC + INSTALLLIB=install-shared +endif + + +all: init $(PREREQ) + +# create dirs +init: + @[ -d $(OBJDIR) ] || /bin/mkdir $(OBJDIR) + @[ -d $(BUILDDIR) ] || /bin/mkdir $(BUILDDIR) + @[ -d $(BUILDDIR)/bin ] || /bin/mkdir $(BUILDDIR)/bin + @[ -d $(BUILDDIR)/lib ] || /bin/mkdir $(BUILDDIR)/lib + +# static library +static: $(PROJECT_LIBS) $(PROJECT_BINS) +# shared library +shared: $(PROJECT_SO) $(PROJECT_BINS) + +libkondo.a: $(LIBKONDO_OBJS) + $(AR) -rc $(BUILDDIR)/lib/$@ $^ + +libkondo.so.$(VERSION): $(LIBKONDO_OBJS) + $(LD) -shared -lm $(LDFLAGS) -o $(BUILDDIR)/lib/$@ $^ + ln -fs ./libkondo.so.$(VERSION) $(BUILDDIR)/lib/libkondo.so + +meankondo: $(MEANKONDO_OBJS) + $(LD) -L$(BUILDDIR)/lib $(LDFLAGS) -o $(BUILDDIR)/bin/$@ $^ $(LIBKONDO_FLAG) -lpthread $(XTRA_LIBS) + +numkondo: $(NUMKONDO_OBJS) + $(LD) -L$(BUILDDIR)/lib $(LDFLAGS) -o $(BUILDDIR)/bin/$@ $^ $(LIBKONDO_FLAG) -lm $(XTRA_LIBS) + +meantools: $(MEANTOOLS_OBJS) + $(LD) -L$(BUILDDIR)/lib $(LDFLAGS) -o $(BUILDDIR)/bin/$@ $^ $(LIBKONDO_FLAG) $(XTRA_LIBS) + +meantools-convert: + cp scripts/meantools-convert $(BUILDDIR)/bin/ + +kondo_preprocess: $(KONDO_PP_OBJS) + $(LD) -L$(BUILDDIR)/lib $(LDFLAGS) -o $(BUILDDIR)/bin/$@ $^ $(LIBKONDO_FLAG) $(XTRA_LIBS) + +%.o : ../$(SRCDIR)/%.c + $(CC) -c $(CFLAGS) $< -o $@ + +install: $(INSTALLLIB) all + mkdir -p $(BINDIR) $(MANDIR) $(SCRIPTDIR) $(DOCDIR) + install -Dm755 $(addprefix $(BUILDDIR)/bin/,$(PROJECT_BINS)) $(BINDIR)/ + install -Dm644 $(addprefix man/,$(PROJECT_MANS)) $(MANDIR)/ + gzip $(addprefix $(MANDIR)/,$(PROJECT_MANS)) + install -Dm644 $(addprefix scripts/,$(PROJECT_SCRIPTS)) $(SCRIPTDIR)/ + install -Dm644 $(addprefix doc/,$(PROJECT_DOCS)) $(DOCDIR)/ + +install-static: all + mkdir -p $(LIBDIR) + install -Dm755 $(addprefix $(BUILDDIR)/lib/,$(PROJECT_LIBS)) $(LIBDIR)/ +install-shared: all + mkdir -p $(LIBDIR) + install -Dm755 $(addprefix $(BUILDDIR)/lib/,$(PROJECT_SO)) $(LIBDIR)/ + for lib in $(PROJECT_SO); do ln -fs ./$$lib $(LIBDIR)/$${lib%.so.*}.so; done + +clean: + @rm -rf $(OBJDIR) + @rm -rf $(BUILDDIR) diff --git a/NOTICE b/NOTICE new file mode 100644 index 0000000..1469df9 --- /dev/null +++ b/NOTICE @@ -0,0 +1,2 @@ +meankondo +Copyright 2015 Ian Jauslin diff --git a/doc/meankondo-doc.html b/doc/meankondo-doc.html new file mode 100644 index 0000000..d1f0694 --- /dev/null +++ b/doc/meankondo-doc.html @@ -0,0 +1,298 @@ + + + + + + + + + + +

meankondo v1.2

+ +

+ This is the official documentation for meankondo, version 1.2. The aim of this document is not to give a technical description of how to use the various programs bundled with meankondo, nor is it to explain where hierarchical models come from and what their meaning is, but rather a conceptual overview of how meankondo approaches the computation of flow equations, and how its programs can be made to interact with one another to compute various quantities. For a more technical description, see the man pages included with the meankondo source code. For a more theoretical discussion of Fermionic hierarchical models, see [G.Benfatto, G.Gallavotti, I.Jauslin, 2015]. +

+ +

Table of contents

+
+ +
+ +

Description

+

+ meankondo is a collection of tools to compute and manipulate the renormalization group flow in Fermionic hierarchical models. The programs included in meankondo are the following: +

+ as well as pre-processors, whose purpose is to help with writing configuration files for specific models: + + In addition, meankondo includes a library, libkondo, which can either be compiled as a shared or a static object, and contains the various structures and functions meankondo is built with. +

+ +

+ Remark: The name "meankondo" comes from the fact that it was originally developed for the hierarchical Kondo model, though the tools in meankondo are quite versatile and can be used for a wide variety of Fermionic hierarchical models. +

+ +

Quickstart

+

+ We first discuss the more elementary commands that can be used to compute and iterate flow equations. The rest of this document is dedicated to what flow equations are and how meankondo represents and manipulates them. +

+ +

+ Given a configuration file 'config', the flow equation can be computed by + + meankondo config + + and it can be iterated for, say, 100 steps starting from \(\ell_0^{[m]}=-0.01\) using + + meankondo -C config | numkondo -N 100 -I "0:-0.01" + +

+ +

Fermionic hierarchical models

+

+ In this section, we discuss how the models that meankondo can process are defined. A model is specified by a collection of fields, and a propagator between pairs of fields. +

+ +

Fields

+

+ Fields are the elementary objects in terms of which a model is defined. Fields can be one of +

+ The fields are used as a basis for a complex algebra, so that we can take products and linear combinations of fields (in other words, the concept of polynomials over the fields is well defined). Some of the fields (Fermions) anti-commute with each other (two fields \(a\) and \(b\) are said to anti-commute if \(ab\equiv-ba\)), and the rest (Bosons) commute. Which fields are Fermions and which are Bosons is specified in the #!fields entry in the configuration file. (Warning: As of version 1.2, all internal fields must be Fermions.) +

+

+ In the configuration file of the meankondo program, the fields are specified in the #!fields entry. +

+ +

Propagator

+

+ Given a collection of fields, a model is specified as a recipe for computing the average of a polynomial of fields. We will use the following notation: given a polynomial \(F(\psi,\Psi,H)\) in the fields, its average will be denoted by \(\langle F(\psi,\Psi,H)\rangle\). The average is a linear operation, and, as indicated by the name "internal", only acts on internal fields, so that the external and super-external fields are viewed as constants for the averaging operation and can be factored out. It therefore suffices to define the average of a monomial of internal fields. +

+

+ First of all, a monomial of internal fields which does not contain the same number of \(\psi^+_i\) as \(\psi^-_j\) is \(0\). The average of a monomial of the form \(\psi_{i_1}^+\psi_{j_1}^-\cdots\psi_{i_n}^+\psi_{j_n}^-\) in which \(n\in\mathbb N\), \((i_1,\cdots,i_n)\in\{1,\cdots,I\}^n\) and \((j_1,\cdots,j_n)\in\{1,\cdots,I\}^n\), is computed using the Wick rule: + $$ + \langle\psi_{i_1}^+\psi_{j_1}^-\cdots\psi_{i_n}^+\psi_{j_n}^-\rangle= + \sum_{\pi\in\mathcal S_n}(-1)^\pi\langle\psi_{i_1}^+\psi_{j_{\pi(1)}}^-\rangle\cdots\langle\psi_{i_n}^+\psi_{j_{\pi(N)}}^-\rangle + $$ + in which \(\mathcal S_n\) denotes the set of permutations of \(\{1,\cdots,n\}\) and \((-1)^\pi\) denotes the signature of \(\pi\). Using the Wick rule, we can specify any average by defining the quadratic moments of the model, similarly to the moments of a 0-mean Gaussian measure. The collection of all quadratic moments of the form + $$\langle\psi_i^+\psi_j^-\rangle$$ + with \((i,j)\in\{1,\cdots,I\}^2\) is called the propagator of the model. +

+

+ In the configuration file of the meankondo program, the propagator is specified in the #!propagator entry. Note that meankondo recognizes numeric propagators as well as symbolic ones. +

+ +

Flow equation

+

+ In this section, we discuss what flow equations are, and how meankondo computes them. +

+ +

Definition

+

+ We first discuss how a flow equation is defined from a renormalization group map. The discussion below is not, in any sense, precise, and is meant as a guiding idea to understand why meankondo does what it does in the way it does. +

+

+ Consider a map, which we will call the renormalization group flow, of the form: + $$ + V^{[m]}(\Phi^{[m]},H)\longmapsto V^{[m-1]}(\Phi^{[m-1]},H)=\frac1{C^{[m]}}\left(\left\langle\prod_{\nu=1}^k\left(1+V^{[m]}(\Phi_\nu^{[m]},H)\right)\right\rangle-1\right) + $$ + where \(V^{[m]}\) anf \(V^{[m-1]}\) are polynomials with no constant term, \(C^{[m]}\in\mathbb R\setminus\{0\}\) is a constant, \(\langle\cdot\rangle\) is the average defined in Propagator, \(H\) is the collection of super-external fields, \(k=I/E\) which we assume to be an integer, \(\{\Phi^{[m]\pm}_{\nu,i}\}_{i\in\{1,\cdots,E\},\ k\in\{1,\cdots,k\}}\) and \(\{\Phi^{[m-1]\pm}_{i}\}_{i\in\{1,\cdots,E\}}\) are collections of combinations of the internal and external fields, of the form + $$ + \Phi^{[m]\pm}_{\nu,i}=\alpha_i\Psi_i^\pm+\psi_{j_\nu^{(i)}}^\pm,\quad + \Phi^{[m-1]\pm}_{i}=\Psi_i^\pm + $$ + in which \(\alpha_i\in\mathbb R\setminus\{0\}\) is some re-scaling factor, and + $$ + \bigcup_{i\in\{1,\cdots,E\}}^\circ\{j_1^{(i)},\cdots,j_k^{(i)}\}=\{1,\cdots,I\} + $$ + in which the \(\circ\) over the union means it is disjoint. +

+

+ If \(V^{[m-1]}\) can be written in the same form as \(V^{[m]}\), then the renormalization group map can be written as a finite system of equations: if there exist \(p\in\mathbb N\), \((\ell_1^{[m]},\cdots,\ell_p^{[m]})\in\mathbb C^p\), \((\ell_1^{[m-1]},\cdots,\ell_p^{[m-1]})\in\mathbb C^p\), and polynomials \((O_1(\Phi,H),\cdots,O_l(\Phi,H))\) such that + $$ + V^{[m]}(\Phi_\nu^{[m]},H)=\sum_{n=1}^p\ell_n^{[m]}O_n(\Phi_\nu^{[m]},H),\quad + V^{[m-1]}(\Phi^{[m-1]},H)=\sum_{n=1}^p\ell_n^{[m-1]}O_n(\Phi^{[m-1]},H) + $$ + then the renormalization group map can be written as a finite system of equations, called the flow equation, which maps + $$ + (\ell_1^{[m]},\cdots,\ell_p^{[m]})\longmapsto(C^{[m]},\ell_1^{[m-1]},\cdots,\ell_1^{[m-1]}) + $$ + (we added the constant \(C^{[m]}\) which plays an important role). The collection \(\underline\ell\) is called the collection of running coupling constants. +

+ +

Computation

+

+ We now discuss how meankondo computes flow equations. The two relevant configuration file entries are #!input_polynomial and #!id_table. +

+

+ The entry #!input_polynomial specifies a polynomial \(F(\underline\ell,\psi,\Psi,H)\) in the fields, as well as in a family of complex symbolic variables \(\underline\ell=(\ell_1,\cdots,\ell_p)\). In the context of the computation in Definition, it would be of the form + $$ + \prod_{\nu=1}^k\left(1+\sum_{n=1}^p\ell_nO_n(\Phi_{\nu}^{[m]},H)\right). + $$ +

+

+ The entry #!id_table specifies a sequence of \(p\) polynomials \((O_1(\Psi,H),\cdots,O_p(\Psi,H))\) in the external fields. In the context of the computation in Definition the \(n\)-th polynomial would be + $$ + O_n(\Phi_{\nu}^{[m]},H) + $$ + which specifies that \(O_n\) corresponds to \(\ell_n\). +

+

+ meankondo first computes the average + $$ + F'(\underline\ell,\Psi,H):=\langle F(\underline\ell,\psi,\Psi,H)\rangle + $$ + and then groups \(F'\) according to the id table, that is it tries to put it in the form + $$ + F'(\underline\ell,\Psi,H)=C(\underline\ell)\left(1+\sum_{n=1}^p\ell'_n(\underline\ell)O_n(\Psi,H)\right) + $$ + in which \(C(\underline\ell)\) is a complex constant term and \(\ell'_n(\underline\ell)\) is some complex coefficient, which are both functions of \(\underline\ell\). Note that for Fermionic hierarchical models, \(\ell'_n\) is a rational function of \(\underline\ell\), and \(C\) is a polynomial in \(\underline\ell\). In addition, \(\ell'_n\) can be expressed as a polynomial in \((\underline\ell,C^{-1}(\underline\ell))\). The flow equation is then defined as + $$ + \mathcal R:\underline\ell\longmapsto(C(\underline\ell),\ell'_1(\underline\ell),\cdots,\ell'_p(\underline\ell)) + $$ + If \(F'\) cannot be grouped, then meankondo exits with error code -1. +

+ +

Operations on flow equations

+

+ In this section we describe the various operations on flow equations that the tools bundled with meankondo support. +

+ +

Numerical evaluation

+

+ Once a flow equation \(\mathcal R\) has been computed, it can be numerically evaluated by passing a vector \(\underline\ell\) of double precision floating point numbers to meantools eval or numkondo which computes \(\mathcal R(\underline\ell)\) in the former case, and \(\mathcal R^m(\underline\ell)\) for any \(m\in\mathbb N\) in the latter. +

+ +

+ Numerical evaluation is handled in a straightforward manner, but for the following consideration. As was mentioned in Computation, \(\mathcal R\) is a polynomial in \((\underline\ell,C^{-1}(\underline\ell))\), and when evaluating \(\mathcal R(\underline\ell)\), meankondo first evaluates \(C\) and the computes \(\ell'_n(\underline\ell)\). +

+ +

Exponentiation

+

+ Oftentimes the renormalization group flow is expressed in terms of an exponential of an effective potential \(\exp(W)\), in which case the exponential must be computed before it can be processed by meankondo: + $$ + \exp(W)=1+V. + $$ + This is handled by meantools exp, which computes the running coupling constants appearing in \(V\) in terms of those in \(W\). +

+ +

Derivation

+

+ This feature was introduced to compute the susceptibility in the hierarchical Kondo model. In that case, some of the running coupling constants depend on the field, \(h\), and the susceptibility is expressed as a derivative of \(C(\underline\ell(h))\) with respect to \(h\). To that end, we wrote meantools derive to compute the derivatives of a flow equation with respect to an external variable. +

+ +

+ The input of meantools derive consists in a flow equation and a collection of variables \(X\subset\{1,\cdots,p\}\). Each running coupling constant \(\ell_i\) for \(i\in X\) is assumed to depend on an external parameter, \(h\). The flow equation is then derived with respect to \(h\): for every \(n\in\{1,\cdots,p\}\), the derivative of \(\ell_n'(\underline\ell)\) with respect to \(h\) in terms of \(\partial_h\ell_i\) for \(i\in X\) is computed. It is then appended to the input flow equation. +

+ +

Comments on the exactness of the computation

+

+ The computation of the flow equation, as well as its exponentiation and derivation, are exact in the sense that they only involve operations on integers and are not subject to truncations. The coefficients appearing in the flow equation are therefore exact. This statement has one major caveat: integer operations are only correct as long as the integers involved are not too large. The precise meaning of "not too large" is system dependent. In the source code, integers relating to flow equation coefficients are declared with the long int type, which, at least using the C library meankondo was tested with (that is glibc 2.21), means integers are encoded on 64 bits on 64-bit systems and 32 bits on 32-bit systems. All operations are therefore exact as long as all integers are in \([-2^{31},2^{31}-1]\) on 64-bit systems and \([-2^{15},2^{15}-1]\) on 32-bit systems. +

+ + +

+ Numerical evaluations are not exact. The numbers manipulated meankondo are "long doubles", which, when compiled for x86 processors, have a precision of 64 bits, which implies they are accurate to 19 decimal places; and the absolute value of doubles is bounded above by \(2^{16384}-2^{16384-64}\) (that is the number whose binary expansion has \(16383\) digits and whose \(64\) left-most digits are \(1\) whereas the others are \(0\)) and below by \(2^{-16382}\). +

+ + +

Authors

+

+ meankondo was written by Ian Jauslin, in the context of a project in collaboration with Giuseppe Benfatto and Giovanni Gallavotti. +

+ + diff --git a/man/kondo_preprocess.1 b/man/kondo_preprocess.1 new file mode 100644 index 0000000..cb37cd3 --- /dev/null +++ b/man/kondo_preprocess.1 @@ -0,0 +1,163 @@ +.Dd $Mdocdate: April 14 2015 $ +.Dt kondo_preprocess 1.2 +.Os +.Sh NAME +.Nm kondo_preprocess +.Nd A pre-processor to generate configuration files for +.Sy meankondo +for the Kondo model +.Sh SYNOPSIS +.Nm +.Op Fl d Ar dimension +.Op Ar config_file +.Pp +.Nm +.Fl v +.Sh DESCRIPTION +.Nm +generates a configuration file to be read by +.Sy meankondo +for the Kondo model. It generates the '#!fields', '#!symbols', '#!identities', '#!groups', '#!propagator', '#!input_polynomial' and '#!id_table' entries from special '#!propagator', '#!input_polynomial' and '#!id_table' entries, which are much more synthetic than those needed for the Kondo model. +.Pp +The quantities in the configuration file are expressed in terms of the observables A and B, which we do not define here, as well as the magnetic field h. +.Pp +.Nm +is part of a set of tools to compute and manipulate Fermionic hierarchical flows: +.Bl -bullet +.It +.Sy meankondo +: computes flow equations for hierarchical Fermionic models +.It +.Sy numkondo +: numerical evaluation of flow equations. +.It +.Sy meantools, meantools-convert +: perform various operations on flow equations (derivation, exponentiation, evaluation and conversion to other formats). +.El +.Pp +as well as the following pre-processors, which generate configuration files for their associated model: +.Bl -bullet +.It +.Sy kondo_proprocess +: Kondo model +.El +.Pp +.Sh COMMAND-LINE ARGUMENTS +.Bl -tag -width Ds +.It Fl d Ar dimension +The dimension of the field theory for the Kondo model (defaults to 2), including imaginary time (2 if the Fermionic chain is not neglected, 1 if it is). This parameter is used to determine how many boxes contribute to each scale in the hierarchical model: id the dimension is 2, then there are 4 boxes, whereas if it is 1, then there are 2. +.It Fl v +Print version information and exit. +.El +.Pp +.Sh CONFIGURATION FILE +.Nm +reads a configuration file, that can either be passed as a command line argument or to stdin, which specifies the model for which to compute the flow equation. +.Pp +A configuration file is a list of entries, separated by a '&' character, each of which has a title (or header), which is preceded by '#!'. Note that '#!' must be at the beginning of a line in order to be read correctly. +.Pp +Whenever the '#' character is encountered, the rest of the line is treated as a comment and ignored (unless it is followed by '!'). +.Pp +As a general rule, spaces and line breaks in the entries of the configuration file are ignored, so they may be used at the user's discretion. The few entries that require that no extra line breaks be inserted are explicitly mentioned below. +.Pp +.Nm +recognizes the following entries (unless explicitly mentioned, the entries below are mandatory) (entries may be provided in any order) (any extra entries in the configuration file are ignored): +.Bl -tag -width Ds +.It Sy #!input_polynomial +The polynomial whose mean we wish to compute in order to calculate the flow equation. +.Pp +The format of the polynomial is similar to that in +.Sx meankondo Ns (1) , +up to the following differences. +.Bl -bullet +.It +The fields can be specified as scalar products of A's and B's. For each n in {1,...,dimension}, +.Nm +defines An and Bn, as well as symbols for scalar products of the form +.D1 [f An.An] +.D1 [f Bn.Bn] +.D1 [f An.Bn] +.D1 [f An.h] +.D1 [f Bn.h] +In addition, a vector product symbol is defined for (AnxBn).h : +.D1 [f AnxBn.h] +.Pp +.It +In addition, +.Nm +defines external fields for A and B, denoted by a and b. They can be used as fields in the input polynomial using the syntax +.D1 +.D1 +.D1 +.D1 +.D1 +.D1 +.Pp +.It +Furthermore, in order to simplify writing products of polynomials over each box index, if the polynomial contains a '%', then +.Nm +multiplies the polynomial by itself as many times as there are boxes (2^dimension times), replacing '%' with the appropriate box index. For example, if dimension=1 +.D1 '[fA%.A%]+[fB%.B%]' +is equivalent to +.D1 '[fA1.A1]+[fB1.B1] * [fA2.A2]+[fB2.B2]'. +.El +.Pp +Example: +.D1 (1) + (1/2)[l1][fA1.A1] + (1/2)[l2][fB1.h] +.D1 * (1) + (1/2)[l2][fA2.A2] +(1/2)[l2][fB2.h] +.Pp +.It Sy #!id_table +The idtable used to identify the running coupling constants. +.Pp +The idtable has the same syntax as that in +.Sx meankondo Ns (1) , +in which the polynomial can use the fields +.D1 +.D1 +.D1 +.D1 +.D1 +.D1 +defined above. +.Pp +Example: +.D1 1:(1/2), 2:(2) +.Pp +.It Sy #!propagator +The propagator of the model. +.Pp +The propagator syntax differs from +.Sx meankondo Ns (1) , +in that the field indices are specified using An and Bn. +.Pp +Example: +.D1 A1;A2: 1 , A2;A1: -1 , B1;B2: s{-1} , B2;B1: (-1)s{-1} +.Pp +.It Sy extra entries +If there is a '#!symbols' or an '#!identities' entry in the configuration file, then they are appended to the end of those entries in the new configuration file. +.Pp +Any other entry is appended to the new configuration file. This can be useful to pipe the output to tools other than +.Sy meankondo +(e.g. +.Sy meantools ) . +.Pp +.Sh OUTPUT +.Nm +prints the configuration file to stdout. +.Pp +The output of +.Nm +can be piped into +.Sy meankondo +directly. +.Pp +.Sh RETURN CODE +.Nm +returns 0 on success and -1 on error. +.Pp +.Sh SEE ALSO +.Sx meankondo Ns (1) , +.Sx numkondo Ns (1) , +.Sx meantools Ns (1) , +.Sx meantools-convert Ns (1) +.Pp diff --git a/man/meankondo.1 b/man/meankondo.1 new file mode 100644 index 0000000..ea4814b --- /dev/null +++ b/man/meankondo.1 @@ -0,0 +1,221 @@ +.Dd $Mdocdate: April 13 2015 $ +.Dt meankondo 1.2 +.Os +.Sh NAME +.Nm meankondo +.Nd A tool to compute renormalization group flows for Fermionic hierarchical models +.Sh SYNOPSIS +.Nm +.Op Fl t Ar threads +.Op Fl C +.Op Ar config_file +.Pp +.Nm +.Fl v +.Sh DESCRIPTION +.Nm +computes the renormalization group flow equations for Fermionic hierarchical models, which should be defined in the configuration file provided on the command line, following the syntax detailed below. +.Pp +The flow equation is computed by calculating the mean of a polynomial of fields and running coupling constants using the Wick rule and a propagator provided in the configuration file. The running coupling constants are then identified in the resulting polynomial using an idtable provided in the configuration file. +.Pp +.Nm +is part of a set of tools to compute and manipulate Fermionic hierarchical flows: +.Bl -bullet +.It +.Sy meankondo +: computes flow equations for hierarchical Fermionic models +.It +.Sy numkondo +: numerical evaluation of flow equations. +.It +.Sy meantools, meantools-convert +: perform various operations on flow equations (derivation, exponentiation, evaluation and conversion to other formats). +.El +.Pp +as well as the following pre-processors, which generate configuration files for their associated model: +.Bl -bullet +.It +.Sy kondo_proprocess +: Kondo model +.El +.Pp +.Sh COMMAND-LINE ARGUMENTS +.Bl -tag -width Ds +.It Fl t Ar threads +The number of threads to use for the computation. +.It Fl C +Format the ouptput so it can be piped to +.Sy numkondo , +that is, instead of printing the flow equation, print a full configuration file containing the flow equation as well as all the other entries of the configuration file that do not pertain to the computation of the flow equation. +.It Fl v +Print version information and exit. +.El +.Pp +.Sh CONFIGURATION FILE +.Nm +reads a configuration file, that can either be passed as a command line argument or to stdin, which specifies the model for which to compute the flow equation. +.Pp +A configuration file is a list of entries, separated by a '&' character, each of which has a title (or header), which is preceded by '#!'. Note that '#!' must be at the beginning of a line in order to be read correctly. +.Pp +Whenever the '#' character is encountered, the rest of the line is treated as a comment and ignored (unless it is followed by '!'). +.Pp +As a general rule, spaces and line breaks in the entries of the configuration file are ignored, so they may be used at the user's discretion. The few entries that require that no extra line breaks be inserted are explicitly mentioned below. +.Pp +.Nm +recognizes the following entries (unless explicitly mentioned, the entries below are mandatory) (entries may be provided in any order) (any extra entries in the configuration file are ignored): +.Bl -tag -width Ds +.It Sy #!fields +A list of the fields of the model. +.Pp +The fields entry contains 4 lines which start with 'i:', 'x:', 'h:' and 'f:'. Each of these is followed by a ',' separated list of field indices, which are positive integers. +.Bl -bullet +.It +The indices following 'i' correspond to internal fields, which are integrated out using the Wick rule and the propagator provided in the '#!propagator' entry. Each internal field is associated a conjugate field, whose index is the opposite of the field's index (e.g. 'i:101' defines a field whose index is -101) +.It +The indices following 'x' correspond to external fields that are associated conjugate field (e.g. 'x:100' defines a field whose index is -100). External indices may not appear as internal indices. +.It +The indices following 'h' correspond to external fields that are not associated a conjugate field. External indices may not appear as internal indices. +.It +The 'f' line specifies which of the internal and external indices are Fermions, i.e. which fields anti-commute. The fields appearing in the 'f' line should also either appear in the 'i' or 'x' line. WARNING: for the moment, only cases in which all of the internal fields are Fermions are supported. +.El +.Pp +.Em Line breaks are not ignored in this entry. +.Pp +Example: +.D1 i:101,102,201,202 +.D1 x:100,200 +.D1 h:301,302,303 +.D1 f:100,101,102 +.It Sy #!propagator +The propagator of the model. +.Pp +The propagator entry is a ',' separated list whose elements are of the form +.D1 index1;index2: polynomial +where index1 and index2 are internal indices, and polynomial is a polynomial (see the POLYNOMIALS section below for information on how to format polynomials). The polynomial must not depend on the internal fields. Note that a number is a special type of polynomial, so propagators with numerical entries are handled by +.Nm +just as easily as propagators with symbolic entries. Such an entry means that +.D1 = polynomial. +.Pp +Example: +.D1 101;102: 1 , 102;101: -1 , 201;202: s{-1} + (-1)[l10] , 202;201: (-1)s{-1} + [l10] +.It Sy #!symbols +Symbolic variables used as shortcuts for more complicated expressions (optional entry). +.Pp +In order to simplify long expressions, symbolic variables can be defined in this entry. Each variable is assigned an index, which is a positive integer that must be different from any of the internel and external indices defined in the '#!fields' entry. +.Pp +The symbols entry is a ',' separated list, whose elements are of the form +.D1 index= polynomial +where index is the index of the variable and polynomial is the expression it stands for (see the POLYNOMIALS section below for information on how to format polynomials). Note that polynomial can contain other symbolic variables. There is no safeguard against self-referencing definitions that may cause infinite loops. +.Pp +Example: +.D1 1001= (-1)[f-100][f100] + (-1)[f-101][f101] , 2001=[f-100][f100] + [f-201][f201] +.Pp +This entry is optional. +.Pp +.It Sy #!identities +Identities satisfied by some of the fields (optional entry). +.Pp +In some cases, some of the quantities involved in a model will satisfy an identity (e.g. a vector may be of unit-norm), which should simplified out from the flow equation. +.Pp +The identities entry is a ',' separated list, whose elements are of the form +.D1 monomial=polynomial +where monomial represents the left side of the identity and is a sequence of field indices of the form '[f index1][f index2]...' and polynomial represents the right side of the identity (see the POLYNOMIALS section below for information on how to format polynomials). +.Pp +Example: +.D1 [f301][f301]=(1)+(-1)[f302][f302]+(-1)[f303][f303] +.Pp +This entry is optional. +.Pp +.It Sy #!input_polynomial +The polynomial whose mean we wish to compute in order to calculate the flow equation. +.Pp +The format of the polynomial is that specified in the POLYNOMIALS section. In addition, the polynomial can be specified as the product of other polynomials: +.D1 polynomial1 * polynomial2 * ... +Note that there are no parentheses, and therefore, products cannot be nested, nor can a product of polynomials be summed with another polynomial. +.Pp +Example: +.D1 (1) + (1/2)[l1][f1001] * (1) + (1/2)[l2][f2001] +.It Sy #!id_table +The idtable used to identify the running coupling constants. +.Pp +Once the mean of the input polynomial has been computed, we are left with a polynomial of the external fields and the running coupling constants that were in the input polynomial. In order to compute a flow equation from this average, +.Nm +uses an idtable to identify which of the monomials of the average contribute to which running coupling constant. +.Pp +The id_table entry is a ',' separated list, whose elements are of the form +.D1 rcc: polynomial +where rcc is the index of the corresponding running coupling constant, which is a non-negative integer, and polynomial is the polynomial to which rcc refers to (which is a polynomial of the external fields). +.Pp +Example: +.D1 1:(-1)[f-100][f100] , 2:[f-200][f200] +.It Sy #!groups +Groups of independent variables (optional entry). +.Pp +In order to speed up the computation of the mean of the input polynomial, groups of independent variables can be specified. When +.Nm +computes the mean of a monomial containing elements of different groups, it factors the monomial into independent factors, computes the mean of each factor, and then takes their product. This way, +.Nm +does not repeatedly try to pair independent fields. +.Pp +The groups entry is a list of collections of fields or symbols of the following form +.D1 (index1,index2,...) +.Pp +Example: +.D1 (1001,1002) (2001,2002) +.Pp +.Em Warning: +.Nm +does not check that the fields in different groups are truly independent, so cases in which fields in different group have a non-vanishing propagator entry may give unexpected results. +.Pp +This entry is optional. +.El +.Pp +.Sh NUMBERS +.Nm +can parse rational numbers and linear combinations of square roots of integers (positive or negative (which is how complex numbers are implemented)) with rational coefficients (i.e. elements of the field extension of Q generated by sqrt(Z)). +.Pp +A number is a '+' separated list whose elements are of the form +.D1 (a/b)s{r} +where a and r are integers and b is a positive integer. s{r} stands for sqrt(r). +.Pp +If a=b, then the number may be written as 's{r}'. If b=1, then it can be '(a)s{r}'. If r=1, then it can be 'a/b'. If b=r=1, then it can be 'a'. +.Pp +Example: +.D1 (1/2)s{2} + (-1)s{-1} + 3/2 +.Pp +.Sh POLYNOMIALS +Polynomials are '+' separated lists of monomials. Each monomial is a sequence of numbers, rccs and fields. +.Bl -bullet +.It +Numbers are enclosed between '(' and ')'. If there are several numbers in a monomial, then they are multiplied. +.It +rccs are non-negative indices enclosed between '[l' and ']'. +.It +Fields are non-vanishing indices enclosed between '[f' and ']'. Fields must either appear in the '#!fields' entry or the '#!symbols' entry. +.El +.Pp +If the numerical factor of a monomial is 1, then it can be dropped. However, even if the numerical factor is a single integer, its '(' and ')' delimiters cannot be omitted. +.Pp +Example: +.D1 (1) + ((3/2)s{2} + (-1)s{-1} + 3)[l1][l2][f100][f1001][f101] + [l1][f101] + (3)[l2] +.Pp +.Sh OUTPUT +.Nm +prints the flow equation to stdout. +.Pp +The rccs in the flow equation are of the form '[% index]' and the constant term in the flow equation is denoted by '[C1]'. The factor '[/C1^power]' stands for (1/C1^power). +.Pp +.Sh KNOWN ISSUES +.Nm +only supports models in which all of the internal fields are Fermions. +.Pp +.Sh RETURN CODE +.Nm +returns 0 on success and -1 on error. +.Pp +.Sh SEE ALSO +.Sx numkondo Ns (1) , +.Sx meantools Ns (1) , +.Sx meantools-convert Ns (1) , +.Sx kondo_preprocess Ns (1) +.Pp diff --git a/man/meantools-convert.1 b/man/meantools-convert.1 new file mode 100644 index 0000000..f29d877 --- /dev/null +++ b/man/meantools-convert.1 @@ -0,0 +1,137 @@ +.Dd $Mdocdate: June 12 2015 $ +.Dt meantools-convert 1.2 +.Os +.Sh NAME +.Nm meantools-convert +.Nd Translate flow equations to various formats +.Sh SYNOPSIS +.Nm +.Op Fl i Ar file +.Sy C +.Op Fl l Ar oldl_symbol +.Op Fl L Ar newl_symbol +.Op Fl C Ar constant_symbol +.Op Fl S +.Pp +.Nm +.Op Fl i Ar file +.Sy javascript +.Op Fl l Ar oldl_symbol +.Op Fl L Ar newl_symbol +.Op Fl C Ar constant_symbol +.Pp +.Nm +.Op Fl i Ar file +.Sy LaTeX +.Op Fl l Ar oldl_symbol +.Op Fl L Ar newl_symbol +.Op Fl C Ar constant_symbol +.Op Fl s +.Pp +.Sh DESCRIPTION +.Nm +is a python script to convert flow equations to C, javascript or LaTeX code. +.Pp +.Nm +is part of a set of tools to compute and manipulate Fermionic hierarchical flows: +.Bl -bullet +.It +.Sy meankondo +: computes flow equations for hierarchical Fermionic models +.It +.Sy numkondo +: numerical evaluation of flow equations. +.It +.Sy meantools, meantools-convert +: perform various operations on flow equations (derivation, exponentiation, evaluation and conversion to other formats). +.El +.Pp +as well as the following pre-processors, which generate configuration files for their associated model: +.Bl -bullet +.It +.Sy kondo_proprocess +: Kondo model +.El +.Pp +.Sh COMMAND-LINE ARGUMENTS +.Bl -tag -width Ds +.It Fl i Ar file +File containing the flow equation to convert. If this flag is not specified, then +.Nm +reads the equation from stdin. +.El +.Pp +.Sh C +.Nm +supports the following flags when converting to C code. +.Bl -tag -width Ds +.It Fl l Ar oldl_symbol +The symbol for the running coupling constants appearing on the right hand side (defaults to 'l'). +.It Fl L Ar newl_symbol +The symbol for the running coupling constants appearing on the left hand side (defaults to 'newL'). +.It Fl C Ar constant_symbol +The symbol for the constant factor (defaults to 'C'). Since there can be multiple constants in a flow equation, the actual symbols will be appended an integer (e.g. C1, C2, etc). +.It Fl S +Write each equation on a single line. +.El +.Pp +.Ar oldl_symbol +and +.Ar newl_symbol +must be arrays of floating point numbers (floats, doubles or long doubles) with enough memory allocated to them to hold all running coupling constants. Similarly, +.Ar constant_symbol Ns n +with n large enough must be allocated as floating point numbers. +.Pp +In addition, if the flow equation contains derivatives (see +.Sx meantools Ns (1) ) , +then +.Pf d Ar oldl_symbol , +.Pf d Ar newl_symbol +and +.Pf d Ar constant_symbol Ns n +must be allocated as well. The same holds for higher derivatives, with 'd' replaced by 'dd', 'ddd' and so forth. +.Sh JAVASCRIPT +.Qo +.Nm +.Sy javascript +.Qc +is a shorthand for +.Qo +.Nm +.Sy C +-S +.Qc +.Sh LATEX +.Nm +supports the following flags when converting to LaTeX code. +.Bl -tag -width Ds +.It Fl l Ar oldl_symbol +The symbol for the running coupling constants appearing on the right hand side (defaults to "\\ell"). +.It Fl L Ar newl_symbol +The symbol for the running coupling constants appearing on the left hand side (defaults to "\\ell'"). +.It Fl C Ar constant_symbol +The symbol for the constant factor (defaults to "C"). Since there can be multiple constants in a flow equation, the actual symbols will be appended an integer subscript (e.g. C_1, C_2, etc). +.It Fl s +Break the line at every '+'. +.El +.Pp +Since the length of the equations can change a lot from one system to another, +.Nm +does not make many formatting decisions (for instance +.Nm +does not insert any alignment characters, and either breaks the line at every '+' sign or at the end of the equation), which are left to the user. +.Pp +.Sh OUTPUT +.Nm +prints the converted flow equation to stdout. +.Pp +.Sh RETURN CODE +.Nm +returns 0 on success and -1 on error. +.Pp +.Sh SEE ALSO +.Sx meankondo Ns (1) +.Sx numkondo Ns (1) , +.Sx meantools Ns (1) , +.Sx kondo_preprocess Ns (1) +.Pp diff --git a/man/meantools.1 b/man/meantools.1 new file mode 100644 index 0000000..1849326 --- /dev/null +++ b/man/meantools.1 @@ -0,0 +1,164 @@ +.Dd $Mdocdate: April 14 2015 $ +.Dt meantools 1.2 +.Os +.Sh NAME +.Nm meantools +.Nd A tool to manipulate flow equations +.Sh SYNOPSIS +.Nm +.Sy exp +.Op Ar config_file +.Pp +.Nm +.Sy derive +.Op Fl d Ar nderivs +.Op Fl V Ar variables +.Op Ar config_file +.Pp +.Nm +.Sy eval +.Op Fl R Ar values +.Op Ar config_file +.Pp +.Sh DESCRIPTION +.Nm +performs various operations on flow equations generated by +.Sy meankondo. +Namely, it can exponentiate, derive and evaluate flow equations. +.Pp +.Nm +is part of a set of tools to compute and manipulate Fermionic hierarchical flows: +.Bl -bullet +.It +.Sy meankondo +: computes flow equations for hierarchical Fermionic models +.It +.Sy numkondo +: numerical evaluation of flow equations. +.It +.Sy meantools, meantools-convert +: perform various operations on flow equations (derivation, exponentiation, evaluation and conversion to other formats). +.El +.Pp +as well as the following pre-processors, which generate configuration files for their associated model: +.Bl -bullet +.It +.Sy kondo_proprocess +: Kondo model +.El +.Pp +.Sh EXP +When run with the 'exp' command, +.Nm +computes the exponential of a flow equation. All the required parameters are set in the configuration file, which it either reads from the file provided on the command line, or from stdin. +.Pp +The syntax for the configuration file is the same as for +.Sx meankondo Ns (1) , +and will not be belaboured here. The supported entries are +.Bl -tag -width Ds +.It Sy #!input_polynomial +The polynomial whose exponential is to be computed. +.Pp +.It Sy #!fields +The fields appearing in the polynomial +.Pp +.It Sy #!symbols +Symbolic variables (optional entry). +.Pp +.It Sy #!identities +identities between fields (optional entry). +.Pp +.It Sy #!id_table +The idtable used to compute a flow equation from the polynomial. +.El +.Pp +The resulting flow equation is written to stdout. +.Pp +.Sh DERIVE +When run with the 'derive' command, +.Nm +computes derivatives of a flow equation provided in the configuration file, which can either be passed as a command-line argument or through stdin. +.Pp +The derivatives are derivatives with respect to an extra virtual parameter, which all of the rccs are assumed to depend on (to override the default behavior, the '-V' flag can be used to pass a list of rccs that depend on the extra parameter, alternatively such a list can be given in the configuration file). The derivative of the flow equation is a new flow equation for the rccs and their derivatives with respect to the virtual parameter. +.Pp +When multiple derivatives are taken, the flow equation becomes a flow equation for the rccs, their derivatives, second derivatives, and so forth... +.Pp +This operation can be useful, for instance, to compute moments in an interacting system, in which the generating functional can be expressed as an effective potential depending on a parameter with respect to which the result of the integration should be derived. The 'derive' command writes the flow equation for the derived rccs, from which the quantities of interest can be computed. +.Pp +.Sy Command-line arguments: +.Bl -tag -width Ds +.It Fl d Ar nderivs +Number of derivatives (defaults to 1) +.It Fl V Ar variables +The variables that depend on the extra virtual parameter (defaults to all) (WARNING: if one of the variables has a negative index, do not put it first in the list, since +.Nm +would interpret the argument as being a flag, for example, write '-V "0,-1"' instead of '-V "-1,0"'). +.Pp +Can either be a ',' separated list if indices or 'all' to derive with respect to all available variables. +.El +.Pp +.Sy Configuration file: +.Pp +The configuration file contains the flow equation to derive, and optionally a list of variables (similar to the '-V' flag). The following entries are supported: +.Bl -tag -width Ds +.It Sy #!flow_equation +The flow equation to derive. +.Pp +The syntax is identical to that in +.Sx numkondo Ns (1) . +.Pp +If this entry is the only one in the configuration file, the '#!flow_equation' header may be omitted. +.Pp +.It Sy #!variables +The variables that depend on the extra virtual parameter (optional entry). +.Pp +The variables entry is a ',' separated list of indices, or 'all' in which case, all of the variables on the right side of the flow equation are assumed to depend on the flow equation. +.Pp +If the '-V' flag is provided on the command-line, this entry is ignored. +.El +.Pp +The resulting flow equation is written to stdout. +.Pp +.Sh EVAL +When run with the 'eval' command, +.Nm +evaluates a flow equation, provided in a configuration file, numerically, using the values provided on the command-line or in the configuration file provided on the command-line or through stdin. +.Pp +.Sy Command-line arguments: +.Bl -tag -width Ds +.It Fl R Ar values +The values of the rccs with which to evaluate the flow equation. +.Ar values +is formatted like an initial_condition (see +.Sx numkondo Ns (1) ) . +.El +.Pp +.Sy Configuration file: +.Pp +The configuration file contains the flow equation to evaluate, and optionally a list of values for the rccs. The following entries are supported: +.Bl -tag -width Ds +.It Sy #!flow_equation +The flow equation to evaluate. +.Pp +The syntax is identical to that in +.Sx numkondo Ns (1) . +.Pp +If this entry is the only one in the configuration file, the '#!flow_equation' header may be omitted. +.Pp +.It Sy #!initial_condition +The value on which to evaluate the flow equation (optional entry). +.Pp +The syntax is identical to that in +.Sx numkondo Ns (1) . +.Pp +If the '-R' flag is provided on the command-line, this entry is ignored. +.El +.Pp +The result of the evaluation is written to stdout, and is formatted is such a way that it can be used as an initial condition for +.Pp +.Sh SEE ALSO +.Sx meankondo Ns (1) , +.Sx numkondo Ns (1) , +.Sx meantools-convert Ns (1) , +.Sx kondo_preprocess Ns (1) +.Pp diff --git a/man/numkondo.1 b/man/numkondo.1 new file mode 100644 index 0000000..7406a80 --- /dev/null +++ b/man/numkondo.1 @@ -0,0 +1,157 @@ +.Dd $Mdocdate: April 14 2015 $ +.Dt numkondo 1.2 +.Os +.Sh NAME +.Nm numkondo +.Nd A tool to iterate a flow equation numerically. +.Sh SYNOPSIS +.Nm +.Op Fl F +.Op Fl N Ar niter +.Op Fl D Ar tolerance +.Op Fl I Ar initial_condition +.Op Ar config_file +.Pp +.Nm +.Fl v +.Sh DESCRIPTION +.Nm +computes the numerical value of the iterates of a flow equation provided in its configuration file, starting from an initial condition provided in its configuration file. +.Nm +is part of a set of tools to compute and manipulate Fermionic hierarchical flows: +.Bl -bullet +.It +.Sy meankondo +: computes flow equations for hierarchical Fermionic models +.It +.Sy numkondo +: numerical evaluation of flow equations. +.It +.Sy meantools, meantools-convert +: perform various operations on flow equations (derivation, exponentiation, evaluation and conversion to other formats). +.El +.Pp +as well as the following pre-processors, which generate configuration files for their associated model: +.Bl -bullet +.It +.Sy kondo_proprocess +: Kondo model +.El +.Pp +.Sh COMMAND-LINE ARGUMENTS +.Bl -tag -width Ds +.It Fl N Ar niter +Number of iterations +.It Fl F +Only print the last step of the computation, with full precision. The output can be used as an initial condition for further iterations. +.It Fl D Ar tolerance +If this option is provided, any number smaller than +.Ar tolerance +is set to 0. +.It Fl I Ar initial_condition +Set the initial condition from the command-line (overrides the initial condition in the configuration file). The format is the same as the '#!initial_configuration' entry, see below. +.It Fl v +Print version information and exit. +.El +.Pp +.Sh CONFIGURATION FILE +.Nm +reads a configuration file, that can either be passed as a command line argument or to stdin, which specifies the flow equation to iterate, the initial condition, and text labels for each running coupling constant. +.Pp +A configuration file is a list of entries, separated by a '&' character, each of which has a title (or header), which is preceded by '#!'. Note that '#!' must be at the beginning of a line in order to be read correctly. +.Pp +Whenever the '#' character is encountered, the rest of the line is treated as a comment and ignored (unless it is followed by '!'). +.Pp +As a general rule, spaces and line breaks in the entries of the configuration file are ignored, so they may be used at the user's discretion. The few entries that require that no extra line breaks be inserted are explicitly mentioned below. +.Pp +.Nm +recognizes the following entries (unless explicitly mentioned, the entries below are mandatory) (entries may be provided in any order) (any extra entries in the configuration file are ignored): +.Bl -tag -width Ds +.It Sy #!flow_equation +The flow equation to be iterated. +.Pp +The flow equation has the same format as the output of +.Sy meankondo. +That is, a flow equation is a ',' separated list whose elements are of the form +.D1 [% index] = equation +where index is a non-negative integer or +.D1 [C index] = equation +where index is is a positive integer; and equation is formatted as explained below. [% index] stands for the running coupling constant corresponding to index, and [C index] is a special type of running coupling constant, which correspond to terms in the effective potential that do not depend on the fields. The difference between [% ...] and [C ...] is that the latter are evaluated first, and may be used in the equation of [% ...]. +.Pp +Equations are '+' separated lists whose elements are monomials. Each monomial is a sequence of numbers, rccs and denominators: +.Bl -bullet +.It +Numbers are enclosed between '(' and ')'. If there are several numbers in a monomial, then they are multiplied. +.It +rccs are non-negative indices enclosed between '[%' and ']'. +.It +Denominators are of the form [/C index^power] and correspond to 1/C_index^power. Note that the denominators [C index] are computed before they are divided, and therefore, the equation corresponding to [C index] may not contain any [/C index^power]. +.El +.Pp +If the numerical factor of a monomial is 1, then it can be dropped. However, even if the numerical factor is a single integer, its '(' and ')' delimiters cannot be omitted. +.Pp +In addition, in order to deal with derivatives of flow equations, extra running coupling constants can be introduced by adding any number of 'd' before '%' and 'C'. For instance, a flow equation may contain +.D1 [dd% index] = ... +or +.D1 [dC index] = ... +and an equation may contain an rcc of the form +.D1 [d% index] +.Pp +In practice, [C index] is represented as [%-index], and derivatives offset an index by 1000000, so [d% index] is equivalent to [%1000000+index], [dC index] to [%-1000000-index], [dd% index] to [%2000000+index], and so forth... Indices must therefore be smaller than 1000000 if 'd' is used. +.Pp +Example: +.D1 [C1] = [%1] + (1/2)[%1][%2] , +.D1 [%1] = [%1][/C1^1] + [%1][%1][%2][/C1^2] , +.D1 [%2] = [%2][/C1^1] +and with a derivative: +.D1 [C1] = [%1] + (1/2)[%1][%2] , +.D1 [%1] = [%1][/C1^1] + [%1][%1][%2][/C1^2] , +.D1 [%2] = [%2][/C1^1] , +.D1 [dC1] = [d%1] + (1/2)[d%1][%2] + (1/2)[%1][d%2], +.D1 [d%1] = [d%1][/C1^1] + (-1)[%1][dC1][/C1^2] + (2)[%1][d%1][%2][/C1^2] + [%1][%1][d%2][/C1^2] + (-2) [%1][%1][%2][dC1][C1^3], +.D1 [d%2] = [d%2][/C1^1] + (-1)[%2][dC1][/C1^2] +.Pp +.It Sy #!initial_condition +The initial condition for the iteration. +.Pp +The initial_condition entry is a ',' separated list whose elements are of the form +.D1 index:value +where index is that of the corresponding rcc and value is a double precision float. The index may contain 'C' and 'd' characters as in the '#!flow_equation' entry. +.Pp +Example: +.D1 1:1.0 , 2:2.3e-6 , d1:2.0 +.Pp +Note that if +.Nm +is called with the '-F' flag, the corresponding output is formatted in such a way that it can be used as an initial condition for other iterations. +.Pp +.It Sy #!labels +Labels for the running coupling constants, used as headers to display the flow of the iteration. +.Pp +The labels entry is a ',' separated list whose elements are of the form +.D1 index:"label" +where index is the non-negative index of the corresponding rcc, and label is a string. +.Pp +Example: +.D1 1:"one" , 2:"two" +.El +.Pp +.Sh OUTPUT +Unless the '-F' flag is provided, +.Nm +prints the result of the iteration at each step to stdout. +.Pp +If the '-F' flag is provided, +.Nm +prints the last step of the iteration to stdout in a format that can be re-used as an initial condition for subsequent iterations. +.Pp +.Sh RETURN CODE +.Nm +returns 0 on success and -1 on error. +.Pp +.Sh SEE ALSO +.Sx meankondo Ns (1) , +.Sx meantools Ns (1) , +.Sx meantools-convert Ns (1) , +.Sx kondo_preprocess Ns (1) +.Pp diff --git a/scripts/config_functions.sh b/scripts/config_functions.sh new file mode 100644 index 0000000..593ca94 --- /dev/null +++ b/scripts/config_functions.sh @@ -0,0 +1,91 @@ +## Copyright 2015 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +################################################################################################ +## ## +## Functions used to manipulate configuration files ## +## ## +################################################################################################ + + +# remove unwanted entries in configuration files (read fromm stdin) +function strip_config { + tostrip="$1" + + # whether to skip the entry + skip=0 + + # whether this is the first entry + first_entry=1 + + # read file + while read aline; do + # if line is a header + if [[ "$aline" =~ ^#! ]]; then + # check whether to strip + for header in $tostrip; do + if [[ "$aline" = "$header" ]]; then + skip=1 + break + fi + done + + # add a '&' separator + if [[ "$skip" = 0 && "$first_entry" = 0 ]]; then + echo "&" + fi + [[ "$skip" = 0 ]] && first_entry=0 + fi + + # change entry + if [[ "$aline" = "&" ]]; then + skip=0 + # write out if not instructed to skip + elif [[ $skip = 0 ]]; then + echo "$aline" + fi + + done +} + +# find a config entry in a file (read from stdin) +function find_config_entry { + header="$1" + + # whether the entry was found + foundit=0 + + # read stdin + while read aline; do + # write out if the header was found + if [[ $foundit = 1 ]]; then + if [[ "$aline" != "&" ]]; then + echo "$aline" + fi + fi + + # check the header + if [[ "$aline" = "$header" ]]; then + foundit=1 + fi + + # new entry + if [[ "$aline" = "&" ]]; then + if [[ "$foundit" = 1 ]]; then + break + fi + fi + done +} + diff --git a/scripts/meantools-convert b/scripts/meantools-convert new file mode 100755 index 0000000..603749e --- /dev/null +++ b/scripts/meantools-convert @@ -0,0 +1,261 @@ +#!/usr/bin/env python + +## Copyright 2015 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + + +########################################################################## +## ## +## convert flow equations to various formats ## +## ## +########################################################################## + +import sys +import re + +def print_usage(): + print("\nusage:",file=sys.stderr) + print(" export_flow_equation [-i infile] [options]",file=sys.stderr) + print("\nwhere [options] is one of",file=sys.stderr) + print(" C [-l oldl_symbol] [-L newl_symbol] [-C Constant_symbol] [-S]",file=sys.stderr) + print(" javascript [-l oldl_symbol] [-L newl_symbol] [-C Constant_symbol]",file=sys.stderr) + print(" LaTeX [-l oldl_symbol] [-L newl_symbol] [-C Constant_symbol] [-s]",file=sys.stderr) + print("",file=sys.stderr) + +# to C +def C_engine(argv,text): + argc=len(argv) + + i=1 + # defaults + lsym="l" + Lsym="newL" + Csym="C" + oneline=0 + while (i=argc-1): + print("error: '"+argv[i]+"' must be followed by a string",file=sys.stderr) + exit(-1) + + # symbol for l + if (argv[i]=="-l"): + lsym=argv[i+1] + i=i+2 + # symbol for new l + elif (argv[i]=="-L"): + Lsym=argv[i+1] + i=i+2 + # symbol for C + elif (argv[i]=="-C"): + Csym=argv[i+1] + i=i+2 + + # write each equation on a single line + elif (argv[i]=="-S"): + oneline=1 + i=i+1 + + return(convert_C(text,lsym,Lsym,Csym,oneline)) + +# to LaTeX +def latex_engine(argv,text): + argc=len(argv) + + i=1 + # defaults + lsym="\\ell" + Lsym="\\ell'" + Csym="C" + oneline=1 + while (i=argc-1): + print("error: '"+argv[i]+"' must be followed by a string",file=sys.stderr) + exit(-1) + + # symbol for l + if (argv[i]=="-l"): + lsym=argv[i+1] + i=i+2 + # symbol for new l + elif (argv[i]=="-L"): + Lsym=argv[i+1] + i=i+2 + # symbol for C + elif (argv[i]=="-C"): + Csym=argv[i+1] + i=i+2 + + # break lines in equation + elif (argv[i]=="-s"): + oneline=0 + i=i+1 + + return(convert_latex(text,lsym,Lsym,Csym,oneline,columns)) + +# convert to C format +def convert_C(text, lsym, Lsym, Csym, oneline): + # make sure quotients are floats + text=re.sub(r'([0-9])/',r'\1./',text) + + # remove newlines + if (oneline==0): + text=text.replace('\n','\\\n') + # end of lines + text=text.replace(',\\',';') + else: + text=text.replace('\n','') + # end of lines + text=text.replace(',',';\n') + + # remove extra space + text=re.sub(r' {2,}',' ',text) + + # replace left hand side variables + text=re.sub(r'\[(d*)C *([0-9]*)\] =',r'\1'+Csym+r'\2 =', text) + text=re.sub(r'\[(d*)% *([0-9]*)\] =',r'\1'+Lsym+r'[\2] =', text) + # replace right hand side variables + text=re.sub(r'\[(d*)% *([0-9]*)\]',r'*\1'+lsym+r'[\2]',text) + text=re.sub(r'\[(d*)C *([0-9]*)\]',r'*\1'+Csym+r'[\2]',text) + # replace constants in right hand side (and expand powers as products) + match=re.search(r'\[/C([0-9]*)\^([0-9]*)\]',text) + while(match): + power=int(match.group(2)) + # product of constants + prod='' + for i in range(power): + prod=prod+'/'+Csym+match.group(1) + text=text.replace('[/C'+match.group(1)+'^'+match.group(2)+']',prod) + match=re.search(r'\[/C([0-9]*)\^([0-9]*)\]',text) + + # replace square roots + text=re.sub(r's\{([0-9]*)\}',r'sqrt(\1)',text) + text=text.replace(')sqrt',')*sqrt') + + # re-number variables in a sequential fashion + variables=[] + for match in re.finditer(r'\[([0-9]*)\]',text): + index=int(match.group(1)) + if index not in variables: + variables.append(index) + variables.sort() + for index in range(len(variables)): + text=text.replace('['+str(variables[index])+']','[new'+str(index)+']') + text=text.replace('[new','[') + + return(text+';') + +# convert to LaTeX format +def convert_latex(text, lsym, Lsym, Csym, oneline, columns): + # remove newlines + if (oneline==0): + text=text.replace('\n','\\\\\n') + else: + # end of lines + text=text.replace(',',',\\\\') + + # remove extra space + text=re.sub(r' {2,}',' ',text) + + # replace left hand side variables + text=re.sub(r'\[(d*)C *([0-9]*)\] =',r'@\1'+Csym+r'_{\2} =', text) + text=re.sub(r'\[(d*)% *([0-9]*)\] =',r'@\1'+Lsym+r'_{\2} =', text) + # replace right hand side variables + text=re.sub(r'\[(d*)% *([0-9]*)\]',r'@\1'+lsym+r'_{\2}',text) + text=re.sub(r'\[(d*)C *([0-9]*)\]',r'@\1'+Csym+r'_{\2}',text) + + # rewrite derivatives + # find the highest order of derivatives + i=0 + dstr="d" + index=0 + while (index>=0): + index=text.find(r'@'+dstr,index) + i=i+1 + dstr=dstr+"d" + #replace derivatives + for j in range(i): + # number of derivs + nrd=i-1-j + dstr="" + for k in range(nrd): + dstr=dstr+"d" + if (nrd>1): + text=text.replace(r'@'+dstr,"\\partial^{"+str(nrd)+"}") + elif (nrd==1): + text=text.replace(r'@'+dstr,"\\partial ") + else: + text=text.replace(r'@'+dstr,"") + + # replace constants in right hand side + text=re.sub(r'\[/C([0-9]*)\^([0-9]*)\]',r'\\frac1{'+Csym+r'_{\1}^{\2}}',text) + + # remove "^{1}" + text=text.replace("^{1}","") + + # replace square roots + text=text.replace('s{','\\sqrt{') + + # fractions + text=re.sub(r'\((-?)([0-9]*)/([0-9]*)\)',r'\1\\frac{\2}{\3}',text) + # numbers + text=re.sub(r'\(([-0-9]*)\)',r'\1',text) + + + return(text) + +# read arguments +argc=len(sys.argv) +if (argc<=1): + print_usage() + exit(-1) +i=1 +infile="" +while i +#include +#include +#include "istring.h" +#include "definitions.cpp" + +// init +int init_Int_Array(Int_Array* array, int memory){ + (*array).values=calloc(memory,sizeof(int)); + (*array).memory=memory; + (*array).length=0; + return(0); +} +int free_Int_Array(Int_Array array){ + free(array.values); + return(0); +} + +// copy +int int_array_cpy(Int_Array input, Int_Array* output){ + init_Int_Array(output,input.length); + int_array_cpy_noinit(input,output); + return(0); +} +int int_array_cpy_noinit(Int_Array input, Int_Array* output){ + int i; + if((*output).memory=(*output).memory){ + int_array_resize(output,2*(*output).memory+1); + } + (*output).values[(*output).length]=val; + (*output).length++; + return(0); +} + +// concatenate +int int_array_concat(Int_Array input, Int_Array* output){ + int i; + int offset=(*output).length; + if((*output).length+input.length>(*output).memory){ + // make it longer than needed by (*output).length (for speed) + int_array_resize(output,2*(*output).length+input.length); + } + for(i=0;iarray2.length){ + return(1); + } + + // compare terms + for(i=0;iarray2.values[i]){ + return(1); + } + } + + // if equal + return(0); +} + +// check whether an array is a sub-array of another +// allows for the elements to be separated by others, but the order must be respected +int int_array_is_subarray_ordered(Int_Array input, Int_Array test_array){ + int i; + int matches=0; + + for(i=0;i=(*output).memory){ + char_array_resize(output,2*(*output).memory+1); + } + (*output).str[(*output).length]=val; + (*output).length++; + return(0); +} +// append a string +int char_array_append_str(char* str, Char_Array* output){ + char* ptr; + for (ptr=str;*ptr!='\0';ptr++){ + char_array_append(*ptr, output); + } + return(0); +} + +// concatenate +int char_array_concat(Char_Array input, Char_Array* output){ + int i; + int offset=(*output).length; + if((*output).length+input.length>(*output).memory){ + // make it longer than needed by (*output).length (for speed) + char_array_resize(output,2*(*output).length+input.length); + } + for(i=0;isize){ + // resize + free(out_str); + // +1 for '\0' + size=extra_size+1; + out_str=calloc(size,sizeof(char)); + // read format again + va_start(vaptr, fmt); + vsnprintf(out_str,size,fmt,vaptr); + va_end(vaptr); + } + + // write to char array + for(ptr=out_str;*ptr!='\0';ptr++){ + char_array_append(*ptr, output); + } + + free(out_str); + + return(0); +} + + +// replace '%' with given character +int replace_star(char c, Char_Array str, Char_Array* out){ + int i; + init_Char_Array(out, str.length); + for(i=0;i=(*output).memory){ + str_array_resize(output,2*(*output).memory+1); + } + char_array_cpy(val, (*output).strs+(*output).length); + (*output).length++; + return(0); +} +int str_array_append_noinit(Char_Array val, Str_Array* output){ + if((*output).length>=(*output).memory){ + str_array_resize(output,2*(*output).memory+1); + } + (*output).strs[(*output).length]=val; + (*output).length++; + return(0); +} + +// concatenate +int str_array_concat(Str_Array input, Str_Array* output){ + int i; + int offset=(*output).length; + if((*output).length+input.length>(*output).memory){ + // make it longer than needed by (*output).length (for speed) + str_array_resize(output,2*(*output).length+input.length); + } + for(i=0;i(*output).memory){ + // make it longer than needed by (*output).length (for speed) + str_array_resize(output,2*(*output).length+input.length); + } + for(i=0;i=(*output).memory){ + resize_labels(output,2*(*output).memory); + } + + // copy and allocate + char_array_cpy(label,(*output).labels+offset); + (*output).indices[offset]=index; + // increment length + (*output).length++; + return(0); +} +// append an element to a labels without allocating memory +int labels_append_noinit(Char_Array label, int index, Labels* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_labels(output,2*(*output).memory); + } + + // copy without allocating + (*output).labels[offset]=label; + (*output).indices[offset]=index; + // increment length + (*output).length++; + return(0); +} + +// concatenate two labels tables +int labels_concat(Labels input, Labels* output){ + int i; + for(i=0;i +#include +#include "definitions.cpp" +#include "array.h" + +// read a config file and write the output to str_args +int read_config_file(Str_Array* str_args, const char* file, int read_from_stdin){ + FILE* infile; + char c; + Char_Array buffer; + + // allocate memory + init_Str_Array(str_args, ARG_COUNT); + init_Char_Array(&buffer, STR_SIZE); + + // read from file + if(read_from_stdin==0){ + infile=fopen(file,"r"); + if(infile==NULL){ + fprintf(stderr,"error: can't open file %s\n",file); + exit(-1); + } + } + else{ + infile=stdin; + } + + while(1){ + c=fgetc(infile); + // add to str_args + if(c=='&'){ + str_array_append_noinit(buffer,str_args); + init_Char_Array(&buffer, STR_SIZE); + } + else if(c==EOF){ + str_array_append_noinit(buffer,str_args); + break; + } + else{ + char_array_append(c,&buffer); + } + } + fclose(infile); + + return(0); +} + + +// get the title from a string argument +int get_str_arg_title(Char_Array str_arg, Char_Array* out){ + int j,k; + + init_Char_Array(out,STR_SIZE); + + // find "#!" at beginning of line + for(j=0;j=str_arg.length-2){ + fprintf(stderr, "error: an entry in the configuration file does not have a title (which should be preceeded by '#!' at the beginning of a line)\n"); + exit(-1); + } + + // get title until end of line + for(k=j+2;k +#include +#include "definitions.cpp" +#include "rational.h" +#include "istring.h" +#include "array.h" +#include "number.h" +#include "tools.h" + + +// allocate memory +int init_Coefficient(Coefficient* coef,int size){ + (*coef).factors=calloc(size,sizeof(Int_Array)); + (*coef).nums=calloc(size,sizeof(Number)); + (*coef).denoms=calloc(size,sizeof(coef_denom)); + (*coef).length=0; + (*coef).memory=size; + + return(0); +} + +// free memory +int free_Coefficient(Coefficient coef){ + int i; + for(i=0;i(*output).memory){ + resize_Coefficient(output,input.length); + } + + for(i=0;i=(*output).memory){ + resize_Coefficient(output,2*(*output).memory+1); + } + // copy and allocate + int_array_cpy(factor,(*output).factors+offset); + number_cpy(num,(*output).nums+offset); + (*output).denoms[offset]=denom; + // increment length + (*output).length++; + return(0); +} +// append an element to a coefficient without allocating memory +int coefficient_append_noinit(Int_Array factor, Number num, coef_denom denom, Coefficient* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_Coefficient(output,2*(*output).memory+1); + } + // copy without allocating + (*output).factors[offset]=factor; + (*output).nums[offset]=num; + (*output).denoms[offset]=denom; + // increment length + (*output).length++; + return(0); +} + +// concatenate coefficients and simplify result +int coefficient_concat(Coefficient input, Coefficient* output){ + int i; + for(i=0;i0){ + at_least_one=1; + coefficient_append_noinit(factor,number_Qprod_ret(quot(match_count,1),input.nums[i]), input.denoms[i],output); + } + else{ + free_Int_Array(factor); + } + } + + if(at_least_one>0){ + coefficient_simplify(output); + } + else{ + // add a trivial element to the coefficient + init_Int_Array(&factor,0); + denom.index=-1; + denom.power=0; + coefficient_append_noinit(factor,number_zero(),denom,output); + } + + return(0); +} + +// derive a monomial with respect to an index +int monomial_deriv(Int_Array factor, int index, Int_Array* out_factor, int* match_count){ + int j; + + init_Int_Array(out_factor,factor.length); + // init match count + *match_count=0; + + // loop over indices + for(j=0;jdenom2.index){ + return(-1); + } + + if(denom1.powerdenom2.power){ + return(1); + } + + return(0); +} + + +// evaluate a coefficient on a vector +int evalcoef(RCC rccs, Coefficient coef, long double* out){ + int i,j; + long double num_factor; + + *out=0; + + // for each monomial + for(i=0;i0){ + num_factor/=dpower(rccs.values[intlist_find_err(rccs.indices,rccs.length,coef.denoms[i].index)],coef.denoms[i].power); + } + *out+=num_factor*number_double_val(coef.nums[i]); + } + return(0); +} diff --git a/src/coefficient.h b/src/coefficient.h new file mode 100644 index 0000000..a7a151c --- /dev/null +++ b/src/coefficient.h @@ -0,0 +1,78 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* +Coefficients represent the collection of terms which appear + on the right hand side of a single flow equation. + +They are used as an elementary building block of grouped polynomials. +*/ + +#ifndef COEFFICIENT_H +#define COEFFICIENT_H + +#include "types.h" + +// allocate memory +int init_Coefficient(Coefficient* coef,int size); +// free memory +int free_Coefficient(Coefficient coef); + +// copy a coefficient +int coefficient_cpy(Coefficient input, Coefficient* output); +int coefficient_cpy_noinit(Coefficient input, Coefficient* output); + +// resize the memory allocated to a coefficient +int resize_Coefficient(Coefficient* coefficient,int new_size); + +// append an element to a coefficient +int coefficient_append(Int_Array factor,Number num, coef_denom denom, Coefficient* output); +int coefficient_append_noinit(Int_Array factor, Number num, coef_denom denom, Coefficient* output); + +// concatenate coefficients and simplify result +int coefficient_concat(Coefficient input, Coefficient* output); +int coefficient_concat_noinit(Coefficient input, Coefficient* output); + +// simplify a Coefficient +int coefficient_simplify(Coefficient* coefficient); +// sort the terms in an equation (quicksort algorithm) +int sort_coefficient(Coefficient* coefficient, int begin, int end); +// exchange two terms (for the sorting algorithm) +int exchange_coefficient_terms(int i, int j, Coefficient* coefficient); + +// derive a coefficient with respect to an index (as a polynomial) (does not derive the 1/(1+C)^p ) +int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output); +int coefficient_deriv(Coefficient input, int index, Coefficient* output); + +// product of two coefficients +int coefficient_prod(Coefficient coef1, Coefficient coef2, Coefficient* output); +// product of coefficients, output replaces the second coefficient +int coefficient_prod_chain(Coefficient in, Coefficient* inout); + +// print coefficient +int coefficient_sprint(Coefficient coef, Char_Array* output, int offset, char index_pre); + +// read from a string +int char_array_to_Coefficient(Char_Array str_coef, Coefficient* output); +int str_to_Coefficient(char* str, Coefficient* output); + +// compare coefficient denominators +int coef_denom_cmp(coef_denom denom1, coef_denom denom2); + +// evaluate a coefficient on a vector +int evalcoef(RCC rccs, Coefficient coef, long double* out); + +#endif diff --git a/src/definitions.cpp b/src/definitions.cpp new file mode 100644 index 0000000..d32537d --- /dev/null +++ b/src/definitions.cpp @@ -0,0 +1,60 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#ifndef DEFINITIONS_GCC +#define DEFINITIONS_GCC + +#define VERSION "1.2" + +// number of entries in a configuration file +#define ARG_COUNT 10 +// size of string representing a monomial +#define MONOMIAL_SIZE 20 +// size of various strings +#define STR_SIZE 100 +// number of terms in coefficients +#define COEF_SIZE 100 +// number of terms in polynomials +#define POLY_SIZE 100 +// number of equations +#define EQUATION_SIZE 20 +// number of fields +#define FIELDS_SIZE 50 +// number of elements in numbers +#define NUMBER_SIZE 5 +// number of elements in a group +#define GROUP_SIZE 5 + + +// display options +#define DISPLAY_EQUATION 1 +#define DISPLAY_NUMERICAL 2 +#define DISPLAY_EXPLOG 3 +#define DISPLAY_FINAL 4 + +// available preprocessors +#define PREPROCESSOR_KONDO 1 + +// offset derivative indices +#define DOFFSET 1000000 + +// types of fields (the order matters) +#define FIELD_PARAMETER 1 +#define FIELD_EXTERNAL 2 +#define FIELD_INTERNAL 3 +#define FIELD_SYMBOL 4 + +#endif diff --git a/src/expansions.c b/src/expansions.c new file mode 100644 index 0000000..c889829 --- /dev/null +++ b/src/expansions.c @@ -0,0 +1,28 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "expansions.h" + +#include +#include +#include "definitions.cpp" +#include "tools.h" +#include "array.h" +#include "polynomial.h" +#include "number.h" +#include "rational.h" + + diff --git a/src/expansions.h b/src/expansions.h new file mode 100644 index 0000000..85d6c2b --- /dev/null +++ b/src/expansions.h @@ -0,0 +1,34 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* +Compute exp(V) and log(1+W) +*/ + +#ifndef EXPANSIONS_H +#define EXPANSIONS_H + +#include "polynomial.h" +#include "fields.h" + +// exp(V) +int expand_exponential(Polynomial input_polynomial,Polynomial* output, Fields_Table fields); + +// log(1+W) +int expand_logarithm(Polynomial input_polynomial, Polynomial* output, Fields_Table fields); + + +#endif diff --git a/src/fields.c b/src/fields.c new file mode 100644 index 0000000..1b221f2 --- /dev/null +++ b/src/fields.c @@ -0,0 +1,489 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "fields.h" +#include "definitions.cpp" +#include +#include +#include "number.h" +#include "tools.h" +#include "polynomial.h" +#include "array.h" +#include "rational.h" + +// init and free for Fields_Table +int init_Fields_Table(Fields_Table* fields){ + init_Int_Array(&((*fields).parameter),FIELDS_SIZE); + init_Int_Array(&((*fields).external),FIELDS_SIZE); + init_Int_Array(&((*fields).internal),FIELDS_SIZE); + init_Identities(&((*fields).ids), FIELDS_SIZE); + init_Symbols(&((*fields).symbols), FIELDS_SIZE); + init_Int_Array(&((*fields).fermions),FIELDS_SIZE); + return(0); +} +int free_Fields_Table(Fields_Table fields){ + free_Int_Array(fields.parameter); + free_Int_Array(fields.external); + free_Int_Array(fields.internal); + free_Identities(fields.ids); + free_Symbols(fields.symbols); + free_Int_Array(fields.fermions); + return(0); +} + +// determine field type +int field_type(int index, Fields_Table fields){ + if(int_array_find(abs(index), fields.parameter)>=0){ + return(FIELD_PARAMETER); + } + else if(int_array_find(abs(index), fields.external)>=0){ + return(FIELD_EXTERNAL); + } + else if(int_array_find(abs(index), fields.internal)>=0){ + return(FIELD_INTERNAL); + } + else if(intlist_find(fields.symbols.indices, fields.symbols.length, index)>=0){ + return(FIELD_SYMBOL); + } + + fprintf(stderr,"error: index %d is neither a parameter nor an external or an internal field, nor a symbol\n",index); + exit(-1); +} + +// check whether a field anticommutes +int is_fermion(int index, Fields_Table fields){ + if(int_array_find(abs(index), fields.fermions)>=0){ + return(1); + } + else{ + return(0); + } +} + + +// ------------------ Identities -------------------- + +// allocate memory +int init_Identities(Identities* identities,int size){ + (*identities).lhs=calloc(size,sizeof(Int_Array)); + (*identities).rhs=calloc(size,sizeof(Polynomial)); + (*identities).length=0; + (*identities).memory=size; + return(0); +} + +// free memory +int free_Identities(Identities identities){ + int i; + for(i=0;i=(*output).memory){ + resize_identities(output,2*(*output).memory+1); + } + + // copy and allocate + int_array_cpy(lhs,(*output).lhs+offset); + polynomial_cpy(rhs,(*output).rhs+offset); + // increment length + (*output).length++; + return(0); +} +// append an element to a identities without allocating memory +int identities_append_noinit(Int_Array lhs, Polynomial rhs, Identities* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_identities(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).lhs[offset]=lhs; + (*output).rhs[offset]=rhs; + // increment length + (*output).length++; + return(0); +} + +// concatenate two identitiess +int identities_concat(Identities input, Identities* output){ + int i; + for(i=0;i0){ + at_least_one=0; + + // prevent infinite loops + security++; + if(security>1000000){ + fprintf(stderr,"error: 1000000 iterations reached when trying to simplify a monomial\n"); + exit(-1); + } + + // loop over ids + for(j=0;j=fields.ids.lhs[j].length || (*polynomial).monomials[i].values[k]!=fields.ids.lhs[j].values[l]){ + int_array_append((*polynomial).monomials[i].values[k],&monomial); + // sign correction + if(fermion_count % 2 ==1 && is_fermion((*polynomial).monomials[i].values[k], fields)){ + sign*=-1; + } + } + else{ + // increment fermion_count + if(is_fermion(fields.ids.lhs[j].values[l],fields)){ + fermion_count++; + } + // increment "current" field in the id + l++; + } + } + + num=number_Qprod_ret(quot(sign,1),(*polynomial).nums[i]); + // add extra monomials (if there are more than 1) + for(k=1;k=(*output).memory){ + resize_symbols(output,2*(*output).memory+1); + } + + // copy and allocate + (*output).indices[offset]=index; + polynomial_cpy(expr,(*output).expr+offset); + // increment length + (*output).length++; + return(0); +} +// append an element to a symbols without allocating memory +int symbols_append_noinit(int index, Polynomial expr, Symbols* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_symbols(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).indices[offset]=index; + (*output).expr[offset]=expr; + // increment length + (*output).length++; + return(0); +} + +// concatenate two symbolss +int symbols_concat(Symbols input, Symbols* output){ + int i; + for(i=0;i=(*output).memory){ + resize_groups(output,2*(*output).memory+1); + } + + // copy and allocate + int_array_cpy(indices,(*output).indices+offset); + // increment length + (*output).length++; + return(0); +} +// append an element to a groups without allocating memory +int groups_append_noinit(Int_Array indices, Groups* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_groups(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).indices[offset]=indices; + // increment length + (*output).length++; + return(0); +} + +// concatenate two groupss +int groups_concat(Groups input, Groups* output){ + int i; + for(i=0;i +#include +#include "tools.h" +#include "math.h" +#include "definitions.cpp" +#include "number.h" +#include "array.h" +#include "coefficient.h" + + + +// compute flow numerically, no exponentials +// inputs: flow_equation +// init, niter, tol (the allowed error at each step), ls (whether to display the results in terms of ls), display_mode (what to print) +int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, int niter, long double tol, int display_mode){ + // running coupling contants + RCC rccs=init; + int i,j; + + if(display_mode==DISPLAY_NUMERICAL){ + // print labels + printf("%5s ","n"); + for(j=0;j=0){ + evalcoef(*rccs, flow_equation.coefs[i], new_rccs+i); + // if the new rcc is too small, then ignore it + if(fabs(new_rccs[i]) +#include +#include "definitions.cpp" +#include "rational.h" +#include "istring.h" +#include "coefficient.h" +#include "polynomial.h" +#include "array.h" +#include "number.h" +#include "tools.h" + + +// allocate memory +int init_Grouped_Polynomial(Grouped_Polynomial* gpolynomial, int size){ + (*gpolynomial).coefs=calloc(size,sizeof(Coefficient)); + (*gpolynomial).indices=calloc(size,sizeof(int)); + (*gpolynomial).length=0; + (*gpolynomial).memory=size; + + return(0); +} + +// free memory +int free_Grouped_Polynomial(Grouped_Polynomial gpolynomial){ + int i; + for(i=0;i=(*output).memory){ + resize_grouped_polynomial(output,2*(*output).memory+1); + } + + // copy and allocate + (*output).indices[offset]=index; + coefficient_cpy(coef, (*output).coefs+offset); + //increment length + (*output).length++; + + return(0); +} +// append an element to a grouped_polynomial without allocating memory +int grouped_polynomial_append_noinit(int index, Coefficient coef, Grouped_Polynomial* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_grouped_polynomial(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).indices[offset]=index; + (*output).coefs[offset]=coef; + // increment length + (*output).length++; + return(0); +} + +// concatenate two grouped_polynomials +int grouped_polynomial_concat(Grouped_Polynomial input, Grouped_Polynomial* output){ + int i; + for(i=0;i0){ + // stop if the number of iterations exceeds 100 times the length of the polynomial + if(security >= 100*polynomial.length){ + fprintf(stderr,"error: polynomial could not be grouped in less than %d groupings\n", 100*polynomial.length); + exit(-1); + } + security++; + + // index of the last element + i=remainder.length-1; + + // find entry + if(remainder.monomials[i].length==0){ + // constant + index=-1; + } + else{ + // loop over entries + for(j=0,index=-2;j=0){ + ratio=number_quot_ret(remainder.nums[i],idtable.polynomials[index].nums[pos]); + factor=remainder.factors[i]; + // add to coefficient + denom.index=-1; + denom.power=1; + coefficient_append(factor, ratio, denom, (*grouped_polynomial).coefs+index+1); + + // remove from remainder + free_Int_Array(remainder.monomials[i]); + // do not free factor yet + free_Number(remainder.nums[i]); + remainder.length--; + + // add terms from idtable with minus sign + for(j=0;j=0){ + // a vector in which to store the indices that were masked + init_Int_Array(&mask_tmp_flips,idtable.polynomials[index].length); + + // loop over all monomials in that entry of the idtable + for(j=0;j=0){ + (*dflow).indices[i]=flow_equation.indices[i]+DOFFSET; + } + else{ + (*dflow).indices[i]=flow_equation.indices[i]-DOFFSET; + } + + init_Coefficient((*dflow).coefs+i, COEF_SIZE); + // for each index + for(j=0;j=0){ + int_array_append(DOFFSET + indices.values[j], tmp_coef.factors+k); + } + // constants are offset with -doffset (so that the derivatives of constants also have a negative index) + else{ + int_array_append(-DOFFSET + indices.values[j], tmp_coef.factors+k); + } + } + } + + // add to output + coefficient_concat_noinit(tmp_coef, (*dflow).coefs+i); + } + } + + (*dflow).length=flow_equation.length; + + return(0); +} + + +/* +// derive a flow equation with respect to an index +int flow_equation_deriv(Grouped_Polynomial flow_equation, int index, Grouped_Polynomial* output){ + int i,k; + // temp list of indices + Int_Array factor; + // number of times index was found + int match_count; + coef_denom denom; + // store the computation of the derivative of the constant + int previous_constant_index=0; + Coefficient dC; + Coefficient tmp_coef; + + init_Grouped_Polynomial(output, flow_equation.length); + + // loop over equations + for(k=0;k0){ + coefficient_append_noinit(factor,number_Qprod_ret(quot(match_count,1),flow_equation.coefs[k].nums[i]), flow_equation.coefs[k].denoms[i], (*output).coefs+k); + } + else{ + free_Int_Array(factor); + } + + // derivative of the denominator + if(flow_equation.coefs[k].denoms[i].power>0){ + // check whether the derivative was already computed + if(flow_equation.coefs[k].denoms[i].index!=previous_constant_index){ + // if not first, then free + if(previous_constant_index!=0){ + free_Coefficient(dC); + previous_constant_index=0; + } + init_Coefficient(&dC,COEF_SIZE); + coefficient_deriv_noinit(flow_equation.coefs[intlist_find_err(flow_equation.indices, flow_equation.length, flow_equation.coefs[k].denoms[i].index)], index, &dC); + previous_constant_index=flow_equation.coefs[k].denoms[i].index; + } + + init_Coefficient(&tmp_coef, dC.length); + coefficient_append(flow_equation.coefs[k].factors[i], number_Qprod_ret(quot(-flow_equation.coefs[k].denoms[i].power,1), flow_equation.coefs[k].nums[i]), flow_equation.coefs[k].denoms[i], &tmp_coef); + (tmp_coef.denoms[0].power)++; + + coefficient_prod_chain(dC, &tmp_coef); + + coefficient_concat_noinit(tmp_coef, (*output).coefs+k); + } + } + + // memory safe + if((*output).coefs[k].length>0){ + coefficient_simplify((*output).coefs+k); + } + else{ + // add a trivial element to the coefficient + init_Int_Array(&factor,0); + denom.index=-1; + denom.power=0; + coefficient_append_noinit(factor,number_zero(),denom,(*output).coefs+k); + } + } + + free_Coefficient(dC); + return(0); +} +*/ + + +// print a grouped polynomial +// prepend the indices on the left side with lhs_pre, and those on the right by rhs_pre +int grouped_polynomial_print(Grouped_Polynomial grouped_polynomial, char lhs_pre, char rhs_pre){ + int i,j; + Char_Array buffer; + int dcount; + + // for each equation + for(i=0;i +#include +#include "array.h" +#include "polynomial.h" + +// allocate memory +int init_Id_Table(Id_Table* idtable,int size){ + (*idtable).indices=calloc(size,sizeof(int)); + (*idtable).polynomials=calloc(size,sizeof(Polynomial)); + (*idtable).length=0; + (*idtable).memory=size; + return(0); +} + +// free memory +int free_Id_Table(Id_Table idtable){ + int i; + for(i=0;i=(*output).memory){ + resize_idtable(output,2*(*output).memory); + } + + // copy and allocate + polynomial_cpy(polynomial,(*output).polynomials+offset); + (*output).indices[offset]=index; + // increment length + (*output).length++; + return(0); +} +// append an element to a idtable without allocating memory +int idtable_append_noinit(int index, Polynomial polynomial, Id_Table* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_idtable(output,2*(*output).memory); + } + + // copy without allocating + (*output).indices[offset]=index; + (*output).polynomials[offset]=polynomial; + // increment length + (*output).length++; + return(0); +} + +// concatenate two idtables +int idtable_concat(Id_Table input, Id_Table* output){ + int i; + for(i=0;i +#include "istring.h" + +char* str_concat(char* ptr, const char* str){ + char* str_ptr=(char*)str; + + while(*str_ptr!='\0'){ + *ptr=*str_ptr; + ptr++; + str_ptr++; + } + *ptr='\0'; + return(ptr); +} + +int str_concat_memorysafe(char** str_out, int pos, const char* str, int* memory){ + char* out_ptr; + + if(str_len((char*)str)+pos>=*memory){ + *memory=*memory+str_len((char*)str)+pos+1; + resize_str(str_out,*memory); + } + out_ptr=*str_out+pos; + str_concat(out_ptr,str); + return(0); +} + +int resize_str(char** out, int memory){ + char* tmp_str=calloc(memory,sizeof(char)); + char* tmp_ptr=tmp_str; + char* out_ptr=*out; + + for(;*out_ptr!='\0';out_ptr++,tmp_ptr++){ + *tmp_ptr=*out_ptr; + } + *tmp_ptr='\0'; + + free(*out); + *out=tmp_str; + return(0); +} + +char* str_addchar(char* ptr, const char c){ + *ptr=c; + ptr++; + *ptr='\0'; + return(ptr); +} + +int str_len(char* str){ + char* ptr=str; + int ret=0; + while(*ptr!='\0'){ret++;ptr++;} + return(ret); +} + +int str_cmp(char* str1, char* str2){ + char* ptr1=str1; + char* ptr2=str2; + while(*ptr1==*ptr2 && *ptr1!='\0' && *ptr2!='\0'){ + ptr1++; + ptr2++; + } + if(*ptr1=='\0' && *ptr2=='\0'){ + return(1); + } + else{ + return(0); + } +} + + +int str_cpy_noalloc(char* input, char* output){ + char* ptr; + char* out_ptr; + + for(ptr=input,out_ptr=output;*ptr!='\0';ptr++,out_ptr++){ + *out_ptr=*ptr; + } + *out_ptr='\0'; + return(0); +} + diff --git a/src/istring.h b/src/istring.h new file mode 100644 index 0000000..5b47b9b --- /dev/null +++ b/src/istring.h @@ -0,0 +1,40 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* +String manipulation +*/ + +#ifndef ISTRING_H +#define ISTRING_H + +// concatenate str at the position pointed to by ptr, and return a pointer +// to the end of the new string +char* str_concat(char* ptr, const char* str); +// concatenate strings and resize them if necessary +int str_concat_memorysafe(char** str_out, int pos, const char* str, int* memory); +// resize a string +int resize_str(char** out, int memory); +// idem with a single character +char* str_addchar(char* ptr, const char c); +// string length +int str_len(char* str); +// compare strings +int str_cmp(char* str1, char* str2); +// copy a string to another without allocating memory +int str_cpy_noalloc(char* input, char* output); + +#endif diff --git a/src/kondo.c b/src/kondo.c new file mode 100644 index 0000000..c86b4b3 --- /dev/null +++ b/src/kondo.c @@ -0,0 +1,1449 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "kondo.h" +#include +#include + +#include "idtable.h" +#include "array.h" +#include "number.h" +#include "istring.h" +#include "cli_parser.h" +#include "fields.h" +#include "parse_file.h" +#include "polynomial.h" +#include "definitions.cpp" +#include "rational.h" + +// dimension of A, B and h (must be <10) +#define KONDO_DIM 3 +// number of spin components +#define KONDO_SPIN 2 + +// offsets for indices of A, B and h +// order matters for symbols table +#define KONDO_A_OFFSET 1 +#define KONDO_B_OFFSET 2 +#define KONDO_H_OFFSET 3 + +// parsing modes (from parse_file.c) +#define PP_NULL_MODE 0 +// when reading a factor +#define PP_FACTOR_MODE 1 +// reading a monomial +#define PP_MONOMIAL_MODE 2 +// reading a numerator and denominator +#define PP_NUMBER_MODE 3 +// types of fields +#define PP_FIELD_MODE 6 +#define PP_PARAMETER_MODE 7 +#define PP_EXTERNAL_MODE 8 +#define PP_INTERNAL_MODE 9 +// indices +#define PP_INDEX_MODE 10 +// factors or monomials +#define PP_BRACKET_MODE 11 +// labels +#define PP_LABEL_MODE 12 +// polynomial +#define PP_POLYNOMIAL_MODE 13 +// field scalar product +#define PP_FIELD_SCALAR_MODE 14 +#define PP_FIELD_VECTOR_PROD_MODE 15 + + + +// generate configuration file +int kondo_generate_conf(Str_Array* str_args, int box_count){ + Str_Array new_args; + Fields_Table fields; + Char_Array tmp_str; + int arg_index; + int i; + Char_Array title; + + init_Str_Array(&new_args,8); + + // fields + kondo_fields_table(box_count, &tmp_str, &fields); + str_array_append_noinit(tmp_str, &new_args); + + // symbols + kondo_symbols(&tmp_str, box_count, &fields); + arg_index=find_str_arg("symbols", *str_args); + if(arg_index>=0){ + if(tmp_str.length>0){ + char_array_snprintf(&tmp_str,",\n"); + } + char_array_concat((*str_args).strs[arg_index], &tmp_str); + } + parse_input_symbols(tmp_str, &fields); + str_array_append_noinit(tmp_str, &new_args); + + // identities + kondo_identities(&tmp_str); + arg_index=find_str_arg("identities", *str_args); + if(arg_index>=0){ + if(tmp_str.length>0){ + char_array_snprintf(&tmp_str,",\n"); + } + char_array_concat((*str_args).strs[arg_index], &tmp_str); + } + parse_input_identities(tmp_str, &fields); + str_array_append_noinit(tmp_str, &new_args); + + // groups + kondo_groups(&tmp_str, box_count); + str_array_append_noinit(tmp_str, &new_args); + + + // propagator + arg_index=find_str_arg("propagator", *str_args); + if(arg_index>=0){ + kondo_propagator((*str_args).strs[arg_index], &tmp_str); + str_array_append_noinit(tmp_str, &new_args); + } + + // input polynomial + arg_index=find_str_arg("input_polynomial", *str_args); + if(arg_index>=0){ + kondo_input_polynomial((*str_args).strs[arg_index], &tmp_str, fields, box_count); + str_array_append_noinit(tmp_str, &new_args); + } + + // id table + arg_index=find_str_arg("id_table", *str_args); + if(arg_index>=0){ + kondo_idtable((*str_args).strs[arg_index], &tmp_str, fields); + str_array_append_noinit(tmp_str, &new_args); + } + + // copy remaining entries + for(i=0;i<(*str_args).length;i++){ + get_str_arg_title((*str_args).strs[i], &title); + if(str_cmp(title.str, "symbols")==0 &&\ + str_cmp(title.str, "identities")==0 &&\ + str_cmp(title.str, "propagator")==0 &&\ + str_cmp(title.str, "input_polynomial")==0 &&\ + str_cmp(title.str, "id_table")==0 ){ + + char_array_cpy((*str_args).strs[i], &tmp_str); + str_array_append_noinit(tmp_str, &new_args); + } + free_Char_Array(title); + } + + free_Fields_Table(fields); + free_Str_Array(*str_args); + *str_args=new_args; + + return(0); +} + + +// generate the Kondo fields table +int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fields){ + int i,j; + + init_Char_Array(str_fields,STR_SIZE); + char_array_snprintf(str_fields, "#!fields\n"); + + // external fields + char_array_append_str("x:",str_fields); + for(i=0;i=0 && offset[1]>=0 && index[0]>=0 && index[1]>=0){ + // write indices and num + for(i=0;i0){ + // write num + polynomial_multiply_scalar(tmp_poly, tmp_num); + // replace factor + for(i=0;i0){ + for(i=0;i0){ + for(i=0;i': + if(mode==PP_MONOMIAL_MODE){ + // resolve scalar product + kondo_resolve_scalar_prod(buffer, &scalar_prod_poly, fields); + // add to tmp_poly + if(tmp_poly.length==0){ + polynomial_concat(scalar_prod_poly,&tmp_poly); + } + else{ + polynomial_prod_chain(scalar_prod_poly,&tmp_poly,fields); + } + free_Polynomial(scalar_prod_poly); + + mode=PP_NULL_MODE; + } + break; + + // characters to ignore + case ' ':break; + case '&':break; + case '\n':break; + + // comments + case '#': + comment=1; + break; + + default: + if(mode!=PP_NULL_MODE){ + // write to buffer + buffer_ptr=str_addchar(buffer_ptr,*polynomial_ptr); + } + break; + } + } + } + + // last term + if(tmp_poly.length>0){ + polynomial_multiply_scalar(tmp_poly,tmp_num); + for(i=0;i=0){ + index=*ptr-'0'; + } + else{ + dim=*ptr-'0'; + } + } + } + + // turn B3 into B2 and B4 into B1 + if(offset==KONDO_B_OFFSET){ + switch(index){ + case 3: + index=2; + break; + case 4: + index=1; + break; + } + } + + // h's + if(offset==KONDO_H_OFFSET){ + // external field + init_Int_Array(&monomial,1); + init_Int_Array(&factor,1); + int_array_append(10*(dim+10*offset), &monomial); + polynomial_append_noinit(monomial, factor, number_one(), output); + } + // psi's + else{ + // construct spin indices + for(i=0;i0){ + init_Int_Array(&monomial,1); + init_Int_Array(&factor,1); + + int_array_append(10*(i+10*offset)+index, &monomial); + polynomial_append_noinit(monomial, factor, number_one(), psi+i); + } + } + + // multiply by Pauli matrices + Pauli_matrix(dim+1,&pauli_mat); + for(a=0;a=0){ + kondo_polynomial_vector(offset, index, &poly_vect1, fields); + operation=K_VECT_PROD; + } + break; + + // index + default: + // char to int + index=*ptr-'0'; + } + } + + // final scalar product + if(operation==K_SCALAR_PROD){ + if(offset>=0){ + kondo_polynomial_vector(offset, index, &poly_vect2, fields); + kondo_polynomial_scalar_product(poly_vect1, poly_vect2, output, fields); + } + } + + // free memory + for(i=0;i0){ + init_Int_Array(&monomial,1); + init_Int_Array(&factor,1); + + int_array_append(10*(i+10*offset)+index, &monomial); + polynomial_append_noinit(monomial, factor, number_one(), psi+i); + } + } + + // multiply by Pauli matrices + for(i=0;i=10 || *ptr-'0'<0) && (*ptr!='-')){ + break; + } + } + if(*ptr=='\0'){ + sscanf(str,"%d",&index); + return(index); + } + + for(ptr=str;*ptr!='\0';ptr++){ + switch(*ptr){ + case 'A': + offset=KONDO_A_OFFSET; + break; + case 'a': + offset=KONDO_A_OFFSET; + break; + case 'B': + offset=KONDO_B_OFFSET; + break; + case 'b': + offset=KONDO_B_OFFSET; + break; + case 'h': + offset=KONDO_H_OFFSET; + break; + default: + // set index if dim was already set + if(dim>=0){ + index=*ptr-'0'; + } + else{ + dim=*ptr-'0'; + } + } + } + + if(offset==-1){ + return(-1); + } + // no symbol for h + if(offset==KONDO_H_OFFSET){ + return(10*(10*offset+dim)); + } + else{ + return(100*(10*offset+dim)+index); + } +} diff --git a/src/kondo.h b/src/kondo.h new file mode 100644 index 0000000..23756ad --- /dev/null +++ b/src/kondo.h @@ -0,0 +1,76 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* Generate configuration files specific to the Kondo model */ + +#ifndef KONDO_H +#define KONDO_H + +#include "types.h" + +// generate configuration file +int kondo_generate_conf(Str_Array* str_args, int box_count); + +// generate the Kondo fields table +int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fields); + +// generate Kondo symbols +int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields); +// generate Kondo symbols (older method: one symbol for each scalar product) +int kondo_symbols_scalarprod(Char_Array* str_symbols, int box_count, Fields_Table* fields); + +// generate Kondo groups (groups of independent variables) +int kondo_groups(Char_Array* str_groups, int box_count); + +// generate Kondo identities +int kondo_identities(Char_Array* str_identities); + +// convert the Kondo propagator +int kondo_propagator(Char_Array str_kondo_propagator, Char_Array* str_propagator); + +// read a product of polynomials +int parse_kondo_polynomial_factors(Char_Array str_polynomial, Polynomial* output, Fields_Table fields); + +// convert Kondo input polynomial +int kondo_input_polynomial(Char_Array str_kondo_polynomial, Char_Array* str_polynomial, Fields_Table fields, int box_count); + +// convert the Kondo idtable +int kondo_idtable(Char_Array str_kondo_idtable, Char_Array* str_idtable, Fields_Table fields); + +// read a kondo polynomial and convert it to a polynomial expressed in terms of the fields in the fields table +int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_Table fields); +int parse_kondo_polynomial(Char_Array kondo_polynomial_str, Polynomial* polynomial, Fields_Table fields); + +// read Aij, Bij, hi where i is a space dimension and j is a box index +int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields); +// read a Kondo scalar product +int kondo_resolve_scalar_prod(char* str, Polynomial* output, Fields_Table fields); +// compute a scalar product of polynomial vectors +int kondo_polynomial_scalar_product(Polynomial poly_vect1[3], Polynomial poly_vect2[3], Polynomial* output, Fields_Table fields); +// compute a vector product of polynomial vectors +int kondo_polynomial_vector_product(Polynomial (*poly_vect1)[3], Polynomial poly_vect2[3], Fields_Table fields); +// compute the 3 components of a kondo vector +int kondo_polynomial_vector(int offset, int index, Polynomial (*polys)[3], Fields_Table fields); +// read a scalar product of symbols +int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output); + +// get the offset and index of a monomial term +int get_offset_index(char* str, int* offset, int* index); +// get the offsets and index of a scalar product +int get_offsets_index(char* str, int* offset1, int* offset2, int* index); +// get the index of the symbol corresponding to a given string +int get_symbol_index(char* str); +#endif diff --git a/src/kondo_preprocess.c b/src/kondo_preprocess.c new file mode 100644 index 0000000..3371014 --- /dev/null +++ b/src/kondo_preprocess.c @@ -0,0 +1,126 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* +kondo_preprocess + +Generate configuration files for the Kondo problem + +*/ + + +#include +#include + +#include "definitions.cpp" +#include "kondo.h" +#include "cli_parser.h" +#include "array.h" + + +// read cli arguments +int read_args_kondo_pp(int argc,const char* argv[], Str_Array* str_args, Kondopp_Options* opts); +// print usage message +int print_usage_kondo_pp(); + +int main (int argc, const char* argv[]){ + int i; + // string arguments + Str_Array str_args; + // options + Kondopp_Options opts; + + // read command-line arguments + read_args_kondo_pp(argc,argv,&str_args,&opts); + + kondo_generate_conf(&str_args, 2*opts.dimension); + for(i=0;i\n\n"); + return(0); +} + diff --git a/src/mean.c b/src/mean.c new file mode 100644 index 0000000..aad7a5a --- /dev/null +++ b/src/mean.c @@ -0,0 +1,793 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* +As of version 1.0, the mean of a monomial is computed directly +*/ + +#include "mean.h" + +#include +#include +#include +#include "definitions.cpp" +#include "tools.h" +#include "polynomial.h" +#include "rational.h" +#include "array.h" +#include "fields.h" +#include "number.h" + +// mean of a monomial +int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Matrix propagator){ + int sign=1; + // +/- internal fields + Int_Array internal_plus; + Int_Array internal_minus; + + // init out + *out=polynomial_one(); + + // sort first + monomial_sort(monomial, 0, monomial.length-1, fields, &sign); + polynomial_multiply_Qscalar(*out, quot(sign,1)); + // get internals + // (*out).monomials is the first element of out but it only has 1 element + // first, free (*out).monomials[0] + free_Int_Array((*out).monomials[0]); + get_internals(monomial, &internal_plus, &internal_minus, (*out).monomials, fields); + + if(internal_plus.length>0 && internal_minus.length>0){ + mean_internal(internal_plus, internal_minus, out, propagator, fields); + } + + free_Int_Array(internal_plus); + free_Int_Array(internal_minus); + return(0); +} + +// compute the mean of a monomial of internal fields (with split + and -) +int mean_internal(Int_Array internal_plus, Int_Array internal_minus, Polynomial* out, Polynomial_Matrix propagator, Fields_Table fields){ + if(internal_plus.length!=internal_minus.length){ + fprintf(stderr,"error: monomial contains unmatched +/- fields\n"); + exit(-1); + } + int n=internal_minus.length; + // pairing as an array of positions + int* pairing=calloc(n,sizeof(int)); + // specifies which indices are available for pairing + int* mask=calloc(n,sizeof(int)); + // signature of the permutation + int pairing_sign; + // sign from mixing - and + together + int mixing_sign; + Polynomial num; + Polynomial num_summed=polynomial_zero(); + // propagator matrix indices + int* indices_minus=calloc(n,sizeof(int)); + int* indices_plus=calloc(n,sizeof(int)); + int i; + // whether the next pairing exists + int exists_next=0; + + // indices + for(i=0;i=0 && move=0 && move=0 && start=0;i--){ + // move i-th + mask[pairing[i]]=0; + // search for next possible position + for(j=pairing[i]+1;j=n){ + break; + } + } + } + } + // if the next position was found, then don't try to move the previous pairings + break; + } + // if no next position is found, store the pairing outside the array (so the algorithm stops when the first pairing is outside the array) + else{ + pairing[i]=n; + } + } + } + + number_prod_chain(num_summed,outnum); + free_Number(num_summed); + free(pairing); + free(mask); + return(0); +} +*/ + + +// get lists of internal fields from a monomial (separate + and -) +// requires the monomial to be sorted (for the sign to be correct) +int get_internals(Int_Array monomial, Int_Array* internal_plus, Int_Array* internal_minus, Int_Array* others, Fields_Table fields){ + int i; + init_Int_Array(internal_plus, monomial.length); + init_Int_Array(internal_minus, monomial.length); + init_Int_Array(others, monomial.length); + for (i=0;i=0){ + // split +/- fields + if(monomial.values[i]>0){ + int_array_append(monomial.values[i],internal_plus); + } + else{ + int_array_append(monomial.values[i],internal_minus); + } + } + else{ + int_array_append(monomial.values[i], others); + } + } + return(0); +} + + +// compute the mean of a monomial containing symbolic expressions +// keep track of which means were already computed +int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed){ + Int_Array symbol_list; + int i; + int power; + int* current_term; + Polynomial mean_num; + Int_Array tmp_monomial; + Number tmp_num; + Int_Array base_monomial; + int sign; + // whether or not the next term exists + int exists_next=0; + // simplify polynomial periodically + int simplify_freq=1; + Polynomial mean_poly; + + init_Polynomial(output, POLY_SIZE); + + // check whether the mean was already computed + for(i=0;i<(*computed).length;i++){ + if(int_array_cmp((*computed).lhs[i], monomial)==0){ + // write polynomial + polynomial_concat((*computed).rhs[i], output); + return(0); + } + } + + init_Int_Array(&symbol_list, monomial.length); + init_Int_Array(&base_monomial, monomial.length); + + // generate symbols list + for(i=0;i=0 && move=0 && move=0 && start0){ + sign=1; + monomial_sort_groups(monomial, 0, monomial.length-1, fields, groups, &sign); + } + // construct groups and take mean + init_Int_Array(&tmp_monomial, MONOMIAL_SIZE); + for(i=0;i<=monomial.length;i++){ + // new group + if(i0 && next_group!=group) || i==monomial.length){ + mean_symbols(tmp_monomial, &tmp_poly, fields, propagator, groups, computed); + // if zero + if(polynomial_is_zero(tmp_poly)==1){ + // set output to 0 + free_Polynomial(*output); + init_Polynomial(output, 1); + free_Polynomial(tmp_poly); + break; + } + // add to output + if(polynomial_is_zero(*output)==1){ + polynomial_concat(tmp_poly, output); + } + else{ + polynomial_prod_chain_nosimplify(tmp_poly, output, fields); + } + free_Polynomial(tmp_poly); + + // reset tmp_monomial + free_Int_Array(tmp_monomial); + init_Int_Array(&tmp_monomial, MONOMIAL_SIZE); + } + + // add to monomial + if(i +#include + +// pre-compiler definitions +#include "definitions.cpp" + +// various arrays +#include "array.h" + +// list of fields +#include "fields.h" +// numbers +#include "number.h" +// polynomials +#include "polynomial.h" +// list of rccs +#include "idtable.h" +// grouped representation of polynomials +#include "grouped_polynomial.h" +// command line parser +#include "cli_parser.h" +// parse input file +#include "parse_file.h" +// means +#include "mean.h" +// various string operations +#include "istring.h" + +// read cli arguments +int read_args_meankondo(int argc,const char* argv[], Str_Array* str_args, Meankondo_Options* opts); +// print usage message +int print_usage_meankondo(); +// compute flow +int compute_flow(Str_Array str_args, Meankondo_Options opts); +// compute the flow equation +int compute_flow_equation(Polynomial init_poly, Id_Table idtable, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, Grouped_Polynomial* flow_equation); + + +int main (int argc, const char* argv[]){ + // string arguments + Str_Array str_args; + // options + Meankondo_Options opts; + + // read command-line arguments + read_args_meankondo(argc,argv,&str_args,&opts); + + // warning message if representing rational numbers as floats +#ifdef RATIONAL_AS_FLOAT + fprintf(stderr,"info: representing rational numbers using floats\n"); +#endif + + compute_flow(str_args, opts); + + //free memory + free_Str_Array(str_args); + return(0); +} + + +// parse command-line arguments +#define CP_FLAG_THREADS 1 +int read_args_meankondo(int argc,const char* argv[], Str_Array* str_args, Meankondo_Options* opts){ + int i; + // pointers + char* ptr; + // file to read the polynomial from in flow mode + const char* file=""; + // flag that indicates what argument is being read + int flag=0; + // whether a file was specified on the command-line + int exists_file=0; + + + // defaults + // single thread + (*opts).threads=1; + // do not chain + (*opts).chain=0; + + // loop over arguments + for(i=1;i\n\n"); + return(0); +} + + +// compute the renormalization group flow +int compute_flow(Str_Array str_args, Meankondo_Options opts){ + int i; + // index of the entry in the input file + int arg_index; + // header of the entry + Char_Array arg_header; + // list of fields + Fields_Table fields; + // their propagator + Polynomial_Matrix propagator; + // initial polynomial + Polynomial init_poly; + // list of rccs + Id_Table idtable; + // groups of independent fields + Groups groups; + // flow equation + Grouped_Polynomial flow_equation; + + + // parse fields + arg_index=find_str_arg("fields", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no fields entry in the configuration file\n"); + exit(-1); + } + else{ + parse_input_fields(str_args.strs[arg_index],&fields); + } + + // parse id table + arg_index=find_str_arg("id_table", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no id table entry in the configuration file\n"); + exit(-1); + } + else{ + parse_input_id_table(str_args.strs[arg_index],&idtable, fields); + } + + // parse symbols + arg_index=find_str_arg("symbols", str_args); + if(arg_index>=0){ + parse_input_symbols(str_args.strs[arg_index],&fields); + } + else{ + init_Symbols(&(fields.symbols),1); + } + + // parse input polynomial + arg_index=find_str_arg("input_polynomial", str_args); + if(arg_index>=0){ + parse_input_polynomial(str_args.strs[arg_index],&init_poly, fields); + } + else{ + fprintf(stderr,"error: no input polynomial entry in the configuration file\n"); + exit(-1); + } + + // propagator + arg_index=find_str_arg("propagator", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no propagator entry in the configuration file\n"); + exit(-1); + } + else{ + parse_input_propagator(str_args.strs[arg_index],&propagator, fields); + } + + // parse identities + arg_index=find_str_arg("identities", str_args); + if(arg_index>=0){ + parse_input_identities(str_args.strs[arg_index],&fields); + } + else{ + init_Identities(&(fields.ids),1); + } + + // parse groups + arg_index=find_str_arg("groups", str_args); + if(arg_index>=0){ + parse_input_groups(str_args.strs[arg_index],&groups); + } + else{ + init_Groups(&groups, 1); + } + + // flow equation + compute_flow_equation(init_poly, idtable, fields, propagator, groups, opts.threads, &flow_equation); + free_Polynomial(init_poly); + free_Polynomial_Matrix(propagator); + free_Fields_Table(fields); + free_Groups(groups); + + // if chain then print config file + if(opts.chain==1){ + for(i=0;i1){ + polynomial_mean_multithread(&exp_poly, fields, propagator, groups, threads); + } + else{ + polynomial_mean(&exp_poly, fields, propagator, groups); + } + + // grouped representation of expanded_poly + group_polynomial(exp_poly,flow_equation,idtable, fields); + free_Polynomial(exp_poly); + + return(0); +} diff --git a/src/meantools.c b/src/meantools.c new file mode 100644 index 0000000..f45c376 --- /dev/null +++ b/src/meantools.c @@ -0,0 +1,116 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* +meantools + +Utility to perform various operations on flow equations + +*/ + + +#include +#include + +// pre-compiler definitions +#include "definitions.cpp" + +// grouped representation of polynomials +#include "grouped_polynomial.h" +// command line parser +#include "cli_parser.h" +// parse input file +#include "parse_file.h" +// arrays +#include "array.h" +// string functions +#include "istring.h" +// tools +#include "meantools_exp.h" +#include "meantools_deriv.h" +#include "meantools_eval.h" + +#define EXP_COMMAND 1 +#define DERIV_COMMAND 2 +#define EVAL_COMMAND 3 + +// read cli arguments +int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meantools_Options* opts); +// print usage message +int print_usage_meantools(); + + +int main (int argc, const char* argv[]){ + // string arguments + Str_Array str_args; + // options + Meantools_Options opts; + + // read command-line arguments + read_args_meantools(argc,argv,&str_args, &opts); + + switch(opts.command){ + case EXP_COMMAND: tool_exp(str_args); + break; + case DERIV_COMMAND: tool_deriv(str_args,opts); + break; + case EVAL_COMMAND: tool_eval(str_args,opts); + break; + } + + //free memory + free_Str_Array(str_args); + return(0); +} + + +// parse command-line arguments +int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meantools_Options* opts){ + + // if there are no arguments + if(argc==1){ + print_usage_meantools(); + exit(-1); + } + + if(str_cmp((char*)argv[1],"exp")==1){ + (*opts).command=EXP_COMMAND; + tool_exp_read_args(argc, argv, str_args); + } + else if(str_cmp((char*)argv[1],"derive")==1){ + (*opts).command=DERIV_COMMAND; + tool_deriv_read_args(argc, argv, str_args, opts); + } + else if(str_cmp((char*)argv[1],"eval")==1){ + (*opts).command=EVAL_COMMAND; + tool_eval_read_args(argc, argv, str_args, opts); + } + else{ + print_usage_meantools(); + exit(-1); + } + + return(0); +} + +// print usage message +int print_usage_meantools(){ + printf("\nusage:\n meantools exp \n meantools derive [-d derivatives] -V \n meantools eval -R \n\n"); + return(0); +} + + + diff --git a/src/meantools_deriv.c b/src/meantools_deriv.c new file mode 100644 index 0000000..28d8641 --- /dev/null +++ b/src/meantools_deriv.c @@ -0,0 +1,195 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "meantools_deriv.h" + +#include +#include +#include "parse_file.h" +#include "cli_parser.h" +#include "istring.h" +#include "definitions.cpp" +#include "array.h" +#include "grouped_polynomial.h" + + +#define CP_FLAG_DERIVS 1 +#define CP_FLAG_VARS 2 +// read command line arguments +int tool_deriv_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts){ + // file to read the polynomial from in flow mode + const char* file=""; + // whether a file was specified on the command-line + int exists_file=0; + // flag + int flag=0; + // buffer in which to read the variables + Char_Array buffer; + int i; + char* ptr; + + // defaults + // derive once + (*opts).deriv_derivs=1; + // derive with respect to all variables + (*opts).deriv_vars.length=-1; + + + // loop over arguments + for(i=2;i=0){ + int_array_read(str_args.strs[arg_index],&(opts.deriv_vars)); + } + } + } + + // if variables length is negative then set the variables to all of the available ones + if(opts.deriv_vars.length<0){ + init_Int_Array(&(opts.deriv_vars), flow_equation.length); + for(i=0;i=0){ + int_array_append((j+1)*DOFFSET+variables.values[i], &indices); + } + // constants have a negative index + else{ + int_array_append(-(j+1)*DOFFSET+variables.values[i], &indices); + } + } + + // add to flow equation + grouped_polynomial_concat(dflow, flow_equation_derivs); + } + + if(n>0){ + free_Grouped_Polynomial(dflow); + } + free_Int_Array(indices); + return(0); +} diff --git a/src/meantools_deriv.h b/src/meantools_deriv.h new file mode 100644 index 0000000..4ea36a9 --- /dev/null +++ b/src/meantools_deriv.h @@ -0,0 +1,30 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#ifndef MEANTOOLS_DERIV_H +#define MEANTOOLS_DERIV_H + +#include "types.h" + +// read arguments +int tool_deriv_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts); +// derive a flow equation +int tool_deriv(Str_Array str_args, Meantools_Options opts); +// n first derivatives of a flow equation wrt to variables +int flow_equation_derivative(int n, Int_Array variables, Grouped_Polynomial flow_equation, Grouped_Polynomial* flow_equation_derivs); + +#endif + diff --git a/src/meantools_eval.c b/src/meantools_eval.c new file mode 100644 index 0000000..50fff5b --- /dev/null +++ b/src/meantools_eval.c @@ -0,0 +1,129 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "meantools_eval.h" + +#include +#include +#include "parse_file.h" +#include "cli_parser.h" +#include "grouped_polynomial.h" +#include "array.h" +#include "rcc.h" + + +#define CP_FLAG_RCCS 1 +// read command line arguments +int tool_eval_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts){ + // file to read the polynomial from in flow mode + const char* file=""; + // whether a file was specified on the command-line + int exists_file=0; + // flag + int flag=0; + int i; + char* ptr; + + // defaults + // mark rccstring so that it can be recognized whether it has been set or not + (*opts).eval_rccstring.length=-1; + + // loop over arguments + for(i=2;i=0){ + char_array_cpy(str_args.strs[arg_index],&(opts.eval_rccstring)); + } + } + } + + // initialize the rccs + prepare_init(flow_equation.indices,flow_equation.length,&rccs); + // read rccs from string + if(opts.eval_rccstring.length!=-1){ + parse_init_cd(opts.eval_rccstring, &rccs); + free_Char_Array(opts.eval_rccstring); + } + + // evaluate + evaleq(&rccs, flow_equation); + + // print + RCC_print(rccs); + + // free memory + free_Grouped_Polynomial(flow_equation); + free_RCC(rccs); + return(0); +} + diff --git a/src/meantools_eval.h b/src/meantools_eval.h new file mode 100644 index 0000000..518049d --- /dev/null +++ b/src/meantools_eval.h @@ -0,0 +1,29 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#ifndef MEANTOOLS_EVAL_H +#define MEANTOOLS_EVAL_H + +#include "types.h" + +// read arguments +int tool_eval_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts); +// evaluate a flow equation +int tool_eval(Str_Array str_args, Meantools_Options opts); + +#endif + + diff --git a/src/meantools_exp.c b/src/meantools_exp.c new file mode 100644 index 0000000..1aca928 --- /dev/null +++ b/src/meantools_exp.c @@ -0,0 +1,130 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "meantools_exp.h" + +#include +#include +#include "parse_file.h" +#include "cli_parser.h" +#include "polynomial.h" +#include "fields.h" +#include "grouped_polynomial.h" +#include "idtable.h" + +// read command line arguments +int tool_exp_read_args(int argc, const char* argv[], Str_Array* str_args){ + // file to read the polynomial from in flow mode + const char* file=""; + // whether a file was specified on the command-line + int exists_file=0; + + if(argc>=3){ + file=argv[2]; + exists_file=1; + } + read_config_file(str_args, file, 1-exists_file); + + return(0); +} + + +// compute the exponential of the input polynomial +int tool_exp(Str_Array str_args){ + // index of the entry in the input file + int arg_index; + // list of fields + Fields_Table fields; + // input polynomial + Polynomial poly; + // exp as a polynomial + Polynomial exp_poly; + // list of rccs + Id_Table idtable; + // exp + Grouped_Polynomial exp; + int i,j; + + // parse fields + arg_index=find_str_arg("fields", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no fields entry in the configuration file\n"); + exit(-1); + } + else{ + parse_input_fields(str_args.strs[arg_index],&fields); + } + + // parse id table + arg_index=find_str_arg("id_table", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no id table entry in the configuration file\n"); + exit(-1); + } + else{ + parse_input_id_table(str_args.strs[arg_index],&idtable, fields); + } + + // parse input polynomial + arg_index=find_str_arg("input_polynomial", str_args); + if(arg_index>=0){ + parse_input_polynomial(str_args.strs[arg_index],&poly, fields); + } + else{ + fprintf(stderr,"error: no input polynomial entry in the configuration file\n"); + exit(-1); + } + + // parse symbols + arg_index=find_str_arg("symbols", str_args); + if(arg_index>=0){ + parse_input_symbols(str_args.strs[arg_index],&fields); + } + else{ + init_Symbols(&(fields.symbols),1); + } + + // parse identities + arg_index=find_str_arg("identities", str_args); + if(arg_index>=0){ + parse_input_identities(str_args.strs[arg_index],&fields); + } + else{ + init_Identities(&(fields.ids),1); + } + + // exp(V) + polynomial_exponential(poly,&exp_poly, fields); + // grouped representation + group_polynomial(exp_poly, &exp, idtable, fields); + free_Polynomial(exp_poly); + free_Polynomial(poly); + + // no denominators + for(i=0;i +#include +#include +#include "istring.h" +#include "definitions.cpp" +#include "tools.h" +#include "rational.h" +#include "array.h" + +// init +int init_Number(Number* number, int memory){ + (*number).scalars=calloc(memory,sizeof(Q)); + (*number).base=calloc(memory,sizeof(int)); + (*number).memory=memory; + (*number).length=0; + return(0); +} +int free_Number(Number number){ + free(number.scalars); + free(number.base); + return(0); +} + +// copy +int number_cpy(Number input, Number* output){ + init_Number(output,input.length); + number_cpy_noinit(input,output); + return(0); +} +int number_cpy_noinit(Number input, Number* output){ + int i; + if((*output).memory=(*output).memory){ + number_resize(output,2*(*output).memory); + } + (*output).scalars[(*output).length]=scalar; + (*output).base[(*output).length]=base; + (*output).length++; + // not optimal + number_sort(*output,0,(*output).length-1); + return(0); +} + +// concatenate +int number_concat(Number input, Number* output){ + int i; + int offset=(*output).length; + if((*output).length+input.length>(*output).memory){ + // make it longer than needed by (*output).length (for speed) + number_resize(output,2*(*output).length+input.length); + } + for(i=0;i=0){ + (*inout).scalars[index]=Q_add((*inout).scalars[index], scalar); + } + else{ + number_append(scalar, base, inout); + } + return(0); +} +// create a new number +int number_add(Number x1, Number x2, Number* out){ + number_cpy(x1,out); + number_add_chain(x2,out); + return(0); +} +// return the number +Number number_add_ret(Number x1, Number x2){ + Number out; + number_add(x1,x2,&out); + return(out); +} + +// multiply +int number_prod(Number x1, Number x2, Number* out){ + int i,j; + int div; + Q new_scalar; + int new_base; + init_Number(out, x1.length); + for(i=0;i0){ + (*inout).scalars[i]=Q_inverse((*inout).scalars[i]); + (*inout).scalars[i].denominator*=(*inout).base[i]; + } + else if((*inout).base[i]<0){ + (*inout).scalars[i]=Q_inverse((*inout).scalars[i]); + (*inout).scalars[i].denominator*=-(*inout).base[i]; + (*inout).scalars[i].numerator*=-1; + } + else{ + fprintf(stderr,"error: attempting to invert 0\n"); + exit(-1); + } + } + return(0); +} +// write to output +int number_inverse(Number input, Number* output){ + number_cpy(input,output); + number_inverse_inplace(output); + return(0); +} +// return result +Number number_inverse_ret(Number x){ + Number ret; + number_inverse(x,&ret); + return(ret); +} + +// quotient +int number_quot(Number x1, Number x2, Number* output){ + Number inv; + number_inverse(x2, &inv); + number_prod(x1, inv, output); + free_Number(inv); + return(0); +} +int number_quot_chain(Number x1, Number* inout){ + number_inverse_inplace(inout); + number_prod_chain(x1, inout); + return(0); +} +Number number_quot_ret(Number x1, Number x2){ + Number ret; + number_quot(x1, x2, &ret); + return(ret); +} + + +// remove 0's +int number_simplify(Number in, Number* out){ + int i; + init_Number(out, in.length); + for(i=0;i0){ + char_array_snprintf(out," + "); + } + if(simp.length>1 || (simp.length==1 && simp.base[0]!=1)){ + char_array_append('(',out); + } + Q_sprint(simp.scalars[i], out); + if(simp.length>1 || (simp.length==1 && simp.base[0]!=1)){ + char_array_append(')',out); + } + + if(simp.base[i]!=1){ + char_array_snprintf(out,"s{%d}",simp.base[i]); + } + } + + free_Number(simp); + return(0); +} + +// print to stdout +int number_print(Number number){ + Char_Array buffer; + init_Char_Array(&buffer,5*number.length); + number_sprint(number, &buffer); + printf("%s",buffer.str); + return(0); +} + +#define PP_NULL_MODE 0 +#define PP_NUM_MODE 1 +#define PP_SQRT_MODE 2 +// read from a string +int str_to_Number(char* str, Number* number){ + char* ptr; + int mode; + char* buffer=calloc(str_len(str)+1,sizeof(char)); + char* buffer_ptr=buffer; + Q num; + int base; + // whether there are parentheses in the string + int exist_parenthesis=0; + + init_Number(number, NUMBER_SIZE); + + // init num and base + // init to 0 so that if str is empty, then the number is set to 0 + num=quot(0,1); + base=1; + + mode=PP_NULL_MODE; + for(ptr=str;*ptr!='\0';ptr++){ + switch(*ptr){ + // read number + case '(': + if(mode==PP_NULL_MODE){ + // init base + base=1; + mode=PP_NUM_MODE; + exist_parenthesis=1; + } + break; + case ')': + if(mode==PP_NUM_MODE){ + str_to_Q(buffer,&num); + buffer_ptr=buffer; + *buffer_ptr='\0'; + mode=PP_NULL_MODE; + } + break; + + // read sqrt + case '{': + // init num + if(num.numerator==0){ + num=quot(1,1); + } + if(mode==PP_NULL_MODE){ + mode=PP_SQRT_MODE; + } + // if there is a square root, then do not read a fraction (see end of loop) + exist_parenthesis=1; + break; + case '}': + if(mode==PP_SQRT_MODE){ + sscanf(buffer,"%d",&base); + buffer_ptr=buffer; + *buffer_ptr='\0'; + mode=PP_NULL_MODE; + } + break; + + // write num + case '+': + if(mode==PP_NULL_MODE){ + number_add_elem(num, base, number); + // re-init num and base + num=quot(0,1); + base=1; + } + break; + + default: + if(mode!=PP_NULL_MODE){ + buffer_ptr=str_addchar(buffer_ptr,*ptr); + } + break; + } + } + + // last step + if(mode==PP_NULL_MODE){ + if(exist_parenthesis==0){ + str_to_Q(str, &num); + } + number_add_elem(num, base, number); + } + + free(buffer); + return(0); +} + +// with Char_Array input +int char_array_to_Number(Char_Array cstr_num,Number* output){ + char* buffer; + char_array_to_str(cstr_num,&buffer); + str_to_Number(buffer, output); + free(buffer); + return(0); +} + + +// -------------------- Number_Matrix --------------------- + +// init +int init_Number_Matrix(Number_Matrix* matrix, int length){ + int i,j; + (*matrix).matrix=calloc(length,sizeof(Number*)); + (*matrix).indices=calloc(length,sizeof(int)); + for(i=0;i +#include + +// pre-compiler definitions +#include "definitions.cpp" + +// rccs +#include "rcc.h" +// grouped representation of polynomials +#include "grouped_polynomial.h" +// command line parser +#include "cli_parser.h" +// parse input file +#include "parse_file.h" +// numerical flow +#include "flow.h" +// arrays +#include "array.h" + +// read cli arguments +int read_args_numkondo(int argc,const char* argv[], Str_Array* str_args, Numkondo_Options* opts); +// print usage message +int print_usage_numkondo(); +// compute flow +int numflow(Str_Array str_args, Numkondo_Options opts); + + +int main (int argc, const char* argv[]){ + // string arguments + Str_Array str_args; + // options + Numkondo_Options opts; + + // read command-line arguments + read_args_numkondo(argc,argv,&str_args,&opts); + + numflow(str_args, opts); + + //free memory + free_Str_Array(str_args); + return(0); +} + + +// parse command-line arguments +#define CP_FLAG_NITER 1 +#define CP_FLAG_TOL 2 +#define CP_FLAG_RCCS 3 +int read_args_numkondo(int argc,const char* argv[], Str_Array* str_args, Numkondo_Options* opts){ + int i; + // pointers + char* ptr; + // file to read the polynomial from in flow mode + const char* file=""; + // flag that indicates what argument is being read + int flag=0; + // whether a file was specified on the command-line + int exists_file=0; + + // if there are no arguments + if(argc==1){ + print_usage_numkondo(); + exit(-1); + } + + // defaults + // display entire flow + (*opts).display_mode=DISPLAY_NUMERICAL; + // default niter + (*opts).niter=100; + // default to 0 tolerance + (*opts).tol=0; + // mark rccstring so that it can be recognized whether it has been set or not + (*opts).eval_rccstring.length=-1; + +// loop over arguments +for(i=1;i\n\n"); + return(0); +} + + +// numerical computation of the flow +int numflow(Str_Array str_args, Numkondo_Options opts){ + // index of the entry in the input file + int arg_index; + // list of rccs + Labels labels; + // initial condition + RCC init_cd; + // flow equation + Grouped_Polynomial flow_equation; + + // parse id table + arg_index=find_str_arg("labels", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no labels entry in the configuration file\n"); + exit(-1); + } + else{ + parse_labels(str_args.strs[arg_index], &labels); + } + + // parse flow equation + arg_index=find_str_arg("flow_equation", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no flow equation entry in the configuration file\n"); + exit(-1); + } + else{ + char_array_to_Grouped_Polynomial(str_args.strs[arg_index], &flow_equation); + } + + // initial conditions + // check they were not specified on the command line + if(opts.eval_rccstring.length==-1){ + arg_index=find_str_arg("initial_condition", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no initial condition in the configuration file or on the command line\n"); + exit(-1); + } + else{ + char_array_cpy(str_args.strs[arg_index],&(opts.eval_rccstring)); + } + } + // initialize the rccs + prepare_init(flow_equation.indices,flow_equation.length,&init_cd); + // read rccs from string + if(opts.eval_rccstring.length!=-1){ + parse_init_cd(opts.eval_rccstring, &init_cd); + free_Char_Array(opts.eval_rccstring); + } + + numerical_flow(flow_equation, init_cd, labels, opts.niter, opts.tol, opts.display_mode); + + free_RCC(init_cd); + + // free memory + free_Labels(labels); + free_Grouped_Polynomial(flow_equation); + return(0); +} + diff --git a/src/parse_file.c b/src/parse_file.c new file mode 100644 index 0000000..6054372 --- /dev/null +++ b/src/parse_file.c @@ -0,0 +1,796 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "parse_file.h" + +#include +#include +#include "array.h" +#include "fields.h" +#include "rational.h" +#include "number.h" +#include "polynomial.h" +#include "rcc.h" +#include "definitions.cpp" +#include "istring.h" +#include "tools.h" +#include "idtable.h" + + +// parsing modes +#define PP_NULL_MODE 0 +// when reading a factor +#define PP_FACTOR_MODE 1 +// reading a monomial +#define PP_MONOMIAL_MODE 2 +// reading a numerator and denominator +#define PP_NUMBER_MODE 3 +// types of fields +#define PP_FIELD_MODE 6 +#define PP_PARAMETER_MODE 7 +#define PP_EXTERNAL_MODE 8 +#define PP_INTERNAL_MODE 9 +#define PP_FERMIONS_MODE 10 +// indices +#define PP_INDEX_MODE 11 +// factors or monomials +#define PP_BRACKET_MODE 12 +// labels +#define PP_LABEL_MODE 13 +// polynomial +#define PP_POLYNOMIAL_MODE 14 +// group +#define PP_GROUP_MODE 15 + + +// parse fields list +int parse_input_fields(Char_Array str_fields, Fields_Table* fields){ + // buffer + char* buffer=calloc(str_fields.length+1,sizeof(char)); + char* buffer_ptr=buffer; + int i,j; + int mode; + int comment=0; + + // allocate memory + init_Fields_Table(fields); + + // loop over input + mode=PP_NULL_MODE; + for(j=0;j0){ + str_to_Polynomial(buffer, &polynomial); + symbols_append_noinit(index, polynomial, &((*fields).symbols)); + } + + // simplify + for(i=0;i<(*fields).symbols.length;i++){ + polynomial_simplify((*fields).symbols.expr+i, *fields); + } + + free(buffer); + return(0); +} + +// parse groups of independent fields +int parse_input_groups(Char_Array str_groups, Groups* groups){ + // buffer + char* buffer=calloc(str_groups.length+1,sizeof(char)); + char* buffer_ptr=buffer; + int index; + int j; + Int_Array group; + int mode; + int comment=0; + + // alloc + init_Groups(groups, GROUP_SIZE); + + // loop over input + mode=PP_NULL_MODE; + for(j=0;j=0 && index2>=0){ + free_Polynomial((*propagator).matrix[index1][index2]); + str_to_Polynomial(buffer,(*propagator).matrix[index1]+index2); + buffer_ptr=buffer; + *buffer_ptr='\0'; + mode=PP_INDEX_MODE; + } + break; + + // comment + case '#': + comment=1; + break; + + default: + buffer_ptr=str_addchar(buffer_ptr,str_propagator.str[j]); + break; + } + } + } + + // last step + if(mode==PP_POLYNOMIAL_MODE){ + free_Polynomial((*propagator).matrix[index1][index2]); + str_to_Polynomial(buffer,(*propagator).matrix[index1]+index2); + } + + free(buffer); + return(0); +} + + +// parse input polynomial +int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields_Table fields){ + int j; + // buffer + char* buffer=calloc(str_polynomial.length+1,sizeof(char)); + char* buffer_ptr=buffer; + Polynomial tmp_poly; + + // allocate memory + init_Polynomial(output,POLY_SIZE); + + for(j=0;j-DOFFSET){ + (*init).values[i]=1.; + } + else{ + (*init).values[i]=0.; + } + + } + return(0); +} diff --git a/src/parse_file.h b/src/parse_file.h new file mode 100644 index 0000000..b51103a --- /dev/null +++ b/src/parse_file.h @@ -0,0 +1,56 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* +Parse the input file +*/ + +#ifndef PARSE_FILE_H +#define PARSE_FILE_H + +#include "types.h" + +// parse fields list +int parse_input_fields(Char_Array str_fields, Fields_Table* fields); + +// parse symbols list +int parse_input_symbols(Char_Array str_symbols, Fields_Table* fields); + +// parse groups of independent fields +int parse_input_groups(Char_Array str_groups, Groups* groups); + +// parse identities between fields +int parse_input_identities(Char_Array str_identities, Fields_Table* fields); + +// parse propagator +int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propagator, Fields_Table fields); + +// parse input polynomial +int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields_Table fields); + +// parse id table +int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table fields); + +// parse a list of labels +int parse_labels(Char_Array str_labels, Labels* labels); + +// parse the initial condition +int parse_init_cd(Char_Array init_cd, RCC* init); + +// set indices and length of init +int prepare_init(int* indices, int length, RCC* init); + +#endif diff --git a/src/polynomial.c b/src/polynomial.c new file mode 100644 index 0000000..639728a --- /dev/null +++ b/src/polynomial.c @@ -0,0 +1,1263 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "polynomial.h" + +#include +#include +#include "definitions.cpp" +#include "rational.h" +#include "tools.h" +#include "mean.h" +#include "coefficient.h" +#include "istring.h" +#include "array.h" +#include "number.h" +#include "fields.h" + + +// allocate memory +int init_Polynomial(Polynomial* polynomial,int size){ + (*polynomial).monomials=calloc(size,sizeof(Int_Array)); + (*polynomial).factors=calloc(size,sizeof(Int_Array)); + (*polynomial).nums=calloc(size,sizeof(Number)); + (*polynomial).length=0; + (*polynomial).memory=size; + return(0); +} + +// free memory +int free_Polynomial(Polynomial polynomial){ + int i; + for(i=0;i=(*output).memory){ + resize_polynomial(output,2*(*output).memory+1); + } + + // copy and allocate + int_array_cpy(monomial,(*output).monomials+offset); + int_array_cpy(factor,(*output).factors+offset); + number_cpy(num,(*output).nums+offset); + //increment length + (*output).length++; + + return(0); +} +// append an element to a polynomial without allocating memory +int polynomial_append_noinit(Int_Array monomial, Int_Array factor, Number num, Polynomial* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_polynomial(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).monomials[offset]=monomial; + (*output).factors[offset]=factor; + (*output).nums[offset]=num; + // increment length + (*output).length++; + return(0); +} +// noinit, and if there already exists an element with the same monomial and factor, then just add numbers +int polynomial_append_noinit_inplace(Int_Array monomial, Int_Array factor, Number num, Polynomial* output){ + int i; + int foundit=0; + for(i=0;i<(*output).length;i++){ + if(int_array_cmp(monomial, (*output).monomials[i])==0 && int_array_cmp(factor, (*output).factors[i])==0){ + number_add_chain(num,(*output).nums+i); + foundit=1; + free_Number(num); + free_Int_Array(monomial); + free_Int_Array(factor); + break; + } + } + if(foundit==0){ + polynomial_append_noinit(monomial, factor, num, output); + } + return(0); +} + +// concatenate two polynomials +int polynomial_concat(Polynomial input, Polynomial* output){ + int i; + for(i=0;i=input_polynomial.length;i--){ + // try to advance pointer + (current_term[i])++; + // if the step fails, keep on looping, if not, set all the following + // pointers to the latest position + if(current_term[i]1){ + if(power>33){ + fprintf(stderr,"error: trying to take a power of a polynomial that is too high (>33)\n"); + exit(-1); + } + // next power + polynomial_prod(input_polynomial,previous_power,&tmp_poly, fields); + + // free + free_Polynomial(previous_power); + } + else{ + polynomial_cpy(input_polynomial,&tmp_poly); + } + + // if the power is high enough that V^p=0, then stop + if(tmp_poly.length==0){ + free_Polynomial(tmp_poly); + break; + } + + // copy for next power + polynomial_cpy(tmp_poly,&previous_power); + + // 1/p! + polynomial_multiply_Qscalar(tmp_poly,quot(1,factorial(power))); + // append tmp to the output + polynomial_concat_noinit(tmp_poly,output); + + // increase power + power++; + } + + return(0); +} + + +// log(1+W) +int polynomial_logarithm(Polynomial input_polynomial,Polynomial* output, Fields_Table fields){ + // a buffer for the result of a given power + Polynomial tmp_poly; + // a buffer for the previous power + Polynomial previous_power; + // power + int power=1; + + // allocate memory + init_Polynomial(output,POLY_SIZE); + + while(1){ + if(power>1){ + // next power + polynomial_prod(input_polynomial,previous_power,&tmp_poly, fields); + + // free + free_Polynomial(previous_power); + } + else{ + polynomial_cpy(input_polynomial,&tmp_poly); + } + + // if the power is high enough that V^p=0, then stop + if(tmp_poly.length==0){ + free_Polynomial(tmp_poly); + break; + } + + // copy for next power + polynomial_cpy(tmp_poly,&previous_power); + + // (-1)^{p-1}/p + polynomial_multiply_Qscalar(tmp_poly,quot(ipower(-1,power-1),power)); + // append tmp to the output + polynomial_concat_noinit(tmp_poly,output); + + // increase power + power++; + } + + return(0); +} + +// check whether a monomial vanishes +int check_monomial(Int_Array monomial, Fields_Table fields){ + int i,j; + for(i=0;i0){ + match++; + } + else if(monomial.values[i]<0){ + match--; + } + } + + // check for repetitions + if(is_fermion(monomial.values[i], fields)==1){ + for(j=i+1;j0){ + match_internals++; + } + else if((*polynomial).monomials[i].values[j]<0){ + match_internals--; + } + } + // don't remove a term containing symbols + else if(type==FIELD_SYMBOL){ + match_internals=0; + break; + } + } + if(match_internals==0){ + polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], &output); + } + } + + free_Polynomial(*polynomial); + *polynomial=output; + return(0); +} + + +// denominator of a multinomal factor: m1!m2!... +// requires terms to be sorted +int multinomial(int power,int* terms){ + int multiple=1; + int ret=1; + int i; + // the number of numbers to be multiplied in the multinomial is + // equal to power-1 (the first is 1) + for(i=1;i=(*polynomial).length-1 || monomial_cmp!=0 || factor_cmp!=0 ){ + // check that the polynomial is not trivial + if(number_is_zero(new_num)!=1){ + polynomial_append((*polynomial).monomials[i],(*polynomial).factors[i],new_num,&output); + } + + // reset new numerical factor + free_Number(new_num); + init_Number(&new_num,NUMBER_SIZE); + } + } + + free_Number(new_num); + free_Polynomial(*polynomial); + *polynomial=output; + return(0); +} + +// sort a polynomial +// requires the monomials and factors to be sorted +int polynomial_sort(Polynomial* polynomial, int begin, int end){ + int i; + int index; + int monomial_cmp; + int factor_cmp; + + // the pivot: middle of the array + int pivot=(begin+end)/2; + // if the array is non trivial + if(begin type2){ + return(1); + } + + if(monomial.values[pos1] < monomial.values[pos2]){ + return(-1); + } + else if (monomial.values[pos1] > monomial.values[pos2]){ + return(1); + } + + return(0); +} + +// exchange two fields (with sign) +int exchange_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table fields, int* sign){ + int tmp=monomial.values[pos1]; + int f1,f2; + int i; + + monomial.values[pos1]=monomial.values[pos2]; + monomial.values[pos2]=tmp; + + // sign change + // only if the exchange is not trivial + if(pos1!=pos2){ + f1=is_fermion(monomial.values[pos1],fields); + f2=is_fermion(monomial.values[pos2],fields); + // if both Fermions then sign + if(f1==1 && f2==1){ + *sign*=-1; + } + // if only one of them is a Fermion, then count the number of Fermions between them + else if(f1==1 || f2==1){ + for(i=min(pos1,pos2)+1;i group2){ + return(1); + } + + int type1=field_type(monomial.values[pos1], fields); + int type2=field_type(monomial.values[pos2],fields); + + if(type1 < type2){ + return(-1); + } + else if(type1 > type2){ + return(1); + } + + if(monomial.values[pos1] < monomial.values[pos2]){ + return(-1); + } + else if (monomial.values[pos1] > monomial.values[pos2]){ + return(1); + } + + return(0); +} + + +// convert and idtable to a polynomial +int idtable_to_polynomial(Id_Table idtable, Polynomial* polynomial){ + int i,j; + int start=0; + + // allocate memory + init_Polynomial(polynomial,POLY_SIZE); + + for(i=0;i0){ + // initialize coef to store the product of factors + init_Coefficient(&coef, POLY_SIZE); + + // first term + index=intlist_find_err(equations.indices, equations.length, (*polynomial).factors[i].values[0]); + if(index>=0){ + coefficient_cpy_noinit(equations.coefs[index],&coef); + } + + // other terms + for(j=1;j<(*polynomial).factors[i].length;j++){ + index=intlist_find_err(equations.indices, equations.length, (*polynomial).factors[i].values[j]); + if(index>=0){ + coefficient_prod_chain(equations.coefs[index],&coef); + } + } + + // new polynomial terms + for(j=0;j +#include +#include "istring.h" +#include "array.h" +#include "math.h" + +Q quot(long double p, long double q){ + Q ret; + if(q==0){ + fprintf(stderr,"error: %Lf/%Lf is ill defined\n",p,q); + exit(-1); + } + ret.numerator=p; + ret.denominator=q; + return(ret); +} + +// add +Q Q_add(Q x1,Q x2){ + Q ret; + ret.denominator=lcm(x1.denominator,x2.denominator); + ret.numerator=x1.numerator*(ret.denominator/x1.denominator)+x2.numerator*(ret.denominator/x2.denominator); + return(Q_simplify(ret)); +} +//multiply +Q Q_prod(Q x1,Q x2){ + return(Q_simplify(quot(x1.numerator*x2.numerator,x1.denominator*x2.denominator))); +} +// inverse +Q Q_inverse(Q x1){ + if(x1.numerator>0){ + return(quot(x1.denominator,x1.numerator)); + } + else if(x1.numerator<0){ + return(quot(-x1.denominator,-x1.numerator)); + } + else{ + fprintf(stderr,"error: attempting to invert %Lf/%Lf\n",x1.numerator, x1.denominator); + exit(-1); + } + +} +// quotient +Q Q_quot(Q x1, Q x2){ + if(x2.numerator>0){ + return(Q_simplify(quot(x1.numerator*x2.denominator,x1.denominator*x2.numerator))); + } + else if(x2.numerator<0){ + return(Q_simplify(quot(-x1.numerator*x2.denominator,-x1.denominator*x2.numerator))); + } + else{ + fprintf(stderr,"error: attempting to invert %Lf/%Lf\n",x2.numerator, x2.denominator); + exit(-1); + } +} + +// compare +int Q_cmp(Q x1, Q x2){ + Q quo=Q_quot(x1,x2); + if(quo.numerator > quo.denominator){ + return(1); + } + else if(quo.numerator < quo.denominator){ + return(-1); + } + else{ + return(0); + } +} + +// simplify +Q Q_simplify(Q x){ + return(quot(x.numerator/gcd(x.numerator,x.denominator),x.denominator/gcd(x.numerator,x.denominator))); +} +//simplify in place +int Q_simplify_inplace(Q* x){ + Q ret=Q_simplify(*x); + *x=ret; + return(0); +} + +// greatest common divisor +long double gcd(long double x, long double y){ + long double ax=fabsl(x); + long double ay=fabsl(y); + int security=0; + while(ax!=0 && ay!=0){ + if(ax>ay){ax=modld(ax,ay);} + else{ay=modld(ay,ax);} + + security++; + if(security>1000000){ + fprintf(stderr,"error: could not compute gcd( %Lf , %Lf ) after %d tries\n",x,y,security); + exit(-1); + } + } + // if both are negative, gcd should be negative (useful for simplification of fractions and product of square roots) + if(x<0 && y<0){ + ax*=-1; + ay*=-1; + } + if(fabsl(ay)>fabsl(ax)){return(ay);} + else{return(ax);} +} + +long double modld(long double x, long double m){ + long double q=floorl(x/m); + return(x-q*m); +} + +// least common multiple +long double lcm(long double x,long double y){ + if(x!=0 || y!=0){ + return((x*y)/gcd(x,y)); + } + else{ + return(0); + } +} + +// approximate value as double +double Q_double_value(Q q){ + return(1.0*q.numerator/q.denominator); +} + + +// print to string +int Q_sprint(Q num, Char_Array* out){ + if(num.denominator!=1){ + char_array_snprintf(out,"%Lf/%Lf", num.numerator,num.denominator); + } + else{ + char_array_snprintf(out,"%Lf",num.numerator); + } + + return(0); +} + +#define PP_NUMERATOR_MODE 1 +#define PP_DENOMINATOR_MODE 2 +// read from a string +int str_to_Q(char* str, Q* num){ + char* ptr; + int mode; + char* buffer=calloc(str_len(str)+1,sizeof(char)); + char* buffer_ptr=buffer; + + *num=quot(0,1); + + mode=PP_NUMERATOR_MODE; + for(ptr=str;*ptr!='\0';ptr++){ + if(*ptr=='/'){ + sscanf(buffer,"%Lf",&((*num).numerator)); + buffer_ptr=buffer; + *buffer_ptr='\0'; + mode=PP_DENOMINATOR_MODE; + } + else{ + buffer_ptr=str_addchar(buffer_ptr,*ptr); + } + } + + // last step + if(mode==PP_NUMERATOR_MODE){ + sscanf(buffer,"%Lf",&((*num).numerator)); + } + else if(mode==PP_DENOMINATOR_MODE){ + sscanf(buffer,"%Lf",&((*num).denominator)); + } + + free(buffer); + return(0); +} + +#else +int null_declaration_so_that_the_compiler_does_not_complain; +#endif diff --git a/src/rational_float.h b/src/rational_float.h new file mode 100644 index 0000000..931b5ec --- /dev/null +++ b/src/rational_float.h @@ -0,0 +1,64 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* + Rational numbers + uses long doubles to represent integers (to avoid overflow) +*/ + +#ifdef RATIONAL_AS_FLOAT + +#ifndef RATIONAL_H +#define RATIONAL_H + +#include "types.h" + +// Q from int/int +Q quot(long double p, long double q); + +// add +Q Q_add(Q x1,Q x2); +//multiply +Q Q_prod(Q x1,Q x2); +// inverse +Q Q_inverse(Q x1); +// quotient +Q Q_quot(Q x1, Q x2); + +// compare +int Q_cmp(Q x1, Q x2); + +// simplify +Q Q_simplify(Q x); +//simplify in place +int Q_simplify_inplace(Q* x); + +// greatest common divisor +long double gcd(long double x,long double y); +long double modld(long double x, long double m); +// least common multiple +long double lcm(long double x,long double y); + +// approximate value as double +double Q_double_value(Q q); + +// print to string +int Q_sprint(Q num, Char_Array* out); +// read from a string +int str_to_Q(char* str, Q* num); + +#endif +#endif diff --git a/src/rational_int.c b/src/rational_int.c new file mode 100644 index 0000000..ec601a8 --- /dev/null +++ b/src/rational_int.c @@ -0,0 +1,190 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#ifndef RATIONAL_AS_FLOAT + +#include "rational_int.h" +#include +#include +#include "istring.h" +#include "array.h" + +Q quot(long int p, long int q){ + Q ret; + if(q==0){ + fprintf(stderr,"error: %ld/%ld is ill defined\n",p,q); + exit(-1); + } + ret.numerator=p; + ret.denominator=q; + return(ret); +} + +// add +Q Q_add(Q x1,Q x2){ + Q ret; + ret.denominator=lcm(x1.denominator,x2.denominator); + ret.numerator=x1.numerator*(ret.denominator/x1.denominator)+x2.numerator*(ret.denominator/x2.denominator); + return(Q_simplify(ret)); +} +//multiply +Q Q_prod(Q x1,Q x2){ + return(Q_simplify(quot(x1.numerator*x2.numerator,x1.denominator*x2.denominator))); +} +// inverse +Q Q_inverse(Q x1){ + if(x1.numerator>0){ + return(quot(x1.denominator,x1.numerator)); + } + else if(x1.numerator<0){ + return(quot(-x1.denominator,-x1.numerator)); + } + else{ + fprintf(stderr,"error: attempting to invert %ld/%ld\n",x1.numerator, x1.denominator); + exit(-1); + } + +} +// quotient +Q Q_quot(Q x1, Q x2){ + if(x2.numerator>0){ + return(Q_simplify(quot(x1.numerator*x2.denominator,x1.denominator*x2.numerator))); + } + else if(x2.numerator<0){ + return(Q_simplify(quot(-x1.numerator*x2.denominator,-x1.denominator*x2.numerator))); + } + else{ + fprintf(stderr,"error: attempting to invert %ld/%ld\n",x2.numerator, x2.denominator); + exit(-1); + } +} + +// compare +int Q_cmp(Q x1, Q x2){ + Q quo=Q_quot(x1,x2); + if(quo.numerator > quo.denominator){ + return(1); + } + else if(quo.numerator < quo.denominator){ + return(-1); + } + else{ + return(0); + } +} + +// simplify +Q Q_simplify(Q x){ + return(quot(x.numerator/gcd(x.numerator,x.denominator),x.denominator/gcd(x.numerator,x.denominator))); +} +//simplify in place +int Q_simplify_inplace(Q* x){ + Q ret=Q_simplify(*x); + *x=ret; + return(0); +} + +// greatest common divisor +long int gcd(long int x, long int y){ + long int ax=labs(x); + long int ay=labs(y); + int security=0; + while(ax!=0 && ay!=0){ + if(ax>ay){ax=ax%ay;} + else{ay=ay%ax;} + + security++; + if(security>1000000){ + fprintf(stderr,"error: could not compute gcd( %ld , %ld ) after %d tries\n",x,y,security); + exit(-1); + } + } + // if both are negative, gcd should be negative (useful for simplification of fractions and product of square roots) + if(x<0 && y<0){ + ax*=-1; + ay*=-1; + } + if(labs(ay)>labs(ax)){return(ay);} + else{return(ax);} +} + +// least common multiple +long int lcm(long int x,long int y){ + if(x!=0 || y!=0){ + return((x*y)/gcd(x,y)); + } + else{ + return(0); + } +} + +// approximate value as double +double Q_double_value(Q q){ + return(1.0*q.numerator/q.denominator); +} + + +// print to string +int Q_sprint(Q num, Char_Array* out){ + if(num.denominator!=1){ + char_array_snprintf(out,"%ld/%ld", num.numerator,num.denominator); + } + else{ + char_array_snprintf(out,"%ld",num.numerator); + } + + return(0); +} + +#define PP_NUMERATOR_MODE 1 +#define PP_DENOMINATOR_MODE 2 +// read from a string +int str_to_Q(char* str, Q* num){ + char* ptr; + int mode; + char* buffer=calloc(str_len(str)+1,sizeof(char)); + char* buffer_ptr=buffer; + + *num=quot(0,1); + + mode=PP_NUMERATOR_MODE; + for(ptr=str;*ptr!='\0';ptr++){ + if(*ptr=='/'){ + sscanf(buffer,"%ld",&((*num).numerator)); + buffer_ptr=buffer; + *buffer_ptr='\0'; + mode=PP_DENOMINATOR_MODE; + } + else{ + buffer_ptr=str_addchar(buffer_ptr,*ptr); + } + } + + // last step + if(mode==PP_NUMERATOR_MODE){ + sscanf(buffer,"%ld",&((*num).numerator)); + } + else if(mode==PP_DENOMINATOR_MODE){ + sscanf(buffer,"%ld",&((*num).denominator)); + } + + free(buffer); + return(0); +} + +#else +int null_declaration_so_that_the_compiler_does_not_complain; +#endif diff --git a/src/rational_int.h b/src/rational_int.h new file mode 100644 index 0000000..a9f164d --- /dev/null +++ b/src/rational_int.h @@ -0,0 +1,59 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* Rational numbers */ + +#ifndef RATIONAL_AS_FLOAT +#ifndef RATIONAL_H +#define RATIONAL_H + +#include "types.h" + +// Q from int/int +Q quot(long int p, long int q); + +// add +Q Q_add(Q x1,Q x2); +//multiply +Q Q_prod(Q x1,Q x2); +// inverse +Q Q_inverse(Q x1); +// quotient +Q Q_quot(Q x1, Q x2); + +// compare +int Q_cmp(Q x1, Q x2); + +// simplify +Q Q_simplify(Q x); +//simplify in place +int Q_simplify_inplace(Q* x); + +// greatest common divisor +long int gcd(long int x,long int y); +// least common multiple +long int lcm(long int x,long int y); + +// approximate value as double +double Q_double_value(Q q); + +// print to string +int Q_sprint(Q num, Char_Array* out); +// read from a string +int str_to_Q(char* str, Q* num); + +#endif +#endif diff --git a/src/rcc.c b/src/rcc.c new file mode 100644 index 0000000..484e20a --- /dev/null +++ b/src/rcc.c @@ -0,0 +1,95 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "rcc.h" +#include +#include +#include "array.h" + +// init +int init_RCC(RCC* rccs, int size){ + (*rccs).values=calloc(size,sizeof(long double)); + (*rccs).indices=calloc(size,sizeof(int)); + (*rccs).length=size; + return(0); +} + +int free_RCC(RCC rccs){ + free(rccs.values); + free(rccs.indices); + return(0); +} + +// set a given element of an rcc +int RCC_set_elem(long double value, int index, RCC* rcc, int pos){ + (*rcc).values[pos]=value; + (*rcc).indices[pos]=index; + return(0); +} + +int RCC_cpy(RCC input,RCC* output){ + int i; + + init_RCC(output,input.length); + for(i=0;i +#include + +int factorial(int n){ + int i; + int res=1; + for(i=2;i<=n;i++){ + res*=i; + } + return(res); +} + +int ipower(int n, int p){ + int i; + int res=1; + for(i=1;i<=p;i++){ + res*=n; + } + return(res); +} + +// power of a double +long double dpower(long double x, int p){ + int i; + long double res=1.; + if(p>=0){ + for(i=1;i<=p;i++){ + res*=x; + } + } + else{ + for(i=1;i<=-p;i++){ + res/=x; + } + } + + return(res); +} + + +// find an element in a list whose first element is its size +int list_find(int* list, int x){ + int i; + for(i=1;i<=*list;i++){ + if(list[i]==x){return(i);} + } + return(0); +} + +// find an element in a list whose first element is not its number of elements (skips first element) +int unlist_find(int* list, int size, int x){ + int i; + for(i=1;i=x2){ + return(x1); + } + else{ + return(x2); + } +} +int min(int x1, int x2){ + if(x1<=x2){ + return(x1); + } + else{ + return(x2); + } +} diff --git a/src/tools.h b/src/tools.h new file mode 100644 index 0000000..e5f319f --- /dev/null +++ b/src/tools.h @@ -0,0 +1,49 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* +Various tools +*/ + +#ifndef TOOLS_H +#define TOOLS_H + +#include "types.h" + +int factorial(int n); +int ipower(int n, int p); +long double dpower(long double x, int p); + +// find an element in a list whose first element is its size +int list_find(int* list, int x); +// compare two lists +int list_cmp(int* list1, int* list2); + +// find an element in a list whose first element is not its number of elements (skips first element) +int unlist_find(int* list, int size, int x); +// no error +int unlist_find_noerr(int* list, int size, int x); + +// find an element in a list (checks first element) +int intlist_find(int* list, int size, int x); +// with error +int intlist_find_err(int* list, int size, int x); + +// max and min +int max(int x1, int x2); +int min(int x1, int x2); + +#endif diff --git a/src/types.h b/src/types.h new file mode 100644 index 0000000..0452ebc --- /dev/null +++ b/src/types.h @@ -0,0 +1,219 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* + typedefs used by meankondo +*/ + +#ifndef TYPES_H +#define TYPES_H + + +// rational number +typedef struct Q{ +#ifndef RATIONAL_AS_FLOAT + long int numerator; + long int denominator; +#else + long double numerator; + long double denominator; +#endif +} Q; + +// numbers +typedef struct Number{ + Q* scalars; + int* base; + int length; + int memory; +} Number; + +// matrix +typedef struct Number_Matrix{ + Number** matrix; + int* indices; + int length; +} Number_Matrix; + +// list of integers +typedef struct Int_Array{ + int* values; + int length; + int memory; +} Int_Array; + +// string +typedef struct Char_Array{ + char* str; + int length; + int memory; +} Char_Array; + +// string array +typedef struct Str_Array{ + Char_Array* strs; + int length; + int memory; +} Str_Array; + +// polynomial +typedef struct Polynomial{ + Int_Array* monomials; + Int_Array* factors; + Number* nums; + int length; + int memory; +} Polynomial; + +// matrix of polynomial +typedef struct Polynomial_Matrix{ + Polynomial** matrix; + int* indices; + int length; +} Polynomial_Matrix; + +// denominators in coefficients (index of the denominator and the power to which it is raised) +typedef struct coef_denom{ + int index; + int power; +} coef_denom; + +// coefficient +typedef struct Coefficient{ + Int_Array* factors; + Number* nums; + coef_denom* denoms; + int length; + int memory; +} Coefficient; + +// grouped polynomial +typedef struct Grouped_Polynomial{ + Coefficient* coefs; + int* indices; + int length; + int memory; +} Grouped_Polynomial; + +// rcc +typedef struct RCC{ + long double* values; + int* indices; + int length; +} RCC; + +// identities between fields +typedef struct Identities{ + // left hand sides + Int_Array* lhs; + // right hand sides + Polynomial* rhs; + int length; + int memory; +} Identities; + +// symbolic expressions +typedef struct Symbols{ + int* indices; + Polynomial* expr; + int length; + int memory; +} Symbols; + +// groups of independent fields +typedef struct Groups{ + Int_Array* indices; + int length; + int memory; +} Groups; + +// collection of fields +typedef struct Fields_Table{ + // constant parameters + Int_Array parameter; + // anticommuting external fields + Int_Array external; + // anticommuting internal fields + Int_Array internal; + // identities between fields + Identities ids; + // symbolic expressions (commuting) + Symbols symbols; + // list of anti-commuting variables (fields or symbols) + Int_Array fermions; +} Fields_Table; + +// index labels +typedef struct Labels{ + Char_Array* labels; + int* indices; + int length; + int memory; +} Labels; + +// id table +typedef struct Id_Table{ + int* indices; + Polynomial* polynomials; + int length; + int memory; +} Id_Table; + +/* +// polynomial scalar and vectors +typedef struct Polynomial_Scalar{ + Coefficient coef; + int* indices; + int length; +} Polynomial_Scalar; +typedef struct Polynomial_Vector{ + Coefficient* coefv; + int* indices; + int length; +} Polynomial_Vector; +typedef struct Polynomial_Matrix{ + Coefficient** coefm; + int* indices; + int length; +} Polynomial_Matrix; +*/ + + +// command line options +typedef struct Meankondo_Options{ + int threads; + int chain; +} Meankondo_Options; + +typedef struct Numkondo_Options{ + int display_mode; + int niter; + long double tol; + Char_Array eval_rccstring; +} Numkondo_Options; + +typedef struct Meantools_Options{ + int command; + int deriv_derivs; + Int_Array deriv_vars; + Char_Array eval_rccstring; +} Meantools_Options; + +typedef struct Kondopp_Options{ + int dimension; +} Kondopp_Options; + +#endif