home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / picalcu.zip / ESCOTT.C next >
Text File  |  2002-09-20  |  7KB  |  330 lines
  1. /* 16 bit version of the  Escott algorithm for Pi */
  2.  
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include <malloc.h>
  6. #include <ctype.h>
  7.  
  8. #define DECMAX 94835L
  9.  
  10. typedef struct
  11. { unsigned int *fraction;
  12.   unsigned int valid;
  13. } LONGFRACT;
  14.  
  15. LONGFRACT pi;
  16. unsigned int digits;
  17.  
  18. void divlongf (LONGFRACT *dividend, unsigned int divisor, unsigned int rest)
  19. {
  20.   _asm
  21.   {
  22.     les bx,dividend
  23.     mov si,es:[bx].valid
  24.     les bx,es:[bx].fraction
  25.     mov cx,divisor
  26.     mov dx,rest
  27.  
  28.     dec si
  29.     shl si,1
  30.     mov ax,es:[bx+si]
  31.     div cx
  32.     and ax,ax
  33.     jnz div_2
  34.     les bx,dividend
  35.     dec word ptr es:[bx].valid
  36.     les bx,es:[bx].fraction
  37.     jmp div_2
  38.  
  39.   div_1:
  40.     mov ax,es:[bx+si]
  41.     div cx
  42.   div_2:
  43.     mov es:[bx+si],ax
  44.     sub si,2
  45.     jnc div_1
  46.   }
  47. }
  48.  
  49.  
  50. void addlongf (LONGFRACT *bezug, LONGFRACT *zahl)
  51. {
  52.   _asm
  53.   {
  54.     push ds
  55.     push bp
  56.     les bx,zahl
  57.     mov cx,es:[bx].valid
  58.     les bx,es:[bx].fraction
  59.     lds bp,bezug
  60.     mov di,ds:[bp].valid
  61.     lds bp,ds:[bp].fraction
  62.     xor si,si
  63.     shl di,1
  64.     clc
  65.  
  66.   add_1:
  67.     mov ax,es:[bx+si]
  68.     adc ds:[bp+si],ax
  69.     inc si
  70.     inc si
  71.     dec cx
  72.     jnz add_1
  73.  
  74.     jnc add_4
  75.   add_2:
  76.     cmp si,di
  77.     jc add_3
  78.     lds bp,bezug
  79.     inc word ptr ds:[bp].valid
  80.     lds bp,ds:[bp].fraction
  81.     inc di
  82.     inc di
  83.   add_3:
  84.     inc word ptr ds:[bp+si]
  85.     jnz add_4
  86.     inc si
  87.     inc si
  88.     jmp add_2
  89.   add_4:
  90.     pop bp
  91.     pop ds
  92.   }
  93. }
  94.  
  95.  
  96. void sublongf (LONGFRACT *bezug, LONGFRACT *zahl)
  97. {
  98.   _asm
  99.   {
  100.     push ds
  101.     push bp
  102.     les bx,zahl
  103.     mov cx,es:[bx].valid
  104.     les bx,es:[bx].fraction
  105.     lds bp,bezug
  106.     lds bp,ds:[bp].fraction
  107.     xor si,si
  108.     clc
  109.  
  110.   sub_1:
  111.     mov ax,es:[bx+si]
  112.     sbb ds:[bp+si],ax
  113.     inc si
  114.     inc si
  115.     dec cx
  116.     jnz sub_1
  117.  
  118.     jnc sub_3
  119.   sub_2:
  120.     sub word ptr ds:[bp+si],1
  121.     jnc sub_3
  122.     inc si
  123.     inc si
  124.     jmp sub_2
  125.   sub_3:
  126.     pop bp
  127.     pop ds
  128.   }
  129. }
  130.  
  131.  
  132. unsigned int pimalzt (void)
  133. {
  134.   _asm
  135.   {
  136.     mov ax,SEG digits
  137.     mov es,ax
  138.     mov di,es:[digits]
  139.     mov bx,SEG pi
  140.     mov es,bx
  141.     mov bx,OFFSET pi
  142.     mov si,es:[bx].valid
  143.     les bx,es:[bx].fraction
  144.     shl si,1
  145.     shl di,1
  146.  
  147.     mov ax,10000
  148.     mul word ptr es:[bx+si]
  149.     and ax,ax
  150.     jnz mul_2
  151.     mov bx,SEG pi
  152.     mov es,bx
  153.     mov bx,OFFSET pi
  154.     inc word ptr es:[bx].valid
  155.     les bx,es:[bx].fraction
  156.     jmp mul_2
  157.  
  158.   mul_1:
  159.     mov ax,10000
  160.     mul word ptr es:[bx+si]
  161.     add ax,cx
  162.     jnc mul_2
  163.     inc dx
  164.   mul_2:
  165.     mov cx,dx
  166.     mov es:[bx+si],ax
  167.     inc si
  168.     inc si
  169.     cmp si,di
  170.     jc  mul_1
  171.     mov ax,cx
  172.   }
  173. }
  174.  
  175. void cpylongf (LONGFRACT *ziel, LONGFRACT *quelle)
  176. {
  177.   register unsigned int *z, *q, i;
  178.   ziel->valid = i = quelle->valid;
  179.   z = ziel->fraction;
  180.   q = quelle->fraction;
  181.   while (i--) *z++ = *q++;
  182. }
  183.  
  184.  
  185. int main (int argc, char *argv[])
  186. {
  187.   LONGFRACT a,b,c,d,s;
  188.   int minus = 0;
  189.   unsigned int denomin = 1, a_high;
  190.   long decimals = 0;
  191.   char ch, last[5], *next = last;
  192.  
  193.   if (argc > 2)
  194.   {
  195.     fprintf (stderr, "too many arguments (only one number allowed)!\n");
  196.     exit (1);
  197.   }
  198.   if (argc == 2)
  199.   {
  200.     while (ch = *argv[1]++)
  201.     {
  202.       if (isdigit (ch))
  203.       {
  204.         ch -= '0';
  205.         decimals = decimals * 10 + ch;
  206.         if (decimals > DECMAX)
  207.         {
  208.           fprintf (stderr, "number of decimals may not exceed %ld!\n",
  209.                    DECMAX);
  210.           exit (2);
  211.         }
  212.       }
  213.       else
  214.       {
  215.         fprintf (stderr, "parameter must be the number of decimals!\n");
  216.         exit (3);
  217.       }
  218.     }
  219.   }
  220.   else decimals = DECMAX;
  221.  
  222.   digits = (decimals * 49 + 354) / 236;
  223.  
  224.   if (a.fraction = calloc (digits, sizeof (unsigned int)))
  225.   {
  226.     if (b.fraction = calloc (digits, sizeof (unsigned int)))
  227.     {
  228.       if (c.fraction = calloc (digits, sizeof (unsigned int)))
  229.       {
  230.         if (d.fraction = calloc (digits, sizeof (unsigned int)))
  231.         {
  232.           if (s.fraction = calloc (digits, sizeof (unsigned int)))
  233.           {
  234.             if (pi.fraction = calloc (digits, sizeof (unsigned int)))
  235.             {
  236.               a.valid = b.valid = c.valid = d.valid = digits;
  237.  
  238.               a_high = 3;
  239.               divlongf (&a, 28, 4);
  240.               divlongf (&b, 443, 8);
  241.               divlongf (&c, 1393, 20);
  242.               divlongf (&d, 11018, 40);
  243.  
  244.               cpylongf (&pi, &a);
  245.               addlongf (&pi, &b);
  246.               sublongf (&pi, &c);
  247.               sublongf (&pi, &d);
  248.  
  249.               while (a.valid)
  250.               {
  251.                 divlongf (&a, 784, a_high);
  252.                 a_high = 0;
  253.                 fprintf (stderr, "  1: %5d", a.valid);
  254.                 if (a.valid)
  255.                 {
  256.                   cpylongf (&s, &a);
  257.                   if (b.valid)
  258.                   {
  259.                     divlongf (&b, 443, 0);
  260.                     if (b.valid)
  261.                       divlongf (&b, 443, 0);
  262.                     fprintf (stderr, "  2: %5d", b.valid);
  263.                     if (b.valid)
  264.                     {
  265.                       addlongf (&s, &b);
  266.                       if (c.valid)
  267.                       {
  268.                         divlongf (&c, 1393, 0);
  269.                         if (c.valid)
  270.                           divlongf (&c, 1393, 0);
  271.                         fprintf (stderr, "  3: %5d", c.valid);
  272.                         if (c.valid)
  273.                         {
  274.                           sublongf (&s, &c);
  275.                           if (d.valid)
  276.                           {
  277.                             divlongf (&d, 11018, 0);
  278.                             if (d.valid)
  279.                               divlongf (&d, 11018, 0);
  280.                             fprintf (stderr, "  4: %5d", d.valid);
  281.                             if (d.valid)
  282.                               sublongf (&s, &d);
  283.                           } /* d.valid */
  284.                         } /* c.valid 2 */
  285.                       } /* c.valid 1 */
  286.                     } /* b.valid 2 */
  287.                   } /* b.valid 1 */
  288.  
  289.                   divlongf (&s, denomin += 2, 0);
  290.                   if (s.valid)
  291.                   { if (minus = ! minus)
  292.                       sublongf (&pi, &s);
  293.                     else
  294.                       addlongf (&pi, &s);
  295.                   }
  296.                   else a.valid = 0;
  297.  
  298.                 } /* a.valid 2 */
  299.  
  300.                 fprintf (stderr, "\r");
  301.               } /* a.valid 1 */
  302.  
  303.               printf ("3.");
  304.               pi.valid = 0;
  305.  
  306.               while (decimals > 4)
  307.               {
  308.                 printf ("%04d", pimalzt ());
  309.                 decimals -= 4;
  310.               }
  311.  
  312.               sprintf (last, "%04d", pimalzt ());
  313.               while (decimals--)
  314.                 fputchar (*next++);
  315.  
  316.               free (pi.fraction);
  317.             }
  318.             free (s.fraction);
  319.           }
  320.           free (d.fraction);
  321.         }
  322.         free (c.fraction);
  323.       }
  324.       free (b.fraction);
  325.     }
  326.     free (a.fraction);
  327.   }
  328.   exit (0);
  329. }
  330.