Orthogonal signal generator using integer arithmetic
I managed to find an answer to my problem.

The problem is with rounding in integer division. For example, -4/3 will be rounded to -1, just as -5/3. Because of this, the integration error is constantly accumulated. Instead of explaining, here I give a code how to fix this.

Code declaration.

// Global variables
int w = 643398L; // (2*PI)*50*2048
int i1u1=0, i1u2=0, i1u3=0;
int i2u1=0, i2u2=0, i2u3=0;
short i1y1=0, i2y1=0;

// Numerical integrator implementation
short NumInt3rd(int u1, int u2, int u3, short y1) {
    int iu = 23*u1-16*u2+5*u3;
    int iy = 240000L*y1;
    short y0 = (((iu+iy)&gt;&gt;12)*2237+65536)&gt;&gt;17;
    return y0;
}


Main function.

// Get voltage measurements (-2048 to +2048)
short vg = (short) Input(0);

// Numerical integrators, 3rd order
short vg_x = NumInt3rd(i1u1,i1u2,i1u3,i1y1);
short vg_y = NumInt3rd(i2u1,i2u2,i2u3,i2y1);

// Downsample frequency to prevent overflow
short wb = w&gt;&gt;5;

// Update integrator #1 states
i1u3 = i1u2;
i1u2 = i1u1;
i1u1 = ((int)((vg-vg_x)-vg_y)*wb)&gt;&gt;6;
i1y1 = vg_x;

// Update integrator #2 states
i2u3 = i2u2;
i2u2 = i2u1;
i2u1 = ((int)vg_x*wb)&gt;&gt;6;
i2y1 = vg_y;


It should be noted that instead of using integer division, which is very expensive in terms of required number of instruction cycles, I rather use bit shift operation. For example, 1/240000 can be well approximated as 2237/2^29 - the approximation accuracy is 0.00169277%.

And here is the sine wave estimation error. As you can see, the results are much better now.