Showing posts with label Lang_C#. Show all posts
Showing posts with label Lang_C#. Show all posts

Thursday, September 17, 2009

XNA quest for centroid of polygon

Now we will talk about XNA Framework and How to find geometrical center of polygon (aka. so called centroid). At first - Why do we need centroid at all ? Well... for example - for some graphics applications we may want drag&drop shape by it`s centroid. Or maybe in some game applications we may need to find shape`s center of gravity, for simulating physics of balancing objects around center of gravity. And in cases when shape material density is uniform (the same over shape) then shape geometrical center (centroid) is also center of gravity (or center of mass). So lets begin search of centroid... For this you will need: - Visual Studio 2008 - XNA 3.1 Framework First - some notes about XNA framework. It`s a pitty that such great 2D/3D .NET library doesn`t have 2D primitives drawing routines. Well, actually it has some, but for being able to draw some 2D primitives like points and lines, you should dig deeply in 3D knowledge - for drawing mentioned primitives in MSDN way you should define 3D vectors and 3D viewports !! What the hell ? If user only wants 2D stuff, it means that all 3D stuff should be abstracted out of the view. So that at the end, framework should have separate levels of abstraction into 2D and 3D. But this is not the case about drawing 2D primtives... So one way of drawing point in non-MSDN way is to generate small rectangle texture in run-time and draw this generated texture on screen. Next thing is drawing a simple line in non-MSDN way. Recipe for that - we take generated 2D rectangle texture earlier and draw it as line by scaling and rotating it as needed. After defining this stuff we have all prerequisites for target problem. Now- How we will find centroid of polygon ? This is not that hard. According to wikipedia centroid formulas can be computed simply, assuming we have polygon vertex coordinates and assuming polygon is not complex. So centroid computation algorithm in C# is this:
public static Vector2 CalculateCentroid(Vector2[] points, 
                                        int lastPointIndex)
{
    float area = 0.0f;
    float Cx = 0.0f;
    float Cy = 0.0f;
    float tmp = 0.0f;
    int k;

    for (int i = 0; i <= lastPointIndex; i++)
    {
        k = (i + 1) % (lastPointIndex + 1);
        tmp = points[i].X * points[k].Y - 
              points[k].X * points[i].Y;
        area += tmp;
        Cx += (points[i].X + points[k].X) * tmp;
        Cy += (points[i].Y + points[k].Y) * tmp;
    }
    area *= 0.5f;
    Cx *= 1.0f / (6.0f * area);
    Cy *= 1.0f / (6.0f * area);

    return new Vector2(Cx, Cy);
}

NOTE: this algorithm doesn`t apply for complex polygons. But How simply check that polygon is not complex, or not self-intersecting ? One way of doing this is brute-force way -> we can check every polygon edge pair for intersection. If there are some edges which intersects, then polygon is complex, otherwise - simple. So by using this brute-force solution, complex polygon check simplifies to lines self-intersection problem. And How do we check if lines self-intersects ? One way is to find-out 2 linear equations of these lines: y = a1*x + b1 ; y = a2*x + b2. Then self-intersection point is simply where these equations equals: a1*x + b1 = a2*x + b2 --> x_intersection = (b2-b1)/(a1-a2). Substituting into any of these two equations we get y_intersection = (a1*b2-a2*b1)/(a1-a2). Next step - we need to check does lines segments intersects, not the whole lines. This is done just comparing does x_intersection is within x_min,x_max of both lines, the similar check apply for y_intersection. So line intersection check simplifies to finding linear equations slope and intercept constants. This is one way of self-intersection test of lines. But there is also other way to do it (if we don`t need exact intersection point) - suppose we have line1 with end-points l1p1,l1p2 and line2 with end-points l2p1,l2p2. Define function DIRECTION(point1,point2,point3) which returns the direction [CLOCKWISE,COUNTER_CLOCKWISE,LINE] in which we traverse path point1->point2->point3. By having such function we can easily define when 2 line segments intersects. Below are one example 
when lines intersects:
    

Another example when lines do not intersects:  

So lines intersects when : direction(l1p1,l1p2,l2p1) != direction(l1p1,l1p2,l2p2) AND direction(l2p1,l2p2,l1p1) != direction(l2p1,l2p2,l1p2) This type of lines intersection test is implemented in polygon self-intersection test. So you can download centroid XNA 3.1 project and try experimenting by drawing various polygons and seeing where centroid is computed. Final part - some results from this project (centroid is draw as black pixel): Some convex polygons:  

 
 
Now, more interesting -> concave polygons:

   
 
From the concave polygon examples (especially from first concave) can be seen that centroid is not the average of (max,min) coordinates in shape. However in some cases it can be that centroid coincides with AVERAGE(Min,Max) of figure (for example for rectangle). But in general - it not. If to talk informally, then centroid is simply the average of ALL points coordinates in shape. Finally just some examples of complex polygons - as evidence that application detects complex (or self-intersecting) polygons :-) : 

 

Actually we can also compute centroids of complex polygon by using info from above mentioned wikipedia site. That is - to get centroid of complex polygon we need: 1. Decompose complex polygon into N simple polygons. 2. Compute area and centroid of each simple polygon by method already explained in this article. 3. Average computed centroids by some formula. The main task is number 1. - How effectively split complex polygon into simple ones. This task is left as exercise to the reader :-) Have fun with centroids, polygons, XNA and 2D/3D programming in general !!

Tuesday, June 10, 2008

Project Euler p96 - solving sudoku

This is interesting one. Problem description:
-----------------------------------------------
Su Doku (Japanese meaning number place) is the name given to a popular puzzle concept. Its origin is unclear, but credit must be attributed to Leonhard Euler who invented a similar, and much more difficult, puzzle idea called Latin Squares. The objective of Su Doku puzzles, however, is to replace the blanks (or zeros) in a 9 by 9 grid in such that each row, column, and 3 by 3 box contains each of the digits 1 to 9. Below is an example of a typical starting puzzle grid and its solution grid.

0 0 3 0 2 0 6 0 0
9 0 0 3 0 5 0 0 1
0 0 1 8 0 6 4 0 0

0 0 8 1 0 2 9 0 0
7 0 0 0 0 0 0 0 8
0 0 6 7 0 8 2 0 0

0 0 2 6 0 9 5 0 0
8 0 0 2 0 3 0 0 9
0 0 5 0 1 0 3 0 0
===================
4 8 3 9 2 1 6 5 7
9 6 7 3 4 5 8 2 1
2 5 1 8 7 6 4 9 3

5 4 8 1 3 2 9 7 6
7 2 9 5 6 4 1 3 8
1 3 6 7 9 8 2 4 5

3 7 2 6 8 9 5 1 4
8 1 4 2 5 3 7 6 9
6 9 5 4 1 7 3 8 2

A well constructed Su Doku puzzle has a unique solution and can be solved by logic, although it may be necessary to employ "guess and test" methods in order to eliminate options (there is much contested opinion over this). The complexity of the search determines the difficulty of the puzzle; the example above is considered easy because it can be solved by straight forward direct deduction.

