# source:src/libcfa/rational.c@9827c7ba

Last change on this file since 9827c7ba was 9827c7ba, checked in by Peter A. Buhr <pabuhr@…>, 8 years ago

forgot to include in previous commit

• Property mode set to `100644`
File size: 5.5 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
14// Update Count     : 18
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
60
61Rational rational( long int n ) {
62    return (Rational){ n, 1 };
63} // rational
64
65Rational rational( long int n, long int d ) {
66    long int t = simplify( &n, &d );                                    // simplify
67    return (Rational){ n / t, d / t };
68} // rational
69
70
71// getter/setter for numerator/denominator
72
73long int numerator( Rational r ) {
74    return r.numerator;
75} // numerator
76
77long int numerator( Rational r, long int n ) {
78    long int prev = r.numerator;
79    long int t = gcd( abs( n ), r.denominator );                // simplify
80    r.numerator = n / t;
81    r.denominator = r.denominator / t;
82    return prev;
83} // numerator
84
85long int denominator( Rational r ) {
86    return r.denominator;
87} // denominator
88
89long int denominator( Rational r, long int d ) {
90    long int prev = r.denominator;
91    long int t = simplify( &r.numerator, &d );                  // simplify
92    r.numerator = r.numerator / t;
93    r.denominator = d / t;
94    return prev;
95} // denominator
96
97
98// comparison
99
100int ?==?( Rational l, Rational r ) {
101    return l.numerator * r.denominator == l.denominator * r.numerator;
102} // ?==?
103
104int ?!=?( Rational l, Rational r ) {
105    return ! ( l == r );
106} // ?!=?
107
108int ?<?( Rational l, Rational r ) {
109    return l.numerator * r.denominator < l.denominator * r.numerator;
110} // ?<?
111
112int ?<=?( Rational l, Rational r ) {
113    return l < r || l == r;
114} // ?<=?
115
116int ?>?( Rational l, Rational r ) {
117    return ! ( l <= r );
118} // ?>?
119
120int ?>=?( Rational l, Rational r ) {
121    return ! ( l < r );
122} // ?>=?
123
124
125// arithmetic
126
127Rational -?( Rational r ) {
128        Rational t = { -r.numerator, r.denominator };
129    return t;
130} // -?
131
132Rational ?+?( Rational l, Rational r ) {
133    if ( l.denominator == r.denominator ) {                             // special case
134                Rational t = { l.numerator + r.numerator, l.denominator };
135                return t;
136    } else {
137                Rational t = { l.numerator * r.denominator + l.denominator * r.numerator, l.denominator * r.denominator };
138                return t;
139    } // if
140} // ?+?
141
142Rational ?-?( Rational l, Rational r ) {
143    if ( l.denominator == r.denominator ) {                             // special case
144                Rational t = { l.numerator - r.numerator, l.denominator };
145                return t;
146    } else {
147                Rational t = { l.numerator * r.denominator - l.denominator * r.numerator, l.denominator * r.denominator };
148                return t;
149    } // if
150} // ?-?
151
152Rational ?*?( Rational l, Rational r ) {
153    Rational t = { l.numerator * r.numerator, l.denominator * r.denominator };
154        return t;
155} // ?*?
156
157Rational ?/?( Rational l, Rational r ) {
158    if ( r.numerator < 0 ) {
159                r.numerator = -r.numerator;
160                r.denominator = -r.denominator;
161        } // if
162        Rational t = { l.numerator * r.denominator, l.denominator * r.numerator };
163    return t;
164} // ?/?
165
166
167// conversion
168
169double widen( Rational r ) {
170        return (double)r.numerator / (double)r.denominator;
171} // widen
172
173// https://rosettacode.org/wiki/Convert_decimal_number_to_rational#C
174Rational narrow( double f, long int md ) {
175        if ( md <= 1 ) {                                                                        // maximum fractional digits too small?
176                Rational t = rational( f, 1 );                                  // truncate fraction
177                return t;
178        } // if
179
180        // continued fraction coefficients
181        long int a, h[3] = { 0, 1, 0 }, k[3] = { 1, 0, 0 };
182        long int x, d, n = 1;
183        int i, neg = 0;
184
185        if ( f < 0 ) { neg = 1; f = -f; }
186        while ( f != floor( f ) ) { n <<= 1; f *= 2; }
187        d = f;
188
189        // continued fraction and check denominator each step
190        for (i = 0; i < 64; i++) {
191                a = n ? d / n : 0;
192          if (i && !a) break;
193                x = d; d = n; n = x % n;
194                x = a;
195                if (k[1] * a + k[0] >= md) {
196                        x = (md - k[0]) / k[1];
197                  if ( ! (x * 2 >= a || k[1] >= md) ) break;
198                        i = 65;
199                } // if
200                h[2] = x * h[1] + h[0]; h[0] = h[1]; h[1] = h[2];
201                k[2] = x * k[1] + k[0]; k[0] = k[1]; k[1] = k[2];
202        } // for
203        Rational t = rational( neg ? -h[1] : h[1], k[1] );
204        return t;
205} // narrow
206
207
208// I/O
209
210forall( dtype istype | istream( istype ) )
211istype * ?|?( istype *is, Rational *r ) {
212        long int t;
213    is | &(r->numerator) | &(r->denominator);
214        t = simplify( &(r->numerator), &(r->denominator) );
215    r->numerator /= t;
216    r->denominator /= t;
217    return is;
218} // ?|?
219
220forall( dtype ostype | ostream( ostype ) )
221ostype * ?|?( ostype *os, Rational r ) {
222    return os | r.numerator | '/' | r.denominator;
223} // ?|?
224
225// Local Variables: //
226// tab-width: 4 //
227// End: //
Note: See TracBrowser for help on using the repository browser.