software/meankondo
1
Fork 0
Code Issues Pull Requests Releases Activity

Initial commit

Browse Source
This commit is contained in:
Ian Jauslin
2015-06-14 00:52:45 +00:00
commit aa0f3ae298
61 changed files with 13380 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

628
src/array.c Normal file
View File

@ -0,0 +1,628 @@
/*
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 "array.h"
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#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<input.length){
fprintf(stderr,"error: trying to copy an array of length %d to another with memory %d\n",input.length, (*output).memory);
exit(-1);
}
for(i=0;i<input.length;i++){
(*output).values[i]=input.values[i];
}
(*output).length=input.length;
return(0);
}
// resize memory
int int_array_resize(Int_Array* array, int newsize){
Int_Array new_array;
init_Int_Array(&new_array, newsize);
int_array_cpy_noinit(*array,&new_array);
free_Int_Array(*array);
*array=new_array;
return(0);
}
// add a value
int int_array_append(int val, Int_Array* output){
if((*output).length>=(*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;i<input.length;i++){
(*output).values[offset+i]=input.values[i];
}
(*output).length=offset+input.length;
return(0);
}
// find (does not assume the array is sorted)
int int_array_find(int val, Int_Array array){
int i;
for(i=0;i<array.length;i++){
if(array.values[i]==val){
return(i);
}
}
return(-1);
}
int int_array_find_err(int val, Int_Array array){
int i;
for(i=0;i<array.length;i++){
if(array.values[i]==val){
return(i);
}
}
fprintf(stderr,"error: value %d not found in array\n",val);
exit(-1);
}
// sort
int int_array_sort(Int_Array array, int begin, int end){
int i;
int tmp;
int index;
// the pivot: middle of the array
int pivot=(begin+end)/2;
// if the array is non trivial
if(begin<end){
// send pivot to the end
tmp=array.values[pivot];
array.values[pivot]=array.values[end];
array.values[end]=tmp;
// loop over the others
for(i=begin, index=begin;i<end;i++){
// compare with pivot
if(array.values[i]<array.values[end]){
// if smaller, exchange with reference index
tmp=array.values[index];
array.values[index]=array.values[i];
array.values[i]=tmp;
// move reference index
index++;
}
}
// put pivot (which we had sent to the end) in the right place
tmp=array.values[end];
array.values[end]=array.values[index];
array.values[index]=tmp;
// recurse
int_array_sort(array, begin, index-1);
int_array_sort(array, index+1, end);
}
return(0);
}
// compare Int_Array's
int int_array_cmp(Int_Array array1, Int_Array array2){
int i;
// compare lengths
if(array1.length<array2.length){
return(-1);
}
if(array1.length>array2.length){
return(1);
}
// compare terms
for(i=0;i<array1.length;i++){
if(array1.values[i]<array2.values[i]){
return(-1);
}
if(array1.values[i]>array2.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<test_array.length && matches<input.length;i++){
if(input.values[matches]==test_array.values[i]){
matches++;
}
}
if(matches==input.length){
return(1);
}
else{
return(0);
}
}
// print array
int int_array_print(Int_Array array){
int i;
printf("(");
for(i=0;i<array.length-1;i++){
printf("%d,",array.values[i]);
}
printf("%d)",array.values[array.length-1]);
return(0);
}
// read array
int int_array_read(Char_Array str, Int_Array* array){
char* buffer=calloc(str.length+1,sizeof(char));
char* buffer_ptr=buffer;
int i,j;
int comment_mode=0;
// alloc
init_Int_Array(array,str.length);
*buffer_ptr='\0';
// loop over the input
for(j=0;j<str.length;j++){
if(comment_mode==1){
if(str.str[j]=='\n'){
comment_mode=0;
}
}
else{
switch(str.str[j]){
// new term
case ',':
// write
sscanf(buffer,"%d",&i);
int_array_append(i,array);
// reset buffer
buffer_ptr=buffer;
*buffer_ptr='\0';
break;
// ignore
case '(':break;
case ')':break;
// comments
case '#':
comment_mode=1;
break;
default:
// write to buffer
buffer_ptr=str_addchar(buffer_ptr,str.str[j]);
break;
}
}
}
// write
sscanf(buffer,"%d",&i);
int_array_append(i,array);
free(buffer);
return(0);
}
// ------------------- CharArray ------------------------
// init
int init_Char_Array(Char_Array* array, int memory){
(*array).str=calloc(memory,sizeof(char));
(*array).memory=memory;
(*array).length=0;
return(0);
}
int free_Char_Array(Char_Array array){
free(array.str);
return(0);
}
// copy
int char_array_cpy(Char_Array input, Char_Array* output){
init_Char_Array(output,input.length+1);
char_array_cpy_noinit(input,output);
return(0);
}
int char_array_cpy_noinit(Char_Array input, Char_Array* output){
int i;
if((*output).memory<input.length){
fprintf(stderr,"error: trying to copy an array of length %d to another with memory %d\n",input.length, (*output).memory);
exit(-1);
}
for(i=0;i<input.length;i++){
(*output).str[i]=input.str[i];
}
(*output).length=input.length;
return(0);
}
// resize memory
int char_array_resize(Char_Array* array, int newsize){
Char_Array new_array;
init_Char_Array(&new_array, newsize);
char_array_cpy_noinit(*array,&new_array);
free_Char_Array(*array);
*array=new_array;
return(0);
}
// add a value
int char_array_append(char val, Char_Array* output){
if((*output).length>=(*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;i<input.length;i++){
(*output).str[offset+i]=input.str[i];
}
(*output).length=offset+input.length;
return(0);
}
// convert to char*
int char_array_to_str(Char_Array input, char** output){
int i;
(*output)=calloc(input.length+1,sizeof(char));
for(i=0;i<input.length;i++){
(*output)[i]=input.str[i];
}
if((*output)[input.length-1]!='\0'){
(*output)[input.length]='\0';
}
return(0);
}
// noinit (changes the size of input if needed)
char* char_array_to_str_noinit(Char_Array* input){
if((*input).str[(*input).length-1]!='\0'){
if((*input).length==(*input).memory){
char_array_resize(input,(*input).length+1);
}
// add final '\0'
(*input).str[(*input).length]='\0';
}
return((*input).str);
}
// convert from char*
int str_to_char_array(char* str, Char_Array* output){
char* ptr;
init_Char_Array(output, str_len(str));
for(ptr=str;*ptr!='\0';ptr++){
char_array_append(*ptr,output);
}
return(0);
}
// format strings
int char_array_snprintf(Char_Array* output, char* fmt, ...){
size_t size=STR_SIZE;
unsigned int extra_size;
char* out_str=calloc(size,sizeof(char));
char* ptr;
va_list vaptr;
// initialize argument list starting after fmt
va_start(vaptr, fmt);
// print format
extra_size=vsnprintf(out_str, size, fmt, vaptr);
va_end(vaptr);
// if too large
if(extra_size>size){
// 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<str.length;i++){
if(str.str[i]!='%'){
char_array_append(str.str[i],out);
}
else{
char_array_append(c,out);
}
}
return(0);
}
// ------------------- Str_Array ------------------------
// init
int init_Str_Array(Str_Array* array, int memory){
(*array).strs=calloc(memory,sizeof(Char_Array));
(*array).memory=memory;
(*array).length=0;
return(0);
}
int free_Str_Array(Str_Array array){
int i;
for(i=0;i<array.length;i++){
free_Char_Array(array.strs[i]);
}
free(array.strs);
return(0);
}
// copy
int str_array_cpy(Str_Array input, Str_Array* output){
init_Str_Array(output,input.length);
str_array_cpy_noinit(input,output);
return(0);
}
int str_array_cpy_noinit(Str_Array input, Str_Array* output){
int i;
if((*output).memory<input.length){
fprintf(stderr,"error: trying to copy an array of length %d to another with memory %d\n",input.length, (*output).memory);
exit(-1);
}
for(i=0;i<input.length;i++){
char_array_cpy(input.strs[i],(*output).strs+i);
}
(*output).length=input.length;
return(0);
}
// resize memory
int str_array_resize(Str_Array* array, int newsize){
Str_Array new_array;
init_Str_Array(&new_array, newsize);
str_array_cpy_noinit(*array,&new_array);
free_Str_Array(*array);
*array=new_array;
return(0);
}
// add a value
int str_array_append(Char_Array val, Str_Array* output){
if((*output).length>=(*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<input.length;i++){
char_array_cpy(input.strs[i],(*output).strs+offset+i);
}
(*output).length=offset+input.length;
return(0);
}
int str_array_concat_noinit(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<input.length;i++){
(*output).strs[offset+i]=input.strs[i];
}
(*output).length=offset+input.length;
// free input arrays
free(input.strs);
return(0);
}
// ------------------- Labels ------------------------
// allocate memory
int init_Labels(Labels* labels,int size){
(*labels).labels=calloc(size,sizeof(Char_Array));
(*labels).indices=calloc(size,sizeof(int));
(*labels).length=0;
(*labels).memory=size;
return(0);
}
// free memory
int free_Labels(Labels labels){
int i;
for(i=0;i<labels.length;i++){
free_Char_Array(labels.labels[i]);
}
free(labels.labels);
free(labels.indices);
return(0);
}
// resize the memory allocated to a labels table
int resize_labels(Labels* labels,int new_size){
Labels new_labels;
int i;
init_Labels(&new_labels,new_size);
for(i=0;i<(*labels).length;i++){
new_labels.labels[i]=(*labels).labels[i];
new_labels.indices[i]=(*labels).indices[i];
}
new_labels.length=(*labels).length;
free((*labels).labels);
free((*labels).indices);
*labels=new_labels;
return(0);
}
// copy a labels table
int labels_cpy(Labels input, Labels* output){
init_Labels(output,input.length);
labels_cpy_noinit(input,output);
return(0);
}
int labels_cpy_noinit(Labels input, Labels* output){
int i;
if((*output).memory<input.length){
fprintf(stderr,"error: trying to copy a labels collection of length %d to another with memory %d\n",input.length,(*output).memory);
exit(-1);
}
for(i=0;i<input.length;i++){
char_array_cpy(input.labels[i],(*output).labels+i);
(*output).indices[i]=input.indices[i];
}
(*output).length=input.length;
return(0);
}
// append an element to a labels
int labels_append(Char_Array label, int index, Labels* output){
int offset=(*output).length;
if((*output).length>=(*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<input.length;i++){
labels_append(input.labels[i],input.indices[i],output);
}
return(0);
}

134
src/array.h Normal file
View File

@ -0,0 +1,134 @@
/*
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.
*/
/* Array structures */
#ifndef ARRAY_H
#define ARRAY_H
#include "types.h"
// init
int init_Int_Array(Int_Array* array, int memory);
int free_Int_Array(Int_Array array);
// copy
int int_array_cpy(Int_Array input, Int_Array* output);
int int_array_cpy_noinit(Int_Array input, Int_Array* output);
// resize memory
int int_array_resize(Int_Array* array, int newsize);
// add a value
int int_array_append(int val, Int_Array* output);
// concatenate
int int_array_concat(Int_Array input, Int_Array* output);
// find (does not assume the array is sorted)
int int_array_find(int val, Int_Array array);
int int_array_find_err(int val, Int_Array array);
// sort
int int_array_sort(Int_Array array, int begin, int end);
// compare Int_Array's
int int_array_cmp(Int_Array array1, Int_Array array2);
// check whether an array is a sub-array of another
int int_array_is_subarray_ordered(Int_Array input, Int_Array test_array);
// print array
int int_array_print(Int_Array array);
// read array
int int_array_read(Char_Array str, Int_Array* array);
// Char_Array
// init
int init_Char_Array(Char_Array* array, int memory);
int free_Char_Array(Char_Array array);
// copy
int char_array_cpy(Char_Array input, Char_Array* output);
int char_array_cpy_noinit(Char_Array input, Char_Array* output);
// resize memory
int char_array_resize(Char_Array* array, int newsize);
// add a value
int char_array_append(char val, Char_Array* output);
int char_array_append_str(char* str, Char_Array* output);
// concatenate
int char_array_concat(Char_Array input, Char_Array* output);
// convert to char*
int char_array_to_str(Char_Array input, char** output);
// noinit (changes the size of input if needed)
char* char_array_to_str_noinit(Char_Array* input);
// convert from char*
int str_to_char_array(char* str, Char_Array* output);
// format strings
int char_array_snprintf(Char_Array* output, char* fmt, ...);
// replace '*' with given character
int replace_star(char c, Char_Array str, Char_Array* out);
// Str_Array
// init
int init_Str_Array(Str_Array* array, int memory);
int free_Str_Array(Str_Array array);
// copy
int str_array_cpy(Str_Array input, Str_Array* output);
int str_array_cpy_noinit(Str_Array input, Str_Array* output);
// resize memory
int str_array_resize(Str_Array* array, int newsize);
// add a value
int str_array_append(Char_Array val, Str_Array* output);
int str_array_append_noinit(Char_Array val, Str_Array* output);
// concatenate
int str_array_concat(Str_Array input, Str_Array* output);
int str_array_concat_noinit(Str_Array input, Str_Array* output);
// Labels
// allocate memory
int init_Labels(Labels* labels,int size);
// free memory
int free_Labels(Labels labels);
// resize the memory allocated to a labels table
int resize_labels(Labels* labels,int new_size);
// copy a labels table
int labels_cpy(Labels input, Labels* output);
int labels_cpy_noinit(Labels input, Labels* output);
// append an element to a labels table
int labels_append(Char_Array label, int index, Labels* output);
int labels_append_noinit(Char_Array label, int index, Labels* output);
// concatenate two labels tables
int labels_concat(Labels input, Labels* output);
#endif

123
src/cli_parser.c Normal file
View File

@ -0,0 +1,123 @@
/*
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 "cli_parser.h"
#include <stdio.h>
#include <stdlib.h>
#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;j++){
if(str_arg.str[j]=='#' && str_arg.str[j+1]=='!' && (j==0 || str_arg.str[j-1]=='\n')){
break;
}
}
// if no title then error
if(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<str_arg.length && str_arg.str[k]!='\n'; k++){
char_array_append(str_arg.str[k],out);
}
return(0);
}
// find a string argument with the specified title
int find_str_arg(char* title, Str_Array str_args){
int i,k;
char* ptr;
Char_Array buffer;
for(i=0;i<str_args.length;i++){
get_str_arg_title(str_args.strs[i], &buffer);
// match title
for(k=0, ptr=title; k<buffer.length; k++, ptr++){
if(buffer.str[k]!=*ptr){
break;
}
}
// if match
if(k==buffer.length){
free_Char_Array(buffer);
break;
}
free_Char_Array(buffer);
}
// if found
if(i<str_args.length){
return(i);
}
else{
return(-1);
}
}

35
src/cli_parser.h Normal file
View File

@ -0,0 +1,35 @@
/*
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 command-line arguments
*/
#ifndef CLI_PARSER_H
#define CLI_PARSER_H
#include "types.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);
// get the title from a string argument
int get_str_arg_title(Char_Array str_arg, Char_Array* out);
// find a string argument with the specified title
int find_str_arg(char* title, Str_Array str_args);
#endif

739
src/coefficient.c Normal file
View File

