LONGFLOAT for HP49

V. 3.86

2006-01-25

Multiple precision library with + - * / sqrt yx exp ln sin cos tan asin acos atan sinh cosh tanh and some auxiliary functions for both real and complex input and interval arithmetic (precision tracking) for real functions.

  • The precision is set in variable DIGITS.
  • Results are accurate to within ± 1 in last digit.
  • Integers, longfloats and long complex may generally be mixed.
  • Most algebraic functions may be automatically evaluated to multiple precision, set userflag 42.
  • Settings for deg/rad/grad is obeyed.
  • There are 4 different rounding modes, round, round down, round up and truncate.
  • Unit functions included thanks to Thomas Rast.
  • The library is to a great extent based on simplified algorithms of ZMLIB ( a fortran package for multiple precision).
  • is calculated using Werner Huysegoms implementation of Chudnowsky.
  • Constants like , ln(2), ln(10) are automatically stored when recalculated to higher precision in a library data.

I have done my best to keep the library bugfree, however as always, I recommend you back up your HP49 before using this library. I take no responsibility for the results.

Thanks to Werner Huysegoms and Mika Heiskanen for inspiration.
And thanks to Carsten Dominik and Peter Geelhoed for emacs.

Things not working, and which are used differently:

  • automatic list processing ,
  • lists are used for interval arithmetic instead.( but DOLIST/DOSUBS may be ok)

Please report any bugs/errors/ strange results to gjermund.skailand@yahoo.no
The library is free for noncommercial use, but I retain the copyright.

If anyone find this useful, please send me an email with programs etc, for possible future amendmends to the library.

Best regards
Gjermund Skailand
2003-08-01


SHORT USER GUIDE LONGFLOAT v. 3.0

All results is calculated within +/- 1 in last digit.
Precision is set in DIGITS variable, which may be global or local: « 23 'DIGITS' STO »

When the library is started from rightshift-lib menu, the current digits and rounding mode is shown in the display.

  • :20 indicates 20 digits and round to “even” (nearest except when 5 is being rounded), system flags 37 clear, system flag 38 clear.
  • <20 indicates 20 digits and floor rounding, system flags 37 set, system flag 38 clear.
  • >20 indicates 20 digits and ceiling rounding, system flags 37 clear, system flag 38 set.
  • x20 indicates 20 digits and truncate mode, system flags 37 set, system flag 38 set.
  • []20 indicates 20 digits and interval arithmetic, see other document, set userflag 42.

The above flags may also be set with FMode.

If any calculations are performed on interval numbers, then the lower bound is floor rounded, the upper bound is ceiling rounded, and then rounding mode is reset to round to nearest.

One may also use local variable to set precision temporarily, e.g. « 155 DIGITS « ... » »
Default precision is 20 digits.

Algebraic functions may be directly evaluated to multiple precision,
e.g. '2*y^2+2.5*SIN(*x)' FNUM. Thanks to the versatile CAS, longfloat numbers can also be included in the functions. Further, for real functions, without units, you can use interval arithmetic and get strict bounds on the result. Set userflag 42.

Units may be included in algebraics, and automatically evaluated to arbitrary precision,
e.g. '2*y^2+2.5*SIN(*x) × 1_cm UNum converting to real number with UR.

Max exponent for most functions are +/- 9999. For basic arithmetic ( + - * / square root, incl. LN), exponent is extended to +/- 999 999 999 (9 digits). Overflow or underflow should always generate an error ( i.e. "OVERFLOW", "NEGATIVE UNDERFLOW" etc.).

Longreal, longcomplex and integers may generally be mixed, however, it is e.g. faster to use integers when that is possible.

Eq. 2.5 2 FY^X is much faster than 2.5 2. FY^X, which is calculated as FEXP(FLN(2.5)*2)), and FSQ(2.5) is even faster.

For interval arithmetic, longreal_interval and integers max be generally mixed, see separate document.

Basic matrix functions have been implemented, + - * inverse / DET, see separate document.

There are 11 library menu pages (64  commands):

Longfloat menu page 1

Longfloat menu page 2

Longfloat menu page 3

Longfloat menu page 4

Longfloat menu page 5

Longfloat menu page 6

Longfloat menu page 7

Longfloat menu page 8

Longfloat menu page 9

 

Algebraic Functions

RF

string  longreal

zint longreal
real longreal
complex longcomplex
{real_lower real_upper}  longreal_interval

longreal real
longcomplex complex

longreal_interval {real_lower real_upper} 

FADD
FSUB
FMULT
FDIV
FNEG
FINV
FABS

longfloat longfloat

FSQRT

integer longreal
longreal longreal
longcomplex longcomplex
Example:
Square root of 2 with 155 digits:
2 « 155 DIGITS « RF FSQRT » »

