Monday, 25 November 2013

BigInteger Program in C++

Once I was searching for BigInteger programs in C++, then I found this one. I really don't remember the source, as I copied it then.
But now as I am using this program for many problems, I thought of adding this one as it is really good.

If you use Java or python, you will never need this program.
















// The BigInteger Ripper

#include <vector>
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <string>
using namespace std;
typedef long long LL;
// base and base_digits must be consistent
const int base = 1000000000;
const int base_digits = 9;

struct bigint {
    vector<int> a;
    int sign;

    bigint() :
        sign(1) {
    }

    bigint(long long v) {
        *this = v;
    }

    bigint(const string &s) {
        read(s);
    }

    void operator=(const bigint &v) {
        sign = v.sign;
        a = v.a;
    }

    void operator=(long long v) {
        sign = 1;
        if (v < 0)
            sign = -1, v = -v;
        for (; v > 0; v = v / base)
            a.push_back(v % base);
    }

    bigint operator+(const bigint &v) const {
        if (sign == v.sign) {
            bigint res = v;

            for (int i = 0, carry = 0; i < (int) max(a.size(), v.a.size()) || carry; ++i) {
                if (i == (int) res.a.size())
                    res.a.push_back(0);
                res.a[i] += carry + (i < (int) a.size() ? a[i] : 0);
                carry = res.a[i] >= base;
                if (carry)
                    res.a[i] -= base;
            }
            return res;
        }
        return *this - (-v);
    }

    bigint operator-(const bigint &v) const {
        if (sign == v.sign) {
            if (abs() >= v.abs()) {
                bigint res = *this;
                for (int i = 0, carry = 0; i < (int) v.a.size() || carry; ++i) {
                    res.a[i] -= carry + (i < (int) v.a.size() ? v.a[i] : 0);
                    carry = res.a[i] < 0;
                    if (carry)
                        res.a[i] += base;
                }
                res.trim();
                return res;
            }
            return -(v - *this);
        }
        return *this + (-v);
    }

    void operator*=(int v) {
        if (v < 0)
            sign = -sign, v = -v;
        for (int i = 0, carry = 0; i < (int) a.size() || carry; ++i) {
            if (i == (int) a.size())
                a.push_back(0);
            long long cur = a[i] * (long long) v + carry;
            carry = (int) (cur / base);
            a[i] = (int) (cur % base);
            //asm("divl %%ecx" : "=a"(carry), "=d"(a[i]) : "A"(cur), "c"(base));
        }
        trim();
    }

    bigint operator*(int v) const {
        bigint res = *this;
        res *= v;
        return res;
    }

    friend pair<bigint, bigint> divmod(const bigint &a1, const bigint &b1) {
        int norm = base / (b1.a.back() + 1);
        bigint a = a1.abs() * norm;
        bigint b = b1.abs() * norm;
        bigint q, r;
        q.a.resize(a.a.size());

        for (int i = a.a.size() - 1; i >= 0; i--) {
            r *= base;
            r += a.a[i];
            int s1 = r.a.size() <= b.a.size() ? 0 : r.a[b.a.size()];
            int s2 = r.a.size() <= b.a.size() - 1 ? 0 : r.a[b.a.size() - 1];
            int d = ((long long) base * s1 + s2) / b.a.back();
            r -= b * d;
            while (r < 0)
                r += b, --d;
            q.a[i] = d;
        }

        q.sign = a1.sign * b1.sign;
        r.sign = a1.sign;
        q.trim();
        r.trim();
        return make_pair(q, r / norm);
    }

    bigint operator/(const bigint &v) const {
        return divmod(*this, v).first;
    }

    bigint operator%(const bigint &v) const {
        return divmod(*this, v).second;
    }

    void operator/=(int v) {
        if (v < 0)
            sign = -sign, v = -v;
        for (int i = (int) a.size() - 1, rem = 0; i >= 0; --i) {
            long long cur = a[i] + rem * (long long) base;
            a[i] = (int) (cur / v);
            rem = (int) (cur % v);
        }
        trim();
    }

    bigint operator/(int v) const {
        bigint res = *this;
        res /= v;
        return res;
    }

