Personal tools

Euler problems/171 to 180

From HaskellWiki

< Euler problems(Difference between revisions)
Jump to: navigation, search
m
 
(12 intermediate revisions by 8 users not shown)
Line 2: Line 2:
 
Finding numbers for which the sum of the squares of the digits is a square.
 
Finding numbers for which the sum of the squares of the digits is a square.
   
Solution:
+
{{sect-stub}}
<haskell>
 
 
This does not seem Haskell code to me.
 
If the argument: Learning Haskell were valid pure Haskell code would have been given.
 
#include <stdio.h>
 
 
static int result = 0;
 
 
#define digits 20
 
static long long fact[digits+1];
 
static const long long precision = 1000000000;
 
static const long long precision_mult = 111111111;
 
 
#define maxsquare 64 /* must be a power of 2 > digits * 9^2 */
 
 
static inline int issquare( int n )
 
{
 
for( int step = maxsquare/2, i = step;;)
 
{
 
if( i*i == n ) return i;
 
if( !( step >>= 1 ) ) return -1;
 
if( i*i > n ) i -= step;
 
else i += step;
 
}
 
}
 
 
static inline void dodigit( int d, int nr, int sum, long long c, int s )
 
{
 
if( d )
 
for( int n = 0; n <= nr; c *= ++n, s += d, sum += d*d )
 
dodigit( d-1, nr - n, sum, c, s );
 
else if( issquare( sum ) > 0 )
 
result = ( s * ( fact[digits] / ( c * fact[nr] ) )
 
/ digits % precision * precision_mult
 
+ result ) % precision;
 
}
 
 
int main( void )
 
{
 
fact[0] = 1;
 
for( int i = 1; i < digits+1; i++ ) fact[i] = fact[i-1]*i;
 
dodigit( 9, digits, 0, 1, 0 );
 
printf( "%d\n", result );
 
return 0;
 
}
 
problem_171 = main
 
</haskell>
 
   
 
