Arithmetic/Rational
You are encouraged to solve this task according to the task description, using any language you may know.
The objective of this task is to create a reasonably complete implementation of rational arithmetic in the particular language using the idioms of the language.
For example: Define an new type called frac with binary operator "//" of two integers that returns a structure made up of the numerator and the denominator (as per a rational number).
Further define the appropriate rational unary operators abs and '-', with the binary operators for addition '+', subtraction '-', multiplication '×', division '/', integer division '÷', modulo division, the comparison operators (e.g. '<', '≤', '>', & '≥') and equality operators (e.g. '=' & '≠').
Define standard coercion operators for casting int to frac etc.
If space allows, define standard increment and decrement operators (e.g. '+:=' & '-:=' etc.).
Finally test the operators: Use the new type frac to find all perfect numbers less then 219 by summing the reciprocal of the factors.
See also
Ada
The generic package specification: <lang ada>generic
type Number is range <>;
package Generic_Rational is
type Rational is private; function "abs" (A : Rational) return Rational; function "+" (A : Rational) return Rational; function "-" (A : Rational) return Rational; function Inverse (A : Rational) return Rational; function "+" (A : Rational; B : Rational) return Rational; function "+" (A : Rational; B : Number ) return Rational; function "+" (A : Number; B : Rational) return Rational;
function "-" (A : Rational; B : Rational) return Rational; function "-" (A : Rational; B : Number ) return Rational; function "-" (A : Number; B : Rational) return Rational;
function "*" (A : Rational; B : Rational) return Rational; function "*" (A : Rational; B : Number ) return Rational; function "*" (A : Number; B : Rational) return Rational;
function "/" (A : Rational; B : Rational) return Rational; function "/" (A : Rational; B : Number ) return Rational; function "/" (A : Number; B : Rational) return Rational; function "/" (A : Number; B : Number) return Rational; function ">" (A : Rational; B : Rational) return Boolean; function ">" (A : Number; B : Rational) return Boolean; function ">" (A : Rational; B : Number) return Boolean;
function "<" (A : Rational; B : Rational) return Boolean; function "<" (A : Number; B : Rational) return Boolean; function "<" (A : Rational; B : Number) return Boolean;
function ">=" (A : Rational; B : Rational) return Boolean; function ">=" (A : Number; B : Rational) return Boolean; function ">=" (A : Rational; B : Number) return Boolean;
function "<=" (A : Rational; B : Rational) return Boolean; function "<=" (A : Number; B : Rational) return Boolean; function "<=" (A : Rational; B : Number) return Boolean;
function "=" (A : Number; B : Rational) return Boolean; function "=" (A : Rational; B : Number) return Boolean;
function Numerator (A : Rational) return Number;
function Denominator (A : Rational) return Number;
Zero : constant Rational;
One : constant Rational;
private
type Rational is record
Numerator : Number;
Denominator : Number;
end record;
Zero : constant Rational := (0, 1); One : constant Rational := (1, 1);
end Generic_Rational;</lang> The package can be instantiated with any integer type. It provides rational numbers represented by a numerator and denominator cleaned from the common divisors. Mixed arithmetic of the base integer type and the rational type is supported. Division to zero raises Constraint_Error. The implementation of the specification above is as follows: <lang ada>package body Generic_Rational is
function GCD (A, B : Number) return Number is
begin
if A = 0 then
return B;
end if;
if B = 0 then
return A;
end if;
if A > B then
return GCD (B, A mod B);
else
return GCD (A, B mod A);
end if;
end GCD;
function Inverse (A : Rational) return Rational is
begin
if A.Numerator > 0 then
return (A.Denominator, A.Numerator);
elsif A.Numerator < 0 then
return (-A.Denominator, -A.Numerator);
else
raise Constraint_Error;
end if;
end Inverse;
function "abs" (A : Rational) return Rational is
begin
return (abs A.Numerator, A.Denominator);
end "abs";
function "+" (A : Rational) return Rational is
begin
return A;
end "+";
function "-" (A : Rational) return Rational is
begin
return (-A.Numerator, A.Denominator);
end "-";
function "+" (A : Rational; B : Rational) return Rational is
Common : constant Number := GCD (A.Denominator, B.Denominator);
A_Denominator : constant Number := A.Denominator / Common;
B_Denominator : constant Number := B.Denominator / Common;
begin
return (A.Numerator * B_Denominator + B.Numerator * A_Denominator) /
(A_Denominator * B.Denominator);
end "+";
function "+" (A : Rational; B : Number) return Rational is
begin
return (A.Numerator + B * A.Denominator) / A.Denominator;
end "+";
function "+" (A : Number; B : Rational) return Rational is
begin
return B + A;
end "+";
function "-" (A : Rational; B : Rational) return Rational is
begin
return A + (-B);
end "-";
function "-" (A : Rational; B : Number) return Rational is
begin
return A + (-B);
end "-";
function "-" (A : Number; B : Rational) return Rational is
begin
return A + (-B);
end "-";
function "*" (A : Rational; B : Rational) return Rational is
begin
return (A.Numerator * B.Numerator) / (A.Denominator * B.Denominator);
end "*";
function "*" (A : Rational; B : Number) return Rational is
Common : constant Number := GCD (A.Denominator, abs B);
begin
return (A.Numerator * B / Common, A.Denominator / Common);
end "*";
function "*" (A : Number; B : Rational) return Rational is
begin
return B * A;
end "*";
function "/" (A : Rational; B : Rational) return Rational is
begin
return A * Inverse (B);
end "/";
function "/" (A : Rational; B : Number) return Rational is
Common : constant Number := GCD (abs A.Numerator, abs B);
begin
if B > 0 then
return (A.Numerator / Common, A.Denominator * (B / Common));
else
return ((-A.Numerator) / Common, A.Denominator * ((-B) / Common));
end if;
end "/";
function "/" (A : Number; B : Rational) return Rational is
begin
return Inverse (B) * A;
end "/";
function "/" (A : Number; B : Number) return Rational is
Common : constant Number := GCD (abs A, abs B);
begin
if B = 0 then
raise Constraint_Error;
elsif A = 0 then
return (0, 1);
elsif A > 0 xor B > 0 then
return (-(abs A / Common), abs B / Common);
else
return (abs A / Common, abs B / Common);
end if;
end "/";
function ">" (A, B : Rational) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator > 0;
end ">";
function ">" (A : Number; B : Rational) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator > 0;
end ">";
function ">" (A : Rational; B : Number) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator > 0;
end ">";
function "<" (A, B : Rational) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator < 0;
end "<";
function "<" (A : Number; B : Rational) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator < 0;
end "<";
function "<" (A : Rational; B : Number) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator < 0;
end "<";
function ">=" (A, B : Rational) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator >= 0;
end ">=";
function ">=" (A : Number; B : Rational) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator >= 0;
end ">=";
function ">=" (A : Rational; B : Number) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator >= 0;
end ">=";
function "<=" (A, B : Rational) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator <= 0;
end "<=";
function "<=" (A : Number; B : Rational) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator <= 0;
end "<=";
function "<=" (A : Rational; B : Number) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator <= 0;
end "<=";
function "=" (A : Number; B : Rational) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator = 0;
end "=";
function "=" (A : Rational; B : Number) return Boolean is
Diff : constant Rational := A - B;
begin
return Diff.Numerator = 0;
end "=";
function Numerator (A : Rational) return Number is
begin
return A.Numerator;
end Numerator;
function Denominator (A : Rational) return Number is
begin
return A.Denominator;
end Denominator;
end Generic_Rational;</lang> The implementation uses solution of the greatest common divisor task. Here is the implementation of the test: <lang ada>with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions; with Ada.Text_IO; use Ada.Text_IO; with Generic_Rational;
procedure Test_Rational is
package Integer_Rational is new Generic_Rational (Integer); use Integer_Rational;
begin
for Candidate in 2..2**15 loop
declare
Sum : Rational := 1 / Candidate;
begin
for Divisor in 2..Integer (Sqrt (Float (Candidate))) loop
if Candidate mod Divisor = 0 then -- Factor is a divisor of Candidate
Sum := Sum + One / Divisor + Rational'(Divisor / Candidate);
end if;
end loop;
if Sum = 1 then
Put_Line (Integer'Image (Candidate) & " is perfect");
end if;
end;
end loop;
end Test_Rational;</lang> The perfect numbers are searched by summing of the reciprocal of each of the divisors of a candidate except 1. This sum must be 1 for a perfect number. Sample output:
6 is perfect 28 is perfect 496 is perfect 8128 is perfect
ALGOL 68
MODE FRAC = STRUCT( INT num #erator#, den #ominator#);
FORMAT frac repr = $g(-0)"//"g(-0)$;
PROC gcd = (INT a, b) INT: # greatest common divisor #
(a = 0 | b |: b = 0 | a |: ABS a > ABS b | gcd(b, a MOD b) | gcd(a, b MOD a));
PROC lcm = (INT a, b)INT: # least common multiple #
a OVER gcd(a, b) * b;
PROC raise not implemented error = ([]STRING args)VOID: (
put(stand error, ("Not implemented error: ",args, newline));
stop
);
PRIO // = 9; # higher then the ** operator #
OP // = (INT num, den)FRAC: ( # initialise and normalise #
INT common = gcd(num, den);
IF den < 0 THEN
( -num OVER common, -den OVER common)
ELSE
( num OVER common, den OVER common)
FI
);
OP + = (FRAC a, b)FRAC: (
INT common = lcm(den OF a, den OF b);
FRAC result := ( common OVER den OF a * num OF a + common OVER den OF b * num OF b, common );
num OF result//den OF result
);
OP - = (FRAC a, b)FRAC: a + -b,
* = (FRAC a, b)FRAC: (
INT num = num OF a * num OF b,
den = den OF a * den OF b;
INT common = gcd(num, den);
(num OVER common) // (den OVER common)
);
OP / = (FRAC a, b)FRAC: a * FRAC(den OF b, num OF b),# real division #
% = (FRAC a, b)INT: ENTIER (a / b), # integer divison #
%* = (FRAC a, b)FRAC: a/b - FRACINIT ENTIER (a/b), # modulo division #
** = (FRAC a, INT exponent)FRAC:
IF exponent >= 0 THEN
(num OF a ** exponent, den OF a ** exponent )
ELSE
(den OF a ** exponent, num OF a ** exponent )
FI;
OP REALINIT = (FRAC frac)REAL: num OF frac / den OF frac,
FRACINIT = (INT num)FRAC: num // 1,
FRACINIT = (REAL num)FRAC: (
# express real number as a fraction # # a future execise! #
raise not implemented error(("Convert a REAL to a FRAC","!"));
SKIP
);
OP < = (FRAC a, b)BOOL: num OF (a - b) < 0,
> = (FRAC a, b)BOOL: num OF (a - b) > 0,
<= = (FRAC a, b)BOOL: NOT ( a > b ),
>= = (FRAC a, b)BOOL: NOT ( a < b ),
= = (FRAC a, b)BOOL: (num OF a, den OF a) = (num OF b, den OF b),
/= = (FRAC a, b)BOOL: (num OF a, den OF a) /= (num OF b, den OF b);
# Unary operators #
OP - = (FRAC frac)FRAC: (-num OF frac, den OF frac),
ABS = (FRAC frac)FRAC: (ABS num OF frac, ABS den OF frac),
ENTIER = (FRAC frac)INT: (num OF frac OVER den OF frac) * den OF frac;
COMMENT Operators for extended characters set, and increment/decrement "commented out" to save space. COMMENT
Example: searching for Perfect Numbers.
FRAC sum:= FRACINIT 0;
FORMAT perfect = $b(" perfect!","")$;
FOR i FROM 2 TO 2**19 DO
INT candidate := i;
FRAC sum := 1 // candidate;
REAL real sum := 1 / candidate;
FOR factor FROM 2 TO ENTIER sqrt(candidate) DO
IF candidate MOD factor = 0 THEN
sum := sum + 1 // factor + 1 // ( candidate OVER factor);
real sum +:= 1 / factor + 1 / ( candidate OVER factor)
FI
OD;
IF den OF sum = 1 THEN
printf(($"Sum of reciprocal factors of "g(-0)" = "g(-0)" exactly, about "g(0,real width) f(perfect)l$,
candidate, ENTIER sum, real sum, ENTIER sum = 1))
FI
OD
Output:
Sum of reciprocal factors of 6 = 1 exactly, about 1.0000000000000000000000000001 perfect! Sum of reciprocal factors of 28 = 1 exactly, about 1.0000000000000000000000000001 perfect! Sum of reciprocal factors of 120 = 2 exactly, about 2.0000000000000000000000000002 Sum of reciprocal factors of 496 = 1 exactly, about 1.0000000000000000000000000001 perfect! Sum of reciprocal factors of 672 = 2 exactly, about 2.0000000000000000000000000001 Sum of reciprocal factors of 8128 = 1 exactly, about 1.0000000000000000000000000001 perfect! Sum of reciprocal factors of 30240 = 3 exactly, about 3.0000000000000000000000000002 Sum of reciprocal factors of 32760 = 3 exactly, about 3.0000000000000000000000000003 Sum of reciprocal factors of 523776 = 2 exactly, about 2.0000000000000000000000000005
C
C doesn't support classes, so a series of functions is provided that implments the arithmetic of a rational class. <lang c>//#include <math.h>
- include <stdio.h>
- include <stdlib.h>
- include <string.h>
// a structure to contain the numerator and denominator values typedef struct sRational {
int numerator, denominator;
} *Rational, sRational;
int gcd( int a, int b) {
int g;
if (a<b) {
g = a; a = b; b = g;
}
g = a % b;
while (g) {
a = b; b = g;
g = a % b;
}
return abs(b);
}
int lcm( int a, int b) {
return (a/gcd(a,b)*b);
}
Rational NewRational( int n, int d) {
Rational r = (Rational)malloc(sizeof(sRational));
int ndgcd = gcd(n,d);
if (r) {
if (n>0) {
r->numerator = n/ndgcd;
r->denominator = d/ndgcd;
}
else {
r->numerator = -n/ndgcd;
r->denominator = -d/ndgcd;
}
}
if ( d == 0) {
printf("divide by zer0 error\n");
exit(1);
}
return r;
}
- define Ratl_Delete(r) \
{ if(r) free(r); \
r = NULL; }
Rational Ratl_Add( Rational l, Rational r, Rational rzlt) {
int denom; denom= lcm(l->denominator, r->denominator); rzlt->numerator = denom/l->denominator *l->numerator + denom/r->denominator *r->numerator ; rzlt->denominator = denom; return rzlt;
}
Rational Ratl_Negate( Rational l, Rational rzlt) {
rzlt->numerator = - l->numerator; rzlt->denominator = l->denominator; return rzlt;
}
Rational Ratl_Subtract(Rational l, Rational r, Rational rzlt) {
return Ratl_Add(l,Ratl_Negate(r, rzlt), rzlt);
}
Rational Ratl_Multiply(Rational l, Rational r, Rational rzlt) {
int g1 = gcd(l->numerator, r->denominator); int g2 = gcd(r->numerator, l->denominator); rzlt->numerator = l->numerator / g1 * r->numerator / g2; rzlt->denominator = l->denominator / g2 * r->denominator / g1; return rzlt;
}
Rational Ratl_Inverse(Rational l, Rational rzlt) {
if (l->numerator == 0) {
printf("divide by zer0 error\n");
exit(1);
}
else if (l->numerator > 0) {
rzlt->numerator = l->denominator;
rzlt->denominator = l->numerator;
}
else {
rzlt->numerator = -l->denominator;
rzlt->denominator = -l->numerator;
}
return rzlt;
}
Rational Ratl_Divide(Rational l, Rational r, Rational rzlt) {
return Ratl_Multiply(l, Ratl_Inverse(r, rzlt), rzlt);
}
int ipow(int base, int power ) {
int v = base, v2 = 1;
if (power < 0) return 0; // shouldn't happen
while (power > 0) {
if (power & 1)
v2 *=v;
v = v*v;
power >>= 1;
}
return v2;
}
Rational Ratl_Pow( Rational l, int power, Rational rzlt) {
if (power >= 0) {
rzlt->numerator = ipow(l->numerator, power);
rzlt->denominator = ipow(l->denominator, power);
}
else {
rzlt->numerator = ipow(l->denominator, -power);
rzlt->denominator = ipow(l->numerator, -power);
}
return rzlt;
}
int Ratl_Compare(Rational l, Rational r) {
int sign = (l->numerator > 0)? 1 : -1; sRational comp; Rational pcomp;
if ( 0 >= l->numerator * r->numerator) // if opposite signs or one is zero
return (l->numerator - r->numerator);
pcomp = Ratl_Divide(l, r, &comp);
return (pcomp->numerator - pcomp->denominator)*sign;
}
typedef enum {
LT, LE, EQ, GE, GT, NE } CompOp;
// boolean comparisons int Ratl_Cpr( Rational l, CompOp compOp, Rational r) {
int v;
int r1 = Ratl_Compare(l, r);
switch(compOp) {
case LT: v = (r1 <0)? 1 : 0; break;
case LE: v = (r1<=0)? 1 : 0; break;
case GT: v = (r1 >0)? 1 : 0; break;
case GE: v = (r1>=0)? 1 : 0; break;
case EQ: v = (r1==0)? 1 : 0; break;
case NE: v = (r1!=0)? 1 : 0; break;
}
return v;
}
double Ratl_2Real(Rational l)
{
return l->numerator *1.0/l->denominator;
}
int Ratl_2Int(Rational l) {
return l->numerator /l->denominator;
}
int Ratl_2Proper(Rational l) {
int ipart = l->numerator/l->denominator; l->numerator %= l->denominator; return ipart;
}
char *Ratl_2String(Rational l, char *buf, int blen) {
char ibuf[40];
sprintf(ibuf,"%d/%d", l->numerator, l->denominator );
if (buf==NULL) return buf;
if ((int)strlen(ibuf) < blen)
strcpy(buf, ibuf);
else {
strncpy(buf, ibuf, blen-4);
strcat(buf, "...");
}
return buf;
}
Rational Ratl_Abs(Rational l, Rational rzlt) {
rzlt->numerator = abs(l->numerator); rzlt->denominator = abs(l->denominator); // should already be nonnegative return rzlt;
}</lang> Example of use: <lang c>int main(int argc, char *argv[]) {
char ratstr[32], rs1[32],rs2[32]; sRational rtemp1, rtemp2; Rational rz; Rational r1 = NewRational( 5,-7 ); Rational r2 = NewRational( 4,5 ); Rational r3 = NewRational( 3,4);
printf("r3 = %s\n", Ratl_2String(r3, ratstr,32));
rz = Ratl_Multiply( r1,r2, &rtemp1);
printf("%s = %s * %s\n", Ratl_2String(rz, ratstr,32),
Ratl_2String(r1, rs1,32), Ratl_2String(r2, rs2, 32));
rz = Ratl_Divide( r1,r3, &rtemp2);
printf("%s = %s / %s\n", Ratl_2String(rz, ratstr,32),
Ratl_2String(r1, rs1,32), Ratl_2String(r3, rs2, 32));
rz = Ratl_Add( r2,r3, &rtemp1);
printf("%s = %s + %s\n", Ratl_2String(rz, ratstr,32),
Ratl_2String(r2, rs1,32), Ratl_2String(r3, rs2, 32));
rz = Ratl_Subtract( r2,r3, &rtemp1);
printf("%s = %s - %s\n", Ratl_2String(rz, ratstr,32),
Ratl_2String(r2, rs1,32), Ratl_2String(r3, rs2, 32));
printf("%d = %s > %s\n", Ratl_Cpr( r2, GT, &rtemp2),
Ratl_2String(r2, rs1,32), Ratl_2String(&rtemp2, rs2, 32));
rz = Ratl_Pow( r2,-3, &rtemp1);
printf("%s = %s ^ -3\n", Ratl_2String(rz, ratstr,32),
Ratl_2String(r2, rs1,32));
printf("%s = %f\n", Ratl_2String( &rtemp2, ratstr, 32), Ratl_2Real( &rtemp2));
return 0;
}</lang>
Common Lisp
Common Lisp has rational numbers built-in and integrated with all other number types. Common Lisp's number system is not extensible so reimplementing rational arithmetic would require all-new operator names.
<lang lisp>(loop for candidate from 2 below (expt 2 19)
for sum = (+ (/ candidate)
(loop for factor from 2 to (isqrt candidate)
when (zerop (mod candidate factor))
sum (+ (/ factor) (/ (floor candidate factor)))))
when (= sum 1)
collect candidate)</lang>
Fortran
This module currently implements only things needed for the perfect test, and few more (assignment from integer and real, the < operator...). For the GCD as function, see GCD
<lang fortran>module FractionModule
implicit none
type rational
integer :: numerator, denominator
end type rational
interface assignment (=)
module procedure assign_fraction_int, assign_fraction_float
end interface
interface operator (+)
module procedure add_fraction !, &
!float_add_fraction, int_add_fraction, &
!fraction_add_int, fraction_add_float
end interface
interface operator (<)
module procedure frac_lt_frac
end interface
interface int
module procedure frac_int
end interface int
interface real
module procedure frac_real
end interface real
real, parameter :: floatprec = 10e6 private floatprec, lcm , gcd
contains
! gcd can be defined here as private,or must be ! "loaded" by a proper use statement
function lcm(a, b) integer :: lcm integer, intent(in) :: a, b lcm = a / gcd(a, b) * b end function lcm
! int (which truncates, does not round) function frac_int(aratio) integer :: frac_int type(rational), intent(in) :: aratio frac_int = aratio%numerator / aratio%denominator end function frac_int
! get the real value of the fraction function frac_real(f) real :: frac_real type(rational), intent(in) :: f
frac_real = real(f%numerator) / real(f%denominator) end function frac_real
! frac = frac + frac
function add_fraction(f1, f2)
type(rational) :: add_fraction
type(rational), intent(in) :: f1, f2
integer :: common
common = lcm(f1%denominator, f2%denominator)
add_fraction = rational(common / f1%denominator * f1%numerator + &
common / f2%denominator * f2%numerator, common)
end function add_fraction
! assignment: frac = integer subroutine assign_fraction_int(ass, i) type(rational), intent(out), volatile :: ass integer, intent(in) :: i
ass = rational(i, 1) end subroutine assign_fraction_int
! assignment: frac = float(real) subroutine assign_fraction_float(ass, f) type(rational), intent(out), volatile :: ass real, intent(in) :: f
ass = rational(int(f*floatprec), floatprec) end subroutine assign_fraction_float
! simplify (must be called explicitly) subroutine simplify_fraction(fr) type(rational), intent(inout) :: fr integer :: d d = gcd(fr%numerator, fr%denominator) fr = rational(fr%numerator / d, fr%denominator / d) end subroutine simplify_fraction
! relational frac < frac
function frac_lt_frac(f1, f2)
logical :: frac_lt_frac
type(rational), intent(in) :: f1, f2
if ( real(f1) < real(f2) ) then
frac_lt_frac = .TRUE.
else
frac_lt_frac = .FALSE.
end if
end function frac_lt_frac
end module FractionModule</lang>
Testing this core code:
<lang fortran>program RatTesting
use fractionmodule
implicit none
integer :: i, candidate, factor type(rational) :: summa character(10) :: pstr
do i = 2, 2**19
candidate = i
summa = rational(1, candidate)
factor = 2
do while ( factor < sqrt(real(candidate)) )
if ( mod(candidate, factor) == 0 ) then
summa = summa + rational(1, factor) + &
rational(1, candidate/factor)
call simplify_fraction(summa)
end if
factor = factor + 1
end do
if ( summa%denominator == 1 ) then
if ( int(summa) == 1 ) then
pstr = 'perfect!'
else
pstr =
end if
write (*, '(A,1X,I7, = ,I1, 1X, A8)') 'Sum of recipr. factors of', &
candidate, int(summa), pstr
end if
end do
end program RatTesting</lang>
Haskell
Haskell provides a Rational type, which is really an alias for Ratio Integer (Ratio being a polymorphic type implementing rational numbers for any Integral type of numerators and denominators). The fraction is constructed using the % operator.
<lang haskell>import Data.Ratio
-- simply prints all the perfect numbers main = mapM_ print [candidate
| candidate <- [2 .. 2^19],
getSum candidate == 1]
where getSum candidate = 1 % candidate +
sum [1 % factor + 1 % (candidate `div` factor)
| factor <- [2 .. floor(sqrt(fromIntegral(candidate)))],
candidate `mod` factor == 0]</lang>
For a sample implementation of Ratio, see the Haskell 98 Report.
J
J implements rational numbers:
<lang j> 3r4*2r5 3r10</lang>
That said, note that J's floating point numbers work just fine for the stated problem: <lang j> is_perfect_rational=: 2 = (1 + i.) +/@:%@([ #~ 0 = |) ]</lang>
faster version (but the problem, as stated, is still tremendously inefficient): <lang j>factors=: */&>@{@((^ i.@>:)&.>/)@q:~&__ is_perfect_rational=: 2= +/@:%@,@factors</lang>
Exhaustive testing would take forever: <lang j> I.is_perfect_rational@"0 i.2^19 6 28 496 8128
I.is_perfect_rational@x:@"0 i.2^19x
6 28 496 8128</lang>
More limited testing takes reasonable amounts of time: <lang j> (#~ is_perfect_rational"0) (* <:@+:) 2^i.10x 6 28 496 8128</lang>
Mathematica
Mathematica has full support for fractions built-in. If one divides two exact numbers it will be left as a fraction if it can't be simplified. Comparison, addition, division, product et cetera are built-in: <lang Mathematica>4/16 3/8 8/4 4Pi/2 16!/10! Sqrt[9/16] Sqrt[3/4] (23/12)^5 2 + 1/(1 + 1/(3 + 1/4))
1/2+1/3+1/5 8/Pi+Pi/8 //Together 13/17 + 7/31 Sum[1/n,{n,1,100}] (*summation of 1/1 + 1/2 + 1/3 + 1/4+ .........+ 1/99 + 1/100*)
1/2-1/3 a=1/3;a+=1/7
1/4==2/8 1/4>3/8 Pi/E >23/20 1/3!=123/370 Sin[3]/Sin[2]>3/20
Numerator[6/9] Denominator[6/9]</lang> gives back: <lang Mathematica>1/4 3/8 2 2 Pi 5765760 3/4 Sqrt[3]/2 6436343 / 248832 47/17
31/30 (64+Pi^2) / (8 Pi) 522 / 527 14466636279520351160221518043104131447711 / 2788815009188499086581352357412492142272
1/6 10/21
True False True True True
2 3</lang> As you can see, Mathematica automatically handles fraction as exact things, it doesn't evaluate the fractions to a float. It only does this when either the numerator or the denominator is not exact. I only showed integers above, but Mathematica can handle symbolic fraction in the same and complete way: <lang Mathematica>c/(2 c) (b^2 - c^2)/(b - c) // Cancel 1/2 + b/c // Together</lang> gives back: <lang Mathematica>1/2 b+c (2 b+c) / (2 c)</lang> Moreover it does simplification like Sin[x]/Cos[x] => Tan[x]. Division, addition, subtraction, powering and multiplication of a list (of any dimension) is automatically threaded over the elements: <lang Mathematica>1+2*{1,2,3}^3</lang> gives back: <lang Mathematica>{3, 17, 55}</lang> To check for perfect numbers in the range 1 to 2^25 we can use: <lang Mathematica>found={}; CheckPerfect[num_Integer]:=If[Total[1/Divisors[num]]==2,AppendTo[found,num]]; Do[CheckPerfect[i],{i,1,2^25}]; found</lang> gives back: <lang Mathematica>{6, 28, 496, 8128, 33550336}</lang> Final note; approximations of fractions to any precision can be found using the function N.
Objective-C
File frac.h
<lang objc>#import <Foundation/Foundation.h>
@interface RCRationalNumber : NSObject {
@private int numerator; int denominator; BOOL autoSimplify; BOOL withSign;
} +(RCRationalNumber *)valueWithNumerator:(int)num andDenominator: (int)den; +(RCRationalNumber *)valueWithDouble: (double)fnum; +(RCRationalNumber *)valueWithInteger: (int)inum; +(RCRationalNumber *)valueWithRational: (RCRationalNumber *)rnum; -(id)init; -(id)initWithNumerator: (int)num andDenominator: (int)den; -(id)initWithDouble: (double)fnum precision: (int)prec; -(id)initWithInteger: (int)inum; -(id)initWithRational: (RCRationalNumber *)rnum; -(NSComparisonResult)compare: (RCRationalNumber *)rnum; -(id)simplify: (BOOL)act; -(void)setAutoSimplify: (BOOL)v; -(void)setWithSign: (BOOL)v; -(BOOL)autoSimplify; -(BOOL)withSign; -(NSString *)description; // ops -(id)multiply: (RCRationalNumber *)rnum; -(id)divide: (RCRationalNumber *)rnum; -(id)add: (RCRationalNumber *)rnum; -(id)sub: (RCRationalNumber *)rnum; -(id)abs; -(id)neg; -(id)mod: (RCRationalNumber *)rnum; -(int)sign; -(BOOL)isNegative; -(id)reciprocal; // getter -(int)numerator; -(int)denominator; //setter -(void)setNumerator: (int)num; -(void)setDenominator: (int)num; // defraction -(double)number; -(int)integer; @end</lang>
File frac.m
<lang objc>#import <Foundation/Foundation.h>
- import <math.h>
- import "frac.h"
// gcd: Greatest common divisor#Recursive_Euclid_algorithm // if built in as "private" function, add static.
static int lcm(int a, int b) {
return a / gcd(a,b) * b;
}
@implementation RCRationalNumber // initializers -(id)init {
NSLog(@"initialized to unity"); return [self initWithInteger: 1];
}
-(id)initWithNumerator: (int)num andDenominator: (int)den {
if ((self = [super init]) != nil) {
if (den == 0) {
NSLog(@"denominator is zero");
return nil;
}
[self setNumerator: num];
[self setDenominator: den];
[self setWithSign: YES];
[self setAutoSimplify: YES];
[self simplify: YES];
}
return self;
}
-(id)initWithInteger:(int)inum {
return [self initWithNumerator: inum andDenominator: 1];
}
-(id)initWithDouble: (double)fnum precision: (int)prec {
if ( prec > 9 ) prec = 9; double p = pow(10.0, (double)prec); int nd = (int)(fnum * p); return [self initWithNumerator: nd andDenominator: (int)p ];
}
-(id)initWithRational: (RCRationalNumber *)rnum {
return [self initWithNumerator: [rnum numerator] andDenominator: [rnum denominator]];
}
// comparing -(NSComparisonResult)compare: (RCRationalNumber *)rnum {
if ( [self number] > [rnum number] ) return NSOrderedDescending; if ( [self number] < [rnum number] ) return NSOrderedAscending; return NSOrderedSame;
}
// string rapresentation of the Q -(NSString *)description {
[self simplify: [self autoSimplify]]; return [NSString stringWithFormat: @"%@%d/%d", [self isNegative] ? @"-" :
( [self withSign] ? @"+" : @"" ), abs([self numerator]), [self denominator]]; }
// setter options -(void)setAutoSimplify: (BOOL)v {
autoSimplify = v; [self simplify: v];
} -(void)setWithSign: (BOOL)v {
withSign = v;
}
// getter for options -(BOOL)autoSimplify {
return autoSimplify;
}
-(BOOL)withSign {
return withSign;
}
// "simplify" the fraction ... -(id)simplify: (BOOL)act {
if ( act ) {
int common = gcd([self numerator], [self denominator]);
[self setNumerator: [self numerator]/common];
[self setDenominator: [self denominator]/common];
}
return self;
}
// diadic operators -(id)multiply: (RCRationalNumber *)rnum {
int newnum = [self numerator] * [rnum numerator]; int newden = [self denominator] * [rnum denominator]; return [RCRationalNumber valueWithNumerator: newnum
andDenominator: newden]; }
-(id)divide: (RCRationalNumber *)rnum {
return [self multiply: [rnum reciprocal]];
}
-(id)add: (RCRationalNumber *)rnum {
int common = lcm([self denominator], [rnum denominator]); int resnum = common / [self denominator] * [self numerator] + common / [rnum denominator] * [rnum numerator]; return [RCRationalNumber valueWithNumerator: resnum andDenominator: common];
}
-(id)sub: (RCRationalNumber *)rnum {
return [self add: [rnum neg]];
}
-(id)mod: (RCRationalNumber *)rnum {
return [[self divide: rnum]
sub: [RCRationalNumber valueWithInteger: [[self divide: rnum] integer]]]; }
// unary operators -(id)neg {
return [RCRationalNumber valueWithNumerator: -1*[self numerator]
andDenominator: [self denominator]]; }
-(id)abs {
return [RCRationalNumber valueWithNumerator: abs([self numerator])
andDenominator: [self denominator]]; }
-(id)reciprocal {
return [RCRationalNumber valueWithNumerator: [self denominator]
andDenominator: [self numerator]]; }
// get the sign -(int)sign {
return ([self numerator] < 0) ? -1 : 1;
}
// or just test if negativ -(BOOL)isNegative {
return ([self numerator] < 0) ? YES : NO;
}
// Q as real floating point -(double)number {
return (double)[self numerator] / (double)[self denominator];
}
// Q as (truncated) integer -(int)integer {
return [self numerator] / [self denominator];
}
// set num and den indipendently, fixing sign accordingly -(void)setNumerator: (int)num {
numerator = num;
}
-(void)setDenominator: (int)num {
if ( num < 0 ) numerator = -numerator; denominator = abs(num);
}
// getter -(int)numerator {
return numerator;
}
-(int)denominator {
return denominator;
}
// class method +(RCRationalNumber *)valueWithNumerator:(int)num andDenominator: (int)den {
return [[[RCRationalNumber alloc] initWithNumerator: num andDenominator: den] autorelease];
}
+(RCRationalNumber *)valueWithDouble: (double)fnum {
return [[[RCRationalNumber alloc] initWithDouble: fnum] autorelease];
}
+(RCRationalNumber *)valueWithInteger: (int)inum {
return [[[RCRationalNumber alloc] initWithInteger: inum] autorelease];
}
+(RCRationalNumber *)valueWithRational: (RCRationalNumber *)rnum {
return [[[RCRationalNumber alloc] initWithRational: rnum] autorelease];
} @end</lang>
Testing it.
<lang objc>#import <Foundation/Foundation.h>
- import "frac.h"
- import <math.h>
int main() {
NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];
int i;
for(i=2; i < 0x80000; i++) {
int candidate = i;
RCRationalNumber *sum = [RCRationalNumber valueWithNumerator: 1
andDenominator: candidate];
int factor;
for(factor=2; factor < sqrt((double)candidate); factor++) {
if ( (candidate % factor) == 0 ) {
sum = [[sum add: [RCRationalNumber valueWithNumerator: 1
andDenominator: factor]] add: [RCRationalNumber valueWithNumerator: 1 andDenominator: (candidate/factor)]];
}
}
if ( [sum denominator] == 1 ) {
printf("Sum of recipr. factors of %d = %d exactly %s\n",
candidate, [sum integer], ([sum integer]==1) ? "perfect!" : "");
} }
[pool release]; return 0;
}</lang>
OCaml
OCaml's Num library implements arbitrary-precision rational numbers:
<lang ocaml>#load "nums.cma";; open Num;;
for candidate = 2 to 1 lsl 19 do
let sum = ref (num_of_int 1 // num_of_int candidate) in
for factor = 2 to truncate (sqrt (float candidate)) do
if candidate mod factor = 0 then
sum := !sum +/ num_of_int 1 // num_of_int factor
+/ num_of_int 1 // num_of_int (candidate / factor)
done;
if is_integer_num !sum then
Printf.printf "Sum of recipr. factors of %d = %d exactly %s\n%!"
candidate (int_of_num !sum) (if int_of_num !sum = 1 then "perfect!" else "")
done;;</lang>
It might be implemented like this:
[insert implementation here]
Perl
Perl's Math::BigRat core module implements arbitrary-precision rational numbers. The bigrat pragma can be used to turn on transparent BigRat support:
<lang perl>use bigrat;
foreach my $candidate (2 .. 2**19) {
my $sum = 1 / $candidate;
foreach my $factor (2 .. sqrt($candidate)+1) {
if ($candidate % $factor == 0) {
$sum += 1 / $factor + 1 / ($candidate / $factor);
}
}
if ($sum->denominator() == 1) {
print "Sum of recipr. factors of $candidate = $sum exactly ", ($sum == 1 ? "perfect!" : ""), "\n";
}
}</lang>
It might be implemented like this:
[insert implementation here]
Python
Python 3's standard library already implements a Fraction class:
<lang python>from fractions import Fraction
for candidate in range(2, 2**19):
sum = Fraction(1, candidate)
for factor in range(2, int(candidate**0.5)+1):
if candidate % factor == 0:
sum += Fraction(1, factor) + Fraction(1, candidate // factor)
if sum.denominator == 1:
print("Sum of recipr. factors of %d = %d exactly %s" %
(candidate, int(sum), "perfect!" if sum == 1 else ""))</lang>
It might be implemented like this:
<lang python>def lcm(a, b):
return a // gcd(a,b) * b
def gcd(u, v):
return gcd(v, u%v) if v else abs(u)
class Fraction:
def __init__(self, numerator, denominator):
common = gcd(numerator, denominator)
self.numerator = numerator//common
self.denominator = denominator//common
def __add__(self, frac):
common = lcm(self.denominator, frac.denominator)
n = common // self.denominator * self.numerator + common // frac.denominator * frac.numerator
return Fraction(n, common)
def __sub__(self, frac):
return self.__add__(-frac)
def __neg__(self):
return Fraction(-self.numerator, self.denominator)
def __abs__(self):
return Fraction(abs(self.numerator), abs(self.denominator))
def __mul__(self, frac):
return Fraction(self.numerator * frac.numerator, self.denominator * frac.denominator)
def __div__(self, frac):
return self.__mul__(frac.reciprocal())
def reciprocal(self):
return Fraction(self.denominator, self.numerator)
def __cmp__(self, n):
return int(float(self) - float(n))
def __float__(self):
return float(self.numerator / self.denominator)
def __int__(self):
return (self.numerator // self.denominator)</lang>
Ruby
Ruby's standard library already implements a Rational class:
<lang ruby>require 'rational'
for candidate in 2 .. 2**19:
sum = Rational(1, candidate)
for factor in 2 ... candidate**0.5
if candidate % factor == 0
sum += Rational(1, factor) + Rational(1, candidate / factor)
end
end
if sum.denominator == 1
puts "Sum of recipr. factors of %d = %d exactly %s" %
[candidate, sum.to_i, sum == 1 ? "perfect!" : ""]
end
end</lang>
It might be implemented like this:
[insert implementation here]
Slate
Slate uses infinite-precision fractions transparently.
<lang slate>54 / 7. 20 reciprocal. (5 / 6) reciprocal. (5 / 6) as: Float.</lang>
Smalltalk
Smalltalk uses naturally and transparently fractions (through the class Fraction):
st> 54/7 54/7 st> 54/7 + 1 61/7 st> 54/7 < 50 true st> 20 reciprocal 1/20 st> (5/6) reciprocal 6/5 st> (5/6) asFloat 0.8333333333333334
<lang smalltalk>| sum | 2 to: (2 raisedTo: 19) do: [ :candidate |
sum := candidate reciprocal.
2 to: (candidate sqrt) do: [ :factor |
( (candidate \\ factor) = 0 )
ifTrue: [
sum := sum + (factor reciprocal) + ((candidate / factor) reciprocal)
]
].
( (sum denominator) = 1 )
ifTrue: [
('Sum of recipr. factors of %1 = %2 exactly %3' %
{ candidate printString .
(sum asInteger) printString .
( sum = 1 ) ifTrue: [ 'perfect!' ]
ifFalse: [ ' ' ] }) displayNl
]
].</lang>
Tcl
code to find factors of a number not shown <lang tcl>namespace eval rat {}
proc rat::new {args} {
if {[llength $args] == 0} {
set args {0}
}
lassign [split {*}$args] n d
if {$d == 0} {
error "divide by zero"
}
if {$d < 0} {
set n [expr {-1 * $n}]
set d [expr {abs($d)}]
}
return [normalize $n $d]
}
proc rat::split {args} {
if {[llength $args] == 1} {
lassign [::split $args /] n d
if {$d eq ""} {
set d 1
}
} else {
lassign $args n d
}
return [list $n $d]
}
proc rat::join {rat} {
lassign $rat n d
if {$n == 0} {
return 0
} elseif {$d == 1} {
return $n
} else {
return $n/$d
}
}
proc rat::normalize {n d} {
set gcd [gcd $n $d]
return [join [list [expr {$n/$gcd}] [expr {$d/$gcd}]]]
}
proc rat::gcd {a b} {
while {$b != 0} {
lassign [list $b [expr {$a % $b}]] a b
}
return $a
}
proc rat::abs {rat} {
lassign [split $rat] n d
return [join [list [expr {abs($n)}] $d]]
}
proc rat::inv {rat} {
lassign [split $rat] n d return [normalize $d $n]
}
proc rat::+ {args} {
set n 0
set d 1
foreach arg $args {
lassign [split $arg] an ad
set n [expr {$n*$ad + $an*$d}]
set d [expr {$d * $ad}]
}
return [normalize $n $d]
}
proc rat::- {args} {
lassign [split [lindex $args 0]] n d
if {[llength $args] == 1} {
return [join [list [expr {-1 * $n}] $d]]
}
foreach arg [lrange $args 1 end] {
lassign [split $arg] an ad
set n [expr {$n*$ad - $an*$d}]
set d [expr {$d * $ad}]
}
return [normalize $n $d]
}
proc rat::* {args} {
set n 1
set d 1
foreach arg $args {
lassign [split $arg] an ad
set n [expr {$n * $an}]
set d [expr {$d * $ad}]
}
return [normalize $n $d]
}
proc rat::/ {a b} {
set r [* $a [inv $b]]
if {[string match */0 $r]} {
error "divide by zero"
}
return $r
}
proc rat::== {a b} {
return [expr {[- $a $b] == 0}]
}
proc rat::!= {a b} {
return [expr { ! [== $a $b]}]
}
proc rat::< {a b} {
lassign [split [- $a $b]] n d
return [expr {$n < 0}]
}
proc rat::> {a b} {
lassign [split [- $a $b]] n d
return [expr {$n > 0}]
}
proc rat::<= {a b} {
return [expr { ! [> $a $b]}]
}
proc rat::>= {a b} {
return [expr { ! [< $a $b]}]
}
proc is_perfect {num} {
set sum [rat::new 0]
foreach factor [all_factors $num] {
set sum [rat::+ $sum [rat::new 1/$factor]]
}
# note, all_factors includes 1, so sum should be 2
return [rat::== $sum 2]
}
proc get_perfect_numbers {} {
set t [clock seconds]
set limit [expr 2**19]
for {set num 2} {$num < $limit} {incr num} {
if {[is_perfect $num]} {
puts "perfect: $num"
}
}
puts "elapsed: [expr {[clock seconds] - $t}] seconds"
set num [expr {2**12 * (2**13 - 1)}] ;# 5th perfect number
if {[is_perfect $num]} {
puts "perfect: $num"
}
}
source primes.tcl get_perfect_numbers</lang>
perfect: 6 perfect: 28 perfect: 496 perfect: 8128 elapsed: 477 seconds perfect: 33550336
TI-89 BASIC
While TI-89 BASIC has built-in rational and symbolic arithmetic, it does not have user-defined data types.