# My Weblog

## 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


March 23, 2013

## 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


August 10, 2012

## 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

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 )
Macintosh-0026bb610428:Programming mukesh$./Prim < network.txt 2153 August 4, 2012 ## Project Euler Problem 381 Project Euler Problem 381 is easy problem. Using Wilson theorem $1\cdot 2\cdots ( p - 1 )\ \equiv\ -1\ \pmod{p}$ $1\cdot 2\cdots ( p - 2 )\ \equiv\ \frac{-1}{p-1} \pmod{p}$ $1\cdot 2\cdots ( p - 3 )\ \equiv\ \frac{-1}{( p-1) (p-2)} \pmod{p}$ $1\cdot 2\cdots ( p - 4 )\ \equiv\ \frac{-1}{( p-1) (p-2)(p-3)} \pmod{p}$ $1\cdot 2\cdots ( p - 5 )\ \equiv\ \frac{-1}{( p-1) (p-2)(p-3) ( p-4 )} \pmod{p}$ 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

April 28, 2012

## 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 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
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;
}


November 7, 2011

## 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 $2^{2*t} - 2^t = K$. Its given that $2^t$ is positive integer so K will be of form $x^2 - x$.When K will be power of 2 [ t will be positive integer ] then we have perfect partition. In short we have to find the $\frac {\log_2 n } { n-1 } < \frac 1 {12345}$. June 22, 2011 ## 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, $n^2=1 \mod 5$ then $n^2+9=0 \mod 5$. Similarly for n = 2 or 3 mod 5 $n^2=4 \mod 5$ then $n^2+1=0 \mod 5$. It takes couple of minutes on my office’s faster computer. 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; }  June 22, 2011 ## SPOJ 6556. N DIV PHI_N SPOJ 6556. N DIV PHI_N is same as Project Euler problem 69.$\frac n {\phi (n) } = \frac { p_0p_1...p_r} { (p_0-1)(p_1-1)...(p_r-1) }$ and this quantity will be maximum when $p_0*p_1....p_r$ is maximum and $(p_0-1)*(p_1-1)...(p_r-1)$ is minimum.$1 -\frac 1 p$ will be minimal for small p hence $\frac p {p-1}$ will maximal so take product of the primes and product of primes should not exceed N. Accepted Haskell code.You can also have look here. 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

Just ( i , _ ) -> i
Nothing -> error " upseperable ints "

main = BS.interact $BS.unlines.map ( solve .readInt ) . BS.lines  June 17, 2011 ## 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] )



May 29, 2011

## 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


May 26, 2011