#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <windows.h>
#include <sysinfoapi.h>
#define D __device__
#define F float
#define I int
struct v {F x,y,z;D v operator+(v r){return v(x+r.x,y+r.y,z+r.z);}D v operator*(F r){return v(x*r,y*r,z*r);}D F operator%(v r){return x*r.x+y*r.y+z*r.z;}D v(){}D v operator^(v r){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r.y-y*r.x);}D v(F a,F b,F c){x=a;y=b;z=c;}D v operator!(){return *th*rsqrtf(*th%*th);}};D I g[]={247570,280596,280600,249748,18578,18577,231184,16,16};__shared__ I gs;D F R(){gs=(214013*gs+2531011);return __fdividef(((gs>>16)&0x7FFF),(66635.0f));}D I TR(v O,v d,F &t,v &n){t=1e9;I m=0;F p=__fdividef((-O.z),(d.z));if(.01<p){t=p;n=v(0,0,1);m=1;}for(I k=19;k--;)for(I j=9;j--;)if(g[j]&1<<k){v p=O+v(-k,0,-j-4);F b=p%d;F c=p%p-1;F q=b*b-c;if(q>0){F s=-b-sqrtf(q);if(s<t&&s>0.01f){t=s;n=!(p+d*t);m=2;}}}return m;}D v Sample(v O,v d){F at=1.0f;v pc(0,0,0);for(I r=0 ;r<4 ;r++,at /= 2){F t;v n;I mt=TR(O,d,t,n);if(!mt){return pc+v(0.7f,0.6f,1.0f)*__powf(1-d.z,4)*at;}v i=O+d*t;v ld=!(v(9+R(),9+R(),16)+i*-1);v hv=d+n*(n%d*-2);F lf=ld%n;if(lf<0||TR(i,ld,t,n)){lf=0;}if(mt&1){i=i*.2f;v c=((I)(ceilf(i.x)+ceilf(i.y))&1?v(3,1,1):v(3,3,3))*(lf*.2f+.1f);return pc+c*at;}F h=__powf(ld%hv*(lf>0),99);pc=pc+v(h,h,h)*at;O=i;d=hv;}return pc;}__global__ void GetColor(unsigned char *im){I x=blockIdx.x;I y=threadIdx.x;v cd=!v(-6,-16,0);v cu=!(v(0,0,1)^cd)*.002f;v cr=!(cd^cu)*.002f;v e=(cu+cr)*-256+cd;v h(13,13,13);for(I r=64;r--;){v de=cu*(R()-.5)*99.0f+cr*(R()-0.5f)*99.0f;h=Sample(v(17,16,8)+de,!(de*-1+(cu*(R()+x)+cr*(y+R())+e)*16.0f))*3.5f+h;}im[512*y*3+x*3+0]=h.x;im[512*y*3+x*3+1]=h.y;im[512*y*3+x*3+2]=h.z;}
#define gpuErrchk(a){gpuAssert((a),__FILE__,__LINE__);}
inline void gpuAssert(cudaError_t co,const char *f,I li,bool abort=true){if(co != cudaSuccess){fprintf(stderr,"GPUassert: %s %s %d\n",cudaGetErrorString(co),f,li);if(abort)exit(co);}}I main(){cudaSetDevice(0);DWORD st=GetTickCount();char* m=new char[512*512*3];unsigned char *db;cudaMalloc((void**)&db,(512*512*3));gpuErrchk(cudaPeekAtLastError());GetColor<<<512,512>>>(db);gpuErrchk(cudaDeviceSynchronize());cudaMemcpy(m,db,(512*512*3),cudaMemcpyDeviceToHost);gpuErrchk(cudaPeekAtLastError());DWORD em=GetTickCount()-st;fprintf(stderr,"Time: %dms\n",em);char* c=m;for(I y=512;y--;){for(I x=512;x--;){c=&m[y*512*3+x*3];c+=3;}}delete m;return 0;}