#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <math.h>

#define NZERI 10000
#define BINS 20

void ordina_zeri(double zeri[2]);
double calc_err(double zeri_th[2], double zeri_num[2]);
void solve_quadratic(double c[2], double zeri[2]);
void solve_quadratic_best(double c[2], double zeri[2]);
void updisto(int isto[BINS], double epsilon);

int main() {
    srand48(time(NULL));
    int i;
    int istoA[BINS] = {0}, istoB[BINS] = {0};
    double zeri_th[2], zeri_num[2], zeri_num_best[2];
    double c[2];
    double epsilon, epsilon_best;
    
    for(i = 0; i < NZERI; i++) {
        zeri_th[0] = 1e30 + (drand48() - 0.5) * 1e20;
        zeri_th[1] = 1e-30 + drand48() * 1e20;
        c[0] = zeri_th[0] * zeri_th[1];
        c[1] = -(zeri_th[0] + zeri_th[1]);
        solve_quadratic(c, zeri_num);
        solve_quadratic_best(c, zeri_num_best);
        
        ordina_zeri(zeri_th);
        ordina_zeri(zeri_num);
        ordina_zeri(zeri_num_best);
        
        epsilon = calc_err(zeri_th, zeri_num);
        updisto(istoA, epsilon);
        epsilon_best = calc_err(zeri_th, zeri_num_best);
        updisto(istoB, epsilon_best);
        if(i < 10) {
             printf("%.15G %.15G %.15G %.15G\n", zeri_th[0], zeri_th[1], epsilon, epsilon_best);
        }
    }
    
    FILE *out = fopen("isto.dat", "w");
    
    for(i = 0; i < BINS; i++) {
        fprintf(out, "%d %lf %lf\n", i, istoA[i] / (double) NZERI, istoB[i] / (double) NZERI);
    }
    
    fclose(out);
    
    return 0;
}

void ordina_zeri(double zeri[2]) {
    if(zeri[0] > zeri[1]) {
        double tmp = zeri[0];
        zeri[0] = zeri[1];
        zeri[1] = tmp;
    }
}

double calc_err(double zeri_th[2], double zeri_num[2]) {
    return fabs((zeri_th[0] - zeri_num[0]) / zeri_th[0]) + fabs((zeri_th[1] - zeri_num[1]) / zeri_th[1]);
}

void solve_quadratic(double c[2], double zeri[2]) {
    double delta = c[1] * c[1] - 4 * c[0];
    zeri[0] = (-c[1] + delta) / 2.0;
    zeri[1] = (-c[1] - delta) / 2.0;
}

void solve_quadratic_best(double c[2], double zeri[2]) {
    double delta = c[1] * c[1] - 4 * c[0];
    if(c[1] >= 0) {
        zeri[0] = -(c[1] + sqrt(delta)) / 2.0;
    }
    else {
        zeri[0] = (-c[1] + sqrt(delta)) / 2.0;
    }
    if(zeri[0] == 0.0) {
        zeri[1] = 0.0;
    }
    else {
        zeri[1] = c[0] / zeri[0];
    }
}

void updisto(int isto[BINS], double epsilon) {
    int k = floor(log10(epsilon) + BINS);
    
    if(k >= 0 && k < BINS) {
        isto[k]++;
    }
}