Efficient python code for printing the product of divisors of a number

Question

I am trying to solve a problem involving printing the product of all divisors of a given number. The number of test cases is a number 1 <= t <= 300000 , and the number itself can range from 1 <= n <= 500000

I wrote the following code, but it always exceeds the time limit of 2 seconds. Are there any ways to speed up the code ?

from math import sqrt

def divisorsProduct(n):
    ProductOfDivisors=1

    for i in range(2,int(round(sqrt(n)))+1):
        if n%i==0: 
            ProductOfDivisors*=i

            if n/i != i:
                ProductOfDivisors*=(n/i)    

    if ProductOfDivisors <= 9999:
        print ProductOfDivisors
    else:
        result = str(ProductOfDivisors)
        print result[len(result)-4:] 

T = int(raw_input())

for i in range(1,T+1):
    num = int(raw_input())
    divisorsProduct(num)

Thank You.

Answer 1

You need to clarify by what you mean by "product of divisors." The code posted in the question doesn't work for any definition yet. This sounds like a homework question. If it is, then perhaps your instructor was expecting you to think outside the code to meet the time goals.

If you mean the product of unique prime divisors, e.g., 72 gives 2*3 = 6, then having a list of primes is the way to go. Just run through the list up to the square root of the number, multiplying present primes into the result. There are not that many, so you could even hard code them into your program.

If you mean the product of all the divisors, prime or not, then it is helpful to think of what the divisors are. You can make serious speed gains over the brute force method suggested in the other answers and yours. I suspect this is what your instructor intended.

If the divisors are ordered in a list, then they occur in pairs that multiply to n -- 1 and n, 2 and n/2, etc. -- except for the case where n is a perfect square, where the square root is a divisor that is not paired with any other.

So the result will be n to the power of half the number of divisors, (regardless of whether or not n is a square).

To compute this, find the prime factorization using your list of primes. That is, find the power of 2 that divides n, then the power of 3, etc. To do this, take out all the 2s, then the 3s, etc.

The number you are taking the factors out of will be getting smaller, so you can do the square root test on the smaller intermediate numbers to see if you need to continue up the list of primes. To gain some speed, test p*p <= m, rather than p <= sqrt(m)

Once you have the prime factorization, it is easy to find the number of divisors. For example, suppose the factorization is 2^i * 3^j * 7^k. Then, since each divisor uses the same prime factors, with exponents less than or equal to those in n including the possibility of 0, the number of divisors is (i+1)(j+1)(k+1).

E.g., 72 = 2^3 * 3^2, so the number of divisors is 4*3 = 12, and their product is 72^6 = 139,314,069,504.

By using math, the algorithm can become much better than O(n). But it is hard to estimate your speed gains ahead of time because of the relatively small size of the n in the input.

Answer 2

I tried this problem but still i don't understand it this is my program in C++ `

#include <iostream>
#include <cmath>
using namespace std;
#include <stdio.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <conio.h>

void long_pow(long long n, int div);

int primes[] = {0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701};
int prod[500001] ={0};
bool done[500001] = {false};

void calc_product(int n)
{
     if(prod[n]==0)
     {
         int index = 1;
         int n1 = n;
         bool flag = false;
         int total = 1;
         int cur = 1;
         while(n!=0 && index<=126 &&(!flag) )
         {
                    //checking if it is prime
                    if(n1 == primes[index]) break;
                    if(n%primes[index]==0)
                    {
                                          cur++;
                                         // total = total*primes[index];
                                          n/=primes[index];
                    }
                    else
                    {
                        if(n==1) flag = true;
                      //  cout<<primes[index]<<" "<<cur<<" "<<n<<endl;
                        total = total*cur;
                        index++;
                        cur = 1;
                    }
         }
         if(n>1) total*=2;
         //cout<<total<<endl;
        // long_pow(n1,total);
        prod[n1] = (int)pow((double)n1,(total-2)/2);
         //now printing
         if(done[n1])
             {
                 printf("%d%d%d%d",prod[n1]/1000,(prod[n1]%1000)/100,(prod[n1]%100)/10,prod[n1]%10);
             }
             else
                 printf("%d",prod[n1]);
     }
     else
     {
         if(done[n])
         {
             printf("%d%d%d%d",prod[n]/1000,(prod[n]%1000)/100,(prod[n]%100)/10,prod[n]%10);
         }
         else
             printf("%d",prod[n]);
     }
     //cout<<<<endl;
}

