aptestm.cpp

任意精度计算的实现

                    cerr << "Memory limiting calculation split at this level" << endl;
                    avproc = 1;
                }
            }

            if (avproc > 1)
            {
                // Processor available, calculate one half in another process in parallel
                startsubprocess (startproc, midproc - 1, N1, Nm, &proc1);

                // Set this process to use the other processors
                setprocessaffinity (midproc, endproc);
            }
            else
            {
                midproc = startproc;
                Nm = (N1 + N2) / 2;

                // No processor available, do the calculation in this process
                r (startproc, endproc, N1, Nm, &LT, &LQ, &LP);
            }

            // Do the calculation of the other half in this process
            r (midproc, endproc, Nm, N2, &RT, &RQ, (P ? &RP : 0));

            if (avproc > 1)
            {
                waitprocresult (&proc1, &LT, &LQ, &LP);

                currk += (size_t) ((Nm - N1) * (log ((double) (Nm - N1)) / log (2.0) + 1));

                print ();
            }
            else
            {
                currk += Nm - N1;
            }

            // Divide tasks by number of available processes
            // and available memory (all numbers must fit in memory for effective parallelization
            // The inverse square root takes about 1.5 - 1.7 times that of a multiplication

            avproc = endproc - startproc + 1;               // Number of processes we can use

            if (avproc > 1)
            {
                avthreads = avproc;                         // Number of threads we can use
                startthread = startproc;                    // Starting thread number

                avmem = availablememory (avproc);           // Memory available for us

                if (MAX_FORKED_PROCESSES < avproc)
                    avproc = MAX_FORKED_PROCESSES;

                if (toplevel && avmem / (2 * prec / Basedigits + 1) > 1)
                {
                    cerr << "Calculating isqrt in parallel" << endl;
                    avproc--;
                    avmem -= 2 * prec / Basedigits + 1;
                    sqrtp = true;
                }

                maxmem = avmem / (RQ.ap->size + LT.ap->size + 2);
                if (maxmem < avproc)
                {
                    cerr << "Memory limiting processes from " << avproc << " to " << maxmem << endl;
                    avproc = maxmem;
                }

                assert (avproc);                            // Must have sufficient memory

                if (sqrtp)
                {
                    avproc++;                               // Will be decreased by getnthreads
                    getnthreads (&avthreads, &avproc, &startthread, &proc3);
                    startsqrt (&proc3);
                }

                // Set number threads to be used in this process
                getnthreads (&avthreads, &avproc, &startthread, &proc1);
                setprocessaffinity (proc1.startproc, proc1.endproc);
            }
            else
            {
                avproc = 0;
            }

            switch (avproc)
            {
                case 0:
                {
                    *T = RQ * LT + LP * RT;
                    *Q = LQ * RQ;
                    if (P) *P = LP * RP;

                    break;
                }
                case 1:
                {
                    getnthreads (&avthreads, &avproc, &startthread, &proc1);

                    tmp = LP;
                    startmultiply (&tmp, &RT, &proc1);
                    *T = RQ * LT;
                    waitprocresult (&proc1, &tmp);
                    *T += tmp;

                    if (P)
                        startmultiply (&LP, &RP, &proc1);

                    *Q = LQ * RQ;

                    if (P)
                        waitprocresult (&proc1, P);

                    break;
                }
                case 2:
                {
                    getnthreads (&avthreads, &avproc, &startthread, &proc1);
                    getnthreads (&avthreads, &avproc, &startthread, &proc2);

                    tmp = LP;
                    startmultiply (&tmp, &RT, &proc1);
                    tmp = RQ;
                    startmultiply (&tmp, &LT, &proc2);

                    *Q = LQ * RQ;
                    if (P) *P = LP * RP;

                    waitprocresult (&proc1, T);
                    waitprocresult (&proc2, &tmp);

                    *T += tmp;


                    break;
                }
                default:
                {
                    getnthreads (&avthreads, &avproc, &startthread, &proc1);
                    getnthreads (&avthreads, &avproc, &startthread, &proc2);
                    if (P)
                        getnthreads (&avthreads, &avproc, &startthread, &proc3);

                    tmp = LP;
                    startmultiply (&tmp, &RT, &proc1);
                    tmp = RQ;
                    startmultiply (&tmp, &LT, &proc2);

                    if (P)
                        startmultiply (&LP, &RP, &proc3);

                    *Q = LQ * RQ;

                    waitprocresult (&proc1, T);
                    waitprocresult (&proc2, &tmp);
                    *T += tmp;

                    if (P)
                        waitprocresult (&proc3, P);

                    break;
                }
            }

            if (sqrtp)
                waitprocresult (&proc3, &isqrt);            // Read the inverse square root of 640320

            if (endproc > startproc)
                setprocessaffinity (startproc, endproc);    // Restore threads

            currk += N2 - Nm;
        }
    }

    print ();
}

