diff options
Diffstat (limited to 'src/number.c')
| -rw-r--r-- | src/number.c | 223 |
1 files changed, 175 insertions, 48 deletions
diff --git a/src/number.c b/src/number.c index e63427f..5bc87de 100644 --- a/src/number.c +++ b/src/number.c @@ -1,5 +1,5 @@ /* -Copyright 2015 Ian Jauslin +Copyright 2015-2022 Ian Jauslin Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -27,6 +27,7 @@ limitations under the License. #include "tools.h" #include "rational.h" #include "array.h" +#include "parse_file.h" // init int init_Number(Number* number, int memory){ @@ -271,50 +272,107 @@ Number number_Qprod_ret(Q q, Number x){ 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]; +// quotient of two numbers +// recursively simplify numerator/denominator until denominator only has one term +// the output is set to 'numerator' +// both numerator and denominator may be changed by this function +// this algorithm is not optimal in cases where denominator has several terms, in particular if their bases are large +// it is optimal if denominator only has one term +int number_quot_inplace(Number* numerator, Number* denominator){ + Number tmp; + int i,factor; + + switch((*denominator).length){ + // error + case 0: + fprintf(stderr,"error: attempting to invert 0\n"); + exit(-1); + break; + + // trivial case + case 1: + // trivial base + if((*denominator).base[0]==1){ + number_Qprod_chain(Q_inverse((*denominator).scalars[0]), numerator); } - 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; + // non-trivial base + else{ + // set tmp=1/denominator + tmp=number_one(); + tmp.base[0]=(*denominator).base[0]; + tmp.scalars[0].denominator=(*denominator).scalars[0].numerator*abs((*denominator).base[0]); + if((*denominator).base[0]>0){ + tmp.scalars[0].numerator=(*denominator).scalars[0].denominator; + } + else if((*denominator).base[0]<0){ + tmp.scalars[0].numerator=-(*denominator).scalars[0].denominator; + } + else{ + fprintf(stderr,"error: attempting to invert 0\n"); + exit(-1); + } + number_prod_chain(tmp, numerator); + } + break; + + default: + // find non-trivial basis + for(i=0;(*denominator).base[i]==1;i++){} + + // smallest prime factor of the base + if((*denominator).base[i]<0){ + factor=-1; + } + else if((*denominator).base[i]==2){ + factor=2; } else{ - fprintf(stderr,"error: attempting to invert 0\n"); - exit(-1); + // COMMENT: this is not optimal, but, provided the basis is not too large, this should not be a problem + for(factor=3;is_factor(factor, (*denominator).base[i])==0;factor=factor+2){} + } + + // tmp is set to ((terms that k does not divide)-(terms that k divides)) + // tmp will be multiplied to numerator and denominator + init_Number(&tmp,(*denominator).length); + + // find all terms whose base is a multiple of k + for(i=0;i<(*denominator).length;i++){ + if(is_factor(factor, (*denominator).base[i])){ + // add to tmp with a - sign + number_add_elem(quot(-(*denominator).scalars[i].numerator,(*denominator).scalars[i].denominator), (*denominator).base[i], &tmp); + } + else{ + // add to tmp + number_add_elem((*denominator).scalars[i], (*denominator).base[i], &tmp); + } } + + number_prod_chain(tmp, numerator); + number_prod_chain(tmp, denominator); + free_Number(tmp); + + // recurse + number_quot_inplace(numerator, denominator); } + 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 +// not inplace int number_quot(Number x1, Number x2, Number* output){ - Number inv; - number_inverse(x2, &inv); - number_prod(x1, inv, output); - free_Number(inv); + Number numerator, denominator; + number_cpy(x1, &numerator); + number_cpy(x2, &denominator); + number_quot_inplace(&numerator, &denominator); + *output=numerator; + free_Number(denominator); return(0); } -int number_quot_chain(Number x1, Number* inout){ - number_inverse_inplace(inout); - number_prod_chain(x1, inout); +int number_quot_chain(Number* inout, Number x2){ + Number tmp; + number_quot(*inout,x2,&tmp); + free_Number(*inout); + *inout=tmp; return(0); } Number number_quot_ret(Number x1, Number x2){ @@ -419,7 +477,7 @@ int number_print(Number number){ Char_Array buffer; init_Char_Array(&buffer,5*number.length); number_sprint(number, &buffer); - printf("%s",buffer.str); + printf("%s",char_array_to_str_noinit(&buffer)); return(0); } @@ -434,18 +492,54 @@ int str_to_Number(char* str, Number* number){ char* buffer_ptr=buffer; Q num; int base; - // whether there are parentheses in the string - int exist_parenthesis=0; + int ret; + int j; + char* aux_str; + int aux_free=0; init_Number(number, NUMBER_SIZE); + // check whether the string is blank (return 0 in that case) + for(ptr=str;*ptr!='\0';ptr++){ + if(*ptr!=' ' && *ptr!='\n'){ + break; + } + } + // blank string + if(*ptr=='\0'){ + number_add_elem(quot(0,1), 1, number); + free(buffer); + return(0); + } + // init num and base - // init to 0 so that if str is empty, then the number is set to 0 - num=quot(0,1); + num=quot(1,1); base=1; + // check whether the str only contains a rational number, and add parentheses + // keep rtack of the length of str + for(j=0,ptr=str;*ptr!='\0';j++,ptr++){ + if((*ptr<'0' || *ptr>'9') && *ptr!='-' && *ptr!='/'){ + break; + } + } + // only rational + if(*ptr=='\0'){ + aux_str=calloc(j+3,sizeof(char)); + aux_str[0]='('; + for(j=0,ptr=str;*ptr!='\0';ptr++,j++){ + aux_str[j+1]=*ptr; + } + aux_str[j+1]=')'; + aux_str[j+2]='\0'; + aux_free=1; + } + else{ + aux_str=str; + } + mode=PP_NULL_MODE; - for(ptr=str;*ptr!='\0';ptr++){ + for(ptr=aux_str;*ptr!='\0';ptr++){ switch(*ptr){ // read number case '(': @@ -453,7 +547,10 @@ int str_to_Number(char* str, Number* number){ // init base base=1; mode=PP_NUM_MODE; - exist_parenthesis=1; + } + else{ + fprintf(stderr,"syntax error: misplaced '(' in number '%s'\n",str); + exit(-1); } break; case ')': @@ -463,6 +560,10 @@ int str_to_Number(char* str, Number* number){ *buffer_ptr='\0'; mode=PP_NULL_MODE; } + else{ + fprintf(stderr,"syntax error: mismatched ')' in number '%s'\n",str); + exit(-1); + } break; // read sqrt @@ -474,16 +575,26 @@ int str_to_Number(char* str, Number* number){ 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; + else{ + fprintf(stderr,"syntax error: misplaced '{' in number '%s'\n",str); + exit(-1); + } break; case '}': if(mode==PP_SQRT_MODE){ - sscanf(buffer,"%d",&base); + ret=read_int(buffer,&base); + if(ret<0){ + fprintf(stderr,"syntax error: number base should be an integer, got '%s' in '%s'",buffer,str); + exit(-1); + } buffer_ptr=buffer; *buffer_ptr='\0'; mode=PP_NULL_MODE; } + else{ + fprintf(stderr,"syntax error: mismatched '}' in number '%s'\n",str); + exit(-1); + } break; // write num @@ -494,24 +605,40 @@ int str_to_Number(char* str, Number* number){ num=quot(0,1); base=1; } + else{ + fprintf(stderr,"syntax error: misplaced '+' in number '%s'\n",str); + exit(-1); + } break; + // ignore 's', ' ' and '\n' + case 's':break; + case ' ':break; + case '\n':break; default: if(mode!=PP_NULL_MODE){ buffer_ptr=str_addchar(buffer_ptr,*ptr); } + else{ + fprintf(stderr,"syntax error: unrecognized character '%c' in number '%s'\n",*ptr,str); + exit(-1); + } break; } } - // last step if(mode==PP_NULL_MODE){ - if(exist_parenthesis==0){ - str_to_Q(str, &num); - } number_add_elem(num, base, number); } + else{ + fprintf(stderr,"syntax error: mismatched '(' in number '%s'\n",str); + exit(-1); + } + + if(aux_free==1){ + free(aux_str); + } free(buffer); return(0); } |