Source code in the icon programming language

procedure main()
   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


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


   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