Форум Точек.нет - общение без границ !

Ахинея

4.01.2010 - 14:21

Выложенные здесь алгоритмы преследуют исключительно учебные цели.
Код неоптимизирован, местами морально устарел и показывает только принцип решения той или иной задачи. Однако он вполне рабочий и может быть использован с соответствующей дорботкой под Ваши собственные нужды.

Буду рада, если не только я буду их пополнять)

Содержание:
Перевод из ... в ...
Сортировка массивов
Графика
Матрицы
Специальные функции
Методы поиска
Решение систем линейных уравнений
Решение нелинейных уравнений
Решение диф. уравнений

Будет еще дополнено.

Спасибо

Ахинея

4.01.2010 - 14:29

Градусы, минуты, секунды - в градусы:

» Кликните сюда для просмотра оффтоп текста.. «

Градусы - в градусы, минуты, секунды:

» Кликните сюда для просмотра оффтоп текста.. «

Комплексных величин - в экспоненциальную форму

» Кликните сюда для просмотра оффтоп текста.. «

Градусы по Фаренгейту - в градусы по Цельсию

» Кликните сюда для просмотра оффтоп текста.. «

Декартовы координаты - в полярные (2D)

» Кликните сюда для просмотра оффтоп текста.. «

Полярные - в декартовы

» Кликните сюда для просмотра оффтоп текста.. «

Декартовы координаты - в сферические (3D)

» Кликните сюда для просмотра оффтоп текста.. «

Сферические - в декартовы

» Кликните сюда для просмотра оффтоп текста.. «

Ахинея

4.01.2010 - 14:37

Сортировка массивов

Binary insertions

» Кликните сюда для просмотра оффтоп текста.. «

Straight insertion

» Кликните сюда для просмотра оффтоп текста.. «

Heap sort

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Heap sort
//  (c) Johna Smith, 1996
//
//  Method description:
//    This algorithm sorts array using pyramids. The pyramid is the
//    sequence of characters h(L), h(L+1),...,h(R) with the following
//    properties:
//          h(i)<=h(2i);  h(i)<=h(2i+1), i=L..R/2
//                        h1
//                     /     \               This is an example of
//                   h2       h3                   pyramid
//                 /    \   /    \
//               h4     h5 h6     h7
//              .....................
//    There are two steps in this algorythm:
//    1) Building pyramid from the given array using Floyd method of
//    building pyramid "on the same place".
//    2) Sorting built pyramid:
//       Swap last pyramid element (x) with the top one and shift x
//       to the right place using Floyd method
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

const N = 10;  // array size

int array[N]={100,15,234,11,63,78,4,200,0,4}; // array of N integers
int swp; // auxulary variable for swapping
int left,right; // indexes for left and right bounds of the sorted part of the array

void show_array(void) // this function displays array
{
  for (int i=0;i<N;i++)
    printf("%d ",array[i]);
}

// this function shifts element in the pyramid to the right place
// it helps to build pyramid "on the same place" (using Floyd method)

void shift(int left, int right)
{
  int i,j;
  int swp;
  int element;

  i=left;
  j=2*left+1;
  element=array[left];
  if (j<right && array[j]<array[j+1]) j++;
  while (j<=right && element<array[j])
  {
    // swapping elements
    swp=array[i];
    array[i]=array[j];
    array[j]=swp;
    // now i-th element is on j-th place:
    i=j;
    element=array[i];
    // and j-th element is lower (see pyramid picture)
    j=2*j+1;
    if (j<right && array[j]<array[j+1]) j++;
  }
}

void main(void)
{
   // Displaying unsorted array
   printf("Unsorted array: ");
   show_array();

   // Sorting
   left=N/2;
   right=N-1;
   // building pyramid from array
   while (left>0)
   {
     left--;
     shift(left,N-1);
   }
   // sorting built pyramid
   while (right>0)
   {
     // swapping elements
     swp=array[0];
     array[0]=array[right];
     array[right]=swp;
     // shifting element to the right place
     right--;
     shift(0,right);
   }

   // Displaying sorted array
   printf("\nSorted array: ");
   show_array();
}

Straight selections

» Кликните сюда для просмотра оффтоп текста.. «

Shaker sort

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Shaker sort
//  (c) Johna Smith, 1996
//
//  Method description:
//    Advanced bubble sort.
//    There are two main steps:
//     1) Move greatest element to the end, moving forward and
//        remembering place of the last elements exchange
//        this place will be right bound of the unsorted part
//     2) Move smallest element to the beginning, moving backward and
//        remembering place of the last elements exchange
//        this place will be left bound of the unsorted part
//    Repeat this two steps while right bound is greater than left
//    This method is useful when almost all elements in the array
//    are sorted.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

const N=     10;  // array size

int array[N]={100,15,234,11,63,78,4,200,0,4}; // array of N integers
int index; // auxulary variable (the place of last swap)
int swp; // auxulary variable for swapping
int left,right; // indexes for left and right bounds of the sorted part of the array

void show_array(void) // this function displays array
{
  for (int i=0;i<N;i++)
     printf("%d ",array[i]);
}

void main(void)
{
   // Displaying unsorted array
   printf("Unsorted array: ");
   show_array();

   // Sorting
   left=0;
   right=N-1;
   do
   {
     // moving smallest element to the beginning (moving backward)
     for (int i=right;i>left;i--)
     {
        if (array[i-1]>array[i])
        {
          // swapping a(i) and a(i+1)
          swp=array[i];
          array[i]=array[i-1];
          array[i-1]=swp;
          index=i; // save index
        }
     }
     left=index;
     // moving greatest element to the end (moving forward)
     for (i=left;i<right;i++)
     {
       if (array[i]>array[i+1])
       {
         // swapping a(i) and a(i+1)
         swp=array[i];
         array[i]=array[i+1];
         array[i+1]=swp;
         index=i; // save index
       }
     }
     right=index;
} while (left<right);

   // Displaying sorted array
   printf("\nSorted array: ");
   show_array();
}

Shell sort

» Кликните сюда для просмотра оффтоп текста.. «

Quick sort (recursive)

» Кликните сюда для просмотра оффтоп текста.. «

Quick sort (non-recursive)

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Quick sort (non-recursive)
//  (c) Johna Smith, 1996
//
//  Method description:
//    1) Split array into two parts and remember middle element
//    2) Scan left part for element greater than middle
//    3) Scan right part for element less than middle
//    4) Swap these elements
//    So we have an array where all left elements are less than right elements
//    Apply these four steps to each part (left and right) of the array
//    until we have parts that contain only one element.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

const N = 10;         // array size
const STACK_SIZE = 4; // stack size must be >=log N

int array[N]={100,15,234,11,63,78,4,200,0,4}; // array of N integers

void show_array(void) // this function displays array
{
  for (int i=0;i<N;i++)
    printf("%d ",array[i]);
}

struct {int left; int right;} stack[STACK_SIZE]; // stack emulation
unsigned int ss; // current stack size
int i,j;
int element; // auxulary variable for middle element in the interval
int swp; // auxulary variable for swapping
int left,right; // interval bounds

void main(void)
{
  // Displaying unsorted array
  printf("Unsorted array: ");
  show_array();

  // Sorting
  ss=1;
  stack[0].left=0;
  stack[0].right=N-1;
  do
  {
    // pushing from stack last request
    left=stack[ss-1].left;
    right=stack[ss-1].right;
    ss--;
    do
    {
      // splitting array[left]..array[right] interval
      i=left;  // index for left part
      j=right;  // index for right part
      element=array[(left+right)/2]; // middle element
      do
      {
        while (array[i]<element) i++; // scanning left part
        while (element<array[j]) j--; // scanning right part
        if (i<=j)
        {
          // swapping elements
          swp=array[i];
          array[i]=array[j];
          array[j]=swp;
          i++; j--;
        }
      } while (i<=j);
      // pushing request into stack
      // (we select longest part and push request for it to reduce stack size)
      if (j-left<right-i)
      {
        if (i<right) // push request for the right part (because it's longer than left)
        {
          ss++;
          stack[ss-1].left=i;
          stack[ss-1].right=right;
        }
        right=j; // continue sorting left part
      } else
      {
        if (left<j) // push request for the left part (because it's longer than right)
        {
          ss++;
          stack[ss-1].left=left;
          stack[ss-1].right=j;
        }
        left=i; // continue sorting right part
      }
    } while(left<right);
  } while(ss!=0);

  // Displaying sorted array
  printf("\nSorted array: ");
  show_array();
}

Exchanges

» Кликните сюда для просмотра оффтоп текста.. «

Ахинея

5.01.2010 - 13:51

Графика

Фрактал (листья папоротника)

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Leafs drawing
//  (c) Johna Smith, 1996
//
//  Method description:
//    We take a closed area, which will be main leaf. We want to generate
//  four fractal leafs from this one. To make it we need to scale, move and
//  rotate main leaf. These transformations described by the following
//  equations:
//          x'=ax+by+c, y'=dx+ey+f
//    So we need 6 coefficients for each leaf (24 coefficients at all).
//  These coefficients are calculated with special mathematical method
//  (here are already calculated coefficients)
//    To draw whole picture we need an iterational process:
//   1) Take any point within main leaf borders
//   2) Call function iterate for one of 4 leafs (its nuber we select
//      randomly but than the larger is leaf's area the higher is probability
//      that we call 'iterate' for this leaf
//   3) Then we take transformed point as (x,y) and repeat iterations
//
//  ! This program is NOT optimized for best performance !
//  To do so don't use putpixel function - change bytes in video memory directly.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <stdlib.h>
#include <dos.h>
#include <graphics.h>
#include <conio.h>

// this function initializes graphics mode
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void init_gr(void)
{
   /* request autodetection */
   int gdriver = DETECT, gmode, errorcode;

   /* initialize graphics mode */
   initgraph(&gdriver, &gmode, "");

   /* read result of initialization */
   errorcode = graphresult();

   if (errorcode != grOk)    /* an error occurred */
   {
      printf("Graphics error: %s\n", grapherrormsg(errorcode));
      printf("Press any key to halt:");
      getch();
      exit(1);               /* return with error code */
   }
}

// this function shuts graphics mode down
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void end_gr(void)
{
closegraph();
}

// this function puts pixel on the screen in (x,y) position using color 'color'
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void PutPixel(int x, int y, int color)
{
putpixel(x,y,color);
}

#define N       200000L // number of points

float x=0, y=0; // iteration variables (are global)

//                       a     b     c    d    e    f        probability
float coeff[4][6] = {{ 0.00, 0.00, 0.00,0.16,0.00,0.00,},   // 01%
                    { 0.85, 0.04,-0.04,0.85,0.00,1.60,},   // 85%
                    { 0.05,-0.26, 0.23,0.22,0.00,1.60,},   // 07%
                    {-0.15, 0.30, 0.26,0.24,0.00,0.44,},}; // 07%

int colors[5]={10,11,2,3,2};

void iterate(char i,int c)
{
float x1,y1;

x1=x*coeff[i][0]+y*coeff[i][1]+coeff[i][4];
y1=x*coeff[i][2]+y*coeff[i][3]+coeff[i][5];
x=x1;
y=y1;
putpixel ((int)(y*64),240-(int)(x*48),c);
}

int main (void)
{
int leaf,color;

// initializing graphics mode
init_gr();

// drawing leafs
for (long int i=0;i<N;i++)
{
    color=colors[random(5)];  // random color
    leaf=random(1000);  // selecting leaf to draw
    if (leaf<=10) iterate(0,color);
    else if (leaf<=860) iterate(1,color);
    else if (leaf<=930) iterate(2,color);
    else iterate(3,color);
}

/* clean up */
getch();
end_gr();

return 0;
}

Рисование линии (по Брезенхэму)

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Bresengham line drawing
//  (c) Johna Smith, 1996
//
//  Method description:
//    Determine  dy=(y2-y1), dx=(x2-x1)
//    We plot first point of the line and increase errors of plotting
//    (xerr and yerr) by dx and dy respectively. If the error is greater
//    than largest from dx and dy then we should correct next point and
//    shift it to 1 pixel by that axe where error was too big.
//
//  ! This program is NOT optimized for best performance !
//  To do so don't use putpixel function - change bytes in video memory directly.
//
//////////////////////////////////////////////////////////////////////////////

#include <graphics.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>

// this function initializes graphics mode
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void init_gr(void)
{
  /* request autodetection */
  int gdriver = DETECT, gmode, errorcode;

  /* initialize graphics mode */
  initgraph(&gdriver, &gmode, "");

  /* read result of initialization */
  errorcode = graphresult();

  if (errorcode != grOk)    /* an error occurred */
  {
      printf("Graphics error: %s\n", grapherrormsg(errorcode));
      printf("Press any key to halt:");
      getch();
      exit(1);               /* return with error code */
  }
}

// this function shuts graphics mode down
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void end_gr(void)
{
closegraph();
}

// this function puts pixel on the screen in (x,y) position using color 'color'
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void PutPixel(int x, int y, int color)
{
putpixel(x,y,color);
}

void BLine(int x1,int y1,int x2, int y2, int color)
{
int delta_x,delta_y,incx,incy,t,distance;
int xerr=0,yerr=0;

// determine dx and dy
delta_x=x2-x1;
delta_y=y2-y1;

// determine steps by x and y axes (it will be +1 if we move in forward
// direction and -1 if we move in backward direction
if (delta_x>0) incx=1;
else if (delta_x==0) incx=0;
else incx=-1;

if (delta_y>0) incy=1;
else if (delta_y==0) incy=0;
else incy=-1;

delta_x=abs(delta_x);
delta_y=abs(delta_y);

// select largest from deltas and use it as a main distance
if (delta_x>delta_y) distance=delta_x;
else distance=delta_y;

for (t=0;t<=distance+1;t++)
{
   PutPixel(x1,y1,color);
   // increasing error
   xerr+=delta_x;
   yerr+=delta_y;
   // if error is too big then we should decrease it by changing
   // coordinates of the next plotting point to make it closer
   // to the true line
   if(xerr>distance)
   {
     xerr-=distance;
     x1+=incx;
   }
   if (yerr>distance)
   {
     yerr-=distance;
     y1+=incy;
   }
}
}

