C++ Poly::degree方法代码示例

Poly divide(Poly const &a, Poly const &b, Poly &r) {
	Poly c;
	r = a; // remainder
	assert(!b.empty());

	const unsigned k = a.degree();
	const unsigned l = b.degree();
	c.resize(k, 0.);

	for(unsigned i = k; i >= l; i--) {
		assert(i >= 0);
		double ci = r.back()/b.back();
		c[i-l] += ci;
		Poly bb = ci*b;
		//std::cout << ci <<"*(" << b.shifted(i-l) << ") = "
		//          << bb.shifted(i-l) << "     r= " << r << std::endl;
		r -= bb.shifted(i-l);
		r.pop_back();
	}
	//std::cout << "r= " << r << std::endl;
	r.normalize();
	c.normalize();

	return c;
}

double laguerre_internal(Poly const & p,
                         Poly const & pp,
                         Poly const & ppp,
                         double x0,
                         double tol,
                         bool & quad_root) {
    double a = 2*tol;
    double xk = x0;
    double n = p.degree();
    quad_root = false;
    while(a*a > (tol*tol)) {
        //std::cout << "xk = " << xk << std::endl;
        double px = p(xk);
        if(px*px < tol*tol)
            return xk;
        double G = pp(xk) / px;
        double H = G*G - ppp(xk) / px;

        //std::cout << "G = " << G << "H = " << H;
        double radicand = (n - 1)*(n*H-G*G);
        assert(radicand > 0);
        //std::cout << "radicand = " << radicand << std::endl;
        if(G < 0) // here we try to maximise the denominator avoiding cancellation
            a = - sqrt(radicand);
        else
            a = sqrt(radicand);
        //std::cout << "a = " << a << std::endl;
        a = n / (a + G);
        //std::cout << "a = " << a << std::endl;
        xk -= a;
    }
    //std::cout << "xk = " << xk << std::endl;
    return xk;
}

cdouble laguerre_internal_complex(Poly const & p,
                                  double x0,
                                  double tol,
                                  bool & quad_root) {
    cdouble a = 2*tol;
    cdouble xk = x0;
    double n = p.degree();
    quad_root = false;
    const unsigned shuffle_rate = 10;
//    static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};
    unsigned shuffle_counter = 0;
    while(std::norm(a) > (tol*tol)) {
        //std::cout << "xk = " << xk << std::endl;
        cdouble b = p.back();
        cdouble d = 0, f = 0;
        double err = abs(b);
        double abx = abs(xk);
        for(int j = p.size()-2; j >= 0; j--) {
            f = xk*f + d;
            d = xk*d + b;
            b = xk*b + p[j];
            err = abs(b) + abx*err;
        }

        err *= 1e-7; // magic epsilon for convergence, should be computed from tol

        cdouble px = b;
        if(abs(b) < err)
            return xk;
        //if(std::norm(px) < tol*tol)
        //    return xk;
        cdouble G = d / px;
        cdouble H = G*G - f / px;

        //std::cout << "G = " << G << "H = " << H;
        cdouble radicand = (n - 1)*(n*H-G*G);
        //assert(radicand.real() > 0);
        if(radicand.real() < 0)
            quad_root = true;
        //std::cout << "radicand = " << radicand << std::endl;
        if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
            a = - sqrt(radicand);
        else
            a = sqrt(radicand);
        //std::cout << "a = " << a << std::endl;
        a = n / (a + G);
        //std::cout << "a = " << a << std::endl;
        if(shuffle_counter % shuffle_rate == 0)
		{
			//a *= shuffle[shuffle_counter / shuffle_rate];
		}
        xk -= a;
        shuffle_counter++;
        if(shuffle_counter >= 90)
            break;
    }
    //std::cout << "xk = " << xk << std::endl;
    return xk;
}

Poly Poly::operator*(const Poly& p) const {
    Poly result; 
    result.resize(degree() +  p.degree()+1);
    
    for(unsigned i = 0; i < size(); i++) {
        for(unsigned j = 0; j < p.size(); j++) {
            result[i+j] += (*this)[i] * p[j];
        }
    }
    return result;
}

void L::initNiederreiter
(GeneratorMatrixGen<typename A::type> &gm, A a, int d,
 const typename PolynomialRing<A>::type &irred)
{
    typedef typename A::type T;
    typedef PolynomialRing<A> PolyRing;
    typedef typename PolyRing::type Poly;

    const int degree = irred.degree ();

    const int vSize
        = std::max (gm.getM() + degree - 1,   // these elements are copied to gm
                    (gm.getPrec() + degree + 1));  // used in the loop

    Array<T> v (vSize);
    PolyRing poly (a);

    // cout << "Using " << irred << " for d=" << d << endl;

    Poly newPoly = poly.one();

    int u = 0;

    for (int j = 0; j < gm.getPrec(); j++)
    {
        // cout << "  j=" << j << endl;
        // Do we need a new v?

        if (u == 0)
        {
            Poly oldPoly = newPoly;
            int oldDegree = oldPoly.degree ();

            // calculate polyK+1 from polyK

            poly.mulBy (newPoly, irred);
            int newDegree = newPoly.degree ();
            // cout << "    newPolynomial: " << newPoly << endl;

            // kj can be set to any value between 0 <= kj < newDegree

            const int kj = oldDegree   // proposed by BFN
                           // newDegree - 1    // standard, bad???
                           // 0
                           // (newDegree > 3) ? 3 : oldDegree
                           ;

            std::fill (&v[0], &v[kj], T());  // Set leading v's to 0

            v[kj] = a.one();                     // the next one is 1

            if (kj < oldDegree)
            {
                T term = oldPoly [kj];

                for (int r = kj + 1; r < oldDegree; ++r)
                {
                    v [r] = a.one (); // 1 is arbitrary. Could be 0, too

                    a.addTo (term, a.mul (oldPoly [r], v [r]));
                }

                // set v[] not equal to -term

                v [oldDegree] = a.sub (a.one(), term);

                for (int r = oldDegree + 1; r < newDegree; ++r) v [r] = a.one(); //or 0
            }
            else
            {
                for (int r = kj + 1; r < newDegree; ++r) v [r] = a.one(); // or 0..
            }

            // All other elements are calculated by a recursion parameterized
            // by polyK

            for (int r = 0; r < vSize - newDegree; ++r)
            {
                T term = T();

                for (int i = 0; i < newDegree; ++i)
                {
                    a.addTo (term, a.mul (newPoly [i], v [r+i]));
                }
                v [newDegree + r] = a.neg (term);
            }
        }

        // Set data in ci

        for (int r = 0; r < gm.getM(); ++r)  gm.setd (d,r,j, v[r+u]);

        if (++u == degree) u = 0;
    }
}

