OT: Calculating advanced math functions?

  • Follow


Sorry about the OT question, but this is for a C64 demo.  I need to
know how to calculate SIN, SQR, EXP, LOGx and any other functions that
would be useful for graphics.  The C64 is only an 8-bit computer and
would have a hard time calculating these with any speed, but I can
calculate tables at run-time.  My knowledge is limited to a high-
school education; otherwise I probably wouldn't be asking.  :(
0
Reply rose.joseph12 (88) 3/27/2012 1:29:45 PM

On Mar 27, 7:29=A0am, Harry Potter <rose.josep...@yahoo.com> wrote:
> Sorry about the OT question, but this is for a C64 demo. =A0I need to
> know how to calculate SIN, SQR, EXP, LOGx and any other functions that
> would be useful for graphics. =A0The C64 is only an 8-bit computer and
> would have a hard time calculating these with any speed, but I can
> calculate tables at run-time. =A0My knowledge is limited to a high-
> school education; otherwise I probably wouldn't be asking. =A0:(

I hope you don't get intimidated easily.  Trignometric or circular
function evaluation is covered in the first year of calculus.  The
input angles for the series approximations of trig functions need to
be in radian measure.  Multiply your angle in degrees by pi/180 for
the angle in radian measure. (pi=3D3.14159265...)  Infinite series
identities can be found on the lower right hand corner of page 74 of
the US government publication HANDBOOK of MATHEMATICAL FUNCTIONS with
FORMULAS, GRAPHS, and MATHEMATICAL TABLES edited by Milton Abramowitz
and Irene Stegun.  If you can't find it in your public library, try
the engineering library at a nearby university.

Unfortunately, it takes a long time for Commodores to evaluate an
infinite series.  The algorythm usually starts out evaluating a
polynomial of the first N terms of the series.  N is determined by the
magtitude of a remainder function.  The smaller the remainder, the
better the polynomial approximation, but the larger the value of N.
The above algorythm, requires larger values of N for larger angles,
You may want to exploit the periodicity of trig functions so that you
input angles of less than 2*pi.  Trig function periodicities are of
2*pi in length.  I would select a remainder term for sin and cosin
functions and the angle x as follows

     R(x,N)=3D[(x)^(N=3D1)]/[(N+1)!]

For other function evaluation,  just the series expressions  and
remainder terms change.  Only trig functions may require conversions
from degrees to radians.
0
Reply r_u_sure (108) 3/27/2012 7:35:33 PM


27.03.2012 16:29, Harry Potter kirjoitti:
> Sorry about the OT question, but this is for a C64 demo.  I need to
> know how to calculate SIN, SQR, EXP, LOGx and any other functions that
> would be useful for graphics.  The C64 is only an 8-bit computer and
> would have a hard time calculating these with any speed, but I can
> calculate tables at run-time.  My knowledge is limited to a high-
> school education; otherwise I probably wouldn't be asking.  :(

In demos pre-calculated tables and integers are used mostly. They are
fastest.

You can calculate pre-defined tables easily in Basic. You do not need to
have /exact/ math to draw dots in low resolution. You only need to have
good enough approximates.

Pre-calc tables need not to be huge. If 256 bytes are enough, the
pre-calc table is wery fast to access. Just a lda $table,x will do just
fine. 4 cycles.

If you want more precision, 16 bits should be enough. It is slower to
access, but WAY faster than calculating. I think that 512 16 bits values
is enough, and it takes only one kilobyte in memory. Learn
self-modifying code.


0
Reply pekka.NOtakala (8) 3/27/2012 9:55:25 PM

I was planning to calculate some tables at start-up using longs--my C compi=
ler handles longs but not floats--and use the tables for calculations at ru=
n-time.  My demos are simple--they look like a child did them.  However, th=
ey *do* demonstrate the usage of another program.