@ -0,0 +1,739 @@
/*
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 "coefficient.h"
#include <stdio.h>
#include <stdlib.h>
#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<coef.length;i++){
free_Int_Array(coef.factors[i]);
free_Number(coef.nums[i]);
}
free(coef.factors);
free(coef.nums);
free(coef.denoms);
return(0);
}
// copy a coefficient
int coefficient_cpy(Coefficient input, Coefficient* output){
init_Coefficient(output,input.length);
coefficient_cpy_noinit(input,output);
return(0);
}
int coefficient_cpy_noinit(Coefficient input, Coefficient* output){
int i;
// if output does not have enough memory allocated to it
if(input.length>(*output).memory){
resize_Coefficient(output,input.length);
}
for(i=0;i<input.length;i++){
int_array_cpy(input.factors[i],(*output).factors+i);
number_cpy(input.nums[i],(*output).nums+i);
(*output).denoms[i]=input.denoms[i];
}
(*output).length=input.length;
return(0);
}
// resize the memory allocated to a coefficient
int resize_Coefficient(Coefficient* coefficient,int new_size){
Coefficient new_coef;
int i;
init_Coefficient(&new_coef,new_size);
for(i=0;i<(*coefficient).length;i++){
new_coef.factors[i]=(*coefficient).factors[i];
new_coef.nums[i]=(*coefficient).nums[i];
new_coef.denoms[i]=(*coefficient).denoms[i];
}
new_coef.length=(*coefficient).length;
free((*coefficient).factors);
free((*coefficient).nums);
free((*coefficient).denoms);
*coefficient=new_coef;
return(0);
}
// append an element to a coefficient
int coefficient_append(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 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;i<input.length;i++){
coefficient_append(input.factors[i],input.nums[i],input.denoms[i],output);
}
coefficient_simplify(output);
return(0);
}
int coefficient_concat_noinit(Coefficient input, Coefficient* output){
int i;
for(i=0;i<input.length;i++){
coefficient_append_noinit(input.factors[i],input.nums[i],input.denoms[i],output);
}
coefficient_simplify(output);
// free input arrays
free(input.factors);
free(input.nums);
free(input.denoms);
return(0);
}
// simplify a Coefficient
int coefficient_simplify(Coefficient* coefficient){
int i;
Coefficient output;
init_Coefficient(&output,(*coefficient).length);
// the combination of numerical factors
Number new_num;
init_Number(&new_num,NUMBER_SIZE);
// sort the factors in the coefficient
for(i=0;i<(*coefficient).length;i++){
int_array_sort((*coefficient).factors[i],0,(*coefficient).factors[i].length-1);
}
// in order to perform a simplification, the list of terms must be
// sorted (so that terms that are proportional are next to each other)
sort_coefficient(coefficient,0,(*coefficient).length-1);
for(i=0;i<(*coefficient).length;i++){
// if the term actually exists
if(number_is_zero((*coefficient).nums[i])!=1){
// combine numerical factors
number_add_chain((*coefficient).nums[i],&new_num);
}
// if the number is 0, the previous terms that may have the same factors should still be added, hence the 'if' ends here
// if the factor is different from the next then add term
if(i==(*coefficient).length-1 || (int_array_cmp((*coefficient).factors[i],(*coefficient).factors[i+1])!=0) || coef_denom_cmp((*coefficient).denoms[i],(*coefficient).denoms[i+1])!=0){
// check that the coefficient is not trivial
if(number_is_zero(new_num)!=1){
coefficient_append((*coefficient).factors[i],new_num,(*coefficient).denoms[i],&output);
}
// reset new numerical factor
free_Number(new_num);
init_Number(&new_num,NUMBER_SIZE);
}
}
free_Number(new_num);
free_Coefficient(*coefficient);
*coefficient=output;
return(0);
}
// sort the terms in an equation (quicksort algorithm)
int sort_coefficient(Coefficient* coefficient, int begin, int end){
int i;
int index;
// the pivot: middle of the array
int pivot=(begin+end)/2;
// if the array is non trivial
if(begin<end){
// send pivot to the end
exchange_coefficient_terms(pivot,end,coefficient);
// loop over the others
for(i=begin, index=begin;i<end;i++){
// compare with pivot
if(coef_denom_cmp((*coefficient).denoms[i],(*coefficient).denoms[end])<0 || ( coef_denom_cmp((*coefficient).denoms[i],(*coefficient).denoms[end])==0 && (int_array_cmp((*coefficient).factors[i],(*coefficient).factors[end])<0)) ){
// if smaller, exchange with reference index
exchange_coefficient_terms(i,index,coefficient);
// move reference index
index++;
}
}
// put pivot (which we had sent to the end) in the right place
exchange_coefficient_terms(index,end,coefficient);
// recurse
sort_coefficient(coefficient, begin, index-1);
sort_coefficient(coefficient, index+1, end);
}
return(0);
}
// exchange two terms (for the sorting algorithm)
int exchange_coefficient_terms(int i, int j, Coefficient* coefficient){
Int_Array ptmp;
Number tmpq;
coef_denom tmpc;
ptmp=(*coefficient).factors[j];
(*coefficient).factors[j]=(*coefficient).factors[i];
(*coefficient).factors[i]=ptmp;
tmpq=(*coefficient).nums[j];
(*coefficient).nums[j]=(*coefficient).nums[i];
(*coefficient).nums[i]=tmpq;
tmpc=(*coefficient).denoms[j];
(*coefficient).denoms[j]=(*coefficient).denoms[i];
(*coefficient).denoms[i]=tmpc;
return(0);
}
// derive a coefficient with respect to an index
int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output){
int i,j;
// temp list of indices
Int_Array factor;
// number of times index was found
int match_count;
// whether the output contains at least one factor
int at_least_one=0;
coef_denom denom;
// loop over monomials
for(i=0;i<input.length;i++){
init_Int_Array(&factor,input.factors[i].length);
// init match count
match_count=0;
// loop over indices
for(j=0;j<input.factors[i].length;j++){
// if found
if(input.factors[i].values[j]==index){
// if it's the first match, don't add it
if(match_count!=0){
int_array_append(index,&factor);
}
match_count++;
}
// add the index
else{
int_array_append(input.factors[i].values[j],&factor);
}
}
denom=input.denoms[i];
// if the index is that of 1/C
if(index==input.denoms[i].index){
// if no C in the numerator (which is normal behavior)
if(match_count==0){
denom.power++;
}
match_count-=input.denoms[i].power;
}
// if the derivative doesn't vanish, add it to the coefficient
if(match_count!=0){
at_least_one=1;
coefficient_append_noinit(factor,number_Qprod_ret(quot(match_count,1),input.nums[i]), denom, output);
}
else{
free_Int_Array(factor);
}
}
if(at_least_one==1){
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);
}
int coefficient_deriv(Coefficient input, int index, Coefficient* output){
init_Coefficient(output, COEF_SIZE);
coefficient_deriv_noinit(input, index, output);
return(0);
}
/*
// 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 i;
// temp list of indices
Int_Array factor;
// number of times index was found
int match_count;
// whether the output contains at least one factor
int at_least_one=0;
coef_denom denom;
// loop over monomials
for(i=0;i<input.length;i++){
// derivative of monomial
monomial_deriv(input.factors[i], index, &factor, &match_count);
// if the derivative doesn't vanish, add it to the coefficient
if(match_count>0){
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;j<factor.length;j++){
// if found
if(factor.values[j]==index){
// if it's the first match, don't add it
if(*match_count!=0){
int_array_append(index,out_factor);
}
(*match_count)++;
}
// add the index
else{
int_array_append(factor.values[j],out_factor);
}
}
return(0);
}
*/
// product of two coefficients
int coefficient_prod(Coefficient coef1, Coefficient coef2, Coefficient* output){
int i,j;
// tmp factor
Int_Array factor;
coef_denom denom;
// init
init_Coefficient(output,COEF_SIZE);
// loop over coef1
for(i=0;i<coef1.length;i++){
// loop over coef2
for(j=0;j<coef2.length;j++){
init_Int_Array(&factor,coef1.factors[i].length+coef2.factors[j].length);
int_array_concat(coef1.factors[i],&factor);
int_array_concat(coef2.factors[j],&factor);
// don't throw an error if the power is 0
if(coef2.denoms[i].power==0){
coef2.denoms[i].index=coef1.denoms[i].index;
}
else if(coef1.denoms[i].power==0){
coef1.denoms[i].index=coef2.denoms[i].index;
}
if(coef1.denoms[i].index!=coef2.denoms[j].index){
fprintf(stderr,"error: cannot multiply flow equations with different constants\n");
exit(-1);
}
denom=coef1.denoms[i];
denom.power+=coef2.denoms[j].power;
coefficient_append_noinit(factor,number_prod_ret(coef1.nums[i],coef2.nums[j]), denom, output);
}
}
// simplify output
coefficient_simplify(output);
return(0);
}
// product of coefficients, output replaces the second coefficient
int coefficient_prod_chain(Coefficient in, Coefficient* inout){
Coefficient tmp_coef;
coefficient_prod(in,*inout,&tmp_coef);
free_Coefficient(*inout);
*inout=tmp_coef;
return(0);
}
// print coefficient
// offset specifies the amount of blank space to the left of the terms after the first
// prepend indices by ind_pre
int coefficient_sprint(Coefficient coef, Char_Array* output, int offset, char index_pre){
Char_Array buffer;
int i,j,k;
int dcount;
if(coef.length==0){
char_array_snprintf(output, " (0)\n");
}
for(i=0;i<coef.length;i++){
if(i==0){
char_array_snprintf(output, " ");
}
else{
for(j=0;j<=offset;j++){
char_array_append(' ',output);
}
char_array_append('+',output);
}
// print numerical coefficient
char_array_append('(',output);
init_Char_Array(&buffer, STR_SIZE);
number_sprint(coef.nums[i], &buffer);
char_array_concat(buffer, output);
free_Char_Array(buffer);
char_array_append(')',output);
// print factors
for(j=0;j<coef.factors[i].length;j++){
// constant indices
if(coef.factors[i].values[j]<0){
// count derivatives
dcount=-coef.factors[i].values[j]/DOFFSET;
char_array_append('[',output);
for(k=0;k<dcount;k++){
char_array_append('d',output);
}
char_array_snprintf(output,"C%d]",-coef.factors[i].values[j]-dcount*DOFFSET);
}
else{
// count derivatives
dcount=coef.factors[i].values[j]/DOFFSET;
char_array_append('[',output);
for(k=0;k<dcount;k++){
char_array_append('d',output);
}
char_array_snprintf(output,"%c%d]",index_pre,coef.factors[i].values[j]-dcount*DOFFSET);
}
}
// print constant denominators
if(coef.denoms[i].power!=0){
char_array_snprintf(output,"[/C%d^%d]",-coef.denoms[i].index,coef.denoms[i].power);
}
char_array_append('\n',output);
}
return(0);
}
// read from a string
#define PP_NULL_MODE 0
#define PP_BRACKET_MODE 1
#define PP_INDICES_MODE 2
#define PP_POWER_MODE 3
#define PP_COMMENT_MODE 4
#define PP_NUMBER_MODE 5
#define PP_CONSTANT_MODE 6
#define PP_CONSTANT_DENOM_MODE 7
int char_array_to_Coefficient(Char_Array str_coef, Coefficient* output){
// buffer
char* buffer=calloc(str_coef.length+1,sizeof(char));
char* buffer_ptr=buffer;
Int_Array indices;
coef_denom denom;
Number num, tmp1_num;
int mode;
int i,j;
int parenthesis_count=0;
int dcount=0;
// allocate memory
init_Coefficient(output,COEF_SIZE);
// init
init_Int_Array(&indices, MONOMIAL_SIZE);
num=number_one();
denom.index=-1;
denom.power=0;
*buffer_ptr='\0';
// loop over the input polynomial
// start in null mode
mode=PP_NULL_MODE;
for(j=0;j<str_coef.length;j++){
if(mode==PP_COMMENT_MODE){
if(str_coef.str[j]=='\n'){
mode=PP_NULL_MODE;
}
}
else{
switch(str_coef.str[j]){
// new indices
case '+':
if(mode==PP_NULL_MODE){
coefficient_append_noinit(indices, num, denom, output);
// reset indices, num
init_Int_Array(&indices, MONOMIAL_SIZE);
num=number_one();
denom.index=-1;
denom.power=0;
}
break;
// enter indices or power mode
case '[':
if(mode==PP_NULL_MODE){
mode=PP_BRACKET_MODE;
// reset derivatives count
dcount=0;
}
break;
// indices mode
case '%':
if(mode==PP_BRACKET_MODE){
mode=PP_INDICES_MODE;
buffer_ptr=buffer;
*buffer_ptr='\0';
}
break;
case 'C':
if(mode==PP_BRACKET_MODE){
mode=PP_CONSTANT_MODE;
buffer_ptr=buffer;
*buffer_ptr='\0';
}
break;
case '/':
if(mode==PP_BRACKET_MODE){
mode=PP_CONSTANT_DENOM_MODE;
buffer_ptr=buffer;
*buffer_ptr='\0';
}
else if(mode!=PP_NULL_MODE){
// write to buffer
buffer_ptr=str_addchar(buffer_ptr,str_coef.str[j]);
}
break;
// derivatives
case 'd':
if(mode==PP_BRACKET_MODE || mode==PP_INDICES_MODE || mode==PP_CONSTANT_MODE){
dcount++;
}
break;
// power mode
case '^':
if(mode==PP_CONSTANT_DENOM_MODE){
sscanf(buffer,"%d",&i);
denom.index=-i;
mode=PP_POWER_MODE;
buffer_ptr=buffer;
*buffer_ptr='\0';
}
else{
buffer_ptr=str_addchar(buffer_ptr,str_coef.str[j]);
}
break;
// read indices or power
case ']':
sscanf(buffer,"%d",&i);
if(mode==PP_INDICES_MODE){
int_array_append(i+dcount*DOFFSET,&indices);
}
else if(mode==PP_CONSTANT_MODE){
int_array_append(-i-dcount*DOFFSET,&indices);
}
else if(mode==PP_POWER_MODE){
denom.power=i;
}
// switch back to null mode
mode=PP_NULL_MODE;
break;
// numerical pre-factor
case '(':
if(mode==PP_NULL_MODE){
mode=PP_NUMBER_MODE;
parenthesis_count=0;
buffer_ptr=buffer;
*buffer_ptr='\0';
}
else if(mode==PP_NUMBER_MODE){
// match parentheses
parenthesis_count++;
buffer_ptr=str_addchar(buffer_ptr,str_coef.str[j]);
}
break;
case ')':
if(mode==PP_NUMBER_MODE){
if(parenthesis_count==0){
// write num
str_to_Number(buffer,&tmp1_num);
number_prod_chain(tmp1_num,&num);
free_Number(tmp1_num);
// back to null mode
mode=PP_NULL_MODE;
}
else{
parenthesis_count--;
buffer_ptr=str_addchar(buffer_ptr,str_coef.str[j]);
}
}
break;
// characters to ignore
case ' ':break;
case '&':break;
case '\n':break;
// comments
case '#':
mode=PP_COMMENT_MODE;
break;
default:
if(mode!=PP_NULL_MODE){
// write to buffer
buffer_ptr=str_addchar(buffer_ptr,str_coef.str[j]);
}
break;
}
}
}
// last term
coefficient_append_noinit(indices, num, denom, output);
free(buffer);
return(0);
}
int str_to_Coefficient(char* str, Coefficient* output){
Char_Array array;
array.length=str_len(str);
array.str=str;
char_array_to_Coefficient(array, output);
return(0);
}
// compare coefficient denominators
int coef_denom_cmp(coef_denom denom1, coef_denom denom2){
if(denom1.index<denom2.index){
return(1);
}
else if(denom1.index>denom2.index){
return(-1);
}
if(denom1.power<denom2.power){
return(-1);
}
else if(denom1.power>denom2.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;i<coef.length;i++){
// product of factors
for(j=0, num_factor=1.;j<coef.factors[i].length;j++){
num_factor*=rccs.values[intlist_find_err(rccs.indices,rccs.length,coef.factors[i].values[j])];
}
// denominator
if(coef.denoms[i].power>0){
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);
}

78
src/coefficient.h Normal file
View File

@ -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

60
src/definitions.cpp Normal file
View File

@ -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

28
src/expansions.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#include "definitions.cpp"
#include "tools.h"
#include "array.h"
#include "polynomial.h"
#include "number.h"
#include "rational.h"

34
src/expansions.h Normal file
View File

@ -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

489
src/fields.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#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<identities.length;i++){
free_Int_Array(identities.lhs[i]);
free_Polynomial(identities.rhs[i]);
}
free(identities.lhs);
free(identities.rhs);
return(0);
}
// resize
int resize_identities(Identities* identities,int new_size){
Identities new_identities;
int i;
init_Identities(&new_identities,new_size);
for(i=0;i<(*identities).length;i++){
new_identities.lhs[i]=(*identities).lhs[i];
new_identities.rhs[i]=(*identities).rhs[i];
}
new_identities.length=(*identities).length;
free((*identities).lhs);
free((*identities).rhs);
*identities=new_identities;
return(0);
}
// copy
int identities_cpy(Identities input, Identities* output){
init_Identities(output,input.length);
identities_cpy_noinit(input,output);
return(0);
}
int identities_cpy_noinit(Identities input, Identities* output){
int i;
if((*output).memory<input.length){
fprintf(stderr,"error: trying to copy an identities collection of length %d to another with memory %d\n",input.length,(*output).memory);
exit(-1);
}
for(i=0;i<input.length;i++){
int_array_cpy(input.lhs[i],(*output).lhs+i);
polynomial_cpy(input.rhs[i],(*output).rhs+i);
}
(*output).length=input.length;
return(0);
}
// append an element to a identities
int identities_append(Int_Array lhs, Polynomial rhs, Identities* output){
int offset=(*output).length;
if((*output).length>=(*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;i<input.length;i++){
identities_append(input.lhs[i],input.rhs[i],output);
}
return(0);
}
// resolve the identities
// requires both the monomials in polynomial and the ids in fields to be sorted
int resolve_ids(Polynomial* polynomial, Fields_Table fields){
int i,j,k,l;
int sign;
int fermion_count;
int at_least_one;
int security;
Int_Array monomial;
Number num;
Number tmp_num;
// loop over monomials
for(i=0;i<(*polynomial).length;i++){
at_least_one=1;
security=0;
// repeat the simplification until the monomial is fully simplified
while(at_least_one>0){
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.length;j++){
// check whether the monomial matches the id
if(int_array_is_subarray_ordered(fields.ids.lhs[j],(*polynomial).monomials[i])==1){
init_Int_Array(&monomial, (*polynomial).monomials[i].length);
// remove lhs from monomial
// sign from moving the fields out of the monomial
sign=1;
// number of Fermions to remove from the monomial
fermion_count=0;
for(k=0,l=0;k<(*polynomial).monomials[i].length;k++){
// check whether the field is identical to the "current" one in the id
// if l is too large, then keep the field
if(l>=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<fields.ids.rhs[j].length;k++){
number_prod(num, fields.ids.rhs[j].nums[k], &tmp_num);
polynomial_append(monomial, (*polynomial).factors[i], tmp_num, polynomial);
free_Number(tmp_num);
int_array_concat(fields.ids.rhs[j].monomials[k],(*polynomial).monomials+(*polynomial).length-1);
// re-sort monomial
sign=1;
monomial_sort((*polynomial).monomials[(*polynomial).length-1],0,(*polynomial).monomials[(*polynomial).length-1].length-1,fields,&sign);
number_Qprod_chain(quot(sign,1),(*polynomial).nums+(*polynomial).length-1);
}
// correct i-th monomial
free_Number((*polynomial).nums[i]);
(*polynomial).nums[i]=number_prod_ret(num,fields.ids.rhs[j].nums[0]);
free_Int_Array((*polynomial).monomials[i]);
(*polynomial).monomials[i]=monomial;
int_array_concat(fields.ids.rhs[j].monomials[0],(*polynomial).monomials+i);
// re-sort monomial
sign=1;
monomial_sort((*polynomial).monomials[i],0,(*polynomial).monomials[i].length-1,fields,&sign);
number_Qprod_chain(quot(sign,1),(*polynomial).nums+i);
// free num
free_Number(num);
// repeat the step (in order to perform all of the replacements if several are necessary)
j--;
if(at_least_one==0){
at_least_one=1;
}
}
}
}
}
return(0);
}
// ------------------ Symbols --------------------
// allocate memory
int init_Symbols(Symbols* symbols,int size){
(*symbols).indices=calloc(size,sizeof(int));
(*symbols).expr=calloc(size,sizeof(Polynomial));
(*symbols).length=0;
(*symbols).memory=size;
return(0);
}
// free memory
int free_Symbols(Symbols symbols){
int i;
for(i=0;i<symbols.length;i++){
free_Polynomial(symbols.expr[i]);
}
free(symbols.indices);
free(symbols.expr);
return(0);
}
// resize
int resize_symbols(Symbols* symbols,int new_size){
Symbols new_symbols;
int i;
init_Symbols(&new_symbols,new_size);
for(i=0;i<(*symbols).length;i++){
new_symbols.indices[i]=(*symbols).indices[i];
new_symbols.expr[i]=(*symbols).expr[i];
}
new_symbols.length=(*symbols).length;
free((*symbols).indices);
free((*symbols).expr);
*symbols=new_symbols;
return(0);
}
// copy
int symbols_cpy(Symbols input, Symbols* output){
init_Symbols(output,input.length);
symbols_cpy_noinit(input,output);
return(0);
}
int symbols_cpy_noinit(Symbols input, Symbols* output){
int i;
if((*output).memory<input.length){
fprintf(stderr,"error: trying to copy a symbols collection of length %d to another with memory %d\n",input.length,(*output).memory);
exit(-1);
}
for(i=0;i<input.length;i++){
(*output).indices[i]=input.indices[i];
polynomial_cpy(input.expr[i],(*output).expr+i);
}
(*output).length=input.length;
return(0);
}
// append an element to a symbols
int symbols_append(int index, Polynomial expr, Symbols* output){
int offset=(*output).length;
if((*output).length>=(*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<input.length;i++){
symbols_append(input.indices[i],input.expr[i],output);
}
return(0);
}
// ------------------ Groups --------------------
// allocate memory
int init_Groups(Groups* groups,int size){
(*groups).indices=calloc(size,sizeof(Int_Array));
(*groups).length=0;
(*groups).memory=size;
return(0);
}
// free memory
int free_Groups(Groups groups){
int i;
for(i=0;i<groups.length;i++){
free_Int_Array(groups.indices[i]);
}
free(groups.indices);
return(0);
}
// resize
int resize_groups(Groups* groups,int new_size){
Groups new_groups;
int i;
init_Groups(&new_groups,new_size);
for(i=0;i<(*groups).length;i++){
new_groups.indices[i]=(*groups).indices[i];
}
new_groups.length=(*groups).length;
free((*groups).indices);
*groups=new_groups;
return(0);
}
// copy
int groups_cpy(Groups input, Groups* output){
init_Groups(output,input.length);
groups_cpy_noinit(input,output);
return(0);
}
int groups_cpy_noinit(Groups input, Groups* output){
int i;
if((*output).memory<input.length){
fprintf(stderr,"error: trying to copy a groups collection of length %d to another with memory %d\n",input.length,(*output).memory);
exit(-1);
}
for(i=0;i<input.length;i++){
int_array_cpy(input.indices[i],(*output).indices+i);
}
(*output).length=input.length;
return(0);
}
// append an element to a groups
int groups_append(Int_Array indices, Groups* output){
int offset=(*output).length;
if((*output).length>=(*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<input.length;i++){
groups_append(input.indices[i],output);
}
return(0);
}
// find which group an index belongs to
int find_group(int index, Groups groups){
int i,j;
for(i=0;i<groups.length;i++){
for(j=0;j<groups.indices[i].length;j++){
if(groups.indices[i].values[j]==index){
return(i);
}
}
}
return(-1);
}

98
src/fields.h Normal file
View File

@ -0,0 +1,98 @@
/*
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.
*/
/* declaring fields */
#ifndef FIELDS_H
#define FIELDS_H
#include "types.h"
// init
int init_Fields_Table(Fields_Table* fields);
int free_Fields_Table(Fields_Table fields);
// determine field type
int field_type(int index, Fields_Table fields);
// check whether a field anticommutes
int is_fermion(int index, Fields_Table fields);
// init
int init_Identities(Identities* identities,int size);
int free_Identities(Identities identities);
// resize
int resize_identities(Identities* identities,int new_size);
// copy
int identities_cpy(Identities input, Identities* output);
int identities_cpy_noinit(Identities input, Identities* output);
// append an element to a identities
int identities_append(Int_Array lhs, Polynomial rhs, Identities* output);
int identities_append_noinit(Int_Array lhs, Polynomial rhs, Identities* output);
// concatenate two identitiess
int identities_concat(Identities input, Identities* output);
// resolve the identities
int resolve_ids(Polynomial* polynomial, Fields_Table fields);
// init
int init_Symbols(Symbols* symbols,int size);
int free_Symbols(Symbols symbols);
// resize
int resize_symbols(Symbols* symbols,int new_size);
// copy
int symbols_cpy(Symbols input, Symbols* output);
int symbols_cpy_noinit(Symbols input, Symbols* output);
// append an element to a symbols
int symbols_append(int index, Polynomial expr, Symbols* output);
int symbols_append_noinit(int index, Polynomial expr, Symbols* output);
// concatenate two symbolss
int symbols_concat(Symbols input, Symbols* output);
// init
int init_Groups(Groups* groups,int size);
int free_Groups(Groups groups);
// resize
int resize_groups(Groups* groups,int new_size);
// copy
int groups_cpy(Groups input, Groups* output);
int groups_cpy_noinit(Groups input, Groups* output);
// append an element to a groups
int groups_append(Int_Array indices, Groups* output);
int groups_append_noinit(Int_Array indices, Groups* output);
// concatenate two groupss
int groups_concat(Groups input, Groups* output);
// find which group an index belongs to
int find_group(int index, Groups groups);
#define FIELDS_H_DONE
#endif

167
src/flow.c Normal file
View File

@ -0,0 +1,167 @@
/*
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 "flow.h"
#include <stdio.h>
#include <stdlib.h>
#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<rccs.length;j++){
print_label(rccs.indices[j], labels);
}
printf("\n\n");
// print initial values
printf("%5d ",0);
for(j=0;j<rccs.length;j++){
printf("% 14.7Le ",rccs.values[j]);
}
printf("\n");
}
for(i=0;i<niter;i++){
// compute a single step
step_flow(&rccs, flow_equation, tol);
// convert ls to alphas
if(display_mode==DISPLAY_NUMERICAL){
// print the result
printf("%5d ",i+1);
for(j=0;j<rccs.length;j++){
printf("% 14.7Le ",rccs.values[j]);
}
printf("\n");
}
}
if(display_mode==DISPLAY_NUMERICAL){
// print labels
printf("\n");
printf("%5s ","n");
for(j=0;j<rccs.length;j++){
print_label(rccs.indices[j], labels);
}
printf("\n\n");
}
if(display_mode==DISPLAY_FINAL){
RCC_print(rccs);
}
return(0);
}
// single step in the flow no exponentials
// inputs: flow_equation, tol
// input/outputs: rccs
int step_flow(RCC* rccs, Grouped_Polynomial flow_equation, long double tol){
int i;
long double* new_rccs=calloc((*rccs).length,sizeof(long double));
Int_Array computed;
init_Int_Array(&computed, (*rccs).length);
// initialize vectors to 0
for(i=0;i<(*rccs).length;i++){
new_rccs[i]=0.;
}
// compute the constants first
for(i=0;i<flow_equation.length;i++){
if(flow_equation.indices[i]<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])<tol){
new_rccs[i]=0.;
}
(*rccs).values[i]=new_rccs[i];
}
}
// for each equation
for(i=0;i<flow_equation.length;i++){
if(flow_equation.indices[i]>=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])<tol){
new_rccs[i]=0.;
}
}
}
// copy results to rccs
for(i=0;i<(*rccs).length;i++){
(*rccs).values[i]=new_rccs[i];
}
// free memory
free_Int_Array(computed);
free(new_rccs);
return(0);
}
// print the label of an rcc (takes constants and derivatives into account)
int print_label(int index, Labels labels){
int i;
int nderivs;
int posin_labels;
Char_Array label;
// constant term
if(index<0){
nderivs=-index/DOFFSET;
for (i=0;i<12-nderivs;i++){
printf(" ");
}
for(i=0;i<nderivs;i++){
printf("d");
}
printf("C%d ",-index-nderivs*DOFFSET);
}
else{
nderivs=index/DOFFSET;
posin_labels=intlist_find_err(labels.indices, labels.length, index-nderivs*DOFFSET);
label=labels.labels[posin_labels];
for (i=0;i<14-label.length-nderivs;i++){
printf(" ");
}
for(i=0;i<nderivs;i++){
printf("d");
}
printf("%s ",char_array_to_str_noinit(labels.labels+posin_labels));
}
return(0);
}

35
src/flow.h Normal file
View File

@ -0,0 +1,35 @@
/*
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 flow numerically
*/
#ifndef NUMERICAL_FLOW_H
#define NUMERICAL_FLOW_H
#include "grouped_polynomial.h"
#include "rcc.h"
// compute flow
int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, int niter, long double tol, int display_mode);
// single step
int step_flow(RCC* rccs, Grouped_Polynomial flow_equation, long double tol);
// print the label of an rcc (takes constants and derivatives into account)
int print_label(int index, Labels labels);
#endif

765
src/grouped_polynomial.c Normal file
View File

@ -0,0 +1,765 @@
/*
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 "grouped_polynomial.h"
#include <stdio.h>
#include <stdlib.h>
#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<gpolynomial.length;i++){
free_Coefficient(gpolynomial.coefs[i]);
}
free(gpolynomial.coefs);
free(gpolynomial.indices);
return(0);
}
// resize the memory allocated to a grouped_polynomial
int resize_grouped_polynomial(Grouped_Polynomial* grouped_polynomial,int new_size){
Grouped_Polynomial new_poly;
int i;
init_Grouped_Polynomial(&new_poly,new_size);
for(i=0;i<(*grouped_polynomial).length;i++){
new_poly.indices[i]=(*grouped_polynomial).indices[i];
new_poly.coefs[i]=(*grouped_polynomial).coefs[i];
}
new_poly.length=(*grouped_polynomial).length;
free((*grouped_polynomial).indices);
free((*grouped_polynomial).coefs);
*grouped_polynomial=new_poly;
return(0);
}
// copy a grouped_polynomial
int grouped_polynomial_cpy(Grouped_Polynomial input, Grouped_Polynomial* output){
init_Grouped_Polynomial(output,input.length);
grouped_polynomial_cpy_noinit(input,output);
return(0);
}
int grouped_polynomial_cpy_noinit(Grouped_Polynomial input, Grouped_Polynomial* output){
int i;
if((*output).memory<input.length){
fprintf(stderr,"error: trying to copy a grouped polynomial of length %d to another with memory %d\n",input.length,(*output).memory);
exit(-1);
}
for(i=0;i<input.length;i++){
(*output).indices[i]=input.indices[i];
coefficient_cpy(input.coefs[i], (*output).coefs+i);
}
(*output).length=input.length;
return(0);
}
// append an element to a grouped_polynomial
int grouped_polynomial_append(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 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;i<input.length;i++){
grouped_polynomial_append(input.indices[i],input.coefs[i],output);
}
return(0);
}
int grouped_polynomial_concat_noinit(Grouped_Polynomial input, Grouped_Polynomial* output){
int i;
for(i=0;i<input.length;i++){
grouped_polynomial_append_noinit(input.indices[i],input.coefs[i],output);
}
// free input arrays
free(input.indices);
free(input.coefs);
return(0);
}
// construct a grouped polynomial from a polynomial, grouping together the terms specified in the id table
// robust algorithm: allows for a term in the polynomial to contribute to several id table entries
int group_polynomial(Polynomial polynomial, Grouped_Polynomial* grouped_polynomial, Id_Table idtable, Fields_Table fields){
int i,j;
// what is left to group
Polynomial remainder;
int index;
Number ratio;
Number num;
Int_Array factor;
int security=0;
int pos;
coef_denom denom;
// init remainder
polynomial_cpy(polynomial, &remainder);
// allocate memory
init_Grouped_Polynomial(grouped_polynomial, idtable.length+1);
// copy indices from idtable and allocate
// the constant term
(*grouped_polynomial).indices[0]=-1;
init_Coefficient((*grouped_polynomial).coefs, COEF_SIZE);
for(i=1;i<=idtable.length;i++){
(*grouped_polynomial).indices[i]=idtable.indices[i-1];
init_Coefficient((*grouped_polynomial).coefs+i, COEF_SIZE);
}
(*grouped_polynomial).length=idtable.length+1;
// keep on going as long as there are terms in the remainder
while(remainder.length>0){
// 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<idtable.length && index==-2;j++){
// loop over terms in the polynomial
for(pos=0;pos<idtable.polynomials[j].length;pos++){
if(int_array_cmp(idtable.polynomials[j].monomials[pos],remainder.monomials[i])==0){
index=j;
break;
}
}
}
}
if(index==-2){
fprintf(stderr,"error: monomial (");
for(j=0;j<polynomial.monomials[i].length;j++){
fprintf(stderr,"%d", polynomial.monomials[i].values[j]);
if(j<polynomial.monomials[i].length-1){
fprintf(stderr,",");
}
}
fprintf(stderr,") not found in idtable\n");
exit(-1);
}
// if not constant
if(index>=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<idtable.polynomials[index].length;j++){
if(j!=pos){
num=number_prod_ret(ratio, idtable.polynomials[index].nums[j]);
number_Qprod_chain(quot(-1,1),&num);
polynomial_append(idtable.polynomials[index].monomials[j], factor, num, &remainder);
free_Number(num);
}
}
free_Int_Array(factor);
free_Number(ratio);
// simplify remainder
polynomial_simplify(&remainder, fields);
}
// constant
else if(index==-1){
// add to coefficient
denom.index=-1;
denom.power=0;
coefficient_append(remainder.factors[i], remainder.nums[i], denom, (*grouped_polynomial).coefs);
// remove from remainder
free_Int_Array(remainder.monomials[i]);
free_Int_Array(remainder.factors[i]);
free_Number(remainder.nums[i]);
remainder.length--;
}
}
// simplify the result
simplify_grouped_polynomial(grouped_polynomial);
free_Polynomial(remainder);
return(0);
}
// construct a grouped polynomial from a polynomial, grouping together the terms specified in the id table.
// identifies sub-polynomials in the polynomial corresponding to the entire rhs of an entry in the id table.
// requires the polynomial and the idtable to be sorted
// can only treat cases in which monomials in different polynomials of the idtable are distinct
int group_polynomial_pickandchoose(Polynomial polynomial, Grouped_Polynomial* grouped_polynomial, Id_Table idtable){
int i,j,k;
// a mask specifying which terms of the polynomial have already been grouped
int* mask=calloc(polynomial.length, sizeof(int));
int index;
Number ratio, ratio_check;
// whether ratio was ever allocated
int alloc_ratio=0;
// whether the correct index was found
int found_index;
int start_index_search;
Int_Array mask_tmp_flips;
coef_denom denom;
// allocate memory
init_Grouped_Polynomial(grouped_polynomial, idtable.length+1);
// copy indices from idtable and allocate
// the constant term
(*grouped_polynomial).indices[0]=-1;
init_Coefficient((*grouped_polynomial).coefs, COEF_SIZE);
for(i=1;i<=idtable.length;i++){
(*grouped_polynomial).indices[i]=idtable.indices[i-1];
init_Coefficient((*grouped_polynomial).coefs+i, COEF_SIZE);
}
(*grouped_polynomial).length=idtable.length+1;
// loop over monomials
for(i=0;i<polynomial.length;i++){
// check that the term hasn't already been added
if(mask[i]==0){
// loop until the correct index is found (the polynomial must contain all the terms in the index and the numerical factors must match)
found_index=0;
start_index_search=0;
while(found_index==0){
found_index=1;
// find entry
index=find_id(polynomial.monomials[i], idtable,start_index_search);
// easier to debug if the error is here instead of inside find_id
if(index==-2){
fprintf(stderr,"error: monomial not found in idtable\n");
exit(-1);
}
// if not constant
if(index>=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<idtable.polynomials[index].length && found_index==1;j++){
// find the monomial in the polynomial
for(k=i;k<polynomial.length;k++){
// only check if mask==0
// only check if the factors are correct
if(mask[k]==0 && int_array_cmp(polynomial.factors[i],polynomial.factors[k])==0 && int_array_cmp(idtable.polynomials[index].monomials[j],polynomial.monomials[k])==0){
ratio_check=number_quot_ret(polynomial.nums[k],idtable.polynomials[index].nums[j]);
// if the factors don't factor
if(alloc_ratio!=0 && number_compare(ratio,ratio_check)==0){
found_index=0;
break;
}
// check that ratio was allocated
if(alloc_ratio!=0){
free_Number(ratio);
}
ratio=ratio_check;
alloc_ratio=1;
// added to polynomial
mask[k]=1;
// keep track of the flips so that they can be undone if the index turns out to be incorrect
int_array_append(k,&mask_tmp_flips);
break;
}
}
// error if the monomial could not be found
if(k==polynomial.length){
found_index=0;
}
}
// if the index was incorrect
if(found_index==0){
// reset mask
for(j=0;j<mask_tmp_flips.length;j++){
mask[mask_tmp_flips.values[j]]=0;
}
// start index search at next item
start_index_search=index+1;
}
else{
// add to grouped polynomial
denom.index=-1;
denom.power=1;
coefficient_append(polynomial.factors[i], ratio, denom, (*grouped_polynomial).coefs+index+1);
}
if(alloc_ratio==1){
free_Number(ratio);
alloc_ratio=0;
}
free_Int_Array(mask_tmp_flips);
}
// constant
else if(index==-1){
mask[i]=1;
denom.index=-1;
denom.power=0;
coefficient_append(polynomial.factors[i], polynomial.nums[i], denom, (*grouped_polynomial).coefs);
}
}
}
}
// check all the terms were grouped
for(i=0;i<polynomial.length;i++){
if(mask[i]==0){
fprintf(stderr,"error: this polynomial could not be grouped: no matches were found for some of the terms\n");
exit(-1);
}
}
free(mask);
return(0);
}
// find the entry in the idtable containing monomial
// start search at the specified index
int find_id(Int_Array monomial, Id_Table idtable, int start){
int i,j;
// constant
if(monomial.length==0){
return(-1);
}
// loop over entries
for(i=start;i<idtable.length;i++){
// loop over terms in the polynomial
for(j=0;j<idtable.polynomials[i].length;j++){
if(int_array_cmp(idtable.polynomials[i].monomials[j],monomial)==0){
return(i);
}
}
}
return(-2);
}
// simplify grouped polynomial
int simplify_grouped_polynomial(Grouped_Polynomial* polynomial){
int i;
for(i=0;i<(*polynomial).length;i++){
coefficient_simplify((*polynomial).coefs+i);
}
return(0);
}
// derive a flow equation with respect to an unknown variable
// equivalent to DB.dl where dl are symbols for the derivatives of the indices in the flow equation with respect to the unknown variable
// indices specifies the list of indices that depend on the variable
int flow_equation_derivx(Grouped_Polynomial flow_equation, Int_Array indices, Grouped_Polynomial* dflow){
int i,j,k;
Coefficient tmp_coef;
// alloc
init_Grouped_Polynomial(dflow, flow_equation.length);
// for each equation
for(i=0;i<flow_equation.length;i++){
// copy indices
if(flow_equation.indices[i]>=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<indices.length;j++){
coefficient_deriv(flow_equation.coefs[i], indices.values[j], &tmp_coef);
// multiply each coefficient by the appropriate dl[j]
for(k=0;k<tmp_coef.length;k++){
// only in non-trivial cases
if(number_is_zero(tmp_coef.nums[k])==0){
// non-constants
if(indices.values[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;k<flow_equation.length;k++){
init_Coefficient((*output).coefs+k, COEF_SIZE);
// loop over monomials
for(i=0;i<flow_equation.coefs[k].length;i++){
// derivative of the numerator
monomial_deriv(flow_equation.coefs[k].factors[i], index, &factor, &match_count);
// if the derivative doesn't vanish, add it to the coefficient
if(match_count>0){
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<grouped_polynomial.length;i++){
//print lhs
// negative indices are constants
if(grouped_polynomial.indices[i]<0){
// count derivatives
dcount=-grouped_polynomial.indices[i]/DOFFSET;
for(j=0;j<3-dcount;j++){
printf(" ");
}
printf("[");
for(j=0;j<dcount;j++){
printf("d");
}
printf("C%d] =",-grouped_polynomial.indices[i]-dcount*DOFFSET);
}
else{
// count derivatives
dcount=grouped_polynomial.indices[i]/DOFFSET;
for(j=0;j<2-dcount;j++){
printf(" ");
}
printf("[");
for(j=0;j<dcount;j++){
printf("d");
}
printf("%c%2d] =",lhs_pre,grouped_polynomial.indices[i]-dcount*DOFFSET);
}
// rhs
init_Char_Array(&buffer, STR_SIZE);
coefficient_sprint(grouped_polynomial.coefs[i],&buffer,9,rhs_pre);
if(buffer.length>0){
printf("%s",buffer.str);
}
free_Char_Array(buffer);
if(i<grouped_polynomial.length-1){
printf(",");
}
// extra \n
printf("\n");
}
return(0);
}
// read from string
#define PP_NULL_MODE 0
#define PP_COEF_MODE 1
#define PP_INDEX_MODE 3
#define PP_COMMENT_MODE 4
#define PP_BRACKET_MODE 5
#define PP_CONSTANT_MODE 6
int char_array_to_Grouped_Polynomial(Char_Array str, Grouped_Polynomial* output){
// buffer
char* buffer=calloc(str.length+1,sizeof(char));
char* buffer_ptr=buffer;
int index=-2;
Coefficient coef;
int i,j;
int mode;
int dcount=0;
init_Grouped_Polynomial(output, EQUATION_SIZE);
// loop over input
mode=PP_NULL_MODE;
for(j=0;j<str.length;j++){
if(mode==PP_COMMENT_MODE){
if(str.str[j]=='\n'){
mode=PP_NULL_MODE;
}
}
// stay in polynomial mode until ','
else if(mode==PP_COEF_MODE){
if(str.str[j]==','){
// parse polynomial
str_to_Coefficient(buffer, &coef);
// write index and polynomial
grouped_polynomial_append_noinit(index, coef, output);
mode=PP_NULL_MODE;
}
else{
buffer_ptr=str_addchar(buffer_ptr,str.str[j]);
}
}
else{
switch(str.str[j]){
// index
case '[':
if(mode==PP_NULL_MODE){
mode=PP_BRACKET_MODE;
buffer_ptr=buffer;
*buffer_ptr='\0';
// reset derivatives count
dcount=0;
}
break;
case '%':
if(mode==PP_BRACKET_MODE){
mode=PP_INDEX_MODE;
}
break;
// constant term
case 'C':
if(mode==PP_BRACKET_MODE){
mode=PP_CONSTANT_MODE;
}
break;
// derivatives
case 'd':
if(mode==PP_BRACKET_MODE || mode==PP_INDEX_MODE || mode==PP_CONSTANT_MODE){
dcount++;
}
break;
// write index
case ']':
sscanf(buffer,"%d",&i);
if(mode==PP_INDEX_MODE){
index=i+dcount*DOFFSET;
}
else if(mode==PP_CONSTANT_MODE){
index=-i-dcount*DOFFSET;
}
mode=PP_NULL_MODE;
break;
// coef mode
case '=':
if(mode==PP_NULL_MODE){
buffer_ptr=buffer;
*buffer_ptr='\0';
mode=PP_COEF_MODE;
}
break;
// comment
case '#':
mode=PP_COMMENT_MODE;
break;
default:
if(mode!=PP_NULL_MODE){
buffer_ptr=str_addchar(buffer_ptr,str.str[j]);
}
break;
}
}
}
// last step
if(mode==PP_COEF_MODE){
str_to_Coefficient(buffer, &coef);
grouped_polynomial_append_noinit(index, coef, output);
}
free(buffer);
return(0);
}
// evaluate an equation on a vector
int evaleq(RCC* rccs, Grouped_Polynomial poly){
int i;
long double* res=calloc((*rccs).length,sizeof(long double));
if((*rccs).length!=poly.length){
fprintf(stderr, "error: trying to evaluate an flow equation with %d components on an rcc with %d\n",poly.length,(*rccs).length);
exit(-1);
}
// initialize vectors to 0
for(i=0;i<(*rccs).length;i++){
res[i]=0.;
}
// for each equation
for(i=0;i<poly.length;i++){
evalcoef(*rccs, poly.coefs[i], res+i);
}
// copy res to rccs
for(i=0;i<(*rccs).length;i++){
(*rccs).values[i]=res[i];
}
// free memory
free(res);
return(0);
}

74
src/grouped_polynomial.h Normal file
View File

@ -0,0 +1,74 @@
/*
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.
*/
/*
Polynomials can be represented in a grouped way, by putting the terms
proportional to the same monomial together.
The main part of a grouped polynomial is the list of coefficients,
each of which stands for the set of terms proportional to the
corresponding monomial, specified in indices and labels.
*/
#ifndef GROUPED_POLYNOMIAL_H
#define GROUPED_POLYNOMIAL_H
#include "types.h"
// memory
int init_Grouped_Polynomial(Grouped_Polynomial* gpolynomial, int size);
int free_Grouped_Polynomial(Grouped_Polynomial gpolynomial);
// resize the memory allocated to a grouped polynomial
int resize_grouped_polynomial(Grouped_Polynomial* grouped_polynomial,int new_size);
// copy a grouped polynomial
int grouped_polynomial_cpy(Grouped_Polynomial input, Grouped_Polynomial* output);
int grouped_polynomial_cpy_noinit(Grouped_Polynomial input, Grouped_Polynomial* output);
// append an element to a polynomial
int grouped_polynomial_append(int index, Coefficient coef, Grouped_Polynomial* output);
int grouped_polynomial_append_noinit(int index, Coefficient coef, Grouped_Polynomial* output);
// concatenate two polynomials
int grouped_polynomial_concat(Grouped_Polynomial input, Grouped_Polynomial* output);
int grouped_polynomial_concat_noinit(Grouped_Polynomial input, Grouped_Polynomial* output);
// construct a grouped polynomial from a polynomial, grouping together the terms specified in the id table.
int group_polynomial(Polynomial polynomial, Grouped_Polynomial* grouped_polynomial, Id_Table idtable, Fields_Table fields);
// more naive and faster version in which the terms in polynomial corresponding to a polynomial in the id table are grouped together (does not allow parts of terms to be grouped together)
int group_polynomial_pickandchoose(Polynomial polynomial, Grouped_Polynomial* grouped_polynomial, Id_Table idtable);
// find the entry in the idtable containing monomial
int find_id(Int_Array monomial, Id_Table idtable, int start);
// simplify grouped polynomial
int simplify_grouped_polynomial(Grouped_Polynomial* polynomial);
// derive a flow equation with respect to an unknown variable
int flow_equation_derivx(Grouped_Polynomial flow_equation, Int_Array indices, Grouped_Polynomial* dflow);
// print a grouped polynomial
int grouped_polynomial_print(Grouped_Polynomial grouped_polynomial, char lhs_pre, char rhs_pre);
// read from string
int char_array_to_Grouped_Polynomial(Char_Array str, Grouped_Polynomial* output);
// evaluate an equation on an RCC
int evaleq(RCC* rccs, Grouped_Polynomial poly);
#endif

122
src/idtable.c Normal file
View File

@ -0,0 +1,122 @@
/*
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 "idtable.h"
#include <stdio.h>
#include <stdlib.h>
#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<idtable.length;i++){
free_Polynomial(idtable.polynomials[i]);
}
free(idtable.indices);
free(idtable.polynomials);
return(0);
}
// resize the memory allocated to a idtable
int resize_idtable(Id_Table* idtable,int new_size){
Id_Table new_idtable;
int i;
init_Id_Table(&new_idtable,new_size);
for(i=0;i<(*idtable).length;i++){
new_idtable.indices[i]=(*idtable).indices[i];
new_idtable.polynomials[i]=(*idtable).polynomials[i];
}
new_idtable.length=(*idtable).length;
free((*idtable).indices);
free((*idtable).polynomials);
*idtable=new_idtable;
return(0);
}
// copy an idtable
int idtable_cpy(Id_Table input, Id_Table* output){
init_Id_Table(output,input.length);
idtable_cpy_noinit(input,output);
return(0);
}
int idtable_cpy_noinit(Id_Table input, Id_Table* output){
int i;
if((*output).memory<input.length){
fprintf(stderr,"error: trying to copy an idtable of length %d to another with memory %d\n",input.length,(*output).memory);
exit(-1);
}
for(i=0;i<input.length;i++){
(*output).indices[i]=input.indices[i];
polynomial_cpy(input.polynomials[i],(*output).polynomials+i);
}
(*output).length=input.length;
return(0);
}
// append an element to a idtable
int idtable_append(int index, Polynomial polynomial, Id_Table* output){
int offset=(*output).length;
if((*output).length>=(*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<input.length;i++){
idtable_append(input.indices[i],input.polynomials[i],output);
}
return(0);
}

44
src/idtable.h Normal file
View File

@ -0,0 +1,44 @@
/*
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.
*/
/* correspondence table to group a polynomial into a flow equation */
#ifndef IDTABLE_H
#define IDTABLE_H
#include "types.h"
// allocate memory
int init_Id_Table(Id_Table* idtable,int size);
// free memory
int free_Id_Table(Id_Table idtable);
// resize the memory allocated to a idtable
int resize_idtable(Id_Table* idtable,int new_size);
// copy an idtable
int idtable_cpy(Id_Table input, Id_Table* output);
int idtable_cpy_noinit(Id_Table input, Id_Table* output);
// append an element to a idtable
int idtable_append(int index, Polynomial polynomial, Id_Table* output);
int idtable_append_noinit(int index, Polynomial polynomial, Id_Table* output);
// concatenate two idtables
int idtable_concat(Id_Table input, Id_Table* output);
#define IDTABLE_H_DONE
#endif

99
src/istring.c Normal file
View File

@ -0,0 +1,99 @@
/*
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 <stdlib.h>
#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);
}

40
src/istring.h Normal file
View File

@ -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

1449
src/kondo.c Normal file
View File

File diff suppressed because it is too large Load Diff

76
src/kondo.h Normal file
View File

@ -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

126
src/kondo_preprocess.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#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<str_args.length;i++){
printf("%s\n",str_args.strs[i].str);
if(i<str_args.length-1){
printf("&\n");
}
}
//free memory
free_Str_Array(str_args);
return(0);
}
// read cli arguments
#define CP_FLAG_DIMENSION 1
int read_args_kondo_pp(int argc,const char* argv[], Str_Array* str_args, Kondopp_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_kondo_pp();
exit(-1);
}
// defaults
// dimensions (including time)
(*opts).dimension=2;
// loop over arguments
for(i=1;i<argc;i++){
// flag
if(argv[i][0]=='-'){
for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){
switch(*ptr){
case 'd':
flag=CP_FLAG_DIMENSION;
break;
// print version
case 'v':
printf("kondo_preprocess " VERSION "\n");
exit(1);
break;
}
}
}
// number of dimensions
else if (flag==CP_FLAG_DIMENSION){
sscanf(argv[i],"%d",&((*opts).dimension));
flag=0;
}
// read file name from command-line
else{
file=argv[i];
exists_file=1;
}
}
read_config_file(str_args, file, 1-exists_file);
return(0);
}
// print usage message
int print_usage_kondo_pp(){
printf("\nusage:\n kondo_preprocess [-d dimension] <filename>\n\n");
return(0);
}

793
src/mean.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#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<n;i++){
indices_plus[i]=intlist_find_err(propagator.indices, propagator.length, internal_plus.values[i]);
indices_minus[i]=intlist_find_err(propagator.indices, propagator.length, -internal_minus.values[i]);
}
// init pairing and mask
exists_next=init_pairing(pairing, mask, n, propagator, indices_plus, indices_minus)+1;
// initial sign
pairing_sign=permutation_signature(pairing,n);
// mixing sign (from ordering psi+psi-): (-1)^{n(n+1)/2}
if((n*(n+1))/2 %2 ==0){
mixing_sign=1;
}
else{
mixing_sign=-1;
}
// loop over pairings
// loop ends when the first pairing leaves the array
while(exists_next==1){
num=polynomial_one();
// propagator product for the current pairing (only simplify after all pairings)
for(i=0;i<n;i++){
polynomial_prod_chain_nosimplify(propagator.matrix[indices_plus[i]][indices_minus[pairing[i]]],&num, fields);
}
polynomial_multiply_Qscalar(num,quot(mixing_sign*pairing_sign,1));
polynomial_concat_noinit_inplace(num,&num_summed);
exists_next=next_pairing(pairing, mask, n, propagator, indices_plus, indices_minus)+1;
pairing_sign=permutation_signature(pairing,n);
}
// only simplify in mean_symbols
polynomial_prod_chain_nosimplify(num_summed,out,fields);
free_Polynomial(num_summed);
free(pairing);
free(mask);
free(indices_plus);
free(indices_minus);
return(0);
}
// first pairing with a non-vanishing propagator
int init_pairing(int* pairing, int* mask, int n, Polynomial_Matrix propagator, int* indices_plus, int* indices_minus){
// index we want to increment
int move=0;
int i;
for(i=0;i<n;i++){
pairing[i]=-1;
mask[i]=0;
}
// loop until move is out of range
while(move>=0 && move<n){
// move
pairing[move]=next_wick(move, pairing[move], mask, n, propagator, indices_plus, indices_minus);
// if the next term does not exist, then move previous index
if(pairing[move]==-1){
move--;
}
// else move next index
else{
move++;
}
}
// if move==-1, then there is no first term, return -1
if(move==-1){
return(-1);
}
// if the first term exists
return(0);
}
// next pairing with a non-vanishing propagator
int next_pairing(int* pairing, int* mask, int n, Polynomial_Matrix propagator, int* indices_plus, int* indices_minus){
// index we want to increment
int move=n-1;
// last index
mask[pairing[n-1]]=0;
// loop until move is out of range
while(move>=0 && move<n){
// move
pairing[move]=next_wick(move, pairing[move], mask, n, propagator, indices_plus, indices_minus);
// if the next term does not exist, then move previous index
if(pairing[move]==-1){
move--;
}
// else move next index
else{
move++;
}
}
// if move==-1, then there is no next term, return -1
if(move==-1){
return(-1);
}
// if the next term exists
return(0);
}
// next term in the Wick expansion
int next_wick(int index, int start, int* mask, int n, Polynomial_Matrix propagator, int* indices_plus, int* indices_minus){
int i;
// unset mask
if(start>=0 && start<n){
mask[start]=0;
}
// find next position
for(i=start+1;i<n;i++){
// if the propagator doesn't vanish
if(mask[i]==0 && polynomial_is_zero(propagator.matrix[indices_plus[index]][indices_minus[i]])==0){
mask[i]=1;
return(i);
}
}
// no next term
return(-1);
}
/* Older function: propagator as number
// compute the mean of a monomial of internal fields (with split + and -)
// compute all contractions
int mean_internal_slow(Int_Array internal_plus, Int_Array internal_minus, Number* outnum, Polynomial_Matrix propagator){
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;
Number num;
Number num_summed=number_zero();
// propagator matrix indices
int index1, index2;
int i,j,k,l;
// init pairing and mask
for(i=0;i<n-1;i++){
pairing[i]=i;
mask[i]=1;
}
pairing[n-1]=n-1;
pairing_sign=1;
// mixing sign: (-1)^{n(n+1)/2}
if((n*(n+1))/2 %2 ==0){
mixing_sign=1;
}
else{
mixing_sign=-1;
}
// loop over pairings
// loop ends when the first pairing leaves the array
while(pairing[0]<n){
num=number_one();
// propagator product for the current pairing
for(i=0;i<n;i++){
// indices within the propagator matrix
index1=intlist_find_err(propagator.indices, propagator.length, internal_plus.values[i]);
index2=intlist_find_err(propagator.indices, propagator.length, -internal_minus.values[pairing[i]]);
number_prod_chain(propagator.matrix[index1][index2],&num);
}
number_Qprod_chain(quot(mixing_sign*pairing_sign,1),&num);
number_add_chain(num,&num_summed);
free_Number(num);
// next pairing
// last element of the pairing that we can move
for(i=n-1;i>=0;i--){
// move i-th
mask[pairing[i]]=0;
// search for next possible position
for(j=pairing[i]+1;j<n;j++){
if(mask[j]==0){
break;
}
}
// if the next position exists
if(j<n){
// sign correction: change sign by (-1)^{1+(n-i)(n-i-1)/2}
// actually (-1)^{1+(n-1-i)(n-1-i-1)/2
if(((n-i-1)*(n-i-2))/2 % 2==0){
pairing_sign*=-1;
}
pairing[i]=j;
mask[j]=1;
// put the remaining pairings at their left-most possible values
if(i<n-1){
k=i+1;
for(l=0;l<n;l++){
if(mask[l]==0){
mask[l]=1;
pairing[k]=l;
k++;
// if exhausted all indices
if(k>=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<monomial.length;i++){
if(int_array_find(abs(monomial.values[i]),fields.internal)>=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<monomial.length;i++){
if(field_type(monomial.values[i], fields)==FIELD_SYMBOL){
int_array_append(intlist_find_err(fields.symbols.indices, fields.symbols.length, monomial.values[i]), &symbol_list);
}
else{
int_array_append(monomial.values[i], &base_monomial);
}
}
power=symbol_list.length;
// trivial case
if(power==0){
mean(monomial, &mean_num, fields, propagator);
polynomial_concat_noinit(mean_num, output);
free_Int_Array(symbol_list);
free_Int_Array(base_monomial);
return(0);
}
else{
// initialize current term to a position that has no repetitions
current_term=calloc(power,sizeof(int));
exists_next=init_prod(current_term, symbol_list, fields, power, base_monomial)+1;
}
// loop over terms; the loop stops when all the pointers are at the end of the first symbol
while(exists_next==1){
// construct monomial
int_array_cpy(base_monomial, &tmp_monomial);
tmp_num=number_one();
for(i=0;i<power;i++){
int_array_concat(fields.symbols.expr[symbol_list.values[i]].monomials[current_term[i]], &tmp_monomial);
number_prod_chain(fields.symbols.expr[symbol_list.values[i]].nums[current_term[i]], &tmp_num);
}
// check whether the monomial vanishes
if(check_monomial_match(tmp_monomial, fields)==1){
// sort monomial
sign=1;
monomial_sort(tmp_monomial, 0, tmp_monomial.length-1, fields, &sign);
number_Qprod_chain(quot(sign,1), &tmp_num);
// mean
mean_groups(tmp_monomial, &mean_poly, fields, propagator, groups, computed);
// write to output
polynomial_multiply_scalar(mean_poly,tmp_num);
polynomial_concat_noinit_inplace(mean_poly, output);
}
free_Number(tmp_num);
free_Int_Array(tmp_monomial);
// next term
exists_next=next_prod(current_term, symbol_list, fields, power, base_monomial)+1;
// simplfiy every 25 steps (improves both memory usage and performance)
if(simplify_freq %25 ==0){
polynomial_simplify(output, fields);
simplify_freq=0;
}
simplify_freq++;
}
// simplify
polynomial_simplify(output, fields);
// write computed
identities_append(monomial, *output, computed);
// free memory
free(current_term);
free_Int_Array(symbol_list);
free_Int_Array(base_monomial);
return(0);
}
// first term in product with no repetitions
int init_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int power, Int_Array base_monomial){
// index we want to increment
int move=0;
// tmp monomial
Int_Array monomial;
int i;
init_Int_Array(&monomial, base_monomial.length+5*power);
int_array_cpy_noinit(base_monomial, &monomial);
// init current term
for(i=0;i<power;i++){
current_term[i]=-1;
}
// loop until move is out of range
while(move>=0 && move<power){
// move
current_term[move]=next_term_norepeat(current_term[move], fields.symbols.expr[symbol_list.values[move]], &monomial, fields);
// if the next term does not exist, then move previous index
if(current_term[move]==-1){
move--;
}
// else move next index
else{
move++;
}
}
free_Int_Array(monomial);
// if move==-1, then there is no first term, return -1
if(move==-1){
return(-1);
}
// if the next term exists
return(0);
}
// next term in product with no repetitions
int next_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int power, Int_Array base_monomial){
// index we want to increment
int move=power-1;
// tmp monomial
Int_Array monomial;
int i;
// init monomial
init_Int_Array(&monomial, base_monomial.length+5*power);
int_array_cpy_noinit(base_monomial, &monomial);
for(i=0;i<=move;i++){
int_array_concat(fields.symbols.expr[symbol_list.values[i]].monomials[current_term[i]],&monomial);
}
// loop until move is out of range
while(move>=0 && move<power){
// move
current_term[move]=next_term_norepeat(current_term[move], fields.symbols.expr[symbol_list.values[move]], &monomial, fields);
// if the next term does not exist, then move previous index
if(current_term[move]==-1){
move--;
}
// else move next index
else{
move++;
}
}
free_Int_Array(monomial);
// if move==-1, then there is no next term, return -1
if(move==-1){
return(-1);
}
// if the next term exists
return(0);
}
// find the next term in a polynomial that can be multiplied to a given monomial and add it to the monomial
int next_term_norepeat(int start, Polynomial polynomial, Int_Array* monomial, Fields_Table fields){
int i;
// remove last term from monomial
if(start>=0 && start<polynomial.length){
(*monomial).length-=polynomial.monomials[start].length;
}
// find next position
for(i=start+1;i<polynomial.length;i++){
// if no repetitions
if(check_monomial_addterm(*monomial,polynomial.monomials[i],fields)==1){
// append to monomial
int_array_concat(polynomial.monomials[i], monomial);
return(i);
}
}
// no next term
return(-1);
}
// signature of a permutation
int permutation_signature(int* permutation, int n){
int* tmp_array=calloc(n,sizeof(int));
int i;
int ret=1;
for(i=0;i<n;i++){
tmp_array[i]=permutation[i];
}
sort_fermions(tmp_array, 0, n-1, &ret);
free(tmp_array);
return(ret);
}
// sort a list of anti-commuting variables
int sort_fermions(int* array, int begin, int end, int* sign){
int i;
int tmp;
int index;
// the pivot: middle of the monomial
int pivot=(begin+end)/2;
// if the monomial is non trivial
if(begin<end){
// send pivot to the end
if(pivot!=end){
tmp=array[end];
array[end]=array[pivot];
array[pivot]=tmp;
*sign*=-1;
}
// loop over the others
for(i=begin, index=begin;i<end;i++){
// compare with pivot
if(array[i]<array[end]){
// if smaller, exchange with reference index
if(i!=index){
tmp=array[i];
array[i]=array[index];
array[index]=tmp;
*sign*=-1;
}
// move reference index
index++;
}
}
// put pivot (which we had sent to the end) in the right place
if(end!=index){
tmp=array[end];
array[end]=array[index];
array[index]=tmp;
*sign*=-1;
}
// recurse
sort_fermions(array, begin, index-1, sign);
sort_fermions(array, index+1, end, sign);
}
return(0);
}
// mean while factoring groups out
int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed){
Polynomial num_mean;
Int_Array tmp_monomial;
int i;
int group=-2;
int next_group=-2;
Polynomial tmp_poly;
int sign;
init_Polynomial(output, MONOMIAL_SIZE);
// check whether there are symbols
// requires the symbols to be at the end of the monomial
if(monomial.length==0 || field_type(monomial.values[monomial.length-1], fields)!=FIELD_SYMBOL){
// mean
mean(monomial, &num_mean, fields, propagator);
// add to output
polynomial_concat_noinit(num_mean,output);
}
else{
// sort into groups
if(groups.length>0){
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(i<monomial.length){
next_group=find_group(monomial.values[i], groups);
}
// if group changes, take mean
if((i>0 && 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<monomial.length){
int_array_append(monomial.values[i], &tmp_monomial);
}
group=next_group;
}
// sign correction
if(sign==-1){
polynomial_multiply_Qscalar(*output,quot(sign,1));
}
free_Int_Array(tmp_monomial);
}
return(0);
}
// mean of a polynomial
// argument struct for multithreaded mean
struct mean_args{
Polynomial* polynomial;
Fields_Table fields;
Polynomial_Matrix propagator;
Groups groups;
};
// multithreaded
int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads){
int i;
Polynomial thread_polys[threads];
pthread_t thread_ids[threads];
struct mean_args args[threads];
int len=(*polynomial).length;
// alloc
for(i=0;i<threads;i++){
init_Polynomial(thread_polys+i,(*polynomial).length/threads+1);
// arguments
args[i].fields=fields;
args[i].propagator=propagator;
args[i].groups=groups;
}
// split polynomial
for(i=0;i<len;i++){
polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], thread_polys+(i % threads));
}
// start threads
for(i=0;i<threads;i++){
args[i].polynomial=thread_polys+i;
pthread_create(thread_ids+i, NULL, polynomial_mean_thread, (void*)(args+i));
}
free_Polynomial(*polynomial);
init_Polynomial(polynomial, len);
// wait for completion and join
for(i=0;i<threads;i++){
pthread_join(thread_ids[i], NULL);
polynomial_concat_noinit(thread_polys[i], polynomial);
}
polynomial_simplify(polynomial, fields);
return(0);
}
// mean for one of the threads
void* polynomial_mean_thread(void* mean_args){
struct mean_args *args=mean_args;
polynomial_mean((*args).polynomial,(*args).fields,(*args).propagator,(*args).groups);
return(NULL);
}
// single threaded version
int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups){
int i,j;
Polynomial output;
Polynomial tmp_poly;
// a list of already computed means
Identities computed;
init_Polynomial(&output, (*polynomial).length);
init_Identities(&computed, EQUATION_SIZE);
remove_unmatched_plusminus(polynomial, fields);
// mean of each monomial
for(i=0;i<(*polynomial).length;i++){
fprintf(stderr,"computing %d of %d means\n",i,(*polynomial).length-1);
mean_groups((*polynomial).monomials[i], &tmp_poly, fields, propagator, groups, &computed);
// write factors
for(j=0;j<tmp_poly.length;j++){
int_array_concat((*polynomial).factors[i], tmp_poly.factors+j);
number_prod_chain((*polynomial).nums[i], tmp_poly.nums+j);
}
// add to output
polynomial_concat_noinit(tmp_poly, &output);
// simplify (simplify here in order to keep memory usage low)
polynomial_simplify(&output, fields);
}
free_Identities(computed);
free_Polynomial(*polynomial);
*polynomial=output;
return(0);
}

70
src/mean.h Normal file
View File

@ -0,0 +1,70 @@
/*
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 the mean of a monomial or a polynomial
*/
#ifndef MEAN_H
#define MEAN_H
#include "types.h"
// mean of a monomial
int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Matrix propagator);
// 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);
// first pairing with a non-vanishing propagator
int init_pairing(int* pairing, int* mask, int n, Polynomial_Matrix propagator, int* indices_plus, int* indices_minus);
// next pairing with a non-vanishing propagator
int next_pairing(int* pairing, int* mask, int n, Polynomial_Matrix propagator, int* indices_plus, int* indices_minus);
// next term in the Wick expansion
int next_wick(int index, int start, int* mask, int n, Polynomial_Matrix propagator, int* indices_plus, int* indices_minus);
/*
int mean_internal_slow(Int_Array internal_plus, Int_Array internal_minus, Number* outnum, Polynomial_Matrix propagator);
*/
// 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);
// compute the mean of a monomial containing symbolic expressions
int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed);
// first term in product with no repetitions
int init_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int power, Int_Array base_monomial);
// next term in product with no repetitions
int next_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int power, Int_Array base_monomial);
// find the next term in a polynomial that can be multiplied to a given monomial and add it to the monomial
int next_term_norepeat(int start, Polynomial polynomial, Int_Array* monomial, Fields_Table fields);
// signature of a permutation
int permutation_signature(int* permutation, int n);
// sort a list of anti-commuting variables
int sort_fermions(int* array, int begin, int end, int* sign);
// mean while factoring groups out
int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed);
// compute the mean of a polynomial
int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups);
// multithreaded
int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads);
// single thread mean
void* polynomial_mean_thread(void* mean_args);
#endif

