source: src/libcfa/rational.c @ 630a82a

aaron-thesisarm-ehcleanup-dtorsctordeferred_resndemanglerenumforall-pointer-decaygc_noraiijacob/cs343-translationjenkins-sandboxmemorynew-astnew-ast-unique-exprnew-envno_listpersistent-indexerresolv-newstringwith_gc
Last change on this file since 630a82a was 630a82a, checked in by Peter A. Buhr <pabuhr@…>, 6 years ago

C99 compound literals now work, comment rational code, clean up hoisting AddVisit?

  • Property mode set to 100644
File size: 5.6 KB
Line 
1//                               -*- Mode: C -*-
2//
3// Cforall Version 1.0.0 Copyright (C) 2016 University of Waterloo
4//
5// The contents of this file are covered under the licence agreement in the
6// file "LICENCE" distributed with Cforall.
7//
8// rational.c --
9//
10// Author           : Peter A. Buhr
11// Created On       : Wed Apr  6 17:54:28 2016
12// Last Modified By : Peter A. Buhr
13// Last Modified On : Fri Apr  8 15:39:17 2016
14// Update Count     : 17
15//
16
17#include "rational"
18#include "fstream"
19#include "stdlib"
20
21extern "C" {
22#include <stdlib.h>                                                                             // exit
23} // extern
24
25
26// constants
27
28struct Rational 0 = {0, 1};
29struct Rational 1 = {1, 1};
30
31
32// helper
33
34// Calculate greatest common denominator of two numbers, the first of which may be negative. Used to reduce rationals.
35static long int gcd( long int a, long int b ) {
36    for ( ;; ) {                                                                                // Euclid's algorithm
37                long int r = a % b;
38          if ( r == 0 ) break;
39                a = b;
40                b = r;
41    } // for
42        return b;
43} // gcd
44
45static long int simplify( long int *n, long int *d ) {
46    if ( *d == 0 ) {
47                serr | "Invalid rational number construction: denominator cannot be equal to 0." | endl;
48                exit( EXIT_FAILURE );
49    } // exit
50    if ( *d < 0 ) { *d = -*d; *n = -*n; }                               // move sign to numerator
51    return gcd( abs( *n ), *d );                                                // simplify
52} // Rationalnumber::simplify
53
54
55// constructors
56
57Rational rational() {
58    return (Rational){ 0, 1 };
59//      Rational t = { 0, 1 };
60//      return t;
61} // rational
62
63Rational rational( long int n ) {
64    return (Rational){ n, 1 };
65//      Rational t = { n, 1 };
66//      return t;
67} // rational
68
69Rational rational( long int n, long int d ) {
70    long int t = simplify( &n, &d );                                    // simplify
71//    r = (Rational){ n / t, d / t };
72        Rational t = { n / t, d / t };
73        return t;
74} // rational
75
76
77// getter/setter for numerator/denominator
78
79long int numerator( Rational r ) {
80    return r.numerator;
81} // numerator
82
83long int numerator( Rational r, long int n ) {
84    long int prev = r.numerator;
85    long int t = gcd( abs( n ), r.denominator );                // simplify
86    r.numerator = n / t;
87    r.denominator = r.denominator / t;
88    return prev;
89} // numerator
90
91long int denominator( Rational r ) {
92    return r.denominator;
93} // denominator
94
95long int denominator( Rational r, long int d ) {
96    long int prev = r.denominator;
97    long int t = simplify( &r.numerator, &d );                  // simplify
98    r.numerator = r.numerator / t;
99    r.denominator = d / t;
100    return prev;
101} // denominator
102
103
104// comparison
105
106int ?==?( Rational l, Rational r ) {
107    return l.numerator * r.denominator == l.denominator * r.numerator;
108} // ?==?
109
110int ?!=?( Rational l, Rational r ) {
111    return ! ( l == r );
112} // ?!=?
113
114int ?<?( Rational l, Rational r ) {
115    return l.numerator * r.denominator < l.denominator * r.numerator;
116} // ?<?
117
118int ?<=?( Rational l, Rational r ) {
119    return l < r || l == r;
120} // ?<=?
121
122int ?>?( Rational l, Rational r ) {
123    return ! ( l <= r );
124} // ?>?
125
126int ?>=?( Rational l, Rational r ) {
127    return ! ( l < r );
128} // ?>=?
129
130
131// arithmetic
132
133Rational -?( Rational r ) {
134        Rational t = { -r.numerator, r.denominator };
135    return t;
136} // -?
137
138Rational ?+?( Rational l, Rational r ) {
139    if ( l.denominator == r.denominator ) {                             // special case
140                Rational t = { l.numerator + r.numerator, l.denominator };
141                return t;
142    } else {
143                Rational t = { l.numerator * r.denominator + l.denominator * r.numerator, l.denominator * r.denominator };
144                return t;
145    } // if
146} // ?+?
147
148Rational ?-?( Rational l, Rational r ) {
149    if ( l.denominator == r.denominator ) {                             // special case
150                Rational t = { l.numerator - r.numerator, l.denominator };
151                return t;
152    } else {
153                Rational t = { l.numerator * r.denominator - l.denominator * r.numerator, l.denominator * r.denominator };
154                return t;
155    } // if
156} // ?-?
157
158Rational ?*?( Rational l, Rational r ) {
159    Rational t = { l.numerator * r.numerator, l.denominator * r.denominator };
160        return t;
161} // ?*?
162
163Rational ?/?( Rational l, Rational r ) {
164    if ( r.numerator < 0 ) {
165                r.numerator = -r.numerator;
166                r.denominator = -r.denominator;
167        } // if
168        Rational t = { l.numerator * r.denominator, l.denominator * r.numerator };
169    return t;
170} // ?/?
171
172
173// conversion
174
175double widen( Rational r ) {
176        return (double)r.numerator / (double)r.denominator;
177} // widen
178
179// https://rosettacode.org/wiki/Convert_decimal_number_to_rational#C
180Rational narrow( double f, long int md ) {
181        if ( md <= 1 ) {                                                                        // maximum fractional digits too small?
182                Rational t = rational( f, 1 );                                  // truncate fraction
183                return t;
184        } // if
185
186        // continued fraction coefficients
187        long int a, h[3] = { 0, 1, 0 }, k[3] = { 1, 0, 0 };
188        long int x, d, n = 1;
189        int i, neg = 0;
190 
191        if ( f < 0 ) { neg = 1; f = -f; }
192        while ( f != floor( f ) ) { n <<= 1; f *= 2; }
193        d = f;
194 
195        // continued fraction and check denominator each step
196        for (i = 0; i < 64; i++) {
197                a = n ? d / n : 0;
198          if (i && !a) break;
199                x = d; d = n; n = x % n;
200                x = a;
201                if (k[1] * a + k[0] >= md) {
202                        x = (md - k[0]) / k[1];
203                  if ( ! (x * 2 >= a || k[1] >= md) ) break;
204                        i = 65;
205                } // if
206                h[2] = x * h[1] + h[0]; h[0] = h[1]; h[1] = h[2];
207                k[2] = x * k[1] + k[0]; k[0] = k[1]; k[1] = k[2];
208        } // for
209        Rational t = rational( neg ? -h[1] : h[1], k[1] );
210        return t;
211} // narrow
212
213
214// I/O
215
216forall( dtype istype | istream( istype ) )
217istype * ?|?( istype *is, Rational *r ) {
218        long int t;
219    is | &(r->numerator) | &(r->denominator);
220        t = simplify( &(r->numerator), &(r->denominator) );
221    r->numerator /= t;
222    r->denominator /= t;
223    return is;
224} // ?|?
225
226forall( dtype ostype | ostream( ostype ) )
227ostype * ?|?( ostype *os, Rational r ) {
228    return os | r.numerator | '/' | r.denominator;
229} // ?|?
230
231// Local Variables: //
232// tab-width: 4 //
233// End: //
Note: See TracBrowser for help on using the repository browser.