Note: if longreal is negative, it is converted to longcomplex

FY^X

longreal integer longreal
longreal longreal longreal
longcomplex integer longcomplex
longcomplex longreal longcomplex
longcomplex longcomplex longcomplex.
If X is integer then fast binary squaring is used.

FXROOT

longreal integer longreal
integer integer longreal
Example:
-8 3 -2000000000000000000000000E-24
( the integer in stacklevel 2 is converted to longfloat prior to calculating the nroot).

Will error if the result is not real. e.g. -8 2 -> “Undefined Result”

FFC

longreal longreal longcomplex
longcomplex longreal longreal

Fi

longcomplex
Return imaginary 1, (0.,1.)

FPI

.
calculated by Werner Huysegoms implementation of Chudnovsky.

FEXP

zint longreal
longreal longreal
longcomplex longcomplex

FLN

zint longreal
longreal longreal
longcomplex longcomplex

FLNP1

longreal longreal
ln(y-1), use this when y is close to 1, and you have the value x = y-1

FIP

longreal integer

ZSqrt

lnteger integer 1. / 0.
Returns square root such that y-1 < x2 <= y
Flag 1 for exact , ie. y = x2

Trigonometric Functions

SIN
COS
TAN
ASIN
ACOS
ATAN

  • obeys deg, rad, grad mode
  • does not change from real to complex mode, enter complex number if you expect result to be complex, e.g. 1.2 FASIN "Bad argument value"

FSIN
FCOS
FTAN
FASIN
FACOS
FATAN

zint longreal
longreal longreal
longcomplex longcomplex

Hyperbolic Functions

HYPSIN
HYPCOS
HYPTAN

  • inverse functions not implemented
  • always assume rad mode
  • does not change from real to complex mode

FSINH
FCOSH
FTANH

zint longreal
Longreal longreal
longcomplex longcomplex

Interval Aritmetic

F<>

Longreal longreal_interval

Converts longreal to longreal_interval. Accuracy assumed to be +/- 1 of last digit in mantissa, regardless of current DIGITS precision.

Ex.

5789.E-3  {5788.E-3  5790.E-3}

578900.E-5  {578899.E-5 578901.E-5}

F+/-<>

Longreal longreal  longreal_interval

Average +/- accuracy converted to longreal_interval

Ex. DIGITS=6

1230.E-3   1570.E-7    { 122984.E-5 123016.E-3 }

Note: You can also make an longreal_interval directly by including two longreals in a list, ensure that it is {lower_boundary  upper_boundary}. There is no check, and if this is not correct then rounding will not be correct either. 

<>F+/-

Longreal_interval  longreal longreal

Returns average and +/- accuracy

Ex. {2300.E-3 2400.E-3 }   2350.E-3 500.E-4

Evaluation of Algebraic or Symbolic Functions

FSymb

Longreal ’longreal’
longcomplex ’longcomplex’

Longreal interval ’longreal interval’

Convert (quote) a longfloat real, complex or a longfloat interval number to symbolic. Lets you directly include longfloat constants in symbolic or algebraic equations, by using the CAS, which may then be evaluated by ->FNUM

FNUM

Evaluate algebraic to multiple precision floating point number. Works for most algebraic functions, but see below.

  • uses same assumptions as for the equivalent multiple precision number (e.g. always RAD mode)
  • The resulting precision for the algebraic is not guaranteed, unless you use interval_mode,
    however the precision for each separate command is within +/-1 in last digit.
  • Functions are converted to equivalent multiple precision functions whenever such exist.
  • subexpressions may contain user programs. Generally these should result in longfloat types (real or complex) .
  • real and complex numbers input are always converted to longreal or longcomplex.
  • integers in input are first evaluated as integer, if this fails, then integers are converted to longreals and a new evaluation is tried.
  • The algebraic with possible sub-expressions are recursively evaluated
  • Formal variables are not implemented, thus expressions like 'f(x)=.....' is not allowed, however, you may use user programs with local variables instead.
  • There is only a simple check for circular references, user may cancel with ON
  • be aware that reals and integers are treated differently, and the runtime and resulting type depends on that.
    e.g. type speed
    '2^2' integer fast
    '2.^2' longfloat fairly fast
    '2.^2.' longfloat slow [ evaluated as exp(2.*LN(2.)) ]
  • the following constants are recognized: and i ( imaginary symbol)
  • Symbolic matrix is implemented to the extent that algebraics may be inside a single symbolic matrix only. Use the CAS to evaluate/simplify first, but take care so that you don't evaluate any real functions to real values using the CAS.
    [[ 'sin(x)' 'COS(x)' ]] is allowed
    [[ 'sin(x)' 'COS(x)' ]] 'Y' STO
    'Y*Y' is not allowed, a warning/error message would be issued. Use the CAS to evaluate, inspect result and then
    'Y*Y' EVAL FNUM is allowable.

