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;
  if(x>0)
    for(double y=0; y<x; y+=0.000002)
      r = r + r*0.000002;
  if(x<0)
    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;
    else
      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) {
  if(x<0)
    return -x;
  else
    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;
    else
      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) {
  if(in<0)
    n = -n;
  while(in>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;
  }
  if(a>PI/2)
    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));
  exit(0);
  */

  /* 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));
  printf("\n");

  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("\n");

  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]));
}