How can I generate random vectors with conditions? (Using Python)

Question

I'm a beginner in programming and I'm looking for a nice idea how to generate ​n random 2D vectors [a1,b1] [a2,b2] ... [an,bn] that satisfy the 3 conditions below.

(1) a1+a2+...+an=0

(2) b1+b2+...+bn=0

(3) (a1**2+b1**2)=(a2**2+b2**2)=...=(an**2+bn**2)=1/10000

example : for n=4

[0.00258819045102518, 0.00965925826289069]

[0.00866025403784439, 0.00499999999999999]

[-0.00258819045102518, -0.00965925826289069]

[-0.00866025403784439, -0.00499999999999999]

These four vectors satisfy all the conditions.

When n=4, I found the 4 vectors geometrically mean rhombus with a length of one side of 1/100. So I could solve it like this. (I used the Python 3)

import random

from math import *

import math


k1=pi*random.random()*2

k2=pi*random.random()*2

a1=math.cos(k1)/100

b1=math.sin(k1)/100

a2=math.cos(k2)/100

b2=math.sin(k2)/100

c1=[a1,b1]

c2=[a2,b2]

c3=[-a1,-b1]

c4=[-a2,-b2]

print(c1)

print(c2)

print(c3)

print(c4)

But I failed to solve this problem for general n other than 4. Is there anyone who can help me? Thanks for reading this.

Answer 1

Interesting puzzle. You did all hard work for even N (2, 4, 6, etc). After that you just need to wrap this logic inside loop:

import random
from math import pi, cos, sin

# insert how many iterations you want
n = int(input('insert N: '))

# is number of iterations odd?
is_odd = bool(n % 2)

# set number of iterations to even
# and devide it to 2
# because we need to set only half of all vectors
n = int(n / 2)

positive_vectors = []

# set N/2 vectors
for _ in range(0, n):
    k = pi*random.random()*2
    a = cos(k)/100
    b = sin(k)/100
    positive_vectors.append((a, b))

negative_vectors = []

# generate negative vectors with
# reversed values to satisfy 1 and 2 conditions (-0.4 + 0.4 = 0)
for el in positive_vectors:
    negative_vectors.append((-el[0], -el[1]))

# merge positive and negative vectors
vectors = positive_vectors + negative_vectors

if is_odd:
    # do something to satisfy conditions for odd N
    # we need to add one more last Nth element
    pass

# print results
sum_a = 0
sum_b = 0
i = 1

print(' ---')
for vector in vectors:
    sum_a = sum_a + vector[0]
    sum_b = sum_b + vector[1]
    print(vector)
    print(f"a{i}**2 + b{i}**2 =", "%.5f" % (vector[0]**2 + vector[1]**2))
    print(' ---')
    i += 1

print('a1+a2+...+an =', round(sum_a, 16))
print('b1+b2+...+bn =', round(sum_b, 16))

Output will be:

insert N: 6
 ---
(0.009997486005014973, -0.00022421770565627118)
a1**2 + b1**2 = 0.00010
 ---
(0.0030141777667107295, -0.00953492172965603)
a2**2 + b2**2 = 0.00010
 ---
(0.0010585510434668018, -0.009943815650361553)
a3**2 + b3**2 = 0.00010
 ---
(-0.009997486005014973, 0.00022421770565627118)
a4**2 + b4**2 = 0.00010
 ---
(-0.0030141777667107295, 0.00953492172965603)
a5**2 + b5**2 = 0.00010
 ---
(-0.0010585510434668018, 0.009943815650361553)
a6**2 + b6**2 = 0.00010
 ---
a1+a2+...+an = 0.0
b1+b2+...+bn = 0.0

However if N is odd (3, 5, 7, etc) this solution will not work. I guess that we need to change last even element in vectors and add one last odd element to satisfy 1 and 2 conditions.

Answer 2

How about this?

import random as rd
import numpy as np
import math as ma
import matplotlib.pyplot as plt
from math import *

n0 = int(input())
n=n0

def R(b):
  a=[0]*b
  c=rd.random()
  for i in range(b):
    a[i]=(i/b+c)*2*pi
  
  return a



def find(data, target):
  a = []
  b = data
  while True:
    try:
      a.append(b.index(target) + (a[-1]+1 if len(a)!=0 else 0))
      b = data[a[-1]+1:]
    except:
      break     
  return a


