Here are implementations of handlers that returns true if a positive integer is a prime, breaks a positive integer into its prime factors, returns the greatest common divisor of two positive integers, and one that returns the least common multiple of two positive integers.
Enjoy.
(* Made by McUsr 2014-8-12 Source: Kenneth H. Rosen Discrete Mathematics and its Applications.
http://macscripter.net/viewtopic.php?pid=174425#p174425
Version 2 : 2014-8-14 More efficient versions of gcd and lcm, added isCoprime
Version 3 : 2014-8-14 Addend div_alg (Number theory divison)
Version 4: 2014-8-14 Added divides handler, and modulus handler. Bugfix div_alg.
Version 5: 2014-8-14 further bugfix div_alg, corrected usage of dividend and divisor.
Version 6: 2014-8-15 Added relprime handler , expansionToBase2, expansionFromBase2
Version 7: 2014-8-15 Added xgcd Extended euclid for solving linear congruences.
Version 8: 2015-5-25 Added modular_expansion, ceil, floor log2 and ln2.
Version 9: 2015-5-27 Added ln2. -Handlers starting with ln returns a real, handlers starting with log, ints.
Version 10: 2015-6-7 Added P(n,r), c(n,r) and Factorial handlers for simple combinatorics.
Version 10.1: 2015-6-7 Removed bug in c(n,r) handler.
*)
to gcd(a, b)
# Returns the greatest common divisor, this is the 'Euclid' algorithm.
set x to a
set y to b
repeat while y ≠0
set r to x mod y
set x to y
set y to r
end repeat
return x
end gcd
on xgcd(a, b)
(*
Extended euclid for solving linear congruences.
returns a triple: {g,s,t}, so that g = sa * tb
Source: Rosens Discrete Mathematics,
Students solutions Manual, p. 107
*)
set {x, y} to {a, b}
set {oldolds, olds} to {1, 0}
set {oldoldt, oldt} to {0, 1}
repeat while y ≠0
set q to x div y
set r to x mod y
set x to y
set y to r
set s to oldolds - q * olds
set t to oldoldt - q * oldt
set oldolds to olds
set oldoldt to oldt
set olds to s
set oldt to t
end repeat
# x is gcd(a,b) oldold *a + oldoldt *b = x
return {x, oldolds, oldoldt}
end xgcd
to lcm(a, b)
(*
Returns the least common denominator. We take advantage of the fact that
a*b=lcm*gcd . . . (proof elsewhere)
*)
return (a * b) / (gcd(a, b)) as integer
end lcm
to isCoprime(a, b)
# Returns true if the greatest common divisor of a and b is 1
if gcd(a, b) = 1 then
return true
else
return false
end if
end isCoprime
on relprime(L)
# Returns true if all the integers in the list are coprime (relative prime) with each other
set lcount to (count L)
repeat with i from 1 to (lcount - 1)
repeat with j from (i + 1) to lcount
if gcd(item i of L, item j of L) ≠1 then return false
end repeat
end repeat
return true
end relprime
to isPrime(anInt)
# Returns true or false if it is a prime . .
return isPrimeNumber(anInt, false, 2)
end isPrime
to isPrimeNumber(anInt, retDivisor, startval)
(*
Returns true if anInt is prime.
False, or the factor if it isn't.
It takes a startval for probing for not prime, 2 is the ultimate startval.
anInt: Should have an integer as argument, we do not test.
retDivisor: Sometimes we want the divisor back if there was one.
startval: the value to start testing from: optimization for factorization.
If n is a composite integer then n has a prime divisor less than or equal to
sqrt(n) Therefore it is a primenumber if we can't find a factorization
below or equal to sqrt(n)
*)
local ceiling
# if class of anInt is not integer then error "Is no integer" number 3128
if startval = 2 then
if anInt mod 2 = 0 then
if not retDivisor then return false
return 2
else
set startval to 3
end if
end if
set ceiling to (anInt ^ 0.5) div 1
repeat with i from startval to ceiling by 2
if anInt mod i = 0 then
if not retDivisor then return false
return i
end if
end repeat
return true
end isPrimeNumber
to primefactors(thenum)
# Returns a list with primefactors for an integer
local isDone, startval, theFactor, factors
set factors to {}
set isDone to false
set startval to 2
repeat while isDone = false
set theFactor to isPrimeNumber(thenum, true, startval)
if theFactor = true then
set isDone to true
set end of factors to thenum
else
set startval to theFactor
set thenum to thenum div theFactor
set end of factors to theFactor
end if
end repeat
return factors
end primefactors
on div_alg(dividend, divisor)
(*
dividend : integer
divisor : positive integer
returns: {q,r} 0<= r<dividend
This is the division algorithm used in number theory.
The remainder is supposed to always be positive so we can use the result to form linear combinations.
Example:
div_alg(-11,2)
returns: {-6,1} => -11 = -6 * 2 + 1
( as a linear combination expressing -11 )
*)
set r to dividend mod divisor
set q to dividend div divisor
if r < 0 then
set q to q - 1
set r to -((divisor * q) - dividend)
end if
return {q, r}
end div_alg
on modulus(dividend, divisor)
# returns always a positive modulus.
return item 2 of div_alg(dividend, divisor)
end modulus
on divides(dividend, divisor)
(*
Utility handler:
returns true if dividend divides divisor, (without a rest )
*)
return (dividend mod divisor) = 0
end divides
(* ====================================================== *)
(* This package for base changing uses binary as intermediary base. *)
(* ====================================================== *)
on expansionToBase2(n, base_b)
(*
Returns a string with the number converted from the base.
log expansionToBase2(241, 10)
--> (*11110001*)
log expansionToBase2("3EBC", 16)
--> (*0011111010111100*)
log expansionToBase2("17", 8)
--> (*001111*)
One usage is to represent a large power as a binary,
because we can then treat the large number in pieces.
For instance 3^241 = 3^128 * 3^64*3^32*3^16*3
*)
if base_b = 10 then
set q to n
set L to {}
repeat while q ≠0
set end of L to q mod 2
set q to q div 2
end repeat
return reverse of L as string
else if base_b = 16 then
set hexdigitstring to "0123456789ABCDEF"
set hs to reverse of (characters of n)
if hs = {} then return -1
set hstringCount to count hs
set L to {}
repeat with i from 1 to hstringCount
set q to (offset of (item i of hs) in hexdigitstring) - 1
if q = 0 then
repeat 4 times
set end of L to 0
end repeat
else
# we pad out missing zeros
set ctr to 0
repeat while q ≠0
set end of L to q mod 2
set ctr to ctr + 1
set q to q div 2
end repeat
repeat while ctr < 4
set end of L to 0
set ctr to ctr + 1
end repeat
end if
end repeat
return reverse of L as string
else if base_b = 8 then
set os to reverse of (characters of n)
if os = {} then return -1
set octCount to count os
set L to {}
repeat with i from 1 to octCount
set q to item i of os
if q = 0 then
repeat 3 times
set end of L to 0
end repeat
else
# we pad out missing zeros
set ctr to 0
repeat while q ≠0
set end of L to q mod 2
set ctr to ctr + 1
set q to q div 2
end repeat
repeat while ctr < 3
set end of L to 0
set ctr to ctr + 1
end repeat
end if
end repeat
return reverse of L as string
else
return -1
end if
end expansionToBase2
# set L to reverse of (characters of "11110001")
on expansionFromBase2(str, base_b)
(*
16, 10 and 8 are valid values for base b.
log expansionFromBase2("11110001", 10) --> 241
log expansionFromBase2("11111010111100", 16) --> 3EBC
log expansionFromBase2("1111", 8) -->17
*)
# we prepare the string.
set L to reverse of (characters of str)
set bstringCount to count L
if base_b = 10 then
set decimalNumber to 0
repeat with i from 1 to bstringCount
set decimalNumber to decimalNumber + (2 ^ (i - 1) * (item i of L))
end repeat
return decimalNumber as integer
else if base_b is 16 then
set hexdigitstring to "0123456789ABCDEF"
repeat while bstringCount mod 4 ≠0
set end of L to 0
set bstringCount to bstringCount + 1
end repeat
set hexnumber to {}
repeat with j from 1 to bstringCount div 4
set hexdigit to 0
repeat with i from 1 to 4
set hexdigit to hexdigit + (2 ^ (i - 1) * (item (i + ((j - 1) * 4)) of L))
end repeat
set end of hexnumber to item (hexdigit + 1) of hexdigitstring
end repeat
return reverse of hexnumber as string
else if base_b is 8 then
repeat while bstringCount mod 3 ≠0
set end of L to 0
set bstringCount to bstringCount + 1
end repeat
set octnumber to {}
repeat with j from 1 to bstringCount div 3
set octdigit to 0
repeat with i from 1 to 3
set octdigit to octdigit + (2 ^ (i - 1) * (item (i + ((j - 1) * 3)) of L))
end repeat
set end of octnumber to (octdigit as integer)
end repeat
return reverse of octnumber as string
else
return -1
end if
end expansionFromBase2
on modularexpansion(b, L, m)
(*
Expands a number on the form b^n mod m, and keeps the remainder small at all times, this way, we can compute with remainders from large numbers, since we reduce them asap. The exponent is represented as a binary string (n)to make this happen. The base is whatever, no need to be 2 or a number base. 245 is an ok base.
3^644 mod 645 shall yield 36
*)
set x to 1
set the_power to b mod m
repeat with i from (count characters of L) to 1 by -1
if character i of L is "1" then
set x to (x * the_power) mod m
end if
set the_power to (the_power * the_power) mod m
end repeat
return x
end modularexpansion
to ceil (n)
-- ceil: "rounds" upwards to closest integer
set i to n div 1
set f to n mod 1
if f > 0 then
set f to 1
else
set f to 0
end if
return (i + f)
end ceil
to floor (n)
-- floor: "truncates" downwards to closest integer
set i to n div 1
if n < 0 then
set f to n mod 1
if f < 0 then
set f to 1
else
set f to 0
end if
return (n - f) div 1
else
return n div 1
end if
end floor
on log2(num)
-- returns the integer of the logarithm of base 2
-- (truncated to closest integer)
-- log2(900) returns 9, and so will log2(512)
if num < 1 then
return -1 -- that is an error that shouldn't happen!
-- it is your responsibility to feed it values ≥ 0.
else
set res to 0
end if
repeat while num ≥ 65536
set num to num div 65536
set res to res + 16
end repeat
repeat while num ≥ 256
set num to num div 256
set res to res + 8
end repeat
repeat while num > 1
set num to num div 2
set res to res + 1
end repeat
return res
end log2
# The two handlers below needs Satimage.osax
on ln2(num)
-- returns the logarithm of 2 as a floating point number
return (ln num) / 0.69314718056
-- Avoids extra ln call, because uses precomputed constant.
end ln2
on lnB(aNumber, base)
return (ln aNumber) / (ln base)
end lnB
-- Combinatorics
on factorial(n)
if n > 170 then error "Factorial: n too large."
# Result greater than 1.797693E+308!
if n < 0 then error "Factorial: n not positive."
set a to 1
repeat with i from 2 to n
set a to a * i
end repeat
return a
end factorial
on P(n, r)
# counts the number of r-permutations in a set of n elements.
set nom to n
repeat with j from (n - 1) to (n - r + 1) by -1
set nom to nom * j
end repeat
return nom
end P
on c(n, r)
# Counts the number of r-combinations in a set of n elements.
# it is also called the binomial coeffecient.
if r = 0 then return 1
set nom to n
repeat with j from (n - 1) to (n - r + 1) by -1
set nom to nom * j
end repeat
set denom to r
repeat with i from (r - 1) to 2 by -1
set denom to denom * i
end repeat
return nom div denom
end c
