Arithmetic/Rational
From Rosetta Code
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
Contents |
[edit] Ada
[edit] ALGOL 68
Works with: ALGOL 68 version Standard - no extensions to language used
Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386
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:
OP +:= = (REF FRAC a, FRAC b)REF FRAC: ( a := a + b ),
+=: = (FRAC a, REF FRAC b)REF FRAC: ( b := a + b ),
-:= = (REF FRAC a, FRAC b)REF FRAC: ( a := a - b ),
*:= = (REF FRAC a, FRAC b)REF FRAC: ( a := a * b ),
/:= = (REF FRAC a, FRAC b)REF FRAC: ( a := a / b ),
%:= = (REF FRAC a, FRAC b)REF FRAC: ( a := FRACINIT (a % b) ),
%*:= = (REF FRAC a, FRAC b)REF FRAC: ( a := a %* b );
# OP aliases for extended character sets (eg: Unicode, APL, ALCOR and GOST 10859) #
OP × = (FRAC a, b)FRAC: a * b,
÷ = (FRAC a, b)INT: a OVER b,
÷× = (FRAC a, b)FRAC: a MOD b,
÷* = (FRAC a, b)FRAC: a MOD b,
%× = (FRAC a, b)FRAC: a MOD b,
≤ = (FRAC a, b)FRAC: a <= b,
≥ = (FRAC a, b)FRAC: a >= b,
≠ = (FRAC a, b)BOOL: a /= b,
↑ = (FRAC frac, INT exponent)FRAC: frac ** exponent,
÷×:= = (REF FRAC a, FRAC b)REF FRAC: ( a := a MOD b ),
%×:= = (REF FRAC a, FRAC b)REF FRAC: ( a := a MOD b ),
÷*:= = (REF FRAC a, FRAC b)REF FRAC: ( a := a MOD b );
# BOLD aliases for CPU that only support uppercase for 6-bit bytes - wrist watches #
OP OVER = (FRAC a, b)INT: a % b,
MOD = (FRAC a, b)FRAC: a %*b,
LT = (FRAC a, b)BOOL: a < b,
GT = (FRAC a, b)BOOL: a > b,
LE = (FRAC a, b)BOOL: a <= b,
GE = (FRAC a, b)BOOL: a >= b,
EQ = (FRAC a, b)BOOL: a = b,
NE = (FRAC a, b)BOOL: a /= b,
UP = (FRAC frac, INT exponent)FRAC: frac**exponent;
# the required standard assignment operators #
OP PLUSAB = (REF FRAC a, FRAC b)REF FRAC: ( a +:= b ), # PLUS #
PLUSTO = (FRAC a, REF FRAC b)REF FRAC: ( a +=: b ), # PRUS #
MINUSAB = (REF FRAC a, FRAC b)REF FRAC: ( a *:= b ),
DIVAB = (REF FRAC a, FRAC b)REF FRAC: ( a /:= b ),
OVERAB = (REF FRAC a, FRAC b)REF FRAC: ( a %:= b ),
MODAB = (REF FRAC a, FRAC b)REF FRAC: ( a %*:= b );
END 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
[edit] C
[edit] C++
Library: Boost
Boost provides a rational number template.
#include <iostream>
#include "math.h"
#include "boost/rational.hpp"
typedef boost::rational<int> frac;
bool is_perfect(int c)
{
frac sum(1, c);
for (int f = 2;f < sqrt(static_cast<float>(c)); ++f){
if (c % f == 0) sum += frac(1,f) + frac(1, c/f);
}
if (sum.denominator() == 1){
return (sum == 1);
}
return false;
}
int main()
{
for (int candidate = 2; candidate < 0x80000; ++candidate){
if (is_perfect(candidate))
std::cout << candidate << " is perfect" << std::endl;
}
return 0;
}
[edit] 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.
(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)
[edit] Forth
\ Rationals can use any double cell operations: 2!, 2@, 2dup, 2swap, etc.
\ Uses the stack convention of the built-in "*/" for int * frac -> int
: numerator drop ;
: denominator nip ;
: s>rat 1 ; \ integer to rational (n/1)
: rat>s / ; \ integer
: rat>frac mod ; \ fractional part
: rat>float swap s>f s>f f/ ;
: rat. swap 1 .r [char] / emit . ;
\ normalize: factors out gcd and puts sign into numerator
: gcd ( a b -- gcd ) begin ?dup while tuck mod repeat ;
: rat-normalize ( rat -- rat ) 2dup gcd tuck / >r / r> ;
: rat-abs swap abs swap ;
: rat-negate swap negate swap ;
: 1/rat over 0< if negate swap negate else swap then ;
: rat+ ( a b c d -- ad+bc bd )
rot 2dup * >r
rot * >r * r> +
r> rat-normalize ;
: rat- rat-negate rat+ ;
: rat* ( a b c d -- ac bd )
rot * >r * r> rat-normalize ;
: rat/ swap rat* ;
: rat-equal d= ;
: rat-less ( a b c d -- ad<bc )
-rot * >r * r> < ;
: rat-more 2swap rat-less ;
: rat-inc tuck + swap ;
: rat-dec tuck - swap ;
[edit] Fortran
Works with: Fortran version 90 and later
module module_rational
implicit none
private
public :: rational
public :: rational_simplify
public :: assignment (=)
public :: operator (//)
public :: operator (+)
public :: operator (-)
public :: operator (*)
public :: operator (/)
public :: operator (<)
public :: operator (<=)
public :: operator (>)
public :: operator (>=)
public :: operator (==)
public :: operator (/=)
public :: abs
public :: int
public :: modulo
type rational
integer :: numerator
integer :: denominator
end type rational
interface assignment (=)
module procedure assign_rational_int, assign_rational_real
end interface
interface operator (//)
module procedure make_rational
end interface
interface operator (+)
module procedure rational_add
end interface
interface operator (-)
module procedure rational_minus, rational_subtract
end interface
interface operator (*)
module procedure rational_multiply
end interface
interface operator (/)
module procedure rational_divide
end interface
interface operator (<)
module procedure rational_lt
end interface
interface operator (<=)
module procedure rational_le
end interface
interface operator (>)
module procedure rational_gt
end interface
interface operator (>=)
module procedure rational_ge
end interface
interface operator (==)
module procedure rational_eq
end interface
interface operator (/=)
module procedure rational_ne
end interface
interface abs
module procedure rational_abs
end interface
interface int
module procedure rational_int
end interface
interface modulo
module procedure rational_modulo
end interface
contains
recursive function gcd (i, j) result (res)
integer, intent (in) :: i
integer, intent (in) :: j
integer :: res
if (j == 0) then
res = i
else
res = gcd (j, modulo (i, j))
end if
end function gcd
function rational_simplify (r) result (res)
type (rational), intent (in) :: r
type (rational) :: res
integer :: g
g = gcd (r % numerator, r % denominator)
res = r % numerator / g // r % denominator / g
end function rational_simplify
function make_rational (numerator, denominator) result (res)
integer, intent (in) :: numerator
integer, intent (in) :: denominator
type (rational) :: res
res = rational (numerator, denominator)
end function make_rational
subroutine assign_rational_int (res, i)
type (rational), intent (out), volatile :: res
integer, intent (in) :: i
res = i // 1
end subroutine assign_rational_int
subroutine assign_rational_real (res, x)
type (rational), intent(out), volatile :: res
real, intent (in) :: x
integer :: x_floor
real :: x_frac
x_floor = floor (x)
x_frac = x - x_floor
if (x_frac == 0) then
res = x_floor // 1
else
res = (x_floor // 1) + (1 // floor (1 / x_frac))
end if
end subroutine assign_rational_real
function rational_add (r, s) result (res)
type (rational), intent (in) :: r
type (rational), intent (in) :: s
type (rational) :: res
res = r % numerator * s % denominator + r % denominator * s % numerator // &
& r % denominator * s % denominator
end function rational_add
function rational_minus (r) result (res)
type (rational), intent (in) :: r
type (rational) :: res
res = - r % numerator // r % denominator
end function rational_minus
function rational_subtract (r, s) result (res)
type (rational), intent (in) :: r
type (rational), intent (in) :: s
type (rational) :: res
res = r % numerator * s % denominator - r % denominator * s % numerator // &
& r % denominator * s % denominator
end function rational_subtract
function rational_multiply (r, s) result (res)
type (rational), intent (in) :: r
type (rational), intent (in) :: s
type (rational) :: res
res = r % numerator * s % numerator // r % denominator * s % denominator
end function rational_multiply
function rational_divide (r, s) result (res)
type (rational), intent (in) :: r
type (rational), intent (in) :: s
type (rational) :: res
res = r % numerator * s % denominator // r % denominator * s % numerator
end function rational_divide
function rational_lt (r, s) result (res)
type (rational), intent (in) :: r
type (rational), intent (in) :: s
type (rational) :: r_simple
type (rational) :: s_simple
logical :: res
r_simple = rational_simplify (r)
s_simple = rational_simplify (s)
res = r_simple % numerator * s_simple % denominator < &
& s_simple % numerator * r_simple % denominator
end function rational_lt
function rational_le (r, s) result (res)
type (rational), intent (in) :: r
type (rational), intent (in) :: s
type (rational) :: r_simple
type (rational) :: s_simple
logical :: res
r_simple = rational_simplify (r)
s_simple = rational_simplify (s)
res = r_simple % numerator * s_simple % denominator <= &
& s_simple % numerator * r_simple % denominator
end function rational_le
function rational_gt (r, s) result (res)
type (rational), intent (in) :: r
type (rational), intent (in) :: s
type (rational) :: r_simple
type (rational) :: s_simple
logical :: res
r_simple = rational_simplify (r)
s_simple = rational_simplify (s)
res = r_simple % numerator * s_simple % denominator > &
& s_simple % numerator * r_simple % denominator
end function rational_gt
function rational_ge (r, s) result (res)
type (rational), intent (in) :: r
type (rational), intent (in) :: s
type (rational) :: r_simple
type (rational) :: s_simple
logical :: res
r_simple = rational_simplify (r)
s_simple = rational_simplify (s)
res = r_simple % numerator * s_simple % denominator >= &
& s_simple % numerator * r_simple % denominator
end function rational_ge
function rational_eq (r, s) result (res)
type (rational), intent (in) :: r
type (rational), intent (in) :: s
logical :: res
res = r % numerator * s % denominator == s % numerator * r % denominator
end function rational_eq
function rational_ne (r, s) result (res)
type (rational), intent (in) :: r
type (rational), intent (in) :: s
logical :: res
res = r % numerator * s % denominator /= s % numerator * r % denominator
end function rational_ne
function rational_abs (r) result (res)
type (rational), intent (in) :: r
type (rational) :: res
res = sign (r % numerator, r % denominator) // r % denominator
end function rational_abs
function rational_int (r) result (res)
type (rational), intent (in) :: r
integer :: res
res = r % numerator / r % denominator
end function rational_int
function rational_modulo (r) result (res)
type (rational), intent (in) :: r
integer :: res
res = modulo (r % numerator, r % denominator)
end function rational_modulo
end module module_rational
Example:
program perfect_numbers
use module_rational
implicit none
integer, parameter :: n_min = 2
integer, parameter :: n_max = 2 ** 19 - 1
integer :: n
integer :: factor
type (rational) :: sum
do n = n_min, n_max
sum = 1 // n
factor = 2
do
if (factor * factor >= n) then
exit
end if
if (modulo (n, factor) == 0) then
sum = rational_simplify (sum + (1 // factor) + (factor // n))
end if
factor = factor + 1
end do
if (sum % numerator == 1 .and. sum % denominator == 1) then
write (*, '(i0)') n
end if
end do
end program perfect_numbers
Output:
6
28
496
8128
[edit] Go
Go's big package implements arbitrary-precision integers and rational numbers.
package main
import (
"fmt"
"big"
"math"
)
func main() {
max := int64(1<<19)
for candidate := int64(2); candidate < max; candidate++ {
sum := big.NewRat(1, candidate)
max2 := int64(math.Sqrt(float64(candidate)))
for factor := int64(2); factor <= max2; factor++ {
if candidate % factor == 0 {
sum.Add(sum, new(big.Rat).Add(big.NewRat(1, factor), big.NewRat(1, candidate / factor)))
}
}
if sum.Denom().Int64() == 1 {
perfectstring := ""
if sum.Num().Int64() == 1 { perfectstring = "perfect!" }
fmt.Printf("Sum of recipr. factors of %d = %d exactly %s\n",
candidate, sum.Num().Int64(), perfectstring)
}
}
}
It might be implemented like this:
[insert implementation here]
[edit] 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.
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]
For a sample implementation of Ratio, see the Haskell 98 Report.
[edit] Icon and Unicon
[edit] Icon
The IPL provides support for rational arithmetic
- The data type is called 'rational' not 'frac'.
- Use the record constructor 'rational' to create a rational. Sign must be 1 or -1.
- Neither Icon nor Unicon supports operator overloading. Augmented assignments make little sense w/o this.
- Procedures include 'negrat' (unary -), 'addrat' (+), 'subrat' (-), 'mpyrat' (*), 'divrat' (modulo /).
Additional procedures are implemented here to complete the task:
- 'makerat' (make), 'absrat' (abs), 'eqrat' (=), 'nerat' (~=), 'ltrat' (<), 'lerat' (<=), 'gerat' (>=), 'gtrat' (>)
procedure main()Sample output:
limit := 2^19
write("Perfect numbers up to ",limit," (using rational arithmetic):")
every write(is_perfect(c := 2 to limit))
write("End of perfect numbers")
# verify the rest of the implementation
zero := makerat(0) # from integer
half := makerat(0.5) # from real
qtr := makerat("1/4") # from strings ...
one := makerat("1")
mone := makerat("-1")
verifyrat("eqrat",zero,zero)
verifyrat("ltrat",zero,half)
verifyrat("ltrat",half,zero)
verifyrat("gtrat",zero,half)
verifyrat("gtrat",half,zero)
verifyrat("nerat",zero,half)
verifyrat("nerat",zero,zero)
verifyrat("absrat",mone,)
end
procedure is_perfect(c) #: test for perfect numbers using rational arithmetic
rsum := rational(1, c, 1)
every f := 2 to sqrt(c) do
if 0 = c % f then
rsum := addrat(rsum,addrat(rational(1,f,1),rational(1,integer(c/f),1)))
if rsum.numer = rsum.denom = 1 then
return c
end
Perfect numbers up to 524288 (using rational arithmetic): 6 28 496 8128 End of perfect numbers Testing eqrat( (0/1), (0/1) ) ==> returned (0/1) Testing ltrat( (0/1), (1/2) ) ==> returned (1/2) Testing ltrat( (1/2), (0/1) ) ==> failed Testing gtrat( (0/1), (1/2) ) ==> failed Testing gtrat( (1/2), (0/1) ) ==> returned (0/1) Testing nerat( (0/1), (1/2) ) ==> returned (1/2) Testing nerat( (0/1), (0/1) ) ==> failed Testing absrat( (-1/1), ) ==> returned (1/1)
The following task functions are missing from the IPL:
procedure verifyrat(p,r1,r2) #: verification tests for rational procedures
return write("Testing ",p,"( ",rat2str(r1),", ",rat2str(\r2) | &null," ) ==> ","returned " || rat2str(p(r1,r2)) | "failed")
end
procedure makerat(x) #: make rational (from integer, real, or strings)
local n,d
static c
initial c := &digits++'+-'
return case type(x) of {
"real" : real2rat(x)
"integer" : ratred(rational(x,1,1))
"string" : if x ? ( n := integer(tab(many(c))), ="/", d := integer(tab(many(c))), pos(0)) then
ratred(rational(n,d,1))
else
makerat(numeric(x))
}
end
procedure absrat(r1) #: abs(rational)
r1 := ratred(r1)
r1.sign := 1
return r1
end
invocable all # for string invocation
procedure xoprat(op,r1,r2) #: support procedure for binary operations that cross denominators
local numer, denom, div
r1 := ratred(r1)
r2 := ratred(r2)
return if op(r1.numer * r2.denom,r2.numer * r1.denom) then r2 # return right argument on success
end
procedure eqrat(r1,r2) #: rational r1 = r2
return xoprat("=",r1,r2)
end
procedure nerat(r1,r2) #: rational r1 ~= r2
return xoprat("~=",r1,r2)
end
procedure ltrat(r1,r2) #: rational r1 < r2
return xoprat("<",r1,r2)
end
procedure lerat(r1,r2) #: rational r1 <= r2
return xoprat("<=",r1,r2)
end
procedure gerat(r1,r2) #: rational r1 >= r2
return xoprat(">=",r1,r2)
end
procedure gtrat(r1,r2) #: rational r1 > r2
return xoprat(">",r1,r2)
end
link rational
The Library: Icon Programming Library provides rational and gcd in numbers. Record definition and usage is shown below:
record rational(numer, denom, sign) # rational type
addrat(r1,r2) # Add rational numbers r1 and r2.
divrat(r1,r2) # Divide rational numbers r1 and r2.
medrat(r1,r2) # Form mediant of r1 and r2.
mpyrat(r1,r2) # Multiply rational numbers r1 and r2.
negrat(r) # Produce negative of rational number r.
rat2real(r) # Produce floating-point approximation of r
rat2str(r) # Convert the rational number r to its string representation.
real2rat(v,p) # Convert real to rational with precision p (default 1e-10). Warning: excessive p gives ugly fractions
reciprat(r) # Produce the reciprocal of rational number r.
str2rat(s) # Convert the string representation (such as "3/2") to a rational number
subrat(r1,r2) # Subtract rational numbers r1 and r2.
gcd(i, j) # returns greatest common divisor of i and j
[edit] Unicon
The Icon solution works in Unicon.
[edit] J
J implements rational numbers:
3r4*2r5
3r10
That said, note that J's floating point numbers work just fine for the stated problem:
is_perfect_rational=: 2 = (1 + i.) +/@:%@([ #~ 0 = |) ]
faster version (but the problem, as stated, is still tremendously inefficient):
factors=: */&>@{@((^ i.@>:)&.>/)@q:~&__
is_perfect_rational=: 2= +/@:%@,@factors
Exhaustive testing would take forever:
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
More limited testing takes reasonable amounts of time:
(#~ is_perfect_rational"0) (* <:@+:) 2^i.10x
6 28 496 8128
[edit] JavaScript
See Rational Arithmetic/JavaScript
[edit] Lua
function gcd(a,b) return a == 0 and b or gcd(b % a, a) end
do
local function coerce(a, b)
if type(a) == "number" then return rational(a, 1), b end
if type(b) == "number" then return a, rational(b, 1) end
return a, b
end
rational = setmetatable({
__add = function(a, b)
local a, b = coerce(a, b)
return rational(a.num * b.den + a.den * b.num, a.den * b.den)
end,
__sub = function(a, b)
local a, b = coerce(a, b)
return rational(a.num * b.den - a.den * b.num, a.den * b.den)
end,
__mul = function(a, b)
local a, b = coerce(a, b)
return rational(a.num * b.num, a.den * b.den)
end,
__div = function(a, b)
local a, b = coerce(a, b)
return rational(a.num * b.den, a.den * b.num)
end,
__pow = function(a, b)
if type(a) == "number" then return a ^ (b.num / b.den) end
return rational(a.num ^ b, a.den ^ b) --runs into a problem if these aren't integers
end,
__concat = function(a, b)
if getmetatable(a) == rational then return a.num .. "/" .. a.den .. b end
return a .. b.num .. "/" .. b.den
end,
__unm = function(a) return rational(-a.num, -a.den) end}, {
__call = function(z, a, b) return setmetatable({num = a / gcd(a, b),den = b / gcd(a, b)}, z) end} )
end
print(rational(2, 3) + rational(3, 5) - rational(1, 10) .. "") --> 7/6
print((rational(4, 5) * rational(5, 9)) ^ rational(1, 2) .. "") --> 2/3
print(rational(45, 60) / rational(5, 2) .. "") --> 3/10
print(5 + rational(1, 3) .. "") --> 16/3
function findperfs(n)
local ret = {}
for i = 1, n do
sum = rational(1, i)
for fac = 2, i^.5 do
if i % fac == 0 then
sum = sum + rational(1, fac) + rational(fac, i)
end
end
if sum.den == sum.num then
ret[#ret + 1] = i
end
end
return table.concat(ret, '\n')
end
print(findperfs(2^19))
[edit] 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:
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]
gives back:
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
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:
c/(2 c)
(b^2 - c^2)/(b - c) // Cancel
1/2 + b/c // Together
gives back:
1/2
b+c
(2 b+c) / (2 c)
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:
1+2*{1,2,3}^3
gives back:
{3, 17, 55}
To check for perfect numbers in the range 1 to 2^25 we can use:
found={};
CheckPerfect[num_Integer]:=If[Total[1/Divisors[num]]==2,AppendTo[found,num]];
Do[CheckPerfect[i],{i,1,2^25}];
found
gives back:
{6, 28, 496, 8128, 33550336}
Final note; approximations of fractions to any precision can be found using the function N.
[edit] Objective-C
See Rational Arithmetic/Objective-C
[edit] OCaml
OCaml's Num library implements arbitrary-precision rational numbers:
#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;;
It might be implemented like this:
[insert implementation here]
[edit] Perl
Perl's Math::BigRat core module implements arbitrary-precision rational numbers. The bigrat pragma can be used to turn on transparent BigRat support:
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";
}
}
It might be implemented like this:
[insert implementation here]
[edit] Perl 6
Perl 6 supports rational arithmetic natively.
for 2..2**19 -> $candidate {
my $sum = 1 / $candidate;
for 2 .. ceiling(sqrt($candidate)) -> $factor {
if $candidate %% $factor {
$sum += 1 / $factor + 1 / ($candidate / $factor);
}
}
if $sum.denominator == 1 {
say "Sum of reciprocal factors of $candidate = $sum exactly", ($sum == 1 ?? ", perfect!" !! ".");
}
}
Note also that ordinary decimal literals are stored as Rats, so the following loop always stops exactly on 10 despite 0.1 not being exactly representable in floating point:
for 1.0, 1.1, 1.2 ... 10 { .say }
The arithmetic is all done in rationals, which are converted to floating-point just before display so that people don't have to puzzle out what 53/10 means.
[edit] PicoLisp
(load "@lib/frac.l")
(for (N 2 (> (** 2 19) N) (inc N))
(let (Sum (frac 1 N) Lim (sqrt N))
(for (F 2 (>= Lim F) (inc F))
(when (=0 (% N F))
(setq Sum
(f+ Sum
(f+ (frac 1 F) (frac 1 (/ N F))) ) ) ) )
(when (= 1 (cdr Sum))
(prinl
"Perfect " N
", sum is " (car Sum)
(and (= 1 (car Sum)) ": perfect") ) ) ) )
Output:
Perfect 6, sum is 1: perfect Perfect 28, sum is 1: perfect Perfect 120, sum is 2 Perfect 496, sum is 1: perfect Perfect 672, sum is 2 Perfect 8128, sum is 1: perfect Perfect 30240, sum is 3 Perfect 32760, sum is 3 Perfect 523776, sum is 2
[edit] Python
Works with: Python version 3.0
Python 3's standard library already implements a Fraction class:
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 ""))
It might be implemented like this:
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)
[edit] Ruby
Ruby's standard library already implements a Rational class:
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
It might be implemented like this:
[insert implementation here]
[edit] Scala
class Rational(n: Long, d:Long) extends Ordered[Rational]
{
require(d!=0)
private val g:Long = gcd(n, d)
val numerator:Long = n/g
val denominator:Long = d/g
def this(n:Long)=this(n,1)
def +(that:Rational):Rational=new Rational(
numerator*that.denominator + that.numerator*denominator,
denominator*that.denominator)
def -(that:Rational):Rational=new Rational(
numerator*that.denominator - that.numerator*denominator,
denominator*that.denominator)
def *(that:Rational):Rational=
new Rational(numerator*that.numerator, denominator*that.denominator)
def /(that:Rational):Rational=
new Rational(numerator*that.denominator, that.numerator*denominator)
def unary_~ :Rational=new Rational(denominator, numerator)
def unary_- :Rational=new Rational(-numerator, denominator)
def abs :Rational=new Rational(Math.abs(numerator), Math.abs(denominator))
override def compare(that:Rational):Int=
(this.numerator*that.denominator-that.numerator*this.denominator).toInt
override def toString()=numerator+"/"+denominator
private def gcd(x:Long, y:Long):Long=
if(y==0) x else gcd(y, x%y)
}
object Rational
{
def apply(n: Long, d:Long)=new Rational(n,d)
def apply(n:Long)=new Rational(n)
implicit def longToRational(i:Long)=new Rational(i)
}
def find_perfects():Unit=
{
for (candidate <- 2 until 1<<19)
{
var sum= ~Rational(candidate)
for (factor <- 2 until (Math.sqrt(candidate)+1).toInt)
{
if (candidate%factor==0)
sum+= ~Rational(factor)+ ~Rational(candidate/factor)
}
if (sum.denominator==1 && sum.numerator==1)
printf("Perfect number %d sum is %s\n", candidate, sum)
}
}
[edit] Slate
Slate uses infinite-precision fractions transparently.
54 / 7.
20 reciprocal.
(5 / 6) reciprocal.
(5 / 6) as: Float.
[edit] 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
Works with: GNU 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
]
].
[edit] Tcl
[edit] TI-89 BASIC
While TI-89 BASIC has built-in rational and symbolic arithmetic, it does not have user-defined data types.
Categories: Programming Tasks | Arithmetic operations | Ada | ALGOL 68 | C | C++ | Boost | Common Lisp | Forth | Fortran | Go | Haskell | Icon | Icon Programming Library | Unicon | J | JavaScript | Lua | Mathematica | Objective-C | OCaml | Perl | Perl 6 | PicoLisp | Python | Ruby | Scala | Slate | Smalltalk | Tcl | TI-89 BASIC
