[Tutor] greater precision?

[John Collins]

Hi Steve,
Thanks. From
 >>>mkpoints.py 32 32.txt

here's a sample of the output

-0.396087591388 -0.781284022758 0.482400140683
-0.967387012461 -0.0838047084421 -0.239037944614
0.0208969821213 -0.489420208746 0.871797668848
0.887250003871 -0.258893773768 -0.38178717178
0.426352071227 -0.457758408728 -0.780180203927
0.612061168992 -0.280383142359 0.739436555016
-0.887250003871 0.258893773768 0.38178717178
-0.0971745158475 -0.994342015264 -0.0429076933695
-0.120898756509 -0.759794654167 -0.638823586113
-0.0208969821213 0.489420208746 -0.871797668848
0.482396765336 -0.834880934468 -0.265079584379
0.66383194755 0.726669941038 0.176855710123


from
 >>>nfaces.py 32.txt 32fa.txt

here's a sample of the output

point vertex_0_2_12 -0.0288974102653 -0.851924110364 0.66948446135
point vertex_13_24_27 -0.122445373457 1.00045419254 0.398644871129
point vertex_0_7_15 -0.482610585141 -0.963539328269 0.1161816686
point vertex_3_10_17 0.848541556491 -0.639691741968 -0.213520247975
point vertex_1_6_14 -1.05498774772 0.248634080415 -0.00106786881656
point vertex_13_24_31 -0.291794484064 0.826881428947 0.637135994637
point vertex_1_6_25 -1.07023716261 -0.0479998647263 0.164643927545
point vertex_3_17_28 1.05498774772 -0.248634080415 0.00106786881656
point vertex_1_15_25 -0.975134617659 -0.4675255693 -0.073154040374
point vertex_4_8_10 0.291794484064 -0.826881428947 -0.637135994637
point vertex_9_19_20 -0.400242088946 0.638705226984 -0.778897359983
normal face_30 -0.617395768043 -0.359222235879 -0.699844161834
face face_31 vertex_21_29_31
face face_31 vertex_13_21_31
face face_31 vertex_13_24_31
face face_31 vertex_6_24_31
face face_31 vertex_6_29_31
normal face_31 -0.426352071227 0.457758408728 0.780180203927


As you can see, there are exactly 12 significant figures.

BTW, I did not write either script, I can post them, if you think
doing so may show where it may be possible to escape to "decimal"
as you say, perhaps?

> "Decimal" has a specific meaning in Python different to how you
> appear to be using it. It looks to me like you are working with
> floats rather than decimals.
Right.
>> from math import pi, asin, atan2, cos, sin, sqrt
>> from math import pi, asin, atan2, cos, sin, sqrt
Two different scripts.
> "crosspoint" appears to be a custom Python program that we know
> nothing about, apart from what you tell us.
I posted it in it's entirety.
> "==" is used for equality testing:
> x == 2
> returns False, because x does not equal 2; but
> x == 1
> returns True, because x does currently equal 1.
Ah! Thanks!
> Almost. The "from math import blah..." line extracts a bunch of maths
> functions from the math library and makes them available to your
> program.
I'm seeing this now. Python does minimal things, extras are called in
only as needed.
> The number of significant figures is more fundamental that that; your
> operating system understands about floating point numbers (so-called
> "C doubles") with 53 *binary* significant figures, or about 17 *decimal*
> figures. So I'm not sure why you think there are only 12.
Sure. But see the above, from Py 2.7.2, it's precisely 12 sig.figs.!
> You need to show us how the output is generated in the first place.
Done, above.
> Not a chance without some major changes to your program.
I was afraid of this.
> I can't say how major without seeing more of your program.
Well, OK, here's the first, mkpoints.py
#!/usr/bin/env python

# Distribute a set of points randomly across a sphere, allow them
# to mutually repel and find equilibrium.

import sys
import string
import random
from math import pi, asin, atan2, cos, sin, sqrt

args = sys.argv[1:]

if len(args) > 0:
     n = string.atoi(sys.argv[1])
     args = args[1:]
else:
     n = 7

if len(args) > 0:
     outfile = open(args[0], "w")
     args = args[1:]
else:
     outfile = sys.stdout

points = []