How it is currently implemented:

  • 1 the algebraic is converted to the equivalent rpn object and put on the runstream, then evaluated one by one
  • 2 for each object from runstream:
    • 21 any zints, longreals, longcomplex put on stack
    • 22 any real and complex numbers converted to longreal or longcomplex.
    • 23 symbolic evaluate recursively with FNUM
    • 24 global object (IDNT)  or local named variable,
      RCL ; error if type Matrix; convert to symbolic evaluate recursively with FNUM , error if undefined or not exist.
    • 25 secondary ( e.g. functions )
      • 251 check if longfloat equivalent function exist
        • 2511 try to evaluate longfloat function, if successful then ok
        • 2512 if not successful, convert object No. 1 on stack to longfloat type, then evualate. Error if not successful
      • 252 try to evaluate secondary (could be user-function). Error if not successful. Continue at 2 .

Download test examples to calc for more info

NRSOLVE

Solves an equation (finds a root) for a single variable by Newton-Raphson iteration.
The equation to be solved has to be derivable with regard to the variable to be solved for.
Dynamic doubling of precision is carried out.
In order to solve in a reasonable time, you should supply a good guess with 10-12 correct digits ( use the inbuilt solver to get a good guess)
If you supply a real guess, then 12 digits is assumed to be initial precision of guess. If you supply a longreal, then all digits of this guess is assumed for initial precision.
If you need to use longreal in formulas, I suggest you simply store it in a variable, and include that variable in the equation.

Input
Algebraic formula
Variable to be solved for ( Global or local, local is faster)
Guess ( real or longreal )

Output
Longreal
:ERR:Longreal

Ex.
25 'DIGITS' STO
'SIN(2*x)-.11'
'X'
.55111524
NRSOLVE

5511152499287331658823563E-26
:ERR:0.

( In about 16 seconds)
( 100 Digits in about 45 seconds )

 


UNITS

By Thomas Rast

UNum
converts reals, integers and algebraics with or without units to longreals and includes unit calculations too.
If the expression does not contain units, then a dummy unit is appended.
e.g.
35_kph UNum
9_s UNum
* UNum

Relativistic time when travelling at higher speeds:
If you like airplanes and happen to travel at 1000 km/h how does the time pass relative to an observerer at rest?
And after 100 years and you were still flying around, how many more seconds would have passed at the airport?

Relativistic time difference reltim(v) = 1/(1 - v2/c2)
Difference in seconds tdiff(v,t) = (1 - reltime(v)) * t

'reltim(v)'='1/(1-v^2/c^2)' DEF
1_c 'c' STO

1000._kph UNum (convert to longreal with units)
reltime UNum (you have a relative measurement (no units), use UNum)
1 UNum
-
100_yr UNum
*
UNum
'2709245019650046649115400.E-27*1_s'

UR
converts a longreal with units to real with units.

UR
2.709E-3_s

Answer: Not much, i.e. for an observer at rest the time would be 2.7 milli seconds more.
However, For GPS satelites, travelling slightly faster,
this kind of calculations are carried out on the GPS units in order to get more accurate position
fixes.

Be aware:
If your intermediate unit calculations end up with a number with NO units,
you need split the expression and use ->UNum otherwise the next units might disappear, and in the event of unconsistant units, you would not get that error message.
Unfortunately masd does not allow library commands to be compiled to functions

Alternative implementation as a userprogram

« SWAP UNum SWAP UNum v t « v reltime UNum 1 UNum - t * UNum » »
'tdiff' STO

1000._kph
100._yr
tdiff
UR
2.709E-3_s

Symbolic matrixes are not implemented.


Utility functions

FlSF

Longreal integer integer
Float to simple fraction
i.e. 12.5 1250000000 1000000000

SFCF

Integer integer list of integers
Simple fraction to continuous fraction
Uses DIGITS as approx precision. A continous fractions is a number described as
z = a0+ 1/( a1 + 1/ (a2+ 1/( a3+ ..)))

CFSF

List of integers simple fraction.

These functions may be used to find fraction approximations for longreals. The list of integers may be manually inspected and integers removed from the right end (generally cut just in front of a big integer)
e.g (DIGITS = 8)
FPI
FlSF SFCF
{ 3 7 15 1 354 2 6 1 }
354 is a big factor, truncate the list to
{3 7 15 1 }
CFSF
355 113 which is equivalent to 355/113 3.141592 which is with 6 decimals.

Square roots turn out to give periodic patterns when represented as continous fractions :
eg. With DIGITS = 100
2 RF FSQRT
FlSF SFCF
{ 1 2 2 2 2 2 2 2 ..... }