Of course, for pedagogical purposes, my algorithm is incomplete and can only check whether the point of the interior of the tetrahedron.

If the point lie on the surface of the tetrahedron, for example on the face  ABC , then verify it:

1) Solve the system of  x * A + y * B + z * C = P .

2) If  x + y + z = 1  and  0 <= x <= 1  and  0 <= y <= 1  and  0 <= z <= 1  then the point  P  belongs to the face ABC .

3) If the point does not lie within the tetrahedron and does not lie on it's surface, it lies outside of the tetrahedron.

1) Find the equations of planes bounding the tetrahedron. Let they be F1, F2, F3, F4

2) Take a point lying exactly inside the tetrahedron. It be  K((x1+x2+x3+x4)/4, (y1+y2+y3+y4)/4, (z1+z2+z3+z4)/4)=K(x5, y5, z5)

3) If all the equalities hold  sign(F1(x0, y0, z0))=sign(F1(x5, y5, z5)),   sign(F2(x0, y0, z0))=sign(F2(x5, y5, z5)),  sign(F3(x0, y0, z0))=sign(F3(x5, y5, z5)),  sign(F4(x0, y0, z0))=sign(F4(x5, y5, z5)) , then the point  P(x0,y0,z0)  lies  inside tetrahedron ABCD .

Your equation  5 = floor(4969/x)+floor(208/x)  can be solved exactly by using  simple estimates with a further solution in Maple. If  x <= 800  or  x> = 1000, then obviously no solution because in the first case    floor(4969/x)>=6  and in the second case  floor(4969/x)<5  and  floor(208/x)=0 . But if  800<x<1000  
then floor(208/x)=0  and the equation   5 = floor(4969/x)+floor(208/x)  is equivalent to equation 5 = floor(4969/x) .

