SPOJ Pell (Mid pelling)
Pell (Mid pelling) is related to Pell’s equation. It is similar to Project Euler 66 and SPOJ EQU2. I just wrote the quick solution from mathworld but I found Chakravala method very interesting. Accepted Haskell code.
import qualified Data.ByteString.Char8 as BS
import Data.List
import Data.Maybe ( fromJust )
continuedFraction :: Integer -> [ Integer ]
continuedFraction n = map ( \ ( a , _ , _ ) -> a ) . iterate fun $ ( d , 0 , 1 )
where
d = truncate . sqrt . fromIntegral $ n
fun ( a0 , p0 , q0 ) = ( a1 , p1 , q1 ) where
p1 = a0 * q0 - p0
q1 = div ( n - p1 ^ 2 ) q0
a1 = div ( d + p1 ) q1
pellSolver :: Integer -> BS.ByteString
pellSolver n
| perfectSqr n = BS.pack. show $ ( -1 )
| otherwise = ( BS.pack . show $ p ) `BS.append` ( BS.pack " " )
`BS.append` ( BS.pack.show $ q ) where
d = truncate . sqrt . fromIntegral $ n
lst = takeWhile ( /= 2 * d ) . continuedFraction $ n
len = length lst
r@( x : y : xs ) = take ( if even len then len else 2 * len ) .
continuedFraction $ n
( p0 , p1 , q0 , q1 ) = ( x , x * y + 1 , 1 , y )
( p , _ ) = foldl' compute ( p1 , p0 ) $ xs
( q , _ ) = foldl' compute ( q1 , q0 ) $ xs
compute :: ( Integer , Integer ) -> Integer -> ( Integer , Integer )
compute ( p1 , p0 ) a = ( a * p1 + p0 , p1 )
perfectSqr :: Integer -> Bool
perfectSqr n = d * d == n where
d = truncate . sqrt . fromIntegral $ n
readI :: BS.ByteString -> Integer
readI = fst . fromJust . BS.readInteger
main = BS.interact $ BS.unlines . map ( pellSolver . readI ) . tail . BS.lines
Project Euler 197
Project Euler 197 was totally based on intuition ( its still working although now a days I am not very active ). After doing bit of math, I saw the pattern and rest of work was done by Haskell ( iterate function ). I just wrote the function f and ran
ghci>take 4 . drop 0 . iterate f $ ( -1 ) [-1.0,0.7100000000000001,1.001242148,0.708777686] ghci>take 4 . drop 10 . iterate f $ ( -1 ) [1.005832407,0.7042658820000001,1.0068861020000002,0.7032313570000001] ghci>take 4 . drop 500 . iterate f $ ( -1 ) [1.029461825,0.6811758910000001,1.029461827,0.681175889]
Haskell source code
import Data.List f :: Double -> Double f x = 1E-9 * ret where a = 30.403243784 y = a - x ^ 2 ret = fromIntegral ( floor $ 2 ** y ) recur :: [ Double ] -> Double recur ( a : b : c : xs ) | abs ( a - c ) <= 1E-10 = a + b | otherwise = recur ( b : c : xs ) solve :: Double solve = recur . iterate f $ ( - 1 ) main :: IO () main = print $ solve
Project Euler Problem 107 ( Minimum Spanning Tree )
Project Euler Problem 107 is related to minimum spanning tree. Trying to solve some easy problems to cross the 200 mark. This is the first time I have implemented the graph algorithm in Haskell so it was kind of fun. Wrote Prim’s algorithm to calculated the minimum spanning tree. Replaced comma ( , ) in input file by space. See modified input file and total sum of inputs from ideone.
import Data.List
import qualified Data.IntMap.Lazy as M
import qualified Data.PQueue.Min as PQ
import qualified Data.Bits as B
buildGraph :: [ ( Int , Int , Int ) ] -> M.IntMap [ ( Int , Int ) ]
buildGraph xs = M.fromListWith ( ++ ) . map f_1 $ xs where
f_1 ( i , j , val ) = ( i , [ ( val , j ) ] )
visited :: Int -> Integer -> ( Bool , Integer )
visited x vis = ( t == 0 , vis' ) where
t = ( B..&. ) ( B.shiftL ( 1 :: Integer ) x ) vis
vis' = ( B..|. ) ( B.shiftL ( 1 :: Integer ) x ) vis
primAlgorithm :: PQ.MinQueue ( Int , Int ) -> Integer -> Int -> M.IntMap [ ( Int , Int ) ] -> String
primAlgorithm queue vis cost m
| PQ.null queue = show cost ++ "\n"
| t == False = primAlgorithm queue' vis' cost m
| otherwise = primAlgorithm queue'' vis' ( cost + c ) m where
( ( c , x ) , queue' ) = PQ.deleteFindMin queue
( t , vis' ) = visited x vis
queue'' = PQ.union queue' ( PQ.fromList . ( M.! ) m $ x )
test :: [ ( Int , Int , Int ) ] -> String
test xs = primAlgorithm ( PQ.fromList [ ( 0 , 0 ) ] ) 0 0 ( buildGraph xs )
parseInput :: [ ( Int , String ) ] -> [ ( Int , Int ) ]
parseInput [] = []
parseInput ( ( k , x ) : xs )
| x == "-" = parseInput xs
| otherwise = ( k , readD x ) : parseInput xs where
readD :: String -> Int
readD = read
fun :: Int -> [ [ ( Int , Int ) ] ] -> [ [ ( Int , Int , Int ) ] ]
fun _ [] = []
fun k ( x : xs ) = ( map ( \( i , val ) -> ( k , i , val ) ) x ) : fun ( succ k ) xs
main = interact $ primAlgorithm ( PQ.fromList [ ( 0 , 0 ) ] ) 0 0 . buildGraph . concat . fun 0 . map ( parseInput . zip [ 0 , 1 .. ] . words ) . lines
Macintosh-0026bb610428:Programming mukesh$ ghc-7.4.1 –make -O2 Prim.hs
[1 of 1] Compiling Main ( Prim.hs, Prim.o )
Linking Prim …
Macintosh-0026bb610428:Programming mukesh$ ./Prim < network.txt
2153
Project Euler Problem 381
Project Euler Problem 381 is easy problem. Using Wilson theorem
Haskell source code using monad par
import Data.List
import Data.Numbers.Primes
import Control.Monad.Par
import GHC.Conc.Sync
--use wilson theorem ( p - 1 ) ! = -1 mod p
prime = drop 2 . takeWhile ( < 10 ^ 8 ) $ primes
extendedGcd :: Int -> Int -> ( Int , Int )
extendedGcd a b
| b == 0 = ( 1 , 0 )
| otherwise = ( t , s - q * t ) where
( q , r ) = quotRem a b
( s , t ) = extendedGcd b r
modInv :: Int -> Int -> Int
modInv a b
| gcd a b /= 1 = error " gcd is not 1 "
| otherwise = d where
d = until ( > 0 ) ( + b ) . fst.extendedGcd a $ b
computeFinal :: Int -> Int
computeFinal p = computeS ( p - 1 ) 0 1 where
computeS inv ret cnt
| cnt >= 5 = mod ( ret + inv ) p
| otherwise = computeS ( mod ( inv * inv' ) p ) ( mod ( ret + inv ) p ) ( cnt + 1 ) where
inv' = modInv ( p - cnt ) p
chunk :: [ Int ] -> Int -> [ [ Int ] ]
chunk xs n = ret where
tot = numCapabilities
brk = div n tot
ret = unfoldr fun ( xs , 1 ) where
fun ( xs' , cnt )
| null xs' = Nothing
| cnt >= tot = Just ( xs' , ( [] , cnt ) )
| otherwise = Just ( a , ( b , cnt + 1 ) ) where
( a , b ) = splitAt brk xs'
ans = ret' where
ret = runPar $ parMap ( map computeFinal ) ( chunk prime len )
len = length prime
ret' = sum . map sum $ ret
main = do
print ans
Compile it “ghc-7.4.1 -O2 -rtsopts -threaded -fforce-recomp Euler_381.hs” and run it with as many cores you have. My run time is
[mukesh@user Euler]$ time ./Euler_381 +RTS -A16M -H2G -N8 -RTS
139602943319822
real 0m28.473s
user 0m58.204s
sys 0m21.713s
Project Euler 357
Project Euler 357 is easy one. I was trying to solve problem 356 and finally ended up with solving easy problem. Nothing new just sieve of eratosthenes and couple of constraints.
Edit: Finally wrote a Haskell solution. See more discussion on Haskell-Cafe. Compile this code with ghc –make -O2 filename.hs
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed
import Control.Monad
import Data.List
prime :: Int -> UArray Int Bool
prime n = runSTUArray $ do
arr <- newArray ( 2 , n ) True :: ST s ( STUArray s Int Bool )
forM_ ( takeWhile ( \x -> x*x <= n ) [ 2 .. n ] ) $ \i -> do
ai <- readArray arr i
when ( ai ) $ forM_ [ i^2 , i^2 + i .. n ] $ \j -> do
writeArray arr j False
return arr
pList :: UArray Int Bool
pList = prime $ 10 ^ 8
divPrime :: Int -> Bool
divPrime n = all ( \d -> if mod n d == 0 then pList ! ( d + div n d ) else True ) $ [ 1 .. truncate . sqrt . fromIntegral $ n ]
main = putStrLn . show . sum $ ( [ if and [ pList ! i , divPrime . pred $ i ] then ( fromIntegral . pred $ i ) else 0 | i <- [ 2 .. 10 ^ 8 ] ] :: [ Integer ] )
#include<cstdio>
#include<iostream>
#include<vector>
#define Lim 100000001
using namespace std;
bool prime [Lim];
vector<int> v ;
void isPrime ()
{
for( int i = 2 ; i * i <= Lim ; i++)
if ( !prime [i]) for ( int j = i * i ; j <= Lim ; j += i ) prime [j] = 1 ;
for( int i = 2 ; i <= Lim ; i++) if ( ! prime[i] ) v.push_back( i ) ;
//cout<<v.size()<<endl;
//for(int i=0;i<10;i++) cout<<v[i]<<" ";cout<<endl;
}
int main()
{
isPrime();
int n = v.size();
long long sum = 0;
for(int i = 0 ; i < n ; i ++)
{
int k = v[i]-1;
bool f = 0;
for(int i = 1 ; i*i<= k ; i++)
if ( k % i == 0 && prime[ i + ( k / i ) ] ) { f=1 ; break ; }
if ( !f ) sum += k;
}
cout<<sum<<endl;
}
Project Euler Problem 207
import Data.List import Data.Ratio import Data.Bits recur::Integer ->Integer recur n | (%) ( truncate.logBase 2.fromIntegral $ n ) ( n - 1 ) < (%) 1 12345 = n | otherwise = recur ( n + 1 ) solve::Integer solve = x^2 - x where x = recur 2 main = print solve
Project Euler Problem 207 is really nice and cryptic problem. Most of time was spend to understand the problem statement. We can write the given equation 
Project Euler 146
Brute forced with bit of analysis. After carefully observation , the numbers should be multiple of 5. For n = 1 or 4 mod 5, 
Edit:: Forum suggests ” n must be congruent to 10, 80, 130 or 200 mod 210 [2,3,5,7]” so this could be make bit more faster using this fact.Changed the loop from for(_int i=5;i<150000000;i+=5) to for(_int i=10;i<150000000;i+=10) makes it bit faster.
#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;
int main()
{
_int sum=0,k;
for(_int i=10;i<150000000;i+=10)
{
k=( i*i+1 );
if(mpz_probab_prime_p(k.get_mpz_t(),5) )
{
_int next , curr = k;
mpz_nextprime(next.get_mpz_t(),curr.get_mpz_t());
if( next == ( i*i + 3 ) )
{
curr= next ;
mpz_nextprime(next.get_mpz_t() , curr.get_mpz_t());
if( next == ( i*i + 7 ) )
{
curr = next;
mpz_nextprime(next.get_mpz_t() , curr.get_mpz_t());
if ( next == ( i*i + 9 ) )
{
curr = next;
mpz_nextprime(next.get_mpz_t() , curr.get_mpz_t());
if(next == ( i*i + 13 ) )
{
curr = next;
mpz_nextprime(next.get_mpz_t() , curr.get_mpz_t());
if( next == ( i*i + 27 )) sum += i;
}
}
}
}
}
}
cout<<sum<<endl;
}
SPOJ 6556. N DIV PHI_N
SPOJ 6556. N DIV PHI_N is same as Project Euler problem 69.
import Data.List
import qualified Data.ByteString.Char8 as BS
isprime :: [Integer]
isprime = 2:3:filter prime [2*k+1 | k<-[2..]]
prime :: Integer -> Bool
prime n = all ((/=0). mod n ).takeWhile ( <= ( truncate.sqrt.fromIntegral $ n) ) $ isprime
helpSolve :: [Integer] -> Integer -> Integer -> Integer
helpSolve (x:xs) n ret
| ret*x > n = ret
| otherwise = helpSolve xs n ( x * ret )
solve :: Integer -> BS.ByteString
solve n = BS.pack.show $ a where
a = helpSolve isprime n 1
readInt :: BS.ByteString -> Integer
readInt x = case BS.readInteger x of
Just ( i , _ ) -> i
Nothing -> error " upseperable ints "
main = BS.interact $ BS.unlines.map ( solve .readInt ) . BS.lines
Project Euler 329
Project Euler 329 is related to probability. Just follow the problem in Haskell because of recursive nature of problem and Haskell both 
. I was trying to solve it using dynamic programming to make it faster but could not succeed so just pure recursion. If some how i could update the Map by ” prod ” [ M.insert ( i , (x:xs) ) prod m ] and use it for further calculations. Any way its take couple of minutes on my slow PC .
import Data.List
import qualified Data.Map as M
import Data.Ratio
isprime :: Integer -> Bool
isprime n | n<=1 = False
| otherwise = all ( (/=0) . mod n ) . takeWhile (\x -> x*x<=n ) $ [2..]
m_1 , m_2 , m :: M.Map ( Integer , String ) ( Ratio Integer )
m_1 = M.fromList[ if isprime i then ( (i,"P") , (%) 2 3 ) else ( (i,"P") , (%) 1 3 ) | i <-[1..500]]
m_2 = M.fromList[ if isprime i then ( (i,"N") , (%) 1 3 ) else ( (i,"N") , (%) 2 3 ) | i <-[1..500]]
m = M.union m_1 m_2
calprob::Integer -> String -> M.Map ( Integer , String ) ( Ratio Integer ) -> ( Ratio Integer )
calprob i [x] m =
case M.lookup ( i , [x] ) m of
Just t -> t
_ -> error " some thing wrong with your preprocessing\n"
calprob i (x:xs) m =
let
Just t = M.lookup ( i , [x] ) m
prod = ( if i == 1 then t * calprob ( i+1 ) xs m
else
if i == 500 then t*calprob ( i-1 ) xs m
else t * ( % ) 1 2 * ( calprob ( i - 1 ) xs m + calprob ( i + 1 ) xs m ) )
in prod
main = putStrLn . show $ ( (%) 1 500 ) * ( sum . map (\x -> calprob x "PPPPNNPPPNPPNPN" m ) $ [1..500] )
Project Euler 105
At least my brute force solution for problem 103 worked for problem 105 . Replaced all the ‘,’ from input file by space. Same code except reading file by standard input . Takes couple of minutes probably 10 or 15.
import Data.List
subset::[a] ->[[a]]
subset [] = [[]]
subset (x:xs) = subset xs ++ map (x:) (subset xs)
fun ::(Num a, Ord a ) => [a] -> [[a]] -> Bool
fun _ [] = True
fun y (x:xs) =
case intersect y x == [] of
True -> if and[ s_a /= s_b , if l_a > l_b then s_a > s_b else if l_a < l_b then s_a < s_b else True] then fun y xs else False where
s_a = sum y
s_b = sum x
l_a = length y
l_b = length x
_ -> fun y xs
checkOptimum ::(Num a , Ord a ) => [a] -> Bool
checkOptimum xs = optimumHelper ( tail.subset $ xs ) where
optimumHelper [] = True
optimumHelper (y:ys) = if fun y ys then optimumHelper ys else False
solve :: [[Integer]] -> Integer
solve [] = 0
solve ( x : xs ) = if checkOptimum x then sum x + solve xs else solve xs
final ::[[Integer]] -> String
final xs = (show.solve $ xs) ++ "\n"
main = interact $ final . map ( map read . words ) . lines
Edit: Now with improved algorithm , its almost instant answer. I am keeping the original post just to remind me the power of algorithm and thought. I sorted the subsets based on their sum. Now the check is linear and keep checking two adjacent subsets for given condition in problem.
import Data.List
subset::[a] ->[[a]]
subset [] = [[]]
subset (x:xs) = subset xs ++ map (x:) (subset xs)
fun ::(Num a, Ord a ) => [a] -> [a] -> Bool
fun y x =
case intersect y x == [] of
True -> if ( if l_a > l_b then s_a > s_b else if l_a < l_b then s_a < s_b else s_a /= s_b ) then True else False where
s_a = sum y
s_b = sum x
l_a = length y
l_b = length x
_ -> True
checkOptimum ::(Num a , Ord a ) => [a] -> Bool
checkOptimum xs = optimumHelper (sortBy (\a b -> compare ( sum a ) ( sum b) ) . tail.subset $ xs ) where
optimumHelper [x,y] = if fun x y then True else False
optimumHelper (x:y:ys) = if fun x y then optimumHelper (y:ys) else False
solve :: [[Integer]] -> Integer
solve [] = 0
solve ( x : xs ) = if checkOptimum x then sum x + solve xs else solve xs
final ::[[Integer]] -> String
final xs = (show.solve $ xs) ++ "\n"
main = interact $ final . map ( map read . words ) . lines
About
Hello All, My name is Mukesh Tiwari and I graduated from IIITM Gwalior in 2009. I am passionate about programming and mathematics. I love solving problems on SPOJ , UVA , Topcoder and Project Euler. I am curious about functional languages specially Haskell and it’s one the excellent language I encountered after python. I am also looking for some challenging functional programming job specially related to Haskell. You can reach me mukeshtiwariDOTiiitmATgmailDOTcom. Replace DOT with . and AT with @
Recent Comments