void calcpart (size_t startproc, size_t endproc, size_t startterm, size_t endterm, apfloat *T, apfloat *Q, apfloat *P, bool toplevel)
{
    size_t n;
    time_t tt;

    if (toplevel)
    {
        cerr << "Using the Chudnovsky brothers' binary splitting algorithm" << endl;
    }

    A = 13591409;
    B = 545140134;
    J = 101568000; J *= 107701824;       // J = 10939058860032000

    one = 1;
    two = 2;
    five = 5;
    six = 6;

    A.prec (prec);
    B.prec (prec);
    J.prec (prec);

    one.prec (prec);
    two.prec (prec);
    five.prec (prec);
    six.prec (prec);

    n = endterm - startterm + 1;
    maxk = (size_t) (n * (log ((double) n) / log (2.0) + 1));
    currk = 0;

    tt = time (0);

    r (startproc, endproc, startterm, endterm, T, Q, P, toplevel);

    if (toplevel)
    {
        cerr.setf (ios::fixed, ios::floatfield);
        cerr << "100% complete, elapsed time " << setprecision (0) << difftime (time (0), tt) << " seconds" << endl;
    }
}

apfloat final (apfloat T, apfloat Q)
{
    time_t tt;
    apfloat p;

    cerr << "Final value";

    tt = time (0);

    if (!isqrt.ap) isqrt = invroot (apfloat (640320, prec), 2);
    p = invroot (isqrt * T, 1) * 53360 * Q;

    cerr.setf (ios::fixed, ios::floatfield);
    cerr << " took " << setprecision (0) << difftime (time (0), tt) << " seconds" << endl;

    return p;
}