    int operator%(int v) const {
        if (v < 0)
            v = -v;
        int m = 0;
        for (int i = a.size() - 1; i >= 0; --i)
            m = (a[i] + m * (long long) base) % v;
        return m * sign;
    }

    void operator+=(const bigint &v) {
        *this = *this + v;
    }
    void operator-=(const bigint &v) {
        *this = *this - v;
    }
    void operator*=(const bigint &v) {
        *this = *this * v;
    }
    void operator/=(const bigint &v) {
        *this = *this / v;
    }

    bool operator<(const bigint &v) const {
        if (sign != v.sign)
            return sign < v.sign;
        if (a.size() != v.a.size())
            return a.size() * sign < v.a.size() * v.sign;
        for (int i = a.size() - 1; i >= 0; i--)
            if (a[i] != v.a[i])
                return a[i] * sign < v.a[i] * sign;
        return false;
    }

    bool operator>(const bigint &v) const {
        return v < *this;
    }
    bool operator<=(const bigint &v) const {
        return !(v < *this);
    }
    bool operator>=(const bigint &v) const {
        return !(*this < v);
    }
    bool operator==(const bigint &v) const {
        return !(*this < v) && !(v < *this);
    }
    bool operator!=(const bigint &v) const {
        return *this < v || v < *this;
    }

    void trim() {
        while (!a.empty() && !a.back())
            a.pop_back();
        if (a.empty())
            sign = 1;
    }

    bool isZero() const {
        return a.empty() || (a.size() == 1 && !a[0]);
    }

    bigint operator-() const {
        bigint res = *this;
        res.sign = -sign;
        return res;
    }

    bigint abs() const {
        bigint res = *this;
        res.sign *= res.sign;
        return res;
    }

    long long longValue() const {
        long long res = 0;
        for (int i = a.size() - 1; i >= 0; i--)
            res = res * base + a[i];
        return res * sign;
    }

    friend bigint gcd(const bigint &a, const bigint &b) {
        return b.isZero() ? a : gcd(b, a % b);
    }
    friend bigint lcm(const bigint &a, const bigint &b) {
        return a / gcd(a, b) * b;
    }

    void read(const string &s) {
        sign = 1;
        a.clear();
        int pos = 0;
        while (pos < (int) s.size() && (s[pos] == '-' || s[pos] == '+')) {
            if (s[pos] == '-')
                sign = -sign;
            ++pos;
        }
        for (int i = s.size() - 1; i >= pos; i -= base_digits) {
            int x = 0;
            for (int j = max(pos, i - base_digits + 1); j <= i; j++)
                x = x * 10 + s[j] - '0';
            a.push_back(x);
        }
        trim();
    }
   
    int length(){
    int l=0,back=a.back();
    while(back){l++;back/=10;}
    l+=((a.size()-1)*base_digits);
    return l;
    }

    friend istream& operator>>(istream &stream, bigint &v) {
        string s;
        stream >> s;
        v.read(s);
        return stream;
    }

    friend ostream& operator<<(ostream &stream, const bigint &v) {
        if (v.sign == -1)
            stream << '-';
        stream << (v.a.empty() ? 0 : v.a.back());
        for (int i = (int) v.a.size() - 2; i >= 0; --i)
            stream << setw(base_digits) << setfill('0') << v.a[i];
        return stream;
    }

    static vector<int> convert_base(const vector<int> &a, int old_digits, int new_digits) {
        vector<long long> p(max(old_digits, new_digits) + 1);
        p[0] = 1;
        for (int i = 1; i < (int) p.size(); i++)
            p[i] = p[i - 1] * 10;
        vector<int> res;
        long long cur = 0;
        int cur_digits = 0;
        for (int i = 0; i < (int) a.size(); i++) {
            cur += a[i] * p[cur_digits];
            cur_digits += old_digits;
            while (cur_digits >= new_digits) {
                res.push_back(int(cur % p[new_digits]));
                cur /= p[new_digits];
                cur_digits -= new_digits;
            }
        }
        res.push_back((int) cur);
        while (!res.empty() && !res.back())
            res.pop_back();
        return res;
    }

