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