The 6K text file, sudoku.txt (right click and 'Save Link/Target As...'), contains fifty different Su Doku puzzles ranging in difficulty, but all with unique solutions (the first puzzle in the file is the example above).

By solving all fifty puzzles find the sum of the 3-digit numbers found in the top left corner of each solution grid; for example, 483 is the 3-digit number found in the top left corner of the solution grid above.
-----------------------------------------------

Well, back to more understandable language like C# :-). Solution is done as semi brute-force, iterative method. Solution is quick enough - fifty puzzles are solved in just 122 ms. So this is it:



class Program
{

static int[][][] ReadSudoku(string file)
{
StreamReader sr = new StreamReader(file);
int[][][] res = new int[50][][];
int row = 0;
int sud = -1;
string line;

while (!sr.EndOfStream)
{
line = sr.ReadLine();
if (line.ToLower().IndexOf("grid") > -1)
{
row = 0;
sud++;
res[sud] = new int[9][];
}
else
{
res[sud][row] = Array.ConvertAll(line.ToCharArray(), c =>
{ return int.Parse(c.ToString()); });
row++;
}
}
sr.Close();

return res;
}

static void ShowSudoku(int[][] sud)
{
string o;
o = sud.Aggregate("", (agg, next) =>
{
return agg +
((agg == "") ? "":"\n") +
next.Aggregate("|", (agg2, next2) => {
return agg2 + next2.ToString(); }) + "|";
});
o = "-----------\n" + o + "\n-----------";
o = o.Replace("0", ".");
Console.Write(o);
}

static bool ValueIsValid(int[][] grid, int val, int row, int col)
{
int regc, regr;
int lc, lr;

if (val == 0)
return false;

// checking row and column
for (int i = 0; i < 9; i++)
if ((grid[row][i] == val && i != col) ||
(grid[i][col] == val && i != row))
return false;

// checking 3x3 region
regr = 3*(row / 3);
regc = 3*(col / 3);
for (int i = 0; i < 3; i++)
{
lc = regc + i;
for (int j = 0; j < 3; j++)
{
lr = regr + j;
if (lc != col || lr != row)
if (grid[lr][lc] == val)
return false;
}
}

return true;
}

static void CellWithMinValid(int[][] grid,
out int rowo,
out int colo)
{
int row = -1, col = -1;
int minvalid = 9;
int valcount;

for (int i = 0; i < grid.Length; i++)
{
for (int j = 0; j < grid[i].Length; j++)
{
if (grid[i][j] == 0)
{
valcount = 0;
for (int k = 1; k < 10; k++)
if (ValueIsValid(grid, k, i, j))
valcount++;
if (valcount < minvalid)
{
minvalid = valcount;
row = i;
col = j;
}
}
if (minvalid == 1) break;
}
if (minvalid == 1) break;
}

rowo = row;
colo = col;
}

static void SolveSudoku(int[][] grid)
{
int row, col;
int curval;
int iter;
bool force;
int[,] steps = new int[81, 2];

CellWithMinValid(grid, out row, out col);
grid[row][col] = 1;
steps[0, 0] = row;
steps[0, 1] = col;
force = false;
iter = 0;

while (row > -1)
{
curval = grid[steps[iter, 0]]
[steps[iter, 1]];
if (ValueIsValid(grid, curval,
steps[iter, 0],
steps[iter, 1]) && !force)
{
CellWithMinValid(grid, out row, out col);
if (row != -1)
{
force = false;
iter++;
grid[row][col] = 1;
steps[iter, 0] = row;
steps[iter, 1] = col;
}
}
else
{
curval++;
if (curval < 10)
{
grid[steps[iter, 0]]
[steps[iter, 1]] = curval;
force = false;
}
else
{
grid[steps[iter, 0]]
[steps[iter, 1]] = 0;
force = true;
iter--;
}
}
}
}

static void SolveAll(string file)
{
int[][][] sud = ReadSudoku(file);
int tot = 0;

for (int i = 0; i < sud.Length; i++)
{
SolveSudoku(sud[i]);
tot += (sud[i][0][0] * 100) +
(sud[i][0][1] * 10) +
(sud[i][0][2]);
}

Console.WriteLine("Answer to projectEuler {0}\nLast sudoku solution:\n",tot);
ShowSudoku(sud[49]);
}

static void Main(string[] args)
{
System.Diagnostics.Stopwatch sw =
new System.Diagnostics.Stopwatch();
sw.Start();
SolveAll("C:\\ProjectEuler\\sudoku.txt");
sw.Stop();
Console.WriteLine("\n50-ty sudoku`s solved in {0} ms",
sw.ElapsedMilliseconds);
Console.ReadLine();
}
}



As always - have fun!.

Wednesday, May 21, 2008

Project Euler p102 - finding origin

Problem description:
-----------------------------------------------------
Three distinct points are plotted at random on a Cartesian plane, for which -1000 ≤ x, y ≤ 1000, such that a triangle is formed.

Consider the following two triangles:

A(-340,495), B(-153,-910), C(835,-947)

X(-175,41), Y(-421,-714), Z(574,-645)

It can be verified that triangle ABC contains the origin, whereas triangle XYZ does not.

Using triangles.txt (right click and 'Save Link/Target As...'), a 27K text file containing the co-ordinates of one thousand "random" triangles, find the number of triangles for which the interior contains the origin.

NOTE: The first two examples in the file represent the triangles in the example given above.
-----------------------------------------------------

I`ve exploited the idea that for origin to be in the triangle interior, triangle lines segments should cross x and y axis 2 times: x+ / x- and y+ / y-.
Basically it is enough to check does triangle cross y+ and y- axis parts, but I didn`t know that and checked both axis. This worked too :-)