301
src/meankondo.c Normal file
View File

@ -0,0 +1,301 @@
/*
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.
*/
/*
meankondo
A simple tool to compute the renormalization group flow for Fermionic hierarchical models
*/
#include <stdio.h>
#include <stdlib.h>
// 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<argc;i++){
// flag
if(argv[i][0]=='-'){
for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){
switch(*ptr){
// threads
case 't':
flag=CP_FLAG_THREADS;
break;
// chain
case 'C':
(*opts).chain=1;
break;
// print version
case 'v':
printf("meankondo " VERSION "\n");
exit(1);
break;
default:
print_usage_meankondo();
exit(-1);
break;
}
}
}
// threads
else if(flag==CP_FLAG_THREADS){
sscanf(argv[i],"%d",&((*opts).threads));
flag=0;
}
// read file name from command-line
else{
file=argv[i];
exists_file=1;
}
}
read_config_file(str_args, file, 1-exists_file);
return(0);
}
// print usage message
int print_usage_meankondo(){
printf("\nusage:\n meankondo [-t threads] [-C] <filename>\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;i<str_args.length;i++){
// check whether to print the str_arg
get_str_arg_title(str_args.strs[i], &arg_header);
if (\
str_cmp(arg_header.str, "symbols")==0 &&\
str_cmp(arg_header.str, "groups")==0 &&\
str_cmp(arg_header.str, "fields")==0 &&\
str_cmp(arg_header.str, "identities")==0 &&\
str_cmp(arg_header.str, "propagator")==0 &&\
str_cmp(arg_header.str, "input_polynomial")==0 &&\
str_cmp(arg_header.str, "id_table")==0 ){
printf("%s\n&\n",str_args.strs[i].str);
}
free_Char_Array(arg_header);
}
// print flow equation
printf("#!flow_equation\n");
}
// print flow equation
grouped_polynomial_print(flow_equation,'%','%');
// free memory
free_Id_Table(idtable);
free_Grouped_Polynomial(flow_equation);
return(0);
}
// 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){
// expectation
Polynomial exp_poly;
polynomial_cpy(init_poly,&exp_poly);
// average
if(threads>1){
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);
}

