How does vic compute these

[Iltanen]

Does VIC use Taylor series to compute trigonometric functions? My calculator uses them to calculate e^x but does VIC?

[nippur72]

I think a simple polynomial interpolation is used (I guess because it's faster).

I remember that in the ROM Kernal there are the polynomial coefficients for the various functions. I'm sure Mike can give additional details.

[pitcalco]

Does VIC use Taylor series to compute trigonometric functions? My calculator uses them to calculate e^x but does VIC?

I know that I did have to resort to Taylor series to get inverse functions, until I discovered the ATN function - silly me

As to what the VIC does inside its own guts, I am inclined to think that for the kind of precision the VIC has, the Taylor series would make most sense.

Who knows?...oh, right, perhaps Mike does

[nippur72]

well, now that you triggered my curiosity I've looked in the ROM kernal and found that:

LOG() is computed with a polynomial of degree 4
EXP() is computed with a polynomial of degree 8
SIN(), COS(), TAN() are computed with a polynomial of degree 6
ATN() is computed with a polynomial of degree 12

all the others are a combination of the above functions.

Talking of Mike, he has written a nice C routine to convert a floating point number into the CBM format (5 bytes). Perhaps he can write the inverse routine so to get the actual coefficients, taking them directly from the ROM. (In that case we could evaluate the approximation error).

[Mike]

I found time for this yesterday evening, so here we go:

Code: Select all

/*
 *  cbm_sine.c
 *
 *  written 2008-08-21 by Michael Kircher
 *
 *  Analysis of Chebyshev approximation of SIN(x) in CBM BASIC
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define PI 3.141592653589793

double horner(double x, double *coef, unsigned int n)
{
  double result=coef[n];
  while(n>0) result=result*x+coef[--n];
  return(result);
}

double cbm_sine(double x)
{
  double u;
  double coef[6]=
  {/* Coefficients for interpolation polynomial, as given by cbm2dbl. */
     6.283185307,
   -41.3417021,
    81.60522369,
   -76.70417026,
    42.00779712,
   -14.38139067,
  };
  
  /* scale argument, and reduce range to [0,1[ */
  u=x/(2*PI);
  u=u-floor(u);

  /* 
   *  Employ symmetries to bring the argument into the range [-0.25,0.25].
   *  A similar, but more complicated routine also is used in the BASIC ROM.
   */
  if(u>0.75) u=u-1.0;
  if(u>0.25) u=0.5-u;

  /*
   *  Finally, do "series" evaluation in u², and
   *  multiply with u to get odd powers of u.
   */
  return(u*horner(u*u,coef,5));
}

/*
 *  Aside: COS(x)=SIN(x+pi/2), TAN(x)=SIN(x)/COS(x)
 */

double cbm2dbl(unsigned char cbm_exponent, unsigned long cbm_mantissa)
{
  double result;
  if(0 == cbm_exponent && 0 == cbm_mantissa) return(0.0);
  result=ldexp(0.5+(cbm_mantissa & 0x7FFFFFFF)/4294967296.0,cbm_exponent-128);
  return((cbm_mantissa & 0x80000000)?-result:result);
}

int main(int argc, char *argv[])
{
  int i;
  double x;
  double cbm_result,ieee_result;
  double abs_diff;

  /*
   *  Convert the coefficients used for the Chebyshev approximation
   *  of the sine function, from CBM float to IEEE double. Only the
   *  first 10 significant digits are shown, as CBM float doesn't
   *  provide more precision anyway.
   */
  printf("%.10G\n",cbm2dbl(0x83,0x490FDAA2));
  printf("%.10G\n",cbm2dbl(0x86,0xA55DE728));
  printf("%.10G\n",cbm2dbl(0x87,0x2335DFE1));
  printf("%.10G\n",cbm2dbl(0x87,0x99688901));
  printf("%.10G\n",cbm2dbl(0x86,0x2807FBF8));
  printf("%.10G\n",cbm2dbl(0x84,0xE61A2D1B));

  /* first test */
  printf("%.9G\n",cbm_sine(0.25*PI)); /* = 0.707106781, O.K. */

  /*
   *  Check accuracy against IEEE sin(x) function,
   *  in full double precision, for some arguments
   *  in the range 0 to pi/2:
   */
  abs_diff=0.0;
  for(i=0; i<=1000; i++)
  {
    x=PI/2.0*(i/1000.0);
     cbm_result=cbm_sine(x);
    ieee_result=    sin (x);
    printf("x=%f,\tSIN(x)=%.15E,\tsin(x)=%.15E\n",x,cbm_result,ieee_result);
    if(abs_diff<fabs(cbm_result-ieee_result)) abs_diff=fabs(cbm_result-ieee_result);
  }

  printf("Maximum absolute difference: %E\n",abs_diff); /* =8.150047E-011,
  which is quite good for a 5th degree polynomial. CBM BASIC however must do
  the same calculation with CBM floats, which results in additional rounding
  errors. Thus the last 1 or 2 bits of the mantissa may not be correct, but
  this is still accurate to 9 decimal digits. */

  exit(EXIT_SUCCESS);
}

[nippur72]

great Mike, I'm amazed!

And thank you also for the cbm2dbl(), which now completes my set of tools for the CBM floating point format.

I still have to understand the tricks done before the horner() function call (e.g. u^2), they go beyond my comphrension

[Iltanen]

Mike, that was impressive.