## How do I expand my simplex in Expanding Polytope Algorithm (EPA) when the origin lies on it?

I’m doing GJK + EPA in 2D for collision detection + resolution.

For GJK, I was using the vector (1,0) as my first support direction. However, if the origin happens to fall exactly on the simplex line ( quite easily happens with circles placed on an axis-aligned plane), inside EPA, it doesn’t know which direction to expand to. That is because the expansion is supposed to happen in the direction of origin. But because the origin is on the line, the vector towards origin has a length of 0. So although collision is detected in GJK, I don’t get any contact normal + penetration depth from EPA.

I’ve ‘solved’ it by randomizing the initial support direction every time. It can still miss a few times in a second but the visual result looks fine.

Is there a standard way of handling this edge case that I don’t know of? Or maybe this shouldn’t be happening at all and I’m doing something wrong?

## Coupled differential equations metropolis algorithm

I would like to simulate the following coupled SDE. $$d\Omega=-\Omega H\left(\mathrm{n}\right)d\tau$$ $$dn=-Ad\tau+\sum_{i}B^{\left(i\right)}dW_{i}+\sum_{i}C^{\left(i\right)}dW’_{i}$$ provided $$A=\frac{1}{2}n\left(T-UM\right)\left(I-n\right)+\frac{1}{2}\left(I-n\right)\left(T-UM\right)n$$

$$B_{xy}^{\left(i\right)}=\sqrt{\frac{U}{2}}\sum_{\sigma,\sigma’}\sigma_{\sigma\sigma’}^{z}n_{x,\left(i\sigma’\right)}\left(\delta_{\left(i\sigma\right),y}-n_{\left(i\sigma\right),y}\right)$$ $$\sigma^{z}=\begin{bmatrix}1 & 0\ 0 & -1 \end{bmatrix}$$

$$C_{xy}^{\left(i\right)}=\sqrt{\frac{U}{2}}\sum_{\sigma,\sigma’}\sigma_{\sigma\sigma’}^{z}\left(\delta_{x,\left(i\sigma’\right)}-n_{x,\left(i\sigma’\right)}\right)n_{\left(i\sigma\right),y}$$

$$M_{\left(i\sigma\right),\left(j\sigma’\right)}=\delta_{ij}\sum_{\eta,\eta’}n_{\left(i\eta\right),\left(i\eta’\right)}\left(\sigma_{\sigma\sigma’}^{z}\sigma_{\eta\eta’}^{z}-\sigma_{\sigma\eta’}^{z}\sigma_{\eta\sigma’}^{z}\right)$$

Initial condition as $$n=\frac{1}{2}I$$ $$I$$ is a 2M X 2M matrix.$$M$$ is the mode number.

$$H(n)=\sum_{x,y=1}^{2M}T_{xy}n_{xy}+\frac{U}{2}\sum_{i=1}^{M}\sum_{\begin{array}{cc} \eta, & \eta’\ \sigma, & \sigma’ \end{array}\begin{array}{c} \ =1 \end{array}}^{2}\left[n_{\left(i\eta\right),\left(i\sigma’\right)}n_{\left(i\sigma\right),\left(i\eta’\right)}-n_{\left(i\eta\right),\left(i\eta’\right)}n_{\left(i\sigma\right),\left(i\sigma’\right)}+n_{\left(i\eta\right),\left(i\sigma’\right)}\delta_{\sigma\eta’}-n_{\left(i\eta\right),\left(i\sigma’\right)}\delta_{\left(i\sigma\right),\left(i\eta’\right)}\right]\sigma_{\eta\eta’}^{z}\sigma_{\sigma\sigma’}^{z}$$

For example if $$M=2$$ initial value of If M=2, then

Initial value of n matrix will be

$$n=.5\begin{bmatrix}1 & 0 & 0 & 0\ 0 & 1 & 0 & 0\ 0 & 0 & 1 & 0\ 0 & 0 & 0 & 1 \end{bmatrix}$$

Also $$T=\begin{bmatrix}-.1 & -2 & 0 & 0\ -2 & -.1 & 0 & 0\ 0 & 0 & -.1 & -2\ 0 & 0 & -2 & -.1 \end{bmatrix}$$

$$U/t=2$$

Is there a way to solve the above using the successive metropolis method as shown below? I found this in the following paper. https://arxiv.org/abs/0704.3792

## About Fermat’s Integer Factoring algorithm

