Showing posts with label computer science. Show all posts
Showing posts with label computer science. Show all posts

Sunday, February 18, 2024

Casio fx-CG 50 CORDIC Simulation: Approximating Sine and Cosine of Angles

Casio fx-CG 50 CORDIC Simulation:  Approximating Sine and Cosine of Angles




Introduction - How Computers Calculate Mathematical Functions



First developed by Jack E. Volder, the Coordinate Rotation Digital Computer, better known as CORDIC, is an algorithm used to calculate many mathematical functions, including trigonometric functions, logarithms, exponentials, and hyperbolic functions.     CORDIC is a fundamental algorithm, which variants of CORDIC are used in computers and calculators.


Today's focus will be a calculating sines and cosines of angles.  The steps will be detailed in the next sections.  



CORDIC1:  Series of Arctangents


Let Θ be the angle.   The first step is to build Θ as additions and subtractions of terms of arctan(1 ÷ (2^n)), starting at n = 0 and stopping at the required accuracy.  Let A be the approximation.  


Examples: 


Θ = 45° requires only one term:

A = 45° = arctan 1


Θ ≈ 71.56505118°

 A = 71.56505118° = arctan 1 + arctan 1/2


Θ ≈ 32.47119229°

A = 32.47119229° = arctan 1 - arctan 1/2 + arctan 1/4


A series to recreate Θ = 30° takes 32 terms of ±arctan (1 ÷ 2^n) from n = 0 to n = 31.   (8 decimal places)  This is known as coordinate rotation.


The program CORDIC1 determines the number of terms need to created to obtain Θ.   Accuracy in this code is set to 5 decimal places, but can be adjusted (see the line in blue).  


Casio fx-CG 50 Program Code:  CORDIC1


Deg

"ANGLE"?->Θ

0->I

0->A


Lbl 0

tan^-1 (2^(-I))->T

If A<Θ

Then 

A+T->A

Else 

A-T->A

IfEnd


I+1->I


Abs (A-Θ)>1×10^-5=>Goto 0


ClrText

Blue Locate 1,3,"Θ: "

Blue Locate 5,3,Θ


Red Locate 1,4,"A: "

Red Locate 5,4,A


Green Locate 1,7,I


Notes:

-> is the store arrow →

=> is the jump command ⇒


For reference:


arctan 1 = 45°

arctan 1/2 ≈ 26.56505118°

arctan 1/4 ≈ 14.03624347°

arctan 1/8 ≈ 7.125016349°

arctan 1/16 ≈ 3.576334375°

arctan 1/32 ≈ 1.789910608°


CORDIC2:  Calculating Sine and Cosine


Imagine the coordinate (cos Θ, sin Θ) on a unit circle.  A unit circle is a circle with radius of length 1 and center located at the origin (0,0).   


Set the initial angle at 0°.  Then x = cos 0° = 1 and y = sin 0° = 0.   The initial vector is set to be [ [ cos 0° ] [ sin 0° ] ] = [ [ 1 ] [ 0 ] ].


Approximate Θ in terms of arctangent (1 ÷ (2^n)), starting at n = 0.  



Direction of Rotation


Set σ as the direction of rotation.    


A rotation is positive if we add arctan(1 ÷ (2^n)) to A.    For a positive rotation, set σ_i = +1.   


A rotation is negative if we subtract arctan(1 ÷ (2^n)) to A.    For a negative rotation, set σ_i = -1.   


For example:


 Θ ≈ 32.47119229°

A = 32.47119229° = arctan 1 - arctan 1/2 + arctan 1/4


Then:  σ_0 = 1 (positive rotation), σ_1 = -1 (negative rotation), σ_2 =  1 (positive rotation).  



Multiplying Factor


The number of iterations is also used to determine the required multiplication factor:


n = Π( 1 ÷ √(1 + 2^(-2 × I)), I = 0 to I = terms needed - 1)


For the example:


 Θ ≈ 32.47119229°

Needed 3 iterations, i_0 to i_2.  


Then:

n = 1 ÷ √(1 + 2^(-2 × 0)) × 1 ÷ √(1 + 2^(-2 × 1)) × 1 ÷ √(1 + 2^(-2 × 2))

= 4 × √170 ÷ 85

≈ 0.6135719911



Calculating the Next Iteration


The next iteration for x and y are:


x_i+1 = x_i - 2^(-i) × σ_i × y_i


y_i+1 = 2^(-i) × σ_i × x_i + y_i


When A is sufficiently near or equal to Θ, the cosine and sine are approximated as:


cos(A) ≈ n × x_final


