Friday, 13 October 2017


-Safety, speed, usability, etc. all generally speaking originate from the same source, the system doing as little as possible in the simplest and understandable way possible while achieving the goal of doing a lot of the same thing it's good at.

-Transparency: the system must be transparent and open. It must not do anything behind the users back. The user must be able to view all that is going on in their system and understand it.

-Non-redundancy: in UI there ought to be only one way of doing things so as not to confuse the user or waste their time by making them learn many ways of doing the same thing.

-Limited abstraction: working close to the hardware while maintaining a level of readability and functionality. Unnecessary abstraction leads to waste of time, less speed, worse security and confusion.

-Stay out of the way of the user, never popup anything unless requested, never run anything in the background unless specifically requested.

Visit Russia, find subs, observe Militsiya, get a cold (after two years of health).
WM design ideas:

Applications are allowed to bring themselves on top only once (when they are started within certain time). After that the user must ALT+TAB to bring them on top.

OS/Compiler stuff:

1) Compile a static busybox and copy the files to a HDD image (after making the partition):
mount -o offset=32256 sda.bin tmp
2) mknod a few devices
3) Compile linux kernel
Tools these days...
4) qemu-system-x86_64 -m 512 -kernel bzImage -append "root=/dev/sda1" -enable-kvm -hda sda.bin

5) Modify tcc so it can be compiled into static binary, improve the assembler by adding syscall, use statically linked uClibc for startup code (plan: write the startup code entirely from scratch).

x86_64-asm.h: DEF_ASM_OP0(syscall, 0x0f05)
tccpp.c: get rid of ldexp and make your own d = d*((1<<((int)(exp_val-frac_bits))));
tccelf.c: get rid of all things related to libtcc1.a and crt?.o
tcc.c: add missing (dummy) function void __assert_fail (const char *assertion, const char *file, unsigned int line, const char *function) { abort(); }

tcc -DCONFIG_TCC_STATIC -o tcc.o -c tcc.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o libtcc.o -c libtcc.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o tccpp.o -c tccpp.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o tccgen.o -c tccgen.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o tccelf.o -c tccelf.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o tccasm.o -c tccasm.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o tccrun.o -c tccrun.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o x86_64-gen.o -c x86_64-gen.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o i386-asm.o -c i386-asm.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o libtcc1.o -c libtcc1.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o alloca86_64.o -c alloca86_64.S -DTCC_TARGET_X86_64 -I./include -static

gcc -nostdlib -static -o tcc libtcc1.o alloca86_64.o tcc.o libtcc.o tccpp.o tccgen.o tccelf.o tccasm.o tccrun.o x86_64-gen.o i386-asm.o ../uClibc-0.9.33/lib/crt1.o ../uClibc-0.9.33/lib/crti.o ../uClibc-0.9.33/lib/crtn.o ../uClibc-0.9.33/lib/libc.a
strip tcc

6) End up with standalone 290 kB static C-compiler/assembler.

Boots in a second.

