Не правильно работают функции , Не могу разобраться, что не так...

Доброго времени суток!

Я тут просил помочь не так давно с распределением Эрланга, но мою тему перенесли в алгоритмы... Сейчас несколько другой случай, но продолжение того. Начитавшись по ссылкам, которые дал trainer и слабовато поняв (уж, простите - обозначения пугают), нашел с горем пополам книжку Илья Труб "Объектно-ориентированное моделирование на С++". В ней как раз даны примеры ГСЧ по Эрлангу, нормальному и равномерному распределениям. Попытался сгенерировать псевдослучайную последовательность. По равномерному (оно и понятно) генерирует, а по Эрлангу и нормальному какие-то чудеса. На Эрланге выскакивает исключение "Floating point overflow". На нормальном распределении тоже что-то не то. Всё, что добавлено и закоментировано в данный момент - это добавлено мной потом во время моих плясок вокруг ошибок. Результат был таким, что исключения пропали, а вывод-то неправильный. Можете проверить и помочь исправить ошибки? Весь код (в именах функций стоят префиксы, определяющие принадлежность к тому или иному распределению) выкладываю ниже. Также выкладываю ссылку на архив со сканами 7-ми страниц из книги, откуда я взял код (там с комментариями автора - посмотрите). Там же в архиве задание, но мне бы хотя бы распределениями разобраться помогли....

//---------------------------------------------------------------------------
 
#include <vcl.h>
#include <cstdlib.h>
#include <time.h>
#include <math.h>
 
#pragma hdrstop
 
#include "Unit1.h"
 
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
const int N = 10; // длина случайной последовательности для распр.
//---------------------------------------------------------------------------
// Функции, необходимые для вычисления распределений
//---------------------------------------------------------------------------
// Функции, необходимые для вычисления распределения Эрланга
float function_erl(float mu, int k, float t)
// Нахождение значений функции в точке t
{
  long double prod, s;
  prod=1.0; s=1.0;
    for (int i = 1; i < k; i++) {
        prod *= mu;
        prod *= t;
    }
  s=(long double)s+(long double)prod;
  s=s*exp(-mu*t);
  return (1-s);
}
 
/* float search_edge2_erl(float mu, float right, int k)
{
    float edge2, edge2_temp, value;
    edge2 = (float)k/mu+10*sqrt((float)k)/mu;
    value=function_erl(mu, k, edge2);
    while (value<right) {
    __try
    {
        value = function_erl(mu, k, value);
        if (value<right) value*=2;
    }
    catch( ... )
    {
      value = right;
    }
    }
    edge2  = value;
    return edge2;
}  */
 
float equ_erl(float mu, float right, int k, float eps)
// Вычисление функции распределения в точке t
{
  float edge1, edge2, middle, value;
  edge1 = 0.0;
  edge2 = (float)k/mu+10*sqrt((float)k)/mu;         // Эти две строки нужно закомментировать
  while (function_erl(mu,k,edge2)<right) edge2*=2;  // если раскомментировать следующую
 
//  edge2 = search_edge2_erl(mu, right, k); // функцию search_edge2_erl раскомментируйте выше
  while ((edge2-edge1)>eps) {
      middle=(edge1+edge2)/2;
      value=function_erl(mu, k, middle);
      if ((value-right)<0) edge1=middle;
      else edge2=middle;
  }
  return ((edge1+edge2)/2);
}
 
float get_erlang(float mu, int k, float eps)
// Рраспределение Эрланга
{
  int r_num=rand(); float root, right;
  right=(float)r_num/RAND_MAX+1;
  root = equ_erl(mu, right, k, eps);
  return (root);
}
//---------------------------------------------------------------------------
// Функции, необходимые для вычисления нормального распределения
float function_norm(float mean, float disp, float x)
//
{
  float result;
  result = (1.0/(disp*sqrt(2*M_PI)))*exp(-0.5*((x-mean)/disp)*((x-mean)/disp));
  return (result);
}
 
float simpson(float A, float B, float mean, float disp)
//
{
  float k1, k2, k3, s, x, h1, h;
  h = 0.01;
  s = 0; h1 = h/1.5;
  k1=function_norm(mean, disp, A);
  for (x=A; (x<B)&&((x+h-B)<=h1); x=x+h) {
     k2=function_norm(mean, disp, x+h/2);
     k3=function_norm(mean, disp, x+h);
     s=s+k1+4*k2+k3;
     k1=k3;
  }
  s=s*h/6;
  return (s);
}
 