    typedef vector<long long> vll;

    static vll karatsubaMultiply(const vll &a, const vll &b) {
        int n = a.size();
        vll res(n + n);
        if (n <= 32) {
            for (int i = 0; i < n; i++)
                for (int j = 0; j < n; j++)
                    res[i + j] += a[i] * b[j];
            return res;
        }

        int k = n >> 1;
        vll a1(a.begin(), a.begin() + k);
        vll a2(a.begin() + k, a.end());
        vll b1(b.begin(), b.begin() + k);
        vll b2(b.begin() + k, b.end());

        vll a1b1 = karatsubaMultiply(a1, b1);
        vll a2b2 = karatsubaMultiply(a2, b2);

        for (int i = 0; i < k; i++)
            a2[i] += a1[i];
        for (int i = 0; i < k; i++)
            b2[i] += b1[i];

        vll r = karatsubaMultiply(a2, b2);
        for (int i = 0; i < (int) a1b1.size(); i++)
            r[i] -= a1b1[i];
        for (int i = 0; i < (int) a2b2.size(); i++)
            r[i] -= a2b2[i];

        for (int i = 0; i < (int) r.size(); i++)
            res[i + k] += r[i];
        for (int i = 0; i < (int) a1b1.size(); i++)
            res[i] += a1b1[i];
        for (int i = 0; i < (int) a2b2.size(); i++)
            res[i + n] += a2b2[i];
        return res;
    }

    bigint operator*(const bigint &v) const {
        vector<int> a6 = convert_base(this->a, base_digits, 6);
        vector<int> b6 = convert_base(v.a, base_digits, 6);
        vll a(a6.begin(), a6.end());
        vll b(b6.begin(), b6.end());
        while (a.size() < b.size())
            a.push_back(0);
        while (b.size() < a.size())
            b.push_back(0);
        while (a.size() & (a.size() - 1))
            a.push_back(0), b.push_back(0);
        vll c = karatsubaMultiply(a, b);
        bigint res;
        res.sign = sign * v.sign;
        for (int i = 0, carry = 0; i < (int) c.size(); i++) {
            long long cur = c[i] + carry;
            res.a.push_back((int) (cur % 1000000));
            carry = (int) (cur / 1000000);
        }
        res.a = convert_base(res.a, 6, base_digits);
        res.trim();
        return res;
    }
};


void recursive_fraction(bigint a,bigint b)
{
  bigint n,i;
  if(a>b)
  {
    n=a/b;
    for(i=1;i<=n;i+=1)
    {
       cout<<"1 ";
    }
    a=(a-b*n);
  }
  while(a!=0)
  {
    n=b/a;
    if(b%a!=0) n=n+1;
    cout<<n<<" ";
    a=n*a-b;
    b=n*b;
  }
}
int main()
{
 while(true)
 {
      bigint num;
      cin>>num;
      if(num == -1)break;
      bigint nine = 9;
      bigint ans = num/nine;
      bigint zero = 0;
      if(num%nine != zero)
             ans = ans + 1;
      cout<<ans<<endl;
 }
  return 0;
}


PS: (My implementation of Biginteger add, multiply and subtraction, just coded it once in the beginning)

// My Implementation

#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
typedef long long int int64;
using namespace std;
char op1[100001],op2[100001],res[200003];

char * convert(int64 a,char * num)
{
     int len = int(log10(double(a)));
     num[len+1] = '\0';
     while(a>0)
     {
      num[len--] = a%10 + '0';
      a = a/10;        
     }
     return num;
}

char * convert(char * a,char * num)
{
     strcpy(num,a);
     return num;
}

char * reverse(char * num)
{
 char x[200003];int lo = 0;
 int len = strlen(num);
 while(num[len-1]=='0' && len>=1)
                       len--;

 if(len==0)
           {num[0]='0';num[1]='\0';return num;}

 for(int xx = len-1;xx>=0;xx--)
 {
         x[lo++] = num[xx];
 }  
 x[lo]='\0';
 strcpy(num,x);
 return num;
}