void long_pow(long long n, int div)
{
     long long n1 = n;
     bool flag = false;
     long long total = 1;
     //not including the number itself
     if(div%2 == 0)
     {
              div = div/2;
              for(int i=1;i<div;i++)
              {
                      total = total*n;
                      if(total>9999)
                               {
                                    flag = true;
                                    total%=10000;
                               }
                      if(flag)
                           total%=10000;
              }
     }
     else
     {
         n = (int)sqrt((double)n);
         for(int i=1;i<div-1;i++)
         {
                 total*=n;
                 if(total>9999)
                               {
                                       flag = true;
                                       total%=10000;
                               }
                 if(flag)
                           total%=10000;
         }
     }
     done[n1] = flag;
     prod[n1] = total;
     /*if(flag)
     {
             cout<<total/1000<<(total%1000)/100<<(total%100)/10<<total%10;
     }
     else
         cout<<total;
     *///return total;
}

int main()
{
    int cases;
    scanf("%d",&cases);
    while(cases)
    {
    int i;
    scanf("%d",&i);
  // for(i=2;i<=500000;i++)
    calc_product(i);
    cases--;
    if(cases) printf("\n");
    }
   //total time
   printf("time=%.3lf sec.\n",
    (double) (clock())/CLOCKS_PER_SEC);
    getch();
    return 0;
}

now the problem is this is still not sufficient to solve the problem.. i have precomputed my list of primes and and decomposing the number into the prime factors still..got problem .. plz tell me where do i go to ..or help me make my method faster ...

Answer 3

Okay, I think this is close to the optimal algorithm. It produces the product_of_divisors for each number in range(500000).

import math

def number_of_divisors(maxval=500001):
    """ Example: the number of divisors of 12 is 6:  1, 2, 3, 4, 6, 12.
        Given a prime factoring of n, the number of divisors of n is the
        product of each factor's multiplicity plus one (mpo in my variables).

        This function works like the Sieve of Eratosthenes, but marks each
        composite n with the multiplicity (plus one) of each prime factor. """
    numdivs = [1] * maxval   # multiplicative identity
    currmpo = [0] * maxval

    # standard logic for 2 < p < sqrt(maxval)
    for p in range(2, int(math.sqrt(maxval))):
        if numdivs[p] == 1:   # if p is prime
            for exp in range(2,50): # assume maxval < 2^50
                pexp = p ** exp
                if pexp > maxval:
                    break
                exppo = exp + 1
                for comp in range(pexp, maxval, pexp):
                    currmpo[comp] = exppo
            for comp in range(p, maxval, p):
                thismpo = currmpo[comp] or 2
                numdivs[comp] *= thismpo
                currmpo[comp] = 0  # reset currmpo array in place

    # abbreviated logic for p > sqrt(maxval)
    for p in range(int(math.sqrt(maxval)), maxval):
        if numdivs[p] == 1:   # if p is prime
            for comp in range(p, maxval, p):
                numdivs[comp] *= 2

    return numdivs

# this initialization times at 7s on my machine
NUMDIV = number_of_divisors()

def product_of_divisors(n):
    if NUMDIV[n] % 2 == 0:
        # each pair of divisors has product equal to n, for example
        # 1*12 * 2*6 * 3*4  =  12**3
        return n ** (NUMDIV[n] / 2)
    else:
        # perfect squares have their square root as an unmatched divisor
        return n ** (NUMDIV[n] / 2) * int(math.sqrt(n))

# this loop times at 13s on my machine
for n in range(500000):
    a = product_of_divisors(n)

On my very slow machine, it takes 7s to compute the numberofdivisors for each number, then 13s to compute the productofdivisors for each. Of course it can be sped up by translating it into C. (@someone with a fast machine: how long does it take on your machine?)