Hope it's ok to send a few lines of pari/gp code: /* functions correct to order = 4^N-1 */ fpp(y)= { local( t, Y ); t = 1; Y=y; for (k=1, N, t *= (+1+Y+Y^2-Y^3); Y=Y^4; ); return( t-1 ); } fmm(y)= { local( t, Y ); t = 1; Y=y; for (k=1, N, t*= (-1+Y+Y^2+Y^3); Y=Y^4; ); return( t*y-Y ); } hhpm(y)= { local( tp, tm ); tp = fpp(y); tm = fmm(y); return( ((tp+tm) + I*(tp-tm))/2 ); } start gp, copy & paste above into session ... ? N=4; ? taylor(hhpm(y), y) y + I*y^2 - y^3 + I*y^4 + I*y^5 + y^6 - I*y^7 + y^8 + I*y^9 + y^10 - I*y^11 - I*y^12 - y^13 - I*y^14 + y^15 + O(y^16) Now start at the origin and for each coefficient draw a one-unit line to the right/left/up/down when the coefficient is +1/-1/+I/-I. ... or use: ? taylor(hhpm(y)/(1-y), y) y + (1 + I)*y^2 + I*y^3 + 2*I*y^4 + 3*I*y^5 + (1 + 3*I)*y^6 + (1 + 2*I)*y^7 + (2 + 2*I)*y^8 + (2 + 3*I)*y^9 + (3 + 3*I)*y^10 + (3 + 2*I)*y^11 + (3 + I)*y^12 + (2 + I)*y^13 + 2*y^14 + 3*y^15 + O(y^16) That is connect the points in the complex plane (given as coefficients). The following functions gives the turns of the above curve: hht(y)= { /* correct to order = 4^N-1 */ local( t, Y, F, s, p ); t = 1; Y=y; p = 1; F = y + y^2; for(k=2, N, Y = Y^4; t = -F + Y*F + Y^2*F - Y^3*F; p *= 4; if ( 0==(k%2), t += y^(1*p-1); t += y^(3*p); t -= (y+1)*y^(2*p-1); , /* else */ t += y^(1*p); t += y^(3*p-1); ); F = t; ); return( F-(N%2) ); } If you see how to make the latter algorithm more regular, please let me now. Trivial? Well-known? Nice? Boring? Fun? -- p=2^q-1 prime <== q>2, cosh(2^(q-2)*log(2+sqrt(3)))%p=0 Life is hard and then you die.