On Tuesday, March 27, 2012 5:55:25 PM UTC-4, Pekka Takala wrote:
> 27.03.2012 16:29, Harry Potter kirjoitti:
> > Sorry about the OT question, but this is for a C64 demo.  I need to
> > know how to calculate SIN, SQR, EXP, LOGx and any other functions that
> > would be useful for graphics.  The C64 is only an 8-bit computer and
> > would have a hard time calculating these with any speed, but I can
> > calculate tables at run-time.  My knowledge is limited to a high-
> > school education; otherwise I probably wouldn't be asking.  :(
>=20
> In demos pre-calculated tables and integers are used mostly. They are
> fastest.
>=20
> You can calculate pre-defined tables easily in Basic. You do not need to
> have /exact/ math to draw dots in low resolution. You only need to have
> good enough approximates.
>=20
> Pre-calc tables need not to be huge. If 256 bytes are enough, the
> pre-calc table is wery fast to access. Just a lda $table,x will do just
> fine. 4 cycles.
>=20
> If you want more precision, 16 bits should be enough. It is slower to
> access, but WAY faster than calculating. I think that 512 16 bits values
> is enough, and it takes only one kilobyte in memory. Learn
> self-modifying code.

0
Reply rose.joseph12 (88) 3/30/2012 1:10:36 PM

On Mar 27, 7:29=A0am, Harry Potter <rose.josep...@yahoo.com> wrote:
> Sorry about the OT question, but this is for a C64 demo. =A0I need to
> know how to calculate SIN, SQR, EXP, LOGx and any other functions that
> would be useful for graphics. =A0The C64 is only an 8-bit computer and
> would have a hard time calculating these with any speed, but I can
> calculate tables at run-time. =A0My knowledge is limited to a high-
> school education; otherwise I probably wouldn't be asking. =A0:(

Please pardon my delay.  My references were incomplete so I had to
improvise.  This took longer than I would have liked.

I selected 3 candidates for approximating log(z) with a series of n
terms.  For z=3D0 or negative, log(z) doesn't exist or is a complex
numerical value.  My code doesn't check for these invalid values of
z.  For the n terms in each series, a recursion relationship exists
for terms m-1 and m.  That is, term m can be determined in a simple
way from term m-1 for m=3D1 to n.  An application of a modified
geometric series theory shows, that, discarding sign, the product of
the first term not included in the sum and the infinite geometric
total is lager than the infnite term total not included for
approximating log(z).

Candidate 1 begins in line 100 of the Commodore BASIC program listed
below. It impliments a series of terms in the form of ((-1)^(m
+1)*(z-1)^(m))/m  Where m takes on values of 1 to n
For z greater than or =3D to 2, -log(1/z) is approximated.
Candidate 2 begins at line 200 and impliments the sum of terms in the
form of ((z-1)/z)^(m))/m for m=3D 1 to n.  -log(1/z) is approximted for
z =3D<1/2.
Finally, Candidate 3 begins in 300.  The sum of terms on the form of
2*(((z_-1)/(z+1))^(2*m-1))/(2*m-1) for m=3D1 to n.  This last series
converges the fastest for the values of z in the data statements.

10 forl=3D1to5:readz:gosub100:next
20 restore:forl=3D1to5:readz:gosub200:next
30 restore:forl=3D1to5:readz:gosub300:next
90 stop
100 t=3Dz-1:r=3D-t: if z>=3D2thent=3D1-1/z:r=3Dt
110 g=3D1/(1-abs(r)):n=3D1:ln=3D0:printz
120 ln=3Dln+t:t=3Dt*n*r/(n+1)
130 u=3Dabs(t*g):println;t;u;n:n=3Dn+1
140 ifu>0.0001then120
150 stop:return
200 t=3D(z-1)/z:if z=3D<0.5thent=3Dz-1
210 r=3Dabs(t):g=3D1/(1-abs(r)):n=3D1:ln=3D0:printz
220 ln=3Dln+t:t=3Dt*n*r/(n+1)
230 u=3Dabs(t*g):println;t;u;n:n=3Dn+1
240 ifu>0.0001then220
250 stop:return
300 t=3D(z-1)/(z+1):r=3Dt*t:t=3Dt+t
310 g=3D1/(1-abs(r)):n=3D-1:ln=3D0:printz
320 ln=3Dln+t:n=3Dn+2:t=3Dt*n*r/(n+2)
330 u=3Dabs(t*g):println;t;u;(n+1)/2
340 ifu>0.0001then320
350 stop:return
400 data4.48168907,1.947734041
500 data .367879441,2,0.5

Treatments for trig and exponenyial will fowllow,
0
Reply r_u_sure (108) 4/2/2012 10:57:34 AM