void printPoly(const Poly& p)
{
    for (unsigned int i = 0; i <= p.degree(); ++i)
        cout << p.coef(i) << " ";
    cout << endl;
}

void bench_multiply() {
    std::default_random_engine generator;
    generator.seed(getNanoseconds());

    std::uniform_int_distribution<int> degreeDistrib(0, 255);

    std::vector<Poly> polys;

    for (int i = 0;  i < 1000; i++) {
        polys.push_back(Poly::random(degreeDistrib(generator), generator));
    }

    // 1 - Check Correctness of karatsubas
    {
        int tries = 0;
        int successes = 0;

        for (Poly p : polys) {
            for (Poly q : polys) {
                if (p.degree() + q.degree() >= 256) {
                    continue;
                }
                tries ++;

                Poly res1 = p.multiplyNaively(q);
                Poly res2 = p.multiplyKaratsuba32(q);

                if ((res1 + res2).size() == 0) {
                    successes ++;
                }
            }
        }

        std::cout << "Karatsuba32 success ratio : (" << successes << "/" << tries << ")" << std::endl;
    }

    // 2 - Bench the naive method
    {
        int forceBench = 0;

        auto start = std::chrono::high_resolution_clock::now();
        for (Poly p : polys) {
            for (Poly q : polys) {
                if (p.degree() + q.degree() >= 256) {
                    continue;
                }
                Poly r = p.multiplyNaively(q);
                forceBench += r.degree();
            }
        }
        auto end = std::chrono::high_resolution_clock::now();

        volatile int forceBench2 = forceBench;
        (void) forceBench2;

        std::cout << "Naive took " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << " ms" << std::endl;
    }

    // 3 - Bench the karatsuba32 method
    {
        int forceBench = 0;

        auto start = std::chrono::high_resolution_clock::now();
        for (Poly p : polys) {
            for (Poly q : polys) {
                if (p.degree() + q.degree() >= 256) {
                    continue;
                }
                Poly r = p.multiplyKaratsuba32(q);
                forceBench += r.degree();
            }
        }
        auto end = std::chrono::high_resolution_clock::now();

        volatile int forceBench2 = forceBench;
        (void) forceBench2;

        std::cout << "Karatsuba32 took " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << " ms" << std::endl;
    }
}

void bench_shifts() {
    std::default_random_engine generator;
    generator.seed(getNanoseconds());

    std::uniform_int_distribution<int> degreeDistrib(0, 255);

    std::vector<Poly> polys;

    for (int i = 0;  i < 10000; i++) {
        polys.push_back(Poly::random(degreeDistrib(generator), generator));
    }

    // 1 - Check Correctness of shifts
    {
        int tries = 0;
        int successes = 0;

        for (Poly p : polys) {
            for (unsigned i = 0; i < 256 - p.size(); i++) {
                tries ++;

                Poly res1 = naiveShiftLeft(p, i);
                Poly res2 = p << i;

                if ((res1 + res2).size() == 0) {
                    successes ++;
                }
            }

            for (unsigned i = 0;  i <= p.size(); i++) {
                tries ++;

                Poly res1 = naiveShiftRight(p, i);
                Poly res2 = p >> i;

                if ((res1 + res2).size() == 0) {
                    successes ++;
                }
            }
        }

        std::cout << "Shifts success ratio : (" << successes << "/" << tries << ")" << std::endl;
    }

    // 2 - Bench the naive method
    {
        int forceBench = 0;

        auto start = std::chrono::high_resolution_clock::now();
        for (Poly p : polys) {
            for (unsigned i = 0; i < 256 - p.size(); i++) {
                Poly r = naiveShiftLeft(p, i);
                forceBench += r.degree();
            }

            for (unsigned i = 0;  i <= p.size(); i++) {
                Poly r = naiveShiftRight(p, i);
                forceBench += r.degree();
            }
        }
        auto end = std::chrono::high_resolution_clock::now();

        volatile int forceBench2 = forceBench;
        (void) forceBench2;

        std::cout << "Naive took " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << " ms" << std::endl;
    }

    // 3 - Bench the fast method
    {
        int forceBench = 0;

        auto start = std::chrono::high_resolution_clock::now();
        for (Poly p : polys) {
            for (unsigned i = 0; i < 256 - p.size(); i++) {
                Poly r = p << i;
                forceBench += r.degree();
            }

            for (unsigned i = 0;  i <= p.size(); i++) {
                Poly r = p >> i;
                forceBench += r.degree();
            }
        }
        auto end = std::chrono::high_resolution_clock::now();

        volatile int forceBench2 = forceBench;
        (void) forceBench2;

        std::cout << "Naive took " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << " ms" << std::endl;
    }
}