Wednesday, 8 October 2014

Trippy plots and codes...

I noticed wikipedia has a new visualization for the Riemann zeta function. However I'm not sure it's quite best way to show it in every case. This is how it looks:

This is a complex function so each point in 2D-plane has phase and amplitude. In the domain coloring scheme the phase is represented by hue (color) and amplitude by value (intensity). This particular image has some modulo coloring of the amplitude etc. and the difference can be most easily seen by comparing it to these 3D-plots I made where the hue independently represents the phase and amplitude is simply represented by the 3rd dimension.

The second image has phase exaggerated quite a lot to make it clear that for these values the true change of phase is quite different from what it appears in the wikipedia image.

While I was at it, I also plotted the Mandelbrot fractal in terms of the complex value after iteration instead of number of iterations as it's usually done. This too resulted in a rather trippy image. Though obviously hue coloring tends to always results in rather trippy images.

Here hue represents phase angle of the complex number and value (intensity) the amplitude. 3D representation is also rather interesting.
Here only hue is varied.

N = 1024;
s = zeros(N, N);
for y = 1:N
    for x = 1:N
        z0 = (6*x/N-3.5) + j*(6*y/N-3);
        z = 0;
        iteration = 0;
        max_iteration = 256;
        while norm(z) < 8 && iteration < max_iteration
            z = z^2 + z0;
            iteration = iteration + 1;
        s(x, y) = log(iteration);
        S(x, y) = z;


h = angle(S);
h = h-min(min(h));
h = h/max(max(h));

v = abs(S);

hsv(:,:,1) = h;
hsv(:,:,2) = 1;
hsv(:,:,3) = 1;
% hsv(:,:,3) = v/max(max(v));
rgb = hsv2rgb(hsv);

x = 6*[1:N]/N-3.5;
y = 6*[1:N]/N-3;
surf(x, y, v, rgb, 'edgecolor', 'none');
shading interp;

image(x, y, rgb);
imwrite(rgb, 'compmand.png', 'png')

