proc ds2; package math / overwrite = yes; method radians(double angle) returns double; return angle * constant('PI') / 180; end; method rational(double x, double maxden, in_out integer p, in_out integer q) returns double; /* ** FROM: https://www.ics.uci.edu/~eppstein/numth/frap.c ** FROM: https://stackoverflow.com/questions/95727/how-to-convert-floats-to-human-readable-fractions ** ** find rational approximation to given real number ** David Eppstein / UC Irvine / 8 Aug 1993 ** ** With corrections from Arno Formella, May 2008 ** ** Modified for Proc DS2, Richard DeVenezia, Jan 2020. ** ** usage: rational(r,d,p,q) ** x is real number to approx ** maxden is the maximum denominator allowed ** p is return for numerator ** q is return for denominator ** returns 0 if no problems ** ** based on the theory of continued fractions ** if x = a1 + 1/(a2 + 1/(a3 + 1/(a4 + ...))) ** then best approximation is found by truncating this series ** (with some adjustments in the last term). ** ** Note the fraction can be recovered as the first column of the matrix ** ( a1 1 ) ( a2 1 ) ( a3 1 ) ... ** ( 1 0 ) ( 1 0 ) ( 1 0 ) ** Instead of keeping the sequence of continued fraction terms, ** we just keep the last partial product of these matrices. */ declare integer m[0:1,0:1]; declare double startx e1 e2; declare integer ai t result p1 q1 p2 q2; startx = x; /* initialize matrix */ m[0,0] = 1; m[1,1] = 1; m[0,1] = 0; m[1,0] = 0; /* loop finding terms until denom gets too big */ do while (1); ai = x; if not ( m[1,0] * ai + m[1,1] < maxden ) then leave; t = m[0,0] * ai + m[0,1]; m[0,1] = m[0,0]; m[0,0] = t; t = m[1,0] * ai + m[1,1]; m[1,1] = m[1,0]; m[1,0] = t; if x = ai then leave; %* AF: division by zero; x = 1 / (x - ai); if x > 2147483647 /*x'7FFFFFFF'*/ then leave; %* AF: representation failure; end; /* now remaining x is between 0 and 1/ai */ /* approx as either 0 or 1/m where m is max that will fit in maxden */ /* first try zero */ p1 = m[0,0]; q1 = m[1,0]; e1 = startx - 1.0 * p1 / q1; /* now try other possibility */ ai = (maxden - m[1,1]) / m[1,0]; m[0,0] = m[0,0] * ai + m[0,1]; m[1,0] = m[1,0] * ai + m[1,1]; p2 = m[0,0]; q2 = m[1,0]; e2 = startx - 1.0 * p2 / q2; if abs(e1) <= abs(e2) then do; p = p1; q = q1; end; else do; p = p2; q = q2; end; return 0; end; endpackage; run; quit; * Example; proc ds2; data _null_; declare package math math(); declare double x; declare int p1 q1 p q; method run(); streaminit(12345); x = 0; do _n_ = 1 to 20; p1 = ceil(rand('uniform',9)); q1 = ceil(rand('uniform',9)); x + 1. * p1 / q1; math.rational (x, 10000, p, q); put 'add' p1 '/' q1 ' ' x=best16. 'is' p '/' q; end; end; enddata; run; quit;