116
src/meantools.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
// 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 <filename>\n meantools derive [-d derivatives] -V <variables> <filename>\n meantools eval -R <rccs> <filename>\n\n");
return(0);
}

195
src/meantools_deriv.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#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<argc;i++){
// flag
if(argv[i][0]=='-'){
for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){
switch(*ptr){
// number of derivatives
case 'd':
flag=CP_FLAG_DERIVS;
break;
case 'V':
flag=CP_FLAG_VARS;
break;
}
}
}
// number of derivatives
else if(flag==CP_FLAG_DERIVS){
sscanf(argv[i],"%d",&((*opts).deriv_derivs));
flag=0;
}
// variables
else if(flag==CP_FLAG_VARS){
// if the argument is "all" then derive wrt all variables
if(str_cmp((char*)argv[i],"all")){
(*opts).deriv_vars.length=-2;
}
else{
str_to_char_array((char*)argv[i], &buffer);
int_array_read(buffer,&((*opts).deriv_vars));
free_Char_Array(buffer);
}
flag=0;
}
// read file name from command-line
else{
file=argv[i];
exists_file=1;
}
}
read_config_file(str_args, file, 1-exists_file);
return(0);
}
// derive a flow equation
int tool_deriv(Str_Array str_args, Meantools_Options opts){
// index of the entry in the input file
int arg_index;
// flow equation
Grouped_Polynomial flow_equation;
// flow equation for the derivatives
Grouped_Polynomial flow_equation_deriv;
int i;
// parse flow equation
// if there is a unique argument, assume it is the flow equation
if(str_args.length==1){
char_array_to_Grouped_Polynomial(str_args.strs[0], &flow_equation);
}
else{
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);
}
// variables
// check they were not specified on the command line
if(opts.deriv_vars.length==-1){
arg_index=find_str_arg("variables", str_args);
if(arg_index>=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<flow_equation.length;i++){
int_array_append(flow_equation.indices[i], &(opts.deriv_vars));
}
}
// compute derivatives
flow_equation_derivative(opts.deriv_derivs, opts.deriv_vars, flow_equation, &flow_equation_deriv);
grouped_polynomial_print(flow_equation_deriv,'%','%');
// free memory
free_Grouped_Polynomial(flow_equation);
free_Grouped_Polynomial(flow_equation_deriv);
free_Int_Array(opts.deriv_vars);
return(0);
}
// 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){
Grouped_Polynomial dflow;
Grouped_Polynomial tmpflow;
Int_Array indices;
int i,j;
int_array_cpy(variables, &indices);
// output polynomial
grouped_polynomial_cpy(flow_equation, flow_equation_derivs);
for(j=0,dflow=flow_equation;j<n;j++){
// tmp flow contains the result of the previous derivative
grouped_polynomial_cpy(dflow, &tmpflow);
// derive
flow_equation_derivx(tmpflow, indices, &dflow);
// free
free_Grouped_Polynomial(tmpflow);
// add the derived indices as variables for the next derivative
for(i=0;i<variables.length;i++){
if(variables.values[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);
}

30
src/meantools_deriv.h Normal file
View File

@ -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

129
src/meantools_eval.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#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<argc;i++){
// flag
if(argv[i][0]=='-'){
for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){
switch(*ptr){
// evaluate string
case 'R':
flag=CP_FLAG_RCCS;
break;
}
}
}
// rccs
else if(flag==CP_FLAG_RCCS){
str_to_char_array((char*)argv[i], &((*opts).eval_rccstring));
flag=0;
}
// read file name from command-line
else{
file=argv[i];
exists_file=1;
}
}
read_config_file(str_args, file, 1-exists_file);
return(0);
}
// evaluate a flow equation on a vector of rccs
int tool_eval(Str_Array str_args, Meantools_Options opts){
// index of the entry in the input file
int arg_index;
// rccs
RCC rccs;
// flow equation
Grouped_Polynomial flow_equation;
// parse flow equation
// if there is a unique argument, assume it is the flow equation
if(str_args.length==1){
char_array_to_Grouped_Polynomial(str_args.strs[0], &flow_equation);
}
else{
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);
}
// rccs
// 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){
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);
}