Solution ( C# as always):



using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;

namespace ConsoleApplication1
{
class Program
{
static int[][] ReadTriangleCoordinates(string file)
{
int[][] res;
StreamReader sr = new StreamReader(file);
res = Array.ConvertAll(sr.ReadToEnd().TrimEnd().Split("\n".ToCharArray()),
s => {return Array.ConvertAll(s.Split(",".ToCharArray()),
n => { return int.Parse(n); }) ;});
sr.Close();
return res;
}

static void CrossesAxisAt(int x1,
int y1,
int x2,
int y2,
out double? xAt,
out double? yAt)
{
int dx = x2 - x1;
int dy = y2 - y1;
int xmin = Math.Min(x1, x2);
int ymin = Math.Min(y1, y2);
int xmax = Math.Max(x1, x2);
int ymax = Math.Max(y1, y2);
double? xcross;
double? ycross;
double a;
double b;

if (dx == 0)
{
xcross = x1;
ycross = null;
}
else if (dy == 0)
{
xcross = null;
ycross = y1;
}
else
{
a = (double)dy / (double)dx;
b = y1 - (a * x1);
xcross = -(b / a);
ycross = b;
}

if (xcross != null)
{
if (xcross < xmin || xcross > xmax
|| ymin > 0 || ymax < 0)
xcross = null;
}

if (ycross != null)
{
if (ycross < ymin || ycross > ymax
|| xmin > 0 || xmax < 0)
ycross = null;
}

xAt = xcross;
yAt = ycross;
}

static bool ContainsOrigin(int[] abc)
{
double? xc1, xc2, xc3;
double? yc1, yc2, yc3;
int xposc = 0, xnegc = 0, yposc = 0, ynegc = 0;
bool res;

CrossesAxisAt(abc[0], abc[1], abc[4], abc[5], out xc1, out yc1);
CrossesAxisAt(abc[0], abc[1], abc[2], abc[3], out xc2, out yc2);
CrossesAxisAt(abc[2], abc[3], abc[4], abc[5], out xc3, out yc3);

xposc = ((xc1 > 0) ? 1 : 0) + ((xc2 > 0) ? 1 : 0) + ((xc3 > 0) ? 1 : 0);
xnegc = ((xc1 < 0) ? 1 : 0) + ((xc2 < 0) ? 1 : 0) + ((xc3 < 0) ? 1 : 0);
yposc = ((yc1 > 0) ? 1 : 0) + ((yc2 > 0) ? 1 : 0) + ((yc3 > 0) ? 1 : 0);
ynegc = ((yc1 < 0) ? 1 : 0) + ((yc2 < 0) ? 1 : 0) + ((yc3 < 0) ? 1 : 0);

if (xposc > 0 && xnegc > 0 && yposc > 0 && ynegc > 0)
res = true;
else
res = false;

return res;
}

static int TrianglesWithOrigin(string file)
{
int[][] tr = ReadTriangleCoordinates(file);
int trcount = 0;

for (int i = 0; i < tr.Length; i++)
{
if (ContainsOrigin(tr[i]))
trcount++;
}

return trcount;
}

static void Main(string[] args)
{
Console.WriteLine(TrianglesWithOrigin("C:\\ProjectEuler\\triangles.txt"));
Console.ReadLine();
}
}
}


Sunday, May 18, 2008

IBM Ponder This, May 2008 - cellular automaton

From the IBM project "ponder this", may 2008:
-------------------------------------------------
This month's puzzle is about Conway's Game of life.
The game is a cellular automaton on a grid where in each generation, every cell determines its state (dead or alive) based on its eight neighbors: a dead cell will change its state to alive if and only if exactly three of its neighbors are alive; A live cell will stay alive if and only if two or three of its neighbors are alive.
All the changes are done simultaneously.
Find a possible state of the universe whose next generation is depicted below.


....................
....................
..###..###...#...#..
...#....#.#..##.##..
...#....##...#.#.#..
...#....#.#..#...#..
..###..###...#...#..
....................
....................


"." denotes dead cells and "#" denotes live cells.
-------------------------------------------------

So after a bit of thinking i chose Selfish Gene Algorithm, which I also presented on my blog article. Please note - I will present solution which will use 4 classes from the given my blog link above. So if you need to test this solution - download the code from that article also and import these classes:
- SGAallele.cs
- SGAgenotype.cs
- SGAlocus.cs
- SGAlocusGroup.cs
Now about solution. In general - solution with the help of selfish gene algorithm, tries to evolve cellular automaton configuration which next generation is similar to the one described in problem definition. Fitness is calculated as wrong cell amount - that is cell, which is not in place of target configuration, amount. Algorithm tries to minimize that wrong cell amount. So solution is this:



using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SGA;

namespace ConsoleApplication18
{
class Program
{
static int counter = 0;
static bool[][] grid = InitialGrid();
static SGAgenotype gen;

static void ShowCells(bool[][] cells)
{
string outp;
outp = cells.Aggregate("", (acc1, ba) =>
acc1 + ba.Aggregate("", (acc2, b) =>
acc2 + ((b) ? "#" : ".")) + "\n");
Console.Write(outp);
}

static bool[][] InitialGrid()
{
bool[][] c = new bool[9][];
c[0] = new bool[20];
c[1] = new bool[20];
c[2] = new bool[] { false, false, true, true, true, false, false,
true, true, true, false, false, false, true,
false, false, false, true, false, false };
c[3] = new bool[] { false, false, false, true, false, false, false,
false, true, false, true, false, false, true,
true, false, true, true, false, false };
c[4] = new bool[] { false, false, false, true, false, false, false,
false, true, true, false, false, false, true,
false, true, false, true, false, false };
c[5] = new bool[] { false, false, false, true, false, false, false,
false, true, false, true, false, false, true,
false, false, false, true, false, false };
c[6] = new bool[] { false, false, true, true, true, false, false, true,
true, true, false, false, false, true, false,
false, false, true, false, false };
c[7] = new bool[20];
c[8] = new bool[20];

return c;
}

static int LiveCells(bool[][] grid, int row, int col)
{
int live = 0;

for (int i = row - 1; i <= row + 1; i++)
{
for (int j = col - 1; j <= col + 1; j++)
{
if (i != row || j != col)
{
if (grid[i][j])
live++;
}
}
}

return live;
}

static bool WillLive(bool IsLive, int LiveNeighbours)
{
bool state;

if (IsLive)
state = (LiveNeighbours < 2 || LiveNeighbours > 3) ? false : true;
else
state = (LiveNeighbours == 3) ? true : false;

return state;
}

static bool[][] ComputeStates(bool[][] grid)
{
bool[][] res = new bool[9][];
res[0] = new bool[20];
res[1] = new bool[20];
res[2] = new bool[20];
res[3] = new bool[20];
res[4] = new bool[20];
res[5] = new bool[20];
res[6] = new bool[20];
res[7] = new bool[20];
res[8] = new bool[20];

for (int i = 1; i < grid.Length - 1; i++)
{
for (int j = 1; j < grid[i].Length - 1; j++)
{
res[i][j] = WillLive(grid[i][j], LiveCells(grid, i, j));
}
}

return res;
}

static bool[][] CellGridFromGenotype(SGAgenotype gen, SGAgenotype.GenomeNo num)
{
int row, col;
bool v;
int gene = -1;
bool[][] res = new bool[9][];
res[0] = new bool[20];
res[1] = new bool[20];
res[2] = new bool[20];
res[3] = new bool[20];
res[4] = new bool[20];
res[5] = new bool[20];
res[6] = new bool[20];
res[7] = new bool[20];
res[8] = new bool[20];

for (int i = 0; i < gen.locusGroups.Length; i++)
{
row = (i) / 20;
col = (i) - (row * 20);

if (num == SGAgenotype.GenomeNo.FirstGenome)
gene = gen.locusGroups[i].loci[0].FirstGeneIndex;
else if (num == SGAgenotype.GenomeNo.SecondGenome)
gene = gen.locusGroups[i].loci[0].SecondGeneIndex;
else if (num == SGAgenotype.GenomeNo.None)
gene = gen.locusGroups[i].loci[0].BestGeneIndex;

v = Convert.ToBoolean(gen.locusGroups[i].loci[0].Alleles[gene].Value);

res[row][col] = v;
}

return res;
}

static int WrongCellCount(bool[][] previousStates)
{
int res = 0;
bool[][] nextStates = ComputeStates(previousStates);

for (int i = 0; i < grid.Length; i++)
{
for (int j = 0; j < grid[i].Length; j++)
{
if (nextStates[i][j] != grid[i][j])
res++;
}
}

return res;
}

static SGAgenotype.CompareResults CompareCellAutomata()
{
SGAgenotype.CompareResults cr;
bool[][] grid1 = CellGridFromGenotype(gen, SGAgenotype.GenomeNo.FirstGenome);
bool[][] grid2 = CellGridFromGenotype(gen, SGAgenotype.GenomeNo.SecondGenome);
bool[][] gridb = CellGridFromGenotype(gen, SGAgenotype.GenomeNo.None);

int wcellsGene1 = WrongCellCount(grid1);
int wcellsGene2 = WrongCellCount(grid2);
int wcellsGeneb = WrongCellCount(gridb);

cr.CompareEachOther = SGAgenotype.GenomeNo.None;
cr.CompareWithBest = SGAgenotype.GenomeNo.None;

if (wcellsGene1 < wcellsGene2)
cr.CompareEachOther = SGAgenotype.GenomeNo.FirstGenome;
else if (wcellsGene2 < wcellsGene1)
cr.CompareEachOther = SGAgenotype.GenomeNo.SecondGenome;

if (wcellsGene1 < wcellsGeneb)
cr.CompareWithBest = SGAgenotype.GenomeNo.FirstGenome;
else if (wcellsGene2 < wcellsGeneb)
cr.CompareWithBest = SGAgenotype.GenomeNo.SecondGenome;

counter++;

if (counter%1000==0)
{
Console.Clear();
ShowCells(gridb);
Console.WriteLine("====================");
ShowCells(ComputeStates(gridb));
Console.WriteLine("====================\nPlease wait...");
}

return cr;
}

static void Main(string[] args)
{
gen = new SGAgenotype(new SGAlocus[] { new SGAlocus(0, 1, 1) },
180, CompareCellAutomata);
gen.AlleleAffectValue *= 0.1;
gen.MutationProbability *= 1.5;
gen.Evolve(400000);
Console.WriteLine("DONE.");
Console.ReadLine();
}
}
}



After execution we get this picture-


...............#....
.#....#.....#.....##
...##...##..........
...#....#....##.##.#
....#....##..#...#.# ANSWER
..#.....#.....#..#..
...#....##..#.....#.
...#..#.....#...#...
....................
====================
....................
....................
..###..###...#...#..
...#....#.#..##.##..
...#....##...#.#.#.. NEXT GENERATION OF CELLS
...#....#.#..#...#..
..##...###...#...#..
....................
....................
====================


First picture is answer and second (below) is next generation of cells.
We see that this is approximate answer because next generation differs by 1 cell from the given situation in problem definition. Well, this is because I`ve used genetic algorithm which by definition is probabilistic and doesn`t guarantee to give exact solution. But very well calculates approximate solutions.
So all in all - it was fun to play with selfish gene algorithm and cellular automaton.

Tuesday, May 13, 2008

Traveler's dilemma and bonus factor

Traveler's dilemma description is following (taken from wikipedia):
-------------------------------------------------------
An airline loses two suitcases belonging to two different travelers. Both suitcases happen to be identical and contain identical antiques. An airline manager tasked to settle the claims of both travelers explains that the airline is liable for a maximum of $100 per suitcase, and in order to determine an honest appraised value of the antiques the manager separates both travelers so they can't confer, and asks them to write down the amount of their value at no less than $2 and no larger than $100. He also tells them that if both write down the same number, he will treat that number as the true dollar value of both suitcases and reimburse both travelers that amount. However, if one writes down a smaller number than the other, this smaller number will be taken as the true dollar value, and both travelers will receive that amount along with a bonus/malus: $2 extra will be paid to the traveler who wrote down the lower value and a $2 deduction will be taken from the person who wrote down the higher amount. The challenge is: what strategy should both travelers follow to decide the value they should write down?
-------------------------------------------------------

Experiment was conducted as following:
- BONUS was chosen between 2 and 20
- Possible travelers request was between BONUS and 100.
- Travelers chooses cost randomly (for no-strategy imitation)
- Cost random selection was repeated 5000000 times
- After that optimal request was calculated as travelers choise which gives the greatest profit.
- And finally these optimal requests are plot in diagram as function of BONUS.

C# code which conducted experiments on random travelers choise is this:




class Program
{
static void Payoff(int Aoffer, int Boffer, int extra, out int Apayoff, out int Bpayoff)
{
int Aextra, Bextra;
int minoffer = Math.Min(Aoffer, Boffer);

if (Aoffer == Boffer)
{
Aextra = Bextra = 0;
}
else
{
Aextra = (Aoffer == minoffer) ? extra : -extra;
Bextra = (Boffer == minoffer) ? extra : -extra;
}
Apayoff = minoffer + Aextra;
Bpayoff = minoffer + Bextra;
}

static long OptimalSelection(int costFrom, int costTo, int goods)
{
long[] payoffA = new long[costTo + 1];
long[] payoffB = new long[costTo + 1];
Random r = new Random();
int Ac, Bc;
int Pa, Pb;
long Aopt = -1, Bopt = -1;
long Amax = 0, Bmax = 0;

for (int i = 1; i <= goods; i++)
{
Ac = r.Next(costFrom, costTo + 1);
Bc = r.Next(costFrom, costTo + 1);
Payoff(Ac, Bc, costFrom, out Pa, out Pb);
payoffA[Ac] += Pa;
payoffB[Bc] += Pb;
}

for (int j = costFrom; j <= costTo; j++)
{
if (payoffA[j] > Amax) { Amax = payoffA[j]; Aopt = j; }
if (payoffB[j] > Bmax) { Bmax = payoffB[j]; Bopt = j; }
}

return (Aopt + Bopt) / 2;
}

static void MeasureBonusInfluence()
{
double[] bonus = Array.ConvertAll(Enumerable.Range(2, 20).ToArray(),
v => (double)(v));
double[] optimal = Array.ConvertAll(bonus,
v => (double)OptimalSelection((int)v, 100, 5000000));

GoogleChart.PlotLineChart(bonus,
optimal,
"Bonus ($)",
"Optimal request ($)",
"Traveler's dilemma - bonus factor",
"400x200");
}

static void Main(string[] args)
{
MeasureBonusInfluence();
}
}

class GoogleChart
{
public static void PlotLineChart(double[] Xvalues,
double[] Yvalues,
string Xlabel,
string Ylabel,
string Title,
string imgsize)
{
string urltemplate = "http://chart.apis.google.com/chart?cht=lxy&chs=@imgsize@&chd=t:@xvalues@%7C@yvalues@&chds=@xmin@,@xmax@,@ymin@,@ymax@&chxr=0,@xmin@,@xmax@%7C1,@ymin@,@ymax@&chxt=x,y,x,y&chxl=2:%7C@xlabel@%7C3:%7C@ylabel@%7C&chxp=2,50%7C3,50&chtt=@title@";
string xvalagg = Xvalues.Aggregate("", (old, next) =>
old + ((old == "") ? "" : ",") + next.ToString(CultureInfo.InvariantCulture));
string yvalagg = Yvalues.Aggregate("", (old, next) =>
old + ((old == "") ? "" : ",") + next.ToString(CultureInfo.InvariantCulture));

if (imgsize == null) imgsize = "300x200";

urltemplate = urltemplate.Replace("@imgsize@", imgsize);
urltemplate = urltemplate.Replace("@xvalues@", xvalagg);
urltemplate = urltemplate.Replace("@yvalues@", yvalagg);
urltemplate = urltemplate.Replace("@xmin@", Xvalues.Min().ToString(CultureInfo.InvariantCulture));
urltemplate = urltemplate.Replace("@xmax@", Xvalues.Max().ToString(CultureInfo.InvariantCulture));
urltemplate = urltemplate.Replace("@ymin@", Yvalues.Min().ToString(CultureInfo.InvariantCulture));
urltemplate = urltemplate.Replace("@ymax@", Yvalues.Max().ToString(CultureInfo.InvariantCulture));
urltemplate = urltemplate.Replace("@xlabel@", Xlabel);
urltemplate = urltemplate.Replace("@ylabel@", Ylabel);
urltemplate = urltemplate.Replace("@title@", Title);

System.Diagnostics.Process.Start(urltemplate);
}
}




Finally after execution of this program, we get this picture of optimal request relation with bonus amount (if travelers chooses randomly):

Wednesday, May 7, 2008

Google chart API in .NET

Yes I know - there are good .NET libraries for using google chart API, but I was interested to dig a little bit deeper into google chart API. And I have to say that google chart is worse than opensource program GNUPLOT (in the sense of functionality), but ok whatever... And second target was - to have as less code (including libraries) as possible for charting VERY simple XY graphs. So after studing a bit google chart documentation- I wrote very primitive class for plotting simple XY graph. For plotting graphs with this code you need:
- Internet connection
- .NET framework 3.5
- Internet browser (whatever you like)
So this is the code:




class Program
{
static void Main(string[] args)
{
double[] xval;
double[] ysin, yexp;

xval = new double[] { 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5 };
ysin = Array.ConvertAll(xval, v => { return Math.Sin(v); });
yexp = Array.ConvertAll(xval, v => { return Math.Exp(v); });

GoogleChart.PlotLineChart(xval, ysin, "x values", "sin(x)", "sin(x) function", null);
GoogleChart.PlotLineChart(xval, yexp, "x values", "exp(x)", "exp(x) function", null);
}
}

class GoogleChart
{
public static void PlotLineChart(double[] Xvalues,
double[] Yvalues,
string Xlabel,
string Ylabel,
string Title,
string imgsize)
{
string urltemplate = "http://chart.apis.google.com/chart?cht=lxy&chs=@imgsize@&chd=t:@xvalues@%7C@yvalues@&chds=@xmin@,@xmax@,@ymin@,@ymax@&chxr=0,@xmin@,@xmax@%7C1,@ymin@,@ymax@&chxt=x,y,x,y&chxl=2:%7C@xlabel@%7C3:%7C@ylabel@%7C&chxp=2,50%7C3,50&chtt=@title@";
string xvalagg = Xvalues.Aggregate("", (old, next) =>
old + ((old == "") ? "" : ",") + next.ToString(CultureInfo.InvariantCulture));
string yvalagg = Yvalues.Aggregate("", (old, next) =>
old + ((old == "") ? "" : ",") + next.ToString(CultureInfo.InvariantCulture));

if (imgsize == null) imgsize = "300x200";

urltemplate = urltemplate.Replace("@imgsize@", imgsize);
urltemplate = urltemplate.Replace("@xvalues@", xvalagg);
urltemplate = urltemplate.Replace("@yvalues@", yvalagg);
urltemplate = urltemplate.Replace("@xmin@", Xvalues.Min().ToString(CultureInfo.InvariantCulture));
urltemplate = urltemplate.Replace("@xmax@", Xvalues.Max().ToString(CultureInfo.InvariantCulture));
urltemplate = urltemplate.Replace("@ymin@", Yvalues.Min().ToString(CultureInfo.InvariantCulture));
urltemplate = urltemplate.Replace("@ymax@", Yvalues.Max().ToString(CultureInfo.InvariantCulture));
urltemplate = urltemplate.Replace("@xlabel@", Xlabel);
urltemplate = urltemplate.Replace("@ylabel@", Ylabel);
urltemplate = urltemplate.Replace("@title@", Title);

System.Diagnostics.Process.Start(urltemplate);
}
}


With this sample code there are generated these nice graphs:


This is it. Have fun with charting !

Tuesday, May 6, 2008

Project Euler p59 - breaking encryption

Problem definition:
------------------------------------------------
Each character on a computer is assigned a unique code and the preferred standard is ASCII (American Standard Code for Information Interchange). For example, uppercase A = 65, asterisk (*) = 42, and lowercase k = 107.

A modern encryption method is to take a text file, convert the bytes to ASCII, then XOR each byte with a given value, taken from a secret key. The advantage with the XOR function is that using the same encryption key on the cipher text, restores the plain text; for example, 65 XOR 42 = 107, then 107 XOR 42 = 65.

For unbreakable encryption, the key is the same length as the plain text message, and the key is made up of random bytes. The user would keep the encrypted message and the encryption key in different locations, and without both "halves", it is impossible to decrypt the message.

Unfortunately, this method is impractical for most users, so the modified method is to use a password as a key. If the password is shorter than the message, which is likely, the key is repeated cyclically throughout the message. The balance for this method is using a sufficiently long password key for security, but short enough to be memorable.

Your task has been made easy, as the encryption key consists of three lower case characters. Using cipher1.txt (right click and 'Save Link/Target As...'), a file containing the encrypted ASCII codes, and the knowledge that the plain text must contain common English words, decrypt the message and find the sum of the ASCII values in the original text.
------------------------------------------------

C# brute-force solution, running in 42 seconds:



static byte[] ReadCipher(string file)
{
StreamReader sr = new StreamReader(file);
string[] buf = sr.ReadToEnd().Trim().Split(",".ToCharArray());
byte[] res = Array.ConvertAll(buf, s => { return byte.Parse(s); });
sr.Close();

return res;
}

static byte[] TryDecipher(byte[] cip, string key)
{
int ind = 0;
byte[] dec = new byte[cip.Length];

for (int i = 0; i < cip.Length; i++)
{
dec[i] = (byte)((int)cip[i] ^ (int)key[ind]);
ind = (ind >= 2) ? 0 : ind + 1;
}

return dec;
}

static int KeywordsInText(string text, string[] keywords)
{
int kc = 0;

for (int i = 0; i < keywords.Length; i++)
{
if (text.ToLower().IndexOf(" " + keywords[i] + " ") > -1) kc++;
}

return (kc*100) / keywords.Length;
}

static void Decipher(string file)
{
string[] keywords = new string[] {"the","be","is","are","was","to","of","and","a","in","into","that","have","has","had","I","it","for","not","on","with","he","as","you","do","did","done","at","this","but","his","by","from","they","we","say","her","she","or","an","will","shall","my","one","all","would","should","there","here","their","what","so","up","out","if","about","who","get","got","which","go","went","gone","me"};
string passchars = "qwertyuioplkjhgfdsazxcvbnm";
string curkey = "";
byte[] cipher = ReadCipher(file);
byte[] decoded = new byte[cipher.Length];
string dectext;
int keywordcount = 0;
int maxkeywords = 0;
string bestkey = "";
byte[] bestbytes = new byte[cipher.Length];
string besttext = "";
long bytesum = 0;

for (int i1 = 0; i1 < passchars.Length; i1++)
{
for (int i2 = 0; i2 < passchars.Length; i2++)
{
for (int i3 = 0; i3 < passchars.Length; i3++)
{
curkey = passchars[i1].ToString() + passchars[i2].ToString() + passchars[i3].ToString();
decoded = TryDecipher(cipher, curkey);
dectext = ASCIIEncoding.ASCII.GetString(decoded);
keywordcount = KeywordsInText(dectext, keywords);

if (keywordcount > maxkeywords)
{
maxkeywords = keywordcount;
bestkey = curkey;
bestbytes = decoded;
besttext = dectext;
}
}
}
}

for (int i = 0; i < bestbytes.Length; i++)
{
bytesum += (long)bestbytes[i];
}

Console.WriteLine("Byte sum: {0}\nPassword: {1}\nDecoded message 50 chars: {2}\nMessage has {3}% common words",
bytesum,bestkey,besttext.Substring(0,50),maxkeywords);
}

static void Main(string[] args)
{
System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch();
sw.Start();
Decipher("C:\\ProjectEuler\\cipher1.txt");
sw.Stop();
Console.WriteLine("Decryption time: {0} seconds",sw.Elapsed.Seconds);
Console.ReadLine();
}


Tuesday, April 29, 2008

Project Euler p79 - cracking passcode

This one is very interesting. Problem description:
-------------------------------------------
A common security method used for online banking is to ask the user for three random characters from a passcode. For example, if the passcode was 531278, they may asked for the 2nd, 3rd, and 5th characters; the expected reply would be: 317.
The text file, keylog.txt, contains fifty successful login attempts.
Given that the three characters are always asked for in order, analyse the file so as to determine the shortest possible secret passcode of unknown length.
-------------------------------------------

So solution in C# is this, not optimized, but works (in a 16 seconds):



static int[] ReadKeyLog(string file)
{
int[] dig;
string buf;
string[] barr;
StreamReader sr = new StreamReader(file);
buf = sr.ReadToEnd().TrimEnd("\n".ToCharArray());
barr = buf.Split("\n".ToCharArray());
dig = Array.ConvertAll(barr, s => { return Int32.Parse(s); });
sr.Close();
return dig;
}

static bool KeyPartExists(long testnumber, int keypart)
{
int n1 = keypart / 100;
int n2 = (keypart / 10) % 10;
int n3 = keypart % 10;
bool n1e = false;
bool n2e = false;
bool n3e = false;

for (long i = testnumber; i > 0; i /= 10)
{
if (!n3e)
n3e = (i % 10) == n3;
else
if (!n2e)
n2e = (i % 10) == n2;
else
if (!n1e)
n1e = (i % 10) == n1;
else
break;
}
if (n1e && n2e && n3e)
return true;
else
return false;
}

static bool HasAllKeyParts(long testnum, int[] keylog)
{
bool res = true;

for (int i = 0; i < keylog.Length; i++)
{
if (!KeyPartExists(testnum, keylog[i]))
{
res = false;
break;
}
}
return res;
}

static long FindPassCode(string file)
{
int[] keylog = ReadKeyLog(file);
long num = 10000;

while (true)
{
if (HasAllKeyParts(num, keylog))
return num;
num++;
}
}

static void Main(string[] args)
{
System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch();
sw.Start();
Console.WriteLine("Bruteforced passcode is {0}",
FindPassCode("C:\\ProjectEuler\\keylog.txt")
);
sw.Stop();
Console.WriteLine("Cracking time {0} seconds", sw.Elapsed.Seconds);
Console.ReadLine();
}


Thursday, April 24, 2008

Project Euler Problem 17

Definition:
---------------------------------------------
If the numbers 1 to 5 are written out in words: one, two, three, four, five, then there are 3 + 3 + 5 + 4 + 4 = 19 letters used in total.
If all the numbers from 1 to 1000 (one thousand) inclusive were written out in words, how many letters would be used?
NOTE: Do not count spaces or hyphens. For example, 342 (three hundred and forty-two) contains 23 letters and 115 (one hundred and fifteen) contains 20 letters. The use of "and" when writing out numbers is in compliance with British usage.

---------------------------------------------

C# code solution:



static long LettersOfNumber(long num, long[,]bases)
{
long count = 0;
long tempnum;
if (num > 1000 num < 1)
count = -1;
else if (num == 1000)
count = 11;
else
{
tempnum = num % 100;
if (tempnum < 20)
count += bases[tempnum, 0];
else
count += bases[num % 10, 0] + bases[(num / 10) % 10, 1];
if (num > 99)
count += bases[(num / 100), 0] + 7 + ((tempnum == 0) ? 0 : 3);
}
return count;
}
static long LettersUntilNumber(int num)
{
long letcount = 0;
long[,] bas = new long[20, 2] { {0, 0}, {3, 0}, {3, 6}, {5, 6}, {4, 5}, {4, 5}, {3, 5}, {5, 7},
{5, 6}, {4, 6}, {3, 0}, {6, 0}, {6, 0}, {8, 0}, {8, 0}, {7, 0},
{7, 0}, {9, 0}, {8, 0}, {8, 0} };
for (int i = 1; i <= num; i++)
letcount += LettersOfNumber(i, bas);
return letcount;
}
static void Main(string[]args)
{
Console.WriteLine(LettersUntilNumber(1000));
Console.ReadLine();
}


Saturday, April 19, 2008

Project Euler Problem 28

Problem definition-
--------------------------------------------------
Starting with the number 1 and moving to the right in a clockwise direction a 5 by 5 spiral is formed as follows:
21 22 23 24 25
20 7 8 9 10
19 6 1 2 11
18 5 4 3 12
17 16 15 14 13
It can be verified that the sum of both diagonals is 101.
What is the sum of both diagonals in a 1001 by 1001 spiral formed in the same way?
--------------------------------------------------
C# brute-force solution:



static long SpiralDiagonalsSum(long spiralsize)
{
long sum = 1;
long cursize = 3;
long number = 1;
while (cursize <= spiralsize)
{
for (int i = 0; i < 4; i++)
{
number += cursize - 1;
sum += number;
}
cursize += 2;
}
return sum;
}
static void Main(string[]args)
{
System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch();
sw.Start();
long diagsum = SpiralDiagonalsSum(1001);
sw.Stop();
Console.WriteLine("Spiral diagonals sum {0} calculated in {1} ms", diagsum, sw.Elapsed.TotalMilliseconds);
Console.ReadLine();
}




Solution on core duo performs in a 0.2 ms.

Friday, April 18, 2008

Selfish Gene Algorithm

This time I will present a branch of evolutionary algorithms - Selfish gene algorithm.

Scope of article
Question arises - Why we need yet another genetic algorithm (GA) ? Well, its always good to have different point of view. It`s because selfish gene algorithm (further- SGA) is a different beast than 'standard' GA. For that reasons I will present SGA and it's uses in 2 sample problems. For the article being short I will assume reader`s basic knowledge of GA.

Biological background

Several concepts will be used here. Organism genetic makeup is called - genotype. All living organisms has DNA molecules. DNA molecule is composed of genes. Gene location in DNA is called locus. Usually for the same locus there can be different versions of gene. These different versions of gene is called alleles. There can be calculated allele frequency (or probability) in population. In general evolution means, that organisms which succeeds- increases its alleles frequency in population at the expense of children. And vise-versa,- organism which fails- decreases its alleles frequency in population.

Main features of SGA
  • There are no explicit population in algorithm. Individuals are created and destroyed
    only temporary for comparision of their genomes. This means that there are no recombination operator at all. Instead new solutions is generated according to their genes frequency in population.
  • Mutation operator is encoded as random gene selection (not accordingly to frequencies)
  • We need not to assign absolute value of fitness to generated solution, but instead we
    need to compare 2 generated solutions between each other. And we need to return back to the algorithm which solution was better.
Pseudo-code of SGA
  1. Generate 2 solutions, according to solution genes frequency/probability in population.
    1a. If mutation happens - take random solution gene, not according to its freqency.
  2. Compare these 2 solutions, and decide which is better.
    2a. Reward better solution - increase its alleles frequency in population.
    2b. Penalyze worse solution - decrease its alleles frequency in population.
  3. Update best solution found if any.
  4. Repeat steps 1-3 until required solution found or other needed criteria is met-
    for example maximum evolve cycles elapsed or solutions genes frequency is above some threshold or whatever you need.
SGA class diagram



So main SGA classes is 4:

SGAgenotype
- main class for executing GA algorithm. This class is container for SGAlocus class. SGAlocus - class is container for class SGAallele, also holds info about temporary generated solution genes.
SGAlocusGroup is just for possibility of grouping SGAlocus classes, because (you will later see it) it is sometimes usefull to do it. That is acctually SGAgenotype holds SGAlocusGroup class, which in turn groups SGAlocus.
SGAallele is class for describing allele, mainly - allele value and its frequency/probability in population.
Thats it, its all you need for SGA to run on your PC.
Practical examples of SGA usage
There was 2 cases of usage of SGA.
First - box grouping problem.
Given amount of boxes of some random sizes try to arrange boxes in some area while trying to minimize total boxes intersection area and maximize free unfragmented space left in room.
Below is video of how evolves box positions in some number of iterations:

In this type of problem SGAlocusGroup groups 2 loci - one locus for box X position and second locus for box Y position. Total SGAlocusGroup count equals to box count. Alleles encodes all possible values of box positions.
Second - sudoku puzzle generation problem
Given grid size, 9x9 or 16x16, try to generate (if possible) sudoku puzzle or at least something that is close to sudoku puzzle enough :-) . Similarity is measured by duplicating cell colors, that is by error count in grid. Different numbers is represented as different colors.
Below is video how grid evolves in some iterations:

In this problem SGAlocusGroup has only one locus, which holds info about cell color.
SGAlocusGroup count equals to grid cell count. Alleles encodes all possible colors of cell.
In the 9x9 grid case often is that there are left 2-4 color duplicates. In 16x16 case there are much more. But we can adjust evolve cycles to bigger number to overcome these errors.
C# project source files
It is for .NET 2 (VS 2005) but should run in .NET 3.5 also. You can download it here and use as you wish (GNU license):
http://coding-experiments.googlecode.com/files/SelfishGeneAlgorithm.zip
Last- useful references
1. http://www.cad.polito.it/pap/db/sac98.pdf
F. Corno, M. Sonza Reorda, G. Squillero, "The Selfish Gene Algorithm: a New Evolutionary Optimization Strategy," SAC98: 13th Annual ACM Symposium on Applied Computing, Atlanta, Georgia (USA), February 1998, pp. 349-355
2. Book "The Selfish Gene", Richard Dawkins, 1976.

Additional helper links-
http://en.wikipedia.org/wiki/Genotype
http://en.wikipedia.org/wiki/Locus_(genetics)
http://en.wikipedia.org/wiki/Allele
Thats all, have fun !!!

Saturday, April 12, 2008

Project Euler Problem 18

This time we will solve euler problem 18 by semi-bruteforce solution.
Problem description:
--------------------------------------------------------------------------------------
By starting at the top of the triangle below and moving to adjacent numbers on the row below, the maximum total from top to bottom is 23.
3
7 5
2 4 6
8 5 9 3

That is, 3 + 7 + 4 + 9 = 23.
Find the maximum total from top to bottom of the triangle below:
75
95 64
17 47 82
18 35 87 10
20 04 82 47 65
19 01 23 75 03 34
88 02 77 73 07 63 67
99 65 04 28 06 16 70 92
41 41 26 56 83 40 80 70 33
41 48 72 33 47 32 37 16 94 29
53 71 44 65 25 43 91 52 97 51 14
70 11 33 28 77 73 17 78 39 68 17 57
91 71 52 38 17 14 91 43 58 50 27 29 48
63 66 04 68 89 53 67 30 73 16 69 87 40 31
04 62 98 27 23 09 70 98 73 93 38 53 60 04 23
----------------------------------------------------------------------------------------
Solution in C# is more or less short and fast (calculates answer in 1,1 ms):



using System;

class Program
{
static long MaxTotalFromTriangle()
{
long[][] triangle = new long[15][];
long max = 0;
triangle[0] = new long[] { 75 };
triangle[1] = new long[] { 95, 64 };
triangle[2] = new long[] { 17, 47, 82 };
triangle[3] = new long[] { 18, 35, 87, 10 };
triangle[4] = new long[] { 20, 04, 82, 47, 65 };
triangle[5] = new long[] { 19, 01, 23, 75, 03, 34 };
triangle[6] = new long[] { 88, 02, 77, 73, 07, 63, 67 };
triangle[7] = new long[] { 99, 65, 04, 28, 06, 16, 70, 92 };
triangle[8] = new long[] { 41, 41, 26, 56, 83, 40, 80, 70, 33 };
triangle[9] = new long[] { 41, 48, 72, 33, 47, 32, 37, 16, 94, 29 };
triangle[10] = new long[] { 53, 71, 44, 65, 25, 43, 91, 52, 97, 51, 14 };
triangle[11] = new long[] { 70, 11, 33, 28, 77, 73, 17, 78, 39, 68, 17, 57 };
triangle[12] = new long[] { 91, 71, 52, 38, 17, 14, 91, 43, 58, 50, 27, 29, 48 };
triangle[13] = new long[] { 63, 66, 04, 68, 89, 53, 67, 30, 73, 16, 69, 87, 40, 31 };
triangle[14] = new long[] { 04, 62, 98, 27, 23, 09, 70, 98, 73, 93, 38, 53, 60, 04, 23 };

for (int i = 1; i < triangle.Length; i++)
{
for (int j = 0; j < triangle[i].Length; j++)
{
// Accumulating maximum total
if (j == 0)
{
triangle[i][j] += triangle[i - 1][j];
}
else if (j == triangle[i].Length - 1)
{
triangle[i][j] += triangle[i - 1][triangle[i - 1].Length - 1];
}
else
{
triangle[i][j] += Math.Max(triangle[i - 1][j], triangle[i - 1][j - 1]);
}
// Finding maximum from last row
if (i == triangle.Length - 1)
{
if (triangle[i][j] > max) max = triangle[i][j];
}
}
}
return max;
}

static void Main(string[] args)
{
System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch();
sw.Start();
long ans = MaxTotalFromTriangle();
sw.Stop();
Console.WriteLine("Total sum {0}; cpu time {1} ms", ans, sw.Elapsed.TotalMilliseconds);
Console.ReadLine();
}
}



Also this solution is universal and can calculate bigger triangles (lets say 100 rows triangles and more) with a reasonable time.

Have fun!!

Wednesday, April 9, 2008

Project Euler Problem 21

Ok, this is it- problem description:
Let d(n) be defined as the sum of proper divisors of n (numbers less than n which divide evenly into n).
If d(a) = b and d(b) = a, where a ≠ b, then a and b are an amicable pair and each of a and b are called amicable numbers.

For example, the proper divisors of 220 are 1, 2, 4, 5, 10, 11, 20, 22, 44, 55 and 110; therefore d(220) = 284. The proper divisors of 284 are 1, 2, 4, 71 and 142; so d(284) = 220.

Evaluate the sum of all the amicable numbers under 10000.

Fun part,- C# solution:



using System;

class Program
{
static long DivSum(long n)
{
long un = 1 + (long)Math.Sqrt(n);
long sum = 1;

for (int i = 2; i <= un; i++)
{
if (n % i == 0)
sum += i + (n / i);
} return sum;
}

static long AmicablePair(long n)
{
long ds = DivSum(n);
long os = DivSum(ds);
long ret;
ret = (os == n && n != ds) ? ds : 0;
return ret;
}

static long AmicableNumbersSum(long until)
{
bool[] amic = new bool[until];
long nextAmic;
long sum = 0;
for (int i = 1; i < until; i++)
{
if (!amic[i])
{
nextAmic = AmicablePair(i);
if (nextAmic > 0)
{
if (nextAmic < until) amic[nextAmic] = true;
sum += i + nextAmic;
}
}
}
return sum;
}

static void Main(string[] args)
{
System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch();
sw.Start();
Console.WriteLine("Amicable numbers sum {0}", AmicableNumbersSum(10000));
sw.Stop();
Console.WriteLine("Calculation time {0} ms", sw.ElapsedMilliseconds);
Console.ReadLine();
}
}



So solution is fast enough,- runs in about 20 ms. A little bit explanations of why this
solution is fast:

  • We need not to search all divisors, but rather until SQRT limit of number.
    This is because if number, lets say x, have divisor d1 above SQRT(x), then,
    that number x MUST have also a divisor d2 = x/d1, below SQRT(x) limit.
  • We need not to verify are all test numbers amicable, because if we find one
    amicable number, lets say x, then we know that there exists second amicable
    number y, which is sum of proper divisors of x. This rule is true, because
    we acctually searching amicable number pairs. So by applying this rule
    we can use memorization - mark second amicable number in bool array,
    and skip that number from test later.

Have fun!

Sunday, March 23, 2008

Projecteuler.net problem 14 - Csharp vs Fsharp

Problem 14 stated as following:

The following iterative sequence is defined for the set of positive integers:
n → n/2 (n is even)
n → 3n + 1 (n is odd)
Using the rule above and starting with 13, we generate the following sequence:
13 → 40 → 20 → 10 → 5 → 16 → 8 → 4 → 2 → 1
It can be seen that this sequence (starting at 13 and finishing at 1) contains 10 terms. Although it has not been proved yet (Collatz Problem), it is thought that all starting numbers finish at 1.

Which starting number, under one million, produces the longest chain?

Ok, so here it is. Lets have fun with F#:




let seq_length n = n > Seq.unfold(fun x -> match x with
x when x = 0L -> None
x when x = 1L -> Some((1L,0L))
x when x%2L=0L -> Some((x,x/2L))
_ -> Some((x,3L*x + 1L))
) > Seq.length
[for i in 1L..1000000L -> (i, seq_length i)]
> List.fold1_left(fun (ai,al) (i,l) -> if l > al then (i,l) else (ai,al) )
> (fun (x,y) -> x ) > print_any

Elegant, slow version of solution-> ~20 lines; ~10 seconds.

Now lets try C#:




static long CountTerms(long x, ref long[] terms)
{
long tc = 1;
long k = x;

while (k!=1)
{
tc++;
k = (k % 2 == 0) ? k / 2 : 3 * k + 1;
if (k <= x)
{
if (terms[k] > 0)
{
terms[x] = terms[k] + tc;
return terms[x];
}
}
}
terms[x] = tc;
return terms[x] ;
}

static long LongestSeq()
{
long max = 1000000;
long[] termcount = new long[max+1];
long tcmax = 0;
long c = 0;
long ix = 0;

for (int i = 1; i <= max; i++)
{
c = CountTerms(i, ref termcount);
if ( c > tcmax)
{
tcmax = c;
ix = i;
}
}
return ix;
}

static void Main(string[] args)
{
Console.WriteLine(LongestSeq());
Console.ReadLine();
}



Ugly, fast version of solution -> ~0.18 second; ~60 lines.

Conclusion:

F# solutions can be done in more elegant/less code way, but in contrary it is slower that C#
counterpart. Also C# solutions can be written in fast way, but needs more code and attention.
F# solution was by factor 55 slower than C# solution, but also code was by factor of 3 shorter than C# counterpart. So if you want elegancy/generics power - take F#; and if you want only perfomance- take C#.
P.S. There is a way in F# to code (using mutable keyword) very very ugly like in C#, but i`ve not
suggest using it on F#,- why to code ugly in F# if you can grab C# instead ?? :-)

Tuesday, March 18, 2008

Reversing number in Csharp

This time i want to compare 2 methods for reversing a number. One method is based on Array.Reverse method and second - on number modulus 10 operator. Code is this for reversing number as string (C#):



static long ReverseString(long number)
{
char[] str = number.ToString().ToCharArray();
Array.Reverse(str);
return long.Parse(new string(str));
}



And code for reversing number by math operations (C#):



static long ReverseNumber(long number)
{
long digit = number % 10;
long rest = (long) number / 10;
long o = digit;

while (rest > 0)
{
digit = rest % 10;
rest /= 10;
o = (o*10)+ digit;
}
return o;
}



And these separate functions are run for 1000000 iterations.
Then execution time is measured.
Situation is this:

Function ReverseString --> 984 ms.
Function ReverseNumber --> 765 ms.

Conclusion:
Reversing by math operations is a bit faster- faster by factor of 1.3
So if for some reason you need to reverse a number - do it with math,
rather than by string/array operations.