solve((5 <= 4969/x and 4969/x < 6);  #  Exact solution 

                     RealRange(Open(4969/6),4969/5)

evalf(%);

                RealRange(Open(828.1666667), 993.8000000)

You can use my procedure  Picture .  See  http://www.mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-

Example (shaded circular segment: center = [2,1], radius= 3, from  -(1/3)*Pi  to Pi radians):

L := [[[2+3*cos(t), 1+3*sin(t)], t = -(1/3)*Pi .. Pi], [[-1, 1], [7/2, 1-3*sqrt(3)*(1/2)]]]:
plots[display](Picture(L, color = green, [color = black, thickness = 2]), scaling = constrained, view = [-1 .. 6, -2 .. 5]);

Another example:

L := [[3, t = -(1/2)*Pi .. (1/2)*Pi, polar], [[-4+5*cos(t), 5*sin(t)], t = arctan(3/4) .. -arctan(3/4)]]:
plots[display](Picture(L, color = yellow, [color = brown, thickness = 3]), scaling = constrained, view = [-1 .. 4, -4 .. 4]);

Your equation corresponds to the geometric image - the surface of the second order. If this surface is bounded on at least one dimension, the procedure  Isolve  finds all integer points on the surface.

Formal arguments of procedure:  Equation  is the equation of the surface,  Range  is the range of variation of a variable in the format  variable = a .. b. The easiest way to check boundness of the surface is the construction of it's plot. If the surface is unbounded, the procedure finds all the solutions in a specified range.

restart;

Isolve:=proc(Equation, Range)

local L, i, Sol;

L:=[]:

   for i from lhs(rhs(Range)) to rhs(rhs(Range)) do

      Sol:=isolve(eval(Equation, lhs(Range)=i));

      L:=[op(L),seq([op(Sol[k]),z=i], k=1..nops([Sol]))];

   od:

L;

end proc:

Your example:

plots[implicitplot3d](x^2 +y^2 +2*z^2 +x*y -y*z=3, x=-3..3, y=-3..3, z=-3..3, axes=normal, numpoints=10000);

+

We have an ellipsoid.

Isolve(x^2 +y^2 +2*z^2 +x*y -y*z=3, z=-2..2);

[[x = -1, y = 0, z = -1], [x = 1, y = -2, z = -1], [x = 1, y = 0, z = -1], [x = -2, y = 1, z = 0], [x = -1, y = -1, z = 0], [x = -1, y = 2, z = 0], [x = 1, y = -2, z = 0], [x = 1, y = 1, z = 0], [x = 2, y = -1, z = 0], [x = -1, y = 0, z = 1], [x = -1, y = 2, z = 1], [x = 1, y = 0, z = 1]]

Finding of the all integer points on unbounded surface  x^2+y^2-z^2=1  (hyperboloid) in the range  z=-100..100 :

L:=Isolve(x^2 +y^2 -z^2=1, z=-100..100):
nops(L);    # The total number of solutions
 L[1000..1020];     # A few solutions

                                                                                  2892

[[x = 31, y = 22, z = -38], [x = 34, y = -17, z = -38], [x = 34, y = 17, z = -38], [x = 38, y = -1, z = -38], [x = 38, y = 1, z = -38], [x = -37, y = -1, z = -37], [x = -37, y = 1, z = -37], [x = -29, y = -23, z = -37], [x = -29, y = 23, z = -37], [x = -23, y = -29, z = -37], [x = -23, y = 29, z = -37], [x = -1, y = -37, z = -37], [x = -1, y = 37, z = -37], [x = 1, y = -37, z = -37], [x = 1, y = 37, z = -37], [x = 23, y = -29, z = -37], [x = 23, y = 29, z = -37], [x = 29, y = -23, z = -37], [x = 29, y = 23, z = -37], [x = 37, y = -1, z = -37], [x = 37, y = 1, z = -37]]

I think that finding  all  integer points on the unbounded surface of the second order, in general case,  is not easy problem. For example,  Mathematica does not cope with this task.

We use the formula from wiki  http://en.wikipedia.org/wiki/Surface_of_revolution

sol:=solve(eval(x^2/20+y^2/13+z^2/20=1, z=0),  x);

2*Pi*int(sol[1]*sqrt(1+diff(sol[1], y)^2), y=-sqrt(13)..sqrt(13));

S:=combine(%);

evalf(S);

You need  an other  software, for example  

http://umsolver.com/usa/mathematics/free-algebra-equation-solver/ 

interface(rtablesize=infinity):

MyHandler := proc(operator,operands,default_value)

   NumericStatus( division_by_zero = false );

   return infinity;

end proc: 

NumericEventHandler(division_by_zero=MyHandler):

Matrix([[Rad, Degree, sin(phi), cos(phi), tan(phi)], [``,``,``,``,``], seq([seq(j, j=[i, i*180/Pi, convert(sin(i), radical), convert(cos(i), radical), convert(tan(i), radical)])], i=[seq(k*Pi/12, k=0..24)])]);

Fragment of output:

Acer's procedure  MyHandler  was used.

Exact periodic extension of function  g  to the whole real axis can be constructed as follows:

g := t->piecewise(0 <= t and t <= 0.2e-1, 122.63, 0.2e-1 < t and t < .287, 0, .287 <= t and t <= .307, 122.63, .307 < t and t < .335, 0):

R := t->(t/0.335-floor(t/0.335))*0.335:

G := unapply((g@R)(t), t);

Example:

plot(G(t), t = -2*.335 .. 3*.335);

If I understand your question, you want the maximum element in the matrix put in the needed position  and use as a pivoting element. This can be done by interchanging of rows and/or columns. The procedure  MaxPivot sets maximum element of matrix  A  at the position, defined by the list  L:

MaxPivot:=proc(A::Matrix, L::list)

local m, n, M, i, j;

uses LinearAlgebra;

m:=RowDimension(A); n:=ColumnDimension(A);

M:=[[0,0], 0];

for i to m do

for j to n do

if abs(A[i,j])>M[2] then M:=[[i,j], abs(A[i,j])] fi;

od; od;

RowOperation(A, [L[1], M[1,1]]);

ColumnOperation(%, [L[2], M[1,2]]);

end:

Example:

A:=<<1 | 2 | 4> , <5 | 6 | 9> , <-1 | 3 | 4>>:
MaxPivot(A, [1,1]);

The basic idea: if  a<b  that of the inequality  f(x)<g(x)  imply that  int(f(x), x=a..b)<int(g(x), x=a..b) . Sufficient to prove the  original limit for x> 0 .

For  0<x<Pi  consider the obvious double inequality  0<sin(x)<x  and integrate  it several times from 0 to x , expressing everytime  sin(x)  or  cos(x) . This can easily be automated in Maple:

L:=[0<sin(x), sin(x)<x]:

for k to 4 do

L:=[seq(op(solve(map(int, L[i], x=0..x), `if`(k::odd, cos(x), sin(x)))), i=1..2)];

od;  

From the last inequality we express  A = (x-sin (x)) / x^3 :

seq(solve(subs(sin(x)=x-A*x^3,L[i]), {A}) , i=1..2) assuming x^3>0;

So we have a double eniquality  -1/120*x^2+1/6 < A <1/6  , from which it obviously follows that

limit(A, x=0)=1/6

int(sqrt(diff(sqrt(20)*cos(t),t)^2+diff(sqrt(13)*sin(t),t)^2), t=0..2*Pi);
evalf(%);

                                     8*5^(1/2)*EllipticE(1/10*35^(1/2))
                                                      25.44987267

Did you remember to put a semicolon or a colon at the end of the first line?

NumberOfDigits:=proc(N)

local X, n;

X[0]:=0; X[1]:=1; X[2]:=2;

for n from 3 to N do

X[n]:=n*(X[n-3]+X[n-2]+X[n-1]);

od;

convert(X[N], string);

length(%);

end proc;

Your example:

NumberOfDigits(2013);

                  5782

Position:=proc(x::anything, L::list)

local p;

if member(x, L, 'p') then p else 0 fi;

end proc:

Examples:

Position(2, [3,8,2,7]);

Position(5, [3,8,2,7]);

                    3

                    0