void* syscall(void* syscall_number, void* param1, void* param2, void* param3) {
  __asm__ __volatile__(
    "mov %0, %%rax\n"
    "mov %1, %%rdi\n"
    "mov %2, %%rsi\n"
    "mov %3, %%rdx\n"
    : "g" (syscall_number), "g" (param1), "g" (param2), "g" (param3)
    : "rax", "rdi", "rsi", "rdx");

typedef unsigned long int uintptr;
typedef long int intptr;

static intptr write(int fd, void const* data, uintptr nbytes) {
  syscall((void*)1, (void*)(intptr)fd, (void*)data, (void*)nbytes);

int main(int argc, char* argv[]) {
    write(1, "Hello, World!\n", 15);

    return 0;

_start() {
  __asm__ __volatile__(
    "call main\n"
    "mov $60, %rax\n"
    "mov 0, %rdi\n"

7) Compile hello world on a few MB linux installation on a 64 bit CPU.

8) Port my fluid/Navier-Stokes simulator into javascript while you're at it, eliminate need for fft by short relaxation method and simply write del, del2 and other vector and matrix operations into for loops. Much faster with SDL, but pretty much realtime in javascript as well. Also runs fine in my cellphone.

<script type="text/javascript">
window.onload = function() {
  var  t = document.getElementById("a").getContext("2d");
  var  g = t.createImageData(128, 128);
  var  N = 128;
  var  f = new Array(N*N);
  var  p = new Array(N*N);
  var  u = new Array(N*N);
  var  v = new Array(N*N);
  var ux = new Array(N*N);
  var vx = new Array(N*N);
  var uy = new Array(N*N);
  var vy = new Array(N*N);
  var u2 = new Array(N*N);
  var v2 = new Array(N*N);
  var us = new Array(N*N);
  var vs = new Array(N*N);
  var px = new Array(N*N);
  var py = new Array(N*N);

  function r(tframe) {
    var a, b, x, y;

    for(y=-10; y<10; y++)
      for(x=-10; x<10; x++)
        u[0.75*N+x+(N/2+y)*N] = 0.15;

    for(y=1; y<N-1; y++)
      for(x=1; x<N-1; x++) {
        ux[x+y*N] = 0.50*(u[(x+1)+y*N] -                u[(x-1)+y*N]);
        vx[x+y*N] = 0.50*(v[(x+1)+y*N] -                v[(x-1)+y*N]);
        uy[x+y*N] = 0.50*(u[x+(y+1)*N] -                u[x+(y-1)*N]);
        vy[x+y*N] = 0.50*(v[x+(y+1)*N] -                v[x+(y-1)*N]);
        u2[x+y*N] = 0.25*(u[(x+1)+y*N] - 2.0*u[x+y*N] + u[(x-1)+y*N]);
        v2[x+y*N] = 0.25*(v[(x+1)+y*N] - 2.0*v[x+y*N] + v[(x-1)+y*N]);
        u2[x+y*N]+= 0.25*(u[x+(y+1)*N] - 2.0*u[x+y*N] + u[x+(y-1)*N]);
        v2[x+y*N]+= 0.25*(v[x+(y+1)*N] - 2.0*v[x+y*N] + v[x+(y-1)*N]);

    for(x=0; x<N*N; x++) {
      us[x] = u[x] + u[x]*ux[x] + v[x]*uy[x] + u2[x]/20.0 - px[x]/20.0 + 0.00000;
      vs[x] = v[x] + u[x]*vx[x] + v[x]*vy[x] + v2[x]/20.0 - py[x]/20.0 - 0.00000;

    for(y=1; y<N-1; y++)
      for(x=1; x<N-1; x++) {
        f[x+y*N] = 0.5*(us[(x+1)+y*N] - us[(x-1)+y*N] + vs[x+(y+1)*N] - vs[x+(y-1)*N]);
        p[x+y*N] = 0.25*(p[(x+1)+y*N] + p[x+(y+1)*N] + p[(x-1)+y*N] + p[x+(y-1)*N] - f[x+y*N]);

    for(y=1; y<N-1; y++)
      for(x=1; x<N-1; x++) {
        px[x+y*N] = 0.5*(p[(x+1)+y*N] - p[(x-1)+y*N]);
        py[x+y*N] = 0.5*(p[x+(y+1)*N] - p[x+(y-1)*N]);
         u[x+y*N] = us[x+y*N] - px[x+y*N];
         v[x+y*N] = vs[x+y*N] - py[x+y*N];

    for(x=0; x<N*N; x++) {
      a = 4*x;
      b = 1500*Math.sqrt(u[x]*u[x] + v[x]*v[x]);
      if(b>255) b = 255;[a] = b;[a+1] = b;[a+2] = b;[a+3] = 255;
    t.putImageData(g, 0, 0);
  for(x=0; x<128*128; x++) {
     f[x] = 0;
     p[x] = 0;
     u[x] = 0;
     v[x] = 0;
    ux[x] = 0;
    uy[x] = 0;
    vx[x] = 0;
    vy[x] = 0;
    u2[x] = 0;
    v2[x] = 0;
    us[x] = 0;
    vs[x] = 0;
    px[x] = 0;
    py[x] = 0;

<canvas id="a" width="128" height="128" style=width:12cm;height:12cm></canvas>

Monday, 2 October 2017

Belonging, purpose, transcendence and storytelling...

Navier-Stokes, good luck making the code shorter...

N  = 192; dt = 0.25; p = zeros(N, N);
u = zeros(N, N); v = zeros(N, N); ustar = zeros(N, N); vstar = zeros(N, N);
k = pi*[0:(N/2-1) (-N/2):(-1)]; [x y]  = meshgrid(k, k); ds = -(x.^2 + y.^2);
ds(1,1) = 1; px = 0.0; py = 0.0; gx = 0.0; gy = 0.0;
for t = 0:1000
    % constant velocity zone
    u(N/2-10:N/2+10, N/2-10:N/2+10) = 0.2;

    % solve intermediate velocity U*
    [ux uy] = gradient(u);
    [vx vy] = gradient(v);
    ustar = u + (u.*ux + v.*uy + del2(u)/20 - px/20 + gy)*dt;
    vstar = v + (u.*vx + v.*vy + del2(v)/20 - py/20 + gx)*dt;

    % solve Poisson equation: del2(p) = f
    fhat = fft2(divergence(ustar, vstar));
    p = N^2*real(ifft2(fhat./ds));

    % enforce del(U) = 0
    [px py] = gradient(p);
    u = ustar - px*dt;
    v = vstar - py*dt;

    cla; surface(sqrt(u.^2 + v.^2), 'edgecolor', 'none'); drawnow;

Espoo is full of furry animals.
Just a basic wave equation...

u  = zeros(192, 192);
du = zeros(192, 192);
[x y] = meshgrid(linspace(-1, 1, 192), linspace(-1, 1, 192));
for t = 0:0.002:0.4
    u = u + exp(-1000*((x-0.7*sin(7*t)).^2+(y-0.7*sin(11*t)-0.2).^2));
    u = u + 0.5*exp(-1000*((x+0.7*sin(5*t)).^2+(y+0.7).^2))*sin(3*31*t);

    du = du + del2(u);   
    u = u + du;
    cla; imagesc(u); drawnow;


Monday, 11 September 2017

I think I might see it slightly differently

I'm not convinced human minds truly understand the meaning of the questions asked or even the meaning of the answers we produce. Perhaps we simply react in a manner which minimizes the error function and this has come to be as a result of our kind of system adapting to the kinds of environments we've evolved in. Maybe the meaning we experience is at least partially illusory and related to emotional mechanisms which in turn are related survival and such.

It seems that for a consciousness memories can be material as much a computer to its user so at least that part of the human being can be identical to a Philosophical zombie. We can lose our computers and still remain conscious even if without our computers we are in some ways lesser begins. However, if we keep continuously losing all our memories, are we ourselves anymore from one moment to another even if some kind of consciousness remains?

Basically, there is no scientific interpretation of the Schrödinger equation or Einstein equation. Science is just the best prediction, interpretations are philosophy. Though, in the end both philosophy and physics are just practical human endeavours which only deal in conditionals.

Objects in relativity always "move" at the "speed" of light. Under normal circumstances, while they are stationary (relative to the observer), their "speed" is fully in the "direction of time". If they begin to move (in the normal sense), part of their available "speed" is "spent in spatial direction" and less remains for the "temporal direction". Similarly, gravity curves spacetime and the aforementioned "direction", so the objects direction of "movement" results in relative time dilation.

I think it is safe to say that nobody knows what the nature of the universe was prior to big bang, but it's plausible that time loses meaning between Planck time and zero. So talking about time before the big bang is likely nonsensical. Some other (speculative) dimensions relevant to certain concerns might still exist, but they wouldn't have the same meaning to us as time does. Time doesn't exist as the kind of independent dimension our intuition makes it appear anyway so it's not so far fetched that the meaning of time might vanish entirely in the beginning.

Sunday, 23 July 2017

Many worlds interpretation and the hard problem of consciousness

According to the Everett many-worlds interpreration of quantum mechanics the huge dimensionality of the many-particle wave function implies existence of parallel universes in which all possibilities occur simultaneously. While this interpretation might appear to be at odds with Occam's razor it has far fewer postulates than its competitors. It is also true that the information related to alternative histories in some form or another is currently required to predict experimental observations.

Time evolution of the wave function in general is always deterministic. Thus when we ask questions such as why we experience one future instead of another the typical answer is that there is nothing special about any particular world since all possible worlds exist and we ask the same question in all of them (at least in all relevant to this question). Quantum jumps are a consequence of short timescale decoherence, but decoherence itself is also governed by deterministic evolution of the wave function.
One can imagine a machine that doubles a person in a way where one half of their brain is original and the other copied. Both of these now full brains (with half copied and half original) are going to ask why they are the one they are instead of the other. For them the outcome appears random, but from the point of view of the operator (of the copy machine) there is nothing random about it at all. This copy process can be considered as branching of the multiverse and it corresponds to natural evolution of the wave function and it is not in a relevant way subject to the no-cloning theorem.

The experience at least in some sense is deterministic even for the consciousness, it's just that something violent is continuously happening without its knowledge, but if the consciousness knew everything there is to know about the world, it could predict what's going to happen to them as much as a split brain patient or anyone about to experience a qualia they've never experienced before can ever be expected to be able to understand what's in store for them (ignoring the fact that they can't know everything there is to know about the world).

The practical expectation values are among other things a result of the number of universes in which particular outcomes happen. One has to calculate the sum of all possible wave functions (from all the universes) to find a priori probability distributions, but that's another story. Anyway, such a theory allows one to make physically meaningful predictions which are accurate and highly consistent with observations.

However, while for an objective observer there is no particular problem as one branch of the multiverse is as real as any other, an individual conscious observer will likely be quite unsatisfied by the lack of practical deterministic predictions, but even if there was no problem from the point of view of the physicist, for a philosopher the hard problem of consciousness remains.

The hard problem of consciousness is basically asking what if any is the difference between a person who experiences consciousness and a system which simply reacts identically to a conscious person, but doesn't truly experience anything? Not everyone believes there is any problem associated with the question, but one has to consider that according to general relativity time isn't unique or fundamental in a way that would appear necessary for the experience of present, i.e. there is no passage of time and no unique present exists, rather time may simply be a kind of dimension no different from space. How then could consciousness ever experience the passage of time as we appear to experience it? The question becomes what sets apart the present moment from the past and the future for the consciousness when physically there doesn't appear to be anything that could do so?

It has been suggested that time is an entanglement phenomenon and arises from quantum mechanics. Such a time is meaningful only in relation to some other quantum variables often considered to be of thermodynamical origin. Thermodynamics itself appears nothing more than unlikely (low entropy) microstates evolving (chaotically) into likely (high entropy) microstates (that greatly outnumber low entropy ones). The standard problem of time is asking why the universe appears to have started in a highly unlikely low entropy state (perhaps anthropic), but that's separate from the problem of time in relation to consciousness. In a static set time evolution is simply a direction. How could this give rise to a meaningful present moment necessary for consciousness is at best unclear.

Friday, 16 June 2017

Why does the universe exists?

Certain interpretation of Gödel's incompleteness theorems states that not all truths associated with particular systems can be seen from within the system (this is basically a consequence of liars paradox). Does this suggest that there must always be something more than simply the system itself? Maybe not, it may simply be that not all truths can be known even if they're there. However, if it had to be the case, then perhaps it would be necessary that there will always be more complexity "somewhere out there". Then maybe we could claim that it is necessary that the cosmos (all of existence) is necessarily infinitely complex. Perhaps knowing that the cosmos is infinitely complex, is one of those things that cannot be known. Never the less, if that were the case, it could solve a lot of foundational problems, because in such a cosmos nothing is fiction. We would expect to find anything somewhere within that set, even something like the birth of our universe and emergence of human consciousness.

"Once you eliminate the impossible, whatever remains, no matter how improbable, must be the truth."

Compared to the previous, questions asked by your everyday physicist seem almost trivial. We can always build scientific models and fit all available observations from our local neighborhood of the cosmos, the observable universe. Sooner or later that is. These models will explain the universe within reason and give rational and justified expectations of how it likely evolved in the past and how it can be expected to evolve in the future. These models are by definition our best guess and they will to the best of our ability serve the goals we have given them.

There remains, however, one deeply mysterious question. What is the nature of consciousness? Since we are the observer, the one who experiences finite rate for the passage of time among other things (generally speaking physics is tenseless as a consequence of general relativity), it is nontrivial to arrive at an understanding of how human consciousness could come to be. We can imagine a machine that externally speaking reacts to stimuli as we do, but it is less obvious what (if any) separates us from such a machine. If consciousness is nothing more than states of information evolving in time at a finite rate (that can be counted relative to some other states such as global entropy), it still isn't entirely obvious how our experience of the present can exist. Is it somehow complementary to something in the universe like certain aspects of quantum mechanics? Or is that too one of those questions that simply cannot be answered from within?

At least it seems quite likely that simply by constructing certain configuration of matter, the matter will acquire consciousness like we do. This is mind the problem might not be so severe. The question then becomes is there something special about the universe that is related to consciousness in some special way? Like nonlocality of quantum mechanics or perhaps the many worlds interpretation? Suppose there isn't anything particularly special about any of this, then in fact the most important postulate waiting for proof (to the extent anything can be proven) might be that the cosmos is necessarily infinitely complex. Given that, everything more or less follows.

The determinism of many worlds interpretation would be less problematic if there was no consciousness, because then there is nothing special about one path over another. What is, is the whole of multiverse. It's consciousness that experiences only one path and one present that is mysterious, because we would like to understand what determines we experience this path instead of some other. Perhaps one could inject free will and consciousness somewhere in there somehow, or not.

Saturday, 10 June 2017

What determines which universe you're going to experience?

1) Munchhausen's trilemma: Nothing can be proven unconditionally (as far as we know).

2) There is at least one consciousness that experiences things change at a finite rate (mine). It's probable that others exist as well.

3) Quantum mechanics is incomplete, because experimental results are best predicted by probabilistic models and probability is related to incomplete information.

4) General relativity implies that space and time are not independent of each other, fundamental or independent of mass-energy or in general information and events related to them.

5) Entropy grows with time.

6) Science cannot explain why the cosmos exists, only how the observable universe evolves. It tells us what future bets are justified based on the past. This is likely the best anyone can ever hope to do.

7) Purpose of postulates is to add to the predictive power of our existing models of reality. Supernatural is not a useful postulate as it does not make predictions that could be tested. Supernatural is not nonsense because it's probably false, it's nonsense because the postulate cannot be justified.

8) Fine-tuning of the observable universe either doesn't exist or the cosmos is infinitely complex and we must by necessity find ourselves somewhere within such set.

9) The observable universe probably hasn't always existed. The cosmos probably has, to the extent that it even has time. The origin of the observable universe is likely atemporal and thus talk of a cause is likely nonsense.

Thursday, 25 May 2017

How to win 1 M€ with 99% probability?

Play no limit roulette and bet 1 M€ on either black or red and always double your bet until you've made 1 M€ profit (then stop and never look back).

The joke is that you need to be able to go 127 M€ into debt to have 99% probability of winning 1 M€ as you need to be prepared to double in case you lose to win back your money and reach your goal. Your expectation value is always negative as you are likely to lose quite a lot when you do, but if you happen to have access to 127 M€, your chances of winning 1 M€ are approx. 99%.

A rational agent should never bet on negative expectation value, but given 99% rate of success, how many of us wouldn't? Besides of course banks and such to whom 127 M€ is nothing and only cumulative long term profits matter.

European roulette has slightly better odds than American since European has only one 0 instead of two (0 and 00). I wonder what's the limit in Monte Carlo?

If you have only 7 M€ limit (and start with 1 M€ bet), your chances of winning 1 M€ are still around 86%. If you have 100 k€ debt limit and start with 100 k€ bet, your chances of winning 1 M€ are 5.6%.

Beware the Barrenness of a Busy Life

clear all
MAXDEBT = 127;
for x = 1:1e6                       % test the strategy
    money = 0;                      % start with no money
    bet = 1;                        % starting bet
    while money < 1                 % play until you made 1 unit
        money = money - bet;        % bet
        if rand < 18/37             % roulette (18 red, 18 black and 0)
            money = money + 2*bet;  % if you win, you win double your bet
        end                         % otherwise lost the bet
        bet = 2*bet;                % double the bet after each round
        if money-bet<-MAXDEBT       % stop if at the limit
    mon(x) = money;

Ex = mean(mon)                      % ~ -0.2
MIN = min(mon)                      % -127
Pwin = sum(mon>0)/x                 % ~ 0.99

Wednesday, 3 May 2017

The simplest acts of kindness are by far more powerful than a thousand heads bowing in prayer.

Our understanding of the arbitrary biological origin of morality does not determine our willingness to act morally. Subscribing to moral nihilism does not pose any problem. We've evolved to become moral creatures, and that's just the way we are even if such is ultimately arbitrary. There's nothing wrong in killing on some cosmic scale. Right and wrong are not properties of the universe, they are properties of living feeling creatures. Things that are wrong to us, are so because that's the way our species is (for the most part). This is due to our common history and to the extent that we don't share some values then that is just an argument against objective morality. This does not mean we should stop feeling, teaching of morality or punishing the guilty (or at least managing crime and human behavior). It just means morality is a functional state of our species (and to an extent our ecosystem).

It is not inconceivable that a species might exist that would accept murder (or something like that) as part of their daily lives (though it would likely need to exist under circumstances that would not be harmful or would even be beneficial for the survival of that species). Who are we to say it's wrong if every single individual of that species thinks it's right and it serves their survival.

We are not such a species and to kill another human is almost universally considered unacceptable and harmful by humans (for practical good reasons in my opinion). It is also a matter of how such things make us feel. Even if feelings are arbitrary, they are a property we possess and they are closely related to our conscious experience. They are also what makes us individuals and in my opinion we should proudly embrace our individuality, because what else is life and consciousness besides embracing ones own complexities and experience (which to some degree includes that which we call the external world)?

I almost always value the autonomy of any creature over their own lives above anything else. Some people possess values that may conflict mine. When it comes down to the basics, there is not much more to it. A particular value in vacuum is not irrational as such, but often poor reasons are given. Rationalization in general is not difficult and people love to do it, but what is less often discussed is that reasons aren't necessarily needed. Having them and thinking that's all there is might in fact make things worse. People try to rationalize arbitrary emotional values when in fact doing so and being ignorant of the arbitrary emotional basis is a sign of irrationality. Rationalization matters only to the extent that conflicting values are held. Sometimes, however, two people might have conflicting values without holding any conflicting values of their own. One might think this would be the only reason for war. To defend ones way of life against someone whose highest priority is to eliminate your way of life and replace it with something you can never accept. However, from what I know that has never been the primary reason for war. Go figure.

Values become interesting (or even problematic) only after you find the continued existence of your way of life is contributing harm to others or in general you find that your way of life is in conflict with your values. You are forced to weigh the value of your particular experience relative to how much you value that of others. You can try to do balancing, give up something that you don't value that much, but eventually you might reach funny conclusions such as that your particular experience or even existence isn't justified even on the scale of your own species, let alone on some cosmic scale. Then again, you are too small to know this to be true, but as you have nothing but your own faculties of reasoning to arrive at a conclusion, you are forced to arrive at one anyway. Unfortunately after you follow all the lines of reasoning you might still find it's all irrelevant and arbitrary. Sometimes all that matters is that which feels right. Thank goodness we appear to have something resembling the idea of conscience and intrinsic will to live, at least some if not most of us.

Monday, 1 May 2017

What if existence of suffering is necessary for consciousness to exist?

Fractals on Arduino Nano with fixed-point arithmetic today...

These small 128x64 monochrome OLED displays were 3.2 euros from ebay.

Serial output from the board as seen by minicom in linux.
Scientific knowledge is that which is useful for making predictions, everything else is just something else.
This 16 MHz board costs 2.2 euros in ebay.
/* avr-gcc -Os mandelbrot.c -mmcu=atmega328p -o program.elf
 * avr-objcopy -O ihex -R .eeprom program.elf program.hex
 * avrdude -c arduino -b 57600 -P /dev/ttyUSB0 -p atmega328p -vv -U flash:w:program.hex
 * (apt-get install avrdude gcc-avr avr-libc)

#define F_CPU 16000000UL
#define BAUD 115200

#include <avr/io.h>
#include <stdio.h>
#include <util/delay.h>
#include <util/setbaud.h>

putch(char a) {
  loop_until_bit_is_set(UCSR0A, UDRE0);
  UDR0 = a;

int main() {
  int x, y, n, F = 1e3;
  long x0, x1, x2, y0, y1, y2;

  UCSR0A |= _BV(U2X0);
  while(1) {
    for(y=0; y<23; y++) {
      for(x=0; x<80; x++) { 
        x0 = F*6*x/80 - 7*F/2;
  y0 = F*6*y/23 - 3*F;
  x1 = 0;
  y1 = 0;
  n = 0;
  while(sqrt(x1*x1+y1*y1) < F*8 && ++n < 64) {
          x2 = x1*x1/F - y1*y1/F;
          y2 = 2*x1*y1/F;
          x1 = x2 + x0;
          y1 = y2 + y0;

With the OLED display...

int mandelbrot(long xx, long yy) {
  long x0, y0, x1, y1, x2, y2;
  int n, F = 1000;
  x0 = F*3*xx/128 - 7*F/2;
  y0 = F*3*yy/64 - 3*F/2;
  x1 = 0;
  y1 = 0;
  n = 0;
  while(sqrt(x1*x1+y1*y1) < F*8 && ++n < 11) {
    x2 = x1*x1/F - y1*y1/F;
    y2 = 2*x1*y1/F;
    x1 = x2 + x0;
    y1 = y2 + y0;
  return n;

void loop() {
  unsigned char x, y, z, c;
  int xx, yy;

  for(y=0; y<8; y++) {
    xx = 0;
    for(xx=0; xx<256; xx++) {
      z = 1;
      c = 0;
      for(x=0; x<8; x++) {
        yy = x+y*8;
        if(mandelbrot(xx, yy)>10)
          c = c + z;
        z = 2*z;



Sunday, 16 April 2017

I realized that there's a big difference between deciding to leave and knowing where to go.

So apparently I was bored so I wrote a program in C that numerically simulates time evolution of the Schrödinger equation and a wave packet tunneling through a potential barrier. Nothing unusual, but to make it a bit harder, I decided that I will also write my own basic functions such as exponential function, square root and sine using nothing but standard additions, subtractions and multiplication. I also didn't replicate any known algorithm (at least not intentionally), but instead just improvised something (inefficient).

Gif compression is a bit funny, but you get the point. Simulation done by a C program, data plotted in matlab.
// #include <math.h>
#include <stdio.h>

#define PI 3.14159265358979323846264338327950288419716939937508

struct mycomplex {

  double real;
  double imag;

/* e^x is a function that grows at a rate proportional to its value */

/* really inefficient, but there's a more efficient one with taylor series below */
double r_exp(double x) {
  double r = 1;
    for(double y=0; y<x; y+=0.000002)
      r = r + r*0.000002;
    for(double y=0; y>x; y-=0.000002)
      r = r - r*0.000002;
  return r;

/* my routine for natural logarithm as inverse of e^x similarly to sqrt */
double my_log(double x) {
  double low = 0;
  double high = x;
  double neg = 1;
  if(x<1) {
    x = 1/x;
    high = x;
    neg = -1;
  double value = (low + high)/2;

  for(int i=0; i<30; i++) {
    if(r_exp(value) < x)
      low = value;
      high = value;
    value = (low + high)/2;
  return neg*value;

/* x^y */
double my_pow(double x, double y) {
  return r_exp(y*my_log(x));

/* abs(x) */

double my_abs(double x) {
    return -x;
    return x;

/* my routine for sqrt */

double my_sqrt(double x) {
  double low = 1;
  double high = x;
  double root = (low + high)/2;
  if(x<1) {
    low = x;
    high = 1;
  for(int i=0; i<30; i++) {
    if(root*root < x)
      low = root;
      high = root;
    root = (low + high)/2;
  return root;


double my_sqrt(double x) {
  return my_pow(x, 0.5);

/* function for taking a modulo for floating point numbers */

double my_fmod(double in, double n) {
    n = -n;
    in += n;
  return in;

/* my routine for sine */

double my_sin(double a) {
  double l = 0, x0 = 0, y0 = 1, x = 0, y = 1, neg = 1.0;

  /* utilize symmetries to only calculate quarter */

  while(a<0) {
    a = a + 2*PI;
  a = my_fmod(a, 2*PI);
  if(a>PI) {
    neg = -neg;
    a = 2*PI-a;
    a = PI-a;

  /* integrate perimeter length until it corresponds to given angle */

  while(l<a) {
    x = x + 0.000005;
    y = my_sqrt(1-x*x);
    l = l + my_sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0));
    x0 = x;
    y0 = y;

  return neg*x;


/* integrate the perimeter length until x/y corresponds to r */

double my_atan(double r) {
  double l = 0, x0 = 0, y0 = 1, x = 0, y = 1;
  while(x/y<r) {
    x = x + 0.000002;
    y = my_sqrt(1-x*x);
    l = l + my_sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0));
    x0 = x;
    y0 = y;
  return l;

/* my routine for multiplying two complex numbers */

struct mycomplex my_mul(struct mycomplex a, struct mycomplex b) {
  struct mycomplex out;

  out.real = a.real * b.real - a.imag * b.imag;

  out.imag = a.imag * b.real + a.real * b.imag;

  return out;


/* my routine for taking exp(z), where z is complex valued vector */

struct mycomplex my_exp(double real, double imag) {
  struct mycomplex out;

  out.real = r_exp(real)*my_sin(imag+PI/2);

  out.imag = r_exp(real)*my_sin(imag);

  return out;


/* my routine for constructing linear real valued vector */

my_linspace(double a, double b, int n, double out[]) {
  double df = (b-a)/(n-1);
  double f = a;
  for(int i=0; i<n; i++) {
    out[i] = f;
    f = f + df;

main() {

  double k[800];
  struct mycomplex diff2[800];
  struct mycomplex tmp;
  struct mycomplex tmpb;
  struct mycomplex psi[800];
  struct mycomplex potential[800];

  /* just for testing of accuracy */

  printf("my_sqrt(2)=%f sqrt(2)=%f\n", my_sqrt(2), sqrt(2));
  printf("my_sin(pi/4)=%f sin(pi/4)=%f\n", my_sin(PI/4), sin(PI/4));
  printf("4*my_atan(1)=%f 4*atan(1)=%f\n", 4*my_atan(1), 4*atan(1));
  printf("my_exp(1)=%f exp(1)=%f\n", r_exp(1), exp(1));
  printf("my_exp(-1)=%f exp(-1)=%f\n", r_exp(-1), exp(-1));

  /* prepare a potential barrier */

  for(int i=0; i<800; i++) {
    potential[i].real = 0;
    potential[i].imag = 0;
  for(int i=390; i<410; i++)
    potential[i].real = 0.05;

  my_linspace(0, 50, 800, &k);

  /* prepare initial wave function */

  for(int i=0; i<800; i++) {
    psi[i].real = r_exp(-0.07*(k[i]-12)*(k[i]-12));
    tmp = my_exp(0, -900*k[i]);
    psi[i] = my_mul(psi[i], tmp);
    psi[i].real = psi[i].real/8.7;
    psi[i].imag = psi[i].imag/8.7;
    // printf("%f\n", i/800.0);

  /* evolve */

  for(int t=0; t<2000000; t++) {
    diff2[0].real = 0;
    diff2[0].imag = 0;
    for(int i=1; i<799; i++) {
      diff2[i].real = (psi[i+1].real-psi[i].real)-(psi[i].real-psi[i-1].real);
      diff2[i].imag = (psi[i+1].imag-psi[i].imag)-(psi[i].imag-psi[i-1].imag);
    diff2[799].real = 0;
    diff2[799].imag = 0;
    for(int i=0; i<800; i++) {
      tmp = my_mul(potential[i], psi[i]);
      tmpb.real = diff2[i].real - tmp.real;
      tmpb.imag = diff2[i].imag - tmp.imag;
      tmp.real = 0;
      tmp.imag = 0.001;
      tmp = my_mul(tmp, tmpb);
      psi[i].real = psi[i].real + tmp.real;
      psi[i].imag = psi[i].imag + tmp.imag;
    if(t%1000==1) {
      for(int i=0; i<800; i++)
printf("%f, %f\n", psi[i].real, psi[i].imag);


Actually, you only need complex valued exponential function and some iterative inversions to calculate pretty much anything (normal).

#include <math.h>
#include <stdio.h>
#include <complex.h>

#define PI 3.14159265358979323846264338327950288419716939937508

/* iteratively invert a growing function with a and b as limits */
double inv(double (*f)(), double x, double a, double b) {
  for(int i=0; i<20; i++)
    if(f((a+b)/2)<x) a = (a+b)/2;
    else b = (a+b)/2;
  return (a+b)/2;

/* integrate from definition: e^z = de^z/dz, z^0 = 1 */
double complex my_exp(double complex z) {
  double complex v = 1;

  for(double a=0; a<1; a+=1/1e5)
    v = v + v*(z/1e5);
  return v;

/* more efficient (taylor series), less intuitive e^z */
double complex alt_exp(double complex z) {
  double complex v = 1, w = 1;
  for(int k=1; k<20; k++)
    v += (w*=z/k);
  return v;

double real_my_exp(double x) { return creal(my_exp(x)); }

double my_log(double x) {
  if(x<1) return -my_log(1/x);
  return invert(real_my_exp, x, 0, x);

double my_sq(double x) { return x*x; }

double my_sqrt(double x) {
  if(x>1) return invert(my_sq, x, 1, x);
  return invert(my_sq, x, x, 1);

double my_sin(double x) { return cimag(my_exp(I*x)); }
double my_cos(double x) { return creal(my_exp(I*x)); }
double my_tan(double x) { return my_sin(x)/my_cos(x); }
double my_cot(double x) { return my_cos(x)/my_sin(x); }
double my_pow(double x, double y) { return creal(my_exp(y*my_log(x))); }
double my_sinh(double x) { return (creal(my_exp(x))-creal(my_exp(-x)))/2; }
double my_cosh(double x) { return (creal(my_exp(x))+creal(my_exp(-x)))/2; }
double my_tanh(double x) { return my_sinh(x)/my_cosh(x); }
double my_coth(double x) { return my_cosh(x)/my_sinh(x); }
double my_log10(double x) { return my_log(x)/my_log(10); }

double my_atan(double x) {
  if(x<0) return -my_atan(-x);
  return invert(my_tan, x, 0, 1.6); /* 1.6 is just some value above pi/2 */

double my_asin(double x) { return my_atan(x/my_sqrt(1-x*x)); }
double my_acos(double x) { return 2*my_atan(1)-my_atan(x/my_sqrt(1-x*x)); }

/* alternative sin by integration of perimeter until y corresponds to a*/
double alt_sin(double a) {
  double c, x = 1, y = 0, xn = 1, yn = 0, l = 0, dl;

  /* utilize symmetries */
  while(a<0) a = a + 2*PI;
  while(a>2*PI) a = a - 2*PI;
  if(a>PI) return -alt_sin(2*PI-a);
  if(a>PI/2) return alt_sin(PI-a);

  /* integrate the perimeter length */
  while(l<a) {
    x = x - 1e-5;
    y = sqrt(1-x*x);
    dl = sqrt((x-xn)*(x-xn)+(y-yn)*(y-yn));
    l = l + dl;
    xn = x;
    yn = y;
  return y;


/* alternative atan by integrating the perimeter length until x/y corresponds to v */
double alt_atan(double v) {
  double l = 0, x0 = 0, y0 = 1, x = 0, y = 1;
  if(v<0) return -alt_atan(-v);
  while(x/y<v) {
    x = x + 1e-5;
    y = sqrt(1-x*x);
    l = l + sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0));
    x0 = x;
    y0 = y;
  return l;

/* discrete Fourier transform: F(w) = sum(f(t)e^(-iwt)) */
int dft(double complex *a, double complex *b, int n) {
  for(int m=0; m<n; m++) {
    b[m] = 0;
    for(int t=0; t<n; t++) /* w = 2*PI*m/n */
      b[m] = b[m] + a[t]*cexp(-I*2*PI*m*t/n); /* (a+bi)*(c+di) */

/* more efficient, less intuitive Fourier transform */
int fft(double complex a[], double complex b[], int n, int step) {
  double complex t;
  if(step==1) for(int i=0; i<n; i++) b[i] = a[i];
  if(step<n) {
    fft(b, a, n, step * 2);
    fft(b + step, a + step, n, step * 2);
    for(int i=0; i<n; i+=2*step) {
      t = cexp(-I*PI*i/n)*a[i+step];
      b[i/2] = a[i]+t;
      b[(i+n)/2] = a[i]-t;
} /* notice that this fft alters the input vector */

int main() {
  double complex a[] = {1, 2, -1, 4};
  double complex b[4];

  printf("Function  \tMy        \tC       \tAlt\n");
  printf("exp(1)    \t%f\t%f\t%f\n", creal(my_exp(1)), exp(1), creal(alt_exp(1)));
  printf("exp(-1)   \t%f\t%f\t%f\n", creal(my_exp(-1)), exp(-1), creal(alt_exp(-1)));
  printf("sqrt(2)   \t%f\t%f\n", my_sqrt(2), sqrt(2));
  printf("sqrt(0.5) \t%f\t%f\n", my_sqrt(0.5), sqrt(0.5));
  printf("log(2)    \t%f\t%f\n", my_log(2), log(2));
  printf("log(0.5)  \t%f\t%f\n", my_log(0.5), log(0.5));
  printf("atan(1)   \t%f\t%f\t%f\n", my_atan(1), atan(1), alt_atan(1));
  printf("atan(-1)  \t%f\t%f\t%f\n", my_atan(-1), atan(-1), alt_atan(-1));

  dft(a, b, 4);
  printf("Discrete Fourier Transform of 1 2 -1 4\n");
  for(int i=0; i<4; i++)
    printf("%f\t%f\n", creal(b[i]), cimag(b[i]));

  printf("Fast Fourier Transform of 1 2 -1 4\n");
  fft(a, b, 4, 1);
  for(int i=0; i<4; i++)
    printf("%f\t%f\n", creal(b[i]), cimag(b[i]));