Thank you.  I saved your post o my flash drive for later use.  :)

On Monday, April 2, 2012 6:57:34 AM UTC-4, rusure wrote:
> On Mar 27, 7:29=A0am, Harry Potter <rose.josep...@yahoo.com> wrote:
> > Sorry about the OT question, but this is for a C64 demo. =A0I need to
> > know how to calculate SIN, SQR, EXP, LOGx and any other functions that
> > would be useful for graphics. =A0The C64 is only an 8-bit computer and
> > would have a hard time calculating these with any speed, but I can
> > calculate tables at run-time. =A0My knowledge is limited to a high-
> > school education; otherwise I probably wouldn't be asking. =A0:(
>=20
> Please pardon my delay.  My references were incomplete so I had to
> improvise.  This took longer than I would have liked.
>=20
> I selected 3 candidates for approximating log(z) with a series of n
> terms.  For z=3D0 or negative, log(z) doesn't exist or is a complex
> numerical value.  My code doesn't check for these invalid values of
> z.  For the n terms in each series, a recursion relationship exists
> for terms m-1 and m.  That is, term m can be determined in a simple
> way from term m-1 for m=3D1 to n.  An application of a modified
> geometric series theory shows, that, discarding sign, the product of
> the first term not included in the sum and the infinite geometric
> total is lager than the infnite term total not included for
> approximating log(z).
>=20
> Candidate 1 begins in line 100 of the Commodore BASIC program listed
> below. It impliments a series of terms in the form of ((-1)^(m
> +1)*(z-1)^(m))/m  Where m takes on values of 1 to n
> For z greater than or =3D to 2, -log(1/z) is approximated.
> Candidate 2 begins at line 200 and impliments the sum of terms in the
> form of ((z-1)/z)^(m))/m for m=3D 1 to n.  -log(1/z) is approximted for
> z =3D<1/2.
> Finally, Candidate 3 begins in 300.  The sum of terms on the form of
> 2*(((z_-1)/(z+1))^(2*m-1))/(2*m-1) for m=3D1 to n.  This last series
> converges the fastest for the values of z in the data statements.
>=20
> 10 forl=3D1to5:readz:gosub100:next
> 20 restore:forl=3D1to5:readz:gosub200:next
> 30 restore:forl=3D1to5:readz:gosub300:next
> 90 stop
> 100 t=3Dz-1:r=3D-t: if z>=3D2thent=3D1-1/z:r=3Dt
> 110 g=3D1/(1-abs(r)):n=3D1:ln=3D0:printz
> 120 ln=3Dln+t:t=3Dt*n*r/(n+1)
> 130 u=3Dabs(t*g):println;t;u;n:n=3Dn+1
> 140 ifu>0.0001then120
> 150 stop:return
> 200 t=3D(z-1)/z:if z=3D<0.5thent=3Dz-1
> 210 r=3Dabs(t):g=3D1/(1-abs(r)):n=3D1:ln=3D0:printz
> 220 ln=3Dln+t:t=3Dt*n*r/(n+1)
> 230 u=3Dabs(t*g):println;t;u;n:n=3Dn+1
> 240 ifu>0.0001then220
> 250 stop:return
> 300 t=3D(z-1)/(z+1):r=3Dt*t:t=3Dt+t
> 310 g=3D1/(1-abs(r)):n=3D-1:ln=3D0:printz
> 320 ln=3Dln+t:n=3Dn+2:t=3Dt*n*r/(n+2)
> 330 u=3Dabs(t*g):println;t;u;(n+1)/2
> 340 ifu>0.0001then320
> 350 stop:return
> 400 data4.48168907,1.947734041
> 500 data .367879441,2,0.5
>=20
> Treatments for trig and exponenyial will fowllow,

0
Reply rose.joseph12 (88) 4/3/2012 12:53:20 PM

On Apr 3, 6:53=A0am, Harry Potter <rose.josep...@yahoo.com> wrote:
> Thank you. =A0I saved your post o my flash drive for later use. =A0:)


>
> > Treatments for trig and exponential will follow,

