Math functions again

Hallenstal · April 6, 2023, 8:00pm

I have amused myself with writing math functions in AppleScript. What is needed more than sin, cos, tan, asin, acos, atan, ln and exp?
here’s some code

on cos(a as real)
	if a < 0 then return cos(-a)
	if a = 0 then return 1
	if a > 2 * pi then return cos(a - ((a / 2 / pi) div 1) * 2 * pi)
	if a = pi / 2 or a = 3 * pi / 2 then return 0
	if a > pi / 2 and a < 3 * pi / 2 then return -(cos(pi - a))
	set x to sin(a)
	return (1 - x * x) ^ 0.5
end cos

on acos(a as real)
	if a > 1.0 then return "ERROR"
	return pi / 2.0 - asin(a)
end acos

on sin(a as real)
	if a = 0 then return 0
	if a < 0 then return -(sin(-a))
	if a > 2 * pi then set a to a - ((a / 2 / pi) div 1) * 2 * pi
	if a > pi then return -(sin(2 * pi - a))
	if a > pi / 2 then set a to pi - a
	set {fact, x_p, f_x} to {1, a, a}
	repeat with n from 3 to 23 by 2
		set fact to (-fact) * n * (n - 1)
		set x_p to x_p * a * a
		set f_x to f_x + x_p / fact
	end repeat
	return f_x
end sin

on asin(a as real)
	if a > 1.0 then return "ERROR"
	if a = 0 then return a
	if a = 1.0 then return pi / 2
	if a < 0 then return asin(-a)
	if a > 0.707106781187 then return 1.570796326795 - asin((1 - a * a) ^ (0.5))
	set {x, f_x, tlj, nmn_1} to {a, a, 1, 1}
	if a ≤ 0.4 then
		set M to 18
	else if a ≤ 0.5 then
		set M to 23
	else if a ≤ 0.6 then
		set M to 30
	else
		set M to 43
	end if
	repeat with n from 1 to M
		set tlj to tlj * n * (2 * n - 1)
		set nmn_1 to nmn_1 * 2 * n * n
		set x to x * a * a
		set f_x to f_x + x * tlj / (nmn_1 * (2 * n + 1))
	end repeat
	return f_x
end asin

on tan(a as real)
	if ((a + pi / 2) / pi mod 1) = 0 then return "ERROR"
	set x to sin(a)
	return x / ((1 - x * x) ^ 0.5)
end tan

on atan(a as real)
	if a < 0 then return -(atan1(-a))
	if a > 1 then return pi / 2 - atan1(1 / a)
	return asin(a / ((a * a + 1) ^ 0.5))
end atan


on exp(a as real)
	return 2.718281828459 ^ a
end exp
on ln(a as real)
	if a ≤ 0 then return "ERROR"
	if a = 1.0 then return 0.0
	if a < 1.0E+4 and a > 1 then
		set exp10_1 to -4.0
		set a to 1.0E+4 * a
		set off to offset of "+" in (a as text)
	else if a < 1 and a ≥ 1.0E-4 then
		set exp10_1 to 4.0
		set a to 1.0E-4 * a
	else
		set exp10_1 to 0
	end if
	set a_t to a as text
	set off to offset of "×10" in a_t
	set exp10_2 to (text (off + 4) thru (count of a_t) of a_t) as real
	set a to (a / (10 ^ exp10_2)) -- alternaltiv use coercing	
	if a ≥ 4.481689070338 then
		set exp to 2.0
		set a to a / 7.389056098931
	else if a > 1.6487212707 then
		set exp to 1.0
		set a to a / 2.718281828459
	else
		set exp to 0.0
	end if
	set x to (a - 1) / (a + 1)
	set f_x to x * 2.0 + x * x * x * (0.666666666667 + x * x * (0.4 + x * x * (0.285714285714 + x * x * (0.222222222222 + x * x * (0.181818181818 + x * x * (0.153846153846 + x * x * (0.133333333333 + x * x * (0.117647058824 + x * x * (0.105263157895 + x * x * (0.095238095238 + x * x * 0.086956521739))))))))))
	return (exp10_2 + exp10_1) * 2.302585092994 + exp + f_x
end ln

enjoy:-)

Mockman · April 6, 2023, 9:50pm

Where is atan1 defined? I’m getting a script error (-1708).

Hallenstal · April 7, 2023, 9:20am

