[Tutor] greater precision?
John Collins
john at netcore.com.au
Mon Oct 29 11:54:27 CET 2012
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.
More information about the Tutor
mailing list