This is the final installment for approximating transcendental
functions by Taylor and
Mclauren polynomial methods in Commodore BASIC.  The sine, cosine, and
exponential computations begin in lines 500, 600, 700 respectively.
The topic was covered in my third semester in Calculus course.  The
polynomial evaluation exploits a recursion relationship between
contiguous terms in the polynomial.  The precise discrepancy between a
function and it's polynomial is difficult or impossible to evaluate.
However, for triganometric and exponential functions, upper bounds of
the absolute discrepancy can be calculated.  The upper bound is
contained in the variable u.  The polynomials are contained in the
variables s, c, and e for sine, cosine, and exponential functions
respectively.  For a given polynomial term t, the product of t and r
results in the next term in the polynomial.

10 forj=3D1to8:readx:z=3Dx*3.14159265/180
20 gosub500:gosub600:nextj
30 forx=3D-10to10:z=3Dx/2:gosub700:nextx
40 stop
500 printx;sin(z):t=3Dz:r=3D-t*t::u=3D1:n=3D-1:s=3D0
510 n=3Dn+2:s=3Ds+t:printn;t;s;:t=3Dt*r/((n+1)*(n+2))
520 u=3Du*r/(n*(n+1))::printu
530 ifabs(u)>0.0001then510
540 stop:return
600 printx;cos(z):t=3D1:r=3D-z*z::u=3Dz:n=3D0:s=3D0
610 n=3Dn+2:s=3Ds+t:printn;t;s;u
620 t=3Dt*r/(n*(n-1)):u=3Du*r/(n*(n+1)):ifabs(u)>0.0001then610
630 stop:return
700 u=3D1:n=3Dabs(z):ifn=3D0then730
710 ifn-int(n)>0 then n=3Dint(n)+1
720 fori=3D1ton:u=3Du*3:next
730 print3^abs(z);u;z;n;:t=3D1:n=3D0:e=3D0
740 e=3De+t:n=3Dn+1:t=3Dt*z/n:u=3Du*z/n
750 if abs(u)>0.0001 then 740
760 printlog(e):stop:return
900 data 30,60,120,150,210,240,300,330


0
Reply r_u_sure (108) 4/5/2012 1:45:21 AM

Correction

The last line in the discussion

For a given polynomial term t, the product of t and r
results in the next term in the polynomial.

Should be replaced by

For a given polynomial term t, the product of t and r and divided
by a function of n results in the next term in the polynomial.
0
Reply r_u_sure (108) 4/5/2012 1:53:05 AM

Please discard my previous 2 posts

This is the final installment for approximating transcendental
functions by Taylor and Mclauren polynomial methods in Commodore
BASIC.  The sine, cosine, and exponential computations begin in
lines 500, 600 and 700 respectively.  The topic was covered in my
third semester in Calculus course.  The polynomial evaluation
exploits a recursion relationship between contiguous terms in the
polynomial.  The precise discrepancy between a function and it's
polynomial is difficult or impossible to evaluate.  However, for
trigonometric and exponential functions, upper bounds of the
absolute discrepancies can be calculated.  The upper bound is
contained in the variable u.  The polynomials are contained in the
variables s, c, and e for sine, cosine, and exponential functions
respectively.  To get the term t to be added in the next iteration,
multiply the current term and r then divide by (n*(n-1)) for the sine
and by (n+2)*(n+1) for the cosine.  The print  statements in the trig
subroutines have been revised.  In addition, the "s"'s with  "c"'s
in the cosine subroutine beginning in line 600

10 forj=1to8:readx:z=x*3.14159265/180
20 gosub500:gosub600:nextj
30 forx=-10to10:z=x/2:gosub700:nextx
40 stop
500 printx;sin(z);:t=z:r=-t*t::u=1:n=-1:s=0
510 n=n+2:s=s+t:t=t*r/((n+1)*(n+2))
520 u=u*r/(n*(n+1))
530 ifabs(u)>0.0001then510
540 prints,,:return
600 printcos(z);:t=1:r=-z*z::u=z:n=0:c=0
610 n=n+2:c=c+t:t=t*r/(n*(n-1))
620 u=u*r/(n*(n+1)):ifabs(u)>0.0001then610
630 printc:stop:return
700 u=1:n=abs(z):ifn=0then730
710 ifn-int(n)>0 then n=int(n)+1
720 fori=1ton:u=u*3:next
730 print3^abs(z);u;z;n;:t=1:n=0:e=0
740 e=e+t:n=n+1:t=t*z/n:u=u*z/n
750 if abs(u)>0.0001 then 740
760 printlog(e):stop:return
900 data 30,60,120,150,210,240,300,330
0
Reply r_u_sure (108) 4/5/2012 9:00:28 PM

