浮動小数点演算のベンチマーク

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

#define PI 3.14159265358979323846264

int cubic_equation( double a, double b, double c, double d, double roots[] );

int main(){

    double a, b, c, d, roots[ 3 ];
    int i, n;
    DWORD started, terminated;

    started = GetCurrentTime(); 
    for ( i = 1; i <= 1000000; i++ ){
      a = rand(); b = rand(); c = rand(); d = rand();
      if ( a != 0 ) n = cubic_equation( a, b, c, d, roots );
    }
    terminated = GetCurrentTime();
    printf( " total %10.3f seconds\n", ( terminated - started ) * 1e-3 );
    return 0;
}

double cubic_root( double x ){

    double root, prev_root, eps;
    int converged, n;
    
    if ( x == 0 ) return 0;
    else {
       converged = 0;
       root = 1e0;
       n = 1;
       do {
          prev_root = root;
          root -= ( root * root * root - x ) / ( 3 * root * root );
          eps = fabs( root / prev_root - 1e0 );
          if ( ( eps < 1e-8 ) || ( fabs( root ) < 1e-8 ) ) converged = 1;
          n++;
       } while ( converged == 0 );
       return root;
    }
}

int cubic_equation( double a, double b, double c, double d, double roots[] ){

    double p, q, t, a3, b3;
    int num_realroots;

    b /= ( 3 * a );  c /= a;  d /= a;
    p = b * b - c / 3;
    q = ( b * ( c - 2 * b * b ) - d ) / 2;
    a = q * q - p * p * p;

    if ( a > 0 ){
        num_realroots = 1;
        if ( q > 0 ) a3 = cubic_root( q + sqrt( a ) );
        else         a3 = cubic_root( q - sqrt( a ) );
        b3 = p / a3;
        roots[ 0 ] = a3 + b3 - b;
        roots[ 1 ] = -0.5 * ( a3 + b3 ) - b;
        roots[ 2 ] = fabs( a3 - b3 ) * sqrt( 3e0 ) / 2;
    }
    else {
        num_realroots = 3;
        a = sqrt( p );
        t = acos( q / ( p * a ) );
        a *= 2;
        roots[ 0 ] = a * cos( t / 3 ) - b;
        roots[ 1 ] = a * cos( ( t + 2 * PI ) / 3 ) - b;
        roots[ 2 ] = a * cos( ( t + 4 * PI ) / 3 ) - b;
    }
    return num_realroots;
}