29
src/meantools_eval.h Normal file
View File

@ -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

130
src/meantools_exp.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#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<exp.length;i++){
for(j=0;j<exp.coefs[i].length;j++){
exp.coefs[i].denoms[j].power=0;
}
}
grouped_polynomial_print(exp,'%','%');
// free memory
free_Fields_Table(fields);
free_Id_Table(idtable);
free_Grouped_Polynomial(exp);
return(0);
}

27
src/meantools_exp.h Normal file
View File

@ -0,0 +1,27 @@
/*
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_EXP_H
#define MEANTOOLS_EXP_H
#include "types.h"
// read arguments
int tool_exp_read_args(int argc, const char* argv[], Str_Array* str_args);
// compute the exponential of the input polynomial
int tool_exp(Str_Array str_args);
#endif

551
src/number.c Normal file
View File

@ -0,0 +1,551 @@
/*
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 "number.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#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<input.length){
fprintf(stderr,"error: trying to copy a number of length %d to another with memory %d\n",input.length, (*output).memory);
exit(-1);
}
for(i=0;i<input.length;i++){
(*output).scalars[i]=input.scalars[i];
(*output).base[i]=input.base[i];
}
(*output).length=input.length;
return(0);
}
// resize memory
int number_resize(Number* number, int newsize){
Number new_number;
init_Number(&new_number, newsize);
number_cpy_noinit(*number,&new_number);
free_Number(*number);
*number=new_number;
return(0);
}
// add a value
int number_append(Q scalar, int base, Number* output){
if((*output).length>=(*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<input.length;i++){
(*output).scalars[offset+i]=input.scalars[i];
(*output).base[offset+i]=input.base[i];
}
(*output).length=offset+input.length;
return(0);
}
// special numbers
Number number_zero(){
Number ret;
// don't allocate 0 memory since zero's will usually be used to add to
init_Number(&ret,NUMBER_SIZE);
return(ret);
}
Number number_one(){
Number ret=number_zero();
number_add_elem(quot(1,1),1,&ret);
return(ret);
}
// find a base element
int number_find_base(int base, Number number){
return(intlist_find(number.base, number.length, base));
}
// sort
int number_sort(Number number, int begin, int end){
int i;
int index;
// the pivot: middle of the array
int pivot=(begin+end)/2;
// if the array is non trivial
if(begin<end){
// send pivot to the end
number_exchange_terms(pivot, end, number);
// loop over the others
for(i=begin, index=begin;i<end;i++){
// compare with pivot
if(number.base[i]<number.base[end]){
// if smaller, exchange with reference index
number_exchange_terms(i, index, number);
// move reference index
index++;
}
}
// put pivot (which we had sent to the end) in the right place
number_exchange_terms(index, end, number);
// recurse
number_sort(number, begin, index-1);
number_sort(number, index+1, end);
}
return(0);
}
// exchange terms (for sorting)
int number_exchange_terms(int index1, int index2, Number number){
Q qtmp;
int tmp;
qtmp=number.scalars[index1];
number.scalars[index1]=number.scalars[index2];
number.scalars[index2]=qtmp;
tmp=number.base[index1];
number.base[index1]=number.base[index2];
number.base[index2]=tmp;
return(0);
}
// checks whether two numbers are equal
// requires both numbers to be sorted
int number_compare(Number x1, Number x2){
int i;
if(x1.length!=x2.length){
return(0);
}
for(i=0;i<x1.length;i++){
if(x1.base[i]!=x2.base[i] || Q_cmp(x1.scalars[i],x2.scalars[i])!=0){
return(0);
}
}
return(1);
}
// add (write result to second element)
int number_add_chain(Number input, Number* inout){
int i;
for(i=0;i<input.length;i++){
number_add_elem(input.scalars[i], input.base[i], inout);
}
return(0);
}
// add a single element
int number_add_elem(Q scalar, int base, Number* inout){
int index;
index=number_find_base(base,*inout);
if(index>=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;i<x1.length;i++){
for(j=0;j<x2.length;j++){
new_scalar=Q_prod(x1.scalars[i], x2.scalars[j]);
// simplify the base
div=gcd(x1.base[i], x2.base[j]);
new_base=(x1.base[i]/div)*(x2.base[j]/div);
new_scalar.numerator*=div;
number_add_elem(new_scalar, new_base, out);
}
}
return(0);
}
// write to second number
int number_prod_chain(Number input, Number* inout){
Number tmp;
number_prod(input,*inout,&tmp);
free_Number(*inout);
*inout=tmp;
return(0);
}
// return result
Number number_prod_ret(Number x1, Number x2){
Number ret;
number_prod(x1,x2,&ret);
return(ret);
}
// multiply by a rational
int number_Qprod_chain(Q q, Number* inout){
int i;
for(i=0;i<(*inout).length;i++){
(*inout).scalars[i]=Q_prod(q,(*inout).scalars[i]);
}
return(0);
}
// write to output
int number_Qprod(Q q, Number x, Number* inout){
number_cpy(x,inout);
number_Qprod_chain(q,inout);
return(0);
}
// return result
Number number_Qprod_ret(Q q, Number x){
Number ret;
number_Qprod(q,x,&ret);
return(ret);
}
// inverse
int number_inverse_inplace(Number* inout){
int i;
for(i=0;i<(*inout).length;i++){
if((*inout).base[i]>0){
(*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;i<in.length;i++){
if(in.scalars[i].numerator!=0){
number_append(in.scalars[i],in.base[i], out);
}
}
return(0);
}
// check whether a number is 0
int number_is_zero(Number x){
int i;
for(i=0;i<x.length;i++){
if(x.scalars[i].numerator!=0 && x.base[i]!=0){
return(0);
}
}
return(1);
}
// approximate numerical expression
long double number_double_val(Number x){
int i;
long double ret=0.;
long double b;
for(i=0;i<x.length;i++){
if(x.scalars[i].numerator!=0){
b=1.0*x.base[i];
ret+=Q_double_value(x.scalars[i])*sqrt(b);
}
}
return(ret);
}
// print to string
int number_sprint(Number number, Char_Array* out){
int i;
Number simp;
number_simplify(number, &simp);
for(i=0;i<simp.length;i++){
if(i>0){
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<length;i++){
(*matrix).matrix[i]=calloc(length,sizeof(Number));
for(j=0;j<length;j++){
(*matrix).matrix[i][j]=number_zero();
}
}
(*matrix).length=length;
return(0);
}
int free_Number_Matrix(Number_Matrix matrix){
int i,j;
for(i=0;i<matrix.length;i++){
for(j=0;j<matrix.length;j++){
free_Number(matrix.matrix[i][j]);
}
free(matrix.matrix[i]);
}
free(matrix.matrix);
free(matrix.indices);
return(0);
}
// Pauli matrices
int Pauli_matrix(int i, Number_Matrix* output){
init_Number_Matrix(output,2);
switch(i){
case 1:
number_add_elem(quot(1,1),1,(*output).matrix[0]+1);
number_add_elem(quot(1,1),1,(*output).matrix[1]+0);
break;
case 2:
number_add_elem(quot(-1,1),-1,(*output).matrix[0]+1);
number_add_elem(quot(1,1),-1,(*output).matrix[1]+0);
break;
case 3:
number_add_elem(quot(1,1),1,(*output).matrix[0]+0);
number_add_elem(quot(-1,1),1,(*output).matrix[1]+1);
break;
default:
fprintf(stderr,"error: requested %d-th pauli matrix\n",i);
exit(-1);
}
return(0);
}

120
src/number.h Normal file
View File

@ -0,0 +1,120 @@
/*
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.
*/
/*
Numerical constants (rational numbers and their square roots (including \sqrt{-1}=i))
*/
#ifndef NUMBER_H
#define NUMBER_H
#include "types.h"
// init
int init_Number(Number* number, int memory);
int free_Number(Number number);
// copy
int number_cpy(Number input, Number* output);
int number_cpy_noinit(Number input, Number* output);
// resize memory
int number_resize(Number* number, int newsize);
// add a value
int number_append(Q scalar, int base, Number* output);
// concatenate
int number_concat(Number input, Number* output);
// special numbers
Number number_zero();
Number number_one();
// find a base element
int number_find_base(int base, Number number);
// sort
int number_sort(Number number, int begin, int end);
// exchange terms (for sorting)
int number_exchange_terms(int index1, int index2, Number number);
// checks whether two numbers are equal
int number_compare(Number x1, Number x2);
// add (write result to second element)
int number_add_chain(Number input, Number* inout);
// add a single element
int number_add_elem(Q scalar, int base, Number* inout);
// create a new number
int number_add(Number x1, Number x2, Number* out);
// return the number
Number number_add_ret(Number x1, Number x_2);
// multiply
int number_prod(Number x1, Number x2, Number* out);
// write to second number
int number_prod_chain(Number input, Number* inout);
// return result
Number number_prod_ret(Number x1, Number x2);
// multiply by a rational
int number_Qprod_chain(Q q, Number* inout);
// write to output
int number_Qprod(Q q, Number x, Number* inout);
// return result
Number number_Qprod_ret(Q q, Number x);
// inverse
int number_inverse_inplace(Number* inout);
// write to output
int number_inverse(Number input, Number* output);
// return result
Number number_inverse_ret(Number x);
// quotient
int number_quot(Number x1, Number x2, Number* output);
int number_quot_chain(Number x1, Number* inout);
Number number_quot_ret(Number x1, Number x2);
// remove 0's
int number_simplify(Number in, Number* out);
// check whether a number is 0
int number_is_zero(Number x);
// approximate numerical expression
long double number_double_val(Number x);
// print to string
int number_sprint(Number number, Char_Array* out);
// print to stdout
int number_print(Number number);
// read from a string
int str_to_Number(char* str, Number* number);
// char array input
int char_array_to_Number(Char_Array cstr_num, Number* number);
//------------------------ Number_Matrix --------------------------
// init
int init_Number_Matrix(Number_Matrix* matrix, int length);
int free_Number_Matrix(Number_Matrix matrix);
// Pauli matrices
int Pauli_matrix(int i, Number_Matrix* output);
#endif