== [http://projecteuler.net/index.php?section=problems&id=172 Problem 172] ==
 
== [http://projecteuler.net/index.php?section=problems&id=172 Problem 172] ==
Line 41: Line 41:
 
Counting the number of "hollow" square laminae that can form one, two, three, ... distinct arrangements.
 
Counting the number of "hollow" square laminae that can form one, two, three, ... distinct arrangements.
   
Solution: This was my C++ code, published here without my permission nor any attribution, shame on whoever put it here. [[user:henk263|henk263]]
+
{{sect-stub}}
   
 
== [http://projecteuler.net/index.php?section=problems&id=175 Problem 175] ==
 
== [http://projecteuler.net/index.php?section=problems&id=175 Problem 175] ==
Line 59: Line 59:
 
l=length k
 
l=length k
 
p175 x y=
 
p175 x y=
init$foldl (++) "" [a++","|
+
concat $ intersperse "," $
a<-map show $reverse $filter (/=0)$findRat x y]
+
map show $reverse $filter (/=0)$findRat x y
 
problems_175=p175 123456789 987654321
 
problems_175=p175 123456789 987654321
 
test=p175 13 17
 
test=p175 13 17
Line 84: Line 84:
 
Integer angled Quadrilaterals.
 
Integer angled Quadrilaterals.
   
Solution:
+
{{sect-stub}}
  +
  +
== [http://projecteuler.net/index.php?section=problems&id=178 Problem 178] ==
  +
Step Numbers
  +
  +
Count pandigital step numbers.
 
<haskell>
 
<haskell>
This does not seem Haskell code to me.
+
import qualified Data.Map as M
If the argument: Learning Haskell were valid pure Haskell code would have been given.
 
   
#include <stdio.h>
+
data StepState a = StepState { minDigit :: a
#include <math.h>
+
, maxDigit :: a
// gcc --std c99 -lm 177.c
+
, lastDigit :: a
int isint(double x);
+
} deriving (Show, Eq, Ord)
double fabs(double x);
 
   
long long int count;
+
isSolution (StepState i a _) = i == 0 && a == 9
double duparray[100][20];
+
neighborStates m s@(StepState i a n) = map (\x -> (x, M.findWithDefault 0 s m)) $
  +
[StepState (min i (n - 1)) a (n - 1), StepState i (max a (n + 1)) (n + 1)]
   
double PIeight=3.14159265358979323846264338327950288419716939937510/180.0;
+
allStates = [StepState i a n | (i, a) <- range ((0, 0), (9, 9)), n <- range (i, a)]
int main()
+
initialState = M.fromDistinctAscList [(StepState i i i, 1) | i <- [1..9]]
{
+
stepState m = M.fromListWith (+) $ allStates >>= neighborStates m
double x,w,m;
+
numSolutionsInMap = sum . map snd . filter (isSolution . fst) . M.toList
int a,b,c,d;
+
numSolutionsOfSize n = sum . map numSolutionsInMap . take n $ iterate stepState initialState
double maxxvalue;
 
   
int iopt;
+
problem_178 = numSolutionsOfSize 40
int I,j;
+
</haskell>
int N;
 
count=0;
 
double sine[200];
 
double cosine[200];
 
for(int i=1;i<=180;i++)
 
{
 
if(i<=90)
 
{
 
sine[i]=sin(PIeight*(double)i);
 
cosine[i]=cos(PIeight*(double)i);
 
}
 
else
 
{
 
sine[i]=sine[180-i];
 
cosine[i]=-cosine[180-i];
 
}
 
}
 
   
for(int alpha=1;alpha<=45;alpha++)
+
== [http://projecteuler.net/index.php?section=problems&id=179 Problem 179] ==
{
+
Consecutive positive divisors.
for(int beta=alpha;beta<180-alpha;beta++)
+
See http://en.wikipedia.org/wiki/Divisor_function for a simple
{
+
explanation of calculating the number of divisors of an integer,
  +
based on its prime factorization. Then, if you have a lot of
  +
time on your hands, run the following program. You need to load
  +
the Factoring library from David Amos' wonderful Maths library.
  +
See: http://www.polyomino.f2s.com/david/haskell/main.html
   
for(int gamma=alpha;gamma<=180-alpha-beta-alpha;gamma++)
+
<haskell>
{
+
import Factoring
w=sine[alpha+beta]/sine[alpha+beta+gamma];
+
import Data.List
int b=180-alpha-beta-gamma;
 
for(int delta=alpha;delta<=180-gamma-beta-alpha;delta++)
 
{
 
x=sine[beta]/sine[beta+gamma+delta];
 
m=sqrt(w*w+x*x-2*w*x*cosine[delta]);
 
if(x*sine[delta]>m)
 
a=90;
 
else
 
a=(int)(round)(1.0/PIeight*asin(x*sine[delta]/m));
 
if(m*m+w*w-x*x<0)
 
a=180-a;
 
d=180-beta-gamma-delta;
 
   
c=360-a-b-alpha-beta-gamma-delta-d;
+
nFactors :: Integer -> Int
  +
nFactors n =
  +
let a = primePowerFactorsL n
  +
in foldl' (\x y -> x * ((snd y)+1) ) 1 a
   
if(a>=alpha && c>=alpha && a+b+c+d+alpha+beta+gamma+delta==360)
+
countConsecutiveInts l = foldl' (\x y -> if y then x+1 else x) 0 a
{
+
where a = zipWith (==) l (tail l)
if(fabs((sine[delta]*sine[c]/(sine[a]*sine[d]))-
 
(sine[gamma]*sine[alpha]/(sine[beta]*sine[b])))<1e-11)
 
{
 
duparray[1][1]=(double)alpha;
 
duparray[1][2]=(double)beta;
 
duparray[1][3]=(double)gamma;
 
duparray[1][4]=(double)delta;
 
duparray[1][5]=(double)d;
 
duparray[1][6]=(double)c;
 
duparray[1][7]=(double)a;
 
duparray[1][8]=(double)b;
 
for(I=1;I<=3;I++)
 
for(j=1;j<=8;j++)
 
duparray[I+1][j]=duparray[1][(j+I*2-1)%8+1];
 
for(j=1;j<=8;j++)
 
duparray[5][9-j]=duparray[1][j];
 
   
for(I=1;I<=3;I++)
+
problem_179 = countConsecutiveInts $ map nFactors [2..(10^7 - 1)]
for(j=1;j<=8;j++)
 
duparray[I+5][j]=duparray[5][(j+I*2-1)%8+1];
 
N=8;
 
maxxvalue=1e22;
 
iopt=1;
 
for(I=1;I<=N;I++)
 
{
 
duparray[I][9]=0;
 
for(j=1;j<=8;j++)
 
duparray[I][9]=duparray[I][9]*180+duparray[I][j];
 
   
if(duparray[I][9]<maxxvalue-1e-7)
+
main = print problem_179
{
 
maxxvalue=duparray[I][9];
 
iopt=I;
 
}
 
}
 
   
  +
</haskell>
   
  +
This is all well and good, but it runs very slowly. (About 4
  +
minutes on my machine) We have to factor every number between 2
  +
and 10^7, which on a non Quantum CPU takes a while. There is
  +
another way!
   
if(iopt==1)
+
We can sieve for the answer. Every number has 1 for a factor.
count++;
+
Every other number has 2 as a factor, and every third number has
}
+
3 as a factor. So we run a sieve that counts (increments)
}
+
itself for every integer. When we are done, we run through the
  +
resulting array and look at neighboring elements. If they are
  +
equal, we increment a counter. This version runs in about 9
  +
seconds on my machine. HenryLaxen May 14, 2008
   
  +
<haskell>
  +
{-# OPTIONS -O2 -optc-O #-}
  +
import Data.Array.ST
  +
import Data.Array.Unboxed
  +
import Control.Monad
  +
import Control.Monad.ST
  +
import Data.List
   
}
+
r1 = (2,(10^7-1))
}
 
}
 
}
 
   
printf("%lld\n",count);
+
type Sieve s = STUArray s Int Int
   
}
+
incN :: Sieve s -> Int -> ST s ()
problem_177 = main
+
incN a n = do
</haskell>
+
x <- readArray a n
  +
writeArray a n (x+1)
   
== [http://projecteuler.net/index.php?section=problems&id=178 Problem 178] ==
+
incEveryN :: Sieve s -> Int -> ST s ()
Step Numbers
+
incEveryN a n = mapM_ (incN a) [n,n+n..snd r1]
Solution:
+
<haskell>
+
sieve :: Int
This does not seem Haskell code to me.
+
sieve = countConsecutiveInts b
If the argument: Learning Haskell were valid pure Haskell code would have been given.
+
where b = runSTUArray $
  +
do a <- newArray r1 1 :: ST s (STUArray s Int Int)
  +
mapM_ (incEveryN a) [fst r1 .. (snd r1) `div` 2]
  +
return a
   
#include <stdio.h>
+
countConsecutiveInts :: UArray Int Int -> Int
#include <math.h>
+
countConsecutiveInts a =
#define N 40
+
let r1 = [fst (bounds a) .. snd (bounds a) - 1]
+
in length $ filter (\y -> a ! y == a ! (y+1)) $ r1
double f[50][11][11][11];
+
+
main = print sieve
int main()
 
{
 
int x,y,z,i,j,k,m;
 
 
for(m=1;m<=N;m++)
 
{
 
for(i=0;i<=9;i++)
 
for(j=0;j<=9;j++)
 
for(k=0;k<=9;k++)
 
{
 
if(i==j && j==k && m==1)
 
f[m][i][j][k]=1;
 
else
 
f[m][i][j][k]=0;
 
}
 
}
 
 
for(m=2;m<=N;m++)
 
{
 
for(x=0;x<=9;x++)
 
{
 
 
for(y=x+1;y<=9;y++)
 
for(z=x;z<=y;z++)
 
{
 
if(z>x && z<y)
 
{
 
f[m][x][y][z]=f[m-1][x][y][z-1]+f[m-1][x][y][z+1];
 
}
 
if(z==x)
 
{
 
f[m][x][y][z]=f[m-1][x][y][x+1]+f[m-1][x+1][y][x+1];
 
}
 
if(z==y)
 
{
 
f[m][x][y][z]=f[m-1][x][y][y-1]+f[m-1][x][y-1][y-1];
 
}
 
}
 
}
 
}
 
double count=0;
 
for(i=1;i<=N;i++)
 
{
 
for(z=1;z<=9;z++)
 
count+=f[i][0][9][z];
 
}
 
printf("%lf\n",count);
 
}
 
problem_178 = main
 
 
</haskell>
 
</haskell>
   
== [http://projecteuler.net/index.php?section=problems&id=179 Problem 179] ==
 
Consecutive positive divisors.
 
 
{{sect-stub}}
 
   
 
== [http://projecteuler.net/index.php?section=problems&id=180 Problem 180] ==
 
== [http://projecteuler.net/index.php?section=problems&id=180 Problem 180] ==

Latest revision as of 01:49, 13 February 2010

Contents

[edit] 1 Problem 171

Finding numbers for which the sum of the squares of the digits is a square.

[edit] 2 Problem 172

Investigating numbers with few repeated digits.

Solution:

factorial n = product [1..toInteger n]
 
fallingFactorial x n = product [x - fromInteger i | i <- [0..toInteger n - 1] ]
 
choose n k = fallingFactorial n k `div` factorial k
 
-- how many numbers can we get having d digits and p positions
p172 0 _ = 0
p172 d p 
    | p < 4 = d^p
    | otherwise = 
        (p172' p) +  p*(p172' (p-1)) + (choose p 2)*(p172' (p-2)) + (choose p 3)*(p172' (p-3))
    where
    p172' = p172 (d-1)
 
problem_172= (p172 10 18) * 9 `div` 10

[edit] 3 Problem 173

Using up to one million tiles how many different "hollow" square laminae can be formed? Solution:

problem_173=
    let c=div (10^6) 4
        xm=floor$sqrt $fromIntegral c
        k=[div c x|x<-[1..xm]]
    in  sum k-(div (xm*(xm+1)) 2)

[edit] 4 Problem 174

Counting the number of "hollow" square laminae that can form one, two, three, ... distinct arrangements.

[edit] 5 Problem 175

Fractions involving the number of different ways a number can be expressed as a sum of powers of 2. Solution:

sternTree x 0=[]
sternTree x y=
    m:sternTree y n  
    where
    (m,n)=divMod x y
findRat x y
    |odd l=take (l-1) k++[last k-1,1]
    |otherwise=k
    where
    k=sternTree x y
    l=length k
p175 x y= 
    concat $ intersperse "," $
    map show $reverse $filter (/=0)$findRat x y
problems_175=p175 123456789 987654321
test=p175 13 17

[edit] 6 Problem 176

Rectangular triangles that share a cathetus. Solution:

--k=47547 
--2*k+1=95095 = 5*7*11*13*19
lst=[5,7,11,13,19]
primes=[2,3,5,7,11]
problem_176 =
    product[a^b|(a,b)<-zip primes (reverse n)]
    where
    la=div (last lst+1) 2
    m=map (\x->div x 2)$init lst
    n=m++[la]

[edit] 7 Problem 177

Integer angled Quadrilaterals.

[edit] 8 Problem 178

Step Numbers

Count pandigital step numbers.

import qualified Data.Map as M
 
data StepState a = StepState { minDigit  :: a
                             , maxDigit  :: a
                             , lastDigit :: a
                             } deriving (Show, Eq, Ord)
 
isSolution (StepState i a _) = i == 0 && a == 9
neighborStates m s@(StepState i a n) = map (\x -> (x, M.findWithDefault 0 s m)) $
    [StepState (min i (n - 1)) a (n - 1), StepState i (max a (n + 1)) (n + 1)]
 
allStates    = [StepState i a n | (i, a) <- range ((0, 0), (9, 9)), n <- range (i, a)]
initialState = M.fromDistinctAscList [(StepState i i i, 1) | i <- [1..9]]
stepState m  = M.fromListWith (+) $ allStates >>= neighborStates m
numSolutionsInMap    = sum . map snd . filter (isSolution . fst) . M.toList
numSolutionsOfSize n = sum . map numSolutionsInMap . take n $ iterate stepState initialState
 
problem_178 = numSolutionsOfSize 40

[edit] 9 Problem 179

Consecutive positive divisors. See http://en.wikipedia.org/wiki/Divisor_function for a simple explanation of calculating the number of divisors of an integer, based on its prime factorization. Then, if you have a lot of time on your hands, run the following program. You need to load the Factoring library from David Amos' wonderful Maths library. See: http://www.polyomino.f2s.com/david/haskell/main.html

import Factoring
import Data.List
 
nFactors :: Integer -> Int
nFactors n =
  let a = primePowerFactorsL n
  in foldl' (\x y -> x * ((snd y)+1) ) 1 a
 
countConsecutiveInts l = foldl' (\x y -> if y then x+1 else x) 0 a
  where a = zipWith (==) l (tail l)
 
problem_179 =  countConsecutiveInts $ map nFactors [2..(10^7 - 1)]
 
main = print problem_179

This is all well and good, but it runs very slowly. (About 4 minutes on my machine) We have to factor every number between 2 and 10^7, which on a non Quantum CPU takes a while. There is another way!

We can sieve for the answer. Every number has 1 for a factor. Every other number has 2 as a factor, and every third number has 3 as a factor. So we run a sieve that counts (increments) itself for every integer. When we are done, we run through the resulting array and look at neighboring elements. If they are equal, we increment a counter. This version runs in about 9 seconds on my machine. HenryLaxen May 14, 2008

{-# OPTIONS -O2 -optc-O #-}
import Data.Array.ST
import Data.Array.Unboxed
import Control.Monad
import Control.Monad.ST
import Data.List
 
r1 = (2,(10^7-1))
 
type Sieve s = STUArray s Int Int
 
incN :: Sieve s -> Int -> ST s ()
incN a n = do
  x <- readArray a n
  writeArray a n (x+1)
 
incEveryN :: Sieve s -> Int -> ST s ()
incEveryN a n = mapM_ (incN a)  [n,n+n..snd r1]
 
sieve :: Int
sieve = countConsecutiveInts b
  where b = runSTUArray $
            do a <- newArray r1 1 :: ST s (STUArray s Int Int)
               mapM_ (incEveryN a) [fst r1 .. (snd r1) `div` 2]
               return a
 
countConsecutiveInts :: UArray Int Int -> Int                        
countConsecutiveInts a =
  let r1 = [fst (bounds a) .. snd (bounds a) - 1]
  in length $ filter (\y -> a ! y == a ! (y+1)) $ r1
 
main = print sieve


[edit] 10 Problem 180

Rational zeros of a function of three variables. Solution:

import Data.Ratio
 
{-
  After some algebra, we find:
   f1 n x y z = x^(n+1) + y^(n+1) - z^(n+1)
   f2 n x y z = (x*y + y*z + z*x) * ( x^(n-1) + y^(n-1) - z^(n-1) )
   f3 n x y z = x*y*z*( x^(n-2) + y^(n-2) - z^(n-2) )
   f n x y z = f1 n x y z + f2 n x y z - f3 n x y z
   f n x y z = (x+y+z) * (x^n+y^n-z^n)
Now the hard part comes in realizing that n can be negative.  
Thanks to Fermat, we only need examine the cases n = [-2, -1, 1, 2]
Which leads to:
 
f(-2) z = xy/sqrt(x^2 + y^2)
f(-1) z = xy/(x+y)
f(1)  z = x+y
f(2)  z = sqrt(x^2 + y^2)
 
-}
 
unique ::  Eq(a) => [a] -> [a]
unique [] = []
unique (x:xs) | elem x xs = unique xs
              | otherwise = x : unique xs
 
-- Not quite correct, but I don't care about the zeros
ratSqrt :: Rational -> Rational
ratSqrt x = 
  let a = floor $ sqrt $ fromIntegral $ numerator x
      b = floor $ sqrt $ fromIntegral $ denominator x
      c = (a%b) * (a%b)
  in if x == c then (a%b) else 0
 
-- Not quite correct, but I don't care about the zeros
reciprocal :: Rational -> Rational
reciprocal x 
  | x == 0 = 0
  | otherwise = denominator x % numerator x
 
problem_180 =
  let order = 35
      range :: [Rational]            
      range = unique [ (a%b) | b <- [1..order], a <- [1..(b-1)] ]
      fm2,fm1,f1,f2 :: [[Rational]]
      fm2 = [[x,y,z] | x<-range, y<-range, 
            let z = x*y * reciprocal (ratSqrt(x*x+y*y)), elem z range]
      fm1 = [[x,y,z] | x<-range, y<-range, 
            let z = x*y * reciprocal (x+y), elem z range]
      f1  = [[x,y,z] | x<-range, y<-range, 
            let z = (x+y), elem z range]
      f2  = [[x,y,z] | x<-range, y<-range, 
            let z = ratSqrt(x*x+y*y), elem z range]            
      result = sum $ unique $ map (\x -> sum x) (fm2++fm1++f1++f2)
  in (numerator result + denominator result)