float equ_norm(float bottom_bound, float top_bound, float mean, float disp, float almost_all, float eps, float right)
//
{
  float edge1, edge2, middle, cover, value;
  edge1 = bottom_bound;
  edge2 = top_bound;
  if (right>almost_all) return top_bound; else
  if (right<(1-almost_all)) return bottom_bound; else
  cover = 0;
  while ((edge2-edge1)>eps) {
      middle=(edge1+edge2)/2;
      value=simpson(edge1, middle, mean, disp);
      if ((cover+value-right)<0) {
         edge1=middle;
         cover = cover+value;
      }
      else edge2=middle;
  }
//  return ((edge1+edge2)/2);    // в книге этой строки нет
}
 
float get_norm(float mean, float disp, float eps)
// Нормальное распределение
{
  int r_num; float root, bottom_bound, top_bound, almost_all, right;
  bottom_bound=mean-disp*sqrt(-log(2*M_PI*eps*eps*disp*disp));
  top_bound=mean+disp*sqrt(-log(2*M_PI*eps*eps*disp*disp));
  almost_all=simpson(bottom_bound, top_bound, mean, disp);
  r_num=rand();
  right=(float)r_num/RAND_MAX+1;
  root = equ_norm(bottom_bound, top_bound, mean, disp, almost_all, eps, right);
  return (root);
}
 
//---------------------------------------------------------------------------
// Функции, необходимые для вычисления равномерного распределения
float get_uniform()
// Равномерное распределение
{
  int number = rand();
  float result = (float)number/(RAND_MAX+1);
  return (result);
}
//---------------------------------------------------------------------------
// Конец блока функций, необходимых для вычисления распределений
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
    : TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Button1Click(TObject *Sender)
{
  srand((unsigned)time(0));
  Memo1->Clear();
  for (int i = 0; i < N; i++) {
     Memo1->Lines->Add(FloatToStr(get_uniform()));
  }
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Button2Click(TObject *Sender)
{
  Memo2->Clear();
  srand((unsigned)time(0));
  float eps = 0.01;
  for (int i = 0; i < N; i++) {
 
     Memo2->Lines->Add(FloatToStr(get_erlang(90, 2, eps)));
  }
}
//---------------------------------------------------------------------------
 
void __fastcall TForm1::Button3Click(TObject *Sender)
{
  Memo3->Clear();
  srand((unsigned)time(0));
  float eps = 0.01;
  for (int i = 0; i < N; i++) {
    Memo3->Lines->Add(FloatToStr(get_norm(12.35, 24.7875, eps)));
  }
}
//---------------------------------------------------------------------------

Скачать архив со сканами из книги


Заранее благодарен!

Цитата =SAPSAN= @ 26.10.09, 05:53

По равномерному (оно и понятно) генерирует, а по Эрлангу и нормальному какие-то чудеса. На Эрланге выскакивает исключение "Floating point overflow".

А тебе для расчетов принципиально именно тип float использовать, а не double ?

Сообщение отредактировано: leo - 26.10.09, 08:28

Что double, что long double, что float - все равно исключение на очередной итерации в строке

  while (function_erl(mu,k,edge2)<right) edge2*=2;  // если раскомментировать следующую

Я всю башку уже сломал - что делать-то? Проверьте, пожалуйста, сами. Может я что-то не так вообще делаю.

Цитата =SAPSAN= @ 26.10.09, 14:48

что делать-то?

Свести сложные выражения к простым, отладчиком пройтись.

Цитата =SAPSAN= @ 26.10.09, 14:48

Что double, что long double, что float - все равно исключение на очередной итерации в строке

 while (function_erl(mu,k,edge2) < right) edge2*=2;  // если раскомментировать следующую

Похоже на зацикливание и как следствие переполнение edge2, т.к. у тебя right >= 1, а при mu=90, k=2, и начальном edge2 ~ 0.18 результат function_erl всегда <= 1

Сообщение отредактировано: leo - 26.10.09, 15:52

Мама дорогая! Ошибка с Эрлангом была в строчке функции get_erlang

right=(float)r_num/RAND_MAX+1;

Скобок-то нет...
А в функции equ_norm нужна строка

  return ((edge1+edge2)/2);    // в книге этой строки нет

помимо той же ошибки со скобками, что и в Эрланге.

Всем спасибо за ответы!

Сообщение отредактировано: =SAPSAN= - 26.10.09, 17:50

Здравствуйте уважаемый=SAPSAN=! На исходниках.ру наткнулся на Ваше сообщение про распределение Эрланга и нормальное распределение. У меня такая же проблема требуется реализовать пример про автоперевозки из книги Ильи Труба на с++билдер. Я всё написал но появляется сообщение об ошибке floating point overflow. Если Вы реализовали эти алгоритмы не могли бы Вы помочь мне?