Some Integer handlers in AppleScript

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

Updated post above

More efficient versions of gcd and lcm, added isCoprime hander.

Updated post #1

Added division algorithm, as used in number therory for expressing numbers as linear combinations of the quotient and divisor, and the rest, which is defined as a positive quantity, and not as implemented by Applescripts modulus operator.

Updated post #1.

Added divides handler, and modulus handler. Bugfix div_alg.

Nice to have you back McUser :cool:

Thanks DJ, I have been quitting smoking, hahaha, I haven’t done what I used to do while I smoked, but now a half year has passed so. . .

I had to rework the div_alg some more, to make it return the correct sign of the rest in all cases, now it should return both the correct value, and the correct sign, which should be positive in all cases. I have also corrected the usage of divisor and dividend . . divisor is what partitions a number, and the dividend is the number that is partitioned.

Added handlers to the script in post #1

relprime handler, which tests if all the positive integers in a list are relatively prime (gdc(int1,int2)=1 ), expansionToBase2, expansionFromBase2, which conversts from base10 to base16,base2 and base 8. Base2 are used as an intermediary base for all conversions.

Added xgcd Extended euclid for solving linear congruences, it is a really brilliant handler, that expresses the relationship between the gcd between two coprimes as a linear combination: 1 = sa + tb, it is great to have when doing modular arithmetic! :slight_smile:

Hello.

I have added some handlers to the script in post #1: first modular_expansion, which I had forgotten to post, then ceil, floor, log2 and ln2. The latter should all be fast, if not as fast as possible. You’ll need Satimage.Osax for the ln2 handler. I see clearly that I should split this collection into a “prime” library, and add the rest to an “integer/number” library. :slight_smile:

Hello.

I added an lnB(base, aNumber) handler, that returns the logarithm of a number with respect to the given base, as a real.

I use the naming standard that handlers that returns logarithms as ints, are prefixed with log, and those that returns floats are prefixed with ln.

Hello.

I have added handlers for the binomial coeffecient, the number of r-permutations, and the factorial handler.


log c(4, 3)
# 4
log p(4,3)
# 24
log factorial(4)
# 24

Hello.

I added treatment of the case that the number of elements in a combination was zero, thereby removing a bug.

I’m sorry for any inconvenience.