sin(A) ≈ n × y_final


For the example:

 Θ ≈ 32.47119229°


The approximated angle, in this case, happens to be exact angle, A = Θ.  And,

cos(A) ≈ 0.8436614877 

sin(A) ≈ 0.5368754922


For more details, please check out resources in the Sources section.   I particularly like Oxford's A Very Short Introduction series.  



Casio fx-CG 50 Program Code:  CORDIC2


This algorithm is adopted from the Python example (see Wikipedia article).  


Deg

"ANGLE"?->Θ


0->I

0->A

1->X

0->Y

1->N


Lbl 0

tan^-1 (2^(-I))->T

If A<Θ

Then 

A+T->A

1->S

Else 

A-T->A

-1->S

IfEnd

N÷√(1+2^(-2*I))->N

X-S*2^(-I)*Y->P

Y+X*S*2^(-I)->Q

P->X

Q->Y


I+1->I


Abs (A-Θ)>1×10^-5=>Goto 0

X*N->X

Y*N->Y


ClrText

Blue Locate 1,3,"Θ: "

Blue Locate 5,3,Θ


Red Locate 1,4,"A: "

Red Locate 5,4,A


Black Locate 1,5,"cos :"

Black Locate 7,5,X

Black Locate 1,6,"sin :"

Black Locate 7,6,Y


Green Locate 1,7,I



CORDIC3:  Doing it without a Arctangent function


In reality, most of the time the arctangent function also has to be approximated.   There are many ways to approximate, which varying accuracy.   I used the fx-CG50's statistics mode to come up with a regression equation with the following lists:


x_list = sequence of 1÷(2^i) from i = 0 to i = 39 


y_list = sequence of arctan(1÷(2^i)) from i = 0 to i = 39


Of the regression models the fx-CG50 offers, the best regression model is quartic regression (4th-order polynomial):


y ≈ 9.43597784 × x^4 - 22.116232 × x^3 + 0.40116676 × x^2 + 57.2790727 × x + 1.4558 × 10^-5


Remember that I am working with degree angle measurement.  





Casio fx-CG 50 Program Code:  CORDIC3


Deg

"ANGLE"?->Θ


0->I

0->A

1->X

0->Y

1->N


Lbl 0

2^(-I)->K

9.43597784917955K^(4)-22.1162323028173K^(3)+

0.401166766193516K^2+57.2790727477367K+

1.45584715586876×10^-5->T


If A<Θ

Then 

A+T->A

1->S

Else 

A-T->A

-1->S

IfEnd

N÷√(1+2^(-2*I))->N

X-S*2^(-I)*Y->P

Y+X*S*2^(-I)->Q

P->X

Q->Y


I+1->I


Abs (A-Θ)>1×10^-5=>Goto 0

X*N->X

Y*N->Y


ClrText

Blue Locate 1,3,"Θ: "

Blue Locate 5,3,Θ


Red Locate 1,4,"A: "

Red Locate 5,4,A


Black Locate 1,5,"cos :"

Black Locate 7,5,X

Black Locate 1,6,"sin :"

Black Locate 7,6,Y


Green Locate 1,7,I


For the example:

 Θ ≈ 32.47119229°


Approximating Θ within 8 decimal places (10^-5) yields these results:

A:  32.47118598

cos A:  0.8436683456

sin A:  0.5368647154

24 terms used



Sources


"CORDIC"  Wikipedia.   Last Edited January 11, 2024.   Accessed January 12, 2024.  https://en.wikipedia.org/wiki/CORDIC


Brummelen, Glen Van.  Trigonometry:  A Very Short Introduction  Oxford University Press: Oxford, United Kingdom.  2020.  ISBN 978-0-19-881431-3



Eddie


All original content copyright, © 2011-2024.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 


Monday, December 25, 2023

TI Calculator App: RPN83P v0.8 by Brian T. Park

TI Calculator App:  RPN83P v0.8 by Brian T. Park








Reverse Polish Notation on the TI-83/TI-84 Plus


The calculator app, the RPN83P, programmed by Brian T. Park, is a calculator app that runs an HP 42S-like interface and gives the calculator a working, feature-filled RPN (Reverse Polish Notation) calculator.   The application works on the TI-83/84 Plus family that have a monochrome screen:  


TI-83 Plus

TI-83 Plus Silver Edition

TI-84 Plus

