extra/rosetta-code/pythagorean-triples/pythagorean-triples.factor

   1 ! Copyright (c) 2012 Anonymous
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: accessors arrays formatting kernel literals math
   4 math.functions math.matrices ranges sequences ;
   5 IN: rosetta-code.pythagorean-triples
   6
   7 ! http://rosettacode.org/wiki/Pythagorean_triples
   8
   9 ! A Pythagorean triple is defined as three positive integers
  10 ! (a,b,c) where a < b < c, and a2 + b2 = c2. They are called
  11 ! primitive triples if a,b,c are coprime, that is, if their
  12 ! pairwise greatest common divisors gcd(a,b) = gcd(a,c) = gcd(b,c)
  13 ! = 1. Because of their relationship through the Pythagorean
  14 ! theorem, a, b, and c are coprime if a and b are coprime
  15 ! (gcd(a,b) = 1). Each triple forms the length of the sides of a
  16 ! right triangle, whose perimeter is P = a + b + c.
  17
  18 ! Task
  19
  20 ! The task is to determine how many Pythagorean triples there
  21 ! are with a perimeter no larger than 100 and the number of these
  22 ! that are primitive.
  23
  24 ! Extra credit: Deal with large values. Can your program handle
  25 ! a max perimeter of 1,000,000? What about 10,000,000?
  26 ! 100,000,000?
  27
  28 ! Note: the extra credit is not for you to demonstrate how fast
  29 ! your language is compared to others; you need a proper algorithm
  30 ! to solve them in a timely manner.
  31
  32 CONSTANT: T1 {
  33   {  1  2  2 }
  34   { -2 -1 -2 }
  35   {  2  2  3 }
  36 }
  37 CONSTANT: T2 {
  38   {  1  2  2 }
  39   {  2  1  2 }
  40   {  2  2  3 }
  41 }
  42 CONSTANT: T3 {
  43   { -1 -2 -2 }
  44   {  2  1  2 }
  45   {  2  2  3 }
  46 }
  47
  48 CONSTANT: base { 3 4 5 }
  49
  50 TUPLE: triplets-count primitives total ;
  51
  52 : <0-triplets-count> ( -- a ) 0 0 \ triplets-count boa ;
  53
  54 : next-triplet ( triplet T -- triplet' )
  55     [ 1array ] [ mdot ] bi* first ;
  56
  57 : candidates-triplets ( seed -- candidates )
  58     ${ T1 T2 T3 } [ next-triplet ] with map ;
  59
  60 : add-triplets ( current-triples limit triplet -- stop )
  61     sum 2dup > [
  62     /i [ + ] curry change-total
  63     [ 1 + ] change-primitives drop t
  64     ] [ 3drop f ] if ;
  65
  66 : all-triplets ( current-triples limit seed -- triplets )
  67     3dup add-triplets [
  68         candidates-triplets [ all-triplets ] with swapd reduce
  69     ] [ 2drop ] if ;
  70
  71 : count-triplets ( limit -- count )
  72     <0-triplets-count> swap base all-triplets ;
  73
  74 : pprint-triplet-count ( limit count -- )
  75     [ total>> ] [ primitives>> ] bi
  76     "Up to %d: %d triples, %d primitives.\n" printf ;
  77
  78 : pyth ( -- )
  79     8 [1..b] [ 10^ dup count-triplets pprint-triplet-count ] each ;