# source:src/libcfa/rational.c@53ba273

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

switch from std=c99 to std=gnu99, update latex macros, refrat and extend user manual, update limits/iostream/fstream, add rational numbers

• Property mode set to `100644`
File size: 5.3 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     : 7
15//
16
17#include "rational"
18#include <fstream>
19#include <stdlib>
20
21extern "C" {
22#include <stdlib.h>                                                                             // exit
23} // extern
24
25struct Rational 0 = {0, 1};
26struct Rational 1 = {1, 1};
27
28// Calculate the greatest common denominator of two numbers, the first of which may be negative.  It is used to reduce
29// rationals.
30
31long int gcd( long int a, long int b ) {
32    for ( ;; ) {                                                                                // Euclid's algorithm
33                long int r = a % b;
34          if ( r == 0 ) break;
35                a = b;
36                b = r;
37    } // for
38        return b;
39} // gcd
40
41long int simplify( long int *n, long int *d ) {
42    if ( *d == 0 ) {
43                serr | "Invalid rational number construction: denominator cannot be equal to 0." | endl;
44                exit( EXIT_FAILURE );
45    } // exit
46    if ( *d < 0 ) { *d = -*d; *n = -*n; }                               // move sign to numerator
47    return gcd( abs( *n ), *d );                                                // simplify
48} // Rationalnumber::simplify
49
50Rational rational() {                                                                   // constructor
51//    r = (Rational){ 0, 1 };
52        Rational t = { 0, 1 };
53        return t;
54} // rational
55
56Rational rational( long int n ) {                                               // constructor
57//    r = (Rational){ n, 1 };
58        Rational t = { n, 1 };
59        return t;
60} // rational
61
62Rational rational( long int n, long int d ) {                   // constructor
63    long int t = simplify( &n, &d );                                    // simplify
64//    r = (Rational){ n / t, d / t };
65        Rational t = { n / t, d / t };
66        return t;
67} // rational
68
69long int numerator( Rational r ) {
70    return r.numerator;
71} // numerator
72
73long int numerator( Rational r, long int n ) {
74    long int prev = r.numerator;
75    long int t = gcd( abs( n ), r.denominator );                // simplify
76    r.numerator = n / t;
77    r.denominator = r.denominator / t;
78    return prev;
79} // numerator
80
81long int denominator( Rational r, long int d ) {
82    long int prev = r.denominator;
83    long int t = simplify( &r.numerator, &d );                  // simplify
84    r.numerator = r.numerator / t;
85    r.denominator = d / t;
86    return prev;
87} // denominator
88
89int ?==?( Rational l, Rational r ) {
90    return l.numerator * r.denominator == l.denominator * r.numerator;
91} // ?==?
92
93int ?!=?( Rational l, Rational r ) {
94    return ! ( l == r );
95} // ?!=?
96
97int ?<?( Rational l, Rational r ) {
98    return l.numerator * r.denominator < l.denominator * r.numerator;
99} // ?<?
100
101int ?<=?( Rational l, Rational r ) {
102    return l < r || l == r;
103} // ?<=?
104
105int ?>?( Rational l, Rational r ) {
106    return ! ( l <= r );
107} // ?>?
108
109int ?>=?( Rational l, Rational r ) {
110    return ! ( l < r );
111} // ?>=?
112
113Rational -?( Rational r ) {
114        Rational t = { -r.numerator, r.denominator };
115    return t;
116} // -?
117
118Rational ?+?( Rational l, Rational r ) {
119    if ( l.denominator == r.denominator ) {                             // special case
120                Rational t = { l.numerator + r.numerator, l.denominator };
121                return t;
122    } else {
123                Rational t = { l.numerator * r.denominator + l.denominator * r.numerator, l.denominator * r.denominator };
124                return t;
125    } // if
126} // ?+?
127
128Rational ?-?( Rational l, Rational r ) {
129    if ( l.denominator == r.denominator ) {                             // special case
130                Rational t = { l.numerator - r.numerator, l.denominator };
131                return t;
132    } else {
133                Rational t = { l.numerator * r.denominator - l.denominator * r.numerator, l.denominator * r.denominator };
134                return t;
135    } // if
136} // ?-?
137
138Rational ?*?( Rational l, Rational r ) {
139    Rational t = { l.numerator * r.numerator, l.denominator * r.denominator };
140        return t;
141} // ?*?
142
143Rational ?/?( Rational l, Rational r ) {
144    if ( r.numerator < 0 ) {
145                r.numerator = -r.numerator;
146                r.denominator = -r.denominator;
147        } // if
148        Rational t = { l.numerator * r.denominator, l.denominator * r.numerator };
149    return t;
150} // ?/?
151
152double widen( Rational r ) {
153        return (double)r.numerator / (double)r.denominator;
154} // widen
155
156// https://rosettacode.org/wiki/Convert_decimal_number_to_rational#C
157Rational narrow( double f, long int md ) {
158        if ( md <= 1 ) {                                                                        // maximum fractional digits too small?
159                Rational t = rational( f, 1 );                                  // truncate fraction
160                return t;
161        } // if
162
163        // continued fraction coefficients
164        long int a, h[3] = { 0, 1, 0 }, k[3] = { 1, 0, 0 };
165        long int x, d, n = 1;
166        int i, neg = 0;
167
168        if ( f < 0 ) { neg = 1; f = -f; }
169        while ( f != floor( f ) ) { n <<= 1; f *= 2; }
170        d = f;
171
172        // continued fraction and check denominator each step
173        for (i = 0; i < 64; i++) {
174                a = n ? d / n : 0;
175          if (i && !a) break;
176                x = d; d = n; n = x % n;
177                x = a;
178                if (k[1] * a + k[0] >= md) {
179                        x = (md - k[0]) / k[1];
180                  if ( ! (x * 2 >= a || k[1] >= md) ) break;
181                        i = 65;
182                } // if
183                h[2] = x * h[1] + h[0]; h[0] = h[1]; h[1] = h[2];
184                k[2] = x * k[1] + k[0]; k[0] = k[1]; k[1] = k[2];
185        } // for
186        Rational t = rational( neg ? -h[1] : h[1], k[1] );
187        return t;
188} // narrow
189
190ifstream * ?|?( ifstream *is, Rational *r ) {
191        long int t;
192    is | &(r->numerator) | &(r->denominator);
193        t = simplify( &(r->numerator), &(r->denominator) );
194    r->numerator /= t;
195    r->denominator /= t;
196    return is;
197} // ?|?
198
199ofstream * ?|?( ofstream *os, Rational r ) {
200    return os | r.numerator | '/' | r.denominator;
201} // ?|?
202
203// Local Variables: //
204// tab-width: 4 //
205// End: //
Note: See TracBrowser for help on using the repository browser.