int main(void)
{
// initializing graphics mode
init_gr();

/* examples */
BLine(0, 0, 640, 480,15);
BLine(320, 0, 320, 480,10);
BLine(0, 240, 640, 240,11);
BLine(0, 480, 640, 0,12);
BLine(320, 11, 10, 5,13);
BLine(320, 11, 630, 5,13);

/* clean up */
getch();
end_gr();

return 0;
}

Рисование окружности (по Брезенхэму)

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Bresengham circle drawing
//  (c) Johna Smith, 1996
//
//  Method description:
//    In this algorithm all floating point math changed to sequences
//  of additions and substractions. The main idea of this algorithm
//  is to increase x and y by the value of the error between them.
//
//  ! This program is NOT optimized for best performance !
//  To do so don't use putpixel function - change bytes in video memory directly.
//
//////////////////////////////////////////////////////////////////////////////

#include <graphics.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>

// this function initializes graphics mode
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void init_gr(void)
{
  /* request autodetection */
  int gdriver = DETECT, gmode, errorcode;

  /* initialize graphics mode */
  initgraph(&gdriver, &gmode, "");

  /* read result of initialization */
  errorcode = graphresult();

  if (errorcode != grOk)    /* an error occurred */
  {
     printf("Graphics error: %s\n", grapherrormsg(errorcode));
     printf("Press any key to halt:");
     getch();
     exit(1);               /* return with error code */
  }
}

// this function shuts graphics mode down
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void end_gr(void)
{
closegraph();
}

// this function puts pixel on the screen in (x,y) position using color 'color'
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void PutPixel(int x, int y, int color)
{
putpixel(x,y,color);
}

float const ratio=1.0; // you can change this to draw ellipses

// This function plots points that belongs to the circle
// It recieves offsets from center for the fist quadrant
// and plots symmetrical points in all four quadrants
void plot_circle(int x,int y, int x_center, int y_center, int color)
{
int x_start,y_start,x_end,y_end,x1,y1;

// the distanse between start and end can be greater than 1 if ratio!=1
y_start=y*ratio;
y_end=(y+1)*ratio;
x_start=x*ratio;
x_end=(x+1)*ratio;

for (x1=x_start;x1<x_end;++x1)
{
   // plot points in all four quadrants
   PutPixel(x1+x_center,y+y_center,color);
   PutPixel(x1+x_center,y_center-y,color);
   PutPixel(x_center-x1,y+y_center,color);
   PutPixel(x_center-x1,y_center-y,color);
}

for (y1=y_start;y1<y_end;++y1)
{
   // plot points in all four quadrants
   PutPixel(y1+x_center,x+y_center,color);
   PutPixel(y1+x_center,y_center-x,color);
   PutPixel(x_center-y1,x+y_center,color);
   PutPixel(x_center-y1,y_center-x,color);
}
}

// This is main function that draws circle using function

void Circle(int x1,int y1,int radius, int color)
{
int x,y,delta;

//     Y    *              we start from * and increase x step by step
//          |              decreasing y when needed
//          |
//          |
// --------------------
//          |         X
//          |
//          |
//          |
y=radius;
delta=3-2*radius; // delta is an error
// calculate values for first quadrant
for (x=0;x<y;x++) // x is a main axe
{
   // plot points symmetrically in all quadrants
   plot_circle(x,y,x1,y1,color);
   if (delta<0) delta+=4*x+6;
   else
   {
     delta+=4*(x-y)+10;
     y--; // it's time to decrease y
   }
}
x=y;
if (y!=0) plot_circle(x,y,x1,y1,color);
}

int main(void)
{
// initializing graphics mode
init_gr();

/* examples */
Circle(200,200,100,14);
Circle(300,200,100,15);
Circle(400,200,100,13);
Circle(250,100,100,12);
Circle(350,100,100,11);
Circle(50,400,25,2);
Circle(500,400,25,2);

/* clean up */
getch();
end_gr();
return 0;
}

Сглаживание кривой В-сплайном

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Curve fitting using B-splines
//  (c) Johna Smith, 1996
//
//  Method description:
//    We are drawing lines between points using the following formulas:
//    x(t)=((a3*t+a2)*t+a1)*t+a0
//    y(t)=((b3*t+b2)*t+b1)*t+b0
//    t=0..1
//    Look program for formulas for coefficients ai,bi
//    These coefficients depends on coordinates of current point,
//    previous one, next and next-next ones.
//
//////////////////////////////////////////////////////////////////////////////

#include <graphics.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>

// this function initializes graphics mode
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void init_gr(void)
{
  /* request autodetection */
  int gdriver = DETECT, gmode, errorcode;

  /* initialize graphics mode */
  initgraph(&gdriver, &gmode, "");

  /* read result of initialization */
  errorcode = graphresult();

  if (errorcode != grOk)    /* an error occurred */
  {
     printf("Graphics error: %s\n", grapherrormsg(errorcode));
     printf("Press any key to halt:");
     getch();
     exit(1);               /* return with error code */
  }
}

// this function shuts graphics mode down
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void end_gr(void)
{
closegraph();
}

// this function moves CP to (x,y) position
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void MoveTo(int x, int y)
{
moveto(x,y);
}

// this function draws a line to (x,y) position
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void LineTo(int x, int y)
{
lineto(x,y);
}

// this function draws a line from (x1,y1) to (x2,y2)
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void Line(int x1, int y1, int x2, int y2)
{
line(x1,y1,x2,y2);
}

const N=21; // number of points
const M=300.0; // number of steps between two points

// coordinates of all given points
int x[N]={50,25,50,125,200,225,275,475,590,615,615,615,600,475,275,225,200,125,50,25
,50};
int y[N]={140,100,60,40,70,90,85,75,80,85,100,115,120,125,115,110,130,160,140,100,60
};

// coefficients
float t,xA,xB,xC,xD,yA,yB,yC,yD,a0,a1,a2,a3,b0,b1,b2,b3;
int n,i,j,first;

int main(void)
{
// initializing graphics mode
init_gr();

  /* mark the given points */
  for (i=0;i<N;i++)
  {
    Line(x[i]-4,y[i]-4,x[i]+4,y[i]+4);
    Line(x[i]+4,y[i]-4,x[i]-4,y[i]+4);
  }

  /* main loop */
  first=1;
  for(i=1;i<N-2;i++)
  {
    // calculating coefficients
    xA=x[i-1]; xB=x[i]; xC=x[i+1]; xD=x[i+2];
    yA=y[i-1]; yB=y[i]; yC=y[i+1]; yD=y[i+2];
    a3=(-xA+3*(xB-xC)+xD)/6.0;  b3=(-yA+3*(yB-yC)+yD)/6.0;
    a2=(xA-2*xB+xC)/2.0;        b2=(yA-2*yB+yC)/2.0;
    a1=(xC-xA)/2.0;             b1=(yC-yA)/2.0;
    a0=(xA+4*xB+xC)/6.0;        b0=(yA+4*yB+yC)/6.0;
    for (j=0;j<M;j++) // drawing curve between two given points
    {
      t=(float)j/(float)M;
      if (first)
      {
        first=0;
        MoveTo(((a3*t+a2)*t+a1)*t+a0,((b3*t+b2)*t+b1)*t+b0);
      } else LineTo(((a3*t+a2)*t+a1)*t+a0,((b3*t+b2)*t+b1)*t+b0);
    }
  }

  /* clean up */
  getch();
  end_gr();

  return 0;
}

Рисование куба

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Generating and drawing cube
//  (c) Johna Smith, 1996
//
//  Method description:
//    Cube is one of simplest figures in stereometry.
//    We just need to generate 'n' squares parallel to X,Y and Z axes
//
//////////////////////////////////////////////////////////////////////////////

#include <graphics.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <dos.h>

#define Pi 3.1415926536

enum Action{move,draw};

struct Point3D
{
int x;
int y;
int z;
Action action;
};

// this function initializes graphics mode
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void init_gr(void)
{
/* request autodetection */
int gdriver = DETECT, gmode, errorcode;

/* initialize graphics mode */
initgraph(&gdriver, &gmode, "");

/* read result of initialization */
errorcode = graphresult();

if (errorcode != grOk)    /* an error occurred */
{
    printf("Graphics error: %s\n", grapherrormsg(errorcode));
    printf("Press any key to halt:");
    getch();
    exit(1);               /* return with error code */
}
}

// this function shuts graphics mode down
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void end_gr(void)
{
closegraph();
}

// this function moves CP to (x,y) position
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void MoveTo(int x, int y)
{
moveto(x,y);
}

// this function draws a line to (x,y) position
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void LineTo(int x, int y)
{
lineto(x,y);
}

void draw3Dobject(Point3D *object, int N, float rho, float theta,
                 float phi, float dist_to_screen, int xshift, int yshift)
{
int x,y;
float xe,ye,ze,costh,sinph,cosph,sinth,v11,v12,v13,v21,v22,v32,v33,v23,v43;

// calculating coefficients
costh=cos(theta);
sinth=sin(theta);
cosph=cos(phi);
sinph=sin(phi);
v11=-sinth; v12=-cosph*costh; v13=-sinph*costh;
v21=costh;  v22=-cosph*sinth; v23=-sinph*sinth;
             v32=sinph;        v33=-cosph;
                               v43=rho;
for (int i=0;i<N;i++)
{
    // calculating eye coordinates
    xe=v11*(object+i)->x+v21*(object+i)->y;
    ye=v12*(object+i)->x+v22*(object+i)->y+v32*(object+i)->z;
    ze=v13*(object+i)->x+v23*(object+i)->y+v33*(object+i)->z+v43;

    // calculating screen coordinates
    x=dist_to_screen*xe/ze+xshift;
    y=dist_to_screen*ye/ze+yshift;

    // drawing
    if((object+i)->action==move)
      MoveTo(x,y);
    else
      LineTo(x,y);
}
}

int main(void)
{
const int n=4; // number of cubes segments +1
Point3D cube[15*n]; // coordinates for cubes points
float rho=1900,theta=Pi/3,phi=2*Pi/3,dist_to_screen=600; // view point
int xshift=300, yshift=140; // picture offset
float side=300; // cubes parameters
float delta;// auxulary variable

// initializing graphics mode
init_gr();

// generating cube
delta=side/(n-1);
for (int i=0;i<n;i++)
{
   for(int j=i*5;j<i*5+5;j++)
   {
     cube[j].x=i*delta;
     cube[j].action=draw;
   }
   cube[i*5].y=0;
   cube[i*5].z=0;
   cube[i*5].action=move;
   cube[i*5+1].y=side;
   cube[i*5+1].z=0;
   cube[i*5+2].y=side;
   cube[i*5+2].z=side;
   cube[i*5+3].y=0;
   cube[i*5+3].z=side;
   cube[i*5+4].y=0;
   cube[i*5+4].z=0;
}
int c=5*n;
for (i=0;i<n;i++)
{
   for(int j=i*5;j<i*5+5;j++)
   {
     cube[c+j].y=i*delta;
     cube[c+j].action=draw;
   }
   cube[c+i*5].x=0;
   cube[c+i*5].z=0;
   cube[c+i*5].action=move;
   cube[c+i*5+1].x=side;
   cube[c+i*5+1].z=0;
   cube[c+i*5+2].x=side;
   cube[c+i*5+2].z=side;
   cube[c+i*5+3].x=0;
   cube[c+i*5+3].z=side;
   cube[c+i*5+4].x=0;
   cube[c+i*5+4].z=0;
}
c=10*n;
for (i=0;i<n;i++)
{
   for(int j=i*5;j<i*5+5;j++)
   {
      cube[c+j].z=i*delta;
      cube[c+j].action=draw;
   }
   cube[c+i*5].y=0;
   cube[c+i*5].x=0;
   cube[c+i*5].action=move;
   cube[c+i*5+1].y=side;
   cube[c+i*5+1].x=0;
   cube[c+i*5+2].y=side;
   cube[c+i*5+2].x=side;
   cube[c+i*5+3].y=0;
   cube[c+i*5+3].x=side;
   cube[c+i*5+4].y=0;
   cube[c+i*5+4].x=0;
}

// drawing
draw3Dobject(cube,15*n,rho,theta,phi,dist_to_screen,xshift,yshift);

/* clean up */
getch();
end_gr();

return 0;
}

Рисование тора

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Generating and drawing torus
//  (c) Johna Smith, 1996
//
//  Method description:
//    Torus has two main circles, which are described by
//   two systems of eqations:
//    x=Rcos(a)     x=R+rcos(b)
//    y=Rsin(a)     y=0
//    z=0           z=rsin(b)
//    By generating this circles (approximating by polygons) we can get torus
//
//////////////////////////////////////////////////////////////////////////////

#include <graphics.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <dos.h>

#define Pi 3.1415926536

enum Action{move,draw};

struct Point3D
{
int x;
int y;
int z;
Action action;
};

// this function initializes graphics mode
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void init_gr(void)
{
  /* request autodetection */
  int gdriver = DETECT, gmode, errorcode;

  /* initialize graphics mode */
  initgraph(&gdriver, &gmode, "");

  /* read result of initialization */
  errorcode = graphresult();

  if (errorcode != grOk)    /* an error occurred */
  {
     printf("Graphics error: %s\n", grapherrormsg(errorcode));
     printf("Press any key to halt:");
     getch();
     exit(1);               /* return with error code */
  }
}