TI-84 Plus Silver Edition (I'm pretty sure that the SE keyboard attached to the original TI-Nspire would work too)

TI-84 Plus Pocket SE

TI-84 Pocket.fr


RPN83P will not work on the TI-84 Plus CE or CE Python.   


The app takes 32,768 bytes of flash memory, which is the size of 2 flash pages.  Two lists are also created, REGS for memory registers, STK for the four level stack.




Features


RPN83P operates with a traditional 4-level stack (X, Y, Z, T), and all the stack levels are displayed at once.  To quickly clear the stack, press [CLEAR] three times.  The display also shows current fix and angle settings.    


You can turn the calculator off from the app, when we turn the calculator on and run the app, it goes back to the screen where we left off.  




Common key mappings in RPN83P 


[ (  ] :  roll down the stack, R↓


[ ) ]:  swap the contents of the X and Y stacks, x<>y


[ 2nd ] ( ANS ):  Last x 


[ 2nd ] ( ENTRY ):  Show the entire precision of the value in the X stack, regardless of mode.   Press [ ENTER ] to exit. 


[ MATH ]:   Shows the MATH menus and the HELP facility.  The HELP file is a nine page document that shows the common keys and some of the functions of RPN83P.  Scroll the pages with the up and down arrow keys.  


[ MODE ]:  Sets the following modes:

Display:  FIX, SCI, ENG.   FIX 10 sets the float mode.

Angel:  DEG/RAD.  Degrees or Radians



The Math Functions



A full listing can is found here:  https://github.com/bxparks/rpn83p


Math-Math:  Cube, Cube Root, Universal Roots, ATAN2, extended logarithmic and exponential functions including E^X- (e^x - 1) and LN1+ (ln(x + 1))


Math-Num:  Percent, Percent Change, Integer, Fraction, Ceiling, Floor, Rounding, and other number facilities



The percent (%) function operates like the RPN calculators:


Y:   whole

X:   percentage


( % )  →

Y:   whole

X:   whole × percentage/100


To add a percent:  y [ ENTER ] x (% ) [ + ]

To subtract a percent:  y [ ENTER ]  x (%) [ - ]



The PRIM (prime) function is the smallest prime number function and used in with the division key  for prime factorization.


Example:  


254  (PRIM)   X:  2  

[ ÷ ]  (PRIM)   X:  1  (last number finishes the factorization)

[ ÷ ]    X :  127


254 = 2 × 127


Math-Prob:  Probability Functions


Math-Conv:  Deg/Rad, Rec/Pol, Hour/Hour-Minutes-Seconds conversions


Math-Help:  Access the Help documentation


Math-Base:  Base conversion, Boolean Logic, and computer science functions such as shift and rotate (HP 16C) 


Math-Hyp:  Hyperbolic functions


Math-Stat:  Statistics and regression fit:


LIN: linear

LOGF:  logarithmic

EXPF:  exponential

PWRF:  power

BEST:  takes the data and matches it to the best fitting regression


Math-Units:  Conversion of SI/English Units.   The conversions are in pairs and are cleverly organized.  For example:


>°C   ([ Y = ]) - °F>°C

>°F  ([WINDOW])  °C>°F


>hPA ([TRACE])  iHg>hPa

>iHg  ([GRAPH])  hPa>iHg


Scroll with the up and down arrow keys.  


There are 12 sets of conversions.  


Math-TVM:  Finance functions (see the Finance section)


Math-CLR:  Clearing functions, X, stack, registers, statistics, and time value of money registers


Math-MODE:  Mode settings


Math-STK:  Stack functions:  R↑, R↓, x<>y




The Finance Menu



The finance menu is accessed by [ MATH ], [ ↓ ], [ GRAPH ] (TVM).


The menus in finance are:


1st Page:  the Time Value of Money variables:


N   ([ Y = ])

I%YR  ([WINDOW])

PV ([ ZOOM ])

PMT ([TRACE])

FV  ([GRAPH])


2nd Page:  Utilities


P/YR   ([ Y = ]):   set payments per year

BEG  ([WINDOW]):  sets beginning of the period mode

END  ([ ZOOM ]):  sets end of period mode


CLTV  ([GRAPH]):   Clear TMV variables


3rd Page:  Solving for Interest Settings


IYR1   ([ Y = ]):   set first guess (default 0%)

IYR2  ([WINDOW]):  set second guess (default 100%)

TMAX ([ ZOOM ]):   iteration maximum (up to 255, default set at 15)


RSTV  ([GRAPH]):  Reset the solver


The standard cash flow convention is used:  negative for outflows (payments), positive for inflows (receipts).  This mode operates like the HP 17B financial calculator.  


To enter a number in the finance variables, enter the number and then press the corresponding soft key.


To solve for a variable, just press the soft key that corresponds to the variable.


Example:  Find the value of a savings fund if monthly deposits of $300 are made for 4 years.  The account pays at 3.1% rate.


[ ↓ ]  

( CLTV )


[ ↑ ]

4 [ ENTER ] 12 [ × ] ( N )

3.1 (I%YR)

300  [(-)] (PMT)

( FV )  


Result:   15309.85772




Where to get RPN83P 



The RPN83P application is a work in progress:  the version I am reviewing at the time of this post is v0.8.   Planned additions for 2024 include complex number support.  The application is subject to testing, change, and possible bug correction.  The RPN83P is an open source project, with the source code on GitHub (see GitHub page link below).


RPN83P GitHub page:

https://github.com/bxparks/rpn83p


Download the program here:

https://github.com/bxparks/rpn83p/releases


Thread on Museum of HP Calculators:

https://www.hpmuseum.org/forum/thread-20867.html


Special thanks and gratitude to Brian T. Park.  I definitely recommend the RPN83P application and am excited to see the future updates that will be added to this excellent app.  


Eddie


All original content copyright, © 2011-2024.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 


Saturday, April 3, 2021

Python on Calculators: Python Files used as Modules

Python on Calculators:  Python Files used as Modules

(and Programming Tutorials by Mike Dane)


Python Files Can Be Transferred Between Calculators


Python scripts that are calculators have the .py extension, allowing the files to be transferred to and from calculators and python computer applications.   The calculators that I have with Python are:


* Casio fx-9750GIII

* Casio fx-CG 50

* Numworks

* TI-NSpire CX II


To my knowledge, other calculators with a Python module are:


* TI-83 Premium CE Edition Python (primarily sold in France)

* Casio fx-9860GIII (the fx-9750GIII model sold outside of the United States)


Importing Python Files - Extended Functions (extf.py)


For a Python file to imported in a python script, the source python file must be included in the calculator.   


Example: 


I would be working on the file usethesource.py and I want to import the functions defined on the source.py to be used.  All is needed is an import statement:


import source


To use any of the functions from the source.py, use the source-dot prefix, like:


source.function()


Remember to type the name of the module, the decimal point, and the function.  The Casio calculators do not store custom functions in a catalog, while the NSpire CX and Numworks do include the custom functions in the variable list but not the module prefix.


The python scripts extf.py and extfdemo.py serves as an illustration:


extf.py


# create a custom module

# 2021-02-03 EWS


import math


# functions code


def radius(a,b):

  return math.sqrt(a**2+b**2)


def angle(a,b):

  t=math.acos(a/math.sqrt(a**2+b**2))

  t=t*180/math.pi

  if b<0:

    t=-t

  return t


def pchg(a,b):

  return (b-a)/a*100


def taxplus(a,b):

  return a*(1+.01*b)


def pvaf(a,b):  

  return (1-(1+.01*b)**(-a))/(.01*b)

  

def fvaf(a,b):

  return ((1+.01*b)**a-1)/(.01*b)


def rsum(a,b):

  return 1/(1/a+1/b)


# constants

# gravity of Earth

g=9.80665  

# speed of light

c=299792458


extfdemo.py


# 2020-02-03 EWS


# import extf.py file

import extf


print("Import any py file.")

print("The source file must")

print("be in memory.")


# \n new line

print("\n")

print("extf.radius(a,b)=abs(a+bi)")

print("extf.angle(a,b=arg(a+bi) in degrees")

print("a=11.25, b=13.36")

print(extf.radius(11.25,13.36))


# other functions

print("\n")

print("extf.pchg: percent change")

print("old,new")

print("pchg(78.95,96.95")

print("a =22.52,b=29")

print(extf.pchg(22.52,29))


print("\n")

print("present value annuity factor")

print("extf.pvaf(n,i%)")

print("pv = pmt*pvaf")

print("pmt=250, n=24, i%=4.5")

print(extf.pvaf(24,4.5)*250)


# more to discover

print("\n")

print("More to discover!")



You can download the two files here: 

https://drive.google.com/file/d/1-g7oIru7SGFvg-1IhM4CDr6lrmMK7Nrv/view?usp=sharing



Programming Tutorials by Mike Dane


Mike Dane, also known as Giraffe Academy, has many tutorials of current popular programming languages, including Python, C++, Javascript, and HTML. It is his channel where I learned about the ability to use any Python script to be imported as a module and other basic skills for Python.  


Dane has both playlists of short tutorials and a comprehensive tutorial of programming languages.  Please check out his YouTube page at:  https://www.youtube.com/c/GiraffeAcademy/videos


His web site:  https://www.mikedane.com


 

Eddie


All original content copyright, © 2011-2021.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 


Saturday, December 28, 2019

HP Prime: Base Conversions

HP Prime:  Base Conversions 

Introduction

The program BASEDEC converts a number in a base representation to decimal (base 10).   BASEDEC accepts any positive number, including non-integers up to base 36.  Each of the letters of alphabet represents a value:

A = 10
B = 11
C = 12
D = 13
E = 14
F = 15
...
Z = 35

BASEDEC takes two arguments: number to be converted, entered as a string;  and the base representation of the number.

Examples: 

BASEDEC("531",8) converts 531_8 to decimal.  Result:  345

BASEDEC("63.95A", 11) converts 63.95A_11 to decimal.  Result:  69.86701728

BASEDEC("4H", 18) converts 4H_18 to decimal.  Result:  89  (H = 17)

HP Prime Program BASEDEC

EXPORT BASEDEC(s,b)
BEGIN
// string, original base
// number can have fractional
// part in decimal form
// convert to base 10
// A=10,B=11,C=12...Z=35
// 2019-12-09 EWS
LOCAL l,pt,sl,sr,ls,lr,d;
LOCAL k,ac,al,dg,sd;
d:=0; // decimal
l:=DIM(s);
// split the string
pt:=INSTRING(s,".");
CASE
IF pt==0 THEN 
sl:=s;
sr:="";
END;
IF pt==1 THEN 
sl:="";
sr:=RIGHT(s,l-1);
END;
DEFAULT
sl:=LEFT(s,pt-1);
sr:=RIGHT(s,l-pt);
END;
// conversion to decimal
// integer part
ls:=DIM(sl);
IF ls>0 THEN
FOR k FROM 1 TO ls DO
sd:=MID(sl,k,1);
ac:=ASC(sd);
al:=ac(1);
dg:=when(al>64,al-55,EXPR(sd));
d:=d+dg*b^(ls-k);
END;
END;
// fractional part
lr:=DIM(sr);
IF lr>0 THEN
FOR k FROM 1 TO lr DO
sd:=MID(sr,k,1);
ac:=ASC(sd);
al:=ac(1);
dg:=when(al>64,al-55,EXPR(sd));
d:=d+dg/(b^k);
END;
END;
// return decimal
RETURN d;
END;

Introduction Part II

The program DECBASE converts a number in decimal form to any base you want.  DECBASE works will all positive real numbers.  There are three arguments to DECBASE:

*  The number in decimal form
*  The desired base
*  Accuracy:  how many decimal points do you want.  Conversions from decimal to other bases are approximations.   If the number is a integer (frac(number) = 0), then this argument is ignored (put 0).

The result is returned as a string.

Examples:

DECBASE(69, 6, 0) converts 69 to base 6.  Result:  "153"

DECBASE(5.875, 8, 4) converts 5.875 to base 8 to four decimals.  Result:  "5.7000"  (exact:  5.7_8)

DECBASE(635.05, 12, 12) converts 635.05 to base 12 to twelve decimals.  Result:  "44B.072497249724"   (exact:  44B.07249..._12  (7249 repeats))

HP Prime Program DECBASE

EXPORT DECBASE(d,b,acc)

BEGIN

// decimal, base, accuracy

// convert from decimal to base

// 2019-12-09 EWS

LOCAL n,f,nl,nd;

LOCAL nf,k,dg,s,fd;

n:=IP(d); // integer

f:=FP(d); // fractional

s:="";

// integer

IF n>0 THEN

nl:=IP(LOG(n)/LOG(b));

nd:=n;

FOR k FROM nl DOWNTO 0 DO

dg:=IP(nd/(b^k));

IF dg<10 THEN

s:=s+STRING(dg);

ELSE

s:=s+CHAR(dg+55);

END;

nd:=nd-dg*b^k;

END;

END;

// numbers < 1

IF n==0 THEN

s:="0";

END;

// fractional part

IF f>0 THEN

s:=s+".";

fd:=f;

FOR k FROM 1 TO acc DO

dg:=IP(fd*b^k);

IF dg<10 THEN

s:=s+STRING(dg);

ELSE

s:=s+CHAR(dg+55);

END;

fd:=fd-dg/(b^k);

END;

END;

RETURN s;

END;


Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

RPN: DM32 and DM42: Stopping Sight Distance (Metric)

RPN: DM32 and DM42: Stopping Sight Distance (Metric) The Stopping Sight Distance Formula – Derivation The stopping sight di...