How to test shaders with minimal effort on OSX (though it'll work with relatively minor changes on linux as well)...

import ctypes
from ctypes import *

d = cdll.LoadLibrary("SceneKit.framework/SceneKit")

libSDL = cdll.LoadLibrary("libSDL.dylib")

a = c_char_p("""
uniform float t;
vec4 orb;
float ss;
float map(vec3 p){
 float scale = 1.0;
 orb = vec4(1000.0);
 for( int i=0; i<8;i++ ){
  p = -1.0 + 2.0*fract(0.5*p+0.5);
  float r2 = dot(p,p);
         orb = min( orb, vec4(abs(p),r2) );
  float k = max(ss/r2,0.1);
  p *= k;
  scale *= k;
 return 0.25*abs(p.y)/scale;
float trace( in vec3 ro, in vec3 rd ){
 float maxd = 100.0;
 float precis = 0.001;
 float h=precis*2.0;
   float t = 0.0;
 for(int i=0; i<200; i++){
  if( abs(h)<precis||t>maxd ) continue;
    t += h;
      h = map( ro+rd*t );
     if( t>maxd ) t=-1.0;
 return t;
vec3 calcNormal( in vec3 pos ){
 vec3  eps = vec3(.0001,0.0,0.0);
 vec3 nor;
 nor.x = map(pos+eps.xyy) - map(pos-eps.xyy);
 nor.y = map(pos+eps.yxy) - map(pos-eps.yxy);
 nor.z = map(pos+eps.yyx) - map(pos-eps.yyx);
 return normalize(nor);
void main(void){
 vec2 iResolution=vec2(512,384);
 vec2 p = -1.0 + 2.0*gl_FragCoord.xy / iResolution.xy;
    p.x *= iResolution.x/iResolution.y;
 float time = t*0.25;
    ss = 1.1 + 0.5*smoothstep( -0.3, 0.3, cos(0.1*t) );
 vec3 ro = vec3( 2.8*cos(0.1+.33*time), 0.4 + 0.30*cos(0.37*time),2.8*cos(0.5+0.35*time) );
 vec3 ta = vec3( 1.9*cos(1.2+.41*time), 0.4 + 0.10*cos(0.27*time), 1.9*cos(2.0+0.38*time) );
 float roll = 0.2*cos(0.1*time);
 vec3 cw = normalize(ta-ro);
 vec3 cp = vec3(sin(roll), cos(roll),0.0);
 vec3 cu = normalize(cross(cw,cp));
 vec3 cv = normalize(cross(cu,cw));
 vec3 rd = normalize( p.x*cu + p.y*cv + 2.0*cw );
 vec3 col = vec3(0.0);
 float t = trace( ro, rd );
 if( t>0.0 ){
  vec4 tra = orb;
  vec3 pos = ro + t*rd;
  vec3 nor = calcNormal( pos );
         vec3  light1 = vec3(  0.577, 0.577, -0.577 );
         vec3  light2 = vec3( -0.707, 0.000,  0.707 );
  float key = clamp( dot( light1, nor ), 0.0, 1.0 );
  float bac = clamp( 0.2 + 0.8*dot( light2, nor ), 0.0, 1.0 );
  float amb = (0.7+0.3*nor.y);
  float ao = pow( clamp(tra.w*2.0,0.0,1.0), 1.2 );
  vec3 brdf  = 1.0*vec3(0.40,0.40,0.40)*amb*ao;
  brdf += 1.0*vec3(1.00,1.00,1.00)*key*ao;
  brdf += 1.0*vec3(0.40,0.40,0.40)*bac*ao;
  vec3 rgb = vec3(1.0);
  rgb = mix( rgb, vec3(1.0,0.80,0.2), clamp(6.0*tra.y,0.0,1.0) );
  rgb = mix( rgb, vec3(1.0,0.55,0.0), pow(clamp(1.0-2.0*tra.z,0.0,1.0),8.0) );
  col = rgb*brdf*exp(-0.2*t);
 col = sqrt(col);
 col = mix( col, smoothstep( 0.0, 1.0, col ), 0.25 );

b = d.glCreateShader(35632)
c = d.glCreateProgram();
d.glShaderSource(b, 1, c_long(ctypes.addressof(a)), 0)
d.glAttachShader(c, b);
b = d.glGetUniformLocation(c, "t");

while 1:
 d.glUniform1f(b, c_float(t));
 d.glRecti(-1, -1, 1, 1);

How about this shit (
[bits 32]
; Use Python to load the SDL library, open an SDL window, mmap us into memory,
; then start the intro.

db `python -c"from ctypes import *;S=CDLL('');`
db `CFUNCTYPE(int)(S.mmap(0,1,7,2,'$0',2),0)+167)`
db `(S.SDL_SetVideoMode(640,480,0,0),S.SDL_PollEvent,S.SDL_Flip);" #`
; Global constants
%define XRES          640        ; X resolution of screen
%define YRES          480        ; Y resolution of screen
%define VOXEL_SIZE    4          ; Size of a voxel (2^VOXEL_SIZE)
%define MAXITS        32         ; Maximum number of iterations when marching
%define ROT_SPEED     120        ; Speed of rotation (larger = slower)

; Pseudonyms for stack addresses within main()

%define SURFACE       ebp + 08 ; Pointer to an SDL_Surface
%define SDL_POLLEVENT ebp + 12 ; Pointer to SDL_Pollevent() function
%define SDL_FLIP      ebp + 16 ; Pointer to SDL_Flip() function

    push ebp
    mov  ebp, esp

This fellow does interesting trick to avoid having to make a proper binary on linux, resulting in 512 byte demo.

One can dynamically (or sometimes even not so dynamically) load function addresses in C as well.

#include <dlfcn.h>
typedef void func(int,int,int,int);
        dlopen("SceneKit.framework/SceneKit", RTLD_LAZY);

        void *a = dlsym(RTLD_DEFAULT, "glCreateShader");
        void *b = dlsym(RTLD_DEFAULT, "glCreateProgram");
        void *c = dlsym(RTLD_DEFAULT, "glShaderSource");
        void *d = dlsym(RTLD_DEFAULT, "glCompileShader");
        void *e = dlsym(RTLD_DEFAULT, "glAttachShader");
        void *f = dlsym(RTLD_DEFAULT, "glLinkProgram");
        void *g = dlsym(RTLD_DEFAULT, "glUseProgram");
        void *h = dlsym(RTLD_DEFAULT, "glGetUniformLocation");
        void *i = dlsym(RTLD_DEFAULT, "glUniform1fv");
        void *j = dlsym(RTLD_DEFAULT, "glRecti");
        void *k = dlsym(RTLD_DEFAULT, "printf");

        printf("glCreateShader: %p\n", a);
        printf("glCreateProgram: %p\n", b);
        printf("glShaderSource: %p\n", c);
        printf("glCompileShader: %p\n", d);
        printf("glAttachShader: %p\n", e);
        printf("glLinkProgram: %p\n", f);
        printf("glUseProgram: %p\n", g);
        printf("glGetUniformLocation: %p\n", h);
        printf("glUniform1fv: %p\n", i);
        printf("glRecti: %p\n", j);
        printf("printf: %p\n", k);

 func* fu = (func*)0x10f7f8619;
 printf("%p\n", &printf);
// f(640,480,0,0);
One can also mmap raw binary from a file quite easily and execute it:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <sys/mman.h>

 char buffer[1000];

 FILE *f;
 f = fopen("h64", "rb");
 fseek(f, 0, SEEK_END);
 int size = ftell(f);
 fseek(f, 0, SEEK_SET);
 fread(buffer, size, 1, f);

 printf("Code Size: %d\n", size);
 printf("Prot: %d\n", PROT_EXEC|PROT_READ|PROT_WRITE);
 printf("Flags: %d\n", MAP_PRIVATE|MAP_ANON);

 void* mem_map = mmap(0, sizeof(buffer), PROT_EXEC|PROT_READ|PROT_WRITE, MAP_PRIVATE|MAP_ANON, -1, 0);
 printf("Address: %p\n", mem_map);
 memcpy(mem_map, buffer, sizeof(buffer));


...and the asm (this is for MacOS so if you want to do it for linux, alter it a bit):

; nasm -f bin binary64.asm
[BITS 64]

 mov rax, 0x2000004  ; write
 mov rdi, 1   ; stdout
 lea rsi, [rel msg]
 mov rdx, msg.len


msg: db 'Hello!',0x0A
.len: equ $ - msg

Native OpenGL on Mac OS X:
python -c"from ctypes import *;from Cocoa import *;a=c_char_p(\"float t=gl_Color.x*999.;int m=6;float c=1./(3.*pow(2.,float(m)));vec2 rot(vec2 u,float a){return vec2(u.x*cos(a)-u.y*sin(a),u.y*cos(a)+u.x*sin(a));}void main(void){vec2 u=vec2(640,480);u=-.5*(u-2.*gl_FragCoord.xy)/u.x;u=rot(u,t);u*=sin(t)*.5+1.5;float s=.3;for(int i=0;i<m;i++){u=abs(u)-s;u=rot(u,t);s=s/2.1;}float c=length(u)>c?0.:1.;gl_FragColor=vec4(c,c,c,1.);}\");b=cdll.LoadLibrary(\"SceneKit.framework/SceneKit\");NSApplicationLoad();w=NSWindow.alloc().initWithContentRect_styleMask_backing_defer_(NSMakeRect(0,0,640,480),1,2,1);w.setContentView_(NSOpenGLView.alloc().initWithFrame_pixelFormat_(w.contentRectForFrameRect_(w.frame()),NSOpenGLPixelFormat.alloc().initWithAttributes_([73,5])));w.orderFrontRegardless();d=b.glCreateShader(35632);c=b.glCreateProgram();b.glShaderSource(d,1,c_long(addressof(a)),0);b.glCompileShader(d);b.glAttachShader(c,d);b.glLinkProgram(c);b.glUseProgram(c);
while 1:c+=1;b.glColor3us(c_int(c),0,0);b.glRecti(-1,-1,1,1);b.glSwapAPPLE()"
Looks like one can save 4 bytes in the gzip header technique by just leaving last 4 out. They are crc-32, so they aren't really needed. It'll say truncated input, but works just fine. Previous code gives 675 byte executable (Mac).

echo 'cp \x240 a;sed 1d \x240|zcat>a;./a' > intro
gzip -cn9 > intro1.tmp
dd if=./intro1.tmp of=./intro0.tmp bs=$$((`ls -l intro1.tmp|awk '{ print $$5 }'`-4)) count=1
cat intro0.tmp >> intro
chmod +x intro
rm *.tmp
...\x24 is used in makefile where $ is reserved.

MMAP hello, world in python (binary with nasm):

[bits 32]

db `python -c"from ctypes import *;S=CDLL('');`
db `CFUNCTYPE(int)(S.mmap(0,1,7,2,'$0',2),0)+109)();" #`

 jmp     msg

        mov     eax, 4
        mov     ebx, 1
        pop     ecx
        mov     edx, 7
        int     0x80

msg:    call    start
        db 'Hello!',0xA

Execute binary data in C:

char* code = "\xeb\x13\xb8\x04\x00\x00\x00\xbb\x01\x00\x00\x00\x59\xba\x07\x00\x00\x00\xcd\x80\xc3\xe8\xe8\xff\xff\xff\x48\x65\x6c\x6c\x6f\x21\x0a";

main() {
...stupid blogspot fails to format the text again, but I can't be bothered to give a shit right now.