226
src/numkondo.c Normal file
View File

@ -0,0 +1,226 @@
/*
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.
*/
/*
numkondo
Compute the flow of a flow equation numerically
*/
#include <stdio.h>
#include <stdlib.h>
// 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<argc;i++){
// flag
if(argv[i][0]=='-'){
for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){
switch(*ptr){
// final step: display the final step of the integration with maximal precision
case 'F':
(*opts).display_mode=DISPLAY_FINAL;
break;
// niter
case 'N':
flag=CP_FLAG_NITER;
break;
// tolerance
case 'D':
flag=CP_FLAG_TOL;
break;
// initial condition
case 'I':
flag=CP_FLAG_RCCS;
break;
// print version
case 'v':
printf("numkondo " VERSION "\n");
exit(1);
break;
}
}
}
// if the niter flag is up
else if (flag==CP_FLAG_NITER){
// read niter
sscanf(argv[i],"%d",&((*opts).niter));
// reset flag
flag=0;
}
// tolerance
else if (flag==CP_FLAG_TOL){
sscanf(argv[i],"%Lf",&((*opts).tol));
flag=0;
}
// init condition
else if(flag==CP_FLAG_RCCS){
str_to_char_array((char*)argv[i], &((*opts).eval_rccstring));
flag=0;
}
// read file name from command-line
else{
file=argv[i];
exists_file=1;
}
}
read_config_file(str_args, file, 1-exists_file);
return(0);
}
// print usage message
int print_usage_numkondo(){
printf("\nusage:\n numkondo [-F] [-N niter] [-D tolerance] [-I initial_condition] <filename>\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);
}

796
src/parse_file.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#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;j<str_fields.length;j++){
if(comment==1){
if(str_fields.str[j]=='\n'){
comment=0;
}
}
else{
switch(str_fields.str[j]){
// parameters
case 'h':
if(mode==PP_NULL_MODE){
mode=PP_PARAMETER_MODE;
}
break;
// external fields
case 'x':
if(mode==PP_NULL_MODE){
mode=PP_EXTERNAL_MODE;
}
break;
// internal fields
case 'i':
if(mode==PP_NULL_MODE){
mode=PP_INTERNAL_MODE;
}
break;
case 'f':
if(mode==PP_NULL_MODE){
mode=PP_FERMIONS_MODE;
}
break;
// reset buffer
case ':':
buffer_ptr=buffer;
*buffer_ptr='\0';
break;
// write to fields
case ',':
sscanf(buffer,"%d",&i);
if(mode==PP_PARAMETER_MODE){
int_array_append(i,&((*fields).parameter));
}
else if(mode==PP_EXTERNAL_MODE){
int_array_append(i,&((*fields).external));
}
else if(mode==PP_INTERNAL_MODE){
int_array_append(i,&((*fields).internal));
}
else if(mode==PP_FERMIONS_MODE){
int_array_append(i,&((*fields).fermions));
}
buffer_ptr=buffer;
*buffer_ptr='\0';
break;
// back to null mode
case '\n':
sscanf(buffer,"%d",&i);
if(mode==PP_PARAMETER_MODE){
int_array_append(i,&((*fields).parameter));
}
else if(mode==PP_EXTERNAL_MODE){
int_array_append(i,&((*fields).external));
}
else if(mode==PP_INTERNAL_MODE){
int_array_append(i,&((*fields).internal));
}
else if(mode==PP_FERMIONS_MODE){
int_array_append(i,&((*fields).fermions));
}
mode=PP_NULL_MODE;
break;
// comment
case '#':
comment=1;
break;
default:
if(mode!=PP_NULL_MODE){
buffer_ptr=str_addchar(buffer_ptr,str_fields.str[j]);
}
break;
}
}
}
free(buffer);
return(0);
}
// parse symbols list
// write result to fields
int parse_input_symbols(Char_Array str_symbols, Fields_Table* fields){
// buffer
char* buffer=calloc(str_symbols.length+1,sizeof(char));
char* buffer_ptr=buffer;
Polynomial polynomial;
int index;
int i,j;
int mode;
int comment=0;
// loop over input
mode=PP_INDEX_MODE;
for(j=0;j<str_symbols.length;j++){
if(comment==1){
if(str_symbols.str[j]=='\n'){
comment=0;
}
}
// stay in polynomial mode until ','
else if(mode==PP_POLYNOMIAL_MODE){
if(str_symbols.str[j]==','){
// parse polynomial
str_to_Polynomial(buffer, &polynomial);
// write index and polynomial
symbols_append_noinit(index, polynomial, &((*fields).symbols));
mode=PP_INDEX_MODE;
buffer_ptr=buffer;
*buffer_ptr='\0';
}
else{
buffer_ptr=str_addchar(buffer_ptr,str_symbols.str[j]);
}
}
else{
switch(str_symbols.str[j]){
// polynomial mode
case '=':
if(mode==PP_INDEX_MODE){
// read index
sscanf(buffer,"%d",&index);
// reset buffer
buffer_ptr=buffer;
*buffer_ptr='\0';
mode=PP_POLYNOMIAL_MODE;
}
break;
// comment
case '#':
comment=1;
break;
default:
buffer_ptr=str_addchar(buffer_ptr,str_symbols.str[j]);
break;
}
}
}
// last step
if(polynomial.length>0){
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<str_groups.length;j++){
if(comment==1){
if(str_groups.str[j]=='\n'){
comment=0;
}
}
else{
switch(str_groups.str[j]){
// group mode
case '(':
if(mode==PP_NULL_MODE){
// reset buffer
buffer_ptr=buffer;
*buffer_ptr='\0';
// init
init_Int_Array(&group, GROUP_SIZE);
mode=PP_GROUP_MODE;
}
break;
case')':
if(mode==PP_GROUP_MODE){
sscanf(buffer,"%d",&index);
int_array_append(index, &group);
groups_append_noinit(group, groups);
mode=PP_NULL_MODE;
}
break;
// read index
case',':
if(mode==PP_GROUP_MODE){
sscanf(buffer,"%d",&index);
int_array_append(index, &group);
buffer_ptr=buffer;
*buffer_ptr='\0';
}
break;
// comment
case '#':
comment=1;
break;
default:
if(mode!=PP_NULL_MODE){
buffer_ptr=str_addchar(buffer_ptr,str_groups.str[j]);
}
break;
}
}
}
free(buffer);
return(0);
}
// parse identities between fields
// write result to fields
int parse_input_identities(Char_Array str_identities, Fields_Table* fields){
// buffer
char* buffer=calloc(str_identities.length+1,sizeof(char));
char* buffer_ptr=buffer;
Int_Array monomial;
Polynomial polynomial;
int i,j;
int sign;
int tmp;
int mode;
int comment=0;
init_Int_Array(&monomial, MONOMIAL_SIZE);
// loop over input
mode=PP_NULL_MODE;
for(j=0;j<str_identities.length;j++){
if(comment==1){
if(str_identities.str[j]=='\n'){
comment=0;
}
}
// stay in polynomial mode until ','
else if(mode==PP_POLYNOMIAL_MODE){
if(str_identities.str[j]==','){
// parse polynomial
str_to_Polynomial(buffer, &polynomial);
// write monomial and polynomial
identities_append_noinit(monomial, polynomial, &((*fields).ids));
// realloc
init_Int_Array(&monomial, MONOMIAL_SIZE);
mode=PP_NULL_MODE;
}
else{
buffer_ptr=str_addchar(buffer_ptr,str_identities.str[j]);
}
}
else{
switch(str_identities.str[j]){
// monomial
case '[':
if(mode==PP_NULL_MODE){
mode=PP_INDEX_MODE;
}
if(mode==PP_INDEX_MODE){
mode=PP_BRACKET_MODE;
buffer_ptr=buffer;
*buffer_ptr='\0';
}
break;
// for notational homogeneity
case 'f':
if(mode==PP_BRACKET_MODE){
mode=PP_FIELD_MODE;
}
break;
// write monomial term
case ']':
if(mode==PP_FIELD_MODE){
sscanf(buffer,"%d",&i);
int_array_append(i,&monomial);
mode=PP_INDEX_MODE;
}
break;
// polynomial mode
case '=':
if(mode==PP_INDEX_MODE){
buffer_ptr=buffer;
*buffer_ptr='\0';
mode=PP_POLYNOMIAL_MODE;
}
break;
// comment
case '#':
comment=1;
break;
default:
if(mode!=PP_NULL_MODE){
buffer_ptr=str_addchar(buffer_ptr,str_identities.str[j]);
}
break;
}
}
}
// last step
if(mode==PP_POLYNOMIAL_MODE){
str_to_Polynomial(buffer, &polynomial);
identities_append_noinit(monomial, polynomial, &((*fields).ids));
}
else{
free_Int_Array(monomial);
}
// sort
// don't use the identities to simplify
tmp=(*fields).ids.length;
(*fields).ids.length=0;
for(i=0;i<tmp;i++){
sign=1;
monomial_sort((*fields).ids.lhs[i], 0, (*fields).ids.lhs[i].length-1, *fields, &sign);
polynomial_simplify((*fields).ids.rhs+i, *fields);
polynomial_multiply_Qscalar((*fields).ids.rhs[i],quot(sign,1));
}
(*fields).ids.length=tmp;
free(buffer);
return(0);
}
// parse propagator
int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propagator, Fields_Table fields){
// buffer
char* buffer=calloc(str_propagator.length+1,sizeof(char));
char* buffer_ptr=buffer;
int i,j;
int index1=-1;
int index2=-1;
int mode;
int comment=0;
// allocate memory
init_Polynomial_Matrix(propagator, fields.internal.length);
// copy indices
for(i=0;i<fields.internal.length;i++){
(*propagator).indices[i]=fields.internal.values[i];
}
// loop over input
mode=PP_INDEX_MODE;
for(j=0;j<str_propagator.length;j++){
if(comment==1){
if(str_propagator.str[j]=='\n'){
comment=0;
}
}
else{
switch(str_propagator.str[j]){
// indices
case ';':
if(mode==PP_INDEX_MODE){
sscanf(buffer,"%d",&i);
index1=intlist_find_err((*propagator).indices, (*propagator).length, i);
buffer_ptr=buffer;
*buffer_ptr='\0';
}
break;
case ':':
if(mode==PP_INDEX_MODE){
sscanf(buffer,"%d",&i);
index2=intlist_find_err((*propagator).indices, (*propagator).length, i);
buffer_ptr=buffer;
*buffer_ptr='\0';
mode=PP_POLYNOMIAL_MODE;
}
break;
// num
case ',':
if(mode==PP_POLYNOMIAL_MODE && index1>=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<str_polynomial.length;j++){
switch(str_polynomial.str[j]){
case '*':
str_to_Polynomial(buffer, &tmp_poly);
if((*output).length==0){
polynomial_concat(tmp_poly, output);
}
else{
polynomial_prod_chain(tmp_poly, output, fields);
}
free_Polynomial(tmp_poly);
buffer_ptr=buffer;
*buffer_ptr='\0';
break;
default:
buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]);
break;
}
}
//last step
str_to_Polynomial(buffer, &tmp_poly);
if((*output).length==0){
polynomial_concat(tmp_poly, output);
}
else{
polynomial_prod_chain(tmp_poly, output, fields);
}
free_Polynomial(tmp_poly);
free(buffer);
return(0);
}
// parse id table
// fields argument for sorting
int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table fields){
// buffer
char* buffer=calloc(str_idtable.length+1,sizeof(char));
char* buffer_ptr=buffer;
int index;
Polynomial polynomial;
int j;
int mode;
int comment=0;
// allocate memory
init_Id_Table(idtable,EQUATION_SIZE);
// loop over input
mode=PP_INDEX_MODE;
for(j=0;j<str_idtable.length;j++){
if(comment==1){
if(str_idtable.str[j]=='\n'){
comment=0;
}
}
else{
switch(str_idtable.str[j]){
// end polynomial mode
case ',':
// write polynomial
if(mode==PP_POLYNOMIAL_MODE){
str_to_Polynomial(buffer,&polynomial);
// add to idtable
idtable_append_noinit(index,polynomial,idtable);
mode=PP_INDEX_MODE;
buffer_ptr=buffer;
*buffer_ptr='\0';
}
break;
case ':':
if(mode==PP_INDEX_MODE){
sscanf(buffer,"%d",&index);
mode=PP_POLYNOMIAL_MODE;
buffer_ptr=buffer;
*buffer_ptr='\0';
}
break;
// comment
case '#':
comment=1;
break;
default:
if(mode!=PP_NULL_MODE){
buffer_ptr=str_addchar(buffer_ptr,str_idtable.str[j]);
}
break;
}
}
}
//last step
if(mode==PP_POLYNOMIAL_MODE){
str_to_Polynomial(buffer,&polynomial);
idtable_append_noinit(index,polynomial,idtable);
}
// sort
for(j=0;j<(*idtable).length;j++){
polynomial_simplify((*idtable).polynomials+j, fields);
}
free(buffer);
return(0);
}
// parse a list of labels
int parse_labels(Char_Array str_labels, Labels* labels){
// buffer
char* buffer=calloc(str_labels.length+1,sizeof(char));
char* buffer_ptr=buffer;
Char_Array label;
int index;
int j;
int mode;
int comment=0;
// allocate memory
init_Labels(labels,EQUATION_SIZE);
// loop over input
mode=PP_INDEX_MODE;
for(j=0;j<str_labels.length;j++){
if(comment==1){
if(str_labels.str[j]=='\n'){
comment=0;
}
}
else{
switch(str_labels.str[j]){
case '"':
if(mode==PP_INDEX_MODE){
mode=PP_LABEL_MODE;
buffer_ptr=buffer;
*buffer_ptr='\0';
}
// write
else if(mode==PP_LABEL_MODE){
str_to_char_array(buffer,&label);
labels_append_noinit(label,index,labels);
mode=PP_INDEX_MODE;
buffer_ptr=buffer;
*buffer_ptr='\0';
}
break;
case ':':
// write
if(mode==PP_INDEX_MODE){
sscanf(buffer,"%d",&index);
}
break;
// characters to ignore
case ' ':break;
case '&':break;
case '\n':break;
case ',':break;
// comment
case '#':
comment=1;
break;
default:
buffer_ptr=str_addchar(buffer_ptr,str_labels.str[j]);
break;
}
}
}
free(buffer);
return(0);
}
// read initial condition for numerical computation
int parse_init_cd(Char_Array init_cd, RCC* init){
char* buffer=calloc(init_cd.length+1,sizeof(char));
char* buffer_ptr=buffer;
int index=0;
int i,j;
int comment_mode=0;
int dcount=0;
*buffer_ptr='\0';
// loop over the input
for(j=0;j<init_cd.length;j++){
if(comment_mode==1){
if(init_cd.str[j]=='\n'){
comment_mode=0;
}
}
else{
switch(init_cd.str[j]){
// new term
case ',':
// write init
sscanf(buffer,"%Lf",(*init).values+intlist_find_err((*init).indices,(*init).length,index));
// reset buffer
buffer_ptr=buffer;
*buffer_ptr='\0';
// reset derivatives counter
dcount=0;
break;
// separator
case ':':
// write index
sscanf(buffer,"%d",&i);
if(i<0){
index=i-dcount*DOFFSET;
}
else{
index=i+dcount*DOFFSET;
}
buffer_ptr=buffer;
*buffer_ptr='\0';
break;
// derivatives
case 'd':
dcount++;
break;
// characters to ignore
case ' ':break;
case '&':break;
case '\n':break;
// comments
case '#':
comment_mode=1;
break;
default:
// write to buffer
buffer_ptr=str_addchar(buffer_ptr,init_cd.str[j]);
break;
}
}
}
// write init
sscanf(buffer,"%Lf",(*init).values+unlist_find((*init).indices,(*init).length,index));
free(buffer);
return(0);
}
// set indices and length of init
int prepare_init(int* indices, int length, RCC* init){
int i;
init_RCC(init, length);
for(i=0;i<length;i++){
(*init).indices[i]=indices[i];
// set constants to 1
if(indices[i]<0 && indices[i]>-DOFFSET){
(*init).values[i]=1.;
}
else{
(*init).values[i]=0.;
}
}
return(0);
}

