diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/definitions.cpp | 2 | ||||
| -rw-r--r-- | src/fields.c | 125 | ||||
| -rw-r--r-- | src/fields.h | 4 | ||||
| -rw-r--r-- | src/kondo.c | 150 | ||||
| -rw-r--r-- | src/kondo.h | 2 | ||||
| -rw-r--r-- | src/mean.c | 8 | ||||
| -rw-r--r-- | src/parse_file.c | 24 | ||||
| -rw-r--r-- | src/polynomial.c | 60 | ||||
| -rw-r--r-- | src/polynomial.h | 8 | ||||
| -rw-r--r-- | src/types.h | 2 |
10 files changed, 266 insertions, 119 deletions
diff --git a/src/definitions.cpp b/src/definitions.cpp index 982f7a7..1884488 100644 --- a/src/definitions.cpp +++ b/src/definitions.cpp @@ -17,7 +17,7 @@ limitations under the License. #ifndef DEFINITIONS_GCC #define DEFINITIONS_GCC -#define VERSION "1.2.1" +#define VERSION "1.3" // number of entries in a configuration file #define ARG_COUNT 10 diff --git a/src/fields.c b/src/fields.c index 1b221f2..bfd1f39 100644 --- a/src/fields.c +++ b/src/fields.c @@ -32,6 +32,7 @@ int init_Fields_Table(Fields_Table* fields){ init_Identities(&((*fields).ids), FIELDS_SIZE); init_Symbols(&((*fields).symbols), FIELDS_SIZE); init_Int_Array(&((*fields).fermions),FIELDS_SIZE); + init_Int_Array(&((*fields).noncommuting),FIELDS_SIZE); return(0); } int free_Fields_Table(Fields_Table fields){ @@ -41,6 +42,7 @@ int free_Fields_Table(Fields_Table fields){ free_Identities(fields.ids); free_Symbols(fields.symbols); free_Int_Array(fields.fermions); + free_Int_Array(fields.noncommuting); return(0); } @@ -73,6 +75,16 @@ int is_fermion(int index, Fields_Table fields){ } } +// check whether a field is non-commuting +int is_noncommuting(int index, Fields_Table fields){ + if(int_array_find(abs(index), fields.noncommuting)>=0){ + return(1); + } + else{ + return(0); + } +} + // ------------------ Identities -------------------- @@ -180,13 +192,16 @@ int identities_concat(Identities input, Identities* output){ // resolve the identities // requires both the monomials in polynomial and the ids in fields to be sorted +// IMPORTANT: the sorting must be such that noncommuting fields must come before the other fields int resolve_ids(Polynomial* polynomial, Fields_Table fields){ int i,j,k,l; int sign; int fermion_count; + int first_field; int at_least_one; int security; - Int_Array monomial; + Int_Array pre_monomial; + Int_Array post_monomial; Number num; Number tmp_num; @@ -207,29 +222,38 @@ int resolve_ids(Polynomial* polynomial, Fields_Table fields){ // 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 + first_field=int_array_is_subarray_noncommuting(fields.ids.lhs[j],(*polynomial).monomials[i],fields); + if(first_field>=0){ + init_Int_Array(&pre_monomial, (*polynomial).monomials[i].length); + init_Int_Array(&post_monomial, (*polynomial).monomials[i].length); + + // add whatever is before the first field to pre + for(k=0;k<first_field;k++){ + int_array_append((*polynomial).monomials[i].values[k],&pre_monomial); + } + + // find others and move them together + // sign from moving the fields sign=1; - // number of Fermions to remove from the monomial + // number of Fermions to jump over fermion_count=0; - for(k=0,l=0;k<(*polynomial).monomials[i].length;k++){ + for(l=1,k=first_field+1;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; + // add to post + int_array_append((*polynomial).monomials[i].values[k],&post_monomial); + // count Fermions to jump + if(is_fermion((*polynomial).monomials[i].values[k],fields)){ + fermion_count++; } } else{ - // increment fermion_count - if(is_fermion(fields.ids.lhs[j].values[l],fields)){ - fermion_count++; + // sign correction + if(is_fermion(fields.ids.lhs[j].values[l],fields) && fermion_count % 2 == 1){ + sign*=-1; } // increment "current" field in the id l++; @@ -240,30 +264,33 @@ int resolve_ids(Polynomial* polynomial, Fields_Table fields){ // 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); + polynomial_append(pre_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); + int_array_concat(post_monomial,(*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); + monomial_sort((*polynomial).monomials[(*polynomial).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; + (*polynomial).monomials[i]=pre_monomial; int_array_concat(fields.ids.rhs[j].monomials[0],(*polynomial).monomials+i); + int_array_concat(post_monomial,(*polynomial).monomials+i); + free_Int_Array(post_monomial); // re-sort monomial sign=1; - monomial_sort((*polynomial).monomials[i],0,(*polynomial).monomials[i].length-1,fields,&sign); + monomial_sort((*polynomial).monomials[i],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--; + // repeat the replacement (in order to perform all of the replacements if several are necessary) + j=0; if(at_least_one==0){ at_least_one=1; } @@ -275,6 +302,62 @@ int resolve_ids(Polynomial* polynomial, Fields_Table fields){ return(0); } +// check whether an array is a sub-array of another +// requires noncommuting elements to be next to each other +// other elements may be separated, but the order must be respected +// returns the first index of the sub-array +// IMPORTANT: the noncommuting elements must precede all others in input and in test_array +int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fields_Table fields){ + int i,j; + int matches=0; + int post_nc=0; + int match_nc; + int first=-1; + + // bound noncommuting elements + while(is_noncommuting(input.values[post_nc], fields)==1){ + post_nc++; + } + + for(i=0,match_nc=0;i<test_array.length;i++){ + if(test_array.values[i]==input.values[0]){ + match_nc=1; + } + for(j=1;j<post_nc;j++){ + if(test_array.values[i+j]!=input.values[j]){ + match_nc=0; + } + } + if(match_nc==1){ + first=i; + break; + } + } + + if(first<0){ + return(-1); + } + + if(post_nc>0){ + matches=post_nc; + } + else{ + matches=1; + } + + for(i=first+1;i<test_array.length && matches<input.length;i++){ + if(input.values[matches]==test_array.values[i]){ + matches++; + } + } + if(matches==input.length){ + return(first); + } + else{ + return(-1); + } +} + // ------------------ Symbols -------------------- diff --git a/src/fields.h b/src/fields.h index 3795d92..b6f05c2 100644 --- a/src/fields.h +++ b/src/fields.h @@ -29,6 +29,8 @@ int free_Fields_Table(Fields_Table fields); int field_type(int index, Fields_Table fields); // check whether a field anticommutes int is_fermion(int index, Fields_Table fields); +// check whether a field is non-commuting +int is_noncommuting(int index, Fields_Table fields); // init @@ -51,6 +53,8 @@ int identities_concat(Identities input, Identities* output); // resolve the identities int resolve_ids(Polynomial* polynomial, Fields_Table fields); +// check whether an array is a sub-array of another, support for noncommuting elements +int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fields_Table fields); // init diff --git a/src/kondo.c b/src/kondo.c index 39d74bb..79b7c56 100644 --- a/src/kondo.c +++ b/src/kondo.c @@ -29,16 +29,17 @@ limitations under the License. #include "definitions.cpp" #include "rational.h" -// dimension of A, B and h (must be <10) +// dimension of A, B, h and t (must be <10) #define KONDO_DIM 3 // number of spin components #define KONDO_SPIN 2 -// offsets for indices of A, B and h +// offsets for indices of A, B, h and t // order matters for symbols table #define KONDO_A_OFFSET 1 #define KONDO_B_OFFSET 2 #define KONDO_H_OFFSET 3 +#define KONDO_T_OFFSET 4 // parsing modes (from parse_file.c) #define PP_NULL_MODE 0 @@ -193,13 +194,19 @@ int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fiel // parameters char_array_append_str("h:",str_fields); + // h for(i=0;i<KONDO_DIM;i++){ - char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_H_OFFSET)); + char_array_snprintf(str_fields, "%d,", 10*(i+10*KONDO_H_OFFSET)); + } + // t + for(i=0;i<KONDO_DIM;i++){ + char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_T_OFFSET)); if(i<KONDO_DIM-1){ - char_array_append(',',str_fields); + char_array_append(',', str_fields); } } - char_array_append('\n',str_fields); + char_array_append('\n', str_fields); + // declare Fermions char_array_append_str("f:",str_fields); @@ -226,6 +233,16 @@ int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fiel } char_array_append('\n',str_fields); + // declare noncommuting + char_array_append_str("a:",str_fields); + for(i=0;i<KONDO_DIM;i++){ + char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_T_OFFSET)); + if(i<KONDO_DIM-1){ + char_array_append(',',str_fields); + } + } + char_array_append('\n', str_fields); + // parse fields table parse_input_fields(*str_fields, fields); @@ -245,7 +262,7 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){ // loop over box index for(i=1;i<=box_count;i++){ - // loop over letters + // loop over letters (A and B) for(j=0;j<2;j++){ // loop over space dimension for(k=0;k<KONDO_DIM;k++){ @@ -255,7 +272,7 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){ init_Char_Array(&tmp_str, 6); char_array_snprintf(&tmp_str, "%c%d%d", letters[j], k, i); // compute corresponding polynomial - kondo_resolve_ABh(tmp_str.str, &poly, *fields); + kondo_resolve_ABht(tmp_str.str, &poly, *fields); free_Char_Array(tmp_str); // write to output polynomial_sprint(poly, str_symbols); @@ -320,45 +337,6 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){ return(0); } -// 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){ - int i,j,k; - Char_Array tmp_str; - Polynomial poly; - char letters[3]={'A','B','h'}; - - init_Char_Array(str_symbols, STR_SIZE); - char_array_snprintf(str_symbols, "#!symbols\n"); - - // loop over box index - for(i=1;i<=box_count;i++){ - // loop over letters - for(j=0;j<3;j++){ - for(k=0;k<3;k++){ - // write index - char_array_snprintf(str_symbols, "%d=", 100*(10*(KONDO_A_OFFSET+j)+KONDO_A_OFFSET+k)+i); - // write the name of the scalar product - init_Char_Array(&tmp_str, 6); - char_array_snprintf(&tmp_str, "%c%d.%c%d", letters[j], i, letters[k], i); - // compute corresponding polynomial - kondo_resolve_scalar_prod(tmp_str.str, &poly, *fields); - free_Char_Array(tmp_str); - // write to output - polynomial_sprint(poly, str_symbols); - free_Polynomial(poly); - // add , - if(i<box_count || j<2 || k<2){ - char_array_snprintf(str_symbols,",\n"); - } - } - } - } - - parse_input_symbols(*str_symbols, fields); - - return(0); -} - // generate Kondo groups (groups of independent variables) int kondo_groups(Char_Array* str_groups, int box_count){ @@ -395,12 +373,21 @@ int kondo_groups(Char_Array* str_groups, int box_count){ // generate Kondo identities // h_3^2=1-h_1^2-h_2^2 +// and Pauli matrices int kondo_identities(Char_Array* str_identities){ int i; init_Char_Array(str_identities,STR_SIZE); char_array_snprintf(str_identities, "#!identities\n"); + // Pauli matrices + for(i=0;i<KONDO_DIM;i++){ + char_array_snprintf(str_identities,"[f%d][f%d]=(1),\n",10*(10*KONDO_T_OFFSET+i),10*(10*KONDO_T_OFFSET+i)); + char_array_snprintf(str_identities,"[f%d][f%d]=(s{-1})[f%d],\n",10*(10*KONDO_T_OFFSET+i),10*(10*KONDO_T_OFFSET+(i+1)%3),10*(10*KONDO_T_OFFSET+(i+2)%3)); + char_array_snprintf(str_identities,"[f%d][f%d]=((-1)s{-1})[f%d],\n",10*(10*KONDO_T_OFFSET+(i+2)%3),10*(10*KONDO_T_OFFSET+(i+1)%3),10*(10*KONDO_T_OFFSET+i)); + } + + // h char_array_snprintf(str_identities, "[f%d][f%d]=(1)",10*(KONDO_DIM-1+10*KONDO_H_OFFSET),10*(KONDO_DIM-1+10*KONDO_H_OFFSET)); for(i=0;i<KONDO_DIM-1;i++){ char_array_snprintf(str_identities, "+(-1)[f%d][f%d]",10*(i+10*KONDO_H_OFFSET),10*(i+10*KONDO_H_OFFSET)); @@ -809,11 +796,10 @@ int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_ if(tmp_poly.length>0){ for(i=0;i<tmp_poly.length;i++){ if(mode==PP_FIELD_SCALAR_MODE){ - if(offset1!=KONDO_H_OFFSET || offset2!=KONDO_H_OFFSET){ - int_array_append(1000*(10*offset1+offset2)+index, tmp_poly.monomials+i); - } + int_array_append(1000*(10*offset1+offset2)+index, tmp_poly.monomials+i); } else{ + // vector product int_array_append(100*(100*KONDO_A_OFFSET+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+index, tmp_poly.monomials+i); } } @@ -822,28 +808,15 @@ int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_ else{ init_Int_Array(&tmp_monomial, MONOMIAL_SIZE); if(mode==PP_FIELD_SCALAR_MODE){ - if(offset1!=KONDO_H_OFFSET || offset2!=KONDO_H_OFFSET){ - int_array_append(1000*(10*offset1+offset2)+index, &tmp_monomial); - } + int_array_append(1000*(10*offset1+offset2)+index, &tmp_monomial); } else{ + // vector product int_array_append(100*(100*KONDO_A_OFFSET+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+index, &tmp_monomial); } init_Int_Array(&dummy_factor, 1); polynomial_append_noinit(tmp_monomial, dummy_factor, number_one(), &tmp_poly); } - /* // older method in which a scalar product was expanded in A, B and h - // resolve scalar product - kondo_resolve_scalar_prod_symbols(buffer, &scalar_prod_poly); - // add to tmp_poly - if(tmp_poly.length==0){ - polynomial_concat(scalar_prod_poly,&tmp_poly); - } - else{ - polynomial_prod_chain(scalar_prod_poly,&tmp_poly,fields); - } - free_Polynomial(scalar_prod_poly); - */ } // switch back to null mode mode=PP_NULL_MODE; @@ -941,10 +914,10 @@ int parse_kondo_polynomial(Char_Array kondo_polynomial_str, Polynomial* polynomi } -// 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 Aij, Bij, hi, ti where i is a space dimension and j is a box index +int kondo_resolve_ABht(char* str, Polynomial* output, Fields_Table fields){ char* ptr; - // offset (A,B or H) + // offset (A,B, H or T) int offset=-1; // dimension int dim=-1; @@ -978,6 +951,9 @@ int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields){ case 'h': offset=KONDO_H_OFFSET; break; + case 't': + offset=KONDO_T_OFFSET; + break; default: // set index if dim was already set if(dim>=0){ @@ -1001,8 +977,8 @@ int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields){ } } - // h's - if(offset==KONDO_H_OFFSET){ + // h's and t's + if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){ // external field init_Int_Array(&monomial,1); init_Int_Array(&factor,1); @@ -1062,7 +1038,7 @@ int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields){ // read a Kondo scalar product (generalized to vector products as well) int kondo_resolve_scalar_prod(char* str, Polynomial* output, Fields_Table fields){ char* ptr; - // offset of each term (A,B or H) + // offset of each term (A,B,H or T) int offset=-1; // index of each term (0,...,box_count) int index=0; @@ -1091,6 +1067,9 @@ int kondo_resolve_scalar_prod(char* str, Polynomial* output, Fields_Table fields case 'h': offset=KONDO_H_OFFSET; break; + case 't': + offset=KONDO_T_OFFSET; + break; // scalar product case '.': @@ -1195,8 +1174,8 @@ int kondo_polynomial_vector(int offset, int index, Polynomial (*polys)[3], Field init_Polynomial((*polys)+i,POLY_SIZE); } - // h's - if(offset==KONDO_H_OFFSET){ + // h's and t's + if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){ // construct every component field for(i=0;i<KONDO_DIM;i++){ // external field @@ -1261,7 +1240,7 @@ int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output){ char* ptr; // first or second term int term=0; - // offset of each term (A,B or H) + // offset of each term (A,B,H or T) int offset[2]; // index of each term (0,...,box_count) int index[2]={0,0}; @@ -1289,6 +1268,9 @@ int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output){ case 'h': offset[term]=KONDO_H_OFFSET; break; + case 't': + offset[term]=KONDO_T_OFFSET; + break; // switch term case '.': term=1-term; @@ -1304,13 +1286,13 @@ int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output){ init_Int_Array(&monomial,2); init_Int_Array(&factor, 1); - if(offset[0]==KONDO_H_OFFSET){ + if(offset[0]==KONDO_H_OFFSET || offset[0]==KONDO_T_OFFSET){ int_array_append(10*(10*offset[0]+i)+index[0], &monomial); } else{ int_array_append(100*(10*offset[0]+i)+index[0], &monomial); } - if(offset[1]==KONDO_H_OFFSET){ + if(offset[1]==KONDO_H_OFFSET || offset[1]==KONDO_T_OFFSET){ int_array_append(10*(10*offset[1]+i)+index[1], &monomial); } else{ @@ -1344,6 +1326,9 @@ int get_offset_index(char* str, int* offset, int* index){ case 'h': *offset=KONDO_H_OFFSET; break; + case 't': + *offset=KONDO_T_OFFSET; + break; default: // char to int *index=*ptr-'0'; @@ -1378,6 +1363,9 @@ int get_offsets_index(char* str, int* offset1, int* offset2, int* index){ case 'h': offset[term]=KONDO_H_OFFSET; break; + case 't': + offset[term]=KONDO_T_OFFSET; + break; // switch term case '.': term=1-term; @@ -1391,6 +1379,11 @@ int get_offsets_index(char* str, int* offset1, int* offset2, int* index){ *offset1=offset[0]; *offset2=offset[1]; + // if no A's or B's, then index=0 + if((offset[0]==KONDO_H_OFFSET || offset[0]==KONDO_T_OFFSET) && (offset[1]==KONDO_H_OFFSET || offset[1]==KONDO_T_OFFSET)){ + *index=0; + } + return(0); } @@ -1429,6 +1422,9 @@ int get_symbol_index(char* str){ case 'h': offset=KONDO_H_OFFSET; break; + case 't': + offset=KONDO_T_OFFSET; + break; default: // set index if dim was already set if(dim>=0){ @@ -1443,8 +1439,8 @@ int get_symbol_index(char* str){ if(offset==-1){ return(-1); } - // no symbol for h - if(offset==KONDO_H_OFFSET){ + // no symbol for h or t + if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){ return(10*(10*offset+dim)); } else{ diff --git a/src/kondo.h b/src/kondo.h index 23756ad..c534145 100644 --- a/src/kondo.h +++ b/src/kondo.h @@ -55,7 +55,7 @@ int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, 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); +int kondo_resolve_ABht(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 @@ -42,7 +42,7 @@ int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Ma *out=polynomial_one(); // sort first - monomial_sort(monomial, 0, monomial.length-1, fields, &sign); + monomial_sort(monomial, 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 @@ -417,7 +417,7 @@ int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Po if(check_monomial_match(tmp_monomial, fields)==1){ // sort monomial sign=1; - monomial_sort(tmp_monomial, 0, tmp_monomial.length-1, fields, &sign); + monomial_sort(tmp_monomial, fields, &sign); number_Qprod_chain(quot(sign,1), &tmp_num); // mean @@ -628,7 +628,7 @@ int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Pol init_Polynomial(output, MONOMIAL_SIZE); // check whether there are symbols - // requires the symbols to be at the end of the monomial + // IMPORTANT: the symbols must 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); @@ -639,7 +639,7 @@ int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Pol // sort into groups if(groups.length>0){ sign=1; - monomial_sort_groups(monomial, 0, monomial.length-1, fields, groups, &sign); + monomial_sort_groups(monomial, fields, groups, &sign); } // construct groups and take mean init_Int_Array(&tmp_monomial, MONOMIAL_SIZE); diff --git a/src/parse_file.c b/src/parse_file.c index 6054372..19b91c6 100644 --- a/src/parse_file.c +++ b/src/parse_file.c @@ -44,16 +44,17 @@ limitations under the License. #define PP_EXTERNAL_MODE 8 #define PP_INTERNAL_MODE 9 #define PP_FERMIONS_MODE 10 +#define PP_NONCOMMUTING_MODE 11 // indices -#define PP_INDEX_MODE 11 +#define PP_INDEX_MODE 12 // factors or monomials -#define PP_BRACKET_MODE 12 +#define PP_BRACKET_MODE 13 // labels -#define PP_LABEL_MODE 13 +#define PP_LABEL_MODE 14 // polynomial -#define PP_POLYNOMIAL_MODE 14 +#define PP_POLYNOMIAL_MODE 15 // group -#define PP_GROUP_MODE 15 +#define PP_GROUP_MODE 16 // parse fields list @@ -101,6 +102,11 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){ mode=PP_FERMIONS_MODE; } break; + case 'a': + if(mode==PP_NULL_MODE){ + mode=PP_NONCOMMUTING_MODE; + } + break; // reset buffer case ':': @@ -123,6 +129,9 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){ else if(mode==PP_FERMIONS_MODE){ int_array_append(i,&((*fields).fermions)); } + else if(mode==PP_NONCOMMUTING_MODE){ + int_array_append(i,&((*fields).noncommuting)); + } buffer_ptr=buffer; *buffer_ptr='\0'; break; @@ -142,6 +151,9 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){ else if(mode==PP_FERMIONS_MODE){ int_array_append(i,&((*fields).fermions)); } + else if(mode==PP_NONCOMMUTING_MODE){ + int_array_append(i,&((*fields).noncommuting)); + } mode=PP_NULL_MODE; break; @@ -417,7 +429,7 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){ (*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); + monomial_sort((*fields).ids.lhs[i], *fields, &sign); polynomial_simplify((*fields).ids.rhs+i, *fields); polynomial_multiply_Qscalar((*fields).ids.rhs[i],quot(sign,1)); } diff --git a/src/polynomial.c b/src/polynomial.c index 639728a..d082b60 100644 --- a/src/polynomial.c +++ b/src/polynomial.c @@ -688,7 +688,7 @@ int polynomial_simplify(Polynomial* polynomial, Fields_Table fields){ // sort monomials and factors for(i=0;i<(*polynomial).length;i++){ sign=1; - monomial_sort((*polynomial).monomials[i],0,(*polynomial).monomials[i].length-1,fields,&sign); + monomial_sort((*polynomial).monomials[i],fields,&sign); number_Qprod_chain(quot(sign,1),(*polynomial).nums+i); int_array_sort((*polynomial).factors[i],0,(*polynomial).factors[i].length-1); } @@ -787,7 +787,30 @@ 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){ +// if the monomial contains noncommuting elements, put them at the beginning of the monomial +int monomial_sort(Int_Array monomial, Fields_Table fields, int* sign){ + int i,j; + int tmp; + // first index after noncommuting indices + int post_nc=0; + + for(i=0;i<monomial.length;i++){ + if(is_noncommuting(monomial.values[i], fields)){ + tmp=monomial.values[i]; + for(j=i;j>post_nc;j--){ + monomial.values[j]=monomial.values[j-1]; + } + monomial.values[post_nc]=tmp; + post_nc++; + } + } + + monomial_sort_nonc(monomial, post_nc, monomial.length-1, fields, sign); + + return(0); +} +// without noncommuting terms +int monomial_sort_nonc(Int_Array monomial, int begin, int end, Fields_Table fields, int* sign){ int i; int index; // the pivot: middle of the monomial @@ -812,8 +835,8 @@ int monomial_sort(Int_Array monomial, int begin, int end, Fields_Table fields, i exchange_monomial_terms(monomial, index, end, fields, sign); // recurse - monomial_sort(monomial, begin, index-1, fields, sign); - monomial_sort(monomial, index+1, end, fields, sign); + monomial_sort_nonc(monomial, begin, index-1, fields, sign); + monomial_sort_nonc(monomial, index+1, end, fields, sign); } return(0); } @@ -872,7 +895,30 @@ int exchange_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table // 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){ +// if the monomial contains noncommuting elements, put them at the beginning of the monomial +int monomial_sort_groups(Int_Array monomial, Fields_Table fields, Groups groups, int* sign){ + int i,j; + int tmp; + // first index after noncommuting indices + int post_nc=0; + + for(i=0;i<monomial.length;i++){ + if(is_noncommuting(monomial.values[i], fields)){ + tmp=monomial.values[i]; + for(j=post_nc;j<i;j++){ + monomial.values[j+1]=monomial.values[j]; + } + monomial.values[post_nc]=tmp; + post_nc++; + } + } + + monomial_sort_groups_nonc(monomial, post_nc, monomial.length-1, fields, groups, sign); + + return(0); +} +// without noncommuting terms +int monomial_sort_groups_nonc(Int_Array monomial, int begin, int end, Fields_Table fields, Groups groups, int* sign){ int i; int index; // the pivot: middle of the monomial @@ -897,8 +943,8 @@ int monomial_sort_groups(Int_Array monomial, int begin, int end, Fields_Table fi exchange_monomial_terms(monomial, index, end, fields, sign); // recurse - monomial_sort(monomial, begin, index-1, fields, sign); - monomial_sort(monomial, index+1, end, fields, sign); + monomial_sort_groups_nonc(monomial, begin, index-1, fields, groups, sign); + monomial_sort_groups_nonc(monomial, index+1, end, fields, groups, sign); } return(0); } diff --git a/src/polynomial.h b/src/polynomial.h index ec50227..6aff76d 100644 --- a/src/polynomial.h +++ b/src/polynomial.h @@ -99,14 +99,18 @@ int polynomial_sort(Polynomial* polynomial, int begin, int end); 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); +int monomial_sort(Int_Array monomial, Fields_Table fields, int* sign); +// without noncommuting terms +int monomial_sort_nonc(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); +int monomial_sort_groups(Int_Array monomial, Fields_Table fields, Groups groups, int* sign); +// without noncommuting terms +int monomial_sort_groups_nonc(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); diff --git a/src/types.h b/src/types.h index 0452ebc..6b029d5 100644 --- a/src/types.h +++ b/src/types.h @@ -154,6 +154,8 @@ typedef struct Fields_Table{ Symbols symbols; // list of anti-commuting variables (fields or symbols) Int_Array fermions; + // list of non-commuting variables (fields or symbols) + Int_Array noncommuting; } Fields_Table; // index labels |