home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / picalcu.zip / EULER.C < prev    next >
Text File  |  2002-09-20  |  4KB  |  237 lines

  1. /* Determining of the  Euler constant e - 16 bit version */
  2.  
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include <malloc.h>
  6. #include <ctype.h>
  7.  
  8. #define DECMAX 157761L
  9. /* 287188 would be allowed for the  WORD format of the divisor */
  10.  
  11. typedef struct
  12. { unsigned int *fraction;
  13.   unsigned int valid;
  14. } LONGFRACT;
  15.  
  16. LONGFRACT e, f;
  17. unsigned int digits;
  18.  
  19. void divlongf (unsigned int d)
  20. {
  21.   _asm
  22.   {
  23.     mov bx,SEG f
  24.     mov es,bx
  25.     mov bx,OFFSET f
  26.     mov si,es:[bx].valid
  27.     les bx,es:[bx].fraction
  28.     mov cx,d
  29.     xor dx,dx
  30.  
  31.     dec si
  32.     shl si,1
  33.     mov ax,es:[bx+si]
  34.     div cx
  35.     and ax,ax
  36.     jnz div_2
  37.     mov bx,SEG f
  38.     mov es,bx
  39.     mov bx,OFFSET f
  40.     dec word ptr es:[bx].valid
  41.     les bx,es:[bx].fraction
  42.     jmp div_2
  43.  
  44.   div_1:
  45.     mov ax,es:[bx+si]
  46.     div cx
  47.   div_2:
  48.     mov es:[bx+si],ax
  49.     sub si,2
  50.     jnc div_1
  51.   }
  52. }
  53.  
  54.  
  55. void addlongf ()
  56. {
  57.   _asm
  58.   {
  59.     push ds
  60.     push bp
  61.     mov bx,SEG e
  62.     mov es,bx
  63.     mov bx,OFFSET e
  64.     les bx,es:[bx].fraction
  65.     mov bp,SEG f
  66.     mov ds,bp
  67.     mov bp,OFFSET f
  68.     mov cx,ds:[bp].valid
  69.     lds bp,ds:[bp].fraction
  70.     xor si,si
  71.     clc
  72.  
  73.   add_1:
  74.     mov ax,ds:[bp+si]
  75.     adc es:[bx+si],ax
  76.     inc si
  77.     inc si
  78.     dec cx
  79.     jnz add_1
  80.  
  81.     jnc add_3
  82.   add_2:
  83.     inc word ptr es:[bx+si]
  84.     jnz add_3
  85.     inc si
  86.     inc si
  87.     jmp add_2
  88.   add_3:
  89.     pop bp
  90.     pop ds
  91.   }
  92. }
  93.  
  94.  
  95. void addlonge (unsigned int v)
  96. {
  97.   _asm
  98.   {
  99.     mov bx,SEG e
  100.     mov es,bx
  101.     mov bx,OFFSET e
  102.     les bx,es:[bx].fraction
  103.     mov ax,v
  104.     add es:[bx],ax
  105.     jnc ade_2
  106.   ade_1:
  107.     inc bx
  108.     inc bx
  109.     inc word ptr es:[bx]
  110.     jz  ade_1
  111.   ade_2:
  112.   }
  113. }
  114.  
  115.  
  116. unsigned int emalzt (void)
  117. {
  118.   _asm
  119.   {
  120.     mov ax,SEG digits
  121.     mov es,ax
  122.     mov di,es:[digits]
  123.     mov bx,SEG e
  124.     mov es,bx
  125.     mov bx,OFFSET e
  126.     mov si,es:[bx].valid
  127.     les bx,es:[bx].fraction
  128.     shl si,1
  129.     shl di,1
  130.  
  131.     mov ax,10000
  132.     mul word ptr es:[bx+si]
  133.     and ax,ax
  134.     jnz mul_2
  135.     mov bx,SEG e
  136.     mov es,bx
  137.     mov bx,OFFSET e
  138.     inc word ptr es:[bx].valid
  139.     les bx,es:[bx].fraction
  140.     jmp mul_2
  141.  
  142.   mul_1:
  143.     mov ax,10000
  144.     mul word ptr es:[bx+si]
  145.     add ax,cx
  146.     jnc mul_2
  147.     inc dx
  148.   mul_2:
  149.     mov cx,dx
  150.     mov es:[bx+si],ax
  151.     inc si
  152.     inc si
  153.     cmp si,di
  154.     jc  mul_1
  155.     mov ax,cx
  156.   }
  157. }
  158.  
  159. int main (int argc, char **argv)
  160. {
  161.   unsigned int d;
  162.   long decimals = 0;
  163.   char c, last[5], *plast = last;
  164.  
  165.   if (argc > 2)
  166.   {
  167.     fprintf (stderr, "too many arguments (only one number allowed)!\n");
  168.     exit (1);
  169.   }
  170.   if (argc == 2)
  171.   {
  172.     while (c = *argv[1]++)
  173.     {
  174.       if (isdigit (c))
  175.       {
  176.         decimals = decimals * 10 + c - '0';
  177.         if (decimals > DECMAX)
  178.         {
  179.           fprintf (stderr, "number of decimals may not exceed %ld!\n",
  180.                    DECMAX);
  181.           exit (2);
  182.         }
  183.       }
  184.       else
  185.       {
  186.         fprintf (stderr, "parameter must be the number of decimals!\n");
  187.         exit (3);
  188.       }
  189.     }
  190.   }
  191.   else decimals = DECMAX;
  192.  
  193.   digits = (decimals * 49 + 354) / 236;
  194.  
  195.   if (e.fraction = malloc (digits * sizeof (unsigned int)))
  196.   {
  197.     if (f.fraction = malloc (digits * sizeof (unsigned int)))
  198.     {
  199.       d = digits - 1;
  200.       e.fraction[d] = 0xAAAA;       /* e =  2/3 */
  201.       f.fraction[d] = 0x0AAA;       /* f = 1/24 */
  202.       while (d--)
  203.       {
  204.         e.fraction[d] = 0xAAAA;
  205.         f.fraction[d] = 0xAAAA;
  206.       }
  207.       f.valid = digits;
  208.       d = 4;
  209.  
  210.       while (f.valid)
  211.       {
  212.         fprintf (stderr, "  f: %5d\r", f.valid);
  213.         addlongf ();
  214.         divlongf (++d);
  215.       }
  216.       addlonge ((d >> 1) - 1);  /* compensate for rounding errors */
  217.  
  218.       printf ("2.");
  219.       e.valid = 0;
  220.  
  221.       while (decimals > 4)
  222.       {
  223.         printf ("%04d", emalzt ());
  224.         decimals -= 4;
  225.       }
  226.  
  227.       sprintf (plast, "%04d", emalzt ());
  228.       while (decimals--)
  229.         fputchar (*plast++);
  230.  
  231.       free (f.fraction);
  232.     }
  233.     free (e.fraction);
  234.   }
  235.   exit (0);
  236. }
  237.