// this function shuts graphics mode down
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void end_gr(void)
{
closegraph();
}

// this function moves CP to (x,y) position
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void MoveTo(int x, int y)
{
moveto(x,y);
}

// this function draws a line to (x,y) position
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void LineTo(int x, int y)
{
lineto(x,y);
}

void draw3Dobject(Point3D *object, int N, float rho, float theta,
                 float phi, float dist_to_screen, int xshift, int yshift)
{
int x,y;
float xe,ye,ze,costh,sinph,cosph,sinth,v11,v12,v13,v21,v22,v32,v33,v23,v43;

// calculating coefficients
costh=cos(theta);
sinth=sin(theta);
cosph=cos(phi);
sinph=sin(phi);
v11=-sinth; v12=-cosph*costh; v13=-sinph*costh;
v21=costh;  v22=-cosph*sinth; v23=-sinph*sinth;
             v32=sinph;        v33=-cosph;
                               v43=rho;
for (int i=0;i<N;i++)
{
     // calculating eye coordinates
     xe=v11*(object+i)->x+v21*(object+i)->y;
     ye=v12*(object+i)->x+v22*(object+i)->y+v32*(object+i)->z;
     ze=v13*(object+i)->x+v23*(object+i)->y+v33*(object+i)->z+v43;

     // calculating screen coordinates
     x=dist_to_screen*xe/ze+xshift;
     y=dist_to_screen*ye/ze+yshift;

     // drawing
     if((object+i)->action==move)
       MoveTo(x,y);
     else
       LineTo(x,y);
}
}

int main(void)
{
const int n=20;  // number of torus' segments
Point3D torus[2*n*(n+1)]; // coordinates for torus' points
float rho=1800,theta=0,phi=3*Pi/4,dist_to_screen=600; // view point
int xshift=300, yshift=250; // picture offset
float delta=2.0*Pi/n, r=75, R=300; // torus' parameters
float alpha,cosa,sina,beta,x; // auxulary variables

// initializing graphics mode
init_gr();

// generating torus
for (int i=0;i<n;i++)
{
    alpha=i*delta;
    cosa=cos(alpha);
    sina=sin(alpha);
    for (int j=0;j<n+1;j++)
    {
      beta=j*delta;
      x=R+r*cos(beta);
      torus[i*(n+1)+j].x=cosa*x;
      torus[i*(n+1)+j].y=sina*x;
      torus[i*(n+1)+j].z=r*sin(beta);
      torus[i*(n+1)+j].action=((i==0 && j==0)?move:draw);
    }
}
int c=n*n+n;
for (i=0;i<n;i++)
{
    beta=i*delta;
    x=R+r*cos(beta);
    for (int j=0;j<n+1;j++)
    {
      alpha=j*delta;
      cosa=cos(alpha);
      sina=sin(alpha);
      torus[c+i*(n+1)+j].x=cosa*x;
      torus[c+i*(n+1)+j].y=sina*x;
      torus[c+i*(n+1)+j].z=r*sin(beta);
      torus[c+i*(n+1)+j].action=draw;
    }
}

// drawing
draw3Dobject(torus,2*n*(n+1),rho,theta,phi,dist_to_screen,xshift,yshift);

/* clean up */
getch();
end_gr();

return 0;
}

Рисование 3D объекта

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Drawing 3D object
//  (c) Johna Smith, 1996
//
//  Method description:
//    Function draw3Dobject recieves the following parameters:
//     object - reference to array of points of the drawn object
//     N - number of points in this array
//     rho     \
//     phi     |- spherical coordinates of point of view
//     theta   /
//     dist_to_screen - distance from viewpoint to screen
//     xshift, yshift - picture will be shifted by this values
//
//////////////////////////////////////////////////////////////////////////////

#include <graphics.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <dos.h>

#define Pi 3.1415926536

enum Action{move,draw};

struct Point3D
{
int x;
int y;
int z;
Action action;
};

// this function initializes graphics mode
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void init_gr(void)
{
       /* request autodetection */
       int gdriver = DETECT, gmode, errorcode;

       /* initialize graphics mode */
       initgraph(&gdriver, &gmode, "");

       /* read result of initialization */
       errorcode = graphresult();

       if (errorcode != grOk)    /* an error occurred */
     {
               printf("Graphics error: %s\n", grapherrormsg(errorcode));
               printf("Press any key to halt:");
                                        getch();
               exit(1);               /* return with error code */
       }
}

// this function shuts graphics mode down
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void end_gr(void)
{
closegraph();
}

// this function moves CP to (x,y) position
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void MoveTo(int x, int y)
{
moveto(x,y);
}

// this function draws a line to (x,y) position
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void LineTo(int x, int y)
{
lineto(x,y);
}

void draw3Dobject(Point3D *object, int N, float rho, float theta,
                 float phi, float dist_to_screen, int xshift, int yshift)
{
int x,y;
float xe,ye,ze,costh,sinph,cosph,sinth,v11,v12,v13,v21,v22,v32,v33,v23,v43;

// calculating coefficients
costh=cos(theta);
sinth=sin(theta);
cosph=cos(phi);
sinph=sin(phi);
v11=-sinth; v12=-cosph*costh; v13=-sinph*costh;
v21=costh;  v22=-cosph*sinth; v23=-sinph*sinth;
                                 v32=sinph;        v33=-cosph;
                               v43=rho;
for (int i=0;i<N;i++)
{
     // calculating eye coordinates
     xe=v11*(object+i)->x+v21*(object+i)->y;
     ye=v12*(object+i)->x+v22*(object+i)->y+v32*(object+i)->z;
     ze=v13*(object+i)->x+v23*(object+i)->y+v33*(object+i)->z+v43;
     // calculating screen coordinates
     x=dist_to_screen*xe/ze+xshift;
     y=dist_to_screen*ye/ze+yshift;
     // drawing
     if((object+i)->action==move)
         MoveTo(x,y);
     else
         LineTo(x,y);
}
}

int main(void)
{
Point3D thetr[]={{100,100,100,move},
                             {100,200,100,draw},
                             {200,100,100,draw},
                             {100,100,100,draw},
                             {100,100,200,draw},
                             {200,100,100,draw},
                             {100,100,200,move},
                             {100,200,100,draw}};
int N=sizeof(thetr)/sizeof(Point3D);
float rho=700,theta=Pi/4,phi=Pi/4,dist_to_screen=300;
int xshift=300, yshift=150;

// initializing graphics mode
init_gr();

/* examples */
while (!kbhit())
{
     theta+=0.05;
     // rotating viewpoint
     if (phi>2*Pi) phi-=2*Pi;
     setcolor(11);
     draw3Dobject(thetr,N,rho,theta,phi,dist_to_screen,xshift,yshift);
     setcolor(0);
     delay(15);
     draw3Dobject(thetr,N,rho,theta,phi,dist_to_screen,xshift,yshift);
}

/* clean up */
getch();
end_gr();

return 0;
}

Вращение 3D объекта

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Rotating 3D object
//  (c) Johna Smith, 1996
//
//  Method description:
//   Given: vector X, angle a
//   To rotate this vector about axe X we should apply the following
//   operation:  X'=X*Rx, where Rx - is matrix of rotation
//        [ 1 0      0     0 ]
//        [ 0 cos a  sin a 0 ]
//     Rx=[ 0 -sin a cos a 0 ]
//        [ 0 0      0     1 ]
//   Applying this operation to every point of the given object we
//   rotate it.
//
//////////////////////////////////////////////////////////////////////////////

#include <graphics.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <dos.h>

#define Pi 3.1415926536

enum Action{move,draw};

struct Point3D
{
float x;
float y;
float z;
Action action;
};

// this function initializes graphics mode
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void init_gr(void)
{
       /* request autodetection */
       int gdriver = DETECT, gmode, errorcode;

       /* initialize graphics mode */
       initgraph(&gdriver, &gmode, "");

       /* read result of initialization */
       errorcode = graphresult();

       if (errorcode != grOk)    /* an error occurred */
       {
              printf("Graphics error: %s\n", grapherrormsg(errorcode));
              printf("Press any key to halt:");
              getch();
              exit(1);               /* return with error code */
       }
}

// this function shuts graphics mode down
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void end_gr(void)
{
closegraph();
}

// this function moves CP to (x,y) position
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void MoveTo(int x, int y)
{
moveto(x,y);
}

// this function draws a line to (x,y) position
// it will work only if you're using Borland C++ compiler & BGI drivers
// if you're using another compiler you should overwrite body of this function

void LineTo(int x, int y)
{
lineto(x,y);
}

void draw3Dobject(Point3D *object, int N, float rho, float theta,
                 float phi, float dist_to_screen, int xshift, int yshift)
{
int x,y;
float xe,ye,ze,costh,sinph,cosph,sinth,v11,v12,v13,v21,v22,v32,v33,v23,v43;

// calculating coefficients
costh=cos(theta);
sinth=sin(theta);
cosph=cos(phi);
sinph=sin(phi);
v11=-sinth; v12=-cosph*costh; v13=-sinph*costh;
v21=costh;  v22=-cosph*sinth; v23=-sinph*sinth;
                     v32=sinph;             v33=-cosph;
                                                    v43=rho;
for (int i=0;i<N;i++)
{
        // calculating eye coordinates
        xe=v11*(object+i)->x+v21*(object+i)->y;
        ye=v12*(object+i)->x+v22*(object+i)->y+v32*(object+i)->z;
        ze=v13*(object+i)->x+v23*(object+i)->y+v33*(object+i)->z+v43;

        // calculating screen coordinates
        x=dist_to_screen*xe/ze+xshift;
        y=dist_to_screen*ye/ze+yshift;

       // drawing
        if((object+i)->action==move)
             MoveTo(x,y);
        else
             LineTo(x,y);
}
}

void RotateX(Point3D *object, int n, float angle)
{
float cosa,sina,y,z;

cosa=cos(angle);
sina=sin(angle);
for(int i=0;i<n;i++)
{
        y=(object+i)->y*cosa-(object+i)->z*sina;
        z=(object+i)->y*sina+(object+i)->z*cosa;
        (object+i)->y=y;
        (object+i)->z=z;
}
}

void RotateY(Point3D *object, int n, float angle)
{
float cosa,sina,x,z;

cosa=cos(angle);
sina=sin(angle);
for(int i=0;i<n;i++)
{
        x=(object+i)->x*cosa+(object+i)->z*sina;
        z=-(object+i)->x*sina+(object+i)->z*cosa;
        (object+i)->x=x;
        (object+i)->z=z;
}
}

void RotateZ(Point3D *object, int n, float angle)
{
float cosa,sina,x,y;

cosa=cos(angle);
sina=sin(angle);
for(int i=0;i<n;i++)
{
        x=(object+i)->x*cosa-(object+i)->y*sina;
        y=(object+i)->x*sina+(object+i)->y*cosa;
        (object+i)->x=x;
        (object+i)->y=y;
}
}

int main(void)
{
Point3D thetr[]={{100,100,100,move},
                             {100,200,100,draw},
                             {200,100,100,draw},
                             {100,100,100,draw},
                             {100,100,200,draw},
                             {200,100,100,draw},
                             {100,100,200,move},
                             {100,200,100,draw}};

Point3D cube[]={{-50,-50,-50,move},
                             {-50,50,-50,draw},
                             {-50,50,50,draw},
                             {50,50,50,draw},
                             {50,50,-50,draw},
                             {50,-50,-50,draw},
                             {50,-50,50,draw},
                             {-50,-50,50,draw},
                             {-50,50,50,draw},
                             {-50,-50,50,move},
                             {-50,-50,-50,draw},
                             {50,-50,-50,draw},
                             {50,50,-50,move},
                             {-50,50,-50,draw},
                             {50,-50,50,move},
                             {50,50,50,draw}};

int N=sizeof(thetr)/sizeof(Point3D);
int M=sizeof(cube)/sizeof(Point3D);
float rho=1500,theta=Pi/4,phi=-3*Pi/4,dist_to_screen=700;
int xshift=300, yshift=150,xshift1=200,yshift1=200;

// initializing graphics mode
init_gr();

/* examples */
while (!kbhit())
{
         RotateX(cube,M,Pi/24);
         RotateY(cube,M,Pi/24);
         RotateZ(cube,M,Pi/24);
         RotateY(thetr,N,Pi/30);
         setcolor(12);
         draw3Dobject(thetr,N,rho,theta,phi,dist_to_screen,xshift,yshift);
         draw3Dobject(cube,M,rho,theta,phi,dist_to_screen,xshift1,yshift1);
         setcolor(0);
         delay(35);
         draw3Dobject(thetr,N,rho,theta,phi,dist_to_screen,xshift,yshift);
         draw3Dobject(cube,M,rho,theta,phi,dist_to_screen,xshift1,yshift1);
}

/* clean up */
getch();
end_gr();

return 0;
}

Ахинея

10.01.2010 - 19:21

Матрицы:

Определитель

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating det A, where A is matrix N x N
//  (c) Johna Smith, 1996
//
//  Method description:
//    This method based on two statements:
//    1) det A = a11, if A is matrix 1 x 1
//    2) |a11 a12 a13 ... a1n|                       |a11 a12 .. a1(i-1) a1(i+1) .. a1n |
//       |a21 a22 a23 ... a2n|    n     i+1          |a21 a22 .. a2(i-1) a2(i+1) .. a2n |
// det A=|a31 a32 a33 ... a3n| = SUM (-1)  a1i * det |a31 a32 .. a3(i-1) a3(i+1) .. a3n |
//       |...................|   i=1                 |................................  |
//       |...................|                       |an1 an2 .. an(i-1) an(i+1) .. ann |
//       |an1 an2 an3 ... ann|
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <alloc.h>

