Модераторы: maxim1000

Поиск:

Закрытая темаСоздание новой темы Создание опроса
> Коллекция алгоритмов от Johna Smith, в качестве учебного материала 
:(
    Опции темы
podval
Дата 22.12.2004, 21:37 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



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


Binary insertions

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Binary insertion (array sorting)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Split our array into two parts - sorted (left) and unsorted (right)
//    Initially sorted part contains only one element - first array element
//    Take an element from the unsorted part and put it in the right place in
//    the sorted part (using binary search to find this right place). So
//    sorted part contains now two elements. Repeat this process until
//    no elements left in the unsorted part.
//
//////////////////////////////////////////////////////////////////////////////

#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 element; // value of the element we want to insert
unsigned int right,left,middle; // auxulary variables for binary search

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

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

  // Sorting
  for (int i=1;i<N;i++) // main loop
  {
    // Searching place for element i (binary search)
    left=0;
    right=i;
    element=array[i];
    while (left<right)
    {
      middle=(left+right)/2;
      if (array[middle]<=element) left=middle+1; else right=middle;
    }

    // Inserting element i into the right place
    for (int j=i;j>right;j--) array[j]=array[j-1];
    array[right]=element;
  }

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



Straight insertion

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Straight insertion (array sorting)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Split our array into two parts - sorted (left) and unsorted (right)
//    Initially sorted part contains only one element - first array element
//    Take an element from the unsorted part and put it in the right place in
//    the sorted part. So sorted part contains now two elements. Repeat
//    this process until no elements left in the unsorted part.
//
//////////////////////////////////////////////////////////////////////////////

#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 element; // value of the element we want to insert
unsigned int index; // index of the right position of the element

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

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

  // Sorting
  for (int i=1;i<N;i++) // main loop
  {
    index=0;
    // Searching place for element i
    while (array[index]<array[i]) index++;
    // Inserting element i into the right place
    element=array[i];
    for (int j=i;j>index;j--) array[j]=array[j-1];
    array[index]=element;
  }

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



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

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Straight selections (array sorting)
//  (б) Johna Smith, 1996
//
//  Method description:
//    searching for smallest element in array, swapping with first element,
//    repeat this operation starting search from second element and so on
//    array become sorted when we'll start search from (n-1)-th element
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

const N = 10;  // number of elements in array

int array[N]={100,15,234,11,63,78,4,200,0,4};
int min; // smallest element in array
int index; // index of the smallest element in array
int swp; // auxulary variable for swapping

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

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

   // Sorting
   for (int i=0;i<N-1;i++) // main loop
   {
     // searching for smallest element from i-th
     min=array[i]; // setting i-th elementas smallest
     index=i;
     for (int j=i+1;j<N;j++)
       if (array[j]<min)  // if there is element less than smallest
       {
          min=array[j];  // then remember its value
          index=j;       // and index
       }
     // swapping i-th and smallest element
     swp=array[index];
     array[index]=array[i];
     array[i]=swp;
   }

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



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

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Shell sort (array sorting)
//  (c) Johna Smith, 1996
//
//  Method description:
//    1) Sort elements (0,4,8,12,...) (step=4) using straight insertions method
//    2) Sort elements (0,2,4,6,8,...) (step=2) using straight insertions method
//    3) Sort all elements (step=1) using straight insertions method
//    Shell sort is faster than simple straight insertions because at each
//    step more and more elements are already sorted (productivity increased
//    by ~300% for random array of 256 elements & about 700% (!) for random array
//    of 2048 elements (N. Wirth, Algoritms and data structures)
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

const T = 4;   // number of distances
const N = 20;  // array size

int array[N]={100,15,234,11,63,78,4,200,0,4,14,32,44,58,1,2,3,9,8,7}; // array of N integers
int interval[T]={9,5,3,1}; // last interval must be 1
int step; // current step
int element; // value of the element we want to insert
unsigned int index; // index of the right position of the element

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

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

  // Sorting
  for (int m=0;m<T;m++)  // main loop (by all distances)
  {
    step=interval[m];
    // sorting (0,step,2*step,3*step,...) elements using straight insertions
    for (int i=0;i<N;i+=step)
    {
      index=0;
      // Searching place for element i
      while (array[index]<array[i]) index+=step;
      // Inserting element i into the right place
      element=array[i];
      for (int j=i;j>index;j-=step) array[j]=array[j-step];
      array[index]=element;
    }
  }

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



Quick sort (recursive)

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Quick sort (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

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

void sort(int left,int right)
{
  int i,j;
  int element; // auxulary variable for middle element in the interval
  int swp; // auxulary variable for swapping

  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);
  if (left<j) sort(left,j); // applying the same procedure to the left part
  if (i<right) sort(i,right); // applying the same procedure to the right part
}

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

  // Sorting
  sort(0,N-1);

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



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

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Exchanges (array sorting)
//  (б) Johna Smith, 1996
//
//  Method description:
//    Using linear search find smallest i with the following property:
//    a(i)>a(i+1), swap a(i) and a(i+1) and repeat this process searching
//    from a(i+1). After one pass greatest number will be placed at right
//    place. We'll decrease amount of acting elements on each pass.
//    When 2 elements will left we'll say that array is sorted
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

const N=     10;  // number of array elements

int array[N]={100,15,234,11,63,78,4,200,0,4};
int swp; // auxulary variable for swapping
int m; // number of active elements on current pass

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

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

   // Sorting
   m=N; // All elements acting at first pass
   while (m>1) // repeat passes while there is more than one element
   {
     for (int i=0;i<m-1;i++) // main loop
     {
       // searching for smallest element
       if (array[i]>array[i+1]) // if a(i)>a(i+1)
       {
         // swapping a(i) and a(i+1)
         swp=array[i];
         array[i]=array[i+1];
         array[i+1]=swp;
       }
     }
     m--; // decreasing amount of acting elements
   }

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


Это сообщение отредактировал(а) maxim1000 - 17.12.2006, 03:43
PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 09:20 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



Сортировка файлов


Straight merge

Код

//////////////////////////////////////////////////////////////////////////////
//
//  File sort (straight merge)
//  (c) Johna Smith, 1996
//
//  Method description:
//   This program reads data from file 1.dat and writes sorted data to 1.srt
//   Format of 1.dat: i1 i2 i3 i4 i5 i6 ..., where i(k) is integer
//   It also creates four temporary files: 1.tmp, 2.tmp, 3.tmp, 4.tmp
//   We use four files to sort given sequence of data.
//     1) Split file into 2 parts (A and B)
//     2) Combine A and B to C, sorting 2-element groups
//     3) Repeat these steps for file C, sorting 4-element groups,
//        8-elements,... while size of group less than file size.
//     We can optimize this method by uniting splitting and combining:
//     we just need to change destination on each step: combine one
//     group to C another to D, next group again to C and so on
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

int copy(FILE *A, FILE *B, FILE *C, int n)
{
 int a,b,q,r,m=0;

 q=r=0;
 // combining to C
 if (!feof(A)) fread(&a,sizeof(int),1,A);
 if (!feof(B)) fread(&b,sizeof(int),1,B);
 while (!feof(A) && !feof(B) && q<n && r<n)
 {
   if (a<b)
   {
     fwrite(&a,sizeof(int),1,C);m++;q++;
     if (q<n) fread(&a,sizeof(int),1,A);
   }
   else
   {
     fwrite(&b,sizeof(int),1,C);m++;r++;
     if (r<n) fread(&b,sizeof(int),1,B);
   }
 }
 // writing remainders
 while (!feof(A) && q<n)
 {
   fwrite(&a,sizeof(int),1,C);m++;q++;
   if (q<n) fread(&a,sizeof(int),1,A);
 }
 while (!feof(B) && r<n)
 {
   fwrite(&b,sizeof(int),1,C);m++;r++;
   if (r<n) fread(&b,sizeof(int),1,B);
 }

 return m;
}

void main(void)
{
 FILE *A,*B,*C,*D,*input,*output;
 int N=0,c,m,p,l;

 // opening needed files (file1.dat must exist)
 input=fopen("1.dat","rb");
 output=fopen("1.srt","wb");
 l=0;
 A=fopen("1.tmp","w+b");
 B=fopen("2.tmp","w+b");
 C=fopen("3.tmp","w+b");
 D=fopen("4.tmp","w+b");


 // reading from source and splitting data into A & B
 while (!feof(input))
 {
   fscanf(input,"%d ",&c);
   fwrite(&c,sizeof(int),1,(N%2==0?A:B));
   N++;
 }
 rewind(A);
 rewind(B);

 p=1;
 while (p<N)
 {
   m=N;
   while (m!=0)
   {
     // uniting splitting and combining
     // here we combine from A and B to C and D
     m-=copy(A,B,C,p);
     m-=copy(A,B,D,p);
   }
   // swap A,B and C,D
   // now we'll combine from C and D and split to A and B
   fclose(A);
   fclose(B);
   fclose(C);
   fclose(D);
   A=fopen((l?"1.tmp":"3.tmp"),"rb");
   B=fopen((l?"2.tmp":"4.tmp"),"rb");
   C=fopen((!l?"1.tmp":"3.tmp"),"wb");
   D=fopen((!l?"2.tmp":"4.tmp"),"wb");
   l=!l;
   // groups will be larger twice
   p*=2;
 }

 fread(&c,sizeof(int),1,A);
 while (!feof(A))
 {
   fprintf(output,"%d ",c);
   fread(&c,sizeof(int),1,A);
 }

 fclose(A);
 fclose(B);
 fclose(C);
 fclose(D);
 fclose(input);
 fclose(output);
}



Natural merge

Код

//////////////////////////////////////////////////////////////////////////////
//
//  File sort (natural merge)
//  (c) Johna Smith, 1996
//
//  Method description:
//   This program reads data from file 1.dat and writes sorted data to 1.srt
//   Format of 1.dat: i1 i2 i3 i4 i5 i6 ..., where i(k) is integer
//   It also creates four temporary files: 1.tmp, 2.tmp, 3.tmp
//   We use three files to sort given sequence of data.
//
//   Series - is a sequence of numbers with the following property: a(i)<=a(i+1)
//   1) Divide copy of the given file C into parts A and B, writing first
//      series from C to A, second one to B, third - to C, and so on
//   2) Combine files A and B into C, sorting distributed series
//   3) Repeat first two steps until there is more than one series
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

void main(void)
{
 FILE *A,*B,*C,*input,*output;
 int sorted,a,b,c,old;

 // opening needed files (file1.dat must exist)
 input=fopen("1.dat","rb");
 output=fopen("1.srt","wb");
 A=fopen("1.tmp","w+b");
 B=fopen("2.tmp","w+b");
 C=fopen("3.tmp","w+b");

 // reading from source and copying data into C
 while (!feof(input))
 {
   fscanf(input,"%d ",&c);
   fwrite(&c,sizeof(int),1,C);
 }
 rewind(C);

 do
 {
   // reopeneing files A and B
   fclose(A);
   fclose(B);
   A=fopen("1.tmp","w+b");
   B=fopen("2.tmp","w+b");
   rewind(C);

   // distributing series from C to A and B
   fread(&c,sizeof(int),1,C);
   sorted=1;
   while (!feof(C))
   {
     // series to A
     do
     {
        old=c;
        fwrite(&c,sizeof(int),1,A);
        fread(&c,sizeof(int),1,C);
     } while (!feof(C) && old<=c);
     old=c;
     // series to B
     while (!feof(C) && old<=c)
     {
       old=c;
       fwrite(&c,sizeof(int),1,B);
       fread(&c,sizeof(int),1,C);
       sorted=0; // if there are more than 1 series then sequence isn't sorted
     }
   }
   // reopening file C
   fclose(C);
   C=fopen("3.tmp","w+b");
   rewind(A);
   rewind(B);
   // combining A and B back to C
   fread(&a,sizeof(int),1,A);
   fread(&b,sizeof(int),1,B);
   while (!feof(A) && !feof(B))
   {
     if (a<b)
     {
       old=a;
       fwrite(&a,sizeof(int),1,C);
       fread(&a,sizeof(int),1,A);
       if (feof(A) || a<old)
       do
       {
         old=b;
         fwrite(&b,sizeof(int),1,C);
         fread(&b,sizeof(int),1,B);
       } while(!feof(B) && b>=old);
     }
     else
     {
       old=b;
       fwrite(&b,sizeof(int),1,C);
       fread(&b,sizeof(int),1,B);
       if (feof(B) || b<old)
       do
       {
         old=a;
         fwrite(&a,sizeof(int),1,C);
         fread(&a,sizeof(int),1,A);
       } while(!feof(A) && a>=old);
     }
   }

   // copying remainder from A (or B) to C
   while (!feof(A))
   {
     fwrite(&a,sizeof(int),1,C);
     fread(&a,sizeof(int),1,A);
   }
   while (!feof(B))
   {
     fwrite(&b,sizeof(int),1,C);
     fread(&b,sizeof(int),1,B);
   }
 } while (!sorted);  // sort until there will be only one series

 // writing sorted sequence
 rewind(C);
 fread(&c,sizeof(int),1,C);
 while (!feof(C))
 {
   fprintf(output,"%d ",c);
   fread(&c,sizeof(int),1,C);
 }

 // closing all opened files
 fclose(A);
 fclose(B);
 fclose(C);
 fclose(input);
 fclose(output);
}



Balanced merge

Код

//////////////////////////////////////////////////////////////////////////////
//
//  File sort (balanced merge)
//  (c) Johna Smith, 1996
//
//  Method description:
//   This program reads data from file 1.dat and writes sorted data to 1.srt
//   Format of 1.dat: i1 i2 i3 i4 i5 i6 ..., where i(k) is integer
//   It also creates four temporary files: 1.tmp, 2.tmp, 3.tmp
//   We use N files to sort given sequence of data.
//
//   Series - is a sequence of numbers with the following property: a(i)<=a(i+1)
//   1) Divide copy of the given file C into N files, writing first
//      series from C to F1, second one to F2, third - to F3, and so on
//   2) Combine files F1,F2,F3,...,Fn, sorting distributed series
//   3) Repeat first two steps until there is more than one series
//
//////////////////////////////////////////////////////////////////////////////

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

#define N 10  // number of temporary files (must be even and greater than 2)
#define Nh N/2 // must be an integer

struct ExtFile
{
 FILE *F;
 int first; // current element of this file
 char eor; // indicates end of series
};

// this function copies data element from x to y, updating (.first) members
void copy(ExtFile *x, ExtFile *y)
{
 y->first=x->first;
 fwrite(&y->first,sizeof(int),1,y->F);
 fread(&x->first,sizeof(int),1,x->F);
 x->eor=x->first<y->first;
}

void main(void)
{
 FILE *input,*output;
 ExtFile F[N];
 int l,old,j,mx,tx,min,x,t[N],ta[N],k1,k2;
 char name[13];

 // opening needed files (file1.dat must exist)
 input=fopen("1.dat","rb");
 output=fopen("1.srt","wb");
 for (int i=0;i<N;i++)
 {
   gcvt(i,8,name);
   strcat(name,".tmp");
   F[i].F=fopen(name,"w+b");
   F[i].eor=0;
 }

 // distributing given data
 j=Nh;
 l=0;
 fscanf(input,"%d ",&x);
 while (!feof(input))
 {
   if (j<Nh-1) j++; else j=0;
   do
   {
     old=x;
     fwrite(&x,sizeof(int),1,F[j].F);
     fscanf(input,"%d ",&x);
   } while (!feof(input) && x>=old);
   l++;
 }
 if (x<old)
 {
   if (j<Nh-1) j++; else j=0;
   l++;
 }
 fwrite(&x,sizeof(int),1,F[j].F);

 for (i=0;i<N;i++) t[i]=i;

 do
 {
   // merging t[0]..t[Nh-1] to t[Nh]..t[N-1]
   if(l<Nh) k1=l-1; else k1=Nh-1;
   for(i=0;i<=k1;i++)
   {
     rewind(F[t[i]].F);
     ta[i]=t[i];
     fread(&F[t[i]].first,sizeof(int),1,F[t[i]].F);
   }
   l=0; // number of input series
   j=Nh; // index of input sequence
   // merging input series to t[j]
   do
   {
     l++;
     k2=k1;
     // selecting minimal element from F1..Fn
     do
     {
       i=1;
       mx=0;
       min=F[ta[0]].first;
       while (i<=k2)
       {
         x=F[ta[i]].first;
         if(x<min)
         {
           min=x;
           mx=i;
         }
         i++;
       }
       copy(&F[ta[mx]],&F[t[j]]);
       if (feof(F[ta[mx]].F))
       {
         // excluding file
         fclose(F[ta[mx]].F);
         gcvt(ta[mx],8,name);
         strcat(name,".tmp");
         F[ta[mx]].F=fopen(name,"w+b");
         F[ta[mx]].eor=0;
         ta[mx]=ta[k2];
         ta[k2]=ta[k1];
         k1--;
         k2--;
       } else
       if (F[ta[mx]].eor)
       {
         // closing series
         tx=ta[mx];
         ta[mx]=ta[k2];
         ta[k2]=tx;
         k2--;
       }
     } while(k2>=0);
     if(j<N-1) j++; else j=Nh;
   } while(k1>=0);
   for(i=0;i<Nh;i++)
   {
     tx=t[i];
     t[i]=t[i+Nh];
     t[i+Nh]=tx;
   }
 } while (l!=1);

 // writing sorted sequence
 rewind(F[0].F);
 fread(&x,sizeof(int),1,F[0].F);
 while (!feof(F[0].F))
 {
   fprintf(output,"%d ",x);
   fread(&x,sizeof(int),1,F[0].F);
 }

 // closing all opened files
 for (i=0;i<N;i++) fclose(F[i].F);
 fclose(input);
 fclose(output);
}



Polyphase merge

Код

//////////////////////////////////////////////////////////////////////////////
//
//  File sort (polyphase merge)
//  (c) Johna Smith, 1996
//
//  Method description:
//   This program reads data from file 1.dat and writes sorted data to 1.srt
//   Format of 1.dat: i1 i2 i3 i4 i5 i6 ..., where i(k) is integer
//   It also creates four temporary files: 1.tmp, 2.tmp, 3.tmp
//   We use N files to sort given sequence of data.
//
//   Series - is a sequence of numbers with the following property: a(i)<=a(i+1)
//   1) Divide copy of the given file C into N-1 files, writing first
//      series from C to F1, second one to F2, third - to F3, and so on
//   2) Combine series from files F1,F2,F3,...,F(n-1) to Fn,
//      sorting distributed series
//   3) When we reach the end of one of F1..F(n-1) files we turn sequences
//      and will combine series to this new empty file. For example if we reach
//      end of F3 we will gather data from F1,F2,F4,..Fn to F3
//   4) Repeat first three steps until there is more than one series
//
//////////////////////////////////////////////////////////////////////////////

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

#define N 6  // number of temporary files (must be greater than 2)

struct ExtFile
{
 FILE *F;
 int first; // current element of this file
 char eor; // indicates end of series
};

 FILE *output;
 ExtFile input;
 ExtFile F[N];  // array of files
 int level,i,j,k,mx,tx,xmin,x,t[N],ta[N],a[N],d[N],z,tn,dn;
 char name[13];

// this functions sets new value of variable j
// on the next step of data distribution we'll write to sequence #j
void select(void)
{
 if (d[j]<d[j+1]) j++;  // d[j] shows number of series in sequence #j
 else
 {
   if (d[j]==0)
   {
     level++;
     z=a[0];
     for (int i=0;i<N-1;i++)
     {
       d[i]=z+a[i+1]-a[i];
       a[i]=z+a[i+1];
     }
   }
   j=0;
 }
 d[j]--;
}

// this function copies data element from x to y, updating (.first) members
void copy(ExtFile *x, ExtFile *y)
{
 y->first=x->first;
 fwrite(&y->first,sizeof(int),1,y->F);
 fread(&x->first,sizeof(int),1,x->F);
 x->eor=(x->first<y->first || feof(x->F));
}

// this function copies series from input to f[j]
void copyrun(void)
{
 do
 {
   F[j].first=input.first;
   fwrite(&F[j].first,sizeof(int),1,F[j].F);
   fscanf(input.F,"%d",&input.first);
   input.eor=input.first<F[j].first;
 } while (!input.eor && !feof(input.F));
 if (feof(input.F))
 {
   if (input.first<F[j].first) select();
   F[j].first=input.first;
   fwrite(&F[j].first,sizeof(int),1,F[j].F);
   d[j]++;
 }
}

void main(void)
{
 // opening required files (file 1.dat must exist)
 input.F=fopen("1.dat","rb");
 fscanf(input.F,"%d ",&input.first);
 F[t[N-1]].eor=0;
 output=fopen("1.srt","wb");
 for (int i=0;i<N-1;i++)
 {
   gcvt(i,8,name);
   strcat(name,".tmp");
   F[i].F=fopen(name,"w+b");
   fread(&F[i].first,sizeof(int),1,F[i].F);
   F[i].eor=0;
   a[i]=1;
   d[i]=1;
 }

 // initializing variables
 level=1;
 j=0;
 a[N-1]=0;
 d[N-1]=0;

 // distributing data from input to first N-1 sequences
 do
 {
   select();
   copyrun();
 } while (!feof(input.F) && j!=N-2);

 while (!feof(input.F))
 {
   select();
   if (F[j].first<=input.first)
   {
     copyrun();
     if (feof(input.F)) d[j]++; else copyrun();
   } else copyrun();
 }

 // preparing first N-1 series for reading
 for (i=0;i<N-1;i++)
 {
   t[i]=i;
   rewind(F[i].F);
   fread(&F[i].first,sizeof(int),1,F[i].F);
   F[i].eor=0;
 }
 t[N-1]=N-1;

 // main loop
 do
 {
   // gathering data from t[0]..t[N-2] into t[N-1]
   z=a[N-2];
   d[N-1]=0;
   // preparing file to write
   fclose(F[t[N-1]].F);
   gcvt(t[N-1],8,name);
   strcat(name,".tmp");
   F[t[N-1]].F=fopen(name,"w+b");
   F[t[N-1]].eor=0;
   do
   {
     // gathering one series
     k=-1;
     for (i=0;i<N-1;i++)
     {
       if (d[i]>0) d[i]--;
       else
       {
         k++;
         ta[k]=t[i];
       }
     }
     if (k==-1) d[N-1]++;
     else
     {
       // gathering one real series from t[0]..t[k] into t[N-1]
       do
       {
         i=0;
         mx=0;
         xmin=F[ta[i]].first;
         // selecting series with minimal element to place it first
         while (i<k)
         {
           i++;
           x=F[ta[i]].first;
           if (x<xmin)
           {
             xmin=x;
             mx=i;
           }
         }
         // copying smallest element to output series
         copy(&F[ta[mx]],&F[t[N-1]]);
         // if we reach the end of the current input sequence
         // or end of series in it then close this source
         if (F[ta[mx]].eor)
         {
           // closing this source
           ta[mx]=ta[k];
           k--;
         }
       } while (k!=-1);
     }  
     z--;
   } while (z!=0);

   // swapping sequences
   rewind(F[t[N-1]].F);
   fread(&F[t[N-1]].first,sizeof(int),1,F[t[N-1]].F);
   F[t[N-1]].eor=0;

   tn=t[N-1];
   dn=d[N-1];
   z=a[N-2];
   for(i=N-1;i>0;i--)
   {
     t[i]=t[i-1];
     d[i]=d[i-1];
     a[i]=a[i-1]-z;
   }
   t[0]=tn;
   d[0]=dn;
   a[0]=z;

   // preparing file F[t[N-1]] to write
   fclose(F[t[N-1]].F);
   gcvt(t[N-1],8,name);
   strcat(name,".tmp");
   F[t[N-1]].F=fopen(name,"w+b");
   F[t[N-1]].eor=0;

   level--;
 } while (level!=0); // while there are more than one series

 // writing sorted sequence from F[t[0]]
 rewind(F[t[0]].F);
 fread(&x,sizeof(int),1,F[t[0]].F);
 while (!feof(F[t[0]].F))
 {
   fprintf(output,"%d ",x);
   fread(&x,sizeof(int),1,F[t[0]].F);
 }

 // closing all opened files
 for (i=0;i<N;i++) fclose(F[i].F);
 fclose(input.F);
 fclose(output);
}

PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 09:36 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



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


Binary search in sorted string

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Binary search in sorted string
//  (c) Johna Smith, 1996
//
//  Method description:
//     1) Split string into two parts
//     2) If the middle symbol of the string is greater than searched one -
//        select left part else select right part
//     3) Split selected part into two parts (see 1)
//   When in the selected part will be only one symbol compare it with
//   searched one.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

char string[]="0245789AAABDFHJKLXYZ"; // sorted ASCIIZ string
char symbol='J'; // symbol to find
unsigned int left,right; // bounds of the area where we search
unsigned int middle; // center of the area where we search
unsigned int length=sizeof(string)/sizeof(char); // string length

void main(void)
{
 printf("String: %s\n",string);

 left=0;
 right=length-1;

 // Here we avoid checking condition 'string[middle]==symbol'
 // because maximum number of steps is log N (where N is length
 // of the string) and some auxulary steps will be more effective
 // than comparing elements each step.
 while (left<right)
 {
   middle=(left+right)/2;
   if (string[middle]<symbol) left=middle+1; else right=middle;
 }

 // We'll stop when left==right => if there is searched symbol in
 // the string both bounds will point at it
 if (string[right]==symbol) printf("Symbol '%c' found at position %d.",symbol,right+1);
 else printf("Symbol '%c' not found.",symbol);
}



Substring search. Boyer-Moore algorithm

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Substring search. Boyer-Moore algorithm
//  (c) Johna Smith, 1996
//
//  Method description:
//    1) Build auxulary array D. D contains 256 elements (it's
//    amount of symbols in ASCII table. D(i) equals to length
//    of the substring if it doesn't contain character with code i.
//    And if substring contains character with code i then D(i)
//    equals to the distance between end of string and position of
//    this character closest to the end of string.
//    2) Scan string from the start, comparing it with the substring
//    If diferrence found we can skip next d(string(i)) positions,
//    where string(i) is the code of i-th character of the string and
//    i is the current position in the string.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

char string[]="Unsorted string. ABCBBCBCA."; // ASCIIZ string
char substring[]="string"; // substring to find
const unsigned int string_length=sizeof(string)/sizeof(char)-1; // string length
const unsigned int substring_length=sizeof(substring)/sizeof(char)-1; // substring length
int i,j,k;
int d[256]; // auxulary array

void main(void)
{
 printf("String: %s\n",string);


 // Analysing substring, building array d
 for (int i=0;i<256;i++) d[i]=substring_length;
 for (i=0;i<substring_length-1;i++) d[substring[i]]=substring_length-i-1;

 // Searching
 i=substring_length;
 do
 {
   j=substring_length; k=i;
   do // search difference between string and substring
   {
     k--;
     j--;
   } while (j>=0 && string[k]==substring[j]);
   i+=d[string[i-1]]; // skip next d(string(i)) positions
 } while (j>=0 && i<string_length);

 // condition j<0 means that substring is found
 if (j<0) printf("Substring \"%s\" found.",substring);
 else printf("Substring \"%s\" not found.",substring);
}



Linear substring search

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Linear substring search
//  (c) Johna Smith, 1996
//
//  Method description:
//    Simple search method.
//    1) Compare string and substring from first symbol of the string
//    2) If difference found compare string and substring from
//       second symbol of the string
//    Repeat these steps until no differencies found or end of string reached
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

char string[]="Unsorted string"; // ASCIIZ string
char substring[]="ring"; // substring to find
unsigned int position; // comparing from this position of the string
unsigned int subposition; // compared position of the substring
unsigned int string_length=sizeof(string)/sizeof(char); // string length
unsigned int substring_length=sizeof(substring)/sizeof(char); // substring length

void main(void)
{
 printf("String: %s\n",string);

 position=0;
 // finish compare when all symbols of substring compared (substring found)
 // or there left less than substring_length uncompared symbols
 // in the string (substring not found)
 while (subposition!=substring_length && position!=string_length-substring_length+1)
 {
   subposition=0; // comparing from the beginning of the substring
   while (subposition<substring_length &&
             string[position+subposition]==substring[subposition])
     subposition++;
   // compare from next position of the string
   position++;
 }

 if (subposition==substring_length) printf("Substring \"%s\" found from position %d.",substring,position);
 else printf("Substring \"%s\" not found.",substring);
}



Substring search. Knuth-Morris-Pratt algorithm

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Substring search. Knuth-Morris-Pratt algorithm
//  (c) Johna Smith, 1996
//
//  Method description:
//    1) Build auxulary array D. D contains "substring_length" elements.
//    D(i) equals to the length of longest sequence of characters
//    of the substring with such properties:
//         ss(0)..ss(d(i)-1)==ss(i-d(i))..p(i-1);  (ss is a substring)
//         p(d(i))!=p(i);
//    2) Scan string from the start, comparing it with the substring
//    If diferrence found we can skip next d(string(i)) positions,
//    where string(i) is the code of i-th character of the string and
//    i is the current position in the string.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

char string[]="Unsorted string. ABCBBCBCA."; // ASCIIZ string
char substring[]="BCB"; // substring to find
const unsigned int string_length=sizeof(string)/sizeof(char)-1; // string length
const unsigned int substring_length=sizeof(substring)/sizeof(char)-1; // substring length
int i,j,k;
int d[substring_length]; // auxulary array

void main(void)
{
 printf("String: %s\n",string);

 // Analysing substring, building array d
 j=0; k=-1; d[0]=-1;
 while (j<substring_length-1)
 {
   while (k>=0 && substring[j]!=substring[k]) k=d[k];
   j++; k++;
   if (substring[j]==substring[k]) d[j]=d[k]; else d[j]=k;
 }

 // Searching
 i=j=k=0;
 while (j<substring_length && i<string_length)
 {
   // here j is the current position in the substring and i is
   // the current position in the string
   while (j>=0 && string[i]!=substring[j]) j=d[j]; //skip next d[j] positions
   i++; j++;
 }
 // Condition j==substring length means that substring is found
 if (j==substring_length) printf("Substring \"%s\" found.",substring);
 else printf("Substring \"%s\" not found.",substring);
}



Linear search in the string

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Linear search
//  (c) Johna Smith, 1996
//
//  Method description:
//    Simple search method. Take first element - if it is what we
//   need - stop searh else take next element.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

char string[]="Unsorted string"; // ASCIIZ string
char symbol='d'; // symbol to find
char *p;
unsigned int position; // position of the found symbol
unsigned int length=sizeof(string)/sizeof(char); // string length

void main(void)
{
 printf("String: %s\n",string);

 // changing last (terminating) symbol of the string
 // to the symbol that we search allows us to avoid
 // checking condition 'end of string' - we always
 // will find symbol in the modified string
 string[length-1]=symbol;
 for (p=string,position=1;*p!=symbol;p++,position++);
 string[length-1]='\0'; //restoring last symbol
 
 if (position!=length) printf("Symbol '%c' found at position %d.",symbol,position);
 else printf("Symbol '%c' not found.",symbol);
}

PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 09:57 (ссылка)    | (голосов:1) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



Матрицы


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

Код

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



Обращение

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Matrix operations (inversion)
//  (c) Johna Smith, 1996
//
//  Given: A (N x N), det A != 0                 -1
//  This algorithm inverts matrix A and returns A  .
//      -1
//   A*A  = E
//  We simply convert matrix A into matrix E, and do the same operations
//  over matrix B (initially B=E). Finally A=E, B=A^(-1).
//
//////////////////////////////////////////////////////////////////////////////

#define N               3

#include <stdio.h>

float a[N][N]={{4,8,0},{8,8,8},{2,0,1}};

void Invert(float *matrix)
{
 float e[N][N];

 // initializing matrix e
 for (int i=0;i<N;i++)
   for (int j=0;j<N;j++)
     e[i][j]=(i==j?1:0);
 // converting matrix to e
 for(i=0;i<N;i++)
 {
   // normalizing row (making first element =1)
   float tmp=*(matrix+i*N+i);
   for(int j=N-1;j>=0;j--)
   {
     e[i][j]/=tmp;
     *(matrix+i*N+j)/=tmp;
   }
   // excluding i-th element from each row except i-th one
   for(j=0;j<N;j++)
     if (j!=i)
     {
       tmp=*(matrix+j*N+i);
       for(int k=N-1;k>=0;k--)
       {
         e[j][k]-=e[i][k]*tmp;
         *(matrix+j*N+k)-=*(matrix+i*N+k)*tmp;
       }
     }
 }

 // now e contains inverted matrix so we need only to copy e to matrix
 for(i=0;i<N;i++)
   for(int j=0;j<N;j++)
     *(matrix+i*N+j)=e[i][j];
}

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

void main(void)
{
 // display matrix A
 show_matrix(a[0]);
 // Invert it
 Invert(a[0]);
 // display the inverted matrix
 show_matrix(a[0]);
 // Invert it again
 Invert(a[0]);
 // display the inversion of the inverted matrix
 show_matrix(a[0]);
}



Сложение

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Matrix operations (addition)
//  (c) Johna Smith, 1996
//
//  Given: A,B - matrixes M x N
//  This algorithm sum these matrixes and display the result.
//
//////////////////////////////////////////////////////////////////////////////

#define M    5
#define N    6

#include <stdio.h>

int a[N][M]={{1,2,3,4,5},{5,4,3,2,1},{1,2,3,4,5},{5,4,3,2,1},{1,2,3,4,5},{0,1,0,1,0}};
int b[N][M]={{5,4,3,2,1},{1,2,3,4,5},{5,4,3,2,1},{1,2,3,4,5},{5,4,3,2,1},{1,0,1,0,1}};
int c[N][M];

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

void main(void)
{
 // display matrixes A and B
 show_matrix(a[0]);
 show_matrix(b[0]);

 // sum these matrixes
 for(int i=0;i<N;i++)
   for(int j=0;j<M;j++)
     c[i][j]=a[i][j]+b[i][j];

 // display the sum
 show_matrix(c[0]);
}



Умножение

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Matrix operations (multiplication)
//  (c) Johna Smith, 1996
//
//  Given: A (M x N) and B (N x P)
//  This algorithm multiplies these matrixes and display the result.
//
//////////////////////////////////////////////////////////////////////////////

#define M     3
#define N     4
#define P     2

#include <stdio.h>

int a[M][N]={{2,0,3,1},{5,1,2,0},{0,0,4,1}};
int b[N][P]={{1,3},{2,1},{4,0},{3,5}};
int c[M][P];

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

void main(void)
{
 // display matrixes A and B
 show_matrix(a[0],M,N);
 show_matrix(b[0],N,P);

 // substract these matrixes
 for(int i=0;i<M;i++)
   for(int j=0;j<P;j++)
   {
     c[i][j]=0;
     for (int k=0;k<N;k++) c[i][j]+=a[i][k]*b[k][j];
   }
 // display the difference
 show_matrix(c[0],M,P);
}



Вычитание

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Matrix operations (substraction)
//  (c) Johna Smith, 1996
//
//  Given: A,B - matrixes M x N
//  This algorithm substract these matrixes and display the result.
//
//////////////////////////////////////////////////////////////////////////////

#define M   5
#define N   6

#include <stdio.h>

int a[N][M]={{4,4,4,4,4},{5,4,3,2,1},{1,2,3,4,5},{5,4,3,2,1},{1,2,3,4,5},{0,1,0,1,0}};
int b[N][M]={{3,4,4,4,4},{5,3,3,2,1},{1,2,2,4,5},{5,4,3,1,1},{1,2,3,4,4},{0,1,0,1,0}};
int c[N][M];

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

void main(void)
{
 // display matrixes A and B
 show_matrix(a[0]);
 show_matrix(b[0]);

 // substract these matrixes
 for(int i=0;i<N;i++)
   for(int j=0;j<M;j++)
     c[i][j]=a[i][j]-b[i][j];

 // display the difference
 show_matrix(c[0]);
}



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

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Matrix operations (transponating)
//  (c) Johna Smith, 1996
//
//  Given: A (N x N),                                 T
//  This algorithm transponates matrix A and returns A  .
//          T
//   A(ji)=A (ij)
//
//////////////////////////////////////////////////////////////////////////////

#define N  3

#include <stdio.h>

float a[N][N]={{4,8,0},{8,8,8},{2,0,1}};

void Transponate(float *matrix)
{
 float swp;

 for (int i=0;i<N;i++)
   for (int j=i+1;j<N;j++)
   {
     swp=*(matrix+i*N+j);
     *(matrix+i*N+j)=*(matrix+j*N+i);
     *(matrix+j*N+i)=swp;
   }
}

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

void main(void)
{
 // display matrix A
 show_matrix(a[0]);
 // Transponate it
 Transponate(a[0]);
 // display transponated matrix
 show_matrix(a[0]);
 // Transponate it again
 Transponate(a[0]);
 // display the transponation of the transponated matrix
 show_matrix(a[0]);
}

PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 10:08 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



Операции с комплексными величинами


Сложение и вычитание в алгебраической форме

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (addition and substraction in a+bi form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=a1+i*b1, z2=a2+i*b2
//  z1+z2=(a1+a2)+i*(b1+b2)
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

struct complex
{
 float re;
 float im;
};

void show_complex(complex c) // this functions displays complex value
{
 printf("%f%+fi",c.re,c.im);
}

complex Add(complex a, complex b)
{
 complex c;

 c.re=a.re+b.re;
 c.im=a.im+b.im;

 return c;
}

complex Substract(complex a, complex b)
{
 complex c;

 c.re=a.re-b.re;
 c.im=a.im-b.im;

 return c;
}

complex a={2,3},b={4,6};

void main(void)
{
 // addition
 show_complex(a);
 printf(" + ");
 show_complex(b);
 printf(" = ");
 show_complex(Add(a,b));

 // substraction
 printf("\n");
 show_complex(b);
 printf(" - ");
 show_complex(a);
 printf(" = ");
 show_complex(Substract(b,a));
}



Деление в алгебраической форме

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (division in a+bi form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=a1+i*b1, z2=a2+i*b2
//  z1/z2=(a1*a2-b1*b2)/(a2^2+b2^2)+i*(b1*a2-a1*b2)/(a2^2+b2^2)
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

struct complex
{
 float re;
 float im;
};

void show_complex(complex c) // this functions displays complex value
{
 printf("%f%+fi",c.re,c.im);
}

complex Div(complex a, complex b)
{
 complex c;

 c.re=(a.re*b.re+a.im*b.im)/(b.re*b.re+b.im*b.im);
 c.im=(a.im*b.re-a.re*b.im)/(b.re*b.re+b.im*b.im);

 return c;
}

complex a={2,3},b={-1,2};

void main(void)
{
 show_complex(a);
 printf(" / ");
 show_complex(b);
 printf(" = ");
 show_complex(Div(a,b));
}



Умножение в алгебраической форме

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (multiplication in a+bi form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=a1+i*b1, z2=a2+i*b2
//  z1*z2=(a1*a2-b1*b2)+i*(a1*b2+a2*b1)
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

struct complex
{
 float re;
 float im;
};

void show_complex(complex c) // this functions displays complex value
{
 printf("%f%+fi",c.re,c.im);
}

complex Mul(complex a, complex b)
{
 complex c;

 c.re=a.re*b.re-a.im*b.im;
 c.im=a.im*b.re+a.re*b.im;

 return c;
}

complex a={2,3},b={-1,2};

void main(void)
{
 show_complex(a);
 printf(" * ");
 show_complex(b);
 printf(" = ");
 show_complex(Mul(a,b));
}



Деление в экспоненциальной форме

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (division in M*exp(i*phi) form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=M1*exp(i*phi1), z2=M2*exp(i*phi2)
//  z1/z2=M1/M2*exp(i*(phi1-phi2))
//
//////////////////////////////////////////////////////////////////////////////

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

#define Pi      3.1415926536

struct exp_complex
{
 float M;
 float phi;
};

void show_exp_complex(exp_complex c) // this functions displays complex value
{
 printf("%f*exp(%f*i)",c.M,c.phi);
}

exp_complex Div(exp_complex a,exp_complex b)
{
 exp_complex c;

 c.M=a.M/b.M;
 c.phi=a.phi-b.phi;

 return c;
}

exp_complex a={1,0}, b={-1,-Pi/2};

void main(void)
{
 show_exp_complex(a);
 printf(" / ");
 show_exp_complex(b);
 printf(" = ");
 show_exp_complex(Div(a,b));
}



Умножение в экспоненциальной форме

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Complex values operations (multiplication in M*exp(i*phi) form)
//  (c) Johna Smith, 1996
//
//  Given: z1,z2 - complex values
//  z1=M1*exp(i*phi1), z2=M2*exp(i*phi2)
//  z1*z2=M1*M2*exp(i*(phi1+phi2))
//
//////////////////////////////////////////////////////////////////////////////

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

#define Pi      3.1415926536

struct exp_complex
{
 float M;
 float phi;
};

void show_exp_complex(exp_complex c) // this functions displays complex value
{
 printf("%f*exp(%f*i)",c.M,c.phi);
}

exp_complex Mul(exp_complex a,exp_complex b)
{
 exp_complex c;

 c.M=a.M*b.M;
 c.phi=a.phi+b.phi;

 return c;
}

exp_complex a={-1,Pi/2}, b={-1,-Pi/2};

void main(void)
{
 show_exp_complex(a);
 printf(" * ");
 show_exp_complex(b);
 printf(" = ");
 show_exp_complex(Mul(a,b));
}

PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 10:28 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



Математическая статистика


Интервальная оценка среднего (рассмотрены случаи значений доверительной вероятности 0.95 и 0.99)

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating <x>, (x-dx,x+dx)
//  (c) Johna Smith, 1996
//
//  Method description:
//    dx=t(alpha,n)*Sn/sqrt(n)
//    <x>=Sum x/n
//
//////////////////////////////////////////////////////////////////////////////

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

#define N   6 // number of experiments
float x[N]={3,5,4,5,3,4}; // experimental data

double Student095(unsigned int n)
{
 double result;

 result=1.96+2.4/(n-1)+5.901610*pow(n-1,-2.372);

 return result;
}

double Student099(unsigned int n)
{
 double result;

 result=2.576+5.0/(n-1)+29.12178*pow(n-1,-2.591843);

 return result;
}

// this function calculates average x
double Average(void)
{
 double sum=0;

 for (int i=0;i<N;i++) sum+=x[i];

 return sum/N;
}

// this function calculates delta x
double DeltaX(void)
{
 double ave,sum=0,Sn;

 ave=Average();
 for (int i=0;i<N;i++) sum+=x[i]*x[i];
 Sn=sqrt(((double)N/(N-1))*((double)sum/N-ave*ave));

 return Student095(N)*Sn/sqrt(N); // you can use Student099 if you want
}

void main(void)
{
 printf("Data: ");
 for (int i=0;i<N;i++) printf("%f ",x[i]);
 printf("\nx=%f +/- %f",Average(),DeltaX());
}



Вычисление коэффициента Стьюдента

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating Student's coefficients t(alpha,n), if alpha=0.95 or 0.99
//  (c) Johna Smith, 1996
//
//  Method description:
//   To calculate Student's coefficients we use the following
//  approximation:
//    t(n)=C0+C1/(n-1)+C2(n-1)^(-C3)
//  Coefficients Ci for alpha=0.95 and alpha=0.99 is given
//  n - number of experiments
//
//////////////////////////////////////////////////////////////////////////////

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

double Student095(unsigned int n)
{
 double result;

 result=1.96+2.4/(n-1)+5.901610*pow(n-1,-2.372);

 return result;
}

double Student099(unsigned int n)
{
 double result;

 result=2.576+5.0/(n-1)+29.12178*pow(n-1,-2.591843);

 return result;
}

void main(void)
{
 printf("alpha=0.95 n=5 t5= %g\n",Student095(5));
 printf("alpha=0.99 n=15 t15= %g",Student099(15));
}



Теоретический коэффициент корреляции для двух величин

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating theoretical corellation coefficient for two values
//  (c) Johna Smith, 1996
//
//  Method description:
//    Corellation coefficient shows dependencies between two values
//
//////////////////////////////////////////////////////////////////////////////

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

#define N       3 // number of experiments

double x[N]={1,2,3};
double y[N]={2,4,5};

double TheoreticalCorellation(void)
{
 double sxy=0,sx2=0,sy2=0;

 for (int i=0;i<N;i++)
 {
   sxy+=x[i]*y[i];
   sx2+=x[i]*x[i];
   sy2+=y[i]*y[i];
 }

 return sxy/sqrt(sx2*sy2);
}

void main(void)
{
 printf("Data:\n");
 for (int i=0;i<N;i++) printf("%f %f\n",x[i],y[i]);
 printf("\nCorellation coefficient r= %f\n",TheoreticalCorellation());
}



Сглаживание экспериментальных данных по 3 точкам

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Smoothing experimental data by 3 points
//  (c) Johna Smith, 1996
//
//  Method description:
//    We smooth experimental data to enclose them to real curve.
//    All formulas use smallest sqares method to approximate data.
//    This program uses smoothing using 3-points groups
//
//////////////////////////////////////////////////////////////////////////////

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

#define N 10      // number of experiments

float data[N]={0.9,2.12,2.92,4.15,4.9,6.1,6.92,8.15,9.05,9.8};
float smoothed[N];

void main(void)
{
 printf("Data:\n");
 for (int i=0;i<N;i++) printf("%.4f ",data[i]);

 //smoothing
 smoothed[0]=(5*data[0]+2*data[1]-data[2])/6;
 for (i=1;i<N-1;i++) smoothed[i]=(data[i-1]+data[i]+data[i+1])/3;
 smoothed[N-1]=(5*data[N-1]+2*data[N-2]-data[N-3])/6;

 printf("\nSmoothed data:\n");
 for (i=0;i<N;i++) printf("%.4f ",smoothed[i]);
}



Сглаживание экспериментальных данных по 5 точкам

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Smoothing experimental data by 5 points
//  (c) Johna Smith, 1996
//
//  Method description:
//    We smooth experimental data to enclose them to real curve.
//    All formulas use smallest sqares method to approximate data.
//    This program uses smoothing using 5-points groups
//
//////////////////////////////////////////////////////////////////////////////

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

#define N  10      // number of experiments

float data[N]={0.9,2.12,2.92,4.15,4.9,6.1,6.92,8.15,9.05,9.8};
float smoothed[N];

void main(void)
{
 printf("Data:\n");
 for (int i=0;i<N;i++) printf("%.4f ",data[i]);

 //smoothing
 smoothed[0]=(3*data[0]+2*data[1]+data[2]-data[4])/5;
 smoothed[1]=(4*data[0]+3*data[1]+2*data[2]+data[3])/10;
 for (i=2;i<N-2;i++) smoothed[i]=(data[i-2]+data[i-1]+data[i]+
                                  data[i+1]+data[i+2])/5;
 smoothed[N-2]=(data[N-4]+2*data[N-3]+3*data[N-2]+4*data[N-1])/10;
 smoothed[N-1]=(3*data[N-1]+2*data[N-2]+data[N-3]-data[N-5])/5;

 printf("\nSmoothed data:\n");
 for (i=0;i<N;i++) printf("%.4f ",smoothed[i]);
}



Сглаживание экспериментальных данных по 7 точкам

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Smoothing experimental data by 7 points
//  (c) Johna Smith, 1996
//
//  Method description:
//    We smooth experimental data to enclose them to real curve.
//    All formulas use smallest sqares method to approximate data.
//    This program uses smoothing using 7-points groups
//
//////////////////////////////////////////////////////////////////////////////

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

#define N 10      // number of experiments

float data[N]={0.7,1.08,1.39,1.64,1.76,1.99,2.04,2.22,2.28,2.42};
float smoothed[N];

void main(void)
{
 printf("Data:\n");
 for (int i=0;i<N;i++) printf("%.4f ",data[i]);

 //smoothing
 smoothed[0]=(39*data[0]+8*data[1]-4*(data[2]+data[3]-data[5]) + data[4]-2*data[6])/42;
 smoothed[1]=(8*data[0]+19*data[1]+16*data[2]+6*data[3]-4*data[4] - 7*data[5]+4*data[7])/42;
 smoothed[2]=(-4*data[0]+16*data[1]+19*data[2]+12*data[3]+2*data[4] - 4*data[5]+data[6])/42;

 for (i=3;i<N-3;i++)
   smoothed[i]=(7*data[i]+6*(data[i-1]+data[i+1]) +
         3*(data[i-2]+data[i+2])-2*(data[i-3]+data[i+3]))/21;

 smoothed[N-3]=(data[N-7]-4*data[N-6]+2*data[N-5]+12*data[N-4] +
         9*data[N-3]+16*data[N-2]-4*data[N-1])/42;
 smoothed[N-2]=(4*data[N-7]-7*data[N-6]-4*data[N-5]+6*data[N-4] +
         16*data[N-3]+19*data[N-2]+8*data[N-1])/42;
 smoothed[N-1]=(-2*data[N-7]+4*data[N-6]+data[N-5]-4*data[N-4] -
         4*data[N-3]+8*data[N-2]+39*data[N-1])/42;

 printf("\nSmoothed data:\n");
 for (i=0;i<N;i++) printf("%.4f ",smoothed[i]);
}



Эмпирический коэффициент корреляции для двух величин

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating empirical corellation coefficient for two values
//  (c) Johna Smith, 1996
//
//  Method description:
//    Corellation coefficient shows dependencies between two values
//
//////////////////////////////////////////////////////////////////////////////

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

#define N 3 // number of experiments

double x[N]={1,2,3};
double y[N]={2,4,5};

double EmpiricalCorellation(void)
{
 double sxy=0,sx2=0,sy2=0,sx=0,sy=0;

 for (int i=0;i<N;i++)
 {
   sxy+=x[i]*y[i];
   sx2+=x[i]*x[i];
   sy2+=y[i]*y[i];
   sx+=x[i];
   sy+=y[i];
 }

 return (sxy-sx*sy/N)/(sqrt(sx2-sx*sx/N)*sqrt(sy2-sy*sy/N));
}

void main(void)
{
 printf("Data:\n");
 for (int i=0;i<N;i++) printf("%f %f\n",x[i],y[i]);
 printf("\nCorellation coefficient r= %f\n",EmpiricalCorellation());
}



Эмпирический коэффициент корреляции для трех величин

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating corellation coefficient for three values
//  (c) Johna Smith, 1996
//
//  Method description:
//    Corellation coefficient shows dependencies between three values
//
//////////////////////////////////////////////////////////////////////////////

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

#define N   3 // number of experiments

double x[N]={1,2,3};
double y[N]={2,4,5};
double z[N]={1,7,0};

double EmpiricalCorellation2D(double *i, double *j)
{
 double sij=0,si2=0,sj2=0,si=0,sj=0;

 for (int a=0;a<N;a++)
 {
   sij+=*(i+a)**(j+a);
   si2+=*(i+a)**(i+a);
   sj2+=*(j+a)**(j+a);
   si+=*(i+a);
   sj+=*(j+a);
 }

 return (sij-si*sj/N)/(sqrt(si2-si*si/N)*sqrt(sj2-sj*sj/N));
}

double Corellation3D(void)
{
 double rxz,ryz,rxy;

 // calculating partial corellation coefficients
 rxy=EmpiricalCorellation2D(x,y);
 rxz=EmpiricalCorellation2D(x,z);
 ryz=EmpiricalCorellation2D(y,z);

 return sqrt((rxz*rxz-ryz*ryz-2*rxy*rxz*ryz)/(1-rxy*rxy));
}

void main(void)
{
 printf("Data:\n");
 for (int i=0;i<N;i++) printf("%f %f %f\n",x[i],y[i],z[i]);
 printf("\nCorellation coefficient R= %f\n",Corellation3D());
}



Частный коэффициент корреляции для трех величин

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating partial corellation coefficient for three values
//  (c) Johna Smith, 1996
//
//  Method description:
//   Partial corellation coefficient shows dependencies between two values
//   using set of three values (we exclude influence of one of these values)
//
//////////////////////////////////////////////////////////////////////////////

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

// excluding influence of z, recieving corellation
// coefficients between x,y  y,z and x,z
double PartialCorellation3D(double rxy,double rxz,double ryz)
{
 return (rxy-rxz*ryz)/sqrt((1-rxz*rxz)*(1-ryz*ryz));
}

void main(void)
{
 double rxy=0.998753,rxz=0.69753,ryz=0.7473;

 printf("Data:\n");
 printf("%f %f %f\n",rxy,rxz,ryz);
 printf("\nPartial corellation coefficient r= %.10f\n",
 PartialCorellation3D(rxy,rxz,ryz));
}

PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 10:38 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



Интегрирование функций


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

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Integrating function (trapecias method)
//  (c) Johna Smith, 1996
//
//  Method description:
//   We approximate area under function with trapecias, using
//   the same intervals by the X axe.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

double f(double x)
{
 // !!! you should specify your own function here

 return x;
}

double Integrate(double x1, double x2, double step)
{
 double x=x1;
 double I=0; // integral value

 while (x<x2-step)
 {
   // calculating trapecia's area as f(x)+f(x+step) * halfheight of it
   I += (f(x)+f(x+step))*step/2;
   x += step;
 }

 return I;
}

void main(void)
{
 printf("I=%f",Integrate(0,1,1e-4));
}



Метод Симпсона

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Integrating function (Simpson's method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    We approximate area under function with parabolas, using
//    three points. We repeat calculation with less step if
//    the diffference between I(h) and I(h/2) is greater than eps
//
//////////////////////////////////////////////////////////////////////////////

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

double f(double x)
{
 // you should specify your own function here
 return x*x*sin(x);
}

double Integrate(double x1, double x2, double step, double precision)
{
 double x=x1;
 double I1=0,I2=0; // integral value
 double s;

 do
 {
   I1=I2;
   s=(f(x1)-f(x2))/2;
   x=x1+step;

   while (x<x2)
   {
     s+=2*f(x)+f(x+step);
     x+=2*step;
   }
   I2=2*step*s/3;

   step/=2.0;  // try once more using less step
 } while (fabs(I1-I2)>precision);

 return I2;
}

void main(void)
{
 printf("I=%f",Integrate(0,1,1e-4,1e-6));
}



Метод прямоугольников

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Integrating function (rectanges method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    We approximate area under function with rectangles, using
//    the same intervals by the X axe.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

double f(double x)
{
 // you should specify your own function here
 return x;
}

double Integrate(double x1, double x2, double step)
{
 double x=x1;
 double I=0; // integral value
 
 while (x<x2-step)
 {
   // calculating function at the halfpoint of interval
   I+=step*f(x+step/2);
   x+=step;
 }
 
 return I;
}

void main(void)
{
 printf("I=%f",Integrate(0,1,1e-4));
}



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

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Integrating function (Monte-Karlo method)
//  (c) Johna Smith, 1996
//
//  Method description:
//   We generate random numbers from interval [a;b] where we want to
//   integrate our function and sum function values in this random points.
//   Here we use approximation of the integral:
//   I=<y>*sigma, where <y> is average value of the function f on [a;b]
//   and sigma is volume of integrating area (if we integrate by interval it
//   is its length)
//
//////////////////////////////////////////////////////////////////////////////

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

double f(double x)
{
 // you should specify your own function here
 return x*x*sin(x);
}

double Integrate(double x1, double x2, long int steps)
{
 double I=0; // integral value

 randomize();  // initializing random numbers generator
 for (long int i=0;i<steps;i++)
   I+=f((random(10000)/10000.0)*(x2-x1)+x1);
 I*=(x2-x1)/steps;

 return I;
}

void main(void)
{
 printf("I=%f",Integrate(0,1,1e5));
}



Календарь


Число дней между двумя датами

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating number of days between two dates
//  (c) Johna Smith, 1996
//
//  Method description:
//    We calculte two Julian dates and substract them
//    Result will be number of days between two dates
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

enum Month {January=1,February,March,April,May,June,July,August,
           September,October,November,December};

struct Date
{
 float day;
 Month month;
 int year;
};

double JulianDate(Date date)
{
 long int A,B,Y,M;

 if (date.month==January || date.month==February)
 {
   Y=date.year-1;
   M=date.month+12;
 } else
 {
   Y=date.year;
   M=date.month;
 }
 if ((date.year>1581) ||
     (date.year==1581 && date.month>October) ||
     (date.year==1581 && date.month==October && date.day>=15))
 {
   A=Y/100;
   B=2-A+A/4;
 } else B=0;

 return B+(long int)(365.25*Y)+(long int)(30.6001*(M+1))+date.day+1720994.5;
}

Date Revolution={17.5,October,1917};
Date Today={11,October,1996};

void main(void)
{
 printf("Days between %f/%d/%d and %f/%d/%d = %f",
        Revolution.day,Revolution.month,Revolution.year,
        Today.day,Today.month,Today.year,
        JulianDate(Today)-JulianDate(Revolution));
}        



Определение дня недели по дате

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating day of week by date
//  (c) Johna Smith, 1996
//
//  Method description:
//  Here we use special formulas different for January & February
//  and other monthes. Use this algorithm only for dates after
//  October, 15, 1528
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

enum Month {January=1,February,March,April,May,June,July,August,
           September,October,November,December};

struct Date
{
 float day;
 Month month;
 int year;
};

int DayOfWeek(Date date)
{
 float F;

 if (date.month<March)
   F=365*date.year+date.day+31*(date.month-1)+
     (int)((date.year-1)/4)-(int)(3*(int)((date.year-1)/100+1)/4);
 else
   F=365*date.year+date.day+31*(date.month-1)-(int)(0.4*date.month+2.3)+
     (int)(date.year/4)-(int)(3*(int)(date.year/100+1)/4);

 return  F-7*(int)(F/7)-1;
}


Date Today={11,October,1996};

void main(void)
{
 printf("%f/%d/%d is ",Today.day,Today.month,Today.year);
 switch (DayOfWeek(Today))
 {
   case -1:printf("Sunday\n");break;
   case 0:printf("Monday\n");break;
   case 1:printf("Tuesday\n");break;
   case 2:printf("Wednesday\n");break;
   case 3:printf("Thursday\n");break;
   case 4:printf("Friday\n");break;
   case 5:printf("Saturday\n");break;
 }
}



Вычисление даты по юлианскому календарю

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Julian date calculating
//  (c) Johna Smith, 1996
//
//  Method description:
//    Julian date is the number of days since afternoon (GMT) of
//    January, 1, 4713 BC.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

enum Month {January=1,February,March,April,May,June,July,August,
           September,October,November,December};

struct Date
{
 float day;
 Month month;
 int year;
};

double JulianDate(Date date)
{
 long int A,B,Y,M;

 if (date.month==January || date.month==February)
 {
   Y=date.year-1;
   M=date.month+12;
 } else
 {
   Y=date.year;
   M=date.month;
 }
 if ((date.year>1581) ||
     (date.year==1581 && date.month>October) ||
     (date.year==1581 && date.month==October && date.day>=15))
 {
   A=Y/100;
   B=2-A+A/4;
 }
 else
   B=0;

 return B+(long int)(365.25*Y)+(long int)(30.6001*(M+1))+date.day+1720994.5;
}

Date Revolution={17.5,October,1917};

void main(void)
{
 printf("%f/%d/%d = %f",Revolution.day,Revolution.month,
                        Revolution.year,JulianDate(Revolution));
}

PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 11:23 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



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


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

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (square root method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    From given matrix A we build two auxulary triangular matrixes S & D
//    Then solving equations Gy=b (where g[I][j]=s[j][I]*d[j]) and Sx=y
//    we obtain vector x.
//
//////////////////////////////////////////////////////////////////////////////

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

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

float matrix[N][N1]=
{{10.9,1.2,2.1,0.9,23.2},
 {1.2,11.2,1.5,2.5,38.1},
 {2.7,1.5,9.8,1.3,40.3},
 {0.9,2.5,1.3,12.1,58.2}
};

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
 register char i,j,k;
 float s[N][N],x[N],y[N],d[N],sum;

 // Printing given matrix
 ShowMatrix();

 // Building matrixes S and D
 for (j=0;j<N;j++)
   for (i=0;i<=j;i++)
   {
     sum=matrix[i][j];
     for (k=0;k<i;k++)
       sum-=s[k][i]*d[k]*s[k][j];
     if (i==j)
     {
       d[i]=(sum>0?1:-1);
       s[i][j]=sqrt(fabs(sum));
     }
     else s[i][j]=sum/(d[i]*s[i][i]);
   }

 // Solving equation Gy=b (G: g[I][j]=s[j][I]*d[j])
 for (i=0;i<N;i++)
 {
   y[i]=matrix[i][N];
   for (j=0;j<i;j++) y[i]-=s[j][i]*d[j]*y[j];
   y[i]/=s[i][i]*d[i];
 }

 // Solving equation Sx=y
 for (i=N-1;i>=0;i--)
 {
   x[i]=y[i];
   for (j=i+1;j<N;j++) x[i]-=s[i][j]*x[j];
   x[i]/=s[i][i];
 }

 // 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 (Nekrasov method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This is an iterative method of solving systems of linear equations.
//    Criteria for end of iterations is ||xn||-||x(n-1)||<epsilon, where
//    ||x|| is norm of vector x.
//
//////////////////////////////////////////////////////////////////////////////

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

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

float matrix[N][N1]=
     {{10.9,1.2,2.1,0.9,23.2},
      {1.2,11.2,1.5,2.5,38.1},
      {2.7,1.5,9.8,1.3,40.3},
      {0.9,2.5,1.3,12.1,58.2}
    };

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

// this function calculates ||x||
float norma(float *x)
{
 float sum=0;

 for (unsigned i=0;i<N;i++) sum+=fabs(*(x+i));
 return sum;
}

void main(void)
{
 // Variables declaration
 register char i,j;
 float x_current[N],x_next[N],sum1,sum2,epsilon;

 epsilon=1e-7;
 // Printing given matrix
 ShowMatrix();

 // Solving
 do
 {
   for (i=0;i<N;i++) x_current[i]=x_next[i];
   for (i=0;i<N;i++)
   {
     sum1=0;sum2=0;
     for(j=0;j<i;j++) sum1+=matrix[i][j]*x_next[j];
     for(j=i+1;j<N;j++) sum2+=matrix[i][j]*x_current[j];
     x_next[i]=(matrix[i][N]-sum1-sum2)/matrix[i][i];
   }
 } while (fabs(norma(x_current)-norma(x_next))>epsilon);

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



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

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (3-diagonal matrix)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This method is developed for 3-diagonal matrixes:
//    (a11 a12 0 ..............)
//    (a21 a22 a23 0 ..........)
//    (0   a32 a33 a34 0 ......) = M
//    (........................)
//    (.......... 0 an(n-1) ann)
//    Cause all elements except three diagonals is 0 we'll store
//    only non zero elements in the A,B,C arrays. Array D stores
//    vector B (M*X=B). We build two auxulary vectors P and Q and
//    then calculate vector X using the following recurrent formula:
//      X(k-1)=P(k)*X(k)+Q(k)
//
//////////////////////////////////////////////////////////////////////////////


#include <stdio.h>

#define N    20
#define N1   N+1

void main(void)
{
 // Variables declaration
 register int i,i1,k;  // using processor registers for better performance
 float a[N]={0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};
 float b[N]={2,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,2};
 float c[N]={-1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0};
 float d[N]={0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0};
 float x[N1],p[N1],q[N1],tmp;

 // Some initializations

 p[0]=0;
 q[0]=0;
 x[N]=0;
 c[N-1]=0; // to avoid overflow

 // Building vectors P and Q
 for (i=0;i<N;i++)
 {
   i1=i+1;
   if ((tmp=b[i]+a[i]*p[i])==0)
   {
     printf("Wrong method for such system...\n");
     return;
   };
   p[i1]=-c[i]/tmp;
   q[i1]=(d[i]-a[i]*q[i])/tmp;
 }

 // Building vector X
 for (k=N;k>0;k--) x[k-1]=p[k]*x[k]+q[k];

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



Метод Гаусса

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Gauss method)
//  (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)
//
//////////////////////////////////////////////////////////////////////////////

#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 tmp,x[N1];
 register short int i,j,k;

 // Printing given matrix
 ShowMatrix();
 // Main loop
 for (i=0;i<N;i++)
 {
   // 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 (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]);
}



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

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Zaidel method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: system of linear equations
//       a11 x1 + a12 x2 + ... + a1n xn = b1
//       a21 x1 + a22 x2 + ... + a2n xn = b2
//       ...
//       an1 x1 + an2 x2 + ... + ann xn = bn
//
//    We use following iterative formula to solve this system:
//
//                     1    i-1                N
//     xi(j+1)=xi(j)- --- ( SUM aik xk(j+1) + SUM aik xk(j) -bi)
//                    aii   k=1               k=i
//
//     where xi(j+1) is j+1-th iteration, xi(j) is j-th iteration
//
//////////////////////////////////////////////////////////////////////////////

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

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

float matrix[N][N1]=
{{10,1,1,12},
 {2,10,1,13},
 {2,2,10,14}
};
float z[N]={0,0,0};  // first iteration

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];
 short int i,j;
 int iterations=0,finish=0;

 // Printing given matrix
 ShowMatrix();

 while (!finish)
 {
   finish=1;
   for (i=0;i<N;i++)
   {
     x[i]=-matrix[i][N];
     for (j=0;j<N;j++) x[i]+=matrix[i][j]*z[j];
     // don't stop iterations until required accuracy is reached
     if (fabs(x[i]/matrix[i][i])>=epsilon) finish=0;
     x[i]=z[i]-x[i]/matrix[i][i];
     // next iteration
     z[i]=x[i];
   }
   iterations++;
 }

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

PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 12:23 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



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


Метод хорд

Код

//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (chords method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    given: f(x)=0, [a;b], f(a)*f(b)<0
//    find x0: f(x0)=0
//    We approximate curve by chord between its borders and get
//    point where the chord crosses X axe as a new border of
//    interval with the root. When interval becomes small enough
//    we can take one of its borders as a solution.
//
//////////////////////////////////////////////////////////////////////////////

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

// this function returns value of f(x)
double f(double x)
{
 // sin 2x - ln x = 0
 return sin(2*x)-log(x);
}

double Solve(double a, double b, double epsilon)
{
 double u,v,x,y;

 u=f(a); v=f(b);
 do
 {
   // x - new border
   x=(a*v-b*u)/(v-u);
   y=f(x);
   // determine whether x is left or right border
   // f(a)*f(b)<0 condition must be true
   if (y*u<0)
   {
     b=x;
     v=y;
   } else
   {
     a=x;
     u=y;
   }
 } while (b-a>epsilon);
 return x;
}

void main(void)
{
 printf("sin 2x - ln x = 0, [a;b]=[1.3;1.5]");
 printf("\nx0=%f",Solve(1.3,1.5,1e-9));
}



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

Код

//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (dichotomia method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    given: f(x)=0, [a;b], f(a)*f(b)<0
//    find x0: f(x0)=0
//    We split interval into two parts and select that part that
//    contains root of the equation by checking condition
//    f(a)*f(b)<0 We take middle point of the interval as a new border of
//    interval with the root. When interval becomes small enough
//    we can take one of its borders as a solution.
//
//////////////////////////////////////////////////////////////////////////////

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

// this function returns value of f(x)
double f(double x)
{
 // sin 2x - ln x = 0
 return sin(2*x)-log(x);
}

double Solve(double a, double b, double epsilon)
{
 double x;

 do
 {
   // x - new border
   x=(b+a)/2;
   // determine whether x is left or right border
   // f(a)*f(b)<0 condition must be true
   if (f(x)*f(a)<0) b=x;
   else a=x;
 } while (b-a>epsilon);
 return x;
}

void main(void)
{
 printf("sin 2x - ln x = 0, [a;b]=[1.3;1.5]");
 printf("\nx0=%f",Solve(1.3,1.5,1e-9));
}



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

Код

//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (Eitken-Steffenson method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This method is used for solving equations like x=f(x)
//    (see also iterations method)
//    Here we use the following iterations:
//     x(n+1)=(x0*x2-x1^2)/(x0-2x1+x2)
//
//////////////////////////////////////////////////////////////////////////////

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


// this function returns value of f(x)
double f(double x)
{
 // x = x + 0.37(sin 2x - ln x)
 return x+0.37*(sin(2*x)-log(x));
}

double Solve(double x, double epsilon)
{
 double x1,x2,xold,tmp;

 do
 {
   xold=x;
   x1=f(x);
   x2=f(x1);
   tmp=x-2*x1+x2;
   x=(x*x2-x1*x1)/tmp;
 } while (fabs(xold-x)>epsilon && tmp!=0);

 return x;
}

void main(void)
{
 printf("sin 2x - ln x = 0, [a;b]=[1.3;1.5]");
 printf("\nx0=%f",Solve(1.4,1e-6));
}



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

Код

//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (iterations method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This method is used for solving equations like x=f(x)
//    We should set first approximation of the root and then
//    make some iterations x=f(x). When the difference between
//    two consequental iterations will be less than epsilon we
//    can say that we found the root.
//    To apply this method to equations like F(x)=0 we should
//    transform it into iterational form x=f(x), where f(x):
//      1) defined on [a;b]
//      2) for every point on [a;b] exists f'(x)
//      3) for every x from [a;b] f(x) is in [a;b]
//      4) exists number q : |f'(x)|<=q<1 for all x from [a;b]
//
//////////////////////////////////////////////////////////////////////////////

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

// this function returns value of f(x)
double f(double x)
{
 // x = x + 0.37(sin 2x - ln x)
 return x+0.37*(sin(2*x)-log(x));
}

double Solve(double x, double epsilon)
{
 double x1;

 do
 {
   x1=x;
   x=f(x);
 } while (fabs(x1-x)>epsilon);

 return x;
}

void main(void)
{
 printf("sin 2x - ln x = 0");
 printf("\nx0=%f",Solve(1.4,1e-6));
}



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

Код

//////////////////////////////////////////////////////////////////////////////
//  Solving nonlinear equations (Newton method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Here we use the following iterations:
//     x(n+1)=xn-f(xn)/f'(xn)
//
//////////////////////////////////////////////////////////////////////////////

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


// this function returns value of f(x)
double f(double x)
{
 // sin 2x - ln x = 0
 return sin(2*x)-log(x);
}

double Solve(double x, double epsilon)
{
 double x1;

 do
 {
   x1=x;
   x=x-f(x)*epsilon/(f(x+epsilon)-f(x));
 } while (fabs(x1-x)>epsilon);

 return x;
}

void main(void)
{
 printf("sin 2x - ln x = 0, [a;b]=[1.3;1.5]");
 printf("\nx0=%f",Solve(1.4,1e-9));
}



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

Код

//////////////////////////////////////////////////////////////////////////////
//  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;
 }
}

PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 13:54 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



Интерполяция


Метод Эйткена

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Interpolation (using Eitken method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: several points of the unknown function
//    We should find values of this function between given points
//    Interpolation polynom isn't calculated directly:
//    To calculate values of function between given points we should
//    use the following formulas of linear interpolation:
//                  1     |y0  x0-x|
//     y(x,x0,x1)=------- |        |
//                (x1-x0) |y1  x1-x|
//
//                  1     |y0  x0-x|
//     y(x,x0,x2)=------- |        |
//                (x2-x0) |y2  x2-x|
//
//                      1   |y(x,x0,x1)  x1-x|
//     y(x,x0,x1,x2)=-------|                | and so on...
//                   (x2-x1)|y(x,x0,x2)  x2-x|
//
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

#define N  4

float x[N]={-0.12, 1.68, 3.41, 5.62};
float y[N]={0.324, -0.6, -0.4, -0.21};

float xvalues[N-1],f[N-1],tmp[N];
float xa;

void main(void)
{
 // initializing variables
 for(int i=0;i<N-1;i++)
 {
   xvalues[i]=(x[i+1]+x[i])/2; // we'll calculate function between given points
   f[i]=0;
 }

 // calculating function values between given points
 for (int z=0;z<N-1;z++)
 {
   for(int j=0;j<N;j++) tmp[j]=y[j];
   for(j=0;j<N-1;j++)
     for(i=j+1;i<N;i++)
       tmp[i]=((xvalues[z]-x[j])*tmp[i]-(xvalues[z]-x[i])*tmp[j])/(x[i]-x[j]);
   f[z]=tmp[N-1];
 }

 for (i=0;i<N-1;i++)
   printf("F(%f)= %f\n",xvalues[i],f[i]);
}



Полином Лагранжа

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Interpolation (using Lagrange polynom)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: several points of the unknown function
//    We should find values of this function between given points and
//    approximate  function by polynom.
//    Interpolation polynom calculated by the following formula:
//             N       (x-x0)...(x-x(i-1))(x-x(i+1))...(x-xn)
//       L(x)=SUM y(i)----------------------------------------
//            i=0     (xi-x0)..(xi-x(i-1))(xi-x(i+1))..(xi-xn)
//    To calculate values of function between given points we should
//    calculate value of interpolation polynom in these points.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

#define N  4

float x[N]={-0.12, 1.68, 3.41, 5.62};
float y[N]={0.324, -0.6, -0.4, -0.21};
float xvalues[N-1],f[N-1];
int i,j,z;
float l,k;
float pow[N],temp[N],tempaux[N];

void main(void)
{
 // initializing variables
 for(z=0;z<N-1;z++)
 {
   xvalues[z]=(x[z+1]+x[z])/2; // calculating function between given points
   f[z]=0;
 }
 for(z=0;z<N;z++) pow[z]=0;

 // calculating function values
 for(z=0;z<N-1;z++)
   for(i=0;i<N;i++)
   {
     l=1;
     for(j=0;j<N;j++)
       if(i!=j) l*=(xvalues[z]-x[j])/(x[i]-x[j]);
     l*=y[i];
     f[z]+=l;
   }

 // Determining interpolation polynom
 for(z=0;z<N;z++)
 {
   k=y[z];
   for(i=0;i<N;i++)
     if(i!=z) k/=(x[z]-x[i]);
   for(i=0;i<N;i++) temp[i]=0;
   temp[0]=k;
   for(i=0;i<N;i++)
     if(i!=z)
     {
       for(j=0;j<N-1;j++) tempaux[j+1]=temp[j];
       tempaux[0]=0;
       for(j=0;j<N;j++) tempaux[j]-=temp[j]*x[i];
       for(j=0;j<N;j++) temp[j]=tempaux[j];
     }
   for(i=0;i<N;i++) pow[i]+=temp[i];
 }

 // printing results
 for(z=0;z<N-1;z++) printf("F(%f)=%f\n",xvalues[z],f[z]);
 printf("\nPolynom: F(X)=");
 for(z=N-1;z>0;z--) printf("%fX^%d+\n",pow[z],z);
 printf("%f\n",pow[0]);
}



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

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Interpolation (using Newton polynom)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: several points of the unknown function
//    We should find values of this function between given points and
//    approximate  function by polynom.
//    Interpolation polynom is searched by the following formulas:
//
// P(x)=y1+(x-x1)f(x1,x2)+(x-x1)(x-x2)f(x1,x2,x3)+...              n
//                       +(x-x1)(x-x2)...(x-x(n-1))f(x1,x2,..,xn)=SUM Ak fi
//     i            i      i-j             k  k                   i=0
// fi=MUL (x-xj) = SUM aj x    ;   ak=(-1)  MUL xj
//    j=1          j=0                       j=1
//
//         1  l      m+1          k         m
//  a(k-l)=- SUM (-1)   a(k-l+m) SUM (-1/xp); a0=1
//         l m=1                 p=1
//
//
//    To calculate values of function between given points we should
//    calculate value of interpolation polynom in these points.
//
//////////////////////////////////////////////////////////////////////////////

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

#define N  4

float x[N]={-0.12, 1.68, 3.41, 5.62};
float Y[N]={0.324, -0.6, -0.4, -0.21};
float xvalues[N-1],f[N],a[N],y[N],pw,w,s;
int i,j,z,k,l,m,p;

void main(void)
{
 // initializing variables
 for(z=0;z<N-1;z++)
 {
   xvalues[z]=(x[z+1]+x[z])/2; // calculating function between given points
   f[z]=0;
 }
 for(z=0;z<N;z++) y[z]=Y[z];
 a[0]=1;
 f[N-1]=y[0];
 //           n-1
 // P   (x) = SUM A[k]
 //  n-1      k=0
 for(k=1;k<N;k++)
 {
   for(i=1;i<=N-k;i++) y[i-1]=(y[i]-y[i-1])/(x[i+k-1]-x[i-1]);
   //          k  k
   // a[k]=(-1)  MUL x[j]
   //            j=1
   pw=1;
   for(j=1;j<=k;j++) pw*=x[j-1];
   a[k]=(k%2==0?1:-1)*pw;

   //        1  l      m+1          k            m
   // a[k-l]=- SUM (-1)   a[k-l+m] SUM  (-1/x[p])
   //        l m=1                 p=1
   if (k!=1)
   {
     for(l=1;l<=k;l++)
     {
       w=0;
       for(m=1;m<=l;m++)
       {
         s=0;
         for(p=1;p<=k;p++) s+=pow(-1/x[p-1],m);
         w+=(m%2==0?-1:1)*a[k-l+m]*s;
       }
       a[k-l]=w/l;
     }
   }

   for(j=N;j>=N-k;j--) f[j-1]+=a[j-N+k]*y[0];
 }

 // printing results
 for(z=0;z<N-1;z++)
 {
   s=f[0];
   for(i=1;i<N;i++) s=s*xvalues[z]+f[i];
   printf("F(%f)=%f\n",xvalues[z],s);
 }
 printf("\nPolynom: F(X)=");
 for(z=0;z<N-1;z++) printf("%fX^%d+\n",f[z],N-z-1);
 printf("%f",f[N-1]);
}




Операции с полиномами


Вычисление производной

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating derivative of the polynom (analitically)
//  (c) Johna Smith, 1996
//
//  Method description:
//    P(x) - polynom,  deg P = N
//    P(x) = a0 x^n + a1 x^(n-1) + ... + an
//    calculating P'(x)
//    P'(x) = n a0 x^(n-1) + (n-1) a1 x^(n-2) +...
//
//////////////////////////////////////////////////////////////////////////////

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

void PrintPolynom(float *p, int n)
{
 for(int i=0;i<n;i++)
   printf("%+f x^%d",*(p+i),n-i-1);
 printf("\n");
}

#define N   4  // polynom degree

float Polynom[N]={1,2,3,4}; // x^3 + 2x^2 + 3x + 4 = P(x)

void main(void)
{
 // printing given polynom
 PrintPolynom(Polynom,N);

 // calculating derivative
 for(int i=0;i<N;i++)
   Polynom[i]*=(N-i-1);

 // printing derivative
 PrintPolynom(Polynom,N-1);
}



Деление полиномов

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Dividing polynoms (analitically)
//  (c) Johna Smith, 1996
//
//  Method description:
//    P(x) - polynom,  deg P = N
//    Q(x) - polynom,  deg Q = M, M<=N
//
//    P(x) = a0 x^n + a1 x^(n-1) + ... + an
//    Q(x) = b0 x^m + b1 x^(m-1) + ... + bm
//
//    calculating P(x)/Q(x)=S(x)+R(x)/Q(x), deg R < deg Q
//
//////////////////////////////////////////////////////////////////////////////

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

void PrintPolynom(float *p, int n)
{
 for(int i=0;i<n;i++)
   printf("%+f x^%d",*(p+i),n-i-1);
 printf("\n");
}

#define N    4  // deg P-1
#define M    2  // deg Q-1

float P[N]={1,1,1,-3}; // x^3 + x^2 + x - 3 = P(x)
float Q[M]={1,-1};     // x - 1 = Q(x)
float S[N-M+1];        // S(x)
               

// 0 here means constant that will be in          

// the integral I(P)=I+C

void main(void)
{
 // printing given polynoms
 printf("P(x):\n");
 PrintPolynom(P,N);
 printf("Q(x):\n");
 PrintPolynom(Q,M);

 for (int i=0;i<N+M-1;i++) S[i]=0;

 // dividing
 for(i=0;i<N-M+1;i++)
 {
   S[i]=P[i]/Q[0];
   for(int j=0;j<M;j++)
     P[i+j]-=S[i]*Q[j];
 }

 // printing result
 printf("\nresult:\n");
 PrintPolynom(S,N-M+1);
 // and remainder
 printf("remainder:\n");
 PrintPolynom(P+(N-M+1),M-1);
}



Интеграл от полинома

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating integral of the polynom (analitically)
//  (c) Johna Smith, 1996
//
//  Method description:
//    P(x) - polynom,  deg P = N
//    P(x) = a0 x^n + a1 x^(n-1) + ... + an
//    calculating Integral(P(x)dx)
//    Integral(P(x)dx) =  a0/(n+1) x^(n+1) + a1/n x^n +...
//
//////////////////////////////////////////////////////////////////////////////

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

void PrintPolynom(float *p, int n)
{
 for(int i=0;i<n;i++)
   printf("%+f x^%d",*(p+i),n-i-1);
 printf("\n");
}

#define N  4  // polynom degree

float Polynom[N+1]={8,6,4,1,0}; // 8x^3 + 6x^2 + 4x + 1 = P(x)

// 0 here means constant that will be in
// the integral I(P)=I+C

void main(void)
{
 // printing given polynom
 PrintPolynom(Polynom,N);

 // calculating integral
 for(int i=0;i<N;i++)
   Polynom[i]/=(N-i);

 // printing integral
 PrintPolynom(Polynom,N+1);
}



Умножение полиномов

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Multiplication of polynoms (analitically)
//  (c) Johna Smith, 1996
//
//  Method description:
//    P(x) - polynom,  deg P = N
//    Q(x) - polynom,  deg Q = M
//
//    P(x) = a0 x^n + a1 x^(n-1) + ... + an
//    Q(x) = b0 x^m + b1 x^(m-1) + ... + bm
//
//    calculating S(x)=P(x)*Q(x)
//
//////////////////////////////////////////////////////////////////////////////

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

void PrintPolynom(float *p, int n)
{
 for(int i=0;i<n;i++)
   printf("%+f x^%d",*(p+i),n-i-1);
 printf("\n");
}

#define N       4  // deg P-1
#define M       2  // deg Q-1

float P[N]={2,0,4,1}; // 2x^3 + 4x + 1 = P(x)
float Q[M]={1,-2};    // x - 2 = Q(x)
float S[M+N-1];       // S(x)

// 0 here means constant that will be in
// the integral I(P)=I+C

void main(void)
{
 // printing given polynoms
 PrintPolynom(P,N);
 PrintPolynom(Q,M);

 for (int i=0;i<N+M-1;i++) S[i]=0;

 // multiplying
 for(i=0;i<N;i++)
   for(int j=0;j<M;j++)
     S[i+j]+=P[i]*Q[j];

 // printing results
 PrintPolynom(S,N+M-1);
}

PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 14:18 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



Сжатие и шифрование


Компрессия/декомпрессия RLE

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Compressing/decompressing file (RLE technology)
//  (c) Johna Smith, 1996
//
//  Method description:
//  This compression method is based on coding sequences like AAA..A
//  We can write this sequence as two bytes N A, where A is
//  character and N is amount of its repetitions. So this method is
//  very useful for *.bmp files. But if there aren't such sequences
//  then each separate character will turn into two - 1 A 1 B ...
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <string.h>
#include <errno.h>

void usage(void)
{
 // show how to use this program if there isn't enough parameters
 printf("Usage: compress.exe filename1 COMPRESS|DECOMPRESS filename2\n");
 printf("filename1 - name of source file\n");
 printf("filename2 - name of destination file\n");
 printf("\nExample: compress.exe compress.bmp COMPRESS compress.rle\n");
}

FILE *input, *output;

void compress(void)
{
 unsigned char c1,c,length;
 unsigned long int rcounter=0,wcounter=0;

 fread(&c1,sizeof(char),1,input);
 do
 {
   c=c1;
   length=0;
   do
   {
     fread(&c1,sizeof(char),1,input);
     length++;
     rcounter++;
   } while (!feof(input) && c1==c && length!=255);
   // putting char and number of its repetitions to output
   fwrite(&c,sizeof(char),1,output);
   fwrite(&length,sizeof(char),1,output);
   wcounter++;
 } while (!feof(input));
 printf("%ld byte(s) compressed\n",rcounter);
 printf("ratio: %d\n",wcounter/rcounter);
}

void decompress(void)
{
 unsigned char c,length;
 unsigned long int rcounter=0,wcounter=0;

 fread(&c,sizeof(char),1,input);
 do
 {
   fread(&length,sizeof(char),1,input);
   for(int i=0;i<length;i++)
   {
     fwrite(&c,sizeof(char),1,output);
     wcounter++;
   }
   fread(&c,sizeof(char),1,input);
   rcounter+=2;
 } while (!feof(input));
 printf("%ld byte(s) decompressed\n",wcounter);
 printf("ratio: %d\n",wcounter/rcounter);
}

void main(int argc, char **argv)
{
 // copyleft;)_
 printf("Compress. copyleft {c} 1996 Johna Smith. freeware\n\n");
 // analysing parameters
 if (argc!=4)
 {
   usage();
   return;
 }
 if ((input=fopen(argv[1],"rb"))==NULL)
 {
   // error opening file
   printf("! ERR03: can't open source file :(\n");
   return;
 };
 if ((output=fopen(argv[3],"wb"))==NULL)
 {
   // error opening file
   fclose(input);
   printf("! ERR04: can't open destination file :(\n");
   return;
 };
 if (stricmp(argv[2],"COMPRESS")==0) compress();else
 if (stricmp(argv[2],"DECOMPRESS")==0) decompress();else
 {
   // must be COMPRESS or DECOMPRESS magic word there to know what to do
   printf("! ERR02: serious bug was found in second parameter :)\n");
   fclose(input);
   fclose(output);
   return;
 }
 fclose(input);
 fclose(output);
}



Шифрование/дешифрование сдвигом

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Encrypting/decrypting file (shift method)
//
//  (c) Johna Smith, 1996
//
//  Method description:
//    1) Set the base for randomizer (this is the main code for decrypting)
//    2) Get random number using the following formula: x=fract(11*x+Pi)
//    3) Crypt character using formula ch'=ch+rnd()*255
//    4) Repeat steps 3 & 4 for all chars in the text
//
//////////////////////////////////////////////////////////////////////////////

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

#define Pi  3.14159265358979323846

void usage(void)
{
 // show how to use this program if there isn't enough parameters
 printf("Usage: crypt.exe filename1 DECRYPT|ENCRYPT filename2 key\n");
 printf("filename1 - name of source file\n");
 printf("filename2 - name of destination file\n");
 printf("key - key code, must be a number\n");
 printf("\nExample: crypt.exe crypt.cpp ENCRYPT crypt.crp 123.45678\n");
}

float rnd(float _base=0)  // generates random number
{
 static float base;
 float tmp;

 // reinitializing randomizer if _base!=0
 if (_base!=0) base=_base;

 // generating random number
 tmp=11*base+Pi;
 base=tmp-(int)tmp;

 return base;
}

FILE *input, *output;

void encrypt(void)
{
 unsigned char ch;
 int tmp;
 unsigned long int counter=0;

 fread(&ch,sizeof(char),1,input);
 do
 {
   // encrypting char
   tmp=rnd()*255;
   if (ch+tmp>255) ch+=(tmp-255);
   else ch+=tmp;
   // putting char to output
   fwrite(&ch,sizeof(char),1,output);
   counter++;
   // getting char from input
   fread(&ch,sizeof(char),1,input);
 } while (!feof(input));
 printf("%ld byte(s) encrypted\n");
}

void decrypt(void)
{
 unsigned char ch;
 int tmp;
 unsigned long int counter=0;

 fread(&ch,sizeof(char),1,input);
 do
 {
   // decrypting char
   tmp=rnd()*255;
   if ((int)ch<tmp) ch-=(tmp-255);
   else ch-=tmp;
   // putting char to output
   fwrite(&ch,sizeof(char),1,output);
   counter++;
   // getting char from input
   fread(&ch,sizeof(char),1,input);
 } while (!feof(input));
 printf("%ld byte(s) decrypted\n");
}

void main(int argc, char **argv)
{
 // copyleft;)_
 printf("Crypt. copyleft {c} 1996 Johna Smith. freeware\n\n");
 // analysing parameters
 if (argc!=5)
 {
   usage();
   return;
 }
 rnd(atof(argv[4])); // setting up randomizer base
 if (errno==ERANGE)
 {
   // wrong key code
   printf("! ERR01: key code is out of range\n");
   return;
 }
 if ((input=fopen(argv[1],"rb"))==NULL)
 {
   // error opening file
   printf("! ERR03: can't open source file :(\n");
   return;
 };
 if ((output=fopen(argv[3],"wb"))==NULL)
 {
   // error opening file
   fclose(input);
   printf("! ERR04: can't open destination file :(\n");
   return;
 };
 if (stricmp(argv[2],"ENCRYPT")==0) encrypt();else
 if (stricmp(argv[2],"DECRYPT")==0) decrypt();else
 {
   // must be ENCRYPT or DECRYPT magic word there to know what to do
   printf("! ERR02: serious bug was found in second parameter :)\n");
   fclose(input);
   fclose(output);
   return;
 }
 fclose(input);
 fclose(output);
}



Шифрование/дешифрование ХОR-ом

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Encrypting/decrypting file (xor method)
//
//  (c) Johna Smith, 1996
//
//  Method description:
//    1) Get the key code (this is the main code for crypting)
//    2) Crypt part of text using formula str'=str xor keycode
//    3) Repeat step 2 for all parts of text
//
//   (a xor b) xor b = a -> encrypting==decrypting
//
//////////////////////////////////////////////////////////////////////////////

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

void usage(void)
{
 // show how to use this program if there isn't enough parameters
 printf("Usage: xorcrypt.exe filename1 filename2 key\n");
 printf("filename1 - name of source file\n");
 printf("filename2 - name of destination file\n");
 printf("key - key code, can be up to 20 characters (without spaces)\n");
 printf("\nExample: xorcrypt.exe xorcrypt.cpp xorcrypt.crp SAMPLEKEYCODE\n");
}

FILE *input, *output;
char keycode[20];
char keycodelength;

void encrypt(void)
{
 unsigned char ch;
 char j=0;
 unsigned long int counter=0;

 fread(&ch,sizeof(char),1,input);
 do
 {
   // processing char
   ch=(ch | keycode[j])-(ch & keycode[j]);  // a xor b = a or b-a and b
   j++;
   if (j==keycodelength) j=0;
   // putting char to output
   fwrite(&ch,sizeof(char),1,output);
   counter++;
   // getting char from input
   fread(&ch,sizeof(char),1,input);
 } while (!feof(input));
 printf("%ld byte(s) processed\n");
}

void main(int argc, char **argv)
{
 // copyleft;)_
 printf("XORCrypt. copyleft {c} 1996 Johna Smith. freeware\n\n");
 // analysing parameters
 if (argc!=4)
 {
   usage();
   return;
 }
 if (strlen(argv[3])>20)
 {
   // key code is too long
   printf("! ERR01: key code is too long\n");
   return;
 }
 strcpy(keycode,argv[3]);
 keycodelength=strlen(keycode);
 if ((input=fopen(argv[1],"rb"))==NULL)
 {
   // error opening file
   printf("! ERR03: can't open source file :(\n");
   return;
 };
 if ((output=fopen(argv[2],"wb"))==NULL)
 {
   // error opening file
   fclose(input);
   printf("! ERR04: can't open destination file :(\n");
   return;
 };
 encrypt();
 fclose(input);
 fclose(output);
}

PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 14:31 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



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


Метод Эйлера

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving differential equations (Eiler method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given differential equation y'=f(x,y) and starting conditions
//    x=x0, y(x0)=y0. We should find solution of this equation
//    at [a;b] interval.
//    y(i+1) calculated as y(i) + delta y(i)
//    delta y=h*f(x,y(i))
//
//    In this example y'=cos(y)+3x, [a;b]=[0;1], x0=0, y0=1.3
//
//////////////////////////////////////////////////////////////////////////////

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

const float a=0,b=1; // bound of the interval
const int num_points=10; // number of point to solve
float x0=0,y0=1.3; // initial conditions
int M=1;

float f(float x, float y)
{
 return (cos(y)+3*x); // y'=cos(y)+3*x
}

void calculate(int m,float *y)
{
 float x,yi,h;

 h=(b-a)/(num_points*m);
 yi=y0; x=x0;
 for (int i=0;i<num_points;i++)
 {
   for (int k=0;k<m;k++)
   {
     yi+=h*f(x,yi);
     x+=h;
   }
   *(y+i)=yi;
 }
}

void main(void)
{

 float yh[num_points],yh2[num_points];

 calculate(M,yh);
 calculate(2*M,yh2);  // doubled step for better accuracy
 // epsilon is difference between solutions with single
 // and double steps
 printf("X\t\tYH\t\tYH2\t\tEPSILON\n");
 for (int i=0;i<num_points;i++)
   printf("%f\t%f\t%f\t%f\n",(x0+((i+1)*(b-a))/num_points),
                                     yh[i],yh2[i],yh2[i]-yh[i]);
}



Метод Адамса

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving differential equations (Adams method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given differential equation y'=f(x,y) and starting conditions
//    x=x0, y(x0)=y0. We should find solution of this equation
//    at [a;b] interval.
//    We use Runge-Kutta method to process first 'num_starting_points'
//    points and then use Adams extrapolation to process other points.
//    In this example y'=x+y, [a;b]=[0;2], x0=0, y0=1
//
//////////////////////////////////////////////////////////////////////////////

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

const float a=0,b=2;             // bounds of the interval
const int num_points=10,         // number of points to solve
         num_starting_points=4; // number of points to solve with Runge-Kutta method
float x0=0,y0=1;                 // starting conditions

float f(float x, float y)
{
 return x+y;  // y'=x+y
}

// this function realises Runge-Kutta method for n starting points
void calculate(float *y)
{
 float k1,k2,k3,k4,x,yi,h;

 h=(b-a)/num_points;  // step
 yi=y0; x=x0;
 for (int i=0;i<num_starting_points;i++)
 {
   k1=h*f(x,yi);
   k2=h*f(x+h/2,yi+k1/2);
   k3=h*f(x+h/2,yi+k2/2);
   k4=h*f(x+h,yi+k3);
   yi+=(k1+2*k2+2*k3+k4)/6;
   x+=h;
   *(y+i+1)=yi;
 }
}

void main(void)
{

 float y[num_points+1],h;

 y[0]=y0;
 // apply Runge-Kutta method
 calculate(y);
 h=(b-a)/num_points;
 // extrapolating
 for (int i=num_starting_points;i<num_points;i++)
   y[i] = y[i-1]+h/24*(55*f(x0+(i-1)*h,y[i-1])-
            59*f(x0+(i-2)*h,y[i-2])+
            37*f(x0+(i-3)*h,y[i-3])-
            9*f(x0+(i-4)*h,y[i-4]));

 printf("X\t\tY\t\tExact solution\n");
 for (i=0;i<num_points;i++)
   printf("%f\t%f\t%f\n",(x0+i*h),y[i],(2*exp(x0+i*h)-(x0+i*h)-1));
}



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

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving differential equations (Eiler-Coshi method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given differential equation y'=f(x,y) and starting conditions
//    x=x0, y(x0)=y0. We should find solution of this equation
//    at [a;b] interval.
//    y(i+1) calculated as y(i) + delta y(i)
//    delta y=h*(f(x,y(i))+f(x,z))/2, z=y(i)+h*(f(x(i),y(i))
//
//    In this example y'=cos(y)+3x, [a;b]=[0;1], x0=0, y0=1.3
//
//////////////////////////////////////////////////////////////////////////////

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

const float a=0,b=1;     // bound of the interval
const int num_points=10; // number of point to solve
float x0=0,y0=1.3;       // initial conditions
int M=1;

float f(float x, float y)
{
 return (cos(y)+3*x); // y'=cos(y)+3*x
}

void calculate(int m,float *y)
{
 float x,yi,h,z;

 h=(b-a)/(num_points*m);
 yi=y0; x=x0;
 for (int i=0;i<num_points;i++)
 {
   for (int k=0;k<m;k++)
   {
     z=yi+h*f(x,yi);
     yi+=h*(f(x,yi)+f(x,z))/2;
     x+=h;
   }
   *(y+i)=yi;
 }
}

void main(void)
{

 float yh[num_points],yh2[num_points];

 calculate(M,yh);
 calculate(2*M,yh2);  // doubled step for better accuracy
 // epsilon is difference between solutions with single
 // and double steps
 printf("X\t\tYH\t\tYH2\t\tEPSILON\n");
 for (int i=0;i<num_points;i++)
   printf("%f\t%f\t%f\t%f\n",(x0+((i+1)*(b-a))/num_points),
                                      yh[i],yh2[i],yh2[i]-yh[i]);
}



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

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving differential equations (Runge-Kutta method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given differential equation y'=f(x,y) and starting conditions
//    x=x0, y(x0)=y0. We should find solution of this equation
//    at [a;b] interval.
//    y(i+1) calculated as y(i) + delta y(i)
//    delta y =1/6 k1 + 1/3 k2 + 1/3 k3 + 1/6 k4
//    k1=hf(x,y)
//    k2=hf(x+h/2,y+k1/2)
//    k3=hf(x+h/2,y+k2/2)
//    k4=hf(x+h,y+k3)
//
//    In this example y'=cos(x+y)+x-y, [a;b]=[0;1], x0=0, y0=0
//
//////////////////////////////////////////////////////////////////////////////

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

const float a=0,b=1;     // bound of the interval
const int num_points=10; // number of point to solve
float x0=0,y0=0;         // initial conditions
int M=1;

float f(float x, float y)
{
 return (cos(x+y)+x-y); // y'=cos(x+y)+x-y
}

void calculate(int m,float *y)
{
 float k1,k2,k3,k4,x,yi,h;

 h=(b-a)/(num_points*m);
 yi=y0; x=x0;
 for (int i=0;i<num_points;i++)
 {
   for (int k=0;k<m;k++)
   {
     // calculatin coefficients k
     k1=h*f(x,yi);
     k2=h*f(x+h/2,yi+k1/2);
     k3=h*f(x+h/2,yi+k2/2);
     k4=h*f(x+h,yi+k3);
     yi+=(k1+2*k2+2*k3+k4)/6;
     x+=h;
  }
  *(y+i)=yi;
 }
}

void main(void)
{
 float yh[num_points],yh2[num_points];

 calculate(M,yh);
 calculate(2*M,yh2);  // doubled step for better accuracy
 // epsilon is difference between solutions with single
 // and double steps
 printf("X\t\tYH\t\tYH2\t\tEPSILON\n");
 for (int i=0;i<num_points;i++)
   printf("%f\t%f\t%f\t%f\n",(x0+((i+1)*(b-a))/num_points),
                                   yh[i],yh2[i],(yh2[i]-yh[i]));
}

PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 14:46 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



Вычисление CRC


CRC-16 DOS

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculacting CRC-16 (DOS) of data block
//  (c) Johna Smith, 1996
//
//  Method description:
//    CRC is used to check data sequences for errors
//    CRC-16 DOS is DOS modification of 16-bit control code (CRC)
//    To calculate CRC-16 DOS
//    1) CRC=0
//    2) Call CRC16DOS for each byte of data block
//    3) Call CRC16DOS(0) twice
//    Function CRC16DOS recieves current CRC value and returns updated value
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

int CRC16DOS(char c, unsigned int crc)
{
 unsigned int CRC_MASK=0xa001;

 asm
 {
   mov     al,c      // register AL contains current byte value (8 bits)
   mov     dx,crc    // register DX contains current CRC value (16 bits)
   push    cx        // store CX register in stack
   mov     cx,8      // put number 8 into CX register (CX=8)
 CRC_Loop:
   ror     al,1
   //  ROR  AL,1 - shifting bits of AL to the right, putting lowest
   //  bit to highest bit and to carry flag
   //  Example: 10010110b -> 01001011b and Carry Flag is set to 0
   rcr     dx,1
   //  RCR  DX,1  -  shifting  bits  of  DX to the right (like >> operation)
   //  and putting value from Carry Flag to the highetst bit
   //  Example: 10010110b (Carry Flag=1) -> 11001011 (Carry Flag=0)
   jnc     CRC2  // Goto CRC2 if Carry Flag is 0
   xor     dx,CRC_MASK // DX=DX xor CRC_MASK=(DX | CRC_MASK) - (DX & CRC_MASK)
 CRC2:
   loop    CRC_Loop // CX--; if (CX!=0) goto CRC_Loop
   pop     cx       // restore CX register from stack
   mov     crc,dx   // CRC=DX
 }

 return crc;
}

void main(void)
{
 int crc=0;   // initially CRC=0
 char data[]="This is data sample.";

 // Call CRC16DOS for each byte in the block
 for (int i=0;i<sizeof(data)/sizeof(char);i++) crc=CRC16DOS(data[i],crc);
 // Calculate CRC(0) twice
 crc=CRC16DOS(0,crc);
 crc=CRC16DOS(0,crc);

 // Print calculated value
 printf("%u\n",crc);
}



CRC-16 CCITT

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculacting CRC-16 (CCITT) of data block
//  (c) Johna Smith, 1996
//
//  Method description:
//    CRC is used to check data sequences for errors
//    CRC-16 CCITT is CCITT modification of 16-bit control code (CRC)
//    To calculate CRC-16 CCITT
//    1) CRC=0
//    2) Call CRC16CCITT for each byte of data block
//    3) Call CRC16CCITT(0) twice
//    Function CRC16CCITT recieves current CRC value and returns updated value
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

int CRC16CCITT(char c, unsigned int crc)
{
 unsigned int CRC_MASK=0x1021;

 asm
 {
   mov     al,c     // register AL contains current byte value (8 bits)
   mov     dx,crc   // register DX contains current CRC value (16 bits)
   push    cx       // store CX register in stack
   mov     cx,8     // put number 8 into CX register (CX=8)
 CRC_Loop:
   rol     al,1
   //  ROL  AL,1 - shifting bits of AL to the left, putting highest
   //  bit to the lowest bit and to carry flag
   //  Example: 11010110b -> 10101101b and Carry Flag is set to 1
   rcl     dx,1
   //  RCL  DX,1  -  shifting  bits  of  DX to the left (like << operation)
   //  and putting value from Carry Flag to the lowest bit
   //  Example: 01010110b (Carry Flag=1) -> 10101101b (Carry Flag=0)
   jnc     CRC2     // Goto CRC2 if Carry Flag is 0
   xor     dx,CRC_MASK // DX=DX xor CRC_MASK=(DX | CRC_MASK) - (DX & CRC_MASK)
 CRC2:
   loop    CRC_Loop // CX--; if (CX!=0) goto CRC_Loop
   pop     cx       // restore CX register from stack
   mov     crc,dx   // CRC=DX
 }
 return crc;
}

void main(void)
{
 int crc=0;   // initially CRC=0
 char data[]="This is data sample.";

 for (int i=0;i<sizeof(data)/sizeof(char);i++) crc=CRC16CCITT(data[i],crc);
 // Calculate CRC(0) twice
 crc=CRC16CCITT(0,crc);
 crc=CRC16CCITT(0,crc);

 // Print calculated value
 printf("%u\n",crc);
}



CRC-32

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculacting CRC-32 of data block
//  (c) Johna Smith, 1996
//
//  Method description:
//    CRC is used to check data sequences for errors
//    CRC-32 is 32-bit control code (CRC)
//    To calculate CRC-32:
//    1) CRC=-1
//    2) Call CRC32 for each byte of data block
//    3) CRC=~CRC
//    Function CRC32 recieves current CRC value and returns updated value
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

long int CRC32(char c, long int crc)
{
 int CRC_MASK=0x8320; // CRC-32 mask lo word
 int CRC_MASK1=0xEDB8; // CRC-32 mask hi word
 char CRC_FF=CRC_MASK & 0xFF;
 int d,b;

 d=(int)crc;  // d is lo word of CRC-32
 b=crc>>16;   // b is hi word of CRC-32
 asm
 {
   mov al,c   // register AL contains current byte value (8 bits)
   mov bx,b   // register BX contains hi word of current CRC value (16 bits)
   mov dx,d   // register BX contains lo word of current CRC value (16 bits)
   push cx    // store CX register in stack
   push ax    // store AX register in stack
   xor al,dl  // AL=AL xor DL=(AL | DL) - (AL & DL)
   mov cx,8   // put number 8 into CX register (CX=8)
 CRC_Loop:
   shr bx,1   // bx>>=1 (shifting bits of BX to the right)
   rcr dx,1
   //  RCR  DX,1  -  shifting  bits  of  DX to the right (like >> operation)
   //  and putting value from Carry Flag to the highetst bit
   //  Example: 10010110b (Carry Flag=1) -> 11001011 (Carry Flag=0)
   shr al,1   // al>>=1 (shifting bits of AL to the right)
   jnc CRC2   // Goto CRC2 if Carry Flag is 0
   xor dx,CRC_MASK  // DX=DX xor CRC_MASK=(DX | CRC_MASK) - (DX & CRC_MASK)
   xor bx,CRC_MASK1 // BX=BX xor CRC_MASK1=(BX | CRC_MASK1) - (BX & CRC_MASK1)
   xor al,CRC_FF // AL=AL xor CRC_MFF=(AL | CRC_FF) - (AL & CRC_FF)
 CRC2:
   loop CRC_Loop // CX--; if (CX!=0) goto CRC_Loop
   pop ax        // restore AX register from stack
   pop cx        // restore CX register from stack
   mov b,bx      // b=BX
   mov d,dx      // d=DX
 }
 // combining lo word and hi word to unsigned long crc
 crc=( (unsigned long)( ((unsigned)(d)) | (((unsigned long)((unsigned)(b)))<<16)) );

 return crc;
}

void main(void)
{
 unsigned long int crc32=-1; // initially CRC=-1
 char data32[]={1,2,3,4,5,6,7,8,9,0};

 // Call CRC32 for each byte of data block
 for (int i=0;i<sizeof(data32)/sizeof(char);i++) crc32=CRC32(data32[i],crc32);
 crc32=~crc32;

 // Print calculated value
 printf("%lu",crc32);
}

PM WWW ICQ   Вверх
podval
Дата 23.12.2004, 15:12 (ссылка) |    (голосов:1) Загрузка ... Загрузка ... Быстрая цитата Цитата


Где я? Кто я?
****


Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб

Репутация: 18
Всего: 62



Разное


Число размещений из M по N

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating number of possible allocations
//  (c) Johna Smith, 1996
//
//  Method description:
//   Number of possible allocations from n elements by m-element groups is
//   equal to   m     n!
//             A  = ------
//              n   (n-m)!
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

double factorial (int n)
{
 double factorial=1;

 for (int i=2;i<=n;i++)
 {
   factorial=factorial*i;
 }

 return factorial;
}

double Anm(int n, int m)
{
 return (factorial(n)/factorial(n-m));
}

void main(void)
{
 printf("Number of allocations from 10 elements by 5 is Anm = %f\n",Anm(10,5));
}




Число сочетаний из N по M

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating number of possible combinations
//  (c) Johna Smith, 1996
//
//  Method description:
//    Number of possible combinations of n elements to m-element groups is
//   equal to   m      n!
//             C  = --------
//              n   m!(n-m)!
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

double factorial (int n)
{
 double factorial=1;

 for (int i=2;i<=n;i++)
 {
   factorial=factorial*i;
 }

 return factorial;
}

double Cnm(int n, int m)
{
 return (factorial(n)/(factorial(m)*factorial(n-m)));
}

void main(void)
{
 printf("Number of combination from 10 elements by 5 is Cnm = %f\n",Cnm(10,5));
}



Найти дробную часть Х

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating {x}, where {x} is the fractional part of x
//  (c) Johna Smith, 1996
//
//  Method description:
//   Calculate integer part of x ([x]):
//    1) Assume that integer part c is 0
//    2) If the given value x is positive then increase c while it is less than x
//       if the given value x is negative then decrease c while it is greater than x
//   Calculate the fractional part of x:
//    3) {x}=x-[x]
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

float Frac(float x)
{
 int c=0;

 // calculating integer part of x, [x]=c
 if (x>=0)
   while (x>=c+1) c++;
 else
   while (x<c-1) c--;
 // {x}=x-[x]

 return x-c;
}

void main(void)
{
 printf("Pi=%f {Pi}=%f\n",3.1415926536,Frac(3.1415926536));
 printf("{%f}=%f",-1.234567,Frac(-1.234567));
}



Округлить Х к ближайшему целому, меньшему, чем Х

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating [x], where [x] is the nearest to x integer that is less than x
//  [3.1415]=3 and [-3.1415]=-4 (not -3, because -3 is greater than -3.1415
//  (c) Johna Smith, 1996
//
//  Method description:
//    1) Assume that integer part c is 0
//    2) If the given value x is positive then increase c while it is less than x
//       if the given value x is negative then decrease c while it is greater than x
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

int Integer(float x)
{
 int c=0;

 if (x>=0)
   while (x>=c+1) c++;
 else
   while (x<c) c--;
   
 return c;
}

void main(void)
{
 printf("Pi=3.1415926536 [Pi]=%d\n",Integer(3.1415926536));
 printf("[-1234.5678]=%d",Integer(-1234.5678));
}



Факториал (решение "в лоб")

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Factorial calculation (exact solution)
//  (c) Johna Smith, 1996
//
//  Method description:
//   We just use definition of factorial n!=1*2*3*...*n
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

double factorial (int n)
{
 double factorial=1;

 for (int i=2;i<=n;i++)
 {
   factorial=factorial*i;
 }

 return factorial;
}

void main(void)
{
 printf("5! = %f\n35! = %f\n100! = %f\n",factorial(5),
                             factorial(35),factorial(100));
}



Факториал (через логарифм)

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Factorial calculation (logarithm solution)
//  (c) Johna Smith, 1996
//
//  Method description:
//   We just use logarithmic method:
//            sum(k=1;n) lg(k)
//     n! = 10
//
//////////////////////////////////////////////////////////////////////////////

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

long double factorial (int n)
{
 long double factorial;
 double sum=0;

 for (int i=1;i<=n;i++) sum+=(log(i)/log(10));
 factorial=powl(10,sum);

 return factorial;
}

void main(void)
{
 printf("5! = %Lf\n35! = %Lf\n100! = %Lf\n",factorial(5),
                               factorial(35),factorial(100));
}



Квадратный корень (метод Ньютона)

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Square root calculating. Newton method
//  (c) Johna Smith, 1996
//
//  Method description:
//    0) Given: A      Find: A^(1/2)
//    1) Find out first approximation x0 of the sqare root of A
//    2) x(n+1)=(x(n)+A/x(n))/2 while x(n+1)-x(n)<=1
//
//////////////////////////////////////////////////////////////////////////////

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

// recieving value of A and first approximation of the sqare root

float NewtonSqrt(float A, float x)
{
 float xold;

 do
 {
   xold=x;
   x=(x+A/x)/2.0;
 } while (x!=xold);
 return x;
}

void main(void)
{
 printf("C++ Sqrt(13)=%f,  Newton Sqrt(13)=%f.",sqrt(13),NewtonSqrt(13,3.5));
}



Наибольший общий делитель (алгоритм Евклида)

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Calculating NOD using Evklid algorithm
//  (c) Johna Smith, 1996
//
//  Method description:
//    NOD(a,b) is greatest integer number with the following property:
//     a/NOD=m, b/NOD=n, where m,n - are integer
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

int NOD(int a,int b)
{
 int tmp;

 // on each step we make following modifications
 // a   -> b
 // b%a -> a, where b%a is remainder of b/a
 while (a!=0)
 {
   tmp=a;
   a=b%a;
   b=tmp;
 }

 return b;
}

void main(void)
{
 printf("NOD(1532,24) = %d\n",NOD(1532,24));
}



Нахождение простых чисел

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Finding prime numbers
//  (c) Johna Smith, 1996
//
//  Method description:
//   We take a number and try to divide it. If we can divide it
//   without remainder - this is not prime number.
//   We can take into account only odd numbers, because we can
//   divide all even number by 2. Also we can store all prime
//   numbers that are already found in an array and try to divide
//   all new numbers only by numbers from this array.
//   If we want to find all prime numbers less than N the size of
//   the array should be sqrt(N)/2
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

#define N   160  // so we can find all prime numbers that are less than 100000
#define M   25   // check all numbers less than 250

int Simple[N];
int k=0;
enum {yes,no} simple;

void main(void)
{
 // it's easy: 2 and 3 are prime
 if (M>=2) printf("2\n"); // 2 is simple 'cause we can divide it only by itself and 1
 Simple[k++]=2;
 if (M>=3) printf("3\n");
 Simple[k++]=3;

 // but what we can say about other numbers:
 for(int i=5;i<=M;i+=2)
 {
   simple=yes;
   for(int j=0;j<k;j++)
   {
     if (Simple[j]*Simple[j]>i) break; // other Simple[j] is too big for i
     if ((i%Simple[j])==0) simple=no; // there's no remainder - not prime
   }
   if (simple==yes)
   {
     printf("%d\n",i);
     Simple[k++]=i;
   }
 }
}



Нахождение простых сомножителей числа

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Searching prime multipliers
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: number A
//    We should find a1,a2,a3,..an with the following properties:
//    1) ai is prime numbers
//    2) a1*a2*...*an=A
//
//    We can do it by consequitive dividing A by 2 and all prime numbers
//    while sqrt(A)+1<=d, where d is prime number
//
//////////////////////////////////////////////////////////////////////////////

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

unsigned long int A=1234567890L;

void main(void)
{
 unsigned long int b;

 while (A%2==0 && A!=1)
 {
   // divide by 2 while it is possible
   printf("2 ");
   A/=2;
 }
 b=3;
 while (b<sqrt(A)+1)
 {
   if (A%b==0)
   {
     // dividing by prime number
     printf("%ld ",b);
     A/=b;
   } else b+=2;
 }
 printf("%ld\n",A);
}



Быстрое возведение в степень x^y

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Quick powering
//  (c) Johna Smith, 1996
//
//  Method description:
//    Calculate A to the power of k, A is real, k is natural, k>0
//    If k is even we can use the following formula:
//                              k/2
//                   k   (   2 )
//              b = A  = ( A   )
//    and assume k=k/2; A=A*A;
//    If k is odd we assume k=k-1, b=b*A and k becomes even
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

void main(void)
{
 long double A = 8.0;
 long int    k = 19;
 long double b=1;

 printf("%Lf^%ld = ",A,k);
 while (k!=0)
 {
   if (k%2!=0) b*=A;
   k/=2;
   A*=A;
 }
 printf("%Lf\n",b);
}



Быстрое вычисление корня n-й степени x^(1/y)

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Quick powering
//  (c) Johna Smith, 1996
//
//  Method description:
//    Calculate A to the power of 1/k, A is real, k is natural, k>0
//    We can use iterations:
//               k-1
//    x    = (a/x   +(k-1)x  )/k
//     n+1       n         n
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

long double Root(long double a, long int k)
{
 long double xn=(a+k-1)/k; // first iteration
 long double x,xk;

 do
 {
   x=xn;
   xk=1; // x^(k-1)
   for (long int i=0;i<k-1;i++) xk*=x;
   xn=(a/xk+(k-1)*x)/k;
 } while (x!=xn); // while there's difference between two consequitive iterations
 return x;
}

void main(void)
{
 long double A=1048576L;
 long int k=20;

 printf("%Lf to the power of %Lf is %Lf",A,(long double)(1.0/k),Root(A,k));
}



Генерация случайных чисел с нормальным распределением

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Generating random numbers with normal distribution
//  (c) Johna Smith, 1996
//
//  Method description:
//    1) Select the base for randomizer
//    2) Generate two random numbers from 0 to 1
//    3) Calculate two randoms with normal distribution, return one of them
//       and remember other to return next time without generating new randoms
//  This method gives us random numbers with dispersion=1 and average=0  
//
//////////////////////////////////////////////////////////////////////////////

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

#define Pi  3.14159265358979323846

float nrnd(float _base=0)
{
 static float base;
 static float tact=1;
 static float second_value;
 float factor,argument,tmp;

 // reinitializing randomizer if _base!=0
 if (_base!=0) base=_base;

 if (tact==1)
 {
   tact=2;
   // generating random from 0 to 1
   tmp=11*base+Pi;
   base=tmp-(int)tmp;
   // calculating factor
   factor=sqrt(-2.0*log(base));
   // generating random from 0 to 1
   tmp=11*base+Pi;
   base=tmp-(int)tmp;
   // calculating argumeent
   argument=2.0*Pi*base;
   // generating final random numbers with normal distribution
   second_value=factor*sin(argument);
   return factor*cos(argument);
 }
 else
 {
   tact=1;
   return second_value;
 }
}

void main(void)
{
 // Setting the base for randiomizer (not 0)
 nrnd(987.654321);
 // Generating random numbers
 for (int i=0;i<10;i++) printf("%d. %f\n",i+1,nrnd());
}



Генерация случайных чисел в интервале 0...1

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Generating random number between 0 and 1
//  (c) Johna Smith, 1996
//
//  Method description:
//    1) Select the base for randomizer
//    2) Get random number using the following formula: x=fract(11*x+Pi)
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

#define Pi  3.14159265358979323846

float rnd(float _base=0)
{
 static float base;
 float tmp;

 // reinitializing randomizer if _base!=0
 if (_base!=0) base=_base;

 // generating random number
 tmp=11*base+Pi;
 base=tmp-(int)tmp;

 return base;
}

void main(void)
{
 // Setting the base for randiomizer (not 0)
 rnd(123.456789);
 // Generating random numbers
 for (int i=0;i<10;i++) printf("%d. %f\n",i+1,rnd());
}



Собственые числа и собственные вектора матрицы

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Searching own number and own vector for matrix A
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: matrix A
//    We need to find max lambda and vector X: A*X=lambda*X
//    In this algorythm we'll use iterations to find vector X.
//                         n
//    On each step lambda=SUM (A*x(i),x(i)), where (..) means scalar multiplication
//                        i=1
//
//////////////////////////////////////////////////////////////////////////////

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

#define N   3   // size of matrix

float matrix[N][N]={
  {10,5,6},
  {5,20,4},
  {6,4,30}
 };
float x_current[N]={1,1}; // first iteration
const unsigned m=300;     // number of iterations

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

void main(void)
{
 // Variables declaration
 int i,j,k;
 float A_x_current[N],lambda,alpha;

 // all iterations
 for(k=0;k<m;k++)
 {
   lambda=0;
   alpha=0;
   // multiplying matrix A by vector x_current
   for (i=0;i<N;i++) A_x_current[i]=0;
   
   for (i=0;i<N;i++)  // A*X(k)
     for (j=0;j<N;j++)
       A_x_current[i]+=matrix[j][i]*x_current[j];
   
   for (i=0;i<N;i++)  // (A*X(k),X(k))
     lambda+=A_x_current[i]*x_current[i];

   // normalizing x_current
   for (i=0;i<N;i++)
     alpha+=A_x_current[i]*A_x_current[i]; // alpha(k)
   alpha=1/sqrt(alpha);

   for (i=0;i<N;i++) x_current[i]=alpha*A_x_current[i]; // X(k+1)
 }

 // Printing given matrix
 printf("Problem:\n");
 ShowMatrix();

 // Printing solution
 printf("\nLambda(max)=%f\n",lambda);
 printf("X=(");
 for (i=0;i<N;i++) printf("%f%s",x_current[i],(i==N-1?")":","));
}



Инвестиции, проценты, прибыли

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Investments, percents, incomes
//  (c) Johna Smith, 1996
//
//  Method description:
//    s0 - initial investment
//    s  - final sum on account
//    P  - sum on account increased by P % per year
//    N  - number of years
//   If we know three of these values we can calculate fourth value
//
//////////////////////////////////////////////////////////////////////////////

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

// this function calculates s0
float StartValue(float end_value, float percents_per_period, float periods)
{
 return end_value/pow(1+percents_per_period/100,periods);
}

// this function calculates s
float EndValue(float start_value, float percents_per_period, float periods)
{
 return start_value*pow(1+percents_per_period/100,periods);
}

// this function calculates P
float PercentsPerPeriod(float start_value, float end_value, float periods)
{
 return 100*(pow(end_value/start_value,1/periods)-1);
}

// this function calculates N
float Periods(float start_value, float end_value, float percents_per_period)
{
 return (log10(end_value)-log10(start_value))/log10(1+percents_per_period/100);
}

void main(void)
{
 // If we invest $5000 for 5 years and bank pays 3% per year
 // the final sum on the account will be...
 printf("%f\n",EndValue(5000,3,5));
 // How much should we invest if we want to get $30000 and we
 // can wait 22 years and bank pays 3% per year
 printf("%f\n",StartValue(30000,3,22));
 // How much years we will wait if we invest $20000 and want to
 // get $30000. Bank pays 3% per year
 printf("%f\n",Periods(20000,30000,3));
 // If we want to invest $15000, wait 18 years and get $30000
 // then we need to find bank which pays 'PercentsPerPeriod' % per year
 printf("%f\n",PercentsPerPeriod(15000,30000,18));
}



Регулярные инвестиции, проценты, прибыли

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Regular investments, percents, incomes
//  (c) Johna Smith, 1996
//
//  Method description:
//    s1 - investment per year (we invest each year 's1' dollars)
//    s  - final sum on account
//    P  - sum on account increased by P % per year
//    N  - number of years
//   If we know three of these values we can calculate fourth value
//
//////////////////////////////////////////////////////////////////////////////

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

// this function calculates s1
float InvestmentValue(float end_value, float percents_per_period, float periods)
{
 float I=percents_per_period/100;

 return end_value*I/(pow(1+I,periods)-1);
}

// this function calculates s
float EndValue(float investment_value, float percents_per_period, float periods)
{
 float I=percents_per_period/100;

 return investment_value*(pow(1+I,periods)-1)/I;
}

// this function calculates P
float PercentsPerPeriod(float investment_value, float end_value, float periods)
{
 float i=0, h=0.01;

 while (fabs(h)>2.5e-4)
 {
   do i+=h; while ((end_value*i/(pow(1+i,periods)-1)-investment_value)/h>=0);

   h=-h/4;
 }

 return i*100;
}

// this function calculates N
float Periods(float investment_value, float end_value, float percents_per_period)
{
 float I=percents_per_period/100;

 return log(end_value*I/investment_value+1)/log(1+I);
}

void main(void)
{
 // If we will invest $1000 every year during 10 years and bank pays 3%
 // per year the final sum on the account will be...
 printf("%f\n",EndValue(1000,3,10));

 // How much years we will wait if we invest $1000 every year and want to
 // get $30000. Bank pays 3% per year
 printf("%f\n",Periods(1000,30000,3));

 // How much should we invest every year if we want to get $30000 and we
 // can wait 20 years and bank pays 3% per year
 printf("%f\n",InvestmentValue(30000,3,20));

 // If we want to invest $1000 every year during 20 years and get $30000
 // then we need to find bank which pays 'PercentsPerPeriod' % per year
 printf("%f\n",PercentsPerPeriod(1000,30000,20));
}

PM WWW ICQ   Вверх
REFRESH123
  Дата 10.4.2010, 18:51 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 1
Регистрация: 8.9.2009

Репутация: нет
Всего: нет




Модератор: Сообщение скрыто.

PM MAIL   Вверх
Google
  Дата 15.11.2019, 02:36 (ссылка)  





  Вверх
Страницы: (3) Все 1 [2] 3 
Закрытая темаСоздание новой темы Создание опроса
Правила форума "Алгоритмы"

maxim1000

Форум "Алгоритмы" предназначен для обсуждения вопросов, связанных только с алгоритмами и структурами данных, без привязки к конкретному языку программирования и/или программному продукту.


Если Вам понравилась атмосфера форума, заходите к нам чаще! С уважением, maxim1000.

 
0 Пользователей читают эту тему (0 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема | Алгоритмы | Следующая тема »


 




[ Время генерации скрипта: 0.3910 ]   [ Использовано запросов: 21 ]   [ GZIP включён ]


Реклама на сайте     Информационное спонсорство

 
По вопросам размещения рекламы пишите на vladimir(sobaka)vingrad.ru
Отказ от ответственности     Powered by Invision Power Board(R) 1.3 © 2003  IPS, Inc.