OS/2 Shareware BBS: Science/Science.zip/picalcu.zip/EULER.C

Text File | 2002-09-20 | 4KB | 237 lines

/* Determining of the  Euler constant e - 16 bit version */

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <ctype.h>

#define DECMAX 157761L
/* 287188 would be allowed for the  WORD format of the divisor */

typedef struct
{ unsigned int *fraction;
  unsigned int valid;
} LONGFRACT;

LONGFRACT e, f;
unsigned int digits;

void divlongf (unsigned int d)
{
  _asm
  {
    mov bx,SEG f
    mov es,bx
    mov bx,OFFSET f
    mov si,es:[bx].valid
    les bx,es:[bx].fraction
    mov cx,d
    xor dx,dx

    dec si
    shl si,1
    mov ax,es:[bx+si]
    div cx
    and ax,ax
    jnz div_2
    mov bx,SEG f
    mov es,bx
    mov bx,OFFSET f
    dec word ptr es:[bx].valid
    les bx,es:[bx].fraction
    jmp div_2

  div_1:
    mov ax,es:[bx+si]
    div cx
  div_2:
    mov es:[bx+si],ax
    sub si,2
    jnc div_1
  }
}


void addlongf ()
{
  _asm
  {
    push ds
    push bp
    mov bx,SEG e
    mov es,bx
    mov bx,OFFSET e
    les bx,es:[bx].fraction
    mov bp,SEG f
    mov ds,bp
    mov bp,OFFSET f
    mov cx,ds:[bp].valid
    lds bp,ds:[bp].fraction
    xor si,si
    clc

  add_1:
    mov ax,ds:[bp+si]
    adc es:[bx+si],ax
    inc si
    inc si
    dec cx
    jnz add_1

    jnc add_3
  add_2:
    inc word ptr es:[bx+si]
    jnz add_3
    inc si
    inc si
    jmp add_2
  add_3:
    pop bp
    pop ds
  }
}


void addlonge (unsigned int v)
{
  _asm
  {
    mov bx,SEG e
    mov es,bx
    mov bx,OFFSET e
    les bx,es:[bx].fraction
    mov ax,v
    add es:[bx],ax
    jnc ade_2
  ade_1:
    inc bx
    inc bx
    inc word ptr es:[bx]
    jz  ade_1
  ade_2:
  }
}


unsigned int emalzt (void)
{
  _asm
  {
    mov ax,SEG digits
    mov es,ax
    mov di,es:[digits]
    mov bx,SEG e
    mov es,bx
    mov bx,OFFSET e
    mov si,es:[bx].valid
    les bx,es:[bx].fraction
    shl si,1
    shl di,1

    mov ax,10000
    mul word ptr es:[bx+si]
    and ax,ax
    jnz mul_2
    mov bx,SEG e
    mov es,bx
    mov bx,OFFSET e
    inc word ptr es:[bx].valid
    les bx,es:[bx].fraction
    jmp mul_2

  mul_1:
    mov ax,10000
    mul word ptr es:[bx+si]
    add ax,cx
    jnc mul_2
    inc dx
  mul_2:
    mov cx,dx
    mov es:[bx+si],ax
    inc si
    inc si
    cmp si,di
    jc  mul_1
    mov ax,cx
  }
}

int main (int argc, char **argv)
{
  unsigned int d;
  long decimals = 0;
  char c, last[5], *plast = last;

  if (argc > 2)
  {
    fprintf (stderr, "too many arguments (only one number allowed)!\n");
    exit (1);
  }
  if (argc == 2)
  {
    while (c = *argv[1]++)
    {
      if (isdigit (c))
      {
        decimals = decimals * 10 + c - '0';
        if (decimals > DECMAX)
        {
          fprintf (stderr, "number of decimals may not exceed %ld!\n",
                   DECMAX);
          exit (2);
        }
      }
      else
      {
        fprintf (stderr, "parameter must be the number of decimals!\n");
        exit (3);
      }
    }
  }
  else decimals = DECMAX;

  digits = (decimals * 49 + 354) / 236;

  if (e.fraction = malloc (digits * sizeof (unsigned int)))
  {
    if (f.fraction = malloc (digits * sizeof (unsigned int)))
    {
      d = digits - 1;
      e.fraction[d] = 0xAAAA;       /* e =  2/3 */
      f.fraction[d] = 0x0AAA;       /* f = 1/24 */
      while (d--)
      {
        e.fraction[d] = 0xAAAA;
        f.fraction[d] = 0xAAAA;
      }
      f.valid = digits;
      d = 4;

      while (f.valid)
      {
        fprintf (stderr, "  f: %5d\r", f.valid);
        addlongf ();
        divlongf (++d);
      }
      addlonge ((d >> 1) - 1);  /* compensate for rounding errors */

      printf ("2.");
      e.valid = 0;

      while (decimals > 4)
      {
        printf ("%04d", emalzt ());
        decimals -= 4;
      }

      sprintf (plast, "%04d", emalzt ());
      while (decimals--)
        fputchar (*plast++);

      free (f.fraction);
    }
    free (e.fraction);
  }
  exit (0);
}