56
src/parse_file.h Normal file
View File

@ -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

1263
src/polynomial.c Normal file
View File

File diff suppressed because it is too large Load Diff

131
src/polynomial.h Normal file
View File

@ -0,0 +1,131 @@
/*
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.
*/
/*
Represent a polynomial as a list of monomials with factors
*/
#ifndef POLYNOMIAL_H
#define POLYNOMIAL_H
#include "types.h"
// allocate memory
int init_Polynomial(Polynomial* polynomial,int size);
// free memory
int free_Polynomial(Polynomial polynomial);
// resize the memory allocated to a polynomial
int resize_polynomial(Polynomial* polynomial,int new_size);
// copy a polynomial
int polynomial_cpy(Polynomial input, Polynomial* output);
int polynomial_cpy_noinit(Polynomial input, Polynomial* output);
// append an element to a polynomial
int polynomial_append(Int_Array monomial, Int_Array factor, Number num, Polynomial* output);
int polynomial_append_noinit(Int_Array monomial, Int_Array factor, Number num, Polynomial* output);
// 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);
// concatenate two polynomials
int polynomial_concat(Polynomial input, Polynomial* output);
int polynomial_concat_noinit(Polynomial input, Polynomial* output);
int polynomial_concat_noinit_inplace(Polynomial input, Polynomial* output);
// add polynomials
int polynomial_add_chain(Polynomial input, Polynomial* inout, Fields_Table fields);
// add polynomials (noinit)
int polynomial_add_chain_noinit(Polynomial input, Polynomial* inout, Fields_Table fields);
// multiply a polynomial by a scalar
int polynomial_multiply_scalar(Polynomial polynomial, Number num);
int polynomial_multiply_Qscalar(Polynomial polynomial, Q q);
// change the sign of the monomials in a polynomial
int polynomial_conjugate(Polynomial polynomial);
// returns an initialized polynomial, equal to 1
Polynomial polynomial_one();
// returns an initialized polynomial, equal to 0
Polynomial polynomial_zero();
// check whether a polynomial is 0
int polynomial_is_zero(Polynomial poly);
// compute V^p
int polynomial_power(Polynomial input_polynomial, Polynomial* output, int power, Fields_Table fields);
// compute V*W
int polynomial_prod(Polynomial input1, Polynomial input2, Polynomial* output, Fields_Table fields);
int polynomial_prod_chain_nosimplify(Polynomial input, Polynomial* inout, Fields_Table fields);
int polynomial_prod_chain(Polynomial input, Polynomial* inout, Fields_Table fields);
// exp(V)
int polynomial_exponential(Polynomial input_polynomial,Polynomial* output, Fields_Table fields);
// log(1+W)
int polynomial_logarithm(Polynomial input_polynomial, Polynomial* output, Fields_Table fields);
// check whether a monomial vanishes
int check_monomial(Int_Array monomial, Fields_Table fields);
// check whether the product of two monomials will vanish
int check_monomial_willnot_vanish(Int_Array monomial1, Int_Array monomial2, Fields_Table fields);
// check whether one can add a term to a monomial without creating repetitions
int check_monomial_addterm(Int_Array monomial, Int_Array term, Fields_Table fields);
// check whether a monomial vanishes or has unmatched +/- fields
int check_monomial_match(Int_Array monomial, Fields_Table fields);
// remove terms with more plus internal fields than minus ones
int remove_unmatched_plusminus(Polynomial* polynomial, Fields_Table fields);
// denominator of a multinomal factor: m1!m2!...
int multinomial(int power,int* terms);
// simplify a Polynomial
int polynomial_simplify(Polynomial* polynomial, Fields_Table fields);
// sort a polynomial
int polynomial_sort(Polynomial* polynomial, int begin, int end);
// exchange two terms (for the sorting algorithm)
int exchange_polynomial_terms(int i, int j, Polynomial* polynomial);
// sort a monomial (with sign coming from exchanging two Fermions)
int monomial_sort(Int_Array monomial, int begin, int end, Fields_Table fields, int* sign);
// order fields: parameter, external, internal
int compare_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table fields);
// exchange two fields (with sign)
int exchange_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table fields, int* sign);
// sort a monomial by putting each group together
int monomial_sort_groups(Int_Array monomial, int begin, int end, Fields_Table fields, Groups groups, int* sign);
// order fields: group, then parameter, external, internal
int compare_monomial_terms_groups(Int_Array monomial, int pos1, int pos2, Fields_Table fields, Groups groups);
// convert and idtable to a polynomial
int idtable_to_polynomial(Id_Table idtable, Polynomial* polynomial);
// replace the factors in a polynomial as prescribed by an equation in the form of a Grouped_Polynomial
int replace_factors(Grouped_Polynomial equations, Polynomial* polynomial);
// print a polynomial
int polynomial_sprint(Polynomial polynomial, Char_Array* output);
int polynomial_print(Polynomial polynomial);
// read a polynomial
int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output);
int str_to_Polynomial(char* str_polynomial,Polynomial* output);
//------------------------ Polynomial_Matrix --------------------------
// init
int init_Polynomial_Matrix(Polynomial_Matrix* matrix, int length);
int free_Polynomial_Matrix(Polynomial_Matrix matrix);
#endif