int main (int argc, char *argv[], char *env[])
{
    size_t startproc, endproc, startterm, endterm;
    int n = 2, m = 0, b = 3, base = 10;
    time_t tt;
    apfloat T, Q, P;

    if (argc > 1)
    {
        if (!strncmp (argv[1], "-p", 2))
        {
            m = 1;
            n = 7;
            b = 7;
        }
        else if (!strncmp (argv[1], "-m", 2))
        {
            m = 2;
            n = 6;
            b = 6;
        }
        else if (!strncmp (argv[1], "-s", 2))
        {
            m = 3;
            n = 5;
            b = 5;
        }
    }

    if (argc < n)
    {
        cerr << "USAGE: aptestm digits [processes] [base]" << endl;
        cerr << "       aptestm -p[rocess] digits startproc endproc startterm endterm [base]" << endl;
        cerr << "       aptestm -m[ultiply] startproc endproc in1 in2 [base]" << endl;
        cerr << "       aptestm -s[qrt] startproc endproc digits [base]" << endl;
        cerr << "    base must be 2...36" << endl;
        return 2;
    }

    if (argc > b)
    {
        istringstream s (argv[b]);
        if (!(s >> base))
        {
            cerr << "Invalid argument base: " << argv[b] << endl;
            return 1;
        }
    }

    progname = argv[0];
    baseenv = env;

    apbase (base);

    if (m == 0)
    {
        apfloat p;

        istringstream s (argv[1]);
        if (!(s >> prec) || !prec)
        {
            cerr << "Invalid argument digits: " << argv[1] << endl;
            return 1;
        }

        if (argc > 2)
        {
            istringstream s (argv[2]);
            if (!(s >> endproc) || !endproc || endproc > 32)
            {
                cerr << "Invalid argument processes: " << argv[2] << endl;
                return 1;
            }
        }
        else
        {
            endproc = NProcessors;          // Assumed to be set or detected appropriately
            if (!endproc) endproc = 1;
            else if (endproc > 32) endproc = 32;
        }

        dump ();

        cerr << "Calculating pi to " << prec << " base-" << Basedigit << " digits";
        if (endproc > 1) cerr << " using up to " << endproc << " processes";
        cerr << endl;

        tt = time (0);

        startterm = 0;
        endterm = (size_t) (prec * log ((double) Basedigit) / 32.65445004177);
        totalproc = endproc;
        startproc = 0;
        endproc--;              // Count starts from zero

        output = true;

        setfno (startproc, endproc);
        setprocessaffinity (startproc, endproc);

        calcpart (startproc, endproc, startterm, endterm + 1, &T, &Q, 0, true);
        p = final (T, Q);

        cout << pretty << p << endl;

        cerr.setf (ios::fixed, ios::floatfield);
        cerr << "Total elapsed time " << setprecision (0) << difftime (time (0), tt) << " seconds" << endl;
    }
    else if (m == 1)
    {
        istringstream s2 (argv[2]);
        if (!(s2 >> prec) || !prec)
        {
            cerr << "Invalid argument digits: " << argv[2] << endl;
            return 1;
        }

        istringstream s3 (argv[3]);
        if (!(s3 >> startproc))
        {
            cerr << "Invalid argument startproc: " << argv[3] << endl;
            return 1;
        }

        istringstream s4 (argv[4]);
        if (!(s4 >> endproc) || endproc < startproc)
        {
            cerr << "Invalid argument endproc: " << argv[4] << endl;
            return 1;
        }

        istringstream s5 (argv[5]);
        if (!(s5 >> startterm))
        {
            cerr << "Invalid argument startterm: " << argv[5] << endl;
            return 1;
        }

        istringstream s6 (argv[6]);
        if (!(s6 >> endterm) || endterm < startterm)
        {
            cerr << "Invalid argument endterm: " << argv[6] << endl;
            return 1;
        }

        output = false;

        totalproc = endproc - startproc + 1;

        setfno (startproc, endproc);
        setprocessaffinity (startproc, endproc);

        calcpart (startproc, endproc, startterm, endterm, &T, &Q, &P, false);

        writenumber (startproc, endproc, &T, &Q, &P);
    }
    else if (m == 2)
    {
        apfloat x, y;

        istringstream s2 (argv[2]);
        if (!(s2 >> startproc))
        {
            cerr << "Invalid argument startproc: " << argv[2] << endl;
            return 1;
        }

        istringstream s3 (argv[3]);
        if (!(s3 >> endproc) || endproc < startproc)
        {
            cerr << "Invalid argument endproc: " << argv[3] << endl;
            return 1;
        }

        totalproc = endproc - startproc + 1;

        setfno (startproc, endproc);
        setprocessaffinity (startproc, endproc);

        readnumber (startproc, endproc, &x, argv[4][0]);
        readnumber (startproc, endproc, &y, argv[5][0]);

        x *= y;

        writenumber (startproc, endproc, &x);
    }
    else if (m == 3)
    {
        apfloat x;

        istringstream s2 (argv[2]);
        if (!(s2 >> startproc))
        {
            cerr << "Invalid argument startproc: " << argv[2] << endl;
            return 1;
        }

        istringstream s3 (argv[3]);
        if (!(s3 >> endproc) || endproc < startproc)
        {
            cerr << "Invalid argument endproc: " << argv[3] << endl;
            return 1;
        }

        istringstream s4 (argv[4]);
        if (!(s4 >> prec) || !prec)
        {
            cerr << "Invalid argument digits: " << argv[4] << endl;
            return 1;
        }

        totalproc = endproc - startproc + 1;

        setfno (startproc, endproc);
        setprocessaffinity (startproc, endproc);

        x = invroot (apfloat (640320, prec), 2);

        writenumber (startproc, endproc, &x);
    }

    return 0;
}