My Weblog

Blog about programming and math

Project Euler problem 133

I did not use too much math for Project Euler problem 133. I just had doubt about upper bound which i again set too high (10^1000 for c program and 10^500 for python). I simply check the fact that 10^(10^n)!=1(mod 9p). Initially i took only p and keep getting wrong answer.Then i put 9p 😀 and got the right answer. Now its time to learn some math about 129 ( which still unsolved) , 132 and 133.
compile it g++ proj_133.cpp -lgmpxx -lgmp

#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <queue>
#include <stack>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <string>
#include <cstring>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include<gmpxx.h>
typedef mpz_class _int;
using namespace std;

vector<_int> P;

void isprime()
	{
		bool b[100000];
		memset(b,0,sizeof b);
		for(int i=2;i*i<100000;i++) if(!b[i]) for(int j=i*i;j<100000;j+=i)b[j]=1;
		for(int i=2;i<100000;i++) if(!b[i])P.push_back(_int(i));
	}

int main()
	{	
		isprime();
		int n=P.size();
		_int s=0,exp=1;
		for(int i=1;i<=10000;i++) exp*=10;
		_int rop,base=10,t;
		for(int i=0;i<n;i++) 
		 {
			t=9*P[i];//10^(10^n)!=1(mod 9p)
			mpz_powm(rop.get_mpz_t(),base.get_mpz_t(),exp.get_mpz_t(),t.get_mpz_t());
			if( rop!=1) s+=P[i]; 
		 }
		cout<<s<<endl;
	}

python code for which i wrote first but was getting wrong answer and increasing the exp variable upto 1000 was taking too much time so i had to switch for c++. This code also works well under one minute rule.

# To change this template, choose Tools | Templates
# and open the template in the editor.

__author__="user"
__date__ ="$Jul 25, 2010 12:35:45 AM$"
import random
def rabin_miller(p):
	if(p<2):
		return False
	if(p!=2 and p%2==0):
		return False
	s=p-1
	while(s%2==0):
		s>>=1
	for i in xrange(10):
		a=random.randrange(p-1)+1
		temp=s
		mod=pow(a,temp,p)
		while(temp!=p-1 and mod!=1 and mod!=p-1):
			mod=(mod*mod)%p
			temp=temp*2
		if(mod!=p-1 and temp%2==0):
			return False
	return True

if __name__ == "__main__":
    sum,exp=0,10**500
    for p in range(2,100000,1):
        if(rabin_miller(p) and  (pow(10,exp,9*p)!=1)):
            sum+=p;
    print sum;

Advertisements

July 24, 2010 - Posted by | Programming

No comments yet.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: