Herhangi Matematiksel Şeyler: Euler Projesi 280. Soru

Karınca ve tohumlar

Çalışkan bir karınca 5x5 bir ızgarada rastgele yürüyor. Yürüyüş merkezdeki kareden başlıyor. Her adımda karınca ızgara dışına çıkmadan rastgele komşu bir kareye hareket ediyor; dolayısıyla karıncanın konumuna bağlı olarak her adımda 2, 3 veya 4 olasılık mevcuttur.

Yürüyüşe başlarken alttaki satırın her bir karesine bir tohum yerleştirilir. Karınca bir tohum taşımıyorken tohum bulunan alt satırdaki bir kareye ulaştığında, tohumu taşımaya başlar. Sonrasında üst satırda bulunan ulaştığı ilk boş kareye tohumu bırakır.

Tüm tohumların sonuçta üst satırda bırakılmış olmasına kadar geçen adım sayısının beklenen değeri kaçtır? Cevabınızı 6 ondalık basamağa yuvarlayarak veriniz.
Cevap: 430.088247

Python:
import collections
import itertools
import math

n = 5
k = 2 ** n - 1
delta = ((1, 0), (-1, 0), (0, -1), (0, 1))
eps = 1e-11

c = collections.defaultdict(float)
c[(0, n // 2, n // 2, 0, k)] = 1

probs = [[0] * n for i in range(n)]

for x in range(n):
    for y in range(n):
        probs[x][y] = 1 / sum(0 <= x + dx < n and 0 <= y + dy < n for dx, dy in delta)

expected = ended = i = 0

v = []
while ended < 1 - eps:
    for x in range(n):
        for y in range(n):
            c[(0, x, y, k, 0)] = 0

cc = collections.defaultdict(float)
    i += 1
    for a in c:
        carrying, x, y, upper, lower = a
        for dx, dy in delta:
            xx, yy = x + dx, y + dy
            if 0 <= xx < n and 0 <= yy < n:
                if carrying and xx == 0 and upper & 1 << yy == 0:
                    cc[(0, xx, yy, upper | 1 << yy, lower)] += c[a] * probs[x][y]
                elif xx == n - 1 and carrying == 0 and lower & 1 << yy:
                    cc[(1, xx, yy, upper, lower & (~(1 << yy)))] += c[a] * probs[x][y]
                else:
                    cc[(carrying, xx, yy, upper, lower)] += c[a] * probs[x][y]
    t = 0
    for x in range(n):
        for y in range(n):
            t += cc[(0, x, y, k, 0)]
    v.append(i * t)
    ended += t
    c = cc

print(round(math.fsum(v), 6))

C#:
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Diagnostics;
using System.Drawing;
using System.IO;

using MathUtilities;

namespace MathProblems
{
    public class Problem280
    {
        /// <summary>
        /// Problem280 - Ant and seeds
        /// </summary>
        public void Run()
        {
            double expectedValue = GetExpectedValue();
            expectedValue = Math.Round(expectedValue, 6);
            Console.WriteLine("The expected number of steps is {0}", expectedValue);

Debug.Assert(expectedValue == 430.088247);
        }

static int[] numAdjacencies = { 2, 3, 3, 3, 2, 3, 4, 4, 4, 3, 3, 4, 4, 4, 3, 3, 4, 4, 4, 3, 2, 3, 3, 3, 2 };

double GetExpectedValue()
        {
            // State mask:  [5 bits, top row seeds] [5 bits, bottom row seeds] [5 bits, ant position as 5*row+col] [1 bit, 1 if ant carrying seed, 0 otherwise]
            // Note that not all possible state mask combinations are valid.
            // Build the list of valid states.
            List<int> validStates = new List<int>();
            for (int stateMask = 0; stateMask < 0xf800; stateMask++)
            {
                int topMask = (stateMask >> 11) & 0x1f;
                int botMask = (stateMask >> 6) & 0x1f;
                int realPos = (stateMask >> 1) & 0x1f;
                int currentRow = realPos / 5;
                int currentCol = realPos % 5;
                bool carrying = (stateMask & 1) != 0;
                int bitCount = Utilities.CountBits(stateMask & 0xffc0);
                if (realPos >= 25)
                    continue;
                if ((bitCount == 5 && !carrying) || (bitCount == 4 && carrying))
                {
                    if (currentRow == 0 && carrying && ((topMask & (0x10 >> currentCol)) == 0))
                        continue;
                    if (currentRow == 4 && !carrying && ((botMask & (0x10 >> currentCol)) != 0))
                        continue;
                    validStates.Add(stateMask);
                }
            }
            // Add the final absorbing state.
            validStates.Add(0xf800);
            // Note the row number for the start state.
            int startState = validStates.IndexOf(0x7d8);

int numTransientStates = validStates.Count - 1;
            double[,] IMinusQ = new double[numTransientStates, numTransientStates];
            for (int i = 0; i < numTransientStates; i++)
                IMinusQ[i, i] = 1;

for (int currentState = 0; currentState < numTransientStates; currentState++)
            {
                int stateMask = validStates[currentState];
                int realPos = (stateMask >> 1) & 0x1f;
                int currentRow = realPos / 5;
                int currentCol = realPos % 5;
                bool carrying = (stateMask & 1) != 0;

// Try up, right, down, left
                if (currentRow > 0)
                    AddTransition(IMinusQ, validStates, currentState, currentRow - 1, currentCol);

if (currentCol < 4)
                    AddTransition(IMinusQ, validStates, currentState, currentRow, currentCol + 1);

if (currentRow < 4)
                    AddTransition(IMinusQ, validStates, currentState, currentRow + 1, currentCol);

if (currentCol > 0)
                    AddTransition(IMinusQ, validStates, currentState, currentRow, currentCol - 1);
            }

double[,] N = MatrixUtilities.Invert(IMinusQ);

double val = 0;
            for (int i = 0; i < numTransientStates; i++)
                val += N[startState, i];

return val;
        }

void AddTransition(double[,] T, List<int> validStates, int currentState, int nextRow, int nextCol)
        {
            int stateMask = validStates[currentState];
            int topMask = (stateMask >> 11) & 0x1f;
            int botMask = (stateMask >> 6) & 0x1f;
            int currentPos = (stateMask >> 1) & 0x1f;
            bool carrying = (stateMask & 1) != 0;

int nextPos = (nextRow * 5 + nextCol) << 1;

int nextStateMask = (stateMask & 0xffc1) | nextPos;

if (nextRow == 0 && carrying && ((topMask & (0x10 >> nextCol)) == 0))
            {
                int nextTopMask = topMask | (0x10 >> nextCol);
                if (nextTopMask == 0x1f)
                    nextStateMask = validStates[validStates.Count - 1];
                else
                    nextStateMask = (nextTopMask << 11) | (botMask << 6) | nextPos;
            }
            else if (nextRow == 4 && !carrying && ((botMask & (0x10 >> nextCol)) != 0))
            {
                nextStateMask = (topMask << 11) | ((botMask & ~(0x10 >> nextCol)) << 6) | nextPos | 1;
            }

int nextState = validStates.IndexOf(nextStateMask);
            Debug.Assert(nextState >= 0);
            if (nextState < validStates.Count - 1)
            {
                T[currentState, nextState] += -1.0 / numAdjacencies[currentPos];
            }
        }
    }
}

Java:
import java.util.ArrayList;
import java.util.Comparator;
import java.util.Enumeration;
import java.util.Hashtable;
import java.util.Iterator;
import java.util.TreeSet;

public class AntSeeds55 {
 public static final Comparator<ArrayList<String>> MY_ORDER = new Comparator<ArrayList<string>>() {
  public int compare(ArrayList<string> list1, ArrayList<string> list2) {
   return Integer.valueOf(list1.get(0)).compareTo(Integer.valueOf(list2.get(0)));
  }
 };
 private Hashtable<String, ArrayList<string>> encodedMap = new Hashtable<String, ArrayList<string>>();
 private Hashtable<String, Integer> codeToInd = new Hashtable<String, Integer>();
 private TreeSet<ArrayList<string>> indToCode = new TreeSet<ArrayList<string>>(MY_ORDER);
 private ArrayList<ArrayList<integer>> map = new ArrayList<ArrayList<integer>>();
 private int highestUnusedInd = 0, finishedInd, garbageInd;
 public String getAnswer() {
  double av = 0;
  double[] probs;
  boolean[] initSeeds = {true, true, true, true, true, false, false, false, false, false};
  String initState = encode(2, 2, initSeeds);
  getStates(initState);
  getMap();
  probs = new double[map.size()];
  probs[codeToInd.get(initState)] = 1;
  int lim = 2500; // needs to be sufficiently high
  for (int i = 1; i <= lim; i++) {
   double[] nextProbs = new double[map.size()];
   for (int j = 0; j < map.size(); j++)
    if (probs[j] != 0) {
     double p = 1 / (double)map.get(j).size();
     for (int k = 0; k < map.get(j).size(); k++)
      nextProbs[map.get(j).get(k)] += p * probs[j];
    }
   probs = nextProbs;
   av += i * probs[finishedInd];
  }
  return Double.toString(av);
 }
 public void getMap() {
  Enumeration<String> e = encodedMap.keys();
  while (e.hasMoreElements()) {
   String currState = e.nextElement();
   ArrayList<string> treeEntry = new ArrayList<string>();
   codeToInd.put(currState, highestUnusedInd);
   treeEntry.add(String.valueOf(highestUnusedInd));
   treeEntry.add(currState);
   indToCode.add(treeEntry);
   highestUnusedInd++; 
  }
  finishedInd = highestUnusedInd++;
  garbageInd = highestUnusedInd++;
  Iterator<ArrayList<string>> iter = indToCode.iterator();
  while (iter.hasNext()) {
   ArrayList<string> currTreeEntry = iter.next();
   String currState = currTreeEntry.get(1);
   ArrayList<integer> mapsTo = new ArrayList<integer>();
   ArrayList<string> encodedMapsTo = encodedMap.get(currState);
   for (int i = 0; i < encodedMapsTo.size(); i++) {
    String nextState = encodedMapsTo.get(i);
    if (isFinished(nextState)) mapsTo.add(finishedInd);
    else mapsTo.add(codeToInd.get(nextState));
   }
   map.add(mapsTo);
  }
  ArrayList<Integer> finalMapsTo = new ArrayList<integer>();
  finalMapsTo.add(garbageInd);
  map.add(finalMapsTo); // finished states map to garbage state
  map.add(new ArrayList<integer>(finalMapsTo)); // garbage state maps to itself
 }
 public void getStates(String state) {
  if (encodedMap.containsKey(state) || isFinished(state)) return;
  ArrayList<string> mapsTo = new ArrayList<string>();
  boolean carrying = extractCarrying(state);
  int x = extractX(state);
  int y = extractY(state);
  boolean[] seeds = extractSeeds(state);
  ArrayList<integer> options = new ArrayList<integer>(); // available moves
  for (int i = 0; i < 4; i++)
   options.add(i);
  if (y == 4) options.remove(options.indexOf(0));
  if (x == 4) options.remove(options.indexOf(1));
  if (y == 0) options.remove(options.indexOf(2));
  if (x == 0) options.remove(options.indexOf(3));
  for (int i = 0; i < options.size(); i++) {
   String nextState;
   int nextX = x;
   int nextY = y;
   boolean[] nextSeeds = seeds.clone();
   switch(options.get(i)) {
    case 0: nextY++; break; // up
    case 1: nextX++; break; // right
    case 2: nextY--; break; // down
    case 3: nextX--; break; // left
   }
   if (carrying && nextY == 4 && !nextSeeds[5 + nextX])
    nextSeeds[5 + nextX] = true; // drop seed on top row
   if (!carrying && nextY == 0 && nextSeeds[nextX])
    nextSeeds[nextX] = false; // pick up seed from bottom row
   nextState = encode(nextX, nextY, nextSeeds);
   mapsTo.add(nextState);
  }
  encodedMap.put(state, mapsTo);
  for (int i = 0; i < mapsTo.size(); i++)
   getStates(mapsTo.get(i));
 }
 public boolean[] extractSeeds(String state) {
  boolean[] seeds = new boolean[10];
  String binary = addLeadingZeroes(Integer.toString(Integer.parseInt(state.substring(2)), 2), 10);
  for (int i = 0; i < 10; i++)
   if (binary.charAt(i) == '1') seeds[i] = true;
  return seeds;
 }
 public int extractY(String state) {
  return Integer.parseInt(state.substring(1, 2));
 }
 public int extractX(String state) {
  return Integer.parseInt(state.substring(0, 1));
 }
 public boolean extractCarrying(String state) {
  int count = 0;
  String binary = Integer.toString(Integer.parseInt(state.substring(2)), 2);
  for (int i = 0; i < binary.length(); i++)
   if (binary.charAt(i) == '1') count++;
  return (count == 4);
 }
 public boolean isFinished(String state) {
  return state.substring(2).equals("0031");
 }
 public String encode(int x, int y, boolean[] seeds) {
  String result = String.valueOf(x) + String.valueOf(y);
  String binary = "";
  for (int i = 0; i < 10; i++)
   binary += (seeds[i]) ? "1" : "0";
  return result + addLeadingZeroes(String.valueOf(Integer.parseInt(binary, 2)), 4);
 }
 public String addLeadingZeroes(String s, int x) {
  if (x <= s.length()) return s;
  return addLeadingZeroes("0" + s, x);
 }
}

Matlab:
function total = solve()    
    global fv;
    global fp;    
    fp = zeros(5, 31, 5); fv = zeros(5, 31, 5);    
    bvec = dec2base(1:31, 2)-double('0');
    [ip, iv]    = collect_prob([3,3], [0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 1 1 1 1 1]);        
    for i = 1:5
        for j = 1:31
            [p,v] = collect_prob([1,i], [repmat(zeros(1,5), 4, 1); bvec(j,:)]);
            fp(i,j,1:length(p)) = transpose(p);
            fv(i,j,1:length(v)) = transpose(v);
        end
    end           
    total = 2*ip(1)*(iv(1)+f(1, [0 1 1 1 1], ones(1,5)));    
    total = total + 2*ip(2)*(iv(2)+f(2, [1 0 1 1 1], ones(1,5)));
    total = total + ip(3)*(iv(3)+f(3, [1 1 0 1 1], ones(1,5)));    
end

function z = f(x, a, b)    
    global fv;    
    global fp;    
    bmul = 2.^(transpose(4:-1:0));
    if(sum(b) == 0)
        z = 0;
    else
        z = 0;              
        bv = b*bmul;
        p = fp(x, bv,:);
        v = fv(x, bv,:);
        k=0;
        for i = 1:5
            if(b(i)==1)
                b(i) = 0;   
                k = k + 1;
                z = z + p(k)*(v(k)+f(i, b, a));
                b(i) = 1;
            end
        end
    end    
end

function [p, E] = collect_prob(x, dest, tol)
    B = zeros(7);    
    B(2:6,2:6) = 1/3;
    B(2,2) = 1/2;
    B(2,6) = 1/2;
    B(6,2) = 1/2;
    B(6,6) = 1/2;
    B(3:5,3:5) = 1/4;    
    B(2:6,2:6) = B(2:6,2:6).*(1-dest);    
    X = zeros(7);
    X(x(1)+1,x(2)+1) =1;
    d = zeros(7);
    d(2:6,2:6) = dest;
    d = (d == 1);
    p = zeros(sum(sum(dest)), 1);
    E = zeros(size(p));
    step = 0;
    while(1) %iterate
        step = step + 1;
        X(2:6,2:6) = X(3:7, 2:6).*B(3:7, 2:6)+X(1:5,2:6).*B(1:5,2:6)+X(2:6,1:5).*B(2:6,1:5)+X(2:6,3:7).*B(2:6,3:7);
        yterms = X(d);                
        if(norm(yterms) ~= 0)
            p = p + yterms;
            E = E + yterms*step;
            if(max(yterms) < tol)
                break;
            end       
        end
    end
    E = E./p;
end

Delphi:
program mathC280;
uses sysutils, windows;

const tries = 2000;

var pxy, pxyold: array[0..6, 0..6] of extended;
var neighbours: array[0..6, 0..6] of extended;

type TRow = array[1..5] of extended;
var exp33, prob33: TRow;

var exps, probs: array[1..5, 0..31, 1..5] of extended; // exps[start, state, stop]

var ones: array[1..5] of integer = (1, 2, 4, 8, 16);      // i.e. ones[2]  = 00010
var zeros: array[1..5] of integer = (30, 29, 27, 23, 15); // i. e. zeros[2] = 11101
(**********************************************)
procedure fillNeighbours;
var x, y: integer;
// the number of neighbours beyond the border: 1 (because we divide)
begin
for x:= 0 to 6 do for y:= 0 to 6 do neighbours[x, y]:= 1;
for x:= 1 to 5 do for y:= 1 to 5 do neighbours[x, y]:= 4;
x:= 1; for y:= 1 to 5 do neighbours[x, y]:= 3;
x:= 5; for y:= 1 to 5 do neighbours[x, y]:= 3;
y:= 1; for x:= 1 to 5 do neighbours[x, y]:= 3;
y:= 5; for x:= 1 to 5 do neighbours[x, y]:= 3;
neighbours[1, 1]:= 2; neighbours[1, 5]:= 2;
neighbours[5, 1]:= 2; neighbours[5, 5]:= 2;
end;
(**********************************************)
procedure build33(var targetExp, targetProb: TRow);
var x, y, n, i: integer; p: extended;
begin // probability and expectation for startboard
fillNeighbours;

for i:= 1 to 5 do neighbours[i, 1]:= 1E100; // first row is target (all 5 are stops)
for x:= 0 to 6 do for y:= 0 to 6 do pxyold[x, y]:= 0;
pxyold[3, 3]:= 1; // starts in (3,3)
for i:= 1 to 5 do
  begin targetExp[i]:= 0; targetProb[i]:= 0;
  end;
for n:= 1 to tries do
  begin
  for x:= 1 to 5 do for y:= 1 to 5 do
  pxy[x, y]:= pxyold[x-1, y]/neighbours[x-1, y]+pxyold[x+1, y]/neighbours[x+1, y]
             +pxyold[x, y-1]/neighbours[x, y-1]+pxyold[x, y+1]/neighbours[x, y+1];
  for x:= 1 to 5 do for y:= 1 to 5 do pxyold[x, y]:= pxy[x, y];
  for i:= 1 to 5 do
    begin targetExp[i]:= targetExp[i]+n*pxy[i, 1];
          targetProb[i]:= targetProb[i]+pxy[i, 1];
    end;
  end;
for i:= 1 to 5 do targetExp[i]:= targetExp[i]/targetProb[i];
p:= 0; // a quick check
for i:= 1 to 5 do p:= p+targetProb[i];
if abs(p-1)>1E-10 then writeln('error p');
end;
(**********************************************)
procedure buildExpsAndProbs(startx: integer; targetState: integer);
// targetState: position of the stops i.e. seeds in row 1 or empty in row 5
var x, y, n, i: integer; var targetExp, targetProb: TRow; p: extended;
begin fillNeighbours;
for i:= 1 to 5 do targetProb[i]:= 0;
for i:= 1 to 5 do targetExp[i]:= 0;

for i:= 1 to 5 do if targetState and ones[i]<>0 then neighbours[i, 5]:= 1E100;
for x:= 0 to 6 do for y:= 0 to 6 do pxyold[x, y]:= 0;
pxyold[startx, 1]:= 1; // starting in (startx, 5) gives the same

for n:= 1 to tries do
  begin
  for x:= 1 to 5 do for y:= 1 to 5 do
  pxy[x, y]:= pxyold[x-1, y]/neighbours[x-1, y]+pxyold[x+1, y]/neighbours[x+1, y]
             +pxyold[x, y-1]/neighbours[x, y-1]+pxyold[x, y+1]/neighbours[x, y+1];
  for x:= 1 to 5 do for y:= 1 to 5 do pxyold[x, y]:= pxy[x, y];
  for i:= 1 to 5 do
  if targetState and ones[i]<>0 then // calculate only for stop fields
    begin targetExp[i]:= targetExp[i]+n*pxy[i, 5];
          targetProb[i]:= targetProb[i]+pxy[i, 5];
    end;
  end;

for i:= 1 to 5 do if targetState and ones[i]<>0
                  then targetExp[i]:= targetExp[i]/targetProb[i];
// storing of values
for i:= 1 to 5 do exps[startx, targetState, i]:= targetExp[i];
for i:= 1 to 5 do probs[startx, targetState, i]:= targetProb[i];
p:= 0; // a quick check
for i:= 1 to 5 do p:= p+targetProb[i];
if abs(p-1)>1E-10 then writeln('error p ',p);
end;
(**********************************************)
var totalProb, totalExp: extended; noTree: integer;
(**********************************************)
procedure doWalk(startx, state1, state5: integer; down: boolean;
               var startExp, startProb: extended; level: integer);
var targetx, theState, newState1, newState5, newStartx: integer;
    newStartProb, newStartExp: extended; newDown: boolean;
begin
if down then theState:= state1 else theState:= state5;
newState1:= state1; newState5:= state5;
for targetx:= 1 to 5 do // x-coordinate of target
if theState and ones[targetx]<>0 then // if stop field
  begin
  newStartProb:= StartProb*probs[startx, theState, targetx];
  newStartExp:= StartExp+exps[startx, theState, targetx];
  if down then newState1:= theState and zeros[targetx]
          else newState5:= theState and zeros[targetx];
  if newState5=0 then
    begin
    totalProb:= totalProb+newStartProb;
    totalExp:= totalExp+newStartExp*newStartProb;                                            
    inc(noTree); // debugging
    exit;
    end;
  newStartx:= targetx; newDown:= not down;
  doWalk(newStartx, newState1, newState5, newDown, newStartExp, newStartProb, level+1);
  end;
end;

var startx, targetx, state1, newState1, state5, targetState: integer;
    startProb, startExp: extended; s: string;
begin
build33(exp33, prob33);
for startx:= 1 to 5 do for targetState:= 1 to 31 do
  begin
  buildExpsAndProbs(startx, targetState);
  end;

noTree:= 0; totalProb:= 0; totalExp:= 0;
state1:= 1+2+4+8+16; state5:= 1+2+4+8+16;
for targetx:= 1 to 5 do
  begin
  startProb:= prob33[targetx]; // first move
  startExp:=  exp33[targetx];
  newState1:= state1 and zeros[targetx];
  doWalk(targetx, newState1, state5, false, startExp, startProb, 1);
  end (* for i1 *);
str(totalExp:6:6,s);
writeln('totalExp ', s);
end.

Rust:
const N: usize = 1 << 16;

fn iterate(good: &mut [bool; N], exp_old: &[f64; N], exp_new: &mut [f64; N]) {
 for i in 0..N {
  if !good[i] { continue }
  
  let x = (i & (0b111 << 13)) >> 13;
  let y = (i & (0b111 << 10)) >> 10;
  if x >= 5 || y >= 5 { good[i] = false; continue }

let seeds = i & ((1 << 10) - 1);
  let no_seeds = seeds.count_ones();
  if !(no_seeds == 4 || no_seeds == 5) { good[i] = false; continue }

let mut seeds2 = seeds;

if y == 0 && no_seeds == 5 && (seeds & (1 << x) != 0) { seeds2 = seeds & !(1 << x) }
  if y == 4 && no_seeds == 4 && (seeds & (1 << x+5) == 0) { seeds2 = seeds | (1 << x+5) }

if seeds2 == 0b1111100000 { exp_new[i] = 0.0; continue }

let mut count = 0;
  let mut sum = 0.0;
  if x > 0 { count += 1; sum += exp_old[x-1 << 13 | y << 10 | seeds2] }
  if x < 4 { count += 1; sum += exp_old[x+1 << 13 | y << 10 | seeds2] }
  if y > 0 { count += 1; sum += exp_old[x << 13 | y-1 << 10 | seeds2] }
  if y < 4 { count += 1; sum += exp_old[x << 13 | y+1 << 10 | seeds2] }
  exp_new[i] = 1.0 + sum / (count as f64);
 }
}

fn main() {
 let exp1 = &mut [0.0; N];
 let exp2 = &mut [0.0; N];
 let good = &mut [true; N];

for j in 0..1_200 {
  iterate(good, exp1, exp2);
  iterate(good, exp2, exp1);
  if j > 1000 { println!("run {}: {:.6}", j, (*exp1)[0b0100100000011111]); }
 }
}

Herhangi Matematiksel Şeyler

Sayfalar

18 Mart 2020

Euler Projesi 280. Soru