Yes, I tried 2 versions and forgot to change:-(
Maclaurin series of asin converges faster than atan so when that solution.
the correct atan should be

on atan(a as real)
	if a < 0 then return -(atan(-a))
	if a > 1 then return pi / 2 - atan(1 / a)
	return asin(a / ((a * a + 1) ^ 0.5))
end atan

Nigel_Garvey · April 7, 2023, 10:27am

Hallenstal:

on asin(a as real)
	if a > 1.0 then return "ERROR"
	if a = 0 then return a
	if a = 1.0 then return pi / 2
	if a < 0 then return asin(-a)
	if a > 0.707106781187 then return 1.570796326795 - asin((1 - a * a) ^ (0.5))
	set {x, f_x, tlj, nmn_1} to {a, a, 1, 1}
	if a ≤ 0.4 then
		set M to 18
	else if a ≤ 0.5 then
		set M to 23
	else if a ≤ 0.6 then
		set M to 30
	else
		set M to 43
	end if
	repeat with n from 1 to M
		set tlj to tlj * n * (2 * n - 1)
		set nmn_1 to nmn_1 * 2 * n * n
		set x to x * a * a
		set f_x to f_x + x * tlj / (nmn_1 * (2 * n + 1))
	end repeat
	return f_x
end asin

Just for amusement and interest, since it only makes any noticeable difference when high-density number crunching’s involved:

It takes very slightly less time to test if something’s greater than something else than to test if it’s less than or equal to (or not greater than) that something else. Similarly less than vs. greater than or equal to. So if it’s sensibly possible to arrange for the corresponding faster test to be done instead, you get a slight speed improvement if the test has to be done thousands of times. For instance, repeat until (x > 1000) rather than repeat while (x ≤ 1000). Or if you’re really desperate:

if (a > b) then
else
	-- Do what has to be done if a ≤ b.
end if

In my tinkering with Hallenstal’s asin() handler below, I’ve reversed the order of tests in the if … else block to allow the use of > instead of ≤. This makes the tests themselves slightly faster, but of course it also changes the order of the numbers tested — and thus the number of tests that have to be done before there’s a match in any particular case.

It takes very slightly less time to add a number to itself than to multiply it by 2.
In the asin() handler, a * a only needs to be calculated once, before the repeat, and 2 * n (or n + n) only once during it.

atan(0.5)

on atan(a as real)
	if a < 0 then return -(atan(-a))
	if a > 1 then return pi / 2 - atan(1 / a)
	return asin(a / ((a * a + 1) ^ 0.5))
end atan

on asin(a as real)
	if a > 1.0 then return "ERROR"
	if a = 0.0 then return a
	if a = 1.0 then return pi / 2
	if a < 0.0 then return asin(-a)
	if a > 0.707106781187 then return 1.570796326795 - asin((1 - a * a) ^ (0.5))
	set {x, f_x, aSquared, tlj, nmn_1} to {a, a, a * a, 1, 1}
	if (a > 0.6) then
		set M to 43
	else if (a > 0.5) then
		set M to 30
	else if (a > 0.4) then
		set M to 23
	else
		set M to 18
	end if
	repeat with n from 1 to M
		tell (n + n)
			set tlj to tlj * n * (it - 1)
			set nmn_1 to nmn_1 * it * n
			set x to x * aSquared
			set f_x to f_x + x * tlj / (nmn_1 * (it + 1))
		end tell
	end repeat
	return f_x
end asin

Hallenstal · April 7, 2023, 12:52pm

of course there were other errors:-( (in tan function)
proves testing is needed. Also using pi since it is better accuracy than a decimal number

on cos(a as real)
	if a < 0 then return cos(-a)
	if a = 0 then return 1
	if a > 2 * pi then return cos(a - ((a / 2 / pi) div 1) * 2 * pi)
	if a = pi / 2 or a = 3 * pi / 2 then return 0
	if a > pi / 2 and a < 3 * pi / 2 then return -(cos(pi - a))
	set x to sin(a)
	return (1 - x * x) ^ 0.5
end cos

on acos(a as real)
	if a > 1.0 then return "ERROR"
	return pi / 2.0 - asin(a)
end acos

on sin(a as real)
	if a = 0 then return 0
	if a < 0 then return -(sin(-a))
	if a > 2 * pi then set a to a - ((a / 2 / pi) div 1) * 2 * pi
	if a > pi then return -(sin(2 * pi - a))
	if a > pi / 2 then set a to pi - a
	set {fact, x_p, f_x} to {1, a, a}
	repeat with n from 3 to 23 by 2
		set fact to (-fact) * n * (n - 1)
		set x_p to x_p * a * a
		set f_x to f_x + x_p / fact
	end repeat
	return f_x
end sin

on asin(a as real)
	if a > 1.0 then return "ERROR"
	if a = 0 then return a
	if a = 1.0 then return pi / 2
	if a < 0 then return -(asin(-a))
	if a > 0.707106781187 then return pi / 2 - asin((1 - a * a) ^ (0.5))
	set {x, f_x, tlj, nmn_1} to {a, a, 1, 1}
	if a > 0.6 then
		set M to 43
	else if a > 0.5 then
		set M to 30
	else if a > 0.4 then
		set M to 23
	else
		set M to 18
	end if
	repeat with n from 1 to M
		set tlj to tlj * n * (2 * n - 1)
		set nmn_1 to nmn_1 * 2 * n * n
		set x to x * a * a
		set f_x to f_x + x * tlj / (nmn_1 * (2 * n + 1))
	end repeat
	return f_x
end asin

on tan(a as real)
	if ((a + pi / 2) / pi mod 1) = 0 then return "ERROR"
	set a to ((a / pi) mod 1) * pi
	if a < -pi / 2 then
		set x to (sin(a + pi))
	else if a > -pi / 2 and a < pi / 2 then
		set x to sin(a)
	else
		set x to sin(a - pi)
	end if
	return x / ((1 - x * x) ^ 0.5)
end tan

on atan(a as real)
	if a < 0 then return -(atan(-a))
	if a > 1 then return pi / 2 - atan(1 / a)
	return asin(a / ((a * a + 1) ^ 0.5))
end atan


on exp1(a as real)
	return 2.718281828459 ^ a
end exp1

on ln(a as real)
	if a ≤ 0 then return "ERROR"
	if a = 1.0 then return 0.0
	if a < 1.0E+4 and a > 1 then
		set exp10_1 to -4.0
		set a to 1.0E+4 * a
		set off to offset of "+" in (a as text)
	else if a < 1 and a ≥ 1.0E-4 then
		set exp10_1 to 4.0
		set a to 1.0E-4 * a
	else
		set exp10_1 to 0
	end if
	set a_t to a as text
	set off to offset of "×10" in a_t
	set exp10_2 to (text (off + 4) thru (count of a_t) of a_t) as real
	set a to (a / (10 ^ exp10_2)) -- alternaltiv use coercing	
	if a ≥ 4.481689070338 then
		set exp to 2.0
		set a to a / 7.389056098931
	else if a > 1.6487212707 then
		set exp to 1.0
		set a to a / 2.718281828459
	else
		set exp to 0.0
	end if
	set x to (a - 1) / (a + 1)
	set f_x to x * 2.0 + x * x * x * (0.666666666667 + x * x * (0.4 + x * x * (0.285714285714 + x * x * (0.222222222222 + x * x * (0.181818181818 + x * x * (0.153846153846 + x * x * (0.133333333333 + x * x * (0.117647058824 + x * x * (0.105263157895 + x * x * (0.095238095238 + x * x * 0.086956521739))))))))))
	return (exp10_2 + exp10_1) * 2.302585092994 + exp + f_x
end ln

Mockman · April 7, 2023, 2:21pm

What about an atan2 function? With «cough» degrees rather than radians.

Hallenstal · April 7, 2023, 3:16pm

my understanding of atan2 is that it would be

on atan2(a as real)
   set x to tan(a)
   if x<0 then
      return x+pi
   else
      return x
   end if
end atan2

You could do own functions rename atan2 to _atan2 and have
atan2(x, unit); if unit=“deg” then return _atan2(x)*180/pi else return _atan2(x)

BR

Mockman · April 7, 2023, 8:09pm

You’re right. And actually… you can simply enter the dimensions on each side of a division operator and let it figure out the ratio.

atan (4 / 3) -- in radians
--> 0.927295218002

(180 / pi) * (atan (4 / 3)) -- in degrees
--> 53.130102354156
(180 / pi) * (atan (3 / 4))
--- 36.869897645844

Hallenstal · April 8, 2023, 8:51am

Ferik71,
I think it has to do with the decimal precision of what can be entered and what AppleScript can use in calculations
Checked against excel and the error is < 10^-12 typically inbyte order of abs(10^-13) which in my mind acceptabelt;-)

Piyomaru · April 8, 2023, 2:02pm

calcLibAS Script Library has 36 math functions.
http://piyocast.com/as/asinyaye

Hallenstal · April 8, 2023, 8:55pm

I assume there are some performance penalties going outside of AppleScript