Fifty Shades of J/Chapter 35
Table of Contents ... Glossary ... Previous Chapter ... Next Chapter
Principal Topics
- ?. (roll) o. (circle functions) : (grade down) ~ (reflex), simulation, random numbers, random angles, accept-reject technique
In the world, that is. Lots, I know, but let me ask rather what proportion of all triangles in the world are obtuse-angled? Since any two of the angles, say a and b, of a triangle can be determined independently, the every possible triangle corresponds to a point within the large right-angled triangle in the diagram below
Acute-angled triangles must have a+b>90o and so their points all lie within the inner triangle, from which it follows that 75% of all triangles are obtuse-angled. Well let’s ask J and see what he/she/it says using statistical Monte-Carlo methods. Assume that a point is represented by a 2-list of numbers such as 4 7. The verb `py` calculates the square of the distance between two such points
py=.(+/@:*:@:-/)"2 NB. square of dist between 2 points py 4 7,:1 3 NB. (3,4,5) triangle 25
A triangle is a 3-list of 2-lists (points). To calculate the squares of the three sides, take the points in pairs and apply `py` :
sides =: 2 [\ 4&$ ]t=: 3 2$_3 _2 _1 6 0 2 _3 _2 _1 6 0 2 <"2 sides t ┌─────┬────┬─────┐ │_3 _2│_1 6│ 0 2│ │_1 6│ 0 2│_3 _2│ └─────┴────┴─────┘
A triangle is obtuse if the square of the largest side exceeds the sum of the squares of the other two, so define ifbigone which uses reflexive grade down to test whether the largest of a set of numbers is greater than the sum of all the others
ifbigone =: ({. > +/@}.)@\:~ NB. 1 if biggest > +/rest
ifbigone 49 43 100
1
This leads naturally to
obtest =: ifbigone@:(py"2@sides)"2 NB. 1 if obtuse obtest t 1
randtri1 generates y random triangles with vertices in the unit square:
randtri1 =: 0 ?@$~ ,&3 2 <"2 randtri1 3 ┌─────────────────┬─────────────────┬─────────────────┐ │0.577101 0.430063│0.356508 0.879054│0.190733 0.629572│ │0.689446 0.610395│0.898293 0.11287│ 0.45149 0.964217│ │0.401262 0.267585│0.184275 0.444405│0.482228 0.170644│ └─────────────────┴─────────────────┴─────────────────┘
Everything is now in place to test a single random triangle with given resolution :
trial1=.obtest@:randtri1 NB. 1 if obtuse trial1 8 0 0 1 1 1 0 0 1
Of the triangles in this particular trial four of eight were obtuse-angled. Here are a thousand trials:
+/ trial1 1000 740
... and a few more :
+/ trial1 1000 739 +/ trial1 1000 708 +/ trial1 1000 736
By now you should be getting a feel for where the true value lies, somewhere between 72 and 75 per cent. Another way to conduct the simulation is to generate random angles rather than random sides. Start by writing a verb which generates a number of random values (angles) between 0 and 2π to use as angles from the origin to the vertex:
rnd =: ?@$&0 NB. y rand uniform nos. in (0,1) TAU =: 2p1 NB. tauday.com randang =: TAU * rnd NB. random angles in radians randang 3 5.66985 3.78605 2.69206
The sines and cosines of each angle generate a point within the unit circle, and so random angles can be transformed to random points within the unit circle by
cs =: 2 1 o.~/~ randang cs 8 0.515802 _0.856708 0.19751 _0.980301 0.590711 0.806883 _0.659526 _0.751682 _0.451618 _0.892211 _0.956942 _0.29028 _0.980289 _0.197569 _0.4976 0.867407
and then to y random triangles in a unit circle by
randtri2 =: [: cs ,&3 <"2 randtri2 4 ┌──────────────────┬──────────────────┬─────────────────────┬────────────────────┐ │_0.718971 0.69504│ 0.503937 _0.86374│0.0327604 _0.999463│_0.998892 _0.0470605│ │ 0.96262 0.270857│_0.795225 0.606314│ _0.91815 _0.396232│_0.554286 0.832326│ │ 0.442591 0.896724│_0.401744 0.915752│ 0.999965 _0.00831741│ 0.727516 _0.686091│ └──────────────────┴──────────────────┴─────────────────────┴────────────────────┘
All is in place for determining whether these triangles are obtuse or not, and for doing some more experiments
([: +/ trial2)@:(1000"_) ^:(>: i. 4) '' 763 710 754 736
A little bit higher than our previous estimate, but close to the theoretical 75%. Maybe sampling from a square frame includes disproportionately ‘skinny’ acute-angled triangles which have one vertex tucked into one of the far corners. To test this conjecture, sample the coordinates within a square, but filter only those triangles all of whose points lie inside the unit circle. The polynomial transforms the triangle vertices to have absolute value <: 1 . All of these must have length <: 1.
ssq =: +/@:*: NB. sum of squares randtri3 =: (#~ ([: (*./"1) 1 >: ssq"1))@:(_1 2 p. randtri1) NB. graphics shows the remaining points are within the unit circle require 'plot' 'point'plot _2 j./\ , randtri3 999
As before everything is in place for the next experiment. Since randtri3 generates fewer triangles than requested, include the average :
(+/ % #) trial3 10000 0.729905
Not much different from `trial1` – ah well, that’s one conjecture gone for a bang!
Another accept-reject technique for random triangles consists of generating the values of three random sides in order, and rejecting the result if the sides fail to form a triangle because one of them is bigger than the sum of the other two :
randtri4 =: (#~ (>./ < +/ - >./)"1)@:([: rnd ,&3) randtri4 9 0.482861 0.631554 0.575322 0.790549 0.74182 0.983299 0.21859 0.814304 0.746637 0.758181 0.448893 0.836001
Since randtri4 returns lengths of sides the obtuseness test is applied to the squares of the sides :
trial4 =: ifbigone"1@:*:@randtri4 (+/ % #) trial4 99999 0.57199
Seems I’ve mislaid nearly 20% of the world’s obtuse-angled triangles (careless of me!) Another possibility based on sides is to choose two sides at random, say these were 6 and 10. The range of acceptable values for a third side is then from 10-6 to 10+6, or more generally from m-n to m+n where m and n are the lengths of the sides in decreasing order. The third side can then be taken as (m-n) plus a number drawn at random from 0 to 2n:
avg =: +/ % # rnd =: ?@$&0 third =: , (>./ + <./ * _1 2 p. ?@:0:) rantri5 =: [: third"1 [: rnd ,&2 trial5 =: ifbigone"1@:*:@rantri5 avg trial5 100000 0.73331
Let’s get back to the original question . How many obtuse-angled triangles did I say there were? ...
Code Summary
trial1=: obtest@randtri1 NB. first random trial
obtest =: ifbigone@:(py"2@sides)"2 NB. 1 if obtuse
ifbigone =: ({. > +/@}.)@\:~ NB. 1 if biggest > +/rest
py =: (+/@:*:@:-/)"2 NB. sq of dist between 2 pts
sides =: 2 [\ 4&$ NB. sides from coords
randtri1 =: 0 ?@$~ ,&3 2 NB. generate y triangles in unit square
trial2 =: obtest@randtri2 NB. second random trial
randtri2 =: [: cs ,&3
cs =: 2 1 o.~/~ randang NB. cos/sin of rand angle
TAU =: 2p1
randang =: TAU * rnd NB. random angles in radians
rnd =: ?@$&0 NB. random values on (0,1) with shape y
trial3 =: obtest@randtri3 NB. third random trial
randtri3 =: (#~ ([: (*./"1) 1 >: ssq"1))@:(_1 2 p. randtri1)
NB. generate at most y triangles in unit square,
NB. discarding those outside the inscribed circle.
ssq =: +/@:*: NB. sum of squares
trial4 =: ifbigone@:*:@randtri4 NB. fourth random trial
randtri4 =: (#~ (>./ < +/ - >./)"1)@:([: rnd ,&3)
NB. generate at most y triangles by random length sides
NB. discarding triples which can not form triangles.
trial5 =: ifbigone@:*:@rantri5 NB. fifth random trial
rantri5 =: [: third"1 [: rnd ,&2 NB. choose 2 sides, 3rd side randomly chosen as plausible
third =: , (>./ + <./ * _1 2 p. ?@:0:) NB. join feasible third side