def P(n):
    p = []
    if n < 2:
        return p
    for i in range(2, n+1):
        a = True
        for j in p:
            if i % j == 0:
                a = False 
                break
            elif j > i**0.5:
                break
        if a:
            p.append(i)
    return p
    


p=P(n)
p.reverse()
print(p)



p0=len(p)
q=[0]*p0
r=list(range(p0))


 

i=0
while np.dot(p,q) != n0 :
  if i == p0:
    k0=find(q,0)
    k1=list(set(r)-set(k0))
    
    if k1==[]:
      i=0
    else:
      
      kk0=rd.choice(k1)
      n=n+p[kk0]
      q[kk0] -=1
      i= 0
      
  if i < p0:
    c1=n//p[i]
    c=rd.randint(0,c1)
    q[i]=q[i]+c
    n=n-p[i]*c
    i+=1
  
  else :
    print("Error")
    print(i)
    
print(q)
print(r)


t0=find(q,0)
t1=list(set(r)-set(t0))
print(t1)
t2=range(len(t1))

e1=[]
e2=[]

for i in t2:
  j=q[t1[i]]
  
  for x in range(j):
    R0=R(p[t1[i]])
    for k in range(p[t1[i]]):
      print([0.01*ma.cos(R0[k]),0.01*ma.sin(R0[k])])
      e1.append(0.01*ma.cos(R0[k]))
      e2.append(0.01*ma.sin(R0[k]))

print(sum(e1))
print(sum(e2))      
plt.scatter(e1, e2)
plt.show()

As I observed,the sum of correct n vectors with a length of 0.01 must be zero vector. And when I matched all the starting points of the vector to the origin, I saw a clue to the solution.

For example, when n is 25, 25 can be expressed by prime numbers as follows.

25=0*23+0*19+0*17+0*13+0*11+1*7+2*5+0*3+4*2

From this, we can see that each endpoint of the correct n position vectors of length 0.01 consists of (a), (b), and (c).

(a) End points of the position vectors corresponding to the vertex of the one regular hexagon. Rotating the roots of equation "z^7=1" at any angle c1 on a complex plane.

(b) End points of the position vectors corresponding to the vertex of the two different regular pentagon. Rotating the roots of equation "z^5=1" at any angle c1 and c2 on a complex plane.

(c) Four end points pair of originally symmetric position vectors. Rotating the roots of equation "z^2=1" at any angle c1, c2, c3, and c4 on a complex plane.

Thank you for reading this.

When input is 21, Output will be :

[19, 17, 13, 11, 7, 5, 3, 2]
[1, 0, 0, 0, 0, 0, 0, 1]
[0, 1, 2, 3, 4, 5, 6, 7]
[0, 7]
[0.009953491621510826, 0.0009633299229826691]
[0.009101391276086286, 0.004143027496840982]
[0.007263013963261475, 0.006873763755721376]
[0.004637576390244042, 0.008859621054235392]
[0.0015095854559317348, 0.0098854009403382]
[-0.0017819924861625445, 0.009839944246755682]
[-0.0048804638919187635, 0.008728176911570776]
[-0.007450061306785754, 0.006670576176398388]
[-0.009212328979450691, 0.0038901150078594093]
[-0.009976297863179998, 0.0006880995168652664]
[-0.009659180075222993, -0.0025884822337453183]
[-0.008295340248494279, -0.005584561769889459]
[-0.006032571590375244, -0.007975467384862]
[-0.003116080195446368, -0.009502107356557645]
[0.00013808683962512758, -0.009999046555783323]
[0.00337729002298516, -0.009412433909496783]
[0.006250511428300667, -0.007805838000156213]
[0.008446392933682103, -0.005353358423442713]
[0.009726976705409207, -0.0023207593956346864]
[-0.009906386205501586, 0.0013651052514175907]
[0.009906386205501588, -0.0013651052514175876]
1.734723475976807e-18
-6.938893903907228e-18

when input is 4, Output will be :

[3, 2]
[0, 2]
[0, 1]
[1]
[-0.006464284713103015, -0.007629745942555537]
[0.006464284713103016, 0.007629745942555536]
[0.0015619505179010943, -0.009877262301854118]
[-0.0015619505179010932, 0.00987726230185412]
1.951563910473908e-18
0.0