An implementation of Fermat’s algorithm that is a slight improvement on code in this book is here.

FermatFactor[n_?OddQ]:=Module[{s=Floor@Sqrt[n+01],r,u,v=1},   If[s^2==n,Return@{s,s}];   {r,u}={(s+1)^2-n,2*s+3};   While[r!=0,     While[r>0,       r-=v;v+=2     ];     If[r<0,r+=u;u+=2]   ];   (u+{-v,v-2})/2 ] 

The main improvement above is the fast implementation of Floor[Sqrt[n]] from MichaelE2. The above implementation will factor the following in less than 0.1 seconds!

SeedRandom[100]; p=10^50; n=NextPrime@RandomInteger[{p-80000,p-60000}]*NextPrime@RandomInteger[{p+80000,p+60000}]; FermatFactor[n] (*{99999999999999999999999999999999999999999999920467,100000000000000000000000000000000000000000000075667}*) 

I have Mathematica Version 12.0.0 and FactorInteger[n] takes a very long time on that example. Can FactorInteger in Mathematica 12.2 do it in a reasonable amount of time?

## Algorithm for finding a function connecting points in a space

I have a set of points in a space and I need to find an algorithm which will find the function which connects these points.

From current analysis I know that no three points lie in the same plane. So how can I use Mathematica to find such a function in this space?

## Overview

I want to apply a specified Chop function to every step in an function call in Mathematica (in my case LegendreP), especially when encountering a machine underflow error. The resulting function should be applicable to a dataframe.

## The Setup

I am trying to calculate a function including an associated Legendre function with complex indices

function[t_] = 1/Sin[t] * LegendreP[-1/2 + V * I, l+1, Cos[t]]

where t is between 0 and Pi, V is of order of magnitude 10 and l is between 10 and 100. Beside the function I need its logarithmic derivative function'/function.

Ideally I want to do this by applying the function to a dataframe and appending the result in a separate column.

