Pagina's

2014/07/28

Sum of all divisors of all positive integers <= n.

/*
Example:  N  D        SD  SSD  A024916
          1  1         1    1
          2  1,2       3    4
          3  1,3       4    8
          4  1,2,4     7   15
          5  1,5       6   21
          6  1,2,3,6  12   33
  

Sum of all divisors of all positive integers <= 10^n.

                      A072692 as a simple table
 
     n                       a(n)                                  time
                                                          SSD               SSD3
        
     0                                           1  00:00:00.0000061
     1                                          87  00:00:00.0000058
     2                                        8299  00:00:00.0000092
     3                                      823081  00:00:00.0000215
     4                                    82256014  00:00:00.0000555
     5                                  8224740835  00:00:00.0002824
     6                                822468118437  00:00:00.0007623
     7                              82246711794796  00:00:00.0022483
     8                            8224670422194237  00:00:00.0075755
     9                          822467034112360628  00:00:00.0247891
    10                        82246703352400266400  00:00:00.1029980  00:00:00.0536221
    11                      8224670334323560419029  00:00:00.3310292
    12                    822467033425357340138978  00:00:01.0792704  00:00:00.5139608
    13                  82246703342420509396897774  00:00:03.3665267
    14                8224670334241228180927002517  00:00:10.7148327  00:00:05.1645348  
    15              822467033424114009326065894639  00:00:34.2567848
    16            82246703342411333689227187822414  00:01:49.7741595  00:00:53.1228885
    17          8224670334241132270081671519064067  00:05:55.3252935
    18        822467033424113219363487627735401433  00:19:50.2083007
    19      82246703342411321831834750379375676781  01:12:36.6691244
    20    8224670334241132182459113152714495956789  04:38:51.4847697
    21  822467033424113218236966661186847524013409  15:19:44.3197506


 AMD Athlon II X4 640, 3 GHz, 2 GB, XP  
  
    N=10^16     time(s)   CPU Usage(%)
      SSD(N)     110        25 
     SSD3(N)      53        60
     FSSD(N)      92        25
    FSSD3(N)      41        60
    FSSD4(N)      33        80
*/

using System;
using System.Diagnostics;
using System.Threading.Tasks;
using Xint = System.Numerics.BigInteger;
class A024916
{
    static Xint SSD(Xint N)
    {
        Xint D = 1, Q = N, S = 0;
        while (D < Q)
        {
            S += D++ * (Q - D);
            S += Q++ * Q / 2;
            Q = N / D;
        }
        S += Q++ * Q / 2;
        S -= D-- * D-- * D / 6;
        return S;
    }

    static void Main()
    {
        var sw = new Stopwatch();
        int p = 0;
        Xint N = 1, S;
    L0: sw.Restart();
        S = SSD(N);
        sw.Stop();
        Console.WriteLine(p + "  " + S + "  " + sw.Elapsed);
        p++;
        N *= 10;
        goto L0;
    }

    private static Xint[] SA, DA, QA;

    static Xint SSD3(Xint N)
    {
        SA = new Xint[3];
        DA = new Xint[3];
        QA = new Xint[3];
        Parallel.For(0, 3, (int d) => SSDP3(N, d));
        Xint S = SA[0] + SA[1] + SA[2],
             D = DA[0],
             Q = QA[0];
        if (DA[1] < D)
        {
            D = DA[1];
            Q = QA[1];
        }
        if (DA[2] < D)
        {
            D = DA[2];
            Q = QA[2];
        }
        S += Q++ * Q / 2;
        S -= D-- * D-- * D / 6;
        return S;
    }

    static void SSDP3(Xint N, int d)
    {
        Xint D = d + 1, Q = N / D, S = 0;
        while (D < Q)
        {
            S += D * (Q - D - 1);
            S += Q * (Q + 1) / 2;
            D += 3;
            Q = N / D;
        }
        SA[d] = S;
        DA[d] = D;
        QA[d] = Q;
    }

    static Xint FSSD(Xint N)
    {
        Xint D = 1, Q = N, S = 0;
        while (D < Q)
        {
            S += Q * (Q + 1 + 2 * D) / 2;
            D += 1;
            Q = N / D;
        }
        S += Q * (Q + 1) / 2;
        S -= D * D * (D - 1) / 2;
        return S;
    }

    static Xint FSSD3(Xint N)
    {
        SA = new Xint[3];
        DA = new Xint[3];
        QA = new Xint[3];
        Parallel.For(0, 3, (int d) => FSSDP3(N, d));
        Xint S = SA[0] + SA[1] + SA[2],
             D = DA[0],
             Q = QA[0];
        if (DA[1] < D)
        {
            D = DA[1];
            Q = QA[1];
        }
        if (DA[2] < D)
        {
            D = DA[2];
            Q = QA[2];
        }
        S += Q * (Q + 1) / 2;
        S -= D * D * (D - 1) / 2;
        return S;
    }

    static void FSSDP3(Xint N, int d)
    {
        Xint D = d + 1, Q = N / D, S = 0;
        while (D < Q)
        {
            S += Q * (Q + 1 + 2 * D) / 2;
            D += 3;
            Q = N / D;
        }
        SA[d] = S;
        DA[d] = D;
        QA[d] = Q;
    }

    static Xint FSSD4(Xint N)
    {
        SA = new Xint[4];
        DA = new Xint[4];
        QA = new Xint[4];
        Parallel.For(0, 4, (int d) => FSSDP4(N, d));
        Xint S = SA[0] + SA[1] + SA[2] + SA[3],
             D = DA[0],
             Q = QA[0];
        if (DA[1] < D)
        {
            D = DA[1];
            Q = QA[1];
        }
        if (DA[2] < D)
        {
            D = DA[2];
            Q = QA[2];
        }
        if (DA[3] < D)
        {
            D = DA[3];
            Q = QA[3];
        }
        S += Q * (Q + 1) / 2;
        S -= D * D * (D - 1) / 2;
        return S;
    }

    static void FSSDP4(Xint N, int d)
    {
        Xint D = d + 1, Q = N / D, S = 0;
        while (D < Q)
        {
            S += Q * (Q + 1 + 2 * D) / 2;
            D += 4;
            Q = N / D;
        }
        SA[d] = S;
        DA[d] = D;
        QA[d] = Q;
    }
}