23
src/rational.h Normal file
View File

@ -0,0 +1,23 @@
/*
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.
*/
/*
represent rational numbers either as a quotient of 64-bit integers, or a quotient of 80-bit floats (with 64-bit precision)
*/
// include both, only one is non empty
#include "rational_float.h"
#include "rational_int.h"

196
src/rational_float.c Normal file
View File

@ -0,0 +1,196 @@
/*
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.
*/
#ifdef RATIONAL_AS_FLOAT
#include "rational_float.h"
#include <stdio.h>
#include <stdlib.h>
#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

64
src/rational_float.h Normal file
View File

@ -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

190
src/rational_int.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#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

59
src/rational_int.h Normal file
View File

@ -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

95
src/rcc.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#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<input.length;i++){
RCC_set_elem(input.values[i], input.indices[i], output, i);
}
return(0);
}
// concatenate rccs
int RCC_concat(RCC rccs1, RCC rccs2, RCC* output){
int i;
init_RCC(output,rccs1.length+rccs2.length);
for(i=0;i<rccs1.length;i++){
RCC_set_elem(rccs1.values[i], rccs1.indices[i], output, i);
}
for(i=0;i<rccs2.length;i++){
RCC_set_elem(rccs2.values[i], rccs2.indices[i], output, i+rccs1.length);
}
return(0);
}
// append an rcc at the end of another
int RCC_append(RCC input, RCC* output){
int i;
for(i=0;i<input.length;i++){
RCC_set_elem(input.values[i], input.indices[i], output, i+(*output).length);
}
(*output).length+=input.length;
return(0);
}
// print an rcc vector with maximal precision
int RCC_print(RCC rccs){
int j;
for(j=0;j<rccs.length;j++){
printf("%d:%.19Le",rccs.indices[j],rccs.values[j]);
if(j<rccs.length-1){
printf(",\n");
}
else{
printf("\n");
}
}
return(0);
}

39
src/rcc.h Normal file
View File

@ -0,0 +1,39 @@
/*
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.
*/
/* RCC struct */
#ifndef RCC_H
#define RCC_H
#include "types.h"
// init
int init_RCC(RCC* rccs, int size);
int free_RCC(RCC rccs);
// set an element of an rcc
int RCC_set_elem(long double value, int index, RCC* rcc, int pos);
// copy
int RCC_cpy(RCC input,RCC* output);
// concatenate 2 rccs
int RCC_concat(RCC rccs1, RCC rccs2, RCC* output);
// append an rcc to another
int RCC_append(RCC input, RCC* output);
// print an rcc vector with maximal precision
int RCC_print(RCC rccs);
#endif

142
src/tools.c Normal file
View File

@ -0,0 +1,142 @@
/*
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 "tools.h"
#include <stdio.h>
#include <stdlib.h>
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<size;i++){
if(list[i]==x){return(i);}
}
fprintf(stderr,"error: element %d not found\n",x);
exit(-1);
}
// find an element in a list whose first element is not its number of elements (skips first element)
// no error reporting
int unlist_find_noerr(int* list, int size, int x){
int i;
for(i=1;i<size;i++){
if(list[i]==x){return(i);}
}
return(-1);
}
// find an element in a list (checks first element)
int intlist_find(int* list, int size, int x){
int i;
for(i=0;i<size;i++){
if(list[i]==x){return(i);}
}
return(-1);
}
// with error
int intlist_find_err(int* list, int size, int x){
int i;
for(i=0;i<size;i++){
if(list[i]==x){return(i);}
}
fprintf(stderr,"error: element %d not found\n",x);
exit(-1);
}
// compare two lists
int list_cmp(int* list1, int* list2){
if(*list1!=*list2){
return(0);
}
int* ptr1=list1+1;
int* ptr2=list2+1;
int i;
for(i=1;i<=*list1;i++){
if(*ptr1!=*ptr2){
return(0);
}
ptr1++;
ptr2++;
}
return(1);
}
// max and min
int max(int x1, int x2){
if(x1>=x2){
return(x1);
}
else{
return(x2);
}
}
int min(int x1, int x2){
if(x1<=x2){
return(x1);
}
else{
return(x2);
}
}

49
src/tools.h Normal file
View File

@ -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

219
src/types.h Normal file
View File

@ -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