#define N       5  // matrix size

int matrix[N][N]={{1,2,3,4,5},{2,1,3,4,5},{3,2,1,4,5},{4,3,2,1,5},{5,4,3,2,1}};

void show_matrix(int *matrix, int size) // this functions displays matrix
{
for(int i=0;i<size;i++)
{
  for(int j=0;j<size;j++)
    printf("%d ",*(matrix+i*size+j));
  printf("\n");
}
printf("\n");
}

int Det(int *matrix, int const size)
{
int *submatrix;
int det=0;

if (size>1)
{
   submatrix=(int *)malloc(sizeof(int)*(size-1)*(size-1));
   for (int i=0;i<size;i++)
   {
     // creating new array (submatrix)
     for (int j=0;j<size-1;j++)
       for (int k=0;k<size-1;k++)
         *(submatrix+j*(size-1)+k)=*(matrix+(j+1)*size+(k<i?k:k+1));
     // calling recursively function Det using submatrix as a parameter
     // and adding the result to final value
     det+=*(matrix+i)*(i/2.0==i/2?1:-1)*Det(submatrix,size-1);
   }
   free(submatrix);
} else det=*matrix;
return det;
}

void main(void)
{
// show given matrix
show_matrix(matrix[0],N);

// calculating determinante and displaying the result
printf("det A = %d",Det(matrix[0],N));

Обращение

» Кликните сюда для просмотра оффтоп текста.. «

Умножение

» Кликните сюда для просмотра оффтоп текста.. «

Транспонирование

» Кликните сюда для просмотра оффтоп текста.. «

Ахинея

12.01.2010 - 13:33

Специальные функции:

arccos x (x - real)

» Кликните сюда для просмотра оффтоп текста.. «

arccos z (z - complex)

» Кликните сюда для просмотра оффтоп текста.. «

arcsin x

» Кликните сюда для просмотра оффтоп текста.. «

arcsin z

» Кликните сюда для просмотра оффтоп текста.. «

arctg x

» Кликните сюда для просмотра оффтоп текста.. «

arctg z

» Кликните сюда для просмотра оффтоп текста.. «

Функция Бесселя

» Кликните сюда для просмотра оффтоп текста.. «

Гамма-функция

» Кликните сюда для просмотра оффтоп текста.. «

Бета-функция

» Кликните сюда для просмотра оффтоп текста.. «

Ахинея

14.01.2010 - 18:16

Методы поиска

Binary search in sorted string

» Кликните сюда для просмотра оффтоп текста.. «

Substring search. Boyer-Moore algorithm

» Кликните сюда для просмотра оффтоп текста.. «

Linear substring search

» Кликните сюда для просмотра оффтоп текста.. «

Substring search. Knuth-Morris-Pratt algorithm

» Кликните сюда для просмотра оффтоп текста.. «

Linear search in the string

» Кликните сюда для просмотра оффтоп текста.. «

Ахинея

22.01.2010 - 17:05

Решение систем линейных уравнений

Метод квадратного корня

» Кликните сюда для просмотра оффтоп текста.. «

Метод Некрасова

» Кликните сюда для просмотра оффтоп текста.. «

Метод 3-диагональных матриц

» Кликните сюда для просмотра оффтоп текста.. «

Метод Гаусса

» Кликните сюда для просмотра оффтоп текста.. «

Метод Гаусса с выбором максимального элемента

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Gauss method with max element selecting)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This methos based on excluding variables from equations.
//    Using first equation we exclude x1 from all other equations
//    Using second equation equation we exclude x2 from all equations
//    (excluding first) and so on. So we have triangular matrix from which
//    we can easily get vector X. (AX=B)
//    'Selecting max element' means that on each step we select equation
//    with max element in current column and use it as next equation to exclude
//    variable.
//
//////////////////////////////////////////////////////////////////////////////

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

#define N      4     // size of matrix
#define N1     N+1

float matrix[N][N1]=
{{2,5,4,1,28},
{1,3,2,1,17},
  {2,10,9,7,77},
{3,8,9,2,54}
       };

void ShowMatrix(void)
{
for (int i=0;i<N;i++)
{
   for (int j=0;j<N;j++)
     printf("%+f*x%d",matrix[i][j],i+1);
   printf("=%f\n",matrix[i][N]);
}
}

void main(void)
{
// Variables declaration
float max,swp,tmp,x[N1];
register short int row_with_max_element,i,j,k;

// Printing given matrix
ShowMatrix();
// Main loop
for (i=0;i<N;i++)
{
   // Searching for max element in the current column (i)
   max=fabs(matrix[i][i]);
   row_with_max_element=i;
   for (j=i+1;j<N;j++)
     if (fabs(matrix[j][i])>max)
     {
       row_with_max_element=j;
       max=fabs(matrix[j][i]);
     }
   // Swapping 2 lines of matrix - row_with_max_element & i
   if (row_with_max_element!=i)
   for (j=0;j<N1;j++)
   {
     swp=matrix[row_with_max_element][j];
     matrix[row_with_max_element][j]=matrix[i][j];
     matrix[i][j]=swp;
   }
   // Excluding variable x[i] from equations
   tmp=matrix[i][i];
   for (j=N;j>=i;j--) matrix[i][j]/=tmp;
   for (j=i+1;j<N;j++)
   {
     tmp=matrix[j][i];
     for (k=N;k>=i;k--)
       matrix[j][k]-=tmp*matrix[i][k];
   }
}

// Calculating vector x
x[N-1]=matrix[N-1][N];
for (i=N-2;i>=0;i--)
{
   x[i]=matrix[i][N];
   for (j=i+1;j<N;j++) x[i]-=matrix[i][j]*x[j];
}

// Printing solution
printf("\nSolution:\n");
for (i=0;i<N;i++)
   printf("x%d=%f\n",i+1,x[i]);
}

Метод итераций

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (iterations method)
//  (c) Johna Smith, 1996
//
//  Method description:
//  This program automatically transforms system of equations to
//  iterative form and solves it.
//  Iterative form:
//     x1=a11*x1+a12*x2+... + a1n*xn+b1
//     x2=a21*x1+a22*x2+... + a2n*xn+b2
//     ...
//     xn=an1*x1+an2*x2+... + ann*xn+bn
//  Setting first iteration of vector X we can get next iteration from
//  such form of system of equations. Iterations will be finished when
//  difference between two consequitive iterations will be less than epsilon
//
//////////////////////////////////////////////////////////////////////////////

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

#define N      3     // size of matrix
#define N1     N+1

float matrix[N][N1]=
{{14.38,-2.41,1.39,5.86},
{1.84,25.36,-3.31,-2.28},
{2.46,-3.49,16.37,4.47}
};
int maxiterations=10;  // maximum number of iterations
float epsilon=0.0001;  // required accuracy

void ShowMatrix(void)
{
for (int i=0;i<N;i++)
{
   for (int j=0;j<N;j++)
     printf("%+f*x%d",matrix[i][j],i+1);
   printf("=%f\n",matrix[i][N]);
}
}

void main(void)
{
// Variables declaration
float x[N],y[N],t;
register short int i,j,k;
int iterations=0;

// Printing given matrix
ShowMatrix();

// setting first iteration of vector X
for (i=0;i<N;i++) x[i]=matrix[i][N];
do
{
   for (i=0;i<N;i++)
   {
     t=-matrix[i][N];
     for (j=0;j<N;j++)
       t+=matrix[i][j]*x[j];
       y[i]=(-t+matrix[i][i]*x[i])/matrix[i][i];
   }
   iterations++;
   k=0;
   // checking solution
   while (fabs(x[k]-y[k])<epsilon && k<N) k++;
   // new iteration becomes old
   for(i=0;i<N;i++) x[i]=y[i];
} while (k!=N && iterations<maxiterations);

if (iterations==maxiterations)
{
   printf("Iterations are very slow...");
} else
{
   // Printing solution
   printf("\nSolution:\n");
   for (i=0;i<N;i++)
     printf("x%d=%f\n",i+1,x[i]);
   printf("%d iterations were made",iterations);
}
}

Метод Монте-Карло

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Monte-Karlo method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    We should translate given system of linear equations to special form:
//      x1=a(11)x1+a(12)x2+...+a(1n)xn+a(1n+1)
//      x2=a(21)x2+a(22)x2+...+a(2n)xn+a(2n+1)
//      ....
//      xn=a(n1)x1+a(n2)x2+...+a(nn)xn+a(nn+1)
//            n
//     where SUM |aij|<1 (i=1, 2, ... ,n)
//           j=1
//    Lets divide interval [0;1] into N+1 smaller intervals.
//    Imagine a particle that is moving randomly at [0;1] interval.
//    Remember its moves until it reachs one of the bounds of the interval.
//    So we'll get a trajectory of the particle s(i),s(j),s(k),s(l),s(m),...
//    If particle moves from interval s(i) to s(j) we'll write v(ij)
//    y[i]=v(ij)*v(jk)*...*v(tm)*w(m), v(ij)=sign(aij), v(jk)=sign(ajk)
//                    n
//   w(m)=a(mn+1)/(1-SUM a(mj)),   x[i]=y/M, where M is number of particle runs
//                   j=1
//
//////////////////////////////////////////////////////////////////////////////

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

#define N     3     // size of matrix
#define N1    N+1
#define M     1000 // number of tries

float matrix[N][N1]=
{{3,1,-1,-2},
{1,4,-2,-3},
{2,-1,5,-15}
};

void ShowMatrix(void)
{
for (int i=0;i<N;i++)
{
   printf("x%d=",i+1);
   for (int j=0;j<N;j++)
     printf("%+f*x%d",matrix[i][j],j+1);
   printf("%+f\n",matrix[i][N]);
}
}

void main(void)
{
// Variables declaration
float b[N][N],w[N],y,c,x[N],l;
register short int i,j,v;
unsigned char finish;
int row,   // currently analysed row
     tries; // number of tries to calculate true root by generating randoms

// Transforming matrix
for (i=0;i<N;i++)
{
   l=0;
   for (j=0;j<N1;j++) l+=fabs(matrix[i][j]);
   v=(matrix[i][i]>0?-1:1);
   for (j=0;j<N1;j++) matrix[i][j]=v*(matrix[i][j])/l+(i==j?1:0);
}
// Printing transformed matrix
ShowMatrix();

// filling matrix B
for (i=0;i<N;i++)
{
   b[i][0]=fabs(matrix[i][0]);
   for(j=1;j<N;j++)
     b[i][j]=b[i][j-1]+fabs(matrix[i][j]);
}
// filling matrix w
for(i=0;i<N;i++) w[i]=matrix[i][N]/(1-b[i][N-1]);

// for each equation in the system
for(i=0;i<N;i++)
{
   tries=0;
   row=i;
   y=0;
   v=1;
   while (tries<M)
   {
     // generating random number
     c=random(32767)/32767.0;
     j=N-1;
     finish=0;
     while(j>=0 && !finish)
     {
       if (c>b[row][j])
       {
         if (j==N-1)
         {
           tries++;
           y+=v*w[row];
           row=i;v=1;
         }
         else
         {
           v*=(matrix[row][j+1]>=0?1:-1);
           row=j+1;
         }
         finish=1;
       }
       j--;
     }
     if (!finish)
     {
       v*=(matrix[row][j+1]>=0?1:-1);
       row=0;
     }
   }
   // calculating root
   x[i]=y/M;
}

// Printing solution
printf("\nSolution:\n");
for (i=0;i<N;i++)
   printf("x%d=%f\n",i+1,x[i]);
}

