fft based polynomial multiplication with biginteger coefficients

I am trying to multiply two polynomial mod n ( n not necessairly a prime) with bigint coefficents I am trying to do that both with python and with boost::multiprecision library.

The problem I am facing is that when multiplying bigint with cos or sin the precision goes down so much that the result is not correct, and when I increase the accuracy using something like sympy or multiprecision float the multiplication becomes very slow. Using NTT has its own problems ( multiplying modular arithmetic modulo two numbers does not make sense in many cases). So, the question is, is there a way to do what I want efficiently?

most of what I can find on FFT based polynomial multiplication does not address this issue

any paper or something else I would be willing to read and implement

This is my python code for fft

def fft( b,  invert=False) :
    a=[mp.mpc(0)]*len(b);
    a[0:len(b)]=b;
    n=len(a);
    if (n == 1):
        return a;
    a0=[mp.mpc(0)]*(n//2);
    a1=[mp.mpc(0)]*(n//2);
    for i in range(n//2) :
        a0[i] = mp.mpc(a[2*i]);
        a1[i] = mp.mpc(a[2*i+1]);
    a0=fft(a0, invert);
    a1=fft(a1, invert);
    factor=[-1,1];
    ang = (mp.mpf(2) * mp.pi / n) * factor[int(invert)];
    w=mp.mpc(1);
    wn=mp.exp(mp.mpc(real=0,imag=ang));
    
    for i in range(n//2) :
        a[i] = a0[i] + w * a1[i];
        a[i + n//2] = a0[i] - w * a1[i];
        if (invert) :
            a[i] /= mp.mpc(2);
            a[i + n//2] /= mp.mpc(2);
        w *= wn;

    
    ret

and my python code for polynomial multiplication

def polyMulFast(p1,p2,n):
    r=len(p1);
    m=2**(r.bit_length());
    p1F=[0]*m;
    p2F=[0]*m;
    p1F[0:len(p1)]=[mp.mpc(x) for x in p1];
    p2F[0:len(p1)]=[mp.mpc(x) for x in p2];
    p1F=fft(p1F);
    p2F=fft(p2F);
    result=[0]*m;
    for i in range(m):
        result[i]=p1F[i]*p2F[i];
    result=fft(result,True);
    res=[0]*r;
    for i in range(len(result)):
        res[i%r]+=int(round(result[i].real));
       # res[i%r]%=n;
    #res[i%len(p1)]%=n;
    return res;

it works well on polynomials with small coefficients,

but once the coefficients are big enough( around 500 bits) it does not work