int numroots(int np, poly sseq[], int at[])
    {
                    int	atposinf, atneginf, s;
                    double	f, lf;

                    atposinf = atneginf = 0;


            /*
             * changes at positive infinity
             */
            lf = sseq[0].coef[sseq[0].ord];

            for ( s = 1; s <= np; s++ ){
                            f = sseq[s].coef[sseq[s].ord];
                            if (lf == 0.0 || lf * f < 0)
                                    atposinf++;
                    lf = f;
            }

            /*
             * changes at negative infinity
             */
            if ( (sseq[0].ord & 1) != 0 )
                            lf = -sseq[0].coef[sseq[0].ord];
            else
                            lf = sseq[0].coef[sseq[0].ord];

            for (s = 1; s <= np; s++ ) {
                            if ( ( sseq[s].ord & 1 ) != 0 )
                                    f = -sseq[s].coef[sseq[s].ord];
                            else
                                    f = sseq[s].coef[sseq[s].ord];
                            if (lf == 0 || lf * f < 0)
                                    atneginf++;
                            lf = f;
            }

            at[0] = atneginf;
            at[1] = atposinf;
            //	printf("atneginf, atposinf = %d %d\n", atneginf, atposinf);
            return(atneginf - atposinf);
    }