On 2012-04-05, rusure <r_u_sure@mybluelight.com> wrote:
> Please discard my previous 2 posts
>
> This is the final installment for approximating transcendental
> functions by Taylor and Mclauren polynomial methods in Commodore
> BASIC.  The sine, cosine, and exponential computations begin in
> lines 500, 600 and 700 respectively.  The topic was covered in my
> third semester in Calculus course.  The polynomial evaluation
> exploits a recursion relationship between contiguous terms in the
> polynomial.  The precise discrepancy between a function and it's
> polynomial is difficult or impossible to evaluate.  However, for
> trigonometric and exponential functions, upper bounds of the
> absolute discrepancies can be calculated.  The upper bound is
> contained in the variable u.  The polynomials are contained in the
> variables s, c, and e for sine, cosine, and exponential functions
> respectively.  To get the term t to be added in the next iteration,
> multiply the current term and r then divide by (n*(n-1)) for the sine
> and by (n+2)*(n+1) for the cosine.  The print  statements in the trig
> subroutines have been revised.  In addition, the "s"'s with  "c"'s
> in the cosine subroutine beginning in line 600
>
> 10 forj=1to8:readx:z=x*3.14159265/180
> 20 gosub500:gosub600:nextj
> 30 forx=-10to10:z=x/2:gosub700:nextx
> 40 stop
> 500 printx;sin(z);:t=z:r=-t*t::u=1:n=-1:s=0
> 510 n=n+2:s=s+t:t=t*r/((n+1)*(n+2))
> 520 u=u*r/(n*(n+1))
> 530 ifabs(u)>0.0001then510
> 540 prints,,:return
> 600 printcos(z);:t=1:r=-z*z::u=z:n=0:c=0
> 610 n=n+2:c=c+t:t=t*r/(n*(n-1))
> 620 u=u*r/(n*(n+1)):ifabs(u)>0.0001then610
> 630 printc:stop:return
> 700 u=1:n=abs(z):ifn=0then730
> 710 ifn-int(n)>0 then n=int(n)+1
> 720 fori=1ton:u=u*3:next
> 730 print3^abs(z);u;z;n;:t=1:n=0:e=0
> 740 e=e+t:n=n+1:t=t*z/n:u=u*z/n
> 750 if abs(u)>0.0001 then 740
> 760 printlog(e):stop:return
> 900 data 30,60,120,150,210,240,300,330


So I'm just curious...why do you want to approximate trig functions, roots,
etc.on a C64? Got a neat project cooking? Sorry if I missed it in an earlier
post.

I remember reading about a monte carlo program to calculate Pi in a Transactor
a while back, but my disability, bone-idleness, kept me from messing around 
with it.
0
Reply ceratoph (5) 4/9/2012 3:00:30 PM

> So I'm just curious...why do you want to approximate trig functions, roots,
> etc.on a C64? Got a neat project cooking? Sorry if I missed it in an earlier
> post.
>
> I remember reading about a monte carlo program to calculate Pi in a Transactor
> a while back, but my disability, bone-idleness, kept me from messing around
> with it.

Harry Potter posted the first message in this thread.  He wanted to
construct tables of transcendental functions for graphical purposes.
A complete thread can be found here:

http://groups.google.com/group/comp.sys.cbm/browse_frm/thread/84c812ab68f8f90b?hl=en#
0
Reply r_u_sure (108) 4/10/2012 3:59:39 PM

On Monday, April 9, 2012 11:00:30 AM UTC-4, Michael Grossman wrote:
> So I'm just curious...why do you want to approximate trig functions, roots,
> etc.on a C64? Got a neat project cooking? Sorry if I missed it in an earlier
> post.
> 

In truth, I don't remember. :(  I think I wanted to do file compression and graphics.
0
Reply rose.joseph12 (88) 4/25/2012 1:34:43 PM

11 Replies
51 Views

(page loaded in 0.457 seconds)

Similiar Articles:


















7/18/2012 1:39:53 AM


Reply: