1# The elementary functions are:
2# polynomial operations, nth roots,
3# exp, log, and everything you can get from these
4# In particular, it contains the trig functions and the hyperbolic functions
5
6
7# These are most relevant here.
8SetHelp("rad2deg","functions","Convert radians to degrees");
9function rad2deg(x) = (
10	if(IsMatrix(x)) then
11		return ApplyOverMatrix(x,rad2deg)
12	else if(not IsValue(x)) then
13		(error("rad2deg: argument not a value");bailout);
14	(x*180)/pi
15);
16
17SetHelp("deg2rad", "functions", "Convert degrees to radians");
18function deg2rad(x) = (
19	if(IsMatrix(x)) then
20		return ApplyOverMatrix(x,deg2rad)
21	else if(not IsValue(x)) then
22		(error("deg2rad: argument not a value");bailout);
23	(x*pi)/180
24);
25
26#FIXME: these may not deal well with zero values. (should treat 0 correctly now)
27#FIXME: these should be further double checked
28
29SetHelp("asin","trigonometry","The arcsin (inverse sin) function");
30function asin(x) = (
31	if(IsMatrix(x)) then
32		return ApplyOverMatrix(x,asin)
33	else if(not IsValue(x)) then
34		(error("asin: argument not a value");bailout);
35
36	if x==1 then
37		pi/2
38	else if x==-1 then
39		-pi/2
40	else if IsReal(x) and -1 < x < 1 then
41		atan(x/sqrt(1-x^2))
42	else
43		(-1i)*ln(1i*x+sqrt(abs(1-x^2))*exp((1i/2)*Arg(1-x^2)))
44);
45arcsin = asin
46SetHelpAlias ("asin", "arcsin");
47
48SetHelp("asinh","trigonometry","The arcsinh (inverse sinh) function");
49function asinh(x) = (
50	if(IsMatrix(x)) then
51		return ApplyOverMatrix(x,asinh)
52	else if(not IsValue(x)) then
53		(error("asinh: argument not a value");bailout);
54	if IsReal(x) then
55		ln(x+sqrt((x^2)+1))
56	else
57		ln(x+sqrt(abs(1+x^2))*exp((1i/2)*Arg(1+x^2)))
58);
59arcsinh = asinh
60SetHelpAlias ("asinh", "arcsinh");
61
62SetHelp("acos","trigonometry","The arccos (inverse cos) function");
63function acos(x) = (
64	if(IsMatrix(x)) then
65		return ApplyOverMatrix(x,acos)
66	else if(not IsValue(x)) then
67		(error("acos: argument not a value");bailout);
68
69	if x==1 then
70		0
71	else if x==-1 then
72		pi
73	else if x==0 then
74		pi/2
75	else if IsReal(x) and -1 < x < 1 then
76		atan(sqrt(1-x^2)/x)+(if x>0 then 0 else pi)
77	else
78		(-1i)*ln(x+1i*sqrt(abs(1-x^2))*exp((1i/2)*Arg(1-x^2)))
79);
80arccos = acos
81SetHelpAlias ("acos", "arccos");
82
83SetHelp("acosh","trigonometry","The arccosh (inverse cosh) function");
84function acosh(x) = (
85	if(IsMatrix(x)) then
86		return ApplyOverMatrix(x,acosh)
87	else if(not IsValue(x)) then
88		(error("acosh: argument not a value");bailout);
89	if IsReal(x) and -1 <= x <= 1 then
90		ln(x+sqrt((x^2)-1))
91	else
92		ln(x+sqrt(abs(x^2-1))*exp((1i/2)*Arg(x^2-1)))
93);
94arccosh = acosh
95SetHelpAlias ("acosh", "arccosh");
96
97SetHelp("cot","trigonometry","The cotangent function");
98function cot(x) = (
99	if(IsMatrix(x)) then
100		return ApplyOverMatrix(x,cot)
101	else if(not IsValue(x)) then
102		(error("cot: argument not a value");bailout);
103	1/tan(x)
104);
105
106SetHelp("coth","trigonometry","The hyperbolic cotangent function");
107function coth(x) = (
108	if(IsMatrix(x)) then
109		return ApplyOverMatrix(x,coth)
110	else if(not IsValue(x)) then
111		(error("coth: argument not a value");bailout);
112	1/tanh(x)
113);
114
115SetHelp("acot","trigonometry","The arccot (inverse cot) function");
116function acot(x) = (
117	if(IsMatrix(x)) then
118		return ApplyOverMatrix(x,acot)
119	else if(not IsValue(x)) then
120		(error("acot: argument not a value");bailout);
121
122	#atan(1/x)
123	if IsReal(x) then
124		Re(ln((x+1i)/(x-1i))/2i)
125	else
126		ln((x+1i)/(x-1i))/2i
127);
128arccot = acot
129SetHelpAlias ("acot", "arccot");
130
131SetHelp("acoth","trigonometry","The arccoth (inverse coth) function");
132function acoth(x) = (
133	if(IsMatrix(x)) then
134		return ApplyOverMatrix(x,acoth)
135	else if(not IsValue(x)) then
136		(error("acoth: argument not a value");bailout);
137	#atanh(1/x)
138	ln((x+1)/(x-1))/2
139);
140arccoth = acoth
141SetHelpAlias ("acoth", "arccoth");
142
143SetHelp("tanh","trigonometry","The hyperbolic tangent function");
144function tanh(x) = (
145	if(IsMatrix(x)) then
146		return ApplyOverMatrix(x,tanh)
147	else if(not IsValue(x)) then
148		(error("tanh: argument not a value");bailout);
149	sinh(x)/cosh(x)
150);
151
152SetHelp("atanh","trigonometry","The arctanh (inverse tanh) function");
153function atanh(x) = (
154	if(IsMatrix(x)) then
155		return ApplyOverMatrix(x,atanh)
156	else if(not IsValue(x)) then
157		(error("atanh: argument not a value");bailout);
158	ln((1+x)/(1-x))/2
159);
160arctanh = atanh
161SetHelpAlias ("atanh", "arctanh");
162
163SetHelp("csc","trigonometry","The cosecant function");
164function csc(x) = (
165	if(IsMatrix(x)) then
166		return ApplyOverMatrix(x,csc)
167	else if(not IsValue(x)) then
168		(error("csc: argument not a value");bailout);
169	1/sin(x)
170);
171
172SetHelp("csch","trigonometry","The hyperbolic cosecant function");
173function csch(x) = (
174	if(IsMatrix(x)) then
175		return ApplyOverMatrix(x,csch)
176	else if(not IsValue(x)) then
177		(error("csch: argument not a value");bailout);
178	1/sinh(x)
179);
180
181SetHelp("acsc","trigonometry","The inverse cosecant function");
182function acsc(x) = (
183	if(IsMatrix(x)) then
184		return ApplyOverMatrix(x,acsc)
185	else if(not IsValue(x)) then
186		(error("acsch: argument not a value");bailout);
187	asin(1/x)
188);
189arccsc = acsc
190SetHelpAlias ("acsc", "arccsc");
191
192SetHelp("acsch","trigonometry","The inverse hyperbolic cosecant function");
193function acsch(x) = (
194	if(IsMatrix(x)) then
195		return ApplyOverMatrix(x,acsch)
196	else if(not IsValue(x)) then
197		(error("acsc: argument not a value");bailout);
198	asinh(1/x)
199);
200arccsch = acsch
201SetHelpAlias ("acsch", "arccsch");
202
203SetHelp("sec","trigonometry","The secant function");
204function sec(x) = (
205	if(IsMatrix(x)) then
206		return ApplyOverMatrix(x,sec)
207	else if(not IsValue(x)) then
208		(error("sec: argument not a value");bailout);
209	1/cos(x)
210);
211
212SetHelp("sech","trigonometry","The hyperbolic secant function");
213function sech(x) = (
214	if(IsMatrix(x)) then
215		return ApplyOverMatrix(x,sech)
216	else if(not IsValue(x)) then
217		(error("sech: argument not a value");bailout);
218	1/cosh(x)
219);
220
221SetHelp("asec","trigonometry","The inverse secant function");
222function asec(x) = (
223	if(IsMatrix(x)) then
224		return ApplyOverMatrix(x,asec)
225	else if(not IsValue(x)) then
226		(error("asec: argument not a value");bailout);
227	acos(1/x)
228);
229arcsec = asec
230SetHelpAlias ("asec", "arcsec");
231
232SetHelp("asech","trigonometry","The inverse hyperbolic secant function");
233function asech(x) = (
234	if(IsMatrix(x)) then
235		return ApplyOverMatrix(x,asech)
236	else if(not IsValue(x)) then
237		(error("asech: argument not a value");bailout);
238	acosh(1/x)
239);
240arcsech = asech
241SetHelpAlias ("asech", "arcsech");
242
243SetHelp("log","numeric","Logarithm of any base (calls DiscreteLog if in modulo mode), if base is not given, e is used");
244function log(x,b...) = (
245	m = GetCurrentModulo ();
246	if not IsNull (m) then (
247		if IsNull (b) or elements(b) > 1 then
248			(error("log (discrete): wrong number of arguments");bailout);
249		return DiscreteLog (x, b@(1), m)
250	);
251
252	if IsNull (b) then
253		return ln(x)
254	else if elements(b) > 1 then
255		(error("log: too many arguments");bailout);
256	base = b@(1);
257
258	if IsMatrix(x) or IsMatrix(base) then
259		return ApplyOverMatrix2(x,base,log)
260	else if(not IsValue(x) or not IsValue(base)) then
261		(error("log: arguments not values");bailout);
262	ln(x)/ln(base)
263);
264
265# This is still used for complex values in the hacky computation in funclib.c
266# although for real values we use MPFR
267parameter ErrorFunctionTolerance = 10.0^(-10);
268SetHelp ("ErrorFunctionTolerance", "parameters", "Tolerance of the ErrorFunction (used for complex values only)")
269
270# This is actually done for complex values inside funclib.c
271# This should as some point be replaced by a proper version of erf
272#SetHelp("ErrorFunction","functions","The error function, 2/sqrt(pi) * int_0^x e^(-t^2) dt")
273#function ErrorFunction (x) = (
274#	if IsMatrix (x) then
275#		return ApplyOverMatrix (x, ErrorFunction)
276#	else if not IsValue (x) then
277#		(error("ErrorFunction: argument not a value");bailout);
278#	twosqrtpi = 2/sqrt(pi);
279#	a = 1;
280#	s = 0;
281#	n = 0;
282#	f = 1;
283#	xx = x;
284#	xsq = x^2;
285#	do (
286#		t = xx * a * twosqrtpi;
287#		s = s + t;
288#		increment n;
289#		f = f * n;
290#		a = ((-1)^n) / (((2*n)+1) * f);
291#		xx = xx * xsq
292#	) while (|t| > ErrorFunctionTolerance);
293#	s
294#);
295#erf = ErrorFunction
296#SetHelpAlias ("ErrorFunction", "erf");
297
298#FIXME: Should probably be in a separate source file
299SetHelp("NewtonsMethodPoly","polynomial","Attempt to find a root of a polynomial using Newton's method, returning after two successive values are within epsilon or after maxn tries (then returns null)")
300function NewtonsMethodPoly(poly,guess,epsilon,maxn) = (
301	pf := PolyToFunction (poly);
302	pdf := PolyToFunction (PolyDerivative (poly));
303	guess := float(guess);
304	for n=1 to maxn do (
305		pdfg := pdf(guess);
306		if pdfg == 0.0 then (
307			error ("NewtonsMethodPoly: division by zero");
308			bailout
309		);
310		guessn := guess - pf(guess)/pdfg;
311		if |guessn-guess| <= epsilon then
312			return guessn;
313		guess := guessn
314	);
315	null
316)
317