char * add(char * a, char * b, char * num)
{
     int len1 = strlen(a)-1;
     int len2 = strlen(b)-1;
     int carry = 0;
     int x = 0;
     while(len1>=0 && len2>=0)
     {
      int temp = a[len1--]-'0'+ b[len2--]-'0' + carry;            
      num[x++] = temp%10 + '0';
      carry = temp/10;
     }   
     while(len1>=0)
     {
      int temp = a[len1--]-'0'+ carry;            
      num[x++] = temp%10 + '0';
      carry = temp/10;
     }
     while(len2>=0)
     {
      int temp = b[len2--]-'0'+ carry;            
      num[x++] = temp%10 + '0';
      carry = temp/10;
     }
     if(carry>0)
                num[x++] = carry+'0';
     num[x] = '\0';   
     return reverse(num);
}

char * sub(char * a, char * b, char * num)
{
     int len1 = strlen(a)-1;
     int len2 = strlen(b)-1;
     int carry = 0;
     int x = 0;int togl = 0;
     while(len2>=0)
     {
      int temp;
      if(a[len1]>=b[len2])
      {
       temp = a[len1] - b[len2];
       num[x++] = temp + '0';                   
       len1--,len2--;
      }
      else
      {
       int flag;togl = 1;
       while(a[len1]<=b[len2])
       {
        if(togl==1){num[x++] = a[len1] - b[len2]+10 +'0';togl=0;}
        else {num[x++] = a[len1] - b[len2]+9 +'0';}
        len1--,len2--;            
        if(len2<0)
                  {flag=1;break;} 
       }  
       while(a[len1]=='0'){num[x++]='9';len1--;}
       num[x++] =    a[len1] - (flag==1?'0':b[len2]) - 1 +'0';
       len2--,len1--;
      }            
     }
     while(len1>=0)
                   num[x++]=a[len1--];
     num[x]= '\0';

     return reverse(num);
}


char * mult(char * a, char * b, char * num)
{
     int len1 = strlen(a)-1;
     int len2 = strlen(b)-1;
     int x=0,carry=0,ll=0,mm=0,prevll=-1;
     for(int i=len1;i>=0;i--)
     {
      int mul1 = a[i]-'0';
      for(int j=len2;j>=0;j--)      
             {
                int mul2 = b[j]-'0';          
                int temp = mul1*mul2 + carry + (mm<=prevll?num[mm]-'0':0);
                num[mm++] = temp%10 +'0';
                carry = temp/10;                                    
             }
             if(carry>0)
                        {num[mm++]=carry+'0';
                         carry = 0;}
             prevll = mm-1;
             ll++;
             if(i==0){num[mm]='\0';break;}
             mm = ll;
     }
     return reverse(num);
}



int main()
{
 char a[100001], b[100001];
 while(cin>>a>>b)
 {
  //cout<<add(convert(a,op1),convert(b,op2),res)<<endl;
  cout<<sub(convert(a,op1),convert(b,op2),res)<<endl;
  //cout<<mult(convert(a,op1),convert(b,op2),res)<<endl; 
 
 }
return 0;
}

Friday, 22 March 2013

Fast Euler Totient Function in C++

Euler Totient Function or phi(n) gives the No. of numbers less than or equal to n which are relatively prime to n.Some Postulates of Euler totient Function-
phi(p) = p-1 ; if p is a prime ;
phi(p^n) = (p-1) * p^(n-1) ; still if p is a prime..
phi(a*b) = phi(a) * phi(b) whenever gcd(a,b) = 1
 Read ETF in detail from http://en.wikipedia.org/wiki/Euler%27s_totient_function

Now an easy but slow implementation will be to prime factorize N and then apply phi(p^n) formula to all the prime factors in prime factorization of N.
i.e
Let N = p1^e1 * p2^e2 *..... pn^en
Then
phi(N) = phi(p1^e1) * phi(p2^e2) * ... phi(pn^en)
phi(N) = ((p1 - 1) * (p1 ^ e1 - 1)) * ((p2 - 1) * (p2 ^ (e2 - 1))) * ((pn - 1) * (pn ^ (en - 1)))

Implementation of this one is -
 
#include<iostream>
#include<cmath>
using namespace std;