Метод ортогонализации

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (orthogonalization method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given system of linear equations:
//      a11 x1 + a12 x2 + ... + a1n xn + a(1n+1) = 0
//      a21 x1 + a22 x2 + ... + a2n xn + a(2n+1) = 0
//
//      an1 x1 + an2 x2 + ... + ann xn + a(nn+1) = 0
//
//    Left part of each equation is a result of scalar multiplication of
//    two vectors: ai=(ai1,ai2,...,ain,a(in+1)) and x=(x1,x2,..,xn,1)
//    So to solve this system we need only to build vector x which will be
//    orthogonal to each of ai vectors.
//                            n+1
//    Let u1=a1,   b1=u1/sqrt(SUM u(1j)^2)
//                            j=1
//
//    Other ui we can get from the following iterative formula:
//                      k
//       u(k+1)=u(k+1)-SUM {u(k+1), bj}bj, where {..} means scalar multiplication
//                     j=1
//                          n+1
//       b(k+1)=u(k+1)/sqrt(SUM u(k+1,j)^2)
//                          j=1
//    And finally we can obtain roots of this system from this formula:
//
//    xi=b(n+1,i)/b(n+1n+1)
//
//////////////////////////////////////////////////////////////////////////////

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

#define N     3     // size of matrix
#define N1    N+1

double matrix[N1][N1]=
{{3,1,-1,-2},
{1,4,-2,-3},
{2,-1,5,-15}
};

void ShowMatrix(void)
{
for (int i=0;i<N;i++)
{
   for (int j=0;j<N;j++)
     printf("%+f*x%d",matrix[i][j],i+1);
   printf("%+f=0\n",matrix[i][N]);
}
}

void main(void)
{
// Variables declaration
double c[N],x[N],s;
int i,j,k,l,m;

ShowMatrix();

// expanding system by vector (0,0,0,...,0,1)
for (i=0;i<N;i++) matrix[N][i]=0;
matrix[N][N]=1;

// for all equations
for (i=0;i<N1;i++)
{
   if (i>0)
   {
     k=0;
     // make some iterations to increase accuracy of calculations
     while (k<=3)
     {
       for (l=0;l<i;l++)
       {
         c[l]=0;
         for(m=0;m<N1;m++) c[l]+=matrix[l][m]*matrix[i][m];
       }
       for (j=0;j<N1;j++)
       {
         s=0;
         for(l=0;l<i;l++) s+=c[l]*matrix[l][j];
         matrix[i][j]-=s;
       }
       k++;
     }
   }
    // normalizing vector
    s=0;
    for (k=0;k<N1;k++) s+=matrix[i][k]*matrix[i][k];
    s=sqrt(s);
    for (j=0;j<N1;j++) matrix[i][j]/=s;
}

s=matrix[N][N];
for (j=0;j<N;j++) x[j]=matrix[N][j]/s;

// Printing solution
printf("\nSolution:\n");
for (i=0;i<N;i++)
   printf("x%d=%f\n",i+1,x[i]);
}

Метод Зэйделя

» Кликните сюда для просмотра оффтоп текста.. «

Ахинея

28.01.2010 - 19:18

Решение нелинейных уравнений

Метод хорд

» Кликните сюда для просмотра оффтоп текста.. «

Метод дихотомии

» Кликните сюда для просмотра оффтоп текста.. «

Метод Эйткена-Стеффенсона

» Кликните сюда для просмотра оффтоп текста.. «

Метод итераций

» Кликните сюда для просмотра оффтоп текста.. «

Метод Ньютона

» Кликните сюда для просмотра оффтоп текста.. «

Метод Лобачевского

» Кликните сюда для просмотра оффтоп текста.. «

Код

//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (Lobachevsky method)
//  (c) Johna Smith, 1996
//
//  Method description:
//   Given: a0+a1x+a2x^2+...+anx^n=0
//   This method allows to find modulus of the greatest root of this equation
//   even if it's complex. But in last case there can appear several messages
//   about impossibilty of calculation root of negative number.
//   The main idea of this method is to change given equation to other
//   equation which roots equals to powered roots of given equation. For example
//   if roots of the given equation are x0,x1,.. xn then roots of new equation
//   will be x0^2, x1^2, ..., xn^2. Repeating this operation we get an equation
//   where one root is much greater than other ones. So we can easily
//   obtain modulus of the greatrest root of the given equation.
//   To obtain other roots of equation we need to divide given equation
//   by (x-x0) (where x0 is found root) and apply this method to result.
//
//////////////////////////////////////////////////////////////////////////////

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

#define N   4
#define N1  N+1
#define Iterations  15  // number of iterations

double a[N1]={24,-50,35,-10,1};

void main(void)
{
double r,b[N1],c[N1],g,bi,d;
int z,k;

// printing given equation
printf("%f",a[0]);
for(int i=1;i<N1;i++) printf("%+fx^%d",a[i],i);
printf("=0\n\n");

// preparing auxiliary arrays b and c
for (i=0;i<N1;i++)
{
   b[i]=a[i]/a[N];
   c[i]=0;
}

// setting required parameters
r=1/2.0;
g=1;

// make all iterations
for(int y=0;y<Iterations;y++)
{
   // calculate coefficients c[i] (coefficients of new equation)
   z=1;
   for(i=0;i<N1;i++)
   {
     bi=z*b[i];
     k=(i+1)/2;
     for(int j=i%2;j<N1;j+=2)
     {
       c[k]+=bi*b[j];
       k++;
     }
     z=-z;
   }
   d=z*c[N-1];
   // check whether we could calculate root of d
   if(d>0)
   {
     // calculating and printing new iteration
     g*=powl(d,r);
     printf("%f\n",g);
     for (i=0;i<N1;i++)
     {
       // preparing data for next iteration
       b[i]=c[i]/powl(d,N-i);
       c[i]=0;
     }
     b[N-1]=z;
     b[N]=-z;
   } else
   {
     // d is negative - can't calculate root
     for(i=0;i<N1;i++)
     {
       // preparing data for next iteration
       b[i]=c[i];
       c[i]=0;
     }
     printf("no iteration (can't calculate root from negative number)\n");
   }
   r/=2.0;
}
}

Ахинея

5.03.2010 - 15:26

Решение дифференциальных уравнений

Метод Эйлера

» Кликните сюда для просмотра оффтоп текста.. «

Метод Адамса

» Кликните сюда для просмотра оффтоп текста.. «

Метод Эйлера-Коши

» Кликните сюда для просмотра оффтоп текста.. «

Метод Рунге-Кутта

» Кликните сюда для просмотра оффтоп текста.. «

Сладкая Смерть

3.12.2012 - 8:56

Алгоритм Дейкстры
Неплохой такой алгоритм, я его в дипломе использовала для 3D навигации по универу. Ищет кратчайший путь до заданного объекта от начальной точки.
Неформальное объяснение

» Кликните сюда для просмотра оффтоп текста.. «

Пример кода, который нашла

» Кликните сюда для просмотра оффтоп текста.. «

Мой пример (писала в Delphi)

» Кликните сюда для просмотра оффтоп текста.. «

Цитата

//* uGraph - модуль для работы с графами - алгоритмы на графах

unit uGraph;

interface

uses
uPath;

type
TVertex = record
id : Integer;
x,y,z : Single;
Tag : Integer;
Marked : Boolean;
DistFromStart : Single;
PrevVertex : Integer;
end;

TGraph = class
public
//* Public declarations

M : array of array of Single;
V : array of TVertex;
N : Integer;
Path : TPath;
//*
StartVertex : Integer;
FinishVertex : Integer;
//*
constructor Create;
procedure LoadGraph (const FileName : ShortString);
procedure FindPath (const Tag : Integer );
end;

implementation

//* TGraph *//

procedure TGraph.LoadGraph(const FileName: ShortString);
var
F : TextFile;
i,j : Integer;

begin
AssignFile (F, FileName);
Reset (F);
ReadLn (F, N);

SetLength (M, N, N);
SetLength (V, N);

//

for i := 0 to N-1 do
begin
ReadLn (F, V[i].x, V[i].y, V[i].z, v[i].tag);
StartVertex := 0;
end;

for i := 0 to N-1 do
for j := 0 to N-1 do
begin
try
Read (F, M[i,j]);
except
M[i,j] := 10000000;
end;
end;

//

CloseFile (F);
end;

procedure TGraph.FindPath(const Tag: Integer);
var
i,j : Integer;
NotMarked : Integer;
vm,pv : Integer;
MinDist : Single;

begin
for i := 0 to N-1 do
begin
v[i].Marked := False;
v[i].DistFromStart := M[StartVertex,i];

if M[StartVertex,i] = -1 then
v[i].DistFromStart := 10000000;

v[i].id := i;

if v[i].Tag = Tag then FinishVertex := i;
end;

//*

v[StartVertex].Marked := True;
v[StartVertex].PrevVertex := -1;
NotMarked := N-1;

while NotMarked <> 0 do
begin
MinDist := 10000000;

for i := 0 to N-1 do
if (not v[i].Marked) and (v[i].DistFromStart < MinDist) and (v[i].DistFromStart <> -1) then
begin
vm := i;
MinDist := v[i].DistFromStart;
end;

v[vm].Marked := True;
NotMarked := NotMarked-1;

for i := 0 to N-1 do
if not v[i].Marked then
if (v[i].DistFromStart > v[vm].DistFromStart+M[vm,i]) and
(v[i].DistFromStart <> -1) and (v[vm].DistFromStart <> -1) and
(M[vm,i] <> -1) then
begin
v[i].DistFromStart := v[vm].DistFromStart+M[vm,i];
v[i].PrevVertex := vm;
end;
end;

//*

pv := v[FinishVertex].PrevVertex;
Path.Clear;
Path.AddVertex (v[FinishVertex].x, v[FinishVertex].y, v[FinishVertex].z, v[FinishVertex].Tag);

while pv <> -1 do
if pv <> -1 then
begin
if v[pv].PrevVertex <> -1 then
Path.AddVertex (v[pv].x, v[pv].y, v[pv].z, v[pv].tag);

pv := v[pv].PrevVertex;
end;

//*
end;

constructor TGraph.Create;
begin
Path := TPath.Create;
end;

end.

Сладкая Смерть

6.12.2012 - 19:51

Волновой поиск

» Кликните сюда для просмотра оффтоп текста.. «

Поиск в глубину

» Кликните сюда для просмотра оффтоп текста.. «

Предыдущий алгоритм ещё называется поиском в ширину – он уровень за уровнем просматривает ближайшие клетки. Поиск в глубину выбирает одно направление и просматривает его вглубь на заданное число клеток, переходит к следующему направлению и т.д. пока не отыщется искомая точка. Чтобы найти кратчайший путь, надо вызвать функцию поиска для каждой из клеток, прилегающих к начальной. В вспомогательном массиве Mark такого же размера, как и основная карта хранится 1, если текущая клетка пройдена алгоритмом, и 0 – в противном случае. В переменных dest_x и dest_y должны храниться координаты точки, куда надо найти путь. В переменной length будет храниться длина текущего пути, чтобы он не стал длиннее максимального пути max_length. Вот функция поиска:

Цитата

int deepsearch(int x, int y){
int i;
if length>max_lenght return 0;
int Mark[x][y];
// Здесь забейте весь массив значениями 1
if (x=dest_x && y=dest+y) {
/* Точка найдена! Исходный путь представлен значениями 1 в массиве Mark. Здесь можно запомнить построенный путь и очистить массив Mark, чтобы продолжить поиск или же остановиться, возвратив 1, если найден хотя бы один путь. */
};
length=length+1;
if Mark[x+1][y]=0 deepsearch(x+1,y);
if Mark[x][y+1]=0 deepsearch(x,y+1);
if Mark[x+1][y+1]=0 deepsearch(x+1,y+1);
if Mark[x-1][y-1]=0 deepsearch(x-1,y-1);
if Mark[x-1][y]=0 deepsearch(x-1,y);
if Mark[x][y-1]=0 deepsearch(x,y-1);
length=length-1
};

Метод Хука-Дживса

» Кликните сюда для просмотра оффтоп текста.. «

Пастор

17.12.2012 - 6:29

Про двумерную упаковку: offline алгоритмы

Сегодня, дорогой Читатель, я расскажу тебе историю о комбинаторной оптимизации.
Издревле (как минимум, с начала прошлого века) математики задавались вопросом, как оптимально разместить некоторое количество ~~пива~~ нужных и полезных предметов в рюкзаке. Была сформулирована задача о ранце и ее подзадачи — тысячи их! — которые заинтересовали информатиков, криптографов и даже лингвистов.

От задачи о ранце отпочковалась задача об упаковке в контейнеры (Bin Packing Problem), одной из разновидностей которых является задача двумерной упаковки (2-Dimensional Bin Packing). Снова отбросив несколько вариаций, мы наконец придем к двумерной упаковке в полуограниченную полосу (2-Dimensional Strip Packing, 2DSP). Чувствуете, сколько интересного уже осталось за кадром? Но мы еще не закончили продираться сквозь классификацию. У 2DSP есть два варианта входных данных: когда набор упаковываемых объектов известен заранее (offline-проблема) и когда данные поступают порциями (online-проблема).

В этой статье рассматриваются алгоритмы решения offline-варианта 2DSP. Под катом немного матчасти и много картинок с цветными квадратиками.

В чем, собственно, проблема?

Как частный случай задачи двумерной упаковки, 2DSP заключается в упаковке объектов конкретной формы в конечное число контейнеров конкретной формы таким способом, чтобы число использованных контейнеров было наименьшим или количество или объем объектов (которые упаковывают) были наибольшими. Разница с двумерной упаковкой состоит в том, что есть только один контейнер, и вместо минимизации количества контейнеров, мы минимизируем высоту заполненности одного. Обеспечиваем максимальную плотность упаковки, если хотите.

Так как проблема NP-трудная, исследования сфокусированы главным образом на разработке приближенных алгоритмов решения. (На Хабре упоминался генетический алгоритм). Приближенные алгоритмы находят оптимальное решение с определенной точностью, но не гарантируют оптимальной упаковки для любого набора данных. При этом критерий оптимальности зависит от того, что вы пытались оптимизировать.
Я расскажу о простейших стратегиях и о том, к чему это все применимо в жизни (кроме рюкзака с пивом).

Итак, имеем набор из n прямоугольников и полуограниченный контейнер-стакан с фиксированной шириной W и бесконечной высотой. Каждый прямоугольник по ширине не превышает W. Задача — уложить прямоугольники в стакан без наложений и пересечений так, чтобы стакан стал как можно менее полон. Условимся, что все прямоугольники должны быть упакованы ортогонально, их нельзя поворачивать.


Исходные данные (начало XX века, кубизм)	Полуограниченная полоса	Вариант упаковки (похуже)	Вариант упаковки (получше)

Видно, что задача имеет идеальное решение, при котором высота упакованных прямоугольников равна их суммарной площади деленной на ширину полосы.

Зоопарк алгоритмов

Для offline варианта 2DSP сразу известен размер всех упаковываемых прямоугольников, поэтому их можно отсортировать, выбрать по определенному критерию, сгруппировать или просто натыкать в подходящие места. На этом и основываются большинство алгоритмов:

уровневые (level)
шельфовые (shelf)
плоские (plane)

Плоские алгоритмы размещают прямоугольники строго вплотную друг к другу, но это не самая удачная стратегия, как может показаться на первый взгляд. Уровневые делят полосу на уровни, равные по высоте одному из выбранных прямоугольников, и помещают все остальные на конкретный уровень по определенному критерию. Шельфовые предопределяют сразу несколько шельфов (полок), и распихивают по ним прямоугольники, такое поведение характерно для online-алгоритмов, а это уже совсем другая история.

Чем распространяться на общие слова, лучше обо всем по порядку.

Next Fit Decreasing High

Алгоритм, что называется, «в лоб». Прямоугольники сортируются по не-возрастанию высоты (Decreasing High намекает), самый высокий располагается в левом нижнем углу полосы, тем самым инициализируя первый уровень, по высоте равный ему. Остальные прямоугольники располагаются слева направо, пока есть место на текущем уровне. Прямоугольник, который не поместился на уровне, помещается сверху, образуя следующий уровень, и так далее.
Иллюстрации для каждого алгоритма будут производиться для следующих двух примеров: для наглядности, в левом форма прямоугольников тяготеет к вытянутой, в правом — более квадратные.

Алгоритм NFDH

Input: The number of rectangles to be packed n,
       the dimensions of the rectangles 
       {w(Li); h(Li)} and the strip width W.
Output: The height of a packing obtained in the strip.
 1: level = 0; h(level) = 0; w(level) = 0; i = 1
 2: Sort the rectangles in order of non-increasing height such that h(L1) ≥ h(L2) ≥ ... ≥ h(Ln)
 3: Pack rectangle Li left-justifed at the bottom of the strip
 4: h(level) = h(Li); w(level) = w(Li)
 5: for i = 2..n do
 6:   if W - w(level) ≥ w(Li) then
 7:     pack rectangle Li to the right of rectangle Li-1
 8:     w(level) += w(Li)
 9:   else [W - w(level) < w(Li)]
10:     create a new level above the previous one and pack rectangle Li on the new level
11:     level++; w(level) = w(Li); h(level) = h(level-1) + h(Li)
12:   end if
13: end for
14: print H = h(level)

First Fit Decreasing High

Похожий на предыдущий алгоритм, с той разницей, что для каждого следующего прямоугольника ищется место не только на последнем уровне, а начиная с самого нижнего. Отсюда и «first fit» — прямоугольник помещается на первый подходящий уровень снизу. Интуитивно понятно, что упаковка будет качественнее. Еще одно значительное отличие — в производительности, так как в худшем случае придется рассматривать на каждом шагу все уровни снизу вверх.

Алгоритм FFDH

Input: The number of rectangles to be packed n,
       the dimensions of the rectangles
       {w(Li); h(Li)} and the strip width W.
Output: The height of a packing obtained in the strip.
 1: level = 0; h(level) = 0; i = 1; LevelNum = 1
 2: Sort the rectangles in order of non-increasing height such that h(L1) ≥ h(L2) ≥ ... ≥ h(Ln)
 3: Pack rectangle Li left-justifed at the bottom of the strip; h(level+1) = h(Li)
 4: for i = 2..n do
 5:   search all levels (from the bottom) for the lowest with sufficient space
 6:   if such a level exist then
 7:     pack rectangle Li left justified on that level
 8:   else [there is insufficient space in all existing levels]
 9:     LevelNum++; level = LevelNum; h(level) = h(level-1) + h(Li)
10:     pack rectangle Li on the new level
11:   end if
12: end for
13: print H = h(level)

Best Fit Decreasing High

Модификация предыдущего алгоритма. Суть ее в том, что из уровней, подходящих для упаковки очередного прямоугольника, выбирается не первый, а лучший. Лучший уровень — это такой, на котором останется минимум места после упаковки текущего прямоугольника. Проще говоря, выбирается минимальное подходящее пространство, что способствует лучшему заполнению уровней.

Алгоритм BFDH

Input: The number of rectangles to be packed n,
       the dimensions of the rectangles
       {w(Li); h(Li)} and the strip width W.
Output: The height of a packing obtained in the strip.
 1: level = 0; h(level) = 0; i = 1; LevelNum = 1
 2: Sort the rectangles in order of non-increasing height such that h(L1) ≥ h(L2) ≥ ... ≥ h(Ln)
 3: Pack rectangle Li left-justifed at the bottom of the strip; h(level+1) = h(Li)
 4: for i = 2..n do
 5:   search all levels for the level with sufficient space and has minimum residual space
 6:   if such a level exist then
 7:     pack rectangle Li left justified on that level
 8:   else [there is insufficient space in all existing levels]
 9:     create a new level above the top-most level and pack rectangle Li
10:     LevelNum++; level = LevelNum; h(level) = h(level-1) + h(Li)
11:   end if
12: end for
13: print H = h(level)

Knapsack 0-1

Здесь стоит остановиться подробнее. Knapsack 0-1 — это частный случай задачи о ранце; примечателен тем, что кроме ответа на основной вопрос (нет, не 42, а полученный объем упаковки) дает еще и решение — список предметов, которые нужно упаковать. Порядок действий таков: прямоугольники сортируются по не-возрастанию высоты; упаковывается первый прямоугольник на новый уровень; для этого уровня находится решение задачи Knapsack 0-1, которое дает нам список прямоугольников, максимизированный по площади. Выбранные прямоугольники пакуются и извлекаются из списка неупакованных, уровень закрывается и открывается новый — инициализированный, как водится, первым (самым высоким) из оставшихся. Повторить, пока есть прямоугольники.

Алгоритм KP01

Input: The number of rectangles to be packed n,
       the dimensions of the rectangles
       {w(Li); h(Li)} and the strip width W.
Output: The height of a packing obtained in the strip.
 1: Sort the rectangles in order of non-increasing height such that h(L1) ≥ h(L2) ≥ ... ≥ h(Ln)
 2: level = 0
 3: while there are unpacked rectangles do
 4:   pack first unpacked rectangle, Li say
 5:   h(level) += h(Li)
 6:   solve KP01 instance
 7:   pack selected rectangles
 8:   level = level + 1
 9: end while
10: print H = h(level)

Split Fit

Алгоритм, призванный улучшить FFDN по принципу «разделяй и властвуй». Для начала отбираются прямоугольники, которые шире чем половина полосы. Они будут упакованы в первую очередь. Из них выбираются еще более широкие — шире, чем 2/3 полосы; они упаковываются алгоритмом FFDH. Над ними, начиная со следующего уровня (назовем его уровень R), упаковываются оставшиеся, те, что все еще шире 1/2, но уже 2/3. Они тоже упаковываются с помощью FFDH. Это разделение создает свободный регион шириной 1/3 справа от только что упакованных, начинающийся на уровне R и заканчивающийся текущей верхней границей упаковки (то есть он не простирается выше прямоугольников 1/2 < width <= 2/3). Все оставшиеся прямоугольники, которые уже, чем 1/2 полосы, упаковываются с помощью того же FFDH в первую очередь в образовавшийся регион, а если не помещаются — то сверху. На словах звучит громоздко, на картинке должно быть понятнее. А для тех, кто уже устал от моих литературных упражнений — псевдокод:

Алгоритм SF

Input: The number of rectangles to be packed n,
       the dimensions of the rectangles
       {w(Li); h(Li)} and the strip width W.
Output: The height of a packing obtained in the strip.
1: Let m ≥ 1 be the largest integer for which all rectangles in L
   have width at most 1=m
2: Partition the list of rectangles L into two sublists L1 and L2
   such that L1 is a list of rectangles of width greater than
   1/(m+1), while L2 is a list of rectangles of width at most 1/(m+1)
3: Pack the L1 rectangles into the strip, using the FFDH algorithm
4: Rearrange the blocks of this packing such that those of width
   greater than (m+1)/(m+2) are below those of width
   at most (m+1)/(m+2)
5: Pack rectangles of width at most 1/(m+2) into the region R
   using FFDH algorithm such that no rectangle overlaps the top
   of R and those failing to fit in R are packed above the packing of L1
6: Output the height of the strip, found by adding the height of each level

Join

Судя по всему, этот алгоритм был заточен под определенного характера входные данные, что ж, любая практическая ситуация имеет право на существование. Сейчас сами все поймете. Отсортированные, как водится, по не-возрастанию высоты прямоугольники объединяются в пары так, чтобы разница высоты в паре не превышала заданной доли, обычно 0-10%. Еще одно условие — чтобы их суммарная ширина помещалась в полосу. Полученные «супер-прямоугольники» упаковываются вместе с оставшимися без пары с помощью NFDH и FFDH, выбирается лучшее решение.
Существует вариация этого алгоритма, когда прямоугольники сортируются по ширине и объединяются вертикально, с тем же условием максимального отклонения ширины в паре на заданное количество процентов.

Алгоритм Join

Input: The number of rectangles to be packed n,
       the dimensions of the rectangles
       {w(Li); h(Li)}, the constant gamma as a
       percentage and the strip width W.
Output: The height of a packing obtained in the strip.
 1: Sort the rectangles in order of non-increasing height such that h(L1) ≥ h(L2) ≥ ... ≥ h(Ln)
 2: j = 1
 3: while j+1 ≤ n do
 4:   if (h(Lj) - h(Lj+1))/h(Lj) * 100 < gamma
         and w(Lj) + w(Lj+1) ≤ W then
 5:     w(Lj) += w(Lj+1)
 6:     j += 2
 7:   else
 8:     j++
 9:   end if
10: end while
11: Execute the NFDH and FFDH algorithms to pack the rectangles
12: From the best solution, construct a feasible packing of the original instance
13: Output the height of the strip, found by adding the height of each level

Floor Сeiling No Rotation

Если вы все еще недоумеваете по поводу остающегося на уровнях пространства, то этот алгоритм для вас. Прямоугольники сортируются по не-возрастанию высоты (неожиданно, да?) и применяется алгоритм BFDH с некоторыми модификациями. У каждого уровня есть «пол» (floor) и «потолок» (ceiling). Пока есть возможность, прямоугольники пакуются на «пол» слева направо. Когда место заканчивается, предпринимается попытка упаковать на «потолок» справа налево; если нет места на потолке, то только тогда начинается новый уровень. В лучших традициях BFDH, на каждом шагу просматриваются все уровни — сначала «пол», затем «потолок» — на наличие наиболее подходящего места. Упаковка, как видите, неплохая. Метод удачно упаковывает по «потолкам» самые мелкие прямоугольники.

Алгоритм FCNR

Input: The number of rectangles to be packed n,
       the dimensions of the rectangles
       {w(Li); h(Li)} and the strip width W.
Output: The height of a packing obtained in the strip.
 1: Sort the rectangles in order of non-increasing height such that h(L1) ≥ h(L2) ≥ ... ≥ h(Ln)
 2: for i = 1..n do
 3:   if Li is ceiling feasible then
 4:     pack Li on ceiling with minimum residual space
 5:   else [Li is not ceiling feasible]
 6:     if Li is floor feasible then
 7:       pack Li on the floor with minimum residual space
 8:     else [Li is not floor feasible]
 9:       level++;
10:     end if
11:   end if
12: end for
13: Output the height H of the strip, found by adding the height of each level

Sleator

Настало время «плоских» алгоритмов, без деления на уровни. Алгоритм Sleator использует интуитивный принцип упаковки рюкзака: сложить на дно самые крупные предметы и засыпать сверху мелкими. Вот как это выглядит. Из прямоугольников выбираются самые широкие, шире половины полосы, как вы уже догадались, и хаотично укладываются друг на друга с выравниванием по левому краю. Оставшиеся сортируются по не-возрастанию высоты и начинают укладываться друг за другом слева направо сверху на уже уложенные. Как только их суммарная ширина превысит половину ширины полосы, оставшиеся раскидываются друг на друга то слева (начиная от левого края полосы), то справа (от середины), по принципу «где на данный момент меньше». Как видно из рисунков, этим методом удобнее укладывать книги, чем коробки.

Алгоритм Sleator

Input: The number of rectangles to be packed n,
       the dimensions of the rectangles
       {w(Li); h(Li)} and the strip width W.
Output: The height of a packing obtained in the strip.
 1: Partition L into two sublists L1 and L2 consisting of rectangles
    of width greater than 1/2 and at most 1/2 respectively.
 2: Stack all the rectangles in L1 left justified on top of
    one another starting at the bottom of the strip. Compute Hstack
 3: Packing will continue above Hstack
 4: Sort the rectangles in L2 according to non-increasing height
    such that h(Li) ≥ h(Li+1) for i < n
 5: Let Htall be the height of the tallest rectangle in list L2.
 6: Pack the rectangles, left justified from the left to the right
    edge of the strip until there is insufficient space to pack
    a rectangle or all of the rectangles have been packed
 7: Partition the strip with a vertical line into two equal segments.
    There is possibly one rectangle whose interior may be intercepted
    by the vertical line.
 8: Let Hright and Hleft be the height of the rectangle on the right
    (resp. left) half of the strip whose left (resp. right) edge
    is adjacent to the vertical line or whose interior is intercepted
    by the vertical line
 9: while there are unpacked rectangles do
10:   Draw horizontal lines of length half across
      the rectangles whose height is Hleft and Hright
11:   All subsequent packing will either be on the left
      or right segment of the strip
12:   Select the segment with minimum height and pack
      rectangles from the edge of the strip to the
      vertical line until all rectangles have been
      packed or there is a rectangles which does not fit
13: end while
14: print H = max{Hleft; Hright}

Burke

Снова безуровневый алгоритм, для которого вводится дополнительная «карта высот»:

|               |    |    _   _   ___|
|               |    |   | |_|_| |   |
|               |    |_  | |   | |   |
|_ _ _ _ _ _ _ _| .. |_|_|_|_ _|_|_ _|
 0 0 0 0 0 0 0 0      1 0 3 2 3 0 3 3

Это массив, по которому по мере заполнения полосы можно отслеживать наименее заполненные области и их ширину. В начале он заполнен нулями. Прямоугольники сортируются по не-возрастанию, внезапно, ширины. Затем, на каждом шаге алгоритма:
1. Вычисляется позиция самой низкой области — индекс минимального значения массива;
2. Выбирается наиболее подходящий прямоугольник — во-первых, помещающийся в эту область, во-вторых, максимально ее заполняющий по ширине;
3. Если подходящий прямоугольник найден, он размещается в этой области одним из способов:
3.1 По левому краю области;
3.2 Ближе к более высокому соседу, если один из соседей — край полосы, то ближе к краю;
3.3 Ближе к менее высокому соседу, если один из соседей — край полосы, то дальше от края. К значениям массива, соответствующим ширине прямоугольника, прибавляется его высота.
4. Если подходящего прямоугольника нет, область «заполняется» путем выравнивания ее высоты до высоты ближайшего края.
Вам уже захотелось написать по этому алгоритму бота, играющего в Тетрис?

Алгоритм Burke

Input: The number of rectangles to be packed n,
       the dimensions of the rectangles
       {w(Li); h(Li)} and the strip width W.
Output: The height of a packing obtained in the strip.
 1: Sort the rectangles according to non-increasing width such that w(Li) ≥ w(Li+1) ≥ .. ≥ W(Ln)
 2: for each placement policy
    (leftmost, tallest neighbour, smallest neighbour) do
 3:   while Rectangles not packed do
 4:     find lowest gap
 5:     if w(Li) ≤ GapWidth then
 6:       place best fitting rectangle using placement policy
 7:       raise elements of array to appropriate height
 8:     else
 9:       raise gap to height of the lowest neighbour
10:     end if
11:   end while
12: end for
13: The elements of the array with greatest entry give the total height of the packing
14: Compare total packing heights obtained by each placement policy and return the best solution

Кто лучше?

Результат работы каждого алгоритма — это уровень заполненности полосы, чем меньше, тем лучше. Его можно оценить отношением к нему оптимального:

Сравним полученные результаты для двух приведенных наборов исходных данных:

	Оптимальное решение	NFDH	FFDH	BFDH	KP01	SF	JOIN	FCNR	Sleator	Bruke
Набор 1	149	0,65	0,71	0,71	0,71	0,75	0,61	0,83	0,68	0,72
Набор 2	140	0,66	0,77	0,77	0,78	0,77	0,70	0,84	0,51	0,71

Видим, что для обоих наборов данных победил алгоритм Floor Ceiling. Можно отметить, что Sleator показал гораздо лучший результат для первого набора, а Join, наоборот, больше подходит для второго. Но все это не более, чем статистика.

Эпилог

За счет своей простой формулировки, задача бинарной упаковки может быть подтянута к огромному количеству практических применений: непосредственно к упаковке объектов или раскройке материала, распределению средств, планировке задач на кластерах, etc. В жизни мало рафинированных задач, но, ограничив некоторые условия, можно свести ситуацию к красивой математической модели, имеющей неплохое решение, достижимое за полиномиальное время. А однажды, в результате эффекта бабочки, отброшенные условия возымеют большой вес и решение полетит ко всем Безымянным. Вот из-за таких допущений и несовершенен мир.

Что-то я отвлеклась. В следующей части надеюсь рассказать об online-версии задачи и шельфовых алгоритмах. Good luck!

Источники вдохновения:
Nthabiseng Ntene An Algorithmic Approach to the 2D Oriented Strip Packing Problem
David Pisinger Knapsack problem
Survey on two-dimensional packing
Level Algorithms and Shelf Algorithms (осторожно, дизайн-вырви-глаз)

Код (Qt):
Алгоритмы packager.h packager.cpp
Гуй window.h window.cpp renderarea.h renderarea.cpp main.cpp

Wildy

Пастор

17.12.2012 - 6:45

Про двумерную упаковку: online алгоритмы

Это продолжение поста про оффлайн алгоритмы упаковки.

Суть задачи: имеем полубесконечную полосу — как в тетрисе, только без game over'а, и конечный набор прямоугольников. Данные о прямоугольниках поступают в режиме реального времени; каждый новый прямоугольник необходимо немедленно разместить и больше не двигать с места. Цель — минимизировать общую высоту упакованных прямоугольников.
Это online-вариация задачи об упаковке прямоугольников в полуограниченную полосу (2 Dimensional Strip Packing, 2DSP).

Чуть больше теоретических сведений можно найти в предыдущей статье, а пока, без лишних слов, перейдем к алгоритмам.

Для иллюстрации будет использоваться вот такая последовательность кирпичей:

(Спонсор исходных данных — random.org).
Будут описаны некоторые эвристические алгоритмы, общая идея которых состоит в том, чтобы предварительно разделить полосу на области, и размещать вновь поступающие прямоугольники в конкретной области по определенным критериям.

Level algorithms

Подход всех уровневых алгоритмов в том, что полоса разделяется на уровни, исходя из высоты имеющихся на данном этапе прямоугольников. Прямоугольники размещаются вдоль основания текущего уровня слева направо, пока это возможно; прямоугольник, который не поместился, упаковывается на следующий уровень. Высота каждого уровня определяется по самому высокому прямоугольнику в нем. Вот три варианта упаковки:
Next Fit Level — когда открывается следующий уровень, предыдущий «закрывается» и его больше не рассматриваем;
First Fit Level — на каждом шагу алгоритма просматривается каждый уровень, начиная с самого нижнего, и прямоугольник упаковывается на первый подходящий, на котором достаточно места;
Best Fit Level — на каждом шагу алгоритма просматриваются все уровни, и выбирается наиболее подходящий — тот, на котором после упаковки останется минимум места.
Все три алгоритма проиллюстрированы выше. На этом этапе несложно догадаться, какой где.

Алгоритмы NFL, FFL, BFL

Next Fit Level

Input: The dimensions of the rectangles {w(Li); h(Li)} and the strip width W.
Output: The height H and the entire packing.
 1: level = 0; h(level + 1) = 0; H = 0; i = 0
 2: while there is an unpacked rectangle do
 3:    i++
 4:    if there is sufficient space on current level then
 5:       pack rectangle left justified
 6:       if h(level + 1) < h(Li) then
 7:          h(level + 1) = h(Li)
 8:       end if
 9:    else [there is insufficient space on current level]
10:       open a new level and pack the rectangle
11:       H += h(level + 1); level++
12:    end if
13: end while
14: print H and entire packing

First Fit Level

Input: The dimensions of the rectangles {w(Li); h(Li)} and the strip width W.
Output: The height H and the entire packing.
 1: level = 0; h(level + 1) = h(L1); H = h(L1); i = 1
 2: while there is an unpacked rectangle do
 3:    i++; level = 0
 4:    search for the lowest level with sufficient space
 5:    if such a level exists then
 6:       pack rectangle left justified
 7:       if this is the top-most level and h(level) < h(Li) then
 8:          h(level) = h(Li)
 9:       end if
10:    else [there is insufficient space in all existing levels]
11:       create a new level above the top-most level
12:       pack rectangle left justified
13:       h(level) = h(Li); H += h(Li)
14:    end if
15: end while
16: print H and entire packing

Best Fit Level

Input: The dimensions of the rectangles {w(Li); h(Li)} and the strip width W.
Output: The height H and the entire packing.
 1: level = 0; h(level + 1) = h(L1); H = h(L1); i = 0
 2: while there is an unpacked rectangle do
 3:    i++; level = 0
 4:    search for lowest level with minimum residual horizontal space
 5:    if such a level exists then
 6:       pack rectangle left justified
 7:    else [such a level does not exist]
 8:       create a new level above the top-most level
 9:       pack rectangle left justified
10:       h(level) = h(Li); H += h(Li)
11:    end if
12: end while
13: print H and entire packing

Bi-Level Next Fit

Хитрая модификация Next Fit Level. Уровни создаются по два, на каждом из которых своя стратегия.
Нижний уровень (или все нечетные): первый пришедший прямоугольник упаковывается по левому краю, остальные — справа налево, пока достаточно места. Затем уровень закрывается и его высота определяется по самому высокому прямоугольнику в нем.
Верхний уровень (или все четные): если на нижнем уровне всего один прямоугольник, уровень упаковывается аналогично нижнему. Если два и больше, то первый прямоугольник упаковывается над более низким из двух прижатых к краям полосы, и также прижимается к краю полосы; все последующие упаковываются справа налево, независимо от того, где был упакован первый — слева или справа.

Звучит путано и не сразу понятно, к чему такие танцы с бубном. Единственное, что очевидно обеспечивает этот алгоритм — это более равномерное распределение прямоугольников по объему полосы, без перекоса к левому краю. Но мы еще вернемся к этой стратегии, когда будем рассматривать алгоритмы компрессии.

Алгоритм BiNFL

Bi-Level Next Fit

Input: The dimensions of the rectangles {w(Li); h(Li)} and the strip width W.
Output: The height H and the entire packing.
 1: level = 1; h(level) = 0; H = 0; i = 1
 2: open a new bi-level
 3: call 'pack on lower level'
 4: print H and entire packing
Procedure 'pack on lower level'
 1: the first rectangle packed is left justified
 2: subsequent rectangles that fit are right justified
 3: if there is insufficient space then
 4:    h(level) is given by the height of tallest rectangle packed
 5:    call 'pack on upper level'
 6: end if
Procedure 'pack on the upper level'
 1: if Li+1 is the first rectangle packed then
 2:    left justify on top of Li
 3: else
 4:    if Li+2 is the first rectangle packed then
 5:       pack above min {h(Li); h(Li+1)}
 6:       justify according to min {h(Li); h(Li+1)}
 7:   else
 8:      if Li+j (j > 3) fits on the upper level then
 9:         left justify all other rectangles
10:      else [rectangle does not fit]
11:         H += h(lower) + h(upper); open new bi-level
12:      end if
13:    end if
14: end if

Shelf algorithms

Шельфовые алгоритмы не привязываются к точной высоте прямоугольников и после упаковки первого на уровне (шельфе) остается немного места на случай если следующий будет немного выше. Это «немного» регулируется параметром r ∈ (0;1). Высота шельфа определяется как r^k, где r^k+1 < h(L_i) ≤ r^k для некоторого целого k и высоты первого упаковываемого прямоугольника h(L_i). Аналогично уровневым алгоритмам, применяются стратегии обхода шельфов Next Fit, First Fit и Best Fit. Оставленный на каждом шельфе «резерв» становится выгодным в двух последних случаях (сравните FFL и FFS, BFL и BFS). В примере p=0,7.

Алгоритмы NFS, FFS, BFS

Next Fit Shelf

Input: The dimensions of the rectangles {w(Li); h(Li)}
       the parameter r and the strip width W.
Output: The height H and the entire packing.
 1: shelf = 1; w(shelf) = 0;
    h(shelf) = r^k for some integer k; i = 0; H = h(shelf)
 2: while there is an unpacked rectangle do
 3:    i++; compute k
 4:    if there is sufficient space and r^(k+1) < h(Li) ≤ r^k then
 5:       pack rectangle to the right of the previously packed rectangle
 6:    else [there is insufficient space]
 7:       create a new shelf on top of the previous one
          and pack rectangle Li on the new shelf.
 8:       shelf++; H += r^k
 9:    end if
10: end while
11: print H and entire packing.

First Fit Shelf

Input: The dimensions of the rectangles {w(Li); h(Li)}
       the parameter r and the strip width W.
Output: The height H and the entire packing.
 1: shelf = 1; h(shelf) = r^k for some integer k;
    i = 0; H = h(shelf); shelfnum = 1
 2: while there is an unpacked rectangle do
 3:    i++; compute k
 4:    search for lowest shelf of appropriate height with sufficient space
 5:    if such a shelf exists then
 6:       pack rectangle to the right of the previously packed rectangle
 7:    else [there is insufficient space in all shelves of
             appropriate height or such a shelf does not exist]
 8:       create a new shelf above the top-most shelf
          and pack rectangle Li on the new shelf.
 9:       shelf++; H += r^k; shelfnum++
10:    end if
11: end while
12: print H and entire packing.

Best Fit Shelf

Input: The dimensions of the rectangles {w(Li); h(Li)}
       the parameter r and the strip width W.
Output: The height H and the entire packing.
 1: shelf = 1; h(shelf) = r^k for some integer k;
    i = 0; H = h(shelf); shelfnum = 1
 2: while there is an unpacked rectangle do
 3:    i++; compute k
 4:    search all shelves (starting with the bottom) for one
       with appropriate height and has minimum residual horizontal space.
 5:    if such a shelf exists then
 6:       pack rectangle to the right of the previously packed rectangle
 7:    else [there is insufficient space or
             there is no shelf of appropriate height]
 8:       create a new shelf above the top-most shelf
          and pack the rectangle Li on the new shelf.
 9:       shelf++; H += r^k; shelfnum++
10:    end if
11: end while
12: print H and entire packing.

Harmonic Shelf

Harmonic Shelf пакует на один уровень прямоугольники не только со схожей высотой, но и со схожей шириной. В результате получится больше уровней, но заполнены они будут плотнее. Вводится дополнительный параметр M для расчета допустимой ширины.

Общая ширина полосы, для простоты возьмем ее за единицу, делится на M интервалов I₁..I_M. Автор утверждает, что разумное значение M находится в [3; 12]. Ширина интервалов рассчитывается по формуле:

;

Для каждого прямоугольника рассчитывается значение p. Параметр k для высоты рассчитывается аналогично шельфовым алгоритмам. Проверяются два условия: есть ли для прямоугольника место на каком-нибудь шельфе высотой r^k и соответствует ли его p уже упакованным там. Если хотя бы одно не выполнено, для него создается свой шельф, с блекджеком и максимально комфортными условиями. В примере p=0,7, M=4.

Алгоритм HS

Harmonic Shelf

Input: The dimensions of the rectangles {w(Li); h(Li)}
       the parameters M, r and the strip width W.
Output: The height H and the entire packing.
 1: shelf = 1; h(shelf) = r^k for some integer k; i = 0
 2: while there is an unpacked rectangle do
 3:    i++; compute k such that r^(k+1) < h(Li) ≤ r^k
 4:    compute the value of p such that 1/(p+1) < w(Li) ≤ 1/p
       where 1 ≤ p < M
 5:    amongst shelves of height r^k find the one
       with packing rectangles whose width fits in the interval Ip
 6:    if there is insufficient space or no rectangle of height r^k then
 7:       a new shelf of appropriate height is created above the top-most shelf
 8:       shelf++
 9:    end if
10:    if rectangle Li is the first on the shelf then
11:       h(shelf) = r^k; H += h(shelf)
12:    end if
13: end while
14: print H and entire packing

Azar_Y

Для деления на уровни вводится некое пороговое значение 0 < Y < 1/2. Ширина прямоугольников масштабируется в соответствии с шириной полосы, которая принимается за единицу. Прямоугольники высотой 2^j-1 < h(L_i) ≤ 2^j и шириной 2^-x-1 < h(L_i) ≤ 2^-x упаковываются на один уровень. То есть уровень характеризуется парой (x, j) для некоторого целого j и натурального x.

Прямоугольники шириной не менее Y нарекаются буферными и пакуются каждый на отдельный уровень. То есть если попался такой прямоугольник, нужно срочно открывать новый уровень, высотой равный ему. Все остальные (не-буферные) пакуются на первый подходящий уровень. Подходящий — в смысле с такой же парой (x, j). Если вдруг места не хватило, не-буферный прямоугольник может стать первым в новом уровне, тогда высота этого уровня будет 2^j.
Какая-то нездоровая иерархия, в общем.
В примере Y=0,4.

Алгоритм Azar

Azar_Y

Input: The dimensions of the rectangles {w(Li); h(Li)}
       the parameter Y and the strip width W.
Output: The height H and the entire packing.
 1: h(level) = 0; w(level) = 0; level = 0
 2: while there is an unpacked rectangle do
 3:    if w(Li) ≥ Y then
 4:       level++; w(level) = w(Li); h(level) = h(Li); H += h(level)
 5:    else [w(Li) < Y ]
 6:       compute x and j such that
          2^(j-1) < h(Li) ≤ 2^j and 2^(-x-1) < w(Li) ≤ 2^(-x)
 7:       search for the lowest (x, j) level
 8:       if no (x, j) level is available or Li is blocked then
 9:          level++; h(level) = 2^j; H += h(level)
10:          w(level) += w(Li)
11:       else [(x, j) level is available and not blocked]
12:          w(level) += w(Li)
13:       end if
14:    end if
15: end while
16: print H and entire packing

Compression algorithms

Компрессионные алгоритмы основаны на BiNFL и по сути представляют собой патчи к методу упаковки верхнего уровня.
Общее отличие — верхний уровень всегда упаковывается слева направо.
Compression Part Fit: Для каждого прямоугольника, упакованного на верхний уровень, проверяется, есть ли под ним место на нижнем уровне. Если можно сдвинуть его на нижний уровень так, что часть его останется на верхнем, он сдвигается (отсюда Part Fit — он частично может быть «вжат» в предыдущий уровень). Таким образом можно снизить общую высоту второго уровня. Больше на упаковку он никак кардинально не влияет. На первом рисунке видно прямоугольник, висящий между уровней — он и был сдвинут.
Compression Full Fit: Для каждого прямоугольника, упакованного на верхний уровень, проверяется, может ли он полностью сместиться вниз на предыдущий уровень. Собственно, если может, то и смещается. Это дает нам стратегическое преимущество: освободившееся таким образом место на уровне может быть использовано для следующего прямоугольника. К сожалению, второй рисунок не содержит такую ситуацию, но зато на нем можно увидеть двоих провалившихся.
Compression Combo: Да-да, это комбинация двух предыдущих и, соответственно, комбинация их плюсов. Третий рисунок иллюстрирует.

Алгоритмы CPF, CFF, CC

Compression

Input: The dimensions of the rectangles {w(Li); h(Li)} and the strip width W.
Output: The height H and the entire packing.
 1: i = 0; H = 0
 2: open new bi-level
 3: packing on lower level is similar to BiNFL; H += h(lowerlevel)
 4: if there is insufficient space to pack a rectangle on the lower level then
 5:    call 'pack on the upper level'
 6: end if
 7: print H and entire packing
Procedure 'pack on the upper level'
 1: while there is an unpacked rectangle do
 2:    if i ≥ 3 and Li is the first rectangle packed then
 3:       justify according to the shorter
          of the left-most and right-most rectangles
 4:       slide rectangle downwards if there is sufficient space
          and go to 12
 5:    else
 6:       if i ≥ 3 and Li is the second rectangle packed then
 7:          right justify and go to 12
 8:       else [rectangle does not fit]
 9:          H += h(upperlevel); open new bi-level
10:       end if
11:    end if
12: end while

Compression Part Fit

Input: The dimensions of the rectangles {w(Li); h(Li)} and the strip width W.
Output: The height H and the entire packing.
 1: i = 0; H = 0
 2: open a new bi-level
 3: packing on lower level is similar to BiNFL; H += h(lowerlevel)
 4: if there is insufficient space to pack a rectangle on the lower level then
 5:    call 'pack on the upper level'
 6: end if
 7: print H and entire packing
Procedure 'pack on the upper level'
 1: while there is an unpacked rectangle do
 2:    if h(Li) > VerticalSpace and w(Li) ≤ HorizontalSpace then
 3:       shift rectangle Li downwards; go to 11
 4:    else
 5:       if h(Li) ≤ VerticalSpace and W - w(level) ≥ w(Li) then
 6:          left justify; go to 11
 7:       else [rectangle does not fit]
 8:          H += h(upperlevel); open a new bi-level
 9:       end if
10:    end if
11: end while

Compression Full Fit

Input: The dimensions of the rectangles {w(Li); h(Li)} and the strip width W.
Output: The height H and the entire packing.
 1: i = 0; H = 0
 2: open a new bi-level
 3: packing on lower level is similar to BiNFL; H += h(lowerlevel)
 4: if there is insufficient space to pack a rectangle on the lower level then
 5:    call 'pack on the upper level'
 6: end if
 7: print H and entire packing
Procedure 'pack on the upper level'
 1: while there is an unpacked rectangle do
 2:    if h(Li) ≤ VerticalSpace and w(Li) ≤ HorizontalSpace then
 3:       slide rectangle Li downwards; go to 11
 4:    else
 5:       if h(Li) ≤ VerticalSpace and W - w(level) ≥ w(Li) then
 6:          left justify; go to 11
 7:       else [rectangle does not fit]
 8:          H += h(upperlevel); open new bi-level
 9:       end if
10:    end if
11: end while

Compression Combo

Input: The dimensions of the rectangles {w(Li); h(Li)} and the strip width W.
Output: The height H and the entire packing.
 1: i = 0; H = 0
 2: open a new bi-level
 3: packing on lower level is similar to BiNFL; H += h(lowerlevel)
 4: if there is insufficient space to pack a rectangle on the lower level then
 5:    call 'pack on the upper level'
 6: end if
 7: print H and entire packing
Procedure 'pack on the upper level'
 1: while there is an unpacked rectangle do
 2:    if w(Li) ≤ HorizontalSpace then
 3:       slide rectangle Li downwards; go to 11
 4:    else
 5:       if w(Li) > HorizontalSpace and W - w(level) ≥ w(Li) then
 6:          left justify, go to 11
 7:       else [rectangle does not fit]
 8:          H += h(upperlevel); open new bi-level
 9:       end if
10:    end if
11: end while

Online fit

Более трудоемкий метод, чем предыдущие, но и более производительный. Основан на алгоритме Burke для оффлайн-варианта задачи. Некоторые аспекты я, пожалуй, нарисую.

В процессе упаковки могут образовываться пустые области (Empty Areas), которые алгоритм будет пытаться заполнить в первую очередь. Как они могут образоваться, видно на фото слева. Заметим, что для этого нужно хранить некое представление заполненности полосы.
После заполнения из пустой области может получиться две (а может и не получиться, you know).

Также одним шагом упаковки могут быть созданы сразу несколько областей — пространство под прямоугольником анализируется по вертикали и делится по количеству «ступенек». На среднем фото еще раз показано разбиение пустой области на две после упаковки туда прямоугольника (кстати, мне кажется, в этом месте оптимальнее было бы разбить ее по горизонтали).
В общем, чем дальше в лес, тем больше пустот. На каждом шаге пристально рассматриваются стремительно мельчающие пустоты. Здесь мне видится еще одна оптимизация: можно ввести меру целесообразной площади и отбрасывать неугодные.

Если ни одна пустая область не подходит, упаковка продолжается методом Burke. Я напомню. Строится карта высот для полосы, и прямоугольник примеривается в самую низкую на данный момент область. Если он туда не входит, низкое место «заполняется» до уровня ближайшей ступеньки, и поиск повторяется. Так потенциальное место может быть вытолкнуто на самый верх, и тогда возникает еще один вариант образования пустого пространства (см. справа).
Та же карта высот используется на каждом шаге для выявления пустых областей.

В общем-то, все уже сказано. Сперва пытаемся упаковать по пустым областям, затем в наиболее низкое место, и в последнюю очередь — сверху, вырожденный случай «самого низкого места».

Алгоритм OF

Online Fit

Input: The dimensions of the rectangles {w(Li); h(Li)} and the strip width W.
Output: The height H and the entire packing.
 1: i = 0; H = 0
 2: Create a linear array whose number of elements are equal to W
 3: while there is an unpacked rectangle do
 4:    i++
 5:    Search the list of empty areas for sufficient space to pack Li
 6:    if w(EmptyArea) ≥ w(Li) and h(EmptyArea) ≥ h(Li) then
 7:       pack rectangle in empty area
 8:    else
 9:       if rectangle does not fit in any of the empty areas then
10:          search the linear array for available packing width
11:       else [there is insufficient space in linear array]
12:          pack rectangle on top left justified
             from the index leading to smaller space
13:       end if
14:    end if
15: end while
16: H = highest entry in the linear array
17: print H and entire packing

Послесловие

Чтобы развеять ощущение бесцельности происходящего, приведу один из классических, но вдохновляющих примеров применения алгоритмов двумерной упаковки в полуограниченную полосу. Это планировщик задач для многопроцессорной системы. В современном софте для кластеров используются более «взрослые» алгоритмы, но сама задача планирования сводится к 2DSP. В условиях работы кластеров поток поступающих задач представляет собой не что иное, как online данные для упаковки: размерность по горизонтали — требуемое количество ядер, по вертикали — время работы задачи. (Все как на заре технологий: встаешь в очередь к компьютеру и заранее предупреждаешь, сколько тебе нужно ресурсов и как долго будет крутиться твоя перфокарта). Планировщик должен выделить для задачи конкретные процессоры и определить момент запуска. В чистом виде, кстати, ни один алгоритм использовать нельзя, так как время капает, пока планировщик ждет задачу или думает над ее размещением. Дополнительно программисту могут связывать руки разные специфические условия: приоритет задач; возможность масштабирования задачи во время выполнения; хардварные задержки и самообслуживание системы, в том числе перераспределение с учетом отказавших ресурсов; сопоставление характера коммуникаций внутри задачи с архитектурой кластера, etc.
Любой алгоритм прекрасен в своей математической строгости, пока не начнешь прикладывать его к реальной жизни.

И напоследок — сравнительная характеристика алгоритмов в приятном обезличенном виде. Здесь приведено отношение оптимальной высоты упаковки (сумма площадей прямоугольников, деленная на ширину полосы) к полученной.

NFL	FFL	BFL	BiNFL	NFS	FFS	BFS	HS	Azar_Y	CPF	CFF	CC	OF
0,56	0,63	0,75	0,56	0,45	0,57	0,57	0,37	0,44	0,60	0,59	0,63	0,72

Код (Qt):
Алгоритмы packager.h packager.cpp
Гуй window.h window.cpp renderarea.h renderarea.cpp main.cpp

Wildy

В чем, собственно, проблема?

Зоопарк алгоритмов

Next Fit Decreasing High

Алгоритм NFDH

First Fit Decreasing High

Алгоритм FFDH

Best Fit Decreasing High

Алгоритм BFDH

Knapsack 0-1

Алгоритм KP01

Split Fit

Алгоритм SF

Join

Алгоритм Join

Floor Сeiling No Rotation

Алгоритм FCNR

Sleator

Алгоритм Sleator

Burke

Алгоритм Burke

Кто лучше?

Эпилог

Level algorithms

Алгоритмы NFL, FFL, BFL

Next Fit Level

First Fit Level

Best Fit Level

Bi-Level Next Fit

Алгоритм BiNFL

Shelf algorithms

Алгоритмы NFS, FFS, BFS

First Fit Shelf

Best Fit Shelf

Harmonic Shelf

Алгоритм HS

AzarY

Алгоритм Azar

Compression algorithms

Алгоритмы CPF, CFF, CC

Compression Part Fit

Compression Full Fit

Compression Combo

Online fit

Алгоритм OF

Послесловие

Azar_Y