Append[#, "function" -> function[#timeframe]]

where timeframe is a column with all the t values.

## The Problem

When I run this code for any l bigger than 12 and very small t~1e-5 values, the LegendreP Algorithm throws a machine underflow error General::munfl because it cannot execute a multiplication of extremly small complex numbers.

While for a single call or a plot, it seems to do some chopping or return Indeterminate, when I write it to a dataframe with

Append[#, "function" -> function[#timeframe]]

it just returns Failure and does not write anything to the dataframe.

## What I have tried so far

I have tried to use Chop and Threshold, but this does not seem to apply to the single steps of the algorithm but only the final result.

The way I "solve" the problem at the moment is to catch the error and return 0 instead of my function. This is not ideal since the real or imaginary part of the step in question and the result might not be negligible while the other one is, or it might diverge instead of converge to 0.

Since the multiplication that raises the error lists numbers ~1e-300 or so, I doubt that the problem is solvable by increasing the precision.

## My Goal

Ideally I’d like to call Chop, whenever Mathematica encounters a machine underflow. The behavior of Chop on complex numbers is exactly what I need. This way I should be able to preserve the real or imaginary part that does not vanish.

Is the error handling different, when applied to a dataframe as it relates to this question (for plots or even single evaluation points I don’t have the same issue) or can an indeterminate/NaN be written to a dataframe?

Is there a way to set a "global chop rule"?

Grateful for any hint 😀

## A* algorithm: need help fixing path that come in contact with obstacle

I am using A* as a pathfinding technique for my AI; it works fine until it gets close to an obstacle (in my case rocks). Right now it just continues the loop if it finds a neighbor thats a rock. It should not actually do this. It should check it.

Allowed Moves: FORWARD,LEFT,RIGHT (LEFT & RIGHT at diagonals) turns are basically two phases: FORWARD, turn face then FORWARD (counts as one move with no additional cost)

The AI should know to turn left or right on the rock in the direction of the goal while also taking other rocks into account.

The line that checks for rocks is: if(at == 1 || at == 2) continue;

I guess you could use the neighborlist to check the sides of the ship.

However it shouldn’t always check it. Only when it comes in contact with a rock

Scenario 1: (ship should turn left (one move) once then continue on path)

Scenario 2: (ship should either turn left or right twice (two moves) to unblock itself)

Scenario 3: (ship should turn left or right depending on which path is shorter: doing two lefts will hit rock twice but distance is shorter than if it went right by 1 tile)

In each of these scenarios the face of the ship is the only thing that changes; unless its a forward move into the rock then there is no change. If right/left were used in any other situation (regular tiles) it would change position also.

public class AStarSearch {          private ServerContext context;     private List<AStarNode> openList;     private List<AStarNode> closedList;          private double ORTHOGONAL_COST = 1.0;     private double DIAGONAL_COST = ORTHOGONAL_COST * Math.sqrt(2.0);          public AStarSearch(ServerContext context) {         this.context = context;     }      private Comparator<AStarNode> nodeSorter = new Comparator<AStarNode>() {          @Override         public int compare(AStarNode n0, AStarNode n1) {             if(n1.fCost < n0.fCost) return 1;             if(n1.fCost > n0.fCost) return -1;             return 0;         }              };      public List<AStarNode> findPath(Player bot, Position goal){         openList = new ArrayList<AStarNode>();         closedList = new ArrayList<AStarNode>();         List<AStarNode> neighbors = new ArrayList<AStarNode>();         AStarNode current = new AStarNode(bot, bot.getFace(), MoveType.NONE, null, 0, bot.distance(goal));         openList.add(current);          while(openList.size() > 0) {             Collections.sort(openList, nodeSorter);             current = openList.get(0);             if(current.position.equals(goal)) {                 List<AStarNode> path = new ArrayList<AStarNode>();                 while(current.parent != null) {                     path.add(current);                     current = current.parent;                 }                 openList.clear();                 closedList.clear();                 Collections.reverse(path);                 return path;             }             openList.remove(current);             closedList.add(current);                          int x = current.position.getX();             int y = current.position.getY();             switch (current.face) {                 case NORTH:                     neighbors.add(new AStarNode(new Position(x, y), VesselFace.NORTH, MoveType.NONE,current,0,0));                     neighbors.add(new AStarNode(new Position(x, y+1), VesselFace.NORTH, MoveType.FORWARD,current,0,0));                     neighbors.add(new AStarNode(new Position(x-1, y+1), VesselFace.WEST, MoveType.LEFT,current,0,0));                     neighbors.add(new AStarNode(new Position(x+1, y+1), VesselFace.EAST, MoveType.RIGHT,current,0,0));                     break;                 case EAST:                     neighbors.add(new AStarNode(new Position(x, y), VesselFace.EAST, MoveType.NONE,current,0,0));                     neighbors.add(new AStarNode(new Position(x+1, y), VesselFace.EAST, MoveType.FORWARD,current,0,0));                     neighbors.add(new AStarNode(new Position(x+1, y+1), VesselFace.NORTH, MoveType.LEFT,current,0,0));                     neighbors.add(new AStarNode(new Position(x+1, y-1), VesselFace.SOUTH, MoveType.RIGHT,current,0,0));                     break;                 case SOUTH:                     neighbors.add(new AStarNode(new Position(x, y), VesselFace.SOUTH, MoveType.NONE,current,0,0));                     neighbors.add(new AStarNode(new Position(x, y-1), VesselFace.SOUTH, MoveType.FORWARD,current,0,0));                     neighbors.add(new AStarNode(new Position(x-1, y-1), VesselFace.WEST, MoveType.RIGHT,current,0,0));                     neighbors.add(new AStarNode(new Position(x+1, y-1), VesselFace.EAST, MoveType.LEFT,current,0,0));                     break;                 case WEST:                     neighbors.add(new AStarNode(new Position(x, y), VesselFace.WEST, MoveType.NONE,current,0,0));                     neighbors.add(new AStarNode(new Position(x-1, y), VesselFace.WEST, MoveType.FORWARD,current,0,0));                     neighbors.add(new AStarNode(new Position(x-1, y+1), VesselFace.NORTH, MoveType.RIGHT,current,0,0));                     neighbors.add(new AStarNode(new Position(x-1, y-1), VesselFace.SOUTH, MoveType.LEFT,current,0,0));                     break;             }             for(AStarNode neighborNode : neighbors) {                 // Compute the cost to get *to* the action tile.                 double costToReach = current.position.distance(neighborNode.position);                  int at = context.getMap().getTile(neighborNode.position.getX(), neighborNode.position.getY());                 if(at == 1 || at == 2) continue; // this is the line where it checks if tile is rock or not                  double gCost = current.gCost + costToReach;                 double hCost = heuristicDistance(neighborNode.position,goal);                 AStarNode node = new AStarNode(neighborNode.position, neighborNode.face,neighborNode.move, current, gCost, hCost);                 if(positionInList(closedList, neighborNode.position) && gCost >= node.gCost) continue;                 if(!positionInList(openList, neighborNode.position) || gCost < node.gCost) openList.add(node);             }         }         closedList.clear();         return null;     }          private double getActionCost(Position node, int currentTile) {         if(currentTile > 3 && currentTile < 11) {             return 0.2;         }else {             return 1;            }     }      private double heuristicDistance(Position current, Position goal) {         int xDifference = Math.abs(goal.getX() - current.getX());         int yDifference = Math.abs(goal.getY() - current.getY());          int diagonal = Math.min(xDifference, yDifference);         int orthogonal = xDifference + yDifference - 2 * diagonal;          return orthogonal * ORTHOGONAL_COST + diagonal * DIAGONAL_COST;     }          private boolean positionInList(List<AStarNode> list, Position position) {         for(AStarNode n : list) {             if(n.position.equals(position)) return true;         }         return false;     }  } 

AStarNode:

public class AStarNode {      public Position position;     public VesselFace face;     public MoveType move;     public AStarNode parent;     public double fCost, gCost, hCost;          public AStarNode(Position position, VesselFace face, MoveType move, AStarNode parent, double gCost, double hCost) {         this.position = position;         this.face = face;         this.move = move;         this.parent = parent;         this.gCost = gCost;         this.hCost = hCost;         this.fCost = this.gCost + this.hCost;     }    }   

There will be no additional cost of running into a rock as long as its a shorter route. Also, if a ship tries to turn left or right from its current position; but there is a rock at that tile it will move up one tile and changes its direction.

The overall question/goal: How do I fix my current code to account for these situations; please provide an implementation or instructions.

## Pathfinding algorithm isn’t finding correct route

I am attempting an online coding challenge wherein I am to implement a pathfinding algorithm that finds the shortest path between two points on a 2D grid. The code that is submitted is tested against a number of test cases that I, unfortunately, am unable to see, but it will however tell me if my answer for shortest distance is correct or not. My implementation of the A* algorithm returns a correct answer on 2/3 test cases and I cannot seem to figure out what scenario might create an incorrect answer on the third?

I have tried several of my own test cases and have gotten correct answers for all of those and at this point am feeling a little bit lost. There must be something small in my code that I am not seeing that is causing this third case to fail.

More details

• The grid is w by h and contains only 1’s (passable) and 0’s (impassable) with every edge having a cost of 1 and the pathway cannot move diagonally It all starts with the FindPath function which is to return the length of the shortest path, or -1 if no path is available
• pOutBuffer is used to contain the path taken from beginning to end (excluding the starting point). If multiple paths are available then any will be accepted. So it isnt looking for one path in particular
• I know the issue is not the result of time or memory inefficiency. I has to be either the distance returned is incorrect, or the values in pOutBuffer are incorrect.

Any help would be greatly appreciated as I am just about out of ideas as to what could possibly be wrong here. Thank you.

#include <set> #include <vector> #include <tuple> #include <queue> #include <unordered_map>  inline int PositionToIndex(const int x, const int y, const int w, const int h) {     return x >= 0 && y >= 0 && x < w  && y < h? x + y * w : -1; }  inline std::pair<int, int> IndexToPosition(const int i, const int w) {     return std::make_pair<int, int>(i % w, i / w); }  inline int Heuristic(const int xa, const int ya, const int xb, const int yb) {     return std::abs(xa - xb) + std::abs(ya - yb); }  class Map { public:     const unsigned char* mapData;     int width, height;      const std::vector<std::pair<int, int>> directions = { {1,0}, {0,1}, {-1,0}, {0,-1} };      Map(const unsigned char* pMap, const int nMapWidth, const int nMapHeight)     {         mapData = pMap;         width = nMapWidth;         height = nMapHeight;     }      inline bool IsWithinBounds(const int x, const int y)      {         return x >= 0 && y >= 0 && x < width && y < height;     }      inline bool IsPassable(const int i)     {         return mapData[i] == char(1);     }       std::vector<int> GetNeighbours(const int i)     {         std::vector<int> ret;          int x, y, neighbourIndex;         std::tie(x, y) = IndexToPosition(i, width);          for (auto pair : directions)         {             neighbourIndex = PositionToIndex(x + pair.first, y + pair.second, width, height);             if (neighbourIndex >= 0 && IsWithinBounds(x + pair.first, y + pair.second) && IsPassable(neighbourIndex))                 ret.push_back(neighbourIndex);         }          return ret;     } };  int FindPath(const int nStartX, const int nStartY,     const int nTargetX, const int nTargetY,     const unsigned char* pMap, const int nMapWidth, const int nMapHeight,     int* pOutBuffer, const int nOutBufferSize) {     int ret = -1;      // create the map     Map map(pMap, nMapWidth, nMapHeight);      // get start and end indecies     int targetIndex = PositionToIndex(nTargetX, nTargetY, nMapWidth, nMapHeight);     int startIndex = PositionToIndex(nStartX, nStartY, nMapWidth, nMapHeight);          // if start and end are same exit     if (targetIndex == startIndex) return 0;          std::unordered_map<int, int> pathway = { {startIndex, startIndex} };     std::unordered_map<int, int> distances = { {startIndex, 0} };      // queue for indecies to process     typedef std::pair<int, int> WeightedLocation;     std::priority_queue<WeightedLocation, std::vector<WeightedLocation>, std::greater<WeightedLocation>> queue;      queue.emplace(0, startIndex);          while (!queue.empty())     {         int currentWeight, currentIndex;         std::tie(currentWeight, currentIndex) = queue.top();         queue.pop();          if (currentIndex == targetIndex)             break;          int newDistance = distances[currentIndex] + 1;         for (int n : map.GetNeighbours(currentIndex))         {             if (distances.find(n) == distances.end() || newDistance < distances[n])             {                 distances[n] = newDistance;                  int weight = newDistance + Heuristic(n % nMapWidth, n / nMapWidth, nTargetX, nTargetY);                 queue.emplace(weight, n);                 pathway[n] = currentIndex;             }         }     }      if (pathway.find(targetIndex) != pathway.end())     {         int current = targetIndex;          while (current != startIndex)         {             int outIndex = distances[current] - 1;             pOutBuffer[distances[current] - 1] = current;             current = pathway[current];         }         ret = distances[targetIndex];     }          return ret; } 

## Function that finds the quotient and remainder of the Division Algorithm applied to two integers a and b

Can someone help guide me on writing the code this problem? I am so lost.

Write a function procedure which finds the quotient and remainder of the Division Algorithm applied to two integers a and b.

## Speeding up the Rummikub algorithm – explanation required

Regarding this question: Rummikub algorithm.

I was reading the first part of the solution in the posted answer (specifically, when there are no jokers involved, all tiles are distinct and only four colours are involved). Then, I reached the part in which the says that the algorithm runs in $$O(ABCD(A+B+C+D))$$ time, which is easy to determine why.

However, he the goes on to saying that we can speed up the algorithm so as to run in $$O(ABCD)$$ time by changing "the recurrence to ensure this occurs only once while maintaining correctness, which leads to $$O(1)$$ time for every ‘cell’ in the DP-table".

My problem is: I do not see how this can be done. I have tried playing around a bit with the recurrence, but I do not see how it can be modified, or what else we should keep track of so that we can speed up the time.

## Efficient algorithm for this combinatorial problem [closed]

$$\newcommand{\argmin}{\mathop{\mathrm{argmin}}\limits}$$

I am working on a combinatorial optimization problem and I need to figure out a way to solve the following equation. It naturally popped up in a method I chose to use in my assignment I was working on.

Given a fixed set $$\Theta$$ with each element $$\in (0,1)$$ and total $$N$$ elements ($$N$$ is about 25), I need to find a permutation of elements in $$\Theta$$ such that $$\vec K = \argmin_{\vec k = Permutation(\Theta)} \sum_{i=1}^N t_i D(\mu_i||k_i)$$ where $$\vec t, \vec \mu$$ are given vectors of length $$N$$ and $$D(p||q)$$ is the KL Divergence of the bernoulli distributions with parameters $$p$$ and $$q$$ respectively. Further, all the $$N$$ elements of $$\vec t$$ sum to 1 and $$\vec \mu$$ has all elements in $$[0,1]$$.

It is just impossible to go through all $$N!$$ permutations. A greedy type of algorithm which does not give exact $$\vec K$$ would also be acceptable to me if there is no other apparent method. Please let me know how to proceed!