Context Navigation

rational.cfa @ eae8b37

Last change on this file since eae8b37 was eae8b37, checked in by JiadaL <j82liang@…>, 2 weeks ago
Move enum.hfa/enum.cfa to prelude
Property mode set to `100644`
File size: 6.9 KB

Line
1	//
2	// Cforall Version 1.0.0 Copyright (C) 2016 University of Waterloo
3	//
4	// The contents of this file are covered under the licence agreement in the
5	// file "LICENCE" distributed with Cforall.
6	//
7	// rational.c --
8	//
9	// Author : Peter A. Buhr
10	// Created On : Wed Apr 6 17:54:28 2016
11	// Last Modified By : Peter A. Buhr
12	// Last Modified On : Wed Nov 27 18:06:43 2024
13	// Update Count : 208
14	//
15
16	#include "rational.hfa"
17	#include "fstream.hfa"
18	#include "stdlib.hfa"
19
20	#pragma GCC visibility push(default)
21
22	// Arithmetic, Relational
23	forall( T \| Simple(T) ) {
24	// helper routines
25	// Calculate greatest common denominator of two numbers, the first of which may be negative. Used to reduce
26	// rationals. alternative: https://en.wikipedia.org/wiki/Binary_GCD_algorithm
27	static T gcd( T a, T b ) {
28	for () { // Euclid's algorithm
29	T r = a % b;
30	if ( r == (T){0} ) break;
31	a = b;
32	b = r;
33	} // for
34	return b;
35	} // gcd
36
37	static T simplify( T & n, T & d ) {
38	if ( d == (T){0} ) {
39	abort \| "Invalid rational number construction: denominator cannot be equal to 0.";
40	} // exit
41	if ( d < (T){0} ) { d = -d; n = -n; } // move sign to numerator
42	return gcd( abs( n ), d ); // simplify
43	} // simplify
44	}
45
46	forall( T \| arithmetic( T ) ) {
47	// constructors
48
49	void ?{}( rational(T) & r, zero_t ) {
50	r{ (T){0}, (T){1} };
51	} // rational
52
53	void ?{}( rational(T) & r, one_t ) {
54	r{ (T){1}, (T){1} };
55	} // rational
56
57	void ?{}( rational(T) & r ) {
58	r{ (T){0}, (T){1} };
59	} // rational
60
61	void ?{}( rational(T) & r, T n ) {
62	r{ n, (T){1} };
63	} // rational
64
65	void ?{}( rational(T) & r, T n, T d ) {
66	T t = simplify( n, d ); // simplify
67	r.[numerator, denominator] = [n / t, d / t];
68	} // rational
69
70	// getter for numerator/denominator
71
72	T numerator( rational(T) r ) {
73	return r.numerator;
74	} // numerator
75
76	T denominator( rational(T) r ) {
77	return r.denominator;
78	} // denominator
79
80	[ T, T ] ?=?( & [ T, T ] dst, rational(T) src ) {
81	return dst = src.[ numerator, denominator ];
82	} // ?=?
83
84	// setter for numerator/denominator
85
86	T numerator( rational(T) r, T n ) {
87	T prev = r.numerator;
88	T t = gcd( abs( n ), r.denominator ); // simplify
89	r.[numerator, denominator] = [n / t, r.denominator / t];
90	return prev;
91	} // numerator
92
93	T denominator( rational(T) r, T d ) {
94	T prev = r.denominator;
95	T t = simplify( r.numerator, d ); // simplify
96	r.[numerator, denominator] = [r.numerator / t, d / t];
97	return prev;
98	} // denominator
99
100	// comparison
101
102	int ?==?( rational(T) l, rational(T) r ) {
103	return l.numerator * r.denominator == l.denominator * r.numerator;
104	} // ?==?
105
106	int ?!=?( rational(T) l, rational(T) r ) {
107	return ! ( l == r );
108	} // ?!=?
109
110	int ?!=?( rational(T) l, zero_t ) {
111	return ! ( l == (rational(T)){ 0 } );
112	} // ?!=?
113
114	int ?<?( rational(T) l, rational(T) r ) {
115	return l.numerator * r.denominator < l.denominator * r.numerator;
116	} // ?<?
117
118	int ?<=?( rational(T) l, rational(T) r ) {
119	return l.numerator * r.denominator <= l.denominator * r.numerator;
120	} // ?<=?
121
122	int ?>?( rational(T) l, rational(T) r ) {
123	return ! ( l <= r );
124	} // ?>?
125
126	int ?>=?( rational(T) l, rational(T) r ) {
127	return ! ( l < r );
128	} // ?>=?
129
130	// arithmetic
131
132	rational(T) +?( rational(T) r ) {
133	return (rational(T)){ r.numerator, r.denominator };
134	} // +?
135
136	rational(T) -?( rational(T) r ) {
137	return (rational(T)){ -r.numerator, r.denominator };
138	} // -?
139
140	rational(T) ?+?( rational(T) l, rational(T) r ) {
141	if ( l.denominator == r.denominator ) { // special case
142	return (rational(T)){ l.numerator + r.numerator, l.denominator };
143	} else {
144	return (rational(T)){ l.numerator * r.denominator + l.denominator * r.numerator, l.denominator * r.denominator };
145	} // if
146	} // ?+?
147
148	rational(T) ?+=?( rational(T) & l, rational(T) r ) {
149	l = l + r;
150	return l;
151	} // ?+?
152
153	rational(T) ?+=?( rational(T) & l, one_t ) {
154	l = l + (rational(T)){ 1 };
155	return l;
156	} // ?+?
157
158	rational(T) ?-?( rational(T) l, rational(T) r ) {
159	if ( l.denominator == r.denominator ) { // special case
160	return (rational(T)){ l.numerator - r.numerator, l.denominator };
161	} else {
162	return (rational(T)){ l.numerator * r.denominator - l.denominator * r.numerator, l.denominator * r.denominator };
163	} // if
164	} // ?-?
165
166	rational(T) ?-=?( rational(T) & l, rational(T) r ) {
167	l = l - r;
168	return l;
169	} // ?-?
170
171	rational(T) ?-=?( rational(T) & l, one_t ) {
172	l = l - (rational(T)){ 1 };
173	return l;
174	} // ?-?
175
176	rational(T) ?*?( rational(T) l, rational(T) r ) {
177	return (rational(T)){ l.numerator * r.numerator, l.denominator * r.denominator };
178	} // ?*?
179
180	rational(T) ?*=?( rational(T) & l, rational(T) r ) {
181	return l = l * r;
182	} // ?*?
183
184	rational(T) ?/?( rational(T) l, rational(T) r ) {
185	if ( r.numerator < (T){0} ) {
186	r.[numerator, denominator] = [-r.numerator, -r.denominator];
187	} // if
188	return (rational(T)){ l.numerator * r.denominator, l.denominator * r.numerator };
189	} // ?/?
190
191	rational(T) ?/=?( rational(T) & l, rational(T) r ) {
192	return l = l / r;
193	} // ?/?
194	} // distribution
195
196	// I/O
197
198	forall( T ) {
199	forall( istype & \| istream( istype ) \| { istype & ?\|?( istype &, T & ); } \| Simple(T) )
200	istype & ?\|?( istype & is, rational(T) & r ) {
201	is \| r.numerator \| r.denominator;
202	T t = simplify( r.numerator, r.denominator );
203	r.numerator /= t;
204	r.denominator /= t;
205	return is;
206	} // ?\|?
207
208	forall( ostype & \| ostream( ostype ) \| { ostype & ?\|?( ostype &, T ); } ) {
209	ostype & ?\|?( ostype & os, rational(T) r ) {
210	return os \| r.numerator \| '/' \| r.denominator;
211	} // ?\|?
212	OSTYPE_VOID_IMPL( os, rational(T) )
213	} // distribution
214	} // distribution
215
216	// Exponentiation
217
218	forall( T \| arithmetic( T ) \| { T ?\?( T, unsigned long ); } ) {
219	rational(T) ?\?( rational(T) x, long int y ) {
220	if ( y < 0 ) {
221	return (rational(T)){ x.denominator \ -y, x.numerator \ -y };
222	} else {
223	return (rational(T)){ x.numerator \ y, x.denominator \ y };
224	} // if
225	} // ?\?
226
227	rational(T) ?\=?( rational(T) & x, long int y ) {
228	return x = x \ y;
229	} // ?\?
230	} // distribution
231
232	// Conversion
233
234	forall( T \| arithmetic( T ) \| { double convert( T ); } )
235	double widen( rational(T) r ) {
236	return convert( r.numerator ) / convert( r.denominator );
237	} // widen
238
239	forall( T \| arithmetic( T ) \| { double convert( T ); T convert( double ); } )
240	rational(T) narrow( double f, T md ) {
241	// http://www.ics.uci.edu/~eppstein/numth/frap.c
242	if ( md <= (T){1} ) { // maximum fractional digits too small?
243	return (rational(T)){ convert( f ), (T){1}}; // truncate fraction
244	} // if
245
246	// continued fraction coefficients
247	T m00 = {1}, m11 = { 1 }, m01 = { 0 }, m10 = { 0 };
248	T ai, t;
249
250	// find terms until denom gets too big
251	for () {
252	ai = convert( f );
253	if ( ! (m10 * ai + m11 <= md) ) break;
254	t = m00 * ai + m01;
255	m01 = m00;
256	m00 = t;
257	t = m10 * ai + m11;
258	m11 = m10;
259	m10 = t;
260	double temp = convert( ai );
261	if ( f == temp ) break; // prevent division by zero
262	f = 1 / (f - temp);
263	if ( f > (double)0x7FFFFFFF ) break; // representation failure
264	} // for
265	return (rational(T)){ m00, m10 };
266	} // narrow
267
268	// Local Variables: //
269	// tab-width: 4 //
270	// End: //

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format