Algorithms employed:

sin(x): 
First, the task is converted to the calculation of the sine for x from 0 to p/2. For x<0.88 the Taylor expansion up to x9, centered at 0 is used. For x³0.88 the expansion up to x8, centered at 1 (shifted cosine expansion) is employed. The boundary 0.88 is chosen so that the errors of both approximations are roughly the same at this point.

cos(x) = sin(x-p/2); 
tan(x) = sin(x) / cos(x).

arcsin(x):
If |x|<0.51 the rational Pade approximation is used directly (in the nominator up to x7, in the denominator up to x6). Else sign(x) × {p/2 - 2 arcsin [ sqrt( (1-x) / 2 )]} is calculated where arcsin is evaluated by the same rational expansion. In fact, the relation cos2(x) = [1 + cos(2x)] / 2 is employed in order to keep the absolute value of the argument of the expansion as low as possible (the smaller argument the smaller error of the approximation).

arccos(x) = p/2 - arcsin(x); 
arctan(x) = arcsin[x / sqrt(1 + x2)].

ln(x):
First, an exponent is handled separately (it is done by string operations) - the determination of the exponent, ex, leads to the "zero" estimation, y0 = ex × ln(10). Mantissa, m, should be between 1 and 9.999... It is multiplied by 0.3 (X = m × 0.3) to transform the range to 0.3 - 3 [the value -ln(0.3) must be added to the estimation]. Then the integral of 1/x from 1 to X is evaluated by a Gauss quadrature [note that ò1/x = ln(x)]. This leads to the first correction, y1 = A / (B-C) + A / (B+C), where A = X-1, B = X+1, C = A / sqrt(3). Then the same integral is calculated for X' = X / exp(y1). This gives the second, more precise correction, y2 (the precision of the used approximation of the integral is "good enough" if 0.3<X<3 and "excellent" if X is very close to 1). Finally, ln(x) = y0 + ln(0.3) + y1 + y2.

log(x) = ln(x) / ln(10).

xy and sqrt(x) are WMLScript-library functions.