long long int phi(long long x)
 {
   long long int ret = 1,i,pow;
   for (i = 2; x != 1; i++) 
   {
     pow = 1;
     if(i>sqrt(x))break;
     while (!(x%i)) 
     {
       x /= i;
       pow *= i;
     }
     ret *= (pow - (pow/i));
    }
    if(x!=1)ret*=(x-1);
    return ret;
}

int main()
{
    int t;
    cin>>t;
    while(t--)
    {
           int n;
           cin>>n;
           cout<<phi(n)<<endl;   
    }
    return 0;
} 
 
Complexity of this Solution is O(T * sqrt(N)) as T can be maximum 10^3 and N <= 10^6
so overall complexity of this solution is O(10^6) which will not time out...
But there exists many better and faster implementation of the Euler Totient Function by using sieve

I will discuss two of them -
First implementation -
In first implementation using sieve we will keep one of the prime factors of the numbers from 1 to N in the array..
i.e prime[10] = 5, prime[14] = 7, prime[31] = 31 (prime array will contain one of the prime factor of the Number.)
Now when we will get a value N, we can factorize it in O(Number of prime factors)..
Now for a range [1,10^6] the maximum value Number of prime factors can get is <20 (as 2^20  = 10^6)
So after performing the sieve we can compute phi[n] in O( T * number of prime factors of N ).

Implementation of the above sieve -

#include<iostream>
#define N 1000000
#include<cmath>
#include<cstdio>
using namespace std;
int prime[N+1];
 
void sieve(){
    prime[0]=1;
    prime[1]=1;
    for(long long i=2;i<=N;i++){
        if(!prime[i]){
                        prime[i]=i;
            for(long long j=i;j*i<=N;j++){
                prime[j*i]=i;
            }
        }
    }
}

long long int primeFactorize(int n){
    long long int yy=1;
    for( ; n > 1 ; ){
        int p = prime[n] , e = 0 ;
        for( ; n % p == 0 ; n /= p , e++ ) ;  
        yy*= ((p-1)*(pow((double)p,e-1)));
    }
    return yy;
}


int main()
{
   sieve();
   int t;
   scanf("%d",&t);
   while(t--)
   {
    int n;
    scanf("%d",&n);
    long long int ll = primeFactorize(n);
    printf("%lld\n",ll);              
   } 
    
   return 0;
}
 
In the second implementation we are going to precompute all the values of phi[n].
It must be implemented when the number of testcases are high. Because after precomputation, the time required to answer the query becomes O(1).
But in precomputation the time required is O(nloglogn).
In this sieve we again perform standard sieve by combining the result of phi of earlier values relatively prime to current value.. and updating phi with the primes stored in an array..

Here is the implementation of the second Sieve -
 
#include<iostream>
#include<vector>
#include<cmath>
#include<cstdio>
using namespace std;

#define N 1000000
#define MAXN 1000000
int phi[MAXN + 1], prime[MAXN/10], sz=0;
vector<bool> mark(MAXN + 1);

int main()
{
    
   phi[1] = 1; 
   for (int i = 2; i <= MAXN; i++ ){
    if(!mark[i]){
        phi[i] = i-1;
        prime[sz++]= i;
    }
    for (int j=0; j<sz && prime[j]*i <= MAXN; j++ ){
        mark[prime[j]*i]=1;
        if(i%prime[j]==0){
            int ll = 0;int xx = i;
            while(xx%prime[j]==0)
            {
                           xx/=prime[j];
                           ll++;         
                       }
            int mm = 1;
            for(int k=0;k<ll;k++)mm*=prime[j];
            phi[i*prime[j]] = phi[xx]*mm*(prime[j]-1);
            break;
    }
        else phi[i*prime[j]] = phi[i]*(prime[j]-1 );
    }
}
   int t;
   scanf("%d",&t);
   while(t--)
   {
    int n;
    scanf("%d",&n);
    printf("%d\n",phi[n]);             
   } 
   return 0;    
} 
 
All of this will get more clear by reading
1. Multiplicative Function
2. Sieve of Erasthothenes

Questions to Try -
https://www.spoj.com/problems/ETF/
https://www.spoj.com/problems/DCEPCA03/