Ian Jauslin
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorIan Jauslin <ian.jauslin@roma1.infn.it>2015-06-14 00:52:45 +0000
committerIan Jauslin <ian.jauslin@roma1.infn.it>2015-06-14 00:52:45 +0000
commitaa0f3ae2988d372b190b9bde2e75a6d17e744e93 (patch)
tree14482245c2fca27fcdad3078e97d0871352d52a7 /src/number.c
Initial commitv1.2
Diffstat (limited to 'src/number.c')
-rw-r--r--src/number.c551
1 files changed, 551 insertions, 0 deletions
diff --git a/src/number.c b/src/number.c
new file mode 100644
index 0000000..5d4cd18
--- /dev/null
+++ b/src/number.c
@@ -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);
+}