Initial commitv1.2

1 files changed, 196 insertions, 0 deletions

diff --git a/src/rational_float.c b/src/rational_float.c
new file mode 100644
index 0000000..eebc4f4
--- /dev/null
+++ b/src/rational_float.c
@@ -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