def realprint(a):
     for i in range(len(a)):
         outfile.write(str(a[i]))
         if i < len(a)-1:
             outfile.write(" ")
         else:
             outfile.write("\n")

def pprint(*a):
     realprint(a)

for i in range(n):
     # Invent a randomly distributed point.
     #
     # To distribute points uniformly over a spherical surface, the
     # easiest thing is to invent its location in polar coordinates.
     # Obviously theta (longitude) must be chosen uniformly from
     # [0,2*pi]; phi (latitude) must be chosen in such a way that
     # the probability of falling within any given band of latitudes
     # must be proportional to the total surface area within that
     # band. In other words, the probability _density_ function at
     # any value of phi must be proportional to the circumference of
     # the circle around the sphere at that latitude. This in turn
     # is proportional to the radius out from the sphere at that
     # latitude, i.e. cos(phi). Hence the cumulative probability
     # should be proportional to the integral of that, i.e. sin(phi)
     # - and since we know the cumulative probability needs to be
     # zero at -pi/2 and 1 at +pi/2, this tells us it has to be
     # (1+sin(phi))/2.
     #
     # Given an arbitrary cumulative probability function, we can
     # select a number from the represented probability distribution
     # by taking a uniform number in [0,1] and applying the inverse
     # of the function. In this case, this means we take a number X
     # in [0,1], scale and translate it to obtain 2X-1, and take the
     # inverse sine. Conveniently, asin() does the Right Thing in
     # that it maps [-1,+1] into [-pi/2,pi/2].

     theta = random.random() * 2*pi
     phi = asin(random.random() * 2 - 1)
     points.append((cos(theta)*cos(phi), sin(theta)*cos(phi), sin(phi)))


# For the moment, my repulsion function will be simple
# inverse-square, followed by a normalisation step in which we pull
# each point back to the surface of the sphere.

while 1:
     # Determine the total force acting on each point.
     forces = []
     for i in range(len(points)):
         p = points[i]
         f = (0,0,0)
         ftotal = 0
         for j in range(len(points)):
             if j == i: continue
             q = points[j]

             # Find the distance vector, and its length.
             dv = (p[0]-q[0], p[1]-q[1], p[2]-q[2])
             dl = sqrt(dv[0]**2 + dv[1]**2 + dv[2]**2)

             # The force vector is dv divided by dl^3. (We divide by
             # dl once to make dv a unit vector, then by dl^2 to
             # make its length correspond to the force.)
             dl3 = dl ** 3
             fv = (dv[0]/dl3, dv[1]/dl3, dv[2]/dl3)

             # Add to the total force on the point p.
             f = (f[0]+fv[0], f[1]+fv[1], f[2]+fv[2])

         # Stick this in the forces array.
         forces.append(f)

         # Add to the running sum of the total forces/distances.
         ftotal = ftotal + sqrt(f[0]**2 + f[1]**2 + f[2]**2)

     # Scale the forces to ensure the points do not move too far in
     # one go. Otherwise there will be chaotic jumping around and
     # never any convergence.
     if ftotal > 0.25:
         fscale = 0.25 / ftotal
     else:
         fscale = 1

     # Move each point, and normalise. While we do this, also track
     # the distance each point ends up moving.
     dist = 0
     for i in range(len(points)):
         p = points[i]
         f = forces[i]
         p2 = (p[0] + f[0]*fscale, p[1] + f[1]*fscale, p[2] + f[2]*fscale)
         pl = sqrt(p2[0]**2 + p2[1]**2 + p2[2]**2)
         p2 = (p2[0] / pl, p2[1] / pl, p2[2] / pl)
         dv = (p[0]-p2[0], p[1]-p2[1], p[2]-p2[2])
         dl = sqrt(dv[0]**2 + dv[1]**2 + dv[2]**2)
         dist = dist + dl
         points[i] = p2

     # Done. Check for convergence and finish.
     sys.stderr.write(str(dist) + "\n")
     if dist < 2.10005e-16:
         break

# Output the points.
for x, y, z in points:
     pprint(x, y, z)


The nfaces.py used on the output from this is a tad longer!


> Well, yes, but only with some significant changes
Well, it seems I may have to learn "decimal" python, and insert it
into the above, and nfaces.py.

John.