home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
OS/2 Shareware BBS: Science
/
Science.zip
/
picalcu.zip
/
EULER.C
< prev
next >
Wrap
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);
}