/** A puzzle from paniq: https://gist.github.com/paniq/c8062db998d573c55cc7e98a1e0f3d48 */ use Fourier.frink a = """71 41 83 70 28 5 3 4 7 71 4 3 4 48 92 26 60 46 15 95 84 93 16 60 18 53 17 15 6 74 94 96 73 27 83 72 16 81 91 3 53 51 92 43 38 15 48 42 71 51 40 37 71 37 73 61 72 45 96 26 6 30 27 71 61 51 48 7 95 83 27 73 6 29 15 16 29 16 27 17 53 16 57 51 26 6 73 95 27 81 53 51 27 6 51 70 84 28 82 92 5 92 27 94 28 65 64 36 71 6 73 73 70 16 29 38 16 35 50 25 64 27 16 17 5 71 7 38 93 4 45 72 16 15 50 27 62 95 80 56 61 83 27 6 29 72 42 92 56 17 83 25 47 38 16 29 4 93 63 18 64 15 6 3 50 26 28 82 73 4 7 42 49 7 17 94 52 40 41 28 83 95 54 6 74 56 26 72 93 72 4 7 16 5 70 49 83 58 82 62 59 53 62 28 50 72 50 44 61 15 49 38 72 71 27 72 30 59 17 61 42 70 62 27 48 15 50 18 27 17 60 82 73 49 3 16 82 49 92 82 39 73 74 42 28 62 93 48 28 39 95 16 93 35 59 35 5 82 46 38 28 34 95 83 55 72 39 84 17 49 84 62 39 94 95 73 26 94 72 4 28 94 46 5 3 15 70 60 72 16 28 82 25 61 64 58 30 48 36 70 59 82 37 48 81 47 72 63 50 95 92 15 93 71 83 5 29 93 94 83 81 29 92 27 69 28 95 5 51 83 4 51 50 94 15 43 95 6 83 81 85 41 62 61 82 72 85 27 15 6 82 81 94 62 47 72 94 81 27 94 85 85 47 72 5 52 36 71 46 45 93 28 59 15 95 41 35 4 52 82 14 18 63 39 70 17 74 93 3 59 62 71 80 45 5 82 3 35 5 58 28 29 27 74 40 5 27 39 4 82 83 29 37 93 82 81 4 17 95 26 59 28 6 94 28 18 15 40 36 52 46 70 6 94 72 94 47 51 70 4 38 39 59 3 48 61 48 28 73 59 29 73 14 64 61 71 48 5 95 72 72 73 42 15 28 95 16 94 46 95 66 81 64 6 40 58 6 16 74 15 16 39 82 40 50 18 94 38 28 17 4 38 92 38 27 93 93 83 60 58 4 18 82 83 65 52 4 70 40 83 82 47 16 94 5 91 82 28 6 83 14 27 94 49 5 49 51 16 72 6 16 28 51 93 27 17 93 64 6 85 38 71 39 84 60 28 83 93 95 92 82 59 84 49 71 54 95 72 38 28 49 72 71 26 46 72 42 62 28 16 28 4 71 34 39 84 49 28 8 16 38 49 62 38 92 26 27 72 5 4 71 80 94 51 52 27 83 37 37 17 92 16 27 6 5 52 37 71 93 62 51 69 42 7 28 82 62 41 93 47 72 81 71 93 28 37 39 17 82 41 46 61 71 4 72 27 95 94 7 85 83 83 92 92 94 95 53 82 49 54 84 29 71 49 4 39 46 70 6 84 93 5 35 94 96 93 3 83 16 27 37 29 60 29 5 27 61 15 5 52 42 17 16 82 95 17 70 82 38 17 47 5 83 63 94 65 81 17 16 81 95 81 5 18 70 83 93 84 95 83 39 71 53 47 40 70 63 14 40 3 63 38 83 59 50 14 14 38 39 28 37 84 5 72 15 39 60 60 70 51 49 3 94 82 25 65 94 17 84 48 93 5 59 63 49 48 46 83 48 5 26 47 39 83 29 95 49 50 5 84 52 72 48 40 94 16 83 63 82 58 36 26 62 7 94 84 4 6 47 4 70 49 96 73 49 29 16 15 16 39 71 26 72 26 54 93 49 3 52 94 3 39 94 49 59 93 50 70 72 3 15 47 49 52 93 83 15 30 7 26 29 70 5 83 49 83 94 16 83 5 38 27 6 40 96 26 63 16 35 35 73 4 4 49 39 94 73 6 5 72 18 73 94 59 94 37 61 27 15 27 81 30 62 48 15 47 36 5 16 37 93 60 29 81 95 48 7 16 84 16 72 27 6 58 4 62 5 83 83 41 49 82 15 5 28 94 63 7 82 57 71 16 73 62 60 41 58 40 92 73 84 84 4 82 40 81 4 36 65 50 15 17 47 5 71 17 17 64 61 17 36 49 82 61 93 6 46 62 5 81 38 39 70 69 94 71 3 71 8 95 94 84 49 17 72 74 39 91 5 70 93 82 6 16 50 4 82 72 54 73 92 40 83 61 17 62 25 5 83 48 60 83 81 49 84 17 45 70 16 93""" b = eval[split[%r/\s+/, a]] println[length[b]] use correlation.frink //c = autocorrelation[b] //println[c] for size = 1 to sqrt[length[b]] { g = new graphics for i = rangeOf[b] { p = b@i g.color[p/100, p/100, p/100] g.fillRectSize[i div size, i mod size, 1, 1] } // g.show[] g.write["wow$size.png", 240, undef] } for size = 1 to 1 { g = new graphics pl = new polyline g.color[0,0,0] for i = rangeOf[b] { p = b@i pl.addPoint[i,p] } g.add[pl] g.show[] } f = DFT[b] //println[f] g = new graphics // fprime will contain just selected frequency components. fprime = new array[] len = length[f] zeroUnit = 0 f@0 threshold = .1 vscale = 20 c = 0 str = "" for y = f { x = (c + len/2) mod len // Draw low-freq components together in middle // Draw phases g.color[0,1,0] phase = arctan[Im[y],Re[y]] g.fillRectSides[x,0,x+1,-5 phase] // Draw magnitudes g.color[0,0,0,.7] mag = abs[y] if (c != 0) // Don't draw DC component because it's huge. g.fillRectSides[x,0,x+1,-vscale mag] if (mag > threshold) // Decide which components to keep. { daysphase = (phase/(2 pi) * len) mag2 = 2 mag println["Used component $c, mag is $mag, phase is $daysphase days"] if c == 0 str = format[mag2/2, 1, 6] if c > 0 && c < len/2 str = str + " +\n " + format[mag2,1,6] + " cos[2 pi $c/$len d - " + format[phase,1,6] + "]" fprime.push[y] } else fprime.push[zeroUnit] c = c + 1 } println[str] g.show[] g.write["wowFourier.png", 1024, undef]