# Random Point from a 2D Donut Distribution

A donut distribution (for lack of a better term) is a random distribution of points in a 2-dimensional plane, forming a donut-like shape. The distribution is defined by two parameters: the radius r and spread s, in which the distance to the origin follows a normal (Gaussian) distribution around r, with a standard deviation s. The angular distribution is uniform in the range [0,2π).

## The challenge

Given a radius r and spread s, your code should yield the Cartesian ((x,y)) coordinates of a single point chosen from this distribution.

## Remarks

• Running your code multiple times with the same input should result in the specified distribution.
• Outputting polar coordinates is too trivial and not allowed.
• You can output Cartesian coordinates in any way allowed by the default I/O rules.
• This includes complex values.

## Valid approaches

Several algorithms can be used to yield the desired distribution, including but not limited to

1. Choose a from the uniform distribution [0,2π) and b from the normal distribution (r,s).
Let x = b*cos(a) and y = b*sin(a).
2. Choose a from the uniform distribution [0,4) and b from the normal distribution (r,s).
Let x+y*i = b*i^a.
3. Choose a,b,c all from the normal distribution (0,1).
Let d = a+b*i and x+y*i = d/abs(d) * (c*s+r).

## Example distributions (N=1000)

Below: r=1, s=0.1

Below: r=3, s=1

Below: r=1, s=0

Below: r=100, s=5

• the distance to the origin follows a normal (Gaussian) distribution This is confusing, because any Gaussian distribution can produce negative numbers, and a distance cannot be negative. Your methods 1 and 2 (I haven't looked at what method 3 does) correspond to taking the absolute value of the Gaussian (and shifting the phase by 180 degrees, which is not significant) Mar 7 at 22:36
• @LuisMendo That's a valid point. The third method does basically the same. Consider the distance to be the displacement in the chosen angular direction. Mar 8 at 7:49
• @LuisMendo method 3 seems to be the same as methods 1 and 2, with the random angle generated as the direction of the vector formed by two iid normal variables. Mar 8 at 20:08
• @Nitrodon Thanks! Mar 8 at 23:07
• re "lack of a better term": in math (complex analysis etc.) that's usually called "annulus". Mar 9 at 3:08

# APL(Dyalog Unicode), 38 bytes SBCS

Takes r on the left and s on the right and returns a complex number.

{(⍺+⍵×(.5*⍨¯2×⍟?0)×1○○2×?0)×¯12○2×○?0}


Try it on APLgolf!

Uses the first approach presented in the question. Dyalog doesn't have a builtin for sampling from a normal distribution, so this uses the Box–Muller transform to convert to random numbers from $$\(0,1)\$$ to a normally distributed value:

(⍺+⍵×(.5*⍨¯2×⍟?0)×1○○2×?0) draws a normally distributed value $$\b \sim N(\alpha, \omega^2)\$$:
?0 random number $$\c \in (0,1)\$$
1○○2×?0: $$\\sin(2\pi c)\$$
?0 random number $$\d \in (0,1)\$$
.5*⍨¯2×⍟?0: $$\\sqrt{-2\ln{d}}\$$
⍺+⍵× Scale from $$\N(0, 1)\$$ to $$\N(\alpha, \omega^2)\$$

?0 generates a random number $$\a \in (0,1)\$$
¯12○2×○?0: $$\e^{i 2\pi a} = \sin(2\pi a)i + cos(2\pi a)\$$

The product of these two values is the result.

Plotting code and images:

'InitCauseway' ⎕CY 'sharpplot'
InitCauseway ⍬
sp←⎕NEW Causeway.SharpPlot(700)
sp.SetTrellis(2 2)
sp.TrellisStyle←4

F ← {(⍺+⍵×(.5*⍨¯2×⍟?0)×1○○2×?0)×¯12○2×○?0}

:For r s :In (1 0.1)(3 1)(1 0)(100 5)
sp.NewCell
sp.Heading←'r = ',(⍕r),'; s = ',⍕s
sp.SetAxesScales(1)
sp.DrawScatterPlot↓9 11∘.○{r F s}¨⍳1000
:EndFor

sp.SaveSvg(⊂'plot.svg')


# MATL, 10 bytes

Xr*+Jr4*^*


Inputs are r, then s. Output is a complex number.

Try it online! Or see the plot for 1000 points at MATL Online! (it takes 10‒15 seconds).

### How it works

Uses method 2 described in the challenge.

Xr   % Push random number with standard Gaussian distribution
*    % Implicit input: r. Multiply
+    % Implicit input: s. Add
J    % Push imaginary unit
r    % Push random number with stantard uniform distribution
4    % Push 4
*    % Multiply
^    % Power
*    % Multiply. Implicit output


# Factor, 56 bytes

[ normal-random-float 2pi random 2dup cos * -rot sin * ]


Try it online!

Verifying correctness:

## Explanation

                      ! 1 0.1
normal-random-float   ! 1.091729295255315
2pi                   ! 1.091729295255315 6.283185307179586
random                ! 1.091729295255315 4.669140230445313
2dup                  ! 1.091729295255315 4.669140230445313 1.091729295255315 4.669140230445313
cos                   ! 1.091729295255315 4.669140230445313 1.091729295255315 -0.04323526873134136
*                     ! 1.091729295255315 4.669140230445313 -0.04720120946224146
-rot                  ! -0.04720120946224146 1.091729295255315 4.669140230445313
sin                   ! -0.04720120946224146 1.091729295255315 -0.9990649185802335
*                     ! -0.04720120946224146 -1.090708439475907


# R, 40 bytes

function(r,s)rnorm(1,r,s)*1i^runif(1,,4)

• You could save a few bytes by using the R 4.1.0+ shorthand function definition notation: \(r,s)rnorm(1,r,s)*1i^runif(1,,4), though TIO doesn't seem to support it. Mar 8 at 10:50
• @user2554330 - Yes, but see comments here. TLDR: changing a single function into \ to save 7 bytes seems boring and I don't usually bother, since it's nice to have a TIO link. Mar 8 at 12:05

# Wolfram Language (Mathematica), 42 bytes

Re[#+#2√Log[16/r^2]I^r]I^r&
r:=4Random[]


Try it online!

RandomVariate@NormalDistribution is costly (and, as noted by Ben Izd, doesn't work with stdev=0), so this uses Box-Muller to generate a normal distribution from two uniform ones.

Sample distributions (N=10000):

# Julia 1.0, 26 bytes

r\s=(randn()s+r)im^4rand()


Try it online!

uses the second formula. output is a complex number. randn gives a random number from a normal distribution (0,1), and rand from a uniform distribtion in [0,1)

1000 points from 10\1:

## Mathematica - 108 bytes

Following the first method, we could have:

fn=AngleVector/@RandomVariate[ProductDistribution[NormalDistribution[#,#2],UniformDistribution[{0,2Pi}]],1000]&;


then visualize it by:

visualize =
Graphics[{PointSize[.01], Point[fn[#, #2]]}, Frame -> True] &;


r=1, s=0.1:

r=3, s=1:

r=100, s=5:

Notes:

• Variance: 0 in NormalDistribution is not supported (could be hacked by having a small number)
• The question asks for one point, not 1000; those 1000 points are generated for illustration. Random[] is a shorter way to access a uniform distribution. You can also save more bytes by using prefix/infix notation and ## instead of #,#2.
– att
Mar 7 at 17:59

# Python 2, 58 bytes

Might not be the shortest, but here is the base-case for python i guess.

Outputs a complex number:

lambda r,s:1j**uniform(0,4)*gauss(r,s)
from random import*


Try it online!

# R, 45 bytes

function(r,s)exp(runif(1)*2i*pi)*rnorm(1,r,s)


Try it online!

Uses the first method in the description, of course using the neat fact that $$\e^{i\theta}=\cos(\theta)+i\sin(\theta)\$$. Longer than Dominic van Essen's answer by 5 bytes, though.

# Java 17, 127 bytes

(r,s)->{double a=new java.util.Random().nextGaussian(r,s),b=Math.random()*Math.PI*2;return new P(a*Math.cos(b),a*Math.sin(b));}


This is a BiFunction<Double, Double, P> where P is a record P(double x, double y) {}

• Try it online needs to update their JDK. Mar 7 at 19:47
• You should include the imports, so the var g=new Random(); should be var g=new java.util.Random();. But you can golf the var g=new java.util.Random();double a=g.nextGaussian(r,s),b=g.nextDouble(Math.PI*2); to double a=new java.util.Random().nextGaussian(r,s),b=Math.random()*Math.PI*2; for -8 bytes. Mar 8 at 7:47
• @KevinCruijssen Thanks. updated. +2 bytes.. length still fits in Java's irritating signed byte 😂 Mar 8 at 14:27
• Two more minor things to golf: (r,s)-> can be r->s-> for -1 byte and return new P(a*Math.cos(b),a*Math.sin(b)); can be return a*Math.cos(b)+","+a*Math.sin(b); for another -3 (with the Bifunction replaced with a currying Function: Function<Double, Function<Double, String>>) Mar 8 at 14:47
• And if you haven't seen them yet, tips for golfing Java and tips for golfing in 'all languages' might be interesting to read through. :) Mar 8 at 14:48

# Desmos, 8070 68 bytes

b=random(1,t)τ
f(r,s,t)=(normaldist(r,s).random(1,t)(cosb,sinb))[1]


Takes an extra argument t as the seed, which is the only way to re-use the function to get different samples without pressing the randomize button.

Try it on Desmos!

-10 bytes thanks to Aiden Chow

-2 bytes thanks to emanrescu A (\tau to τ)

All functions in Desmos are pure, so they can't return a different value when evaluated at different times, even in a list comprehension. This causes an issue with the CGSE policy of functions being re-usable.

There's a randomize button to re-seed all of the random seed-dependent function calls: . This doesn't vibe with me because it requires user interaction to re-seed, but it would allow the following 53-byte submission:

b=random()τ
f(r,s)=(cosb,sinb)normaldist(r,s).random


In this submission, I opted to take the random seed as an extra argument, which is a common design decision in Desmos if a program needs to avoid user action when re-seeding. This is the only way to get different outputs from random functions without the user pressing the randomize button.

• τ is two bytes shorter than \tau. Sorry about the edit, my mistake\ May 3 at 4:04

# 05AB1E, 36 bytes

9°©D(ŸDÄ®>αÅ1*˜Ω®/*+žq·®*ÝΩ®/DÅ¾sÅ½‚*


Inputs in the order $$\r,s\$$.

Try it online. (With the 9 replaced with 3 so it won't time out.)

Explanation:

Uses the first method described in the challenge description.

However, 05AB1E lacks a Gaussian distribution random builtin, as well as a builtin to get a random decimal number given a range. Both of those are therefore done manually.

9°©D(ŸDÄ®>αÅ1*˜Ω®/ # Push a random value with Gaussian distribution within the
# range [-1,1]:
9°                 #  Push 1,000,000,000 (10**9)
©                #  Store it in variable ® (without popping)
D               #  Duplicate it
(              #  Negate the copy
Ÿ             #  Pop both and push an integer-list list in the range
#  [-1000000000,1000000000]
D            #  Duplicate this list
Ä           #  Get the absolute value of each in the copy
®>α        #  Get the absolute difference between each and ®+1
Å1      #  Map each inner value to a list of that many 1s
*     #  Multiply each to the values at the same positions in the
#  remaining list
˜    #  Flatten this list of lists
Ω   #  Pop and push a random integer
®/ #  Divide it by ®
*+                 # Use the inputs to transform it into s*random+r:
*                  #  Multiply it to the (implicit) input s
+                 #  Add the (implicit) input r
žq·®*ÝΩ®/          # Push a random value with uniform distribution within the
# range [0,2π):
žq                 #  Push builtin PI: 3.141592653589793
·                #  Double it to get tau: 6.283185307179586
®*              #  Multiply it by ®
Ý             #  Pop and push an integer list in the range [0,2π®]
Ω            #  Pop and push a random integer from this list
®/          #  Divide it by ®
DÅ¾sÅ½‚*           # Calculate the resulting [x,y] pair using the two random
# values:
D                  #  Duplicate it
Å¾                #  Pop and push its cosine
s               #  Swap so the random value is at the top again
Å½             #  Pop and push its sine
‚            #  Pair them together
*           #  Multiply both to the earlier random value
# (after which the result is output implicitly)