: # *-*-perl-*-* eval 'exec perl -S $0 ${1+"$@"}' if 0; # if running under some shell $format_debug = 0; q@ Hypercalc, Copyright (C) 1998-2008 Robert P. Munafo. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License (immediately following this text) for more details. To receive a copy of the GNU General Public License, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA For more perl goodness, go to mrob.com/pub/perl @; $GPL_VERSION_2 = < Copyright (C) This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA Also add information on how to contact you by electronic and paper mail. If the program is interactive, make it output a short notice like this when it starts in an interactive mode: Gnomovision version 69, Copyright (C) year name of author Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. This is free software, and you are welcome to redistribute it under certain conditions; type `show c' for details. The hypothetical commands `show w' and `show c' should show the appropriate parts of the General Public License. Of course, the commands you use may be called something other than `show w' and `show c'; they could even be mouse-clicks or menu items--whatever suits your program. You should also get your employer (if you work as a programmer) or your school, if any, to sign a "copyright disclaimer" for the program, if necessary. Here is a sample; alter the names: Yoyodyne, Inc., hereby disclaims all copyright interest in the program `Gnomovision' (which makes passes at compilers) written by James Hacker. , 1 April 1989 Ty Coon, President of Vice This General Public License does not permit incorporating your program into proprietary programs. If your program is a subroutine library, you may consider it more useful to permit linking proprietary applications with the library. If this is what you want to do, use the GNU Library General Public License instead of this License. endgpl $GPLNOTICE = " Hypercalc is distributed under the terms and conditions of the GNU General Public License, version 2, June 1991. Type 'help gpl' at the Hypercalc prompt for details. "; # # hypercalc -- interpreted calculator with special provisions to avoid # overflow # $tagline = " Go ahead -- just TRY to make me overflow! _ _ |_| . . ._ _ ._ _ ._ | _ | | | | | ) (-` | ( ,-| | ( ~ ~ 7 |~ ~' ~ ~ `~` ~ ~ -' Enter expressions, or type 'help' for help. "; # # Which is bigger: 27 ^ (86 !) or (27 ^ 86) ! # Most calculators can't even give the value of 27^86 or of 86!. # # With HyperCalc you can see that 27^86 is 1.251076x10^123, # and 86! is 2.422709x10^130. Some calculators can handle that -- # the current record-holder is AlCalc for the Pilot, which goes # as high as 10^32767 and can handle 9274! (9274 factorial) # # But no other calculator can tell you that # # (27 ^ 86) ! = 10 ^ (1.534607.. x 10 ^ 125) # # or that # # 27 ^ (86!) = 10 ^ (3.467778.. x 10 ^ 130) # # (in other words, the first has over 10^125 digits and the second, # with over 10^130 digits is "just a little bit" larger.) # Project notes follow <<'endquote'; Hypercalc is an open-source interpreted calculator program designed to calculate extremely large numbers (such as your phone number raised to the power of the factorial of the Federal budget deficit) without overflowing. It does this by using a modified form of the level-index number system with a radix of 1.0e300 Performance stats for other calculators Year Make and model Overflow Digits sin/cos (pi, e) radians 1973 TI SR-50 1.0e100 589,459 ? 1980 Sharp EL-5100 1.0e100 12? ? 1989 Casio fx-7500G 1.0e100 59,459 1.57e8 ? Casio fx-115D 1.0e100 6,45 1.57e8 = 1/2 * 1e8 * pi 1995 Casio CFX-9800G 1.0e100 5898,45904 1.57e8 = 1/2 * 1e8 * pi 1997 Pilot AlCalc 1.0e32768 32? n/a 1998 Casio fx-260 1.0e100 6,45 1.57e8 = 1/2 * 1e8 * pi 1998 Casio 9850 1.0e100 5898,45904 1.57e8 = 1/2 * 1e8 * pi 1998 Casio 7400 1.0e100 5898,45904 1.57e8 = 1/2 * 1e8 * pi 1998 Sharp EL-531L 1.0e100 58,44 1.7453e8 = 5/9 * 1e8 * pi 1998 TI-85 1.0e1000 5898.459 1.0e12 1998 TI-92 1.0e1000 5898,459 1.0e13 1999 TI-89 1.0e1000 5898,459 ? Mathematica 1.44e323228010 n/a ? 1998 Hypercalc 32768pt300 58979,45905 1.35e10 = 2^32 * pi (for Palm) 1999 Hypercalc (10^10)pt300 300 1.35e10 = 2^32 * pi (for UNIX) The overflow value for Hypercalc is so large it can't be represented in the standard way. If we use Hypercalc's internal "PT" format it's easy. Hypercalc handles numbers with absolute value greater than the range supported by the floating point library by storing the numbers in many different formats. When the numbers are within normal floating-point range (less than 10^300) they are stored in the normal floating-point format. Between 10^300 and 10^(10^300) they are stored as logarithms, and Logarithmic Number System (LNS) algorithms are used. When the logarithm gets too big to store as a floating point number, the log is taken again, and so on. An integer field is used to keep track of how many times the log has been taken. Here are some examples: pt-notation pt val represents 0pt1.0 0 1.0 1.0 0pt3.45e10 0 3.45e10 3.45 x 10^10 0pt1.0e299 0 1.0e299 1.0 x 10^299 0pt9.9e299 0 9.9e299 9.9 x 10^299 1pt300.0 1 300.0 1.0 x 10^300 1pt300.301 1 300.301 2.0 x 10^300 1pt301.0 1 301.0 1.0 x 10^301 1pt834.173 1 834.173 1.489 x 10^834 2pt79.0 2 79.0 10^(10^79) 3pt34.0 3 34.0 10^(10^(10^34)) 254pt1.0e10 254 1.0e10 10^(10^(10^(10^ ... ^(10^10) ... ))), where there are 256 10's 32767pt1.0e300 32767 1.0e300 10^(10^(... ^(10^300) ...)), with 32768 10's; the highest number Hypercalc can handle (on the Palm platform; UNIX Hypercalc can go to about (10^10)PT1.0e300) (To read about even larger numbers, go to www.mrob.com and click on "Large Numbers".) Each time we transition from the top of one pt range to the bottom of the next, about 2 1/2 digits of precision are lost as the information formerly stored in the exponent has to be absorbed by the mantissa. Then, as we proceed up the range digits are gradually gained back until we reach the top of the range and we once again have a 2.5 digit exponent. So, for example at the top of the pt=0 range the values are things like 1.23456789012345 x 10^299, and there are 53 binary digits of precision in the mantissa, or almost 16 decimal digits. Then we cross over into the pt1 range and store the log instead, which becomes a value like 301.456789012345 -- we still have 15+ digits to work with, but the first three correspond to the exponent of the number and there are only 12+ or 13 digits left for expressing the mantissa. Of course as we keep going up we get to values like 123456.789012345 (which represents 6.15 x 10^123456) we lose even more mantissa digits to exponent, but eventually we'll get to values like 123456789012345000000 = 1.2345... x 10^20, which represents 10^(1.2345... x 10^20) and as we go on up to even bigger numbers we see that since the exponent needs to be printed it once again holds information equivalent to 2 1/2 digits. This entire issue of variable numer of digits and the associated problems it causes with nonintuitive roundoff performance would be avoided if one used a "natural" PT storage format, where e is the base and the representation is such that the floating point value is always in the range [1.0, e]. So, for example, the number 143 would be represented as {2, 1.601979...} because e^(e^1.601979) is 143. Such a format would be unwieldy for normal calculations, however, because you'd have to keep doing e^x and ln(x) all over the place when doing simple calculations like 25 + 2. Revision history for PalmPilot version of Hypercalc: 1998101x Start project from "SampleCalc" example. 19981018 Fairly complete scientific calculator, except trig functions 19981021 Start implementing PT functions, get pt_exp and pt_mul working. 19981022 Implement add/subtract, pow, pow10/log10, and gamma. 19981024 Pretty much complete on the PT functions; they even handle infinities properly. Also, add a "tiny" font to print exponents when using the stdFont. 19981025 Refine the formatting code for PT1's and higher so it computes exactly how many digits of mantissa can be shown. Add some more buttons, but most not implemented yet. Implement rounding (incredibly complex!). Add arc, hyp, sto, rcl, 2nd, and 1/x keys (but only implement 1/x) 19981026 Add the same formatting refinements to PT0's, so it can print contents of memories (which have fewer pixels available). Implement Sto and Rcl keys. 19981028 Add hyperbolic functions and inverse trig (but not inverse hyperbolic trig) 19981030 Add inverse hyperbolic functions. 19981031 Put f_xxx and pt_xxx routines in their own files. Implement floating-point square root based on the grammar school algorithm (greatly increases speed of inverse trig functions). Revision history for Perl implementation of HyperCalc: 19990610 Start writing a simple Perl calculator program using a new concept: expression evaluation via regular expressions (I got the idea while writing the #top100# movie statistics program). Right now it just does + and @ (kind oflike *). 19990701 Break '+' operator out into a separate subroutine add1 (eventually all operators will be done this way) 19990720? Add all the code from the Pilot Hypercalc, to eventually merge and translate into Perl. 19970721 Parsing routine is fairly complete and now includes nested loop to handle parentheses. Subroutines for all four operators +, -, * and /. "e" and "%" in an expression represent 2.71828... and previous result, respectively. 19990725 Add split() and start writing first operator that handles PT types: p_add, pt_add, pt_addpos. 19990727.2125 Do lots of porting work: put all routines in "proper" (Pascal) order; lots of global replaces to change things like "x.pt" to "$x_pt"; replace Taylor and Newton algorithms with builtin functions where available; minimum work to get pt_addpos working. It now properly adds 1e300 + 1e300 (and gets 1pt300.3010299...) 19990727.2154 pt_add fully works; pt_div works. 19990727.2235 pt_sub and pt_mul work now. Output formatting handles some of the special cases to print values like 1p2345.6789 as 4.77 * 10^2345 rather than as "1 PT 2345.6789" 19990727.2722 pt_ln works; parser handles ln() and log(). 19990728.1333 It now handles exp() and pow(), so I can compute really big values without lots of repetitious keystrokes. 19990728.25xx eval2() now stores all operator results into an array, and stores the array index into the expression string. This is to avoid numbers getting converted from strings into floating point and back again, and that dramatically reduces roundoff error. 19990729 Start editing all the f_xxx routines so the primitive floating-point type can be changed easily later. This involves implementing a minimal set of "primitives" like f_int, f_le, f_neg, f_mul, etc. and making all the other f_xxx routines do all their operations by calling these primitives. Also, inline constants like "10" are replaced with globals. 19990801 &pt_root and &pt_log_n work. All of the &f_xxx routines are 'primitivized', but &pt_xxx routines still need some work. 19990801.2529 Add "debug" command. Put most of &f_xxx primitives inside $f64_prim so they can be defined and redefined via &exec. Create $g_pt_inf to distinguish uses of infinity in PT field from its uses in VAL field. A few other changes to support switching VAL primitive precision. Make it auto-promote inlines like "23e456". 19990802 Pretty much finished making the &pt_xxx routines call &f_xxx routines explicitly. 1999080x Use open2() to launch #bc#. Write bce. 19990805 Write &fbc_fix2sci, &fbc_split, &f_cmp, comparison primitives, &f_neg, &me_magcompare, &m_truncround, &me_addpos, &f_add, &f_sub. &fbc_encode renamed to &fbc_sci2fix. Redirect stderr when launching #bc#. 1999080x Write &me_subpos 19990811 Add HC_LOG debug log, lots of calls to dbg1. Fix lots of bugs. Write bc version of &f_mul and &f_div. 19991015 Fix bug that caused small PT1's to be printed as e.g. 10 ^ 301.30103. Make dbg1flag a bitmask to allow debugging functions, expression parsing, or both. 19991017 Add variables (currently limited to all-alphabetic starting with "v") 19991018 Change single-letter function abbreviations and special letters like 'e', 'p' etc. to uppercase, to clear the lowercase namespace for use by user variables. 19991019 Fix some bugs relating to infinity handling and conversion in fbc routines. Four basic functions almost work (subtraction still seems to have problems) 19991117 Variables no longer need to start with 'v'. Add 'sqrt' function. 19991124 Combine parsing of e/pi/phi with the variable and function parsing; add error-check for undefined variables. 20000120 Write fbc versions of f_ln and f_exp; fix bugs in fix2sci and sci2fix; it now correctly computes 2^100 in scale 30. Fix bugs in switching back and forth between f64 and fbc. 20000206 Fix bug that prevented "sqrt(1+2)" from working 20000304 sqrt now goes through f_sqrt. Fix bugs that made bc hi_init not compute g_pi properly. 20000728 Remove dependency on "rpmlib.pl" 20010102 Add ERASE_BS test. 20010103 Clean up internals of eval_2. Fix "right-to-left precedence bug": "4-3-2" used to give "3", and "4/3/2" used to give "2.66667". I am deliberately leaving exponents that way: 4^3^2 still gives 262144. 20010107 Fix bugs: 2+2/(1+1) gave 2; 7^-1 didn't parse; scale 50,27^27 printed in scientific notation. Write &pt_roundup. Fix &prnt1 handling of high PT1's. fbc-based PT calculation is actually usable now! 20010108 Add history array and define_hist. Conversion across scale changes works, at least in the cases I checked. Fix bug in eval_2: sqrt() and other functions had become broken as a result of yesterday's fixes. Clean up fbc version of f_gamma a little, but it still suffers from a fundamental limit of the Stirling formula method, which basically requires that the number being factorialed must be at least as big as the 15th root of 10^curscale. Combined with the current limit of 1.0e300 for the fbc float data type, that means we can't get more than 33 digits of accuracy out of the f_gamma function. Increasing the exponent limit would fix it, but that poses another problem with the scaling loop -- for 50 digits of accuracy, the scaling loop has to loop 2154 times (because 2154 = 10^(50/15)). 2001010x Finish implementing #format# command. 20010109 Fix bug that made history list usable only for first 9 items. 20010115 Write init_pi_2, which calculates g_pi much more quickly. Decrease gammalim. 20010116 Add input history. 20010117 Change letters I/H for input and output history to C/R (commands and results). 20010210 Fix "c2" in case where c2 is a variable assignment, and add ';' symbol to separate commands. 20010216 Add ability to take "1e9" as input (used to require "1.e9") 20010521 Make 'x' a synonym for '*'. This works pretty well, in fact you can even define a variable 'x', and the expressions '2 x 4', '2 x x' and 'x x x' all do the right thing! But, that's not recommended. Also, change default output format to format 1, and make it print multiplication as 'x' because it looks better. 20010521 Map [] in input to (). This almost solves the problem of having output and input formats match -- the one missing piece is allowing the user to type 'PT', such as '3 PT 1.2 x 10^45'. 20010530 Almost fix the ambiguity of '!!': You can now type "4!!" and it will give you (4!)!, rather than '4' followed by the previous typed line. 20010601 When ';' is present in input, print each of the commands with its 'C# =' label as they're being added to input_history. 20010610 Detect presence of UNIX and doesn't try to run #bc# if not on UNIX. 20010613 Fix some of the bugs in handling of '-'. Add &pt_negate 20011026 Fix some bugs in command history expansion. 20011104 Add autodetect of ^H and call stty erase if they type it (Unix only) 20020129 Move automatic stty erase fix to subroutine fixerase. 20020301 Read first expression from command line. 20020305 Fix some bugs in rounding and prnt2 -- but it still has the problem that "scale=15" prints the same number of digits as the default "scale=14" 20020306 Now can put multiple commands including 'scale=' and 'quit' on command line 20020711 Convert tabs to spaces in input 20051228 Begin to implement BASIC interpreter 20051229 Add OLD, SAVE and indenting 20050102 Disambiguate input like "10 x = 2" and "10 x 2" (the first is a BASIC line, the second is a calculation) 20060716 Add &format4 and &format6. 20060719 Add &format5. 20061112 In canon_sci, remove trailing zeros from mantissa. (This affects format5, and might not be what we want in all cases, but it fixes cases like 10^27) 20061126 FOR..NEXT works enough to compute, e.g. the value of Mega. Also add CAT and rudimentary version of PRINT. 20061129 Try to fix broken "x = x + 1"; this might cause other bugs. Make "old" and "save" conscious of current file name. 20070513 Add #help debug#; slight changes aimed at an eventual merge of f64_prim and fbc_prim 20070621 Fix a flaw in format5 20070624 Add mantissa precision to range [1e5..1e7) in format5 20070629 Add BASIC help 20070809 Make format5's output a little more useful 20071224 Add ability to input time in HH:MM:SS.frac format, and format7 to print out in that format. 20080204 Fix bug that prevented nested loops from working. Make BASIC statements not print if they are followed by a semicolon. 20080215 Accept abbreviated PT notation (e.g. "2P100" for Googolplex) in input. 20080508 A couple fixes intended to make negative exponents in scientific notation (like "1.23e-4") work again. 20081001 Slight improvements to #help var# 20100113 Add #help format#. Add a missing pair of {} in prnt1 case p2_e. Remove extra () and [] from most formats. 20100204 Add a check for Cygwin and treat it as UNIX (for now; we really need to also check for #bc# on all "Unix" systems) 20100221 Fix a couple bugs related to creating and using directories. Eliminate reference to my personal environment variable FROMWORK. 20100318 Upgrade PRINT statement to support literal strings and multiple arguments separated either by comma or semicolon (the latter suppresses the trailing newline). Also, a trailing semicolon on ordinary statements, LET and FOR suppresses printing of the value that has been computed. 20100329 Add $pt_str; handle $apt greater than 9 in &format5. 20100601 &docmd3: Don't allow user to "set" the constants pi, e or phi 20100603 &docmd3: setvar errors now stop BASIC; detect setting r1, etc. 20100611 Add $format_debug and several changes to &split, &prnt2, &canon_sci, &format5, and &prnt1 which collectively make format5 labels more accurate, correctly handle denormalized scientific notation input like "13.72e9", and properly deal with infinity (e.g. from inputting "1/0"). Also, &BASIC_old now allows whitespace blank lines anywhere in its input. 20100612 Allow '_' in variable names and program names. 20100614 Turn most commands to lowercase (except lines inserted into BASIC programs, and quoted strings in PRINT statements). 20100705 Bug fix in mantissa rounding in &canon_sci. TO DO and FUTURE IMPROVEMENTS Presently, 10^10^7 prints "1.00000012{x}10^{10000000}". The erroneous low-order digits can be eliminated by the following method: Choose an epsilon e. Convert (x-e), x, and (x+e) to strings Compare "(x-e)" to "x" and keep all digits that are identical. Compare "x" to "(x+e)" and keep all digits that are identical. Keep whichever of the preceding strings yields more digits. When a value "rolls over" from one PT level to the next (for example, when multiplying 10^299 by 100) the number of mantissa bits diminishes (or increases) by about 10 bits (because the rollover cutoff value is approximately 2^2^10). Error introduced during such an operation can propagate from intermediate results into any final result, so th epsilon should be chosen as if the mantissa precision is 43 bits, rather than the actual precision of 53 bits. If the above doesn't always work, use a slightly bigger epsilon. Regardless of what is displayed, retain all the bits of each result for use in subsequent calculations (I think hypercalc already does this) For more precision in f_gamma and pt_gamma, perhaps I can use the Lanczos series given in A090674 and A090675. Details and analysis are in "Pugh 2004 Analysis.pdf". There is a good practical explanation at www.rskey.org/lanczos.htm See also en.wikipedia.org/wiki/Lanczos_approximation endquote $hd = $ENV{"HOME"}; $p_var = "[a-z][a-z_0-9]*"; $p_nonvar = "[^a-z_0-9]"; $p_arg = "[0-9.EINPR]+"; $p_nonarg = "[^0-9.EINPR]"; # The special letters are: # # A - arctan # B - cos() # C - input history string (during input only) # E - for scientific notation # I - intermediate result # K - ln() # L - log() # M - logn() # N - a minus sign within a number # P - PT # Q - sqrt() # R - history result # S - sin() # T - tan() # V - root() # X - exp() # x - Synonym for * sub seterase { # Note: These environment variables are nonstandard. I really need to find # a standard way to test which stty erase needs to be set. It depends # on your operating system and what type of terminal connection you're # running in (xterm, telnet, ssh, etc.) Suggestions are welcome! $erase_bs = $ENV{"ERASE_BS"}; $erase_del = $ENV{"ERASE_DEL"}; if ($erase_del) { system("stty erase '^?'"); } elsif ($erase_bs) { system("stty erase '^H'"); } else { system("stty erase '^?'"); } } # This is from #mira#; it automatically detects and handle improper # setting of backspace or delete erase character sub fixerase { local($l) = @_; if ($l =~ m/\010/) { while ($l =~ m/.\010/) { $l =~ s/.\010//; } print "stty erase ^H\n"; system "stty erase "; } elsif ($l =~ m/\177/) { while ($l =~ m/.\177/) { $l =~ s/.\177//; } print "stty erase ^?\n"; system "stty erase "; } $l =~ s/\011/ /g; return $l; } sub dbg1 { local($a, $mask) = @_; # $a =~ s/\|//g; $a .= "\n"; if ($dbg1flag & $mask) { print $a; } print HC_LOG $a; } #--------------------------------------------------- F64 primitives # # These routines, originally from "hcfloat.c", perform primitive # functions on floating-point values. In most cases the function # is built-in to perl, but we do lots of special checks for # things like infinity and arguments that generate imaginary or # complex answers. # # perl provides: abs, atan2, cos, exp, int, log, sin, sqrt $f64_prim = <<'endquote'; # These are the full primitives sub f_int { local($a) = @_; return int($a); } sub f_lt { local($a, $b) = @_; return ($a < $b); } sub f_le { local($a, $b) = @_; return ($a <= $b); } sub f_eq { local($a, $b) = @_; return ($a == $b); } sub f_ge { local($a, $b) = @_; return ($a >= $b); } sub f_gt { local($a, $b) = @_; return ($a > $b); } sub f_ne { local($a, $b) = @_; return ($a != $b); } sub f_neg { local($a) = @_; return (- $a); } sub f_add { local($a, $b) = @_; return ($a + $b); } sub f_mul { local($a, $b) = @_; return ($a * $b); } # This divide routine catches the divide by zero case. sub f_div { local($n, $d) = @_; if (&f_eq($d, $g_0)) { return $g_inf; } return ($n / $d); } # following are the partial primitives. sub f_sub { local($a, $b) = @_; return ($a - $b); # return &f_add($a, &f_neg($b)); } # Following are the non-primitives # Our exp routine tests the infinity and overflow cases itself, just in # case the library function doesn't. sub f_exp { local($x) = @_; local($v); &dbg1("$curprim f_exp $x", 1); # Handle infinities first if (&f_eq($x, $g_inf)) { return $g_inf; } elsif (&f_eq($x, $g_ninf)) { return $g_0; } # Handle the extreme overflow cases. if (&f_lt($x, $g_n2000)) { return $g_0; } elsif (&f_gt($x, $g_2000)) { return $g_inf; } $v = exp($x); return $v; } # Our ln routine is strange -- it accepts negative arguments and returns # the magnitude of the answer. It also special-cases infinities and zero. sub f_ln { local($x) = @_; local($v); &dbg1("$curprim f_ln $x", 1); # Handle 0 and infinities first if (&f_eq($x, $g_0)) { return $g_ninf; } elsif (&f_eq($x, $g_inf)) { return $g_inf; } elsif (&f_eq($x, $g_ninf)) { return $g_0; } # If the argument is negative we make it positive, because # we want to return the real component of the logarithm of the number, # and for any complex number, the real component of the log is the log # of the number's magnitude. if (&f_le($x, $g_0)) { $x = &f_neg($x); } $v = log($x); return $v; } sub f_pow10 { local($n) = @_; &dbg1("$curprim f_pow10 $n", 1); return &f_exp(&f_mul($g_ln_10, $n)); } sub f_log_n { local($n, $b) = @_; &dbg1("$curprim f_log_n v{ $b } $n", 1); return &f_div(&f_ln($n), &f_ln($b)); } sub f_log10 { local($n) = @_; &dbg1("$curprim f_log10 $n", 1); return &f_mul(&f_ln($n), $g_log10_e); } # sqrt(2 pi n) * n^n * e^(-n) * e^(1/12n) # 1/2 ln(2 pi n) + n ln(n) - n + 1/12n - 1/1260 n^3 + ... # t = n - 1; # l = 1/2 ln(2 pi) + (n + 1/2) ln(n) - n + 1/(12 n) - 1/(360 n^3) + ... # g = e^l sub f_gamma { local($n) = @_; local($l, $scal1, $n2, $np); &dbg1("$curprim f_gamma $n", 1); # For very low negative arguments, the gamma function is so # close to zero that we treat it as zero. However, for negative # integers it's infinite, so we handle that case explicitly. if (&f_lt($n, $g_n50)) { if (&f_eq($n, &f_int($n))) { return $g_inf; } return $g_0; } # Our approximation formula doesn't work well for small values, so # we use the recurrance relation gamma(n+1) = n gamma(n). $scal1 = $g_1; while ($n < $g_10) { # "10" depends on precision $scal1 = &f_mul($scal1, $n); $n = &f_add($n, $g_1); } # Since we're using Stirling's series for factorials, we have # to subtract 1 to make it be a gamma function series. $n = &f_sub($n, $g_1); $l = $g_05_ln_2pi; $l = &f_add($l, &f_add($n, $g_0_5) * &f_ln($n)); $l = &f_sub($l, $n); $n2 = &f_mul($n, $n); $np = $n; $l = &f_add($l, &f_div($g_1, &f_mul($g_12, $np))); $np = &f_mul($np, $n2); $l = &f_sub($l, &f_div($g_1, &f_mul($g_360, $np))); $np = &f_mul($np, $n2); $l = &f_add($l, &f_div($g_1, &f_mul($g_1260, $np))); $np = &f_mul($np, $n2); $l = &f_sub($l, &f_div($g_1, &f_mul($g_1680, $np))); $np = &f_mul($np, $n2); $l = &f_add($l, &f_div($g_1, &f_mul($g_1188, $np))); $np = &f_mul($np, $n2); $l = &f_sub($l, &f_div($g_691, &f_mul($g_360360, $np))); $np = &f_mul($np, $n2); $l = &f_add($l, &f_div($g_7, &f_mul($g_1092, $np))); $np = &f_mul($np, $n2); $l = &f_sub($l, &f_div($g_3617, &f_mul($g_122400, $np))); return (&f_div(&f_exp($l), $scal1)); } # square root routine handles infinities and negative arguments properly sub f_sqrt { local($x) = @_; &dbg1("$curprim f_sqrt $x", 1); # If the argument is negative we return 0, because 0 is # the real component of the square root of a negative number. if (&f_lt($x, $g_0)) { return $g_0; } elsif (&f_eq($x, $g_inf)) { return $g_inf; } elsif (&f_eq($x, $g_ninf)) { return $g_inf; } return sqrt($x); } # The sine algorithm uses the identity # sin(a+b) = sin a cos b + cos a sin b # The value being sin'd is broken into two pieces one of which is a multiple # of pi/16 and the other is between -pi/32 and pi/32. Then we can use the # Taylor series to get the sin of the latter part and a table for the # first part sub f_sin { local($x) = @_; &dbg1("$curprim f_sin $x", 1); # Handle arguments out of range # Note: The best way to handle things like this is to extend the # numeric representation to include an "error" term, expressing the # standard deviation of the computed value vs. the actual value. # Then, the error of the sin function is the error of the argument # times the cosine, but converging on +-1.0 when the error of the # argument is 2 pi or larger. In the limit case the computed value # will be expressed as 0.0 and the error as 1.0. # # Just for fun, I went to an office supply store and tested 5 # current-model scientific calculators (covering all price ranges) # to see what they did when you take sin or cos of larger and larger # numbers. Each calculator had a different cutoff value, but all # seemed to choose a value that was related to the calculator's # internal precision. The five cutoff values were: # 5*10^7 pi # 4.45*10^9 # 1*10^10 # 1*10^12 # 1*10^13 if (&f_gt($x, $g_sin_limit)) { return $g_0; } elsif (&f_lt($x, &f_neg($g_sin_limit))) { return $g_0; } return sin($x); } sub f_cos { local($x) = @_; &dbg1("$curprim f_cos $x", 1); # Handle out-of-range if (&f_gt($x, $g_sin_limit)) { return $g_1; } elsif (&f_lt($x, &f_neg($g_sin_limit))) { return $g_1; } return cos($x); } sub f_tan { local($x) = @_; &dbg1("$curprim f_tan $x", 1); # Handle out-of-range if (&f_gt($x, $g_sin_limit)) { return $g_inf; } elsif (&f_lt($x, &f_neg($g_sin_limit))) { return $g_inf; } return &f_div(sin($x),cos($x)); } # Arctangent function. # # Av{0} = 1 / sqrt(1 + x^{2}) , Bv{0} = 1 # Av{n+1} = (Av{n} + Bv{n}) / 2 , # Bv{n+1} = sqrt(Av{n+1} * Bv{n}) # arctan(x) = LIMIT [x * Av{0} / Av(n)] sub f_arctan { local($x) = @_; local($a, $b, $a0, $ao, $i); &dbg1("$curprim f_arctan $x", 1); # A few special cases, the x==0 case is just for speed if (&f_eq($x, $g_0)) { return $g_0; } elsif (&f_eq($x, $g_inf)) { return &f_div($g_pi, $g_2); } elsif (&f_eq($x, $g_ninf)) { return &f_div($g_pi, $g_n2); } $a0 = &f_div($g_1, &f_sqrt(&f_add($g_1, &f_mul($x, $x)))); $a = $a0; $b = $g_1; # Set up ao to be different from a so the loop has a chance to start $ao = &f_sub($a, $g_1); for($i=0; ($i<25) && (&f_ne($a, $ao)); $i++) { $ao = $a; $a = &f_div(&f_add($a, $b), $g_2); $b = &f_sqrt(&f_mul($a, $b)); } return &f_div(&f_mul($x, $a0), $a); } # arcsin(x) = arctan(x / sqrt(1 - x*x)) sub f_arcsin { local($x) = @_; local($a); &dbg1("$curprim f_arcsin $x", 1); $a = &f_mul($x, $x); if (&f_gt($a, $g_1)) { return $g_0; } return &f_arctan(&f_div($x, &f_sqrt(&f_sub(1, $a)))); } # arccos(x) = pi/2 - arctan(x / sqrt(1 - x*x)) sub f_arccos { local($x) = @_; local($a); &dbg1("$curprim f_arccos $x", 1); $a = &f_mul($x, $x); if (&f_gt($a, $g_1)) { return $g_0; } return &f_div($g_pi, $g_2) - &f_arctan(&f_div($x, &f_sqrt(&f_sub($g_1, $a)))); } # pow function is more complex than exp and pow10, because # it has to deal with any base. Specifically, if the base is negative # it computes the angle and magnitude of the complex answer # so it can calculate the real component. sub f_pow { local($b, $e) = @_; local($mag, $scal2); &dbg1("$curprim f_pow $b $e", 1); # Handle negative bases if (&f_lt($b, $g_0)) { # The angle of the base is pi, so the angle of the # exponent will be pi * e. Take cosine of this to # determine how much of the answer will be in the # real component. $b = &f_neg($b); $scal2 = &f_cos(&f_mul($g_pi, $e)); } else { $scal2 = $g_1; } # Compute magnitude of answer $mag = &f_exp(&f_mul(&f_ln($b), $e)); # If we got an overflow, we might be able to salvage # the computation by using logarithms. This happens when we compute # -20.0 to the power of 238.49999999999 (exactly 10 9's there) # The magnitude is 1.97e310, which overflows, but the scale # is 3.1416e-10 so the answer is actually well within # range (and even within PT0 range) if (&f_gt($mag, $g_ovr1)) { if (&f_gt($scal2, $g_0)) { $mag = &f_add(&f_mul(&f_ln($b), $e), &f_ln($scal2)); return &f_exp($mag); } else { $mag = &f_add(&f_mul(&f_ln($b), $e), &f_ln(&f_neg($scal2))); return &f_neg(&f_exp($mag)); } } else { # Return real component of answer return &f_mul($scal2, $mag); } } # Root function contains the same complexity (pun intended) as # the pow function. sub f_root { local($b, $r) = @_; local($mag, $scal3); &dbg1("$curprim f_root $r v/ $b", 1); # Handle negative bases if (&f_lt($b, $g_0)) { # The angle of the base is pi, so the angle of the # answer will be pi / r. Take cosine of this to # determine how much of the answer will be in the # real component. $b = &f_neg($b); $scal3 = &f_cos(&f_div($g_pi, $r)); } else { $scal3 = $g_1; } $mag = &f_exp(&f_div(&f_ln($b), $r)); return &f_mul($scal3, $mag); } endquote sub f64_init { eval($f64_prim); # These constants are mostly used by the various transcendental functions $g_0 = 0.0; $g_0_5 = 0.5; $g_1 = 1.0; $g_n1 = -1.0; $g_2 = 2.0; $g_n2 = -2.0; $g_4 = 4.0; $g_5 = 5.0; $g_6 = 6.0; $g_7 = 7.0; $g_8 = 8.0; $g_10 = 10.0; $g_12 = 12.0; $g_15 = 15.0; $g_16 = 16.0; $g_n50 = -50.0; $g_360 = 360.0; $g_691 = 691.0; $g_1092 = 1092.0; $g_1188 = 1188.0; $g_1260 = 1260.0; $g_1680 = 1680.0; $g_n2000 = -2000.0; $g_2000 = 2000.0; $g_3617 = 3617.0; $g_122400 = 122400.0; $g_360360 = 360360.0; $g_2p32m1 = 4294967295.0; $g_pi = 3.1415926535897932; } sub f64_close { # undef(&f_int); # undef(&f_lt); } # --------------------------------------------------- FBC primitives # bce is the "engine" that feeds problems to the running #bc# process # and gets results. sub bce { local($expr, $expect_response) = @_; local($result, $l); print $bc_out ($expr . "\n"); $l = ($expect_response ? "\\" : ""); while ($l =~ m|\\$|) { $l = <$bc_in>; chomp $l; $result .= $l; $result =~ s|\\$||; } &dbg1("bce : $expr => $result", 4); return $result; } # fbc_norm normalizes a number in scientific format. For example, it converts # "0.02e3" to "2.0e1", and "23e45" to "2.3e46". sub fbc_norm { local($x) = @_; local($xs, $xm, $xe); local($s, $xi, $xf); local($mi, $mf, $e); &dbg1("fbc_norm $x", 1); if ($x =~ m/inf/i) { &dbg1("fbc_norm d: $x", 1); return $x; } # make sure there is a sign if (!($x =~ m/^[-+]/)) { $x = "+" . $x; } # extract the three parts of the number if ($x =~ m|^([+-])([0-9.]+)E([-+0-9]+)$|) { $xs = $1; $xm = $2; $xe = $3; } elsif ($x =~ m|^([+-])([0-9.]+)$|) { $xs = $1; $xm = $2; $xe = 0; } else { &dbg1("fbc_norm illegal input |$x|", 1); } # Make sure there is a decimal point if (!($xm =~ m|\.|)) { $xm .= ".0"; } # make sure it doesn't start with a bare decimal point if ($xm =~ m/^\./) { $xm = "0" . $xm; } # extract integer and fraction portions if ($xm =~ m|^([0-9]+)\.([0-9]+)$|) { $xi = $1; $xf = $2; &dbg1("fbc_norm xi $xi xf $xf", 1); } else { print "fbc_norm illegal mantissa |$x|\n"; } # handle cases where mantissa was less than 1 (normalize by decreasing exponent) if ($xi == 0) { $x = $xi . $xf; if ($x =~ m|^(0+[1-9])([0-9]*)$|) { $mi = $1; $mf = $2; $e = 1 - length($mi); $mi =~ s|^0+||; $x = $xs . $mi . "." . $mf . "E" . ($e + $xe); &dbg1("fbc_norm a: $x", 1); return $x; } # else, it was zero $x = "0.0E0"; &dbg1("fbc_norm b: $x", 1); return $x; } # else $e = length($xi) - 1; $x = $xi . $xf; $x =~ s|^([0-9])|$1\.|; $e = $e + $xe; # check for overflow if ($e > 300) { $x = $xs . "inf"; &dbg1("fbc_norm e: $x", 1); return $x; } $x = $xs . $x . "E" . $e; &dbg1("fbc_norm c: $x", 1); return $x; } # fix2sci converts from fixed-point to scientific notation. It's very simple # because all the hard work is done by fbc_norm sub fbc_fix2sci { local($x) = @_; &dbg1("fbc_fix2sci $x", 1); $x =~ s/e/E/g; $x =~ s/ //g; if ($x =~ m/inf/i) { &dbg1("fbc_fix2sci a: $x", 1); return $x; } if (!($x =~ m/E/)) { $x = $x . "E+0"; } $x = &fbc_norm($x); &dbg1("fbc_fix2sci b: $x", 1); return $x; } # fbc_split takes a number in any format and returns sign, mantissa and # exponent. It is used by most of the nontrivial primitives. sub fbc_split { local($x) = @_; local($s, $m, $e); &dbg1("fbc_split $x", 1); $x =~ s/ //g; # map infinity %%% might want to just return "inf" in the exponent field # and make callers check for infinity if ($x =~ m/^ninf$/i) { &dbg1("fbc_split a: ( - , 1 , $x )", 1); return("-", 1, $x); } elsif ($x =~ m/^inf$/i) { &dbg1("fbc_split b: ( + , 1 , $x )", 1); return("+", 1, $x); } # Convert to scientific if appropriate if (!($x =~ m/E/)) { $x = &fbc_fix2sci($x); } # Make sure it has a sign if (!($x =~ m|^[+-]|)) { $x = "+" . $x; } # get sign, mantissa and exponent if ($x =~ m|^([+-])([0-9.]+)E([-+0-9]+)$|) { $s = $1; $m = $2; $e = $3; } else { &dbg1("fbc_split illegal input |$x|", 1); } &dbg1("fbc_split c: ( $s , $m , $e )", 1); return($s, $m, $e); } # sci2fix converts from 1.23e10 format to something like # 12300000000.0 # # It takes scientific or fixed input. sub fbc_sci2fix { local($x, $intonly) = @_; local($s, $m, $mi, $mf, $e); local($ri, $rf, $rv); &dbg1("fbc_sci2fix $x $intonly", 1); # get sign, mantissa and exponent if ($x =~ m|^([+-])([.0-9]+)[eE]([+-][0-9]+)$|) { $s = $1; $m = $2; $e = $3; } elsif ($x =~ m|^([+-])([.0-9]+)[eE]([0-9]+)$|) { $s = $1; $m = $2; $e = $3; } elsif ($x =~ m|^([+-])([.0-9]+)$|) { $s = $1; $m = $2; $e = 0; } elsif ($x =~ m|^([.0-9]+)$|) { $s = "+"; $m = $1; $e = 0; } else { print "fbc_sci2fix illegal input |$x|\n"; } # Get integer and fractional parts if ($m =~ m|^([0-9]+)\.([0-9]*)$|) { $mi = $1; $mf = $2; } else { $mi = $m; $mf = ""; } # check sign of exponent if ($e > 0) { # positive exponent - append some digits if (length($mf) < $e) { $ri = $mi . $mf . ("0" x ($e - length($mf))); $rf = ""; } elsif (length($mf) > $e) { $ri = $mi . substr($mf, 0, $e); $rf = substr($mf, $e); } else { $ri = $mi . $mf; $rf = ""; } } elsif ($e < 0) { # negative exponent - take digits from mantissa integer and transfer over # into fraction $e = 0 - $e; if (length($mi) < $e) { $ri = "0"; $rf = ("0" x ($e - length($mi))) . $mi . $mf; } elsif (length($mi) > $e) { $ri = substr($mi, 0 - $e); $rf = substr($mi, 0, 0 - $e) . $mf; } else { $ri = "0"; $rf = $mi . $mf; } } else { # zero exponent - the mantissa is the result $ri = $mi; $rf = $mf; } if ($intonly) { $rv = $s . $ri; } else { $rv = $s . $ri . "." . $rf; } &dbg1("fbc_sci2fix: $rv", 1); return $rv; } # fbc_prim defines the primitives for BC-based flating point. $fbc_prim = <<'endquote'; # These are the full primitives sub f_int { local($a) = @_; return &fbc_fix2sci(&fbc_sci2fix($a, 1)); } # This function does what the <=> operator does for normal Perl numbers. # All the comparative routines (f_lt, f_ge, f_ne etc) call this and compare # the result to 0. sub f_cmp { local($a, $b) = @_; local($as, $am, $ae, $bs, $bm, $be); &dbg1("bc f_cmp $a <=> $b", 1); ($as, $am, $ae) = &fbc_split($a); ($bs, $bm, $be) = &fbc_split($b); $as .= "1"; $bs .= "1"; # handle the easy cases if ($as < $bs) { &dbg1("f_cmp a: -1", 1); return -1; } elsif ($as > $bs) { &dbg1("f_cmp b: 1", 1); return 1; } # zero has to be handled as a special case if (($am == 0) && ($bm == 0)) { &dbg1("f_cmp h: 0", 1); return 0; } elsif ($am == 0) { &dbg1("f_cmp f: " . (-1 * $bs), 1); return (-1 * $bs); } elsif ($bm == 0) { &dbg1("f_cmp g: $as", 1); return ($as); } # They're the same sign. From here on we treat them both as positive # but multiply our return value by the sign # Infinity is the simplest case at this point if (($ae =~ m/inf/i) && ($be =~ m/inf/i)) { &dbg1("f_cmp i: 0", 1); return 0; } elsif ($ae =~ m/inf/i) { &dbg1("f_cmp j: $as", 1); return ($as); } elsif ($be =~ m/inf/i) { &dbg1("f_cmp k: " . (-1 * $as), 1); return (-1 * $as); } # Exponents are another easy compare if ($ae < $be) { &dbg1("f_cmp c: " . (-1 * $as), 1); return (-1 * $as); } elsif ($ae > $be) { &dbg1("f_cmp d: " . $as, 1); return $as; } # Finally we have to compare the mantissas. We use a string order compare # for this. # It won't catch equals if trailing zeros differ $am =~ s/0+$//; $bm =~ s/0+$//; $as = $as * ($am cmp $bm); &dbg1("f_cmp e: " . $as, 1); return $as; } sub f_lt { local($a, $b) = @_; &dbg1("$curprim f_lt $a < $b", 1); return (&f_cmp($a, $b) < 0); } sub f_le { local($a, $b) = @_; &dbg1("$curprim f_le $a <= $b", 1); return (&f_cmp($a, $b) <= 0); } sub f_eq { local($a, $b) = @_; &dbg1("$curprim f_eq $a == $b", 1); return (&f_cmp($a, $b) == 0); } sub f_ge { local($a, $b) = @_; &dbg1("$curprim f_ge $a >= $b", 1); return (&f_cmp($a, $b) >= 0); } sub f_gt { local($a, $b) = @_; &dbg1("$curprim f_gt $a > $b", 1); return (&f_cmp($a, $b) > 0); } sub f_ne { local($a, $b) = @_; &dbg1("$curprim f_ne $a != $b", 1); return (&f_cmp($a, $b) != 0); } # Negate is about as easy as they get. The sign is just a character at the # beginning of the string. sub f_neg { local($a) = @_; &dbg1("$curprim f_neg $a", 1); if ($a =~ m/ninf/i) { $a = "inf"; } elsif ($a =~ m/inf/i) { $a = "Ninf"; } elsif ($a =~ m/^-/) { $a =~ s/^-//; } else { $a = "-" . $a; $a =~ s/-\+/-/; } &dbg1("$curprim f_neg a: $a", 1); return ($a); } # me_magcompare sub me_magcompare { local($am, $ae, $bm, $be) = @_; if ($ae > $be) { return 1; } elsif ($ae < $be) { return -1; } return ($am cmp $bm); } # m_truncround is currently just a placeholder in case I want to explicitly # do roundoff at all stages in my computations. For now, I just rely on the # fact that bc rounds for me. sub m_truncround { local($x) = @_; return $x; } # me_addpos adds two positive floating-point numbers, already split into # mantissa and exponent. sub me_addpos { local($am, $ae, $bm, $be) = @_; local($ed); # put the bigger one first if (&me_magcompare($am, $ae, $bm, $be) < 0) { ($am, $ae, $bm, $be) = ($bm, $be, $am, $ae); } # eliminate any difference in exponents $ed = $ae - $be; if ($ed > 0) { $bm =~ s/\.//; $bm = "0." . ("0" x ($ed - 1)) . $bm; $be = $ae; } # time to add $am = &bce("$am + $bm", 1); # truncate and round $am = &m_truncround($am); # append exponent and normalize $am = $am . "E" . $ae; $am = &fbc_norm($am); return $am; } # me_subpos takes two positive floating-point numbers, already split into # mantissa and exponent. The first one is known to be larger. It subtracts # the second from the first and returns the answer. sub me_subpos { local($am, $ae, $bm, $be) = @_; local($ed, $zp, $rest); # eliminate any difference in exponents $ed = $ae - $be; if ($ed > 0) { $bm =~ s/\.//; $bm = "0." . ("0" x ($ed - 1)) . $bm; $be = $ae; } # time to subtract $am = &bce("$am - $bm", 1); # truncate and round $am = &m_truncround($am); # append exponent and normalize $am = $am . "E" . $ae; $am = &fbc_norm($am); return $am; } sub f_add { local($a, $b) = @_; local($as, $am, $ae, $bs, $bm, $be); local($mc); &dbg1("$curprim f_add $a + $b", 1); ($as, $am, $ae) = &fbc_split($a); ($bs, $bm, $be) = &fbc_split($b); $as .= "1"; $bs .= "1"; # if signs are equal, call addpos. Else use subpos. if ($as == $bs) { if ($as < 0) { return &f_neg(&me_addpos($am, $ae, $bm, $be)); } else { return &me_addpos($am, $ae, $bm, $be); } } # Signs are opposite, so it's a - b or b - a. We have to figure out # which is bigger so we can put it first in the call to me_subpos $mc = &me_magcompare($am, $ae, $bm, $be); if ($mc < 0) { # a smaller than b if ($bs > 0) { # big positive, small negative; answer is positive return &me_subpos($bm, $be, $am, $ae); } else { # big one is negative, answer will be negative return &f_neg(&me_subpos($bm, $be, $am, $ae)); } } else { # b smaller than a if ($as > 0) { # big positive, small negative; answer is positive return &me_subpos($am, $ae, $bm, $be); } else { # big one is negative, answer will be negative return &f_neg(&me_subpos($am, $ae, $bm, $be)); } } } # fbc version of f_mul. sub f_mul { local($a, $b) = @_; local($as, $am, $ae, $bs, $bm, $be); &dbg1("$curprim f_mul $a * $b", 1); ($as, $am, $ae) = &fbc_split($a); ($bs, $bm, $be) = &fbc_split($b); $as .= "1"; $bs .= "1"; if (($am == 0) || ($bm == 0)) { &dbg1("$curprim f_mul a: $g_0", 1); return $g_0; } # handle infinities if (($ae > 300) || ($be > 300)) { if ($as * $bs > 0) { &dbg1("$curprim f_mul b: inf", 1); return "inf"; } &dbg1("$curprim f_mul c: Ninf", 1); return "Ninf"; } $am = &bce("$am * $bm", 1); $ae = $ae + $be; if ($as * $bs > 0) { $as = "+"; } else { $as = "-"; } $a = $as . $am . "E" . $ae; $a = &fbc_norm($a); &dbg1("$curprim f_mul d: " . $a, 1); return $a; } # This divide routine catches the divide by zero case. sub f_div { local($a, $b) = @_; local($as, $am, $ae, $bs, $bm, $be); &dbg1("$curprim f_div $a / $b", 1); ($as, $am, $ae) = &fbc_split($a); ($bs, $bm, $be) = &fbc_split($b); $as .= "1"; $bs .= "1"; # %%% should -2/0 equal Ninf? if (&f_eq($b, $g_0)) { return "inf"; } $am = &bce("$am / $bm", 1); $ae = $ae - $be; if ($as * $bs > 0) { $as = "+"; } else { $as = "-"; } $a = $as . $am . "E" . $ae; $a = &fbc_norm($a); &dbg1("$curprim f_div b: " . $a, 1); return $a; } # Our exp routine tests the infinity and overflow cases itself, just in # case the library function doesn't. sub f_exp { local($x) = @_; &dbg1("$curprim f_exp $x", 1); # Handle infinities first if ($xe =~ m/-inf/) { return $g_0; } elsif ($xe =~ m/inf/) { return $g_inf; } # Handle the extreme overflow cases. if (&f_lt($x, $g_n2000)) { return $g_0; } elsif (&f_gt($x, $g_2000)) { return $g_inf; } $x = &fbc_sci2fix($x); $x =~ s/\+//g; $x = &bce("e( $x )", 1); $x = &fbc_fix2sci($x); &dbg1("$curprim f_exp e: " . $x, 1); return $x; } # Our ln routine is strange -- it accepts negative arguments and returns # the magnitude of the answer. It also special-cases infinities and zero. sub f_ln { local($x) = @_; local($xs, $xm, $xe); &dbg1("$curprim f_ln $x", 1); ($xs, $xm, $xe) = &fbc_split($x); $xs .= "1"; # If the argument is negative we make it positive, because # we want to return the real component of the logarithm of the number, # and for any complex number, the real component of the log is the log # of the number's magnitude. if ($xs < 0) { $x = &f_neg($x); } # Handle 0 and infinities first if ($xm == 0.0) { &dbg1("$curprim f_ln a: " . $g_ninf, 1); return $g_ninf; } elsif ($xe =~ m/inf/) { &dbg1("$curprim f_ln b: " . $g_inf, 1); return $g_inf; } $x = &fbc_sci2fix($x); $x =~ s/\+//g; $x = &bce("l( $x )", 1); $x = &fbc_fix2sci($x); &dbg1("$curprim f_ln d: " . $x, 1); return $x; } # following are the partial primitives. sub f_sub { local($a, $b) = @_; &dbg1("$curprim f_sub $a - $b", 1); $a = &f_add($a, &f_neg($b)); &dbg1("$curprim f_sub: " . $a, 1); return $a; } # Following are the non-primitives sub f_pow10 { local($n) = @_; &dbg1("$curprim f_pow10 $n", 1); $n = &f_exp(&f_mul($g_ln_10, $n)); &dbg1("$curprim f_pow10: " . $n, 1); return $n; } sub f_log_n { local($n, $b) = @_; &dbg1("$curprim f_log_n v{ $b } $n", 1); return &f_div(&f_ln($n), &f_ln($b)); } sub f_log10 { local($n) = @_; &dbg1("$curprim f_log10 $n", 1); return &f_mul(&f_ln($n), $g_log10_e); } # sqrt(2 pi n) * n^n * e^(-n) * e^(1/12n) # 1/2 ln(2 pi n) + n ln(n) - n + 1/12n - 1/1260 n^3 + ... # t = n - 1; # l = 1/2 ln(2 pi) + (n + 1/2) ln(n) - n + 1/(12 n) - 1/(360 n^3) + ... # g = e^l sub f_gamma { local($n) = @_; local($l, $scal4, $n2, $np); &dbg1("$curprim f_gamma $n", 1); # For very low negative arguments, the gamma function is so # close to zero that we treat it as zero. However, for negative # integers it's infinite, so we handle that case explicitly. if (&f_lt($n, $g_n50)) { if (&f_eq($n, &f_int($n))) { return $g_inf; } return $g_0; } # Our approximation formula doesn't work well for small values, so # we use the recurrance relation gamma(n+1) = n gamma(n). $scal4 = $g_1; while (&f_lt($n, $gammaprec)) { # scale up to our current precision $scal4 = &f_mul($scal4, $n); $n = &f_add($n, $g_1); } # Since we're using Stirling's series for factorials, we have # to subtract 1 to make it be a gamma function series. $n = &f_sub($n, $g_1); $l = $g_05_ln_2pi; # %%% $l = &f_add($l, &f_mul(&f_add($n, $g_0_5), &f_ln($n))); $l = &f_sub($l, $n); $n2 = &f_mul($n, $n); $np = $n; $l = &f_add($l, &f_div($g_1, &f_mul($g_12, $np))); $np = &f_mul($np, $n2); $l = &f_sub($l, &f_div($g_1, &f_mul($g_360, $np))); $np = &f_mul($np, $n2); $l = &f_add($l, &f_div($g_1, &f_mul($g_1260, $np))); $np = &f_mul($np, $n2); $l = &f_sub($l, &f_div($g_1, &f_mul($g_1680, $np))); $np = &f_mul($np, $n2); $l = &f_add($l, &f_div($g_1, &f_mul($g_1188, $np))); $np = &f_mul($np, $n2); $l = &f_sub($l, &f_div($g_691, &f_mul($g_360360, $np))); $np = &f_mul($np, $n2); $l = &f_add($l, &f_div($g_7, &f_mul($g_1092, $np))); $np = &f_mul($np, $n2); $l = &f_sub($l, &f_div($g_3617, &f_mul($g_122400, $np))); # The last term is 3617/(122400 N^15), where N is the argument of # the factorial function being approximated. For scale=50, 27! # produces N=60, so the last term is about 6e-29. That's only enough # for 29 digits of accuracy! return (&f_div(&f_exp($l), $scal4)); } # The following method for gamma will work better: # # Determine how many digits are needed # Set x based on how many digits are needed # Determine how many terms are needed (based on z) # Gamma(z) = x^z e^(-x) SIGMA [ x^n/(z(z+1)(z+2)...(z+n)) ] # + R # where R, the "residual", is an integral that is always less than x e^(-x) # square root routine handles infinities and negative arguments properly sub f_sqrt { local($x) = @_; local($xs, $xm, $xe); &dbg1("$curprim f_sqrt $x", 1); ($xs, $xm, $xe) = &fbc_split($x); $xs .= "1"; # If the argument is negative we return 0, because 0 is # the real component of the square root of a negative number. if ($xs < 0) { return $g_0; } elsif (&f_eq($x, $g_inf)) { return $g_inf; } elsif (&f_eq($x, $g_ninf)) { return $g_0; } elsif ($xm == 0.0) { &dbg1("$curprim f_sqrt a: " . $g_0, 1); return $g_0; } $x = &fbc_sci2fix($x); $x =~ s/\+//g; $x = &bce("e(l( $x )/2)", 1); $x = &fbc_fix2sci($x); &dbg1("$curprim f_sqrt d: " . $x, 1); return $x; } # The sine algorithm uses the identity # sin(a+b) = sin a cos b + cos a sin b # The value being sin'd is broken into two pieces one of which is a multiple # of pi/16 and the other is between -pi/32 and pi/32. Then we can use the # Taylor series to get the sin of the latter part and a table for the # first part sub f_sin { local($x) = @_; # Handle arguments out of range # Note: The best way to handle things like this is to extend the # numeric representation to include an "error" term, expressing the # standard deviation of the computed value vs. the actual value. # Then, the error of the sin function is the error of the argument # times the cosine, but converging on +-1.0 when the error of the # argument is 2 pi or larger. In the limit case the computed value # will be expressed as 0.0 and the error as 1.0. # # Just for fun, I went to an office supply store and tested 5 # current-model scientific calculators (covering all price ranges) # to see what they did when you take sin or cos of larger and larger # numbers. Each calculator had a different cutoff value, but all # seemed to choose a value that was related to the calculator's # internal precision. The five cutoff values were: # 5*10^7 pi ~= 1.57*10^8 # 4.45*10^9 # 1*10^10 # 1*10^12 # 1*10^13 if (&f_gt($x, $g_sin_limit)) { return $g_0; } elsif (&f_lt($x, &f_neg($g_sin_limit))) { return $g_0; } return sin($x); } sub f_cos { local($x) = @_; # Handle out-of-range if (&f_gt($x, $g_sin_limit)) { return $g_1; } elsif (&f_lt($x, &f_neg($g_sin_limit))) { return $g_1; } return cos($x); } sub f_tan { local($x) = @_; # Handle out-of-range if (&f_gt($x, $g_sin_limit)) { return $g_inf; } elsif (&f_lt($x, &f_neg($g_sin_limit))) { return $g_inf; } return &f_div(sin($x),cos($x)); } # Arctangent function. # # Av{0} = 1 / sqrt(1 + x^{2}) , Bv{0} = 1 # Av{n+1} = (Av{n} + Bv{n}) / 2 , # Bv{n+1} = sqrt(Av{n+1} * Bv{n}) # arctan(x) = LIMIT [x * Av{0} / Av(n)] sub f_arctan { local($x) = @_; local($a, $b, $a0, $ao, $i); # A few special cases, the x==0 case is just for speed if (&f_eq($x, $g_0)) { return $g_0; } elsif (&f_eq($x, $g_inf)) { return &f_div($g_pi, $g_2); } elsif (&f_eq($x, $g_ninf)) { return &f_div($g_pi, $g_n2); } $a0 = &f_div($g_1, &f_sqrt(&f_add($g_1, &f_mul($x, $x)))); $a = $a0; $b = $g_1; # Set up ao to be different from a so the loop has a chance to start $ao = &f_sub($a, $g_1); for($i=0; ($i<25) && (&f_ne($a, $ao)); $i++) { $ao = $a; $a = &f_div(&f_add($a, $b), $g_2); $b = &f_sqrt(&f_mul($a, $b)); } return &f_div(&f_mul($x, $a0), $a); } # arcsin(x) = arctan(x / sqrt(1 - x*x)) sub f_arcsin { local($x) = @_; local($a); $a = &f_mul($x, $x); if (&f_gt($a, $g_1)) { return $g_0; } return &f_arctan(&f_div($x, &f_sqrt(&f_sub(1, $a)))); } # arccos(x) = pi/2 - arctan(x / sqrt(1 - x*x)) sub f_arccos { local($x) = @_; local($a); $a = &f_mul($x, $x); if (&f_gt($a, $g_1)) { return $g_0; } return &f_div($g_pi, $g_2) - &f_arctan(&f_div($x, &f_sqrt(&f_sub($g_1, $a)))); } # pow function is more complex than exp and pow10, because # it has to deal with any base. Specifically, if the base is negative # it computes the angle and magnitude of the complex answer # so it can calculate the real component. sub f_pow { local($b, $e) = @_; local($mag, $scal5); # Handle negative bases if (&f_lt($b, $g_0)) { # The angle of the base is pi, so the angle of the # exponent will be pi * e. Take cosine of this to # determine how much of the answer will be in the # real component. $b = &f_neg($b); $scal5 = &f_cos(&f_mul($g_pi, $e)); } else { $scal5 = $g_1; } # Compute magnitude of answer $mag = &f_exp(&f_mul(&f_ln($b), $e)); # If we got an overflow, we might be able to salvage # the computation by using logarithms. This happens when we compute # -20.0 to the power of 238.49999999999 (exactly 10 9's there) # The magnitude is 1.97e310, which overflows, but the scale # is 3.1416e-10 so the answer is actually well within # range (and even within PT0 range) if (&f_gt($mag, $g_ovr1)) { if (&f_gt($scal5, $g_0)) { $mag = &f_add(&f_mul(&f_ln($b), $e), &f_ln($scal5)); return &f_exp($mag); } else { $mag = &f_add(&f_mul(&f_ln($b), $e), &f_ln(&f_neg($scal5))); return &f_neg(&f_exp($mag)); } } else { # Return real component of answer return &f_mul($scal5, $mag); } } # Root function contains the same complexity (pun intended) as # the pow function. sub f_root { local($b, $r) = @_; local($mag, $scal6); # Handle negative bases if (&f_lt($b, $g_0)) { # The angle of the base is pi, so the angle of the # answer will be pi / r. Take cosine of this to # determine how much of the answer will be in the # real component. $b = &f_neg($b); $scal6 = &f_cos(&f_div($g_pi, $r)); } else { $scal6 = $g_1; } $mag = &f_exp(&f_div(&f_ln($b), $r)); return &f_mul($scal6, $mag); } endquote sub fbc_init { local($init_scale) = @_; if (!($bc_running)) { use IO::File; local $BCIN = new IO::File; local $BCOUT = new IO::File; use IPC::Open2; $pid = open2($BCIN, $BCOUT, 'bc -l 2>&1'); die "Can't open2 bc -l" unless ($pid); $bc_in = $BCIN; $bc_out = $BCOUT; &dbg1("open2 gave: $pid $bc_in $bc_out\n", 4); $init_scale += 2; &bce("scale=$init_scale", 0); &bce("3*4", 1); &bce("27^143", 1); $bc_running = 1; } # Define the primitives eval($fbc_prim); # These constants are mostly used by the various transcendental functions $g_0 = "0"; $g_0_5 = "0.5"; $g_1 = "1"; $g_n1 = "-1"; $g_2 = "2"; $g_n2 = "-2"; $g_4 = "4"; $g_5 = "5"; $g_6 = "6"; $g_7 = "7"; $g_8 = "8"; $g_10 = "10"; $g_12 = "12"; $g_16 = "15"; $g_16 = "16"; $g_n50 = "-50"; $g_360 = "360"; $g_691 = "691"; $g_1092 = "1092"; $g_1188 = "1188"; $g_1260 = "1260"; $g_1680 = "1680"; $g_n2000 = "-2000"; $g_2000 = "2000"; $g_3617 = "3617"; $g_122400 = "122400"; $g_360360 = "360360"; $g_2p32m1 = "4294967295"; } sub fbc_rescale { local($newscale) = @_; $newscale += 2; &bce("scale=$newscale", 0); } sub fbc_close { &dbg1("fbc_close\n", 4); close($bc_out); while(<$bc_in>) { } close($bc_in); $bc_running = 0; } #--------------------------------------------------- hybrids # # These routines operate on top of both sets of primitives, but aren't # full PT routines. They're used for conversion and formatting. # normalize works on string-format floating-point numbers of the form # [+-]1.23e[+-]456, and works over either type of primitive operators. # Its main function is to promote values like "23e456" to a PT1, and as # such it returns a PT vector as a result. sub normalize { local($x) = @_; local($m, $m1, $m2, $m3, $e, $s, $t1, $t2, $t3); if ($x =~ m|^([-+0-9.]+)E([-+0-9]+)$|) { $m = $1; $e = $2; } elsif ($x =~ m|E|) { print STDERR "normalize (s1): Cannot parse '$x'\n"; exit(-1); } else { # This is a PT or a plain real number return(0, $x); } $m1 = $m2 = ""; # Now extract the first mantissa digit and eliminate any extra zeros if (0) { } elsif ($m =~ m|^([-+]?)0*[.](0+)([1-9])([0-9]*)$|) { # 000.01, .01, 0.00123, or .00123 $s = $1; $t1 = $2; $t2 = $3; $t3 = $4; $m1 = "$s$t2"; $e = $e - (1 + length($t1)); $m2 = $t3 . "0"; } elsif ($m =~ m|^([-+]?)0*[.]([1-9])([0-9]*)$|) { # 0.1, .1, 0.123, or .123 $s = $1; $t1 = $2; $t2 = $3; $m1 = "$s$t1"; $e--; $m2 = $t2 . "0"; } elsif ($m =~ m|^([-+]?0*[1-9])[.]([0-9]*)$|) { # 01., 1., 01.23, or 1.23 $m1 = $1; $m2 = $2 . "0"; } elsif ($m =~ m|^([-+]?0*[1-9])([0-9]+)[.]([0-9]*)$|) { # 012., 12., 012.34, or 12.34 $m1 = $1; $m2 = $2; $m3 = $3; $e += length($m2); $m2 = $m2 . $m3; } elsif ($m =~ m|^([-+]?0*[1-9])([0-9]*)$|) { # 01, 1, 012, or 12 $m1 = $1; $m2 = $2; $e += length($m2); $m2 = $m2 . "0"; } else { print STDERR "normalize (s2): Cannot parse '$x'\n"; exit(-1); } # Rebuild the mantissa in a normal 1.23 format $m = "$m1.$m2"; $x = $m . "E" . $e; if ($e < -300) { return(0, $g_0); } elsif ($e > 300) { # Need to promote. First, extract sign if ($m < 0) { $s = -1; $m = &f_neg($m); } else { $s = 1; } # take log. Note, we will call whichever f_xxx routines are # current and that's the right thing to do. $e = &f_add($e, &f_log10($m)); # put sign back in $e = &f_mul($s, $e); return (1, $e); } # else # no promote necessary return (0, $x); } # split takes string numbers of the format "2P[+-]3.4E[+-]56", or simple # floating-point numbers like "[+-]3.4E[+-]56", or IR inlines like "I17", # performs the appropriate conversion or value retrieval and returns a PT # vector result. sub split { local($x) = @_; local($x_pt, $x_val, $i); local($xv_pt, $xv_val); # for normalizig, if string is for example # "1.2e345" or "2pt1.2e345" if ($x =~ m|^I([0-9]+)|) { # retrieve a stored intermediate result $i = $1; return($ir_pt[$i], $ir_val[$i]); } elsif ($x =~ m|^R([0-9]+)|) { # retrieve a stored history result $i = $1; return($h_pt[$i], $h_val[$i]); } elsif ($x =~ m|^(.+)P(.+)$|) { # convert an inline. # first, assume $x_pt = $1; $x_val = $2; $x_val =~ s|N|-|g; ($xv_pt, $xv_val) = &normalize($x_val); $x_pt += $xv_pt; $x_val = $xv_val; return ($x_pt, $x_val); } # else, this is a normal number "1.23" or "-23.45", or scientific notation # "1.23E45" print "split1 $x : " if ($format_debug); $x =~ s|N|-|g; ($xv_pt, $xv_val) = &normalize($x); print "$xv_pt $xv_val\n" if ($format_debug); return ($xv_pt, $xv_val); } #--------------------------------------------------- hcfloat.h # # This is the data type used for user calculations. The first field "pt" # is an integer representing how many times we've taken the logarithm # (base 10) of the value being represented. The "val" field is # the value, or the log of the value, or the log of the log of the # value, etc. If the value being represented is negative, the "val" # field is negative. # # typedef struct { # unsigned int pt; # double val; # } pt_real; #------------------------------------------------------ hcpt.c # # These routines, originally from "hcpt.c", do all the calculation # functions and operate on "PT" format values. # Addition and subtraction are done through two intermediate routines, # addpos and subpos. addpos adds any two positive values (not necessarily # with the largest first) and subpos subtracts a positive value from another # positive value, where the result is known to be positive or zero. These # two routines combined with appropriate sign manipulation cover all the # cases of addition and subtraction. sub pt_addpos { local($x_pt, $x_val, $y_pt, $y_val) = @_; local($xpt, $xv, $ypt, $yv); local($s1); local($ans_pt, $ans_val); # Put bigger argument first, so we have fewer cases to handle if ($x_pt > $y_pt) { $xpt = $x_pt; $ypt = y_pt; $xv = $x_val; $yv = $y_val; } elsif (($x_pt == $y_pt) && &f_gt($x_val, $y_val)) { $xpt = $x_pt; $ypt = $y_pt; $xv = $x_val; $yv = $y_val; } else { $xpt = $y_pt; $ypt = $x_pt; $xv = $y_val; $yv = $x_val; } # Deal with infinities. if (&f_eq($xv, $g_inf)) { return ($g_pt_inf, $g_inf); } if ($xpt == 0) { # both are simple numbers, add and test for promote $s1 = &f_add($xv, $yv); if (&f_gt($s1, $g_ovr1)) { $ans_pt = 1; $ans_val = &f_log10($s1); return ($ans_pt, $ans_val); } else { $ans_pt = 0; $ans_val = $s1; return ($ans_pt, $ans_val); } } elsif ($xpt == 1) { # Adding a PT0 to a PT1, answer will usually be x, except in cases like # {1, 301.0} + (0, 1.0e299} => {1, 301.004321} # # Adding two PT1's is similar, and usually reduces to just returning X. # The only difference is that we don't have to take the log of y. # # To compute the log of the answer, we take the log of the two inputs # and subtract, then take 10 to the power of that value, take the # antilog, add to 1 and take the log again: # inputs: {1, 301.0} + {0, 1.0e299} # logs: 301.0, 299.0 # subtract: -2.0 # antilog: 0.01 # add to 1: 1.01 # log: 0.004321 # add: 301.004321 # final answer: {1, 301.004321} $ans_pt = 1; if ($ypt == 0) { $ans_val = &f_log10($yv); } else { $ans_val = $yv; } $ans_val = &f_sub($ans_val, $xv); $ans_val = &f_add($g_1, &f_pow10($ans_val)); $ans_val = &f_add($xv, &f_log10($ans_val)); return ($ans_pt, $ans_val); } else { # X is PT2 or larger. Answer is always uncomputably larger than X. return ($xpt, $xv); } $x_val = &f_add($x_val, $y_val); return ($x_pt, $x_val); } # pt_negate is trivial, but is implemented as a separate routine to try to # simplify the work that would be entailed if I change the method of # representing sign in PT's. sub pt_negate { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); $ans_pt = $x_pt; $ans_val = &f_neg($x_val); return ($ans_pt, $ans_val); } sub pt_subpos { local($x_pt, $x_val, $m_pt, $m_val) = @_; local($xpt, $xv, $mpt, $mv, $s1); local($ans_pt, $ans_val); $xpt = $x_pt; $mpt = $m_pt; $xv = $x_val; $mv = $m_val; # Deal with infinities. If it's +inf + -inf, +inf wins. if (&f_eq($xv, $g_inf)) { return ($g_pt_inf, $g_inf); } if (($xpt == $mpt) && (&f_eq($xv, $mv))) { return (0, $g_0); } if (($mpt == 0) && (&f_eq($mv, $g_0))) { return ($x_pt, $x_val); } if ($xpt == 0) { # both are simple numbers, subtract -- no demote is possible $s1 = &f_sub($xv, $mv); return (0, $s1); } elsif (($xpt == 1) && ($mpt == 0)) { # Subtracting a PT0 from a PT1. Answer will usually be x, but in a few cases like # {1, 301.0} - {0, 1.0e299} => {1, 300.995678} # To compute the log of the answer, we take the log of the two inputs and subtract, # then take 10 to the power of that value, take the antilog, add to 1 and take the log again # inputs: {1, 301.0} - {0, 1.0e299} # logs: 301.0, 299.0 # subtract: -2.0 # antilog: 0.01 # sub from 1: 0.99 # log: -0.004321 # add: 300.995678 # final answer: {1, 300.995678} # In addition, we check for demote to a PT0 answer. $s1 = &f_sub(&f_log10($mv), $xv); $s1 = &f_sub($g_1, &f_pow10($s1)); $s1 = &f_add($xv, &f_log10($s1)); if (&f_gt($s1, $g_log10_ovr1)) { return (1, $s1); } # else return (0, &f_pow10($s1)); } elsif ($xpt == 1) { # This is the case where both are PT1's. We actually do the same thing as # in the previous case (although it usually reduces to just returning X) # the only difference being that we don't have to take the log of y. $s1 = &f_sub($mv, $xv); $s1 = &f_sub($g_1, &f_pow10($s1)); $s1 = &f_add($xv, &f_log10($s1)); if (&f_gt($s1, $g_log10_ovr1)) { return (1, $s1); } # else return(0, &f_pow10($s1)); } else { # X is PT2 or larger. No matter what Y is, the answer is X because X is bigger. return ($xpt, $xv); } } q@ # Compare routine -- right now it isn't used anywhere. int pt_compare(pt_real x, pt_real y) { int neg; if (($x_val < 0) && ($y_val > 0)) { return -1; } elsif (($x_val > 0) && ($y_val < 0)) { return 1; } # Both are of same sign. Force to both positive if ($x_val < 0) { $x_val = - $x_val; $y_val = - $y_val; neg = -1; } else { neg = 1; } # check difference in PTs if ($x_pt > $y_pt) { return neg * 1; } elsif ($x_pt < $y_pt) { return neg * -1; } # check difference in values if ($x_val > $y_val) { return neg * 1; } elsif ($x_val < $y_val) { return neg * -1; } # If they made it through all this, they're equal! return 0; } @; # Magnitude compare -- tells you which is "larger" as in further from 0. # Use pt_compare for standard compare. sub pt_magcompare { local($x_pt, $x_val, $y_pt, $y_val) = @_; local($xv, $yv); $xv = $x_val; $yv = $y_val; if (&f_lt($xv, $g_0)) { $xv = &f_neg($xv); } if (&f_lt($yv, $g_0)) { $yv = &f_neg($yv); } # check difference in PTs if ($x_pt > $y_pt) { return 1; } elsif ($x_pt < $y_pt) { return -1; } # check difference in values if (&f_gt($xv, $yv)) { return 1; } elsif (&f_lt($xv, $yv)) { return -1; } # If they made it through all this, they're equal! return 0; } # Now, finally, we can do addition. There are four cases of positive and negative # arguments, and in addition if the signs are different there are two cases # corresponding to which is of greater magnitude. This routine dispatches each # of the 6 cases using the two compare routines to figure out which case is # present and the two arithmetic routines to generate answers, with pt_negate # to massage the signs of arguments and results where necessary. sub pt_add { local($x_pt, $x_val, $y_pt, $y_val) = @_; # Deal with infinities. If it's +inf and -inf, +inf wins. if (&f_eq($x_val, $g_inf)) { return ($g_pt_inf, $g_inf); } elsif (&f_eq($y_val, $g_inf)) { return ($g_pt_inf, $g_inf); } elsif (&f_eq($x_val, $g_ninf)) { return ($g_pt_inf, $g_ninf); } elsif (&f_eq($y_val, $g_ninf)) { return ($g_pt_inf, $g_ninf); } if (&f_ge($x_val, $g_0)) { if (&f_ge($y_val, $g_0)) { # Both positive return &pt_addpos($x_pt, $x_val, $y_pt, $y_val); } else { # +X, -Y if (&pt_magcompare($x_pt, $x_val, $y_pt, $y_val) > 0) { # X is of greater magnitude return &pt_subpos($x_pt, $x_val, &pt_negate($y_pt, $y_val)); } else { return &pt_negate(&pt_subpos(&pt_negate($y_pt, $y_val), $x_pt, $x_val)); } } } elsif (&f_ge($y_val, $g_0)) { # -X, +Y if (&pt_magcompare($x_pt, $x_val, $y_pt, $y_val) > 0) { # X is of greater magnitude return &pt_negate(&pt_subpos(&pt_negate($x_pt, $x_val), $y_pt, $y_val)); } else { return &pt_subpos($y_pt, $y_val, &pt_negate($x_pt, $x_val)); } } else { # Both negative return &pt_negate(&pt_addpos(&pt_negate($x_pt, $x_val),&pt_negate($y_pt, $y_val))); } } # After going to all that work to do addition, subtraction is easy! sub pt_sub { local($x_pt, $x_val, $m_pt, $m_val) = @_; return &pt_add($x_pt, $x_val, &pt_negate($m_pt, $m_val)); } # We define PT versions of pow10 and log10 next, because # the multiply and divide routines use them to reduce certain cases to adds. # pt_log10 is very simple, because it's a single function argument and # there are no strange boundary conditions where we demote or promote. sub pt_log10 { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); # Handle 0 and infinities first if (&f_eq($x_val, $g_0)) { return ($g_pt_inf, $g_ninf); } elsif (&f_eq($x_val, $g_inf)) { return ($g_pt_inf, $g_inf); } elsif (&f_eq($x_val, $g_ninf)) { return (0, $g_0); } if (&f_lt($x_val, $g_0)) { # Real component of log of negative number is the log of the magnitude $x_val = &f_neg($x_val); } if ($x_pt == 0) { # Call normal log10 function $ans_pt = 0; $ans_val = &f_log10($x_val); return ($ans_pt, $ans_val); } # else # Just subtract 1 from PT $ans_pt = $x_pt - 1; $ans_val = $x_val; return ($ans_pt, $ans_val); } # pt_pow10 is almost as simple as pt_log10, and for similar reasons. sub pt_pow10 { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); # Handle infinities first if (&f_eq($x_val, $g_inf)) { return ($g_pt_inf, $g_inf); } elsif (&f_eq($x_val, $g_ninf)) { return (0, $g_0); } if ($x_pt == 0) { # Handle all negative PT0's and any positive PT0's for which # the answer is in PT0 range. if (&f_lt($x_val, $g_log10_ovr1)) { $ans_pt = 0; $ans_val = &f_pow10($x_val); return ($ans_pt, $ans_val); } } # Now handle negative PT1's, etc. We only have to test $x_val here, # because the negative PT0 case was already handled. if (&f_lt($x_val, $g_0)) { return (0, $g_0); } # All other cases: simple promote $ans_pt = $x_pt + 1; $ans_val = $x_val; return ($ans_pt, $ans_val); } # The multiply routine is implemented fully from primitives, with # 12 separate cases for all the combinations of PT class of the # two arguments and the result. It's this way because the # multiply was implemented before the add, and it is left as-is # for accuracy and speed. sub pt_mul { local($px_pt, $px_val, $py_pt, $py_val) = @_; local($x_pt, $x_val, $y_pt, $y_val); local($xpt, $ypt, $yval); local($p1, $d1, $sign); local($ans_pt, $ans_val); $x_pt = $px_pt; $y_pt = $py_pt; $x_val = $px_val; $y_val = $py_val; # Enforce positive arguments and remember sign. $sign = $g_1; if (&f_lt($x_val, $g_0)) { $sign = &f_neg($sign); $x_val = &f_neg($x_val); } if (&f_lt($y_val, $g_0)) { $sign = &f_neg($sign); $y_val = &f_neg($y_val); } # Put bigger argument first, so we have fewer cases to handle if (&pt_magcompare($x_pt, $x_val, $y_pt, $y_val) > 0) { $xpt = $x_pt; $ypt = $y_pt; $xv = $x_val; $yv = $y_val; } else { $xpt = $y_pt; $ypt = $x_pt; $xv = $y_val; $yv = $x_val; } if (&f_eq($yv, $g_0)) { # don't need to test $ypt because $yv==0 only happens when $ypt==0 # trivial case return (0, $g_0); } elsif (&f_eq($xv, $g_inf)) { # don't need to test pt field, val==inf only occurs when pt==MAXINT if (&f_gt($sign, 0)) { return ($g_pt_inf, $g_inf); } else { return ($g_pt_inf, $g_ninf); } } elsif ($xpt == 0) { # %%% it would be good to avoid the overflow before it happens # Simple case: both are plain numbers $p1 = &f_mul($xv, $yv); if (&f_ge($p1, $g_ovr1)) { # Overflowed range of pt0: recompute product as a pt1 $p1 = &f_add(&f_log10($xv), &f_log10($yv)); $ans_pt = 1; $ans_val = &f_mul($sign, $p1); return ($ans_pt, $ans_val); } else { $ans_pt = 0; $ans_val = &f_mul($sign, $p1); return ($ans_pt, $ans_val); } } elsif (($xpt == 1) && ($ypt == 0)) { # X is PT1 and Y is PT0. # We don't worry about promoting to # pt2 because that can only happen if X is near {1, 1e300} # and if it were, we're calculating log(answer) = 1e100 # plus log(y), and since log(y) can be no bigger than 300 # the answer would always be 1e300. # However, we do have to worry about demoting to PT0. $p1 = &f_add($xv, &f_log10($yv)); if (&f_gt($p1, $g_log10_ovr1)) { $ans_pt = 1; $ans_val = &f_mul($sign, $p1); return ($ans_pt, $ans_val); } else { $ans_pt = 0; $ans_val = &f_mul($sign, &f_pow10($p1)); return ($ans_pt, $ans_val); } } elsif ($xpt == 1) { # Both are PT1's -- just add the logs. $d1 = &f_add($xv, $yv); if (&f_gt($d1, $g_ovr1)) { # Promote to PT2 $ans_pt = 2; $ans_val = &f_mul($sign, &f_log10($d1)); return ($ans_pt, $ans_val); } elsif (&f_gt($d1, $g_log10_ovr1)) { $ans_pt = 1; $ans_val = &f_mul($sign, $d1); return ($ans_pt, $ans_val); } else { # Demote to PT0 $ans_pt = 0; $ans_val = &f_mul($sign, &f_pow10($d1)); return ($ans_pt, $ans_val); } } elsif (($xpt == 2) && ($ypt < 2)) { # X is PT2 and Y is PT0 or PT1. The answer is always equal to the larger, which is x. $ans_pt = $xpt; $ans_val = &f_mul($sign, $xv); return ($ans_pt, $ans_val); } elsif ($xpt == 2) { # Both are PT2's. In this case, the product is equal to the larger # except in cases like {2, 312.0} * {2, 308.0} where the two # values are within +- scale or so, where we have to do the # log(sum(exp(diff))) thing. # inputs: {2, 312.0} * {2, 308.0} # logs: 312.0, 308.0 # subtract: -4.0 # antilog: 0.0001 = 1.0e-4 # add to 1: 1.0001 # log: 0.000043 # add: 312.000043 # final answer: {2, 312.000043} # $d1 = &f_sub($xv, $yv); if (&f_lt($d1, $curscale)) { # x is a little bigger $ans_pt = 2; $ans_val = &f_pow10(&f_sub($yv, $xv)); $ans_val = &f_add($g_1, $ans_val); $ans_val = &f_add($xv, &f_log10($ans_val)); $ans_val = &f_mul($sign, $ans_val); return ($ans_pt, $ans_val); } else { # x is way bigger $ans_pt = 2; $ans_val = &f_mul($sign, $xv); return ($ans_pt, $ans_val); } } elsif ($xpt > 2) { # For anything with X bigger than a PT2, the product is the # larger of the two. $ans_pt = $xpt; $ans_val = &f_mul($sign, $xv); return ($ans_pt, $ans_val); } # default fall-through return (0, $g_0); } sub pt_div(pt_real an, pt_real ad) { local($an_pt, $an_val, $ad_pt, $ad_val) = @_; local($n_pt, $n_val, $d_pt, $d_val); local($ans_pt, $ans_val); local($sign, $p1); $n_pt = $an_pt; $d_pt = $ad_pt; $n_val = $an_val; $d_val = $ad_val; # Handle divide by zero case first, it's easy if ($d_val == 0) { return ($g_pt_inf, $g_inf); } # Enforce positive arguments and remember sign. $sign = $g_1; if (&f_lt($n_val, $g_0)) { $sign = &f_neg($sign); $n_val = &f_neg($n_val); } if (&f_lt($d_val, $g_0)) { $sign = &f_neg($sign); $d_val = &f_neg($d_val); } # Deal with infinities. if (&f_eq($n_val, $g_inf)) { if (&f_gt($sign, $g_0)) { return ($g_pt_inf, $g_inf); } else { return ($g_pt_inf, $g_ninf); } } if (($n_pt == 0) && ($d_pt == 0)) { # Simple case: both are plain numbers. Divide and check # for promote (which happens when dividing 1e200 by 1e-200, # for example) $p1 = &f_div($n_val, $d_val); if (&f_ge($p1, $g_ovr1)) { # Overflowed range of pt0: promote to pt1 $p1 = &f_sub(&f_log10($n_val), &f_log10($d_val)); $ans_pt = 1; $ans_val = &f_mul($sign, $p1); return ($ans_pt, $ans_val); } else { $ans_pt = 0; $ans_val = &f_mul($sign, $p1); return ($ans_pt, $ans_val); } } elsif (($n_pt == 1) && ($d_pt == 0)) { # Dividing a PT1 by a PT0. Need to check for demote, but promote # isn't possible. $p1 = &f_sub($n_val, &f_log10($d_val)); if (&f_lt($p1, $g_log10_ovr1)) { $ans_pt = 0; $ans_val = &f_mul($sign, &f_pow10($p1)); return ($ans_pt, $ans_val); } else { $ans_pt = 1; $ans_val = &f_mul($sign, $p1); return ($ans_pt, $ans_val); } } elsif (($n_pt == 0) && ($d_pt == 1)) { # Dividing a PT0 by a PT1. If the PT0 is something like 1e200 # and the PT1 is say 1e305, the answer is 1e-105 which is a # valid nonzero PT0. So, we have to compute an answer. $p1 = &f_sub(&f_log10($n_val), $d_val); $ans_pt = 0; $ans_val = &f_mul($sign, &f_pow10($p1)); return ($ans_pt, $ans_val); } elsif (($n_pt == 1) && ($d_pt == 1)) { # Both arguments are PT1. Subtract the logs and check for demote. $p1 = &f_sub($n_val, $d_val); if (&f_lt($p1, $g_log10_ovr1)) { $ans_pt = 0; $ans_val = &f_mul($sign, &f_pow10($p1)); return ($ans_pt, $ans_val); } else { $ans_pt = 1; $ans_val = &f_mul($sign, $p1); return ($ans_pt, $ans_val); } } elsif ($n_pt < $d_pt) { # In all the other cases, if the arguments are of dirrerent # PT types the answer is simple. return (0, $g_0); } elsif ($n_pt > $d_pt) { # In these cases the answer is the numerator. $ans_pt = $n_pt; $ans_val = &f_mul($sign, $n_val); return ($ans_pt, $ans_val); } else { # Numerator and denominator of same PT type, PT2 or greater. # Ratio uncomputably differs from: # 0 if D is of greater magnitude, # N if N is of greater magnitude # 1.0 or -1.0 if they are equal. if (&f_lt($n_val, $d_val)) { return (0, $g_0); } elsif (&f_eq($n_val, $d_val)) { $ans_pt = 0; $ans_val = $sign; return ($ans_pt, $ans_val); } else { $ans_pt = $n_pt; $ans_val = &f_mul($sign, $n_val); return ($ans_pt, $ans_val); } } } q% pt_real pt_percent(pt_real p, pt_real x) { local($ans_pt, $ans_val); pt_real ans; $ans_pt = 0; $ans_val = 100.0; ans = &pt_div(p, ans); ans = &pt_mul(ans, x); return ($ans_pt, $ans_val); } %; sub pt_exp { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); local($l1); # Handle infinities first if (&f_eq($x_val, $g_inf)) { return ($g_pt_inf, $g_inf); } elsif (&f_eq($x_val, $g_ninf)) { return (0, $g_0); } if ($x_pt == 0) { # There are three cases; the third (involving promotion to PT2) # only happens when we're called from pt_gamma with an argument # like {0, 7e302}, which is about as high as it will ever get. if (&f_lt($x_val, $g_ln_ovr1)) { # Simple case: answer is within PT0 range $ans_pt = 0; $ans_val = &f_exp($x_val); return ($ans_pt, $ans_val); } elsif (&f_lt(&f_div($x_val, $g_ln_10), $g_ovr1)) { # Promote to PT1 $ans_pt = 1; $ans_val = &f_div($x_val, $g_ln_10); return ($ans_pt, $ans_val); } else { # Promote to PT2 $ans_pt = 2; $ans_val = &f_log10(&f_div($x_val, $g_ln_10)); return ($ans_pt, $ans_val); } } # Now handle the rest of the negative cases if (f_lt($x_val, $g_0)) { return (0, $g_0); } # Now handle positive PT1 if ($x_pt == 1) { # e^(10^X) = 10^(log10(e) * 10^X) = 10^(10^(log10(log10(e)) + X )) # X is from: # 10^300 -- answer is 10^(10^(300+logloge)) = 10^(10^299.637) = 10^(4.342e299) # to # 10^(1e300) -- answer is 10^(10^(10^300+logloge)). # Since log10(log10(e)) is negative, we sometimes have to produce # a PT1 answer. $l1 = &f_add($x_val, $g_log10_log10_e); if (&f_lt($l1, $g_log10_ovr1)) { $ans_pt = 1; $ans_val = &f_mul($g_log10_e, &f_pow10($x_val)); return ($ans_pt, $ans_val); } else { return (2, $l1); } } else { # X is PT2 or greater. e^X is uncomputably larger than 10^x, so we # just promote to the next higher PT with the same val. return ($x_pt + 1, $x_val); } # default fall-through (never reached) return (0, $g_0); } sub pt_ln { local($x_pt, $x_val) = @_; local($l1); local($ans_pt, $ans_val); # Handle 0 and infinities first if (&f_eq($x_val, $g_0)) { return ($g_pt_inf, $g_ninf); } elsif (&f_eq($x_val, $g_inf)) { return ($g_pt_inf, $g_inf); } elsif (&f_eq($x_val, $g_ninf)) { return (0, $g_0); } if (&f_lt($x_val, $g_0)) { # Real component of log of negative number is the log of the magnitude $x_val = &f_neg($x_val); } if ($x_pt == 0) { # Call normal ln function $ans_pt = 0; $ans_val = &f_ln($x_val); return ($ans_pt, $ans_val); } elsif ($x_pt == 1) { # ln(x) = log10(x) / log10(e) = $x_val / log10(e) = $x_val * $g_ln_10 $l1 = &f_mul($x_val, $g_ln_10); if (&f_lt($l1, $g_ovr1)) { return (0, $l1); } else { $ans_pt = 1; $ans_val = &f_log10($l1); return ($ans_pt, $ans_val); } } elsif ($x_pt == 2) { # Have to subtract log10(log10(e) $l1 = &f_sub($x_val, $g_log10_log10_e); if (&f_lt($l1, $g_log10_ovr1)) { $ans_pt = 0; $ans_val = &f_pow10($l1); return ($ans_pt, $ans_val); } else { return (1, $l1); } } else { $ans_pt = $x_pt - 1; $ans_val = $x_val; return ($ans_pt, $ans_val); } } # A routine that breaks out all the cases explicitly would be a little # more accurate and a little faster -- but it's not enough to be of # immediate concern for the usability of the calculator. sub pt_pow { local($ab_pt, $ab_val, $ae_pt, $ae_val) = @_; local($b_pt, $b_val, $e_pt, $e_val); local($p1, $scal7); local($ans_pt, $ans_val); $b_pt = $ab_pt; $e_pt = $ae_pt; $b_val = $ab_val; $e_val = $ae_val; # Handle 0's, 1's, and infinities first. Some of these are handled # quickly as a speed optimization, and others because they're special # cases. # Here is a table illustrating the special cases; a total of 19 # special cases are handled. Any entry that says "." is computed # the normal way, but the others are special cases. For the "C*i" # cases we actually compute the cosine of the exponent to determine # the sign of infinity that should result. The rows/columns labeled # -5, -0.5, 0.5 and 5 represent any typical numbers (there is nothing # special about -5 -0.5, 0.5 or 5) # +--------------EXPONENT--------------- # pow |ninf -5 -1 -0.5 0 0.5 1 5 inf # +-----+------------------------------------- # | inf| 0 0 0 0 1 inf inf inf inf # | 5| 0 . . . 1 . 5 . inf # B 1| 1 1 1 1 1 1 1 1 1 # A 0.5| inf . . . 1 . 0.5 . 0 # S 0| inf inf inf inf 1 0 0 0 0 # E -0.5| inf . . . 1 . -0.5 . 0 # | -1| 1 . . . 1 . -1 . 1 # | -5| 0 . . . 1 . -5 . inf # | ninf| 0 0 0 0 1 C*i ninf C*i inf if (&f_eq($e_val, $g_0)) { return (0, 1); } elsif (&f_eq($b_val, $g_1)) { return (0, 1); } elsif (&f_eq($e_val, $g_1)) { return ($ab_pt, $ab_val); } elsif (&f_eq($b_val, $g_inf)) { # Don't have to check exponent == 0 because it's already been checked if (&f_gt($e_val, $g_0)) { return ($g_pt_inf, $g_inf); } else { return (0, $g_0); } } elsif (&f_eq($b_val, $g_0)) { # Don't have to check exponent == 0 because it's already been checked if (&f_gt($e_val, $g_0)) { return (0, $g_0); } else { return ($g_pt_inf, $g_inf); } } $p1 = &f_mul($b_val, $b_val); if (&f_eq($e_val, $g_inf)) { if (&f_gt($p1, $g_1)) { return ($g_pt_inf, $g_inf); } elsif (&f_eq($p1, $g_1)) { return (0, $g_1); } else { return (0, $g_0); } } elsif (&f_eq($e_val, $g_ninf)) { if (&f_gt($p1, $g_1)) { return (0, $g_0); } elsif (&f_eq($p1, $g_1)) { return (0, $g_1); } else { return ($g_pt_inf, $g_inf); } } elsif (&f_eq($b_val, $g_ninf)) { if (&f_lt($e_val, $g_0)) { return (0, $g_0); } elsif (&f_eq($e_val, $g_1)) { return ($g_pt_inf, $g_ninf); } elsif (&f_eq($e_val, $g_inf)) { return ($g_pt_inf, $g_inf); } $p1 = &f_cos(&f_mul($e_val, $g_pi)); if (&f_gt($p1, $g_0)) { return ($g_pt_inf, $g_inf); } elsif (&f_eq($p1, $g_0)) { return (0, $g_0); } else { return ($g_pt_inf, $g_ninf); } } # extract sign and scale information if (&f_lt($b_val, $g_0)) { if ($e_pt == 0) { # Compute the real component of the angle $scal7 = &f_cos(&f_mul($e_val, $g_pi)); } else { # It'll round off to a multiple of 2pi, so just scale # by 1. $scal7 = $g_1; } $b_val = &f_neg($b_val); } else { $scal7 = $g_1; } # Now $b_val is positive, so it's okay to take the log of it if we need to. if (($b_pt == 0) && ($e_pt == 0)) { $p1 = f_pow($b_val, $e_val); if (&f_lt($p1, $g_ovr1)) { $ans_pt = 0; $ans_val = &f_mul($scal7, $p1); return ($ans_pt, $ans_val); } # Overflowed on direct computation, so compute it by logs instead. $p1 = &f_mul(&f_log10($b_val), $e_val); if (&f_lt($p1, $g_ovr1)) { $ans_pt = 1; # %%% bug here: if scale is small it can bring answer back into PT0 range if (&f_eq($scal7, $g_1)) { $ans_val = $p1; } elsif (&f_gt($scal7, $g_0)) { $ans_val = &f_add($p1, &f_log10($scal7)); } else { $ans_val = &f_neg(&f_add($p1, &f_log10(&f_neg($scal7)))); } return ($ans_pt, $ans_val); } else { # This case happens when raising e.g. {0, 1e10} to the power of {0, 5e299} # The answer is {1, 5e300} which must be promoted to {2, 300.698} # The contribution of scale is uncomputable, except for its sign. $ans_pt = 2; $p1 = &f_log10($p1); if (&f_gt($scal7, $g_0)) { $ans_val = $p1; } else { $ans_val = &f_neg($p1); } return ($ans_pt, $ans_val); } } # For the remaining cases, it seems that we should be able to # just use powers and logs. ($ans_pt, $ans_val) = &pt_pow10( &pt_mul(&pt_log10($b_pt, $b_val), $e_pt, $e_val)); # Factor in the scale. # So I can take -1e2000 to the power of 0.25 and it will get 7.07e500 $b_pt = 0; $b_val = $scal7; ($ans_pt, $ans_val) = &pt_mul($ans_pt, $ans_val, $b_pt, $b_val); return ($ans_pt, $ans_val); } # There is a lot of similarity with the pt_pow routine because pt_root # is just dividing the log by the exponent rather than multiplying. # # r-th root of b = pow(b, 1.0 / r) # = exp(ln(b) / r) # = exp(exp(ln(ln(b)) - ln(r))) sub pt_root { local($ab_pt, $ab_val, $ar_pt, $ar_val) = @_; local($ans_pt, $ans_val); local($b_pt, $b_val, $r_pt, $r_val); local($p1, $scal8); ($b_pt, $b_val, $r_pt, $r_val) = ($ab_pt, $ab_val, $ar_pt, $ar_val); # Deal with 0th root, which is like raising to the power of infinity if (&f_eq($r_val, $g_0)) { if (&f_eq($b_val, $g_0)) { # 0th root of 0 is 0 return (0, $g_0); } elsif (&f_gt($b_val, $g_n1) && &f_lt($b_val, $g_1)) { # 0th root of a number less than 1 return (0, $g_1); } else { return ($g_pt_inf, $g_inf); } } # extract sign and scale information if (&f_lt($b_val, $g_0)) { if ($r_pt == 0) { # We can calculate sin $scal8 = &f_cos(&f_div($g_pi, $r_val)); } else { # We're raising to a power that's very very close to 0, # so close that the cosine is essentially 1.0 $scal8 = $g_1; } $b_val = &f_neg($b_val); } else { $scal8 = $g_1; } # Now $b_val is non-negative, so it's okay to take the log of it if we need to. if (($b_pt == 0) && ($r_pt == 0)) { # Compute the root directly. This will overflow if we're taking say # the 0.1th root of 1e100, which is like taking 1e100 to the 10th # power. $p1 = &f_root($b_val, $r_val); if (&f_lt($p1, $g_ovr1)) { $ans_pt = 0; $ans_val = &f_mul($scal8, $p1); return ($ans_pt, $ans_val); } # Overflowed on direct computation, so compute it by logs instead. $p1 = &f_div(&f_log10($b_val), $r_val); if (&f_lt($p1, $g_ovr1)) { $ans_pt = 1; # %%% bug here: if scale is small it can bring answer back into PT0 range if (&f_eq($scal8, 1.0)) { $ans_val = $p1; } elsif (&f_gt($scal8, $g_0)) { $ans_val = &f_add($p1, &f_log10($scal8)); } else { $ans_val = &f_neg(&f_add($p1, &f_log10(&f_neg($scal8)))); } return ($ans_pt, $ans_val); } else { # This case happens when taking e.g. the {0, 2e-300} root of {0, 1e10}, # which is equivalent to raising {0, 1e10} to the power of {0, 5e299} # The answer is {1, 5e300} which must be promoted to {2, 300.698} # The contribution of scale is uncomputable, except for its sign. $ans_pt = 2; $p1 = &f_log10($p1); if (&f_gt($scal8, $g_0)) { $ans_val = $p1; } else { $ans_val = &f_neg($p1); } return ($ans_pt, $ans_val); } } # For the remaining cases, it seems that we should be able to # just use powers and logs. ($ans_pt, $ans_val) = &pt_pow10( &pt_div(&pt_log10($b_pt, $b_val), $r_pt, $r_val)); # Factor in the scale. # So I can take the 4th root of -1e2000 and it will get 7.07e500 ($ans_pt, $ans_val) = &pt_mul($ans_pt, $ans_val, 0, $scal8); return ($ans_pt, $ans_val); } # log base b of n = log(n) / log(b) # == pow10(log(log(n)) - log(log(b)) # # By the way, this function is a better match for being the analogue of the # divide function because of the symmetry of the formula for computing it. sub pt_log_n { local($n_pt, $n_val, $b_pt, $b_val) = @_; local($ans_pt, $ans_val); # Handle 0 and infinities first if (&f_eq($n_val, $g_0)) { return ($g_pt_inf, $g_ninf); } elsif (&f_eq($n_val, $g_inf)) { return ($g_pt_inf, $g_inf); } elsif (&f_eq($n_val, $g_ninf)) { return (0, $g_0); } # Real component of log of a negative number is the log of the magnitude if (&f_lt($n_val, $g_0)) { $n_val = &f_neg($n_val); } if ($n_pt == 0) { if ($b_pt == 0) { # Both arguments are PT0. The biggest the answer can possibly be is # somewhere around 1.0e19 (a wild guess based on just playing around # on the calculator), definitely well within PT0 bounds. $ans_pt = 0; $ans_val = &f_log_n($n_val, $b_val); return ($ans_pt, $ans_val); } } # Otherwise we just defer to the primitives. # %%% This needs to be changed; right now it only works if b is greater than 0 return &pt_div(&pt_log10($n_pt, $n_val), &pt_log10($b_pt, $b_val)); } # t = n - 1 # l = 1/2 ln(2 pi) + (t + 1/2) ln(t) - t + 1/(12 t) - 1/(360 t^3) + ... # gamma(x) = e^l sub pt_gamma { local($n_pt, $n_val) = @_; local($ans_pt, $ans_val); local($t, $l, $adj, $n2, $np); # Handle infinities first if (&f_eq($n_val, $g_inf)) { return ($g_pt_inf, $g_inf); } elsif (&f_eq($n_val, $g_ninf)) { return (0, $g_0); } # Now handle negatives if (&f_lt($n_val, $g_0)) { if ($n_pt > 0) { return ($g_pt_inf, $g_inf); } else { $ans_pt = 0; $ans_val = &f_gamma($n_val); return ($ans_pt, $ans_val); } } if ($n_pt == 0) { # Dispatch all values that are handled perfectly well by f_gamma # Note that n=167 corresponds to the factorial of 166. if (&f_lt($n_val, $g_gammalim)) { $ans_pt = 0; $ans_val = &f_gamma($n_val); return ($ans_pt, $ans_val); } # %%% Note: The following tests for whether we've computed enough # terms are only accurate when curscale is 14 or less! # Start computing terms of log(gamma) # lower upper $t = &f_sub($n_val, $g_1); # 166.0 1.0e300 $l = $g_05_ln_2pi; # 0.918938533... $l = &f_add($l, &f_mul(&f_add($t, $g_0_5), &f_ln($t))); # 852.0649... 6.907e302 $l = &f_sub($l, $t); # 686.0649... 6.897e302 $n2 = &f_mul($t, $t); $np = $t; $adj = &f_div($g_1, &f_mul($g_12, $np)); # 5.02e-4 0.0 if (&f_eq(&f_add($l, $adj), $l)) { # We're so big that the first 3 terms are enough. So, we're ready to # call exp! (This happens when n is bigger than about 1.0e7) return &pt_exp(0, $l); } # No, keep going for a while... $l = &f_add($l, $adj); $np = &f_mul($np, $n2); $adj = &f_div($g_1, &f_mul($g_360, $np)); # 6.0e-10 0.0 if (&f_eq(&f_sub($l, $adj), $l)) { # For values of n greater than about 1e4, this is as # far as we need to go. return &pt_exp(0, $l); } # Alright, add in 3 more terms (the most it would ever take) # and then call exp. $l = &f_sub($l, $adj); $np = &f_mul($np, $n2); $l = &f_add($l, &f_div($g_1, &f_mul($g_1260, $np))); $np = &f_mul($np, $n2); $l = &f_sub($l, &f_div($g_1, &f_mul($g_1680, $np))); return &pt_exp(0, $l); } elsif ($n_pt == 1) { # N is PT1. # The t = n+1 is irrelevant now, as are other similar additive # elements of the formula, so the formula reduces to # l = n ln(n) - n, gamma(x) = e^l # We convert this to l = n (ln(n) - 1) because otherwise the # error in the subtract would dominate. # There's enough robustness in pt_exp and the other pt_ops # that we can do the calculations directly. # lower upper $ans_pt = $n_pt; # 1.0e300 10^(1.0e300) $ans_val = $n_val; ($ans_pt, $ans_val) = &pt_ln($ans_pt, $ans_val); # 690.775 2.302e300 # Notice: n ln(n) would be 6.90775e302 10^(1.0e300) # and n ln(n) - n would be 6.89775e302 0.0! ($ans_pt, $ans_val) = &pt_sub($ans_pt, $ans_val, 0, $g_1); # 689.775 2.302e300 ($ans_pt, $ans_val) = &pt_mul($n_pt, $n_val, $ans_pt, $ans_val); # 6.89775e302 10^(1.0e300) return &pt_exp($ans_pt, $ans_val); } else { # N is PT2 or higher. # As seen in the upper limit case of the PT1 code above, # for sufficiently large values gamma(n) reduces to pt_exp(n). return &pt_exp($n_pt, $n_val); } } # sqrt is done differently from pow(x, 0.5) because there is a builtin # sqrt routine (which is faster) and because sqrt of a negative number # should have a real component of exactly 0. sub pt_sqrt { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); local($v1); # First handle negatives (including negative infinity), 'cause they're easy if (&f_le($x_val, $g_0)) { return (0, $g_0); } # Handle infinity if (&f_eq($x_val, $g_inf)) { return ($g_pt_inf, $g_inf); } if ($x_pt == 0) { # X is PT0, just call normal square root routine. return (0, &f_sqrt($x_val)); } elsif ($x_pt == 1) { # X is PT1. Divide value by 2 and check for demote. $v1 = &f_div($x_val, $g_2); if (&f_lt($v1, $g_log10_ovr1)) { # demote to PT0 return (0, &f_pow10($v1)); } else { return (1, $v1); } } elsif ($x_pt == 2) { # X is PT2. Subtract log of 2 and check for demote. $v1 = &f_sub($x_val, $g_log10_2); if (&f_lt($v1, $g_log10_ovr1)) { # demote to PT1 return (1, &f_pow10($v1)); } else { return (2, $v1); } } # else # X is PT3 or higher. Square root is uncomputably smaller. return ($x_pt, $x_val); } q| pt_real pt_sin(pt_real x) { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); if ($x_pt == 0) { $ans_pt = 0; $ans_val = f_sin($x_val); return ($ans_pt, $ans_val); } else { return (0, 0); } } pt_real pt_cos(pt_real x) { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); if ($x_pt == 0) { $ans_pt = 0; $ans_val = &f_cos($x_val); return ($ans_pt, $ans_val); } else { return (0, 0); } } pt_real pt_tan(pt_real x) { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); if ($x_pt == 0) { $ans_pt = 0; $ans_val = f_tan($x_val); return ($ans_pt, $ans_val); } else { return ($g_pt_inf, $g_inf); } } # sinh(x) = (e^x - e^-x) / 2 # outside [-40..40] this reduces to - e^-x / 2 or e^x / 2 # and for PT1 and higher arguments it's just - e^-x (negative argument) or e^x pt_real pt_sinh(pt_real x) { local($x_pt, $x_val) = @_; return &pt_div(&pt_sub(&pt_exp(x), &pt_exp(&pt_negate(x))), (0, 2)); } # cosh(x) = (e^x + e^-x) / 2 pt_real pt_cosh(pt_real x) { local($x_pt, $x_val) = @_; return &pt_div(&pt_add(&pt_exp(x), &pt_exp(&pt_negate(x))), (0, 2)); } # tanh(x) = (e^2x - 1) / (e^2x + 1). Outside the PT0 range it's always # 1.0 or -1.0 sub pt_tanh { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); local($e2x); # All PT1's and higher and PT0's outside the range [-20..20] have a tanh # of -1 or 1. if (&f_lt($x_val, -20.0)) { return (0, $g_n1); } elsif (&f_gt($x_val, 20.0)) { return (0, $g_1); } # For the remaining values, we compute the function. The e^2x part is # used twice in the formula for tanh. The answer is always PT0. $e2x = &f_exp(&f_mul($g_2, $x_val)); $ans_pt = 0; $ans_val = ($e2x - 1.0) / ($e2x + 1.0); return ($ans_pt, $ans_val); } # arcsinh(x) = ln(x + sqrt(X^2 + 1)) pt_real pt_arcsinh(pt_real x) { local($x_pt, $x_val) = @_; # Since we'd have to handle all the PT's separately and since arcsinh # isn't used that much, we just don't bother to implement this function # with discrete code. return &pt_ln(&pt_add(x, &pt_sqrt(&pt_add(&pt_mul(x, x), (0, 1))))); } # arccosh(x) = ln(x + sqrt(X^2 - 1)) # undefined if argument is less than 1 pt_real pt_arccosh(pt_real x) { local($x_pt, $x_val) = @_; if ($x_val < 1.0) { # catches all negative PTs as well return (0, 0); } # We don't use discrete code for same reason as with arcsinh. return &pt_ln(&pt_add(x, &pt_sqrt(&pt_sub(&pt_mul(x, x), (0, 1))))); } # arctanh(x) = 1/2 ln((-x-1) / (x-1)) # undefined outside [-1, 1] sub pt_arctanh { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); if ($x_val < -1.0) { # automatically catches negative PT1's etc. return (0, 0); } elsif ($x_val > 1.0) { # automatically catches PT1's etc. return (0, 0); } # The answer is always a PT0, because the biggest ratio we can have # is approximately 1::2^54 (i.e. the precision ratio) $ans_pt = 0; $ans_val = &f_ln((-$x_val - 1.0) / ($x_val - 1.0)) / 2.0; return ($ans_pt, $ans_val); } pt_real pt_arcsin(pt_real x) { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); if ($x_pt > 0) { return (0, 0); } $ans_pt = 0; $ans_val = f_arcsin($x_val); return ($ans_pt, $ans_val); } pt_real pt_arccos(pt_real x) { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); if ($x_pt > 0) { return (0, 0); } $ans_pt = 0; $ans_val = f_arccos($x_val); return ($ans_pt, $ans_val); } |; sub pt_arctan { local($x_pt, $x_val) = @_; local($ans_pt, $ans_val); if (&f_eq($x_val, $g_0)) { return (0, $g_0); } elsif (&f_gt($x_pt, $g_0)) { $ans_pt = 0; if (&f_gt($x_val, $g_0)) { $ans_val = &f_div($g_pi, $g_2); } else { $ans_val = &f_div($g_pi, $g_n2); } return ($ans_pt, $ans_val); } $ans_pt = 0; $ans_val = &f_arctan($x_val); return ($ans_pt, $ans_val); } sub pt_roundup { local($x_pt, $x_val) = @_; local($a_val, $ptsave); if ($doround) { # Temporarily coerce to PT0 $ptsave = $x_pt; $x_pt = 0; # Calculate adjust amount. This is a cheat -- it correctly accounts for # sign but does not set proper rounding amount for round-to-nearest. # It should actually be sign(x_val) * 0.5 * 10^(x_exponent - curscale) $a_val = &f_mul($x_val, $roundprec); ($x_pt, $x_val) = &pt_add($x_pt, $x_val, 0, $a_val); # Add the saved PT back in (because the roundup might have increased # the PT) $x_pt += $ptsave; } return($x_pt, $x_val); } #------------------------------------------------- hypercalc.c # comment for the Pilot version: # f8_decomp1 takes a double and returns its sign, mantissa and exponent # as separate quantities (integer, double, integer respectively). It is # used for output formatting. # test patterns q@ -20 ^ 238.4999999999 -- currently gives 6.2e300 as a PT1. Don't know what it *should* give, there's a lot of error because the complex vector is close to the imaginary axis. 0.0001 * (10^151 * 10^151) -- Should give 1e298, as a PT0 not PT1 -- for testing formatting of small numbers in stack display E ^ 30 E ^ 300 E ^ 3e3 E ^ 3e4 E ^ 3e5 E^(E^30) E^(E^300) E^(E^3e3) E^(E^3e4) E^(E^3e5) E^(E^(E^30)) E^(E^(E^300)) E^(E^(E^3e3)) E^(E^(E^3e4)) E^(E^(E^3e5)) E^(E^(E^(E^30))) (etc...) A simple "iterative" way to keep generating larger values (for display testing) is to repeat "%^ln(%)" over and over again. -1*10^1e200 -- should give -1e2000 4 root % -- should give 7.07e500 1e10 ^ 5e299 -- should give 10^(5e300) as a PT2 20 ^ 238.4999999999 -- should give 1.975e310 -20 ^ 238.49999999999 -- should give 6.2e299, as a PT0 (NOT a PT1) -- for testing infinities 1 ent 0 / -- should give Infinity 2 ent 1 logxY -- should give Infinity 0 ent 1 logxY -- could give -Infinity or Infinity, I don't care which 1 ent 0 / 1 x -- Infinity times 1 should be infinity 20 ent E299 x eX -- should give approx. 10^(10^300) lnx -- should give 2e300 e200 eX +- eX -- should give 0 scale = 120 27^(10^1e100) -- should give 2 PT 1.000...000155e100, with about 103 zeros scale=50 5^2*4 - gives 9.9999999, etc. Should round up and increment exponent #Non-Intuitive Results When Working With Huge Numbers# If you spend a while exploring the ranges of huge numbers HyperCalc can handle, you will probably start noticing some paradoxical results and might even start to think the calculator is giving wrong answers. For example, try calculating 27 to the power of *googolplex* (a googolplex is 10 to the power of *googol* and a googol is 10^{100}). Key in: 27^(10^(10^100)) and it prints 10^(10^(1.00 x 10^100)) so the calculator thinks that 27^(10^(10^100)) = 10^(10^(10^100)) This is clearly wrong -- and it doesn't even seem to be a good approximation. What's going on? Let's try calculating the correct answer ourselves. We need to express the answer as 10 to the power of 10 to the power of something, because that's the standard format the calculator is using, and we're going to see how much of an error it made. So, we want to compute 27^{10^{10^{100}}} as a tower of powers of 10. The first step is express the power of 27 as a power of 10 with a product in the exponent, using the formula x^{y} = 10^{(log(x) ^{.} y)} : 27^{10^{10^{100}}} = 10^{(logv{10}27 ^{.} 10^{10^{100}})} logv{10}27 is about 1.43, so we have 27^{10^{10^{100}}} = 10^{1.43 ^{.} 10^{10^{100}})} Now we have a base of 10 but the exponent still needs work. The next step is to express the product as a sum in the next-higher exponent; this time the formula we use is x ^{.} y = 10^{(log(x) + log(y))} : 10^{1.43 ^{.} 10^{10^{100}})} = 10^{10^{(logv{10}1.43 + 10^{100})}} logv{10}1.43 is about 0.155, and if we add this to 10^{100} we get 10^{10^{(0.155 + 10^{100})}} = 10^{10^{1000...000.155}} = 10^{10^{(1.000...000155 ^{.} 10^{100})}} where there are 94 more 0's in place of each of the "...". So our final answer is: 27^{10^{10^{100}}} = 10^{10^{(1.000...000155 ^{.} 10^{100})}} Now that we've expressed the value of 27^*googolplex* precisely enough to see the calculator's error -- look how small the error is! The calculator would need to have at least 104 digits of precision to be able to handle the value "1.000...000155" accurately -- but it only has 16 digits of accuracy. Those 16 digits are taken up by the 1 and the first 15 0's -- so when the calculator gets to the step where we're adding 0.155 to 1.0^{.}10^{100}, it just rounds off the answer to 1.0^{.}10^{100} -- and produces the answer we saw when we performed the calculation: 10^(10^(1.00 x 10^{100})) Even if it did have the precision, it wouldn't have room to print the whole 104 digits on the screen, so the answer you *see* would look the same. And no matter how many digits of accuracy we try to give the calculator, there's always another even bigger number it wouldn't be able to handle. For example, the calculator would need slightly over a *million* digits of accuracy to distinguish 27^{10^{10^{1000000}}} from 10^{10^{10^{1000000}}} and if we just add one more #10# to that tower of exponents, all hope of avoiding roundoff is lost. # 2^(2^(2+27)) = 2.048696 * 10^161614248 */ more tests: To check each of the special cases in output formatting: 143^143 %! %! %! 143! %! %! %! to check promotion of inlines: 10^23e456 @; #----------------------------------------------- "p_" routines # # each of these routines does the same thing -- converts one # or two arguments in string PT format to separate PT and VAL values, # then calls the corresponding "pt_" routine, then converts # the answer back to string PT format. # sub define_ir { local($pt, $val) = @_; $ir_n++; $ir_pt[$ir_n] = $pt; $ir_val[$ir_n] = $val; &dbg1("I$ir_n := $pt PT $val", 2); return "I$ir_n"; } sub define_R_hist { local($pt, $val) = @_; $h_n++; $h_pt[$h_n] = $pt; $h_val[$h_n] = $val; &dbg1("R$h_n := $pt PT $val", 2); return "R$h_n"; } sub p_add { local($x, $y) = @_; local($x_pt, $x_val, $y_pt, $y_val); ($x_pt, $x_val) = &split($x); ($y_pt, $y_val) = &split($y); ($x_pt, $x_val) = &pt_add($x_pt, $x_val, $y_pt, $y_val); return &define_ir($x_pt, $x_val); } sub p_sub { local($x, $y) = @_; local($x_pt, $x_val, $y_pt, $y_val); # print "p_sub $x $y\n"; ($x_pt, $x_val) = &split($x); ($y_pt, $y_val) = &split($y); ($x_pt, $x_val) = &pt_sub($x_pt, $x_val, $y_pt, $y_val); return &define_ir($x_pt, $x_val); } sub p_mul { local($x, $y) = @_; local($x_pt, $x_val, $y_pt, $y_val); ($x_pt, $x_val) = &split($x); ($y_pt, $y_val) = &split($y); ($x_pt, $x_val) = &pt_mul($x_pt, $x_val, $y_pt, $y_val); return &define_ir($x_pt, $x_val); } sub p_div { local($x, $y) = @_; local($x_pt, $x_val, $y_pt, $y_val); ($x_pt, $x_val) = &split($x); ($y_pt, $y_val) = &split($y); ($x_pt, $x_val) = &pt_div($x_pt, $x_val, $y_pt, $y_val); return &define_ir($x_pt, $x_val); } sub p_negate { local($x) = @_; local($x_pt, $x_val); ($x_pt, $x_val) = &split($x); ($x_pt, $x_val) = &pt_negate($x_pt, $x_val); return &define_ir($x_pt, $x_val); } sub p_exp { local($x) = @_; local($x_pt, $x_val); ($x_pt, $x_val) = &split($x); ($x_pt, $x_val) = &pt_exp($x_pt, $x_val); return &define_ir($x_pt, $x_val); } sub p_ln { local($x) = @_; local($x_pt, $x_val); ($x_pt, $x_val) = &split($x); ($x_pt, $x_val) = &pt_ln($x_pt, $x_val); return &define_ir($x_pt, $x_val); } sub p_arctan { local($x) = @_; local($x_pt, $x_val); ($x_pt, $x_val) = &split($x); ($x_pt, $x_val) = &pt_arctan($x_pt, $x_val); return &define_ir($x_pt, $x_val); } sub p_log10 { local($x) = @_; local($x_pt, $x_val); ($x_pt, $x_val) = &split($x); ($x_pt, $x_val) = &pt_log10($x_pt, $x_val); return &define_ir($x_pt, $x_val); } sub p_pow { local($x, $y) = @_; local($x_pt, $x_val, $y_pt, $y_val); ($x_pt, $x_val) = &split($x); ($y_pt, $y_val) = &split($y); ($x_pt, $x_val) = &pt_pow($x_pt, $x_val, $y_pt, $y_val); return &define_ir($x_pt, $x_val); } sub p_fact { local($x) = @_; local($x_pt, $x_val); ($x_pt, $x_val) = &split($x); ($x_pt, $x_val) = &pt_gamma(&pt_add($x_pt, $x_val, 0, $g_1)); return &define_ir($x_pt, $x_val); } sub p_sqrt { local($n) = @_; local($n_pt, $n_val); ($n_pt, $n_val) = &split($n); ($n_pt, $n_val) = &pt_sqrt($n_pt, $n_val); return &define_ir($n_pt, $n_val); } sub p_root { local($r, $n) = @_; local($r_pt, $r_val, $n_pt, $n_val); ($r_pt, $r_val) = &split($r); ($n_pt, $n_val) = &split($n); ($r_pt, $r_val) = &pt_root($n_pt, $n_val, $r_pt, $r_val); return &define_ir($r_pt, $r_val); } sub p_logn { local($b, $n) = @_; local($b_pt, $b_val, $n_pt, $n_val); ($b_pt, $b_val) = &split($b); ($n_pt, $n_val) = &split($n); ($n_pt, $n_val) = &pt_log_n($n_pt, $n_val, $b_pt, $b_val); return &define_ir($n_pt, $n_val); } < "; $s = ($min * 60) + $sec; $e = $pre . $s . $post; # print "$e\n"; } # now it's safe to eliminate spaces $e =~ s| ||g; # put numbers in standard form $e =~ s|([-\(\|\/*xe+,^])-|$1N|g; $e =~ s|([0-9.])e([N+0-9.])|$1E$2|g; $e =~ s|([0-9])p([N+0-9.])|$1P$2|g; &dbg1("ev1c $e", 1); # do the actual parsing and evaluation $e = &eval_2($e); &dbg1("ev1d $e", 1); $e =~ s'\|''g; &dbg1("ev1e $e", 1); if ($e =~ m|^I([0-9]+)$|) { # result is a single IR $e = $1; $a2_pt = $ir_pt[$e]; $a2_val = $ir_val[$e]; dbg1("ev1f res I$e = $a2_pt" . "P$a2_val", 2); } else { # result is a literal ($a2_pt, $a2_val) = &split($e); dbg1("ev1g lit $a2_pt" . "P$a2_val", 2); } dbg1("ev1h $a2_pt" . "P$a2_val", 2); return ($a2_pt, $a2_val); } #--------------------------------------------------- main # # initialization, main command loop, help, etc. # Compute pi using the Baily-Plouffe hexadecimal series, which # converges quite quickly sub init_pi_1 { local($f1, $i, $d); &dbg1("init_pi_1", 2); $g_pi = $g_0; $f1 = $g_1; $i = $g_0; $d = $g_1; while (&f_gt($f1, $prec)) { $i_f = $i; $i_f = &f_mul($g_8, $i_f); $f1 = &f_div($g_4, &f_add($i_f, $g_1)); $f1 = &f_sub($f1, &f_div($g_2, &f_add($i_f, $g_4))); $f1 = &f_sub($f1, &f_div($g_1, &f_add($i_f, $g_5))); $f1 = &f_sub($f1, &f_div($g_1, &f_add($i_f, $g_6))); $f1 = &f_div($f1, $d); $d = &f_mul($d, $g_16); &dbg1("hi_init g_pi term", 1); $g_pi = &f_add($g_pi, $f1); $i = &f_add($i, $g_1); } &dbg1("init_pi_1 done: " . &prnt2($g_pi), 2); } # Compute pi using a formula based on the Gauss arithmetic-geometric mean # method. Number of digits doubles with each loop. sub init_pi_2 { local($x, $y, $a, $b, $c, $ct, $i); &dbg1("init_pi_2", 2); $a = $x = $g_1; $b = &f_div($g_1, &f_sqrt($g_2)); $c = &f_div($g_1, $g_4); for($i=1; $i<$curscale; $i*=2) { $y = $a; $a = &f_div(&f_add($a, $b), $g_2); $b = &f_sqrt(&f_mul($b, $y)); $ct = &f_sub($a, $y); $c = &f_sub($c, &f_mul($x, &f_mul($ct, $ct))); $x = &f_mul($x, $g_2); } $ct = &f_add($a, $b); $g_pi = &f_div(&f_mul($ct, $ct), &f_mul($g_4, $c)); &dbg1("init_pi_2 done: " . &prnt2($g_pi), 2); } # This part of initialization is the same for both types of primitives. sub hi_init { local($i, $f); local($prec); local($f1, $i_f, $d); local($savedbg); $g_ln_10 = &f_ln($g_10); $g_log10_e = &f_div($g_1, $g_ln_10); &dbg1("hi_init prec", 1); $prec = &f_div($g_1, &f_pow10(&f_add($curscale, $g_2))); $roundprec = &f_div($g_0_5, &f_pow10($curscale)); $gammaprec = &f_pow10(&f_div($curscale, $g_15)); # Constants used for handling roundoff, overflow, boundary # conditions, etc. $g_ovr1 = "+1.0E300"; $g_log10_ovr1 = "300.0"; # ln(10^300) = ln(10) * 300 $g_ln_ovr1 = &f_mul($g_ln_10, $g_log10_ovr1); $g_gammalim = "24.0"; # Check if they're setting the scale higher than the precision of our Gamma # function. %%% note: a better Gamma function is implemented in PARI, # functions mpgamma() and mplngamma() in src/basemath/trans2.c if (&f_gt($gammaprec, $g_gammalim)) { $gammaprec = $g_gammalim; $gammadig = &f_int(&f_log10(&f_pow($gammaprec, $g_15))); $msg1 = "Note: Factorial will only give " . &prnt2($gammadig) . " digits of accuracy.\n"; print $msg1; } # high-level constants # compute constant e $g_e = $g_1; $f = 1; for($i=1; ($i < 172) && (&f_gt($f, $prec)); $i++) { $f = &f_div($f, $i); $g_e = &f_add($g_e, $f); } # Compute pi using the Baily-Plouffe hexadecimal series, which # converges quite quickly &init_pi_2(); $g_sin_limit = &f_mul($g_2p32m1, $g_pi); $g_phi = &f_div(&f_add(&f_sqrt($g_5), $g_1), $g_2); &dbg1("hi_init g_05_ln_2pi", 1); $g_05_ln_2pi = &f_mul($g_0_5, &f_ln(&f_mul($g_2, $g_pi))); $g_log10_log10_e = &f_log10($g_log10_e); $g_log10_2 = &f_log10($g_2); # native floating point infinities. These count on f_mul generating # the infinite value. $g_inf = 1.0e200; $g_ninf = &f_mul($g_inf, -1.0e200); $g_inf = &f_mul($g_inf, $g_inf); # PT Infinities have a native (double) infinity in the PT field $g_pt_inf = 1.0e200; $g_pt_inf = &f_mul($g_pt_inf, $g_pt_inf); # constants used for output formatting $g_ee_b4 = &f_pow10($curscale - 1); } # Print a value that is in floating-point or scientific notation, converting # if necessary to floating-point and adding or truncating mantissa digits # to match the current scale. sub prnt2 { local($x, $paren) = @_; local($m1, $e, $s, $m2, $m3, $s2, $output, $e1, $e2); print "prnt2($x, $paren)..." if ($format_debug); $output = ""; # Special cases if ($x eq 'inf') { print "prnt2: Infinity\n" if ($format_debug); return("Infinity"); } # Change "1.23e+21" to "1.23e21" $x =~ s|\+||g; # If it's in fixed format, but has more digits than the current precision # allows, convert to floating-point if (!($x =~ m|[eE]|)) { # extract sign if ($x =~ m/\-/) { $s = "-"; $x =~ s/\-//; } else { $s = ""; } if (&f_gt($x, $g_ee_b4)) { if ("$x.0" =~ m/^([0-9])([0-9]*)\.([0-9]+)/) { $m1 = $1; $m2 = $2; $m3 = $3; $e = length($m2); if ($e > 0) { $x = $m1 . "." . $m2 . $m3 . "e" . $e; } } else { print "prnt2 error converting fixed to EE: |$x|\n"; } } $x = $s . $x; } if ($x =~ m|[eE]|) { # Translate scientific notation to something prettier if ($x =~ m|^(-?)([0-9]+)\.([0-9]*)[eE]([-0-9]+)$|) { $s = $1; $m1 = $2; $m2 = $3; $e = $4; } elsif ($x =~ m|^(-?)([0-9]+)[eE]([-0-9]+)$|) { $s = $1; $m1 = $2; $m2 = "0"; $e = $3; } elsif ($x =~ m|[eE]|) { print "error-0-01 Unhandled format: $x\n"; exit 0; } $s2 = $curscale - 1; if ($m2 eq "") { $m2 = "0"; } # We should never have "12.34e56" if (length($m1) > 1) { print("error-0-02 Badformat mantissa: $s$m1.$m2" . "e$e\n"); exit 0; } # truncate mantissa if (length($m2) > $s2) { $m2 =~ m|^([0-9]{$s2})|; $m2 = $1; } $output .= $p1_lstr if ($paren); # Decide whether to print as a normal number if (($e >= 0) && ($e <= $s2)) { # Get appropriate number of digits $takedig = $e; while(length($m2) < $takedig) { $m2 .= "0"; } $m2 =~ m|^([0-9]{$takedig})([0-9]*)$|; $m1 = $m1 . $1; $m2 = $2; $x = "$s$m1.$m2"; $x =~ s|0+$||; $x =~ s|\.$||; $output .= $x; } else { $x = "$s$m1.$m2"; $x =~ s|0+$||; $x =~ s|\.$||; $output .= "$x$multstr$k10$powstr$e_lstr$e$e_rstr"; } $output .= $p1_rstr if ($paren); } else { # Just a plain number # Extract sign if ($x =~ m/\-/) { $s = "-"; $x =~ s/\-//; } else { $s = ""; } # Truncate digits to scale if ("$x.0" =~ m/^([0-9]+)\.([0-9]+)/) { $m1 = $1; $m2 = $2; $e1 = length($m1); $e2 = length($m2); $e = $curscale - $e1; if ($e2 > $e) { if ($m2 =~ m/^(.{$e})/) { $m2 = $1; } else { print "prnt2 error-2 truncating fixed to scale: |$x|\n"; } $x = $m1 . "." . $m2; } } else { print "prnt2 error-1 truncating fixed to scale: |$x|\n"; } # Discard any trailing ".000" if ($x =~ m|\.|) { $x =~ s|0+$||; $x =~ s|\.$||; } $output = $s . $x; } print "prnt2: $output\n" if ($format_debug); return $output; } # Constants and functions used by the label formatting routines $k_1e5 = 100000.0; $k_1e6 = $k_million = 1000000.0; $k_1e7 = 10000000.0; $k_1e10 = 1.0e10; # Truncates a string to a given number of characters. Before doing so, # pads the end with some character if desired. sub trunc2 { local($x, $n, $extend) = @_; if ($extend ne "") { $x = $x . ($extend x $n); } if ($x =~ m/^(.{$n})/) { $x = $1; } return $x; } # This routine generates the old-style label for numbers.rhtf # # Examples: # # label Internal form Value # l5390en44 0-pt-5.39e-44 5.390 x 10^-44 # l0_2078 0-pt-0.2078 0.2078 # l1_4142 0-pt-1.4142 1.4142 # l15_154 0-pt-15.154 15.154 # l127 0-pt-127 127 # l720720 0-pt-720720 720720 # l1419_c 0-pt-1419869 1419869 # l_1234e10 0-pt-1.234e10 12345654321 or 1.24 x 10^10 # l_1010e38 0-pt-1.010e38 1.010 x 10^38 # l_p1_340_792 1-pt-340.792 6.2 x 10^340 # l_p1_2098e6 1-pt-2.0989e6 4 x 10^2098959 # l_p1_1000e166 1-pt-1.0e166 10 ^ 10 ^ 166 # l_p2_4342e12 2-pt-4.342e12 10 ^ 10 ^ 4342944819032 # l_p3_36305 4-pt-36305.315 10 ^ 10 ^ 10 ^ 2.069 x 10^36305 sub format4 { local($apt, $av) = @_; local($rv, $as, $am, $ae); ($as, $am, $ae) = &fbc_split($av); if ($apt > 0) { if ($av >= $k_million) { # l_p1_1000e166 $rv = "l_p" . $apt . "_" . &trunc2($am, 5, "0") . "e" . $ae; $rv =~ s/\.//; } else { # l_p1_340_792 $rv = "l_p" . $apt . "_" . &trunc2($av, 7, ""); $rv =~ s/\./_/; } } elsif ($av >= $k_1e10) { # l_1701e38 $rv = $am; $rv =~ s/\.//; $rv = &trunc2($rv, 4, "0"); $rv = "l_" . $rv . "e" . $ae; } elsif ($av >= $k_million) { # l1419_c $am =~ s/\.//; $rv = &trunc2($am, 4, ""); $ae =~ tr/6789/cdef/; $rv = "l_" . $rv . "_" . $ae; } elsif ($av >= 0.1) { # l0_2078, l1_4142, l15_154, l127, l720720 $rv = "l" . &trunc2($av, 6, ""); if ($rv =~ m/\./) { $rv =~ s/\./_/; $rv = &trunc2($rv, 7, "0"); } } else { $rv = $am; $rv =~ s/\.//; $rv = &trunc2($rv, 4, "0"); $ae =~ s/-/n/; $rv = "l" . $rv . "e" . $ae; } return $rv; } # Formats a real number with a given number of significant figures, # with no scientific notation; then substitutes "." (if any) with "_". sub canon_real { local($v, $p) = @_; if ($v =~ m/^\./) { $v = "0" . $v; } # Ensure digit before decimal point $p++; if ($v =~ m/^([.0-9]{$p})/) { $v = $1; } # If too long, truncate $v =~ s/\./_/; $v =~ s/_0+$//; # Remove needless trailing 0's $t = $v . "_"; if ($t =~ m/^[0-9][0-9]_/) { $v = "a" . $v; } if ($t =~ m/^[0-9]{3}_/) { $v = "b" . $v; } if ($t =~ m/^[0-9]{4}_/) { $v = "c" . $v; } if ($t =~ m/^[0-9]{5}_/) { $v = "d" . $v; } return $v; } # Formats a scientific-notation number with a given number of # exponent and mantissa digits for use within format5 labels. # Arguments: # # $e exponent # $pe precision for exponent # $m mantissa # $pm precision for mantissa # # For example, if you want to format 6.02x10^23 with a precision of 3 in # both parts, you would call &canon_sci("23", 3, "6.020", 3), which # becomes 023_602 sub canon_sci { local($e, $pe, $m, $pm) = @_; local($ra); print "canon_sci($e, $pe, $m, $pm)...\n" if ($format_debug); # Check for special cases if ($e eq 'inf') { print "$e\n" if ($format_debug); return($e); } if (!($e =~ m/[0-9]{$pe}/)) { $e = ("0" x $pe) . $e; if ($e =~ m/([0-9]{$pe})$/) { $e = $1; } } # truncate precision of mantissa to one leading digit plus $pm trailing # digits. (The extra digit will be removed during rounding) $m = $m . ("0" x $pm); if ($m =~ m/^([1-9][.][0-9]{$pm})/) { $m = $1; } else { print STDERR "canon_sci (e1): Cannot parse mantissa '$m'\n"; exit(-1); } # Perform rounding $ra = "0." . ("0" x ($pm-1)) . "5"; # If $pm is 3, this will be "0.005" $ra += 0; if ($m + $ra < 10.0) { $m = $m + $ra; } else { # We leave "9.99998" alone so it is distinguishable from "10.00002" } # Truncate again to remove the rounding digit if ($m =~ m/^([1-9][0-9.]{$pm})/) { $m = $1; } elsif (($pm == 3) && ($m =~ m/^([1-9]\.[0-9])/)) { $m = $1; } else { print STDERR ( "canon_sci (e2): Cannot re-truncate mantissa '$m'\n" . " by pattern '^([1-9][0-9.]{$pm})'; consider adding another\n" . " special case.\n"); exit(-1); } # remove mantissa embedded '.' $m =~ s/\.//; # remove mantissa trailing 0's $m =~ s/0+$//; print ($e. "_$m\n") if ($format_debug); return ($e . "_" . $m); } # format5 implements the new-style labels that are completely alphabetizable: # # l_234em34 2.34 x 10^-34 # l0_1 0.1 # l9_34 9.34 # la10_2 10.2 # lb256 256 # lc1729 1729 # ld19683 19683 # le005_1 100000 = 1x10^5 # le023_602 6.02 x 10^23 # le299_9 9 x 10^299 # lp1_b300_0 1 PT 300 = 1x10^300 # lp1_c2345_6 1 PT 2345.6 {~=} 4x10^2345 # lp1_e009_345 1 PT 3.45e9 = 10^(3.45*10^9) # lp1_e299_999 1 PT 9.99e299 = 10^(9.99x10^299) # lp2_b300_0 2 PT 300 = 10^(10^300) # lpa10_xxx 10 PT xxx # lpa99_xxx 99 PT xxx # lpb100_xxx 100 PT xxx sub format5 { local($apt, $av) = @_; local($rv, $as, $am, $ae); print "format5($apt, $av)...\n" if ($format_debug); if ($av eq 'inf') { print "format5: = lz_inf\n" if ($format_debug); return("lz_inf"); } ($as, $am, $ae) = &fbc_split($av); if ($apt > 9) { if ($av >= $k_1e5) { # l_pa12_e009_345 $rv = "lp" . &canon_real($apt, 6) . "_e" . &canon_sci($ae, 3, $am, 3); } else { # l_pa12_c2345_6 $rv = "lp" . &canon_real($apt, 6) . "_" . &canon_real($av, 6); } } elsif ($apt > 0) { # for now I'll skip the test for PT > 9 if ($av >= $k_1e5) { # l_p1_e009_345 $rv = "lp" . $apt . "_e" . &canon_sci($ae, 3, $am, 3); } else { # l_p1_c2345_6 $rv = "lp" . $apt . "_" . &canon_real($av, 6); } } elsif ($av >= $k_1e7) { # le032_602 $rv = "le" . &canon_sci($ae, 3, $am, 3); } elsif ($av >= $k_1e6) { # Use 4 digits of mantissa for 7-digit numbers like 1048576, giving # le006_1048 $rv = "le" . &canon_sci($ae, 3, $am, 4); } elsif ($av >= $k_1e5) { # Use 5 digits of mantissa for 6-digit numbers like 196560, giving # le005_19656 $rv = "le" . &canon_sci($ae, 3, $am, 5); } elsif ($av >= 0.1) { # l0_1, l9_34, la10_2, lb256, ld19683 $rv = "l" . &canon_real($av, 6); } else { # Smaller than 0.1 -- not handling alphabetization yet $am =~ s/\.//; $ae =~ s/-//; $rv = "l_" . &trunc2($am, 3, "") . "em" . $ae; } print "format5: = $rv\n" if ($format_debug); return $rv; } # format 6 just prints the internal representation of a number (PT and Val) sub format6 { local($apt, $av) = @_; local($rv); $rv = "|PT" . $apt . "V" . $av . "|"; return $rv; } # format 7 outputs small numbers in HH:MM:SS format sub format7 { local($apt, $av) = @_; local($rv, $h, $m, $s, $f); if (($apt > 0) || ($av >= 60*3600) || ($av < 0)) { return(&format6($apt, $av)); } $f = $av - int($av); $av = int($av); $m = int($av / 60); $s = $av - (60 * $m); $h = int($m / 60); $m = $m - (60 * $h); $m = sprintf("%02d", $m); $s = sprintf("%02d", $s); $rv = "$h:$m:$s"; if ($f > 0) { $f =~ s/^ *0*\.//; $rv .= ".$f"; } return $rv; } # print a PT, given its pt and val as separate arguments. sub prnt1 { local($apt, $av, $newline) = @_; local($avi, $avf, $avl, $avl2, $output); local($f5_lbl); print "prnt1($apt, $av, $newline)...\n" if ($format_debug); if ($fmt_callproc eq "format4") { print(&format4($apt, $av) . ":: "); } elsif ($fmt_callproc eq "format5") { $f5_lbl = &format5($apt, $av); print($f5_lbl . "::\n#"); } elsif ($fmt_callproc eq "format6") { print(&format6($apt, $av) . " "); } elsif ($fmt_callproc eq "format7") { print(&format7($apt, $av) . " "); } if ($doround) { # round up ($apt, $av) = &pt_roundup($apt, $av); } $output = ""; if ($av eq 'inf') { $output = "Infinity"; } elsif ($apt == 0) { # Format a PT0. # print "prnt1 case la-le\n"; $output = &prnt2($av, 0); } else { # We need to print either "$apt PT $av" or "($apt-1) PT (10^$av)" # depending on whether or not $av is less than 10^scale if (&f_lt($av, $g_ee_b4)) { # $av is less than 10^scale -- we have a value like 1P2345.67 or # 1P2.34567e3. We have to subtract 1 from pt and compute mantissa # of 10 ^ val # get exponent $avi = &f_int($av); # figure out how many digits of mantissa we can print $avl = &f_int(&f_log10($avi)); $avl = $curscale - $avl; $avl =~ s| ||g; $avl2 = $avl + 2; # 2 extra for sign and decimal point $avl2 =~ s| ||g; # compute mantissa $avf = &f_sub($av, $avi); $avf = &f_pow10($avf); # discard insignificant digits of mantissa. We know $avf is # of the form D.DDDDDDDDD where D is a digit. if (length($avf) > $avl2) { $avf =~ m|^(.{$avl2})|; $avf = $1; } # Discard any trailing ".000" if ($avf =~ m|\.|) { $avf =~ s|0+$||; $avf =~ s|\.$||; } # Discard leading '+' $avf =~ s|^\+||; # Handle first three PT classes separately. if ($apt == 1) { # print "prnt1 case p1_a\n"; $output .= &prnt2($avf, 0) . "$multstr$k10$powstr$e_lstr" . &prnt2($avi, 0) . $e_rstr; } elsif ($apt == 2) { # print "prnt1 case p2_a\n"; $output .= "$k10$powstr$e_lstr$p1_lstr$avf$multstr$k10$powstr$e_lstr" . &prnt2($avi, 0) . "$e_rstr$p1_rstr$e_rstr"; } elsif ($apt == 3) { # print "prnt1 case p3_a\n"; $output .= "$k10$powstr$e_lstr$p2_lstr$k10$powstr$e_lstr$p1_lstr" . "$avf$multstr$k10$powstr$e_lstr" . &prnt2($avi, 0) . "$e_rstr$p1_rstr$e_rstr$p2_rstr$e_rstr"; } else { # print "prnt1 case p4_a\n"; $apt -= 1; $output .= "$apt$pt_str$p1_lstr$avf$multstr$k10$powstr$e_lstr" . &prnt2($avi, 0) . "$e_rstr$p1_rstr"; } } else { if ($apt == 1) { # print "prnt1 case p1_e\n"; $output .= "$k10$powstr$e_lstr" . &prnt2($av, 1) . $e_rstr; } elsif ($apt == 2) { # print "prnt1 case p2_e\n"; $output .= "$k10$powstr$e_lstr$p2_lstr$k10$powstr$e_lstr" . &prnt2($av, 1) . "$e_rstr$p2_rstr$e_rstr"; } else { # print "prnt1 case p3_e\n"; $output .= "$apt$pt_str" . &prnt2($av, 1); } } } $output =~ s| +| |g; while(length($output) > $lleft) { $output =~ m/^(.{$lleft})(.*)$/; print $1; print "\\\n"; $lleft = 78; $output = $2; } print $output; if ($fmt_callproc eq "format5") { print("#\n See also [$output|#$f5_lbl]."); } if ($newline) { print "\n"; $lleft = 78; } } # Sets the strings used to format numbers sub format1 { local($fmt) = @_; $k10 = "10"; # Handle the special types if ($fmt == 4) { $fmt_callproc = "format4"; $fmt = 2; } elsif ($fmt == 5) { $fmt_callproc = "format5"; $fmt = 2; } elsif ($fmt == 6) { $fmt_callproc = "format6"; $fmt = 1; } elsif ($fmt == 7) { $fmt_callproc = "format7"; $fmt = 1; } else { $fmt_callproc = ""; } # 20100113: I decided to make p1_lstr, etc. null because () and [] # seem superfluous when there is or {} if ($fmt == 3) { $multstr = "×"; $powstr = ""; $e_lstr = ""; $e_rstr = ""; $p1_lstr = ""; # "("; $p1_rstr = ""; # ")"; $p2_lstr = ""; # "["; $p2_rstr = ""; # "]"; $pt_str = " PT "; } elsif ($fmt == 2) { $multstr = "{x}"; $powstr = "^"; $e_lstr = "{"; $e_rstr = "}"; $p1_lstr = ""; # "("; $p1_rstr = ""; # ")"; $p2_lstr = ""; # "["; $p2_rstr = ""; # "]"; $pt_str = "\@pt\@"; } else { $multstr = " x "; $powstr = " ^ "; $e_lstr = ""; $e_rstr = ""; $p1_lstr = " ( "; $p1_rstr = " ) "; $p2_lstr = " [ "; $p2_rstr = " ] "; $pt_str = " PT "; } } # sub expand.ihist # { # local($cmd) = @_; # local($rv); # # $rv = $cmd; # while(" $rv " =~ m|^(.+)[cC]([0-9]+)([^0-9].*)$|) { # $rv = $1 . "(" . $input_history[$2] . ")" . $3; # } # return $rv; # } sub expand_C_hist { local($e) = @_; local(%h); local($pre, $post, $tn, $te, $tv, $rv); $rv = $e; while(" $rv " =~ m|^(.*[^a-zA-Z_0-9])[cC]([0-9]+)([^0-9].*)$|) { $pre = $1; $tn = $2; $post = $3; if ($h{$tn} ne "") { $tv = $h{$tn}; } else { $te = $input_history[$tn]; $tv = &expand_C_hist($te); # Here we would want to actually evaluate $tv and get a string # numeric result, but that's too hard to do right now. $h{$tn} = $tv; # print "tn $tn te $te tv $tv\n"; } print "C$tn: $tv\n"; $lleft = 78; $rv = $pre . "(" . $tv . ")" . $post; } return $rv; } # expand_ih2 interpolates commands like "c3" into the input string, # and changes 'x' to '*' when appropriate. sub expand_ih2 { local($e) = @_; local($rv, $setvar); $rv = &expand_C_hist($e); $setvar = ""; # Detect parenthesized variable assignment, and convert # this is disabled for now because I don't know why it was present. # It changes input "(x=2*3)" to "x=(2+3)" if (0) { if ($rv =~ m|^\(($p_var) *=(.+)\)$|) { $setvar = $1; $e = $2; $rv = "$setvar=($e)"; &dbg1("setvar1 |$setvar|$e|\n", 1); } } # Pre-translate 'x' for multiplication when next to parentheses. These # are the cases where it is clearly not a variable 'x'. $rv =~ s|\)x |\)* |g; $rv =~ s| x\(| *\(|g; $rv =~ s|\)x\(|\)*\(|g; # Allow them to use 'x' to multiply if padded with spaces and next to # a variable. $rv =~ s|($p_var) x |\1 * |g; $rv =~ s| x ($p_var)| * \2|g; $rv =~ s| +||g; return $rv; } sub define_Chist { local ($cmd, $showeach) = @_; if ($cmd =~ m/^\((.*)\)$/) { $cmd = $1; } if ($showeach) { print "C$ihnum = $cmd\n"; $lleft = 78; } $input_history[$ihnum] = $cmd; $ihnum++; } # docmd3 handles 'C[0-9]' substitution and '='; the rest # is done by eval_1 # # If $BASIC_running $prflag controls printing: 0 for none, 1 for just the # number, 2 for surrounding formatting sub docmd3 { local($cmd, $showeach, $prflag) = @_; local($pr, $sverr); # print "docmd3 cmd '$cmd' showeach $showeach prflag $prflag\n"; if ($cmd eq "") { return; } if ($BASIC_running) { $pr = $prflag; } else { $pr = 2; } &dbg1("dc3.a: $cmd", 1); # First, save input string into input history # Expand first, to prevent infinite recursion $cmd = &expand_ih2($cmd); &dbg1("dc3.b: $cmd", 1); if ($pr) { &define_Chist($cmd, $showeach && ($pr>1)); } &dbg1("cmd $cmd\n", 1); # Parse out variable assignment syntax, if any $setvar = ""; if ($cmd =~ m|^($p_var) *=(.+)$|) { $setvar = $1; $cmd = $2; &dbg1("setvar2 $setvar $cmd\n", 1); } elsif ($cmd =~ m|^\(($p_var) *=(.+)\)$|) { $setvar = $1; $cmd = $2; &dbg1("setvar3 $setvar $cmd\n", 1); } $sverr = ""; if (($setvar eq 'pi') || ($setvar eq 'e') || ($setvar eq 'phi')) { $sverr = "Attempt to set builtin constant"; } elsif ($setvar =~ m/^[cr][0-9]+$/) { # %%% NOTE: command variables like "c1" have just been parsed out # by the call to expand_ih2 above, so this currently only catches # result vars like "r1" $sverr = "Attempt to set history value"; } if ($sverr ne "") { if ($BASIC_running) { print "*** "; } print "ERROR: $sverr '$setvar'\n"; $setvar = ""; $BASIC_running = 0; } ($lv_pt, $lv_val) = &eval_1($cmd); # evaluate if ($pr>1) { print "\n"; $lleft = 78; } if ($pr>1) { $hname = &define_R_hist($lv_pt, $lv_val); print "$hname = "; $lleft -= length("$hname = "); } if ($setvar ne "") { $vars_pt{$setvar} = $lv_pt; $vars_val{$setvar} = $lv_val; if ($pr>1) { print "$setvar: "; $lleft -= length("$setvar: "); } } if ($pr) { &prnt1($lv_pt, $lv_val, ($pr>1)); # 20080204: The following seems to produce lots of superfluous # blank lines, so I removed it and added a newline just before # the command prompt in the main loop. # print "\n"; $lleft = 78; } } # docmd4 handle special-variable setting (scale, debug, format) and passes # everything else along to docmd3 # # If $BASIC_running and $prflag is 0, printing is suppressed. sub docmd4 { local($cmd, $showeach, $prflag) = @_; local($cmdre); $cmdre = $cmd; $cmdre =~ s|[+*?\(\)\|\[\]\{\}\\]||g; if ($cmdre eq "") { # no-op; ignore } elsif ($cmd =~ m|^scale *\= *([0-9]+)|) { if ($inunix) { $newscale = $1; # Limit the scale. This is important -- lots of fundamental # assumptions in the PT algorithms fail if the number of # significant digits in the mantissa is greater than the # exponent value at which promotion occurs. Since promotion # happens at 1.0e300, we have to limit the scale to a bit less # than this. if ($newscale > 295) { print "I am setting scale to 295, the maximum allowed.\n"; $newscale = 295; } $newprim = ($newscale > 14) ? "bc" : "nat"; if ($newprim ne $curprim) { # close old functions if ($curprim eq "bc") { &fbc_close(); } else { &f64_close(); } # and open new ones if ($newprim eq "bc") { &fbc_init($newscale); } else { &f64_init(); } } elsif ($newprim eq "bc") { # re-scale &fbc_rescale($newscale); } $curprim = $newprim; $curscale = $newscale; &hi_init(); } else { print "Sorry precision (scale) can only be changed in UNIX.\n"; } } elsif ("quit" =~ m|^$cmdre|) { $going = 0; } else { &docmd3($cmd, $showeach, $prflag); } } # docmd6 handles '!!', '[ ]', and ';' and passes the rest to docmd4 sub docmd6 { local($cmd) = @_; local($c1, $c2, $c3, $showeach); # Map [] onto () (Note: [] are used internally within the eval # routines) $cmd =~ tr|\[\]|\(\)|; $cmd =~ tr/A-Z/a-z/; # Replace '!!' with previous typed command &dbg1("dc6a: $cmd", 1); while($cmd =~ m/ \!/) { $cmd =~ s/ \!/\!/; } $cmd = " " . $cmd; # %%% There's a problem here: if you type "4 x !!" it does not treat the 'x' # as multiplication $didit = 0; while($cmd =~ m/^(.*)([ \+\-\*\/\^\(\;])\!\!(.*)$/) { $c1 = $1; $c2 = $2; $c3 = $3; $cmd = "$c1$c2$lastc6in$c3"; if ($didit == 0) { print "!!: $lastc6in\n"; $didit = 1; } } $cmd =~ s/^ +//; $lastc6in = $cmd; &dbg1("dc6b: $cmd", 1); $showeach = ($cmd =~ m|\;|); $c2 = $cmd; while($c2 =~ m|^([^\;]*)\;(.*)$|) { $c1 = $1; $c2 = $2; &docmd4($c1, $showeach, 1); } if ($c2 ne "") { &docmd4($c2, $showeach, 1); } } # This routine handles the runtime BASIC commands, e.g. the no-ops like # PRINT and LET and flow control ops like FOR..NEXT sub docmd10 { local($cmd, $pc) = @_; local($c2, $v, $f1, $f2, $prflag); if (0) { } elsif ($cmd =~ m/^ *$/) { # no-op; skip } elsif ($cmd =~ m/^ *end *$/i) { $BASIC_running = 0; } elsif ($cmd =~ m/^for ($p_var) *= *([0-9]+) +to +([0-9]+)(;?)/i) { $v = $1; $f1 = $2; $f2 = $3; $prflag = $4; $v = lc($v); $prflag = ($prflag eq ';') ? 0 : 2; # print "FOR $v $f1 $f2\n"; $next_indent++; $next_var[$next_indent] = $v; $next_limit[$next_indent] = $f2; $next_loopto[$next_indent] = $pc; $next_prflag[$next_indent] = $prflag; # convert to LET $c2 = "$v = $f1"; &docmd4($c2,0, $prflag); } elsif ($cmd =~ m/^ *let +(.*)(;?)$/i) { $c2 = $1; $c2 = lc($c2); $prflag = ($c2 =~ m/;$/) ? 0 : 2; &docmd4($c2,0,$prflag); } elsif ($cmd =~ m/^ *next/i) { # convert to LET $v = $next_var[$next_indent]; $prflag = $next_prflag[$next_indent]; $c2 = "$v = $v + 1"; &docmd4($c2, 0, $prflag); # print "\n" if ($prflag); # test for overflow if ($vars_val{$v} > $next_limit[$next_indent]) { # we're done $next_indent--; } else { # jump back to loop return $NEXT{$next_loopto[$next_indent]}; } } elsif ($cmd =~ m/^ *print +(.*)$/i) { # PRINT: Loop on arguments, evaluate and print $c2 = $1; $c2 .= ","; # Ensure last piece is parsed while($c2 ne "") { # print "p1 c2 $c2\n"; $v = $f2; if ($c2 =~ m|^([^,;]*)([,;])(.*)$|) { $f1 = $1; $f2 = $2; $c2 = $3; # print "p2 f1 $f1\n"; $f1 =~ s|^ +||; $f1 =~ s| +$||; if ($f1 =~ m|^"([^"]+)"$|) { $f1 = $1; print $f1; } elsif ($f1 ne "") { $f1 =~ s/=//g; # Remove '=' if any to prevent assignment $f1 = lc($f1); &docmd3($f1, 0, 1); } if ($f2 eq ',') { print "\t"; } } else { $c2 = ""; } } if ($v eq ",") { print "\t"; } elsif ($v ne ";") { print "\n"; } } elsif ($cmd =~ m/^ *rem +/i) { # remark -- ignore } else { # Anything else is taken as being a normal command. $cmd = lc($cmd); $prflag = ($cmd =~ m/; *$/) ? 0 : 2; &docmd4($cmd, 0, $prflag); # print "\n" if ($prflag); } return 0; } # This routine scans through the BASIC program, setting up the pointers # used for the flow control operators, indenting for LIST command, and a # few simple error checks (like NEXT without FOR) sub BASIC_parse { local($ln, $prev, $cmd, $indent, $next_indent); local($v, $f1, $f2); # First we have to set up the next pointers $lleft = 78; $prev = ""; $next_indent = 0; $IF_level = 0; undef %IF_next; foreach $ln (sort {$a <=> $b} (keys %BASIC)) { $indent = $next_indent; if ($prev eq "") { $BASIC_first = $ln; } else { $NEXT{$prev} = $ln; } $prev = $ln; $cmd = $BASIC{$ln}; $cmd =~ s/^ *//; $BASIC{$ln} = $cmd; $wentup = 0; if ($cmd eq "") { # no-op; skip } elsif ($cmd =~ m/^for ($p_var) *= *([0-9]+) +to +([0-9]+)/i) { $next_indent++; } elsif ($cmd =~ m/^for /i) { print "*** syntax error in FOR statement, line $ln\n"; return -1; } elsif ($cmd =~ m/^ *next /i) { if ($indent > 0) { $indent--; } else { print "*** NEXT without FOR in line $ln\n"; return -1; } $next_indent = $indent; } elsif ($cmd =~ m/^if +/i) { $next_indent++; $IF_level++; $IF_ln = $ln; } elsif ($cmd =~ m/^else +/i) { if ($IF_level > 0) { $indent--; } else { print "*** ELSE without IF in line $ln\n"; return -1; } $IF_next{$IF_ln} = $ln; $IF_ln = $ln; } elsif ($cmd =~ m/^end +if *$/i) { if ($IF_level > 0) { $indent--; } else { print "*** END IF without IF in line $ln\n"; return -1; } $next_indent = $indent; $IF_next{$IF_ln} = $ln; $IF_next{$ln} = 0; $IF_level--; } elsif ($cmd =~ m/^ *rem +/i) { # remark -- ignore } $BASIC_indent{$ln} = $indent; } return 0; } sub BASIC_run { local($ln, $prev); local($pc, $cmd, $err, $next); $err = &BASIC_parse(); if ($err) { return; } $pc = $BASIC_first; $BASIC_running = 1; while($BASIC_running) { # Fetch instruction $cmd = $BASIC{$pc}; # Execute command $next = &docmd10($cmd, $pc); # If we're still running, test for end of file if ($BASIC_running) { if ($next != 0) { $pc = $next; } else { $pc = $NEXT{$pc}; } if ($pc eq "") { print "*** REACHED END OF FILE WITHOUT 'END' STATEMENT\n"; $BASIC_running = 0; } } } } sub make_dir { local($dir) = @_; local($d1, $d2); # print "make_dir : '$dir'\n"; if (-d $dir) { # all set # print "'$dir' exists.\n"; } elsif (-e $dir) { print STDERR "ERR03: Can't create directory $dir, file already there\n"; exit -1; } else { if ($dir =~ m|^(.+)/([^/]+)$|) { $d1 = $1; $d2 = $2; # print "recursive md '$d1'\n"; &make_dir($d1); $dir = $d2; } # print "mkdir $dir\n"; system("mkdir $dir"); } # At this point either the directory already exists or we just created it. # print "chdir '$dir'\n"; chdir($dir); } sub getprogname { local($pn, $cmdname) = @_; $pn =~ s/^ *//; $pn =~ s/ *$//; if ($pn eq "") { print $cmdname; $pn = <>; chop $pn; $pn =~ s/^ *//; $pn =~ s/ *$//; if ($pn =~ m/^$p_var$/i) { # okay } else { print "*** Program name must be alphanumeric\n"; return ""; } } return $pn; } sub go_basic_dir { local($bd, $ans); $bd = "shared/proj/hypercalc"; chdir; if (-d $bd) { chdir($bd); } else { print "In order to use this feature I need to create a directory at '$hd/$bd'.\n"; print "Is this okay? (Y/N) : "; $ans = <>; chomp $ans; $ans =~ s/^[ \t]*//; if ($ans =~ m/^[yY]/) { &make_dir($bd); } else { print STDERR "Exiting.\n"; exit(1); } } } sub BASIC_save { local($pn) = @_; local($ln); if ($pn eq "") { $pn = $PROGNAME; } $pn = &getprogname($pn, "Name to save as: "); if ($pn eq "") { return; } $PROGNAME = $pn; &go_basic_dir(); open(SAVE, "> $pn.basic"); foreach $ln (sort {$a <=> $b} (keys %BASIC)) { print SAVE "$ln $BASIC{$ln}\n"; } close SAVE; } sub BASIC_cat { local($l); &go_basic_dir(); open(CAT, "ls |"); while ($l = ) { chop $l; if ($l =~ m/\.basic$/) { print "$l\n"; } } close CAT; } sub BASIC_old { local($pn) = @_; local($l, $ln, $cmd); if ($pn eq "") { $pn = $PROGNAME; } $pn = &getprogname($pn, "Name of program to load: "); if ($pn eq "") { return; } $PROGNAME = $pn; &go_basic_dir(); if (-e "$pn.basic") { undef %BASIC; open(OLD, "$pn.basic"); while($l = ) { chomp $l; $l =~ s/\t/ /g; $l =~ s/^ +//; $l =~ s/ +$//; if ($l eq '') { } elsif ($l =~ m/^([0-9]+) (.*)$/) { $ln = $1; $cmd = $2; $BASIC{$ln} = $cmd; } else { print "*** syntax error: '$l'\n"; } } close OLD; } else { print "*** OLD: '$pn' not found.\n"; } } sub mainloop { local($lv_pt, $lv_val); $lleft = 78; $lv_pt = 0; $lv_val = 0; print "$GPLNOTICE$tagline\n"; $dbg1flag = 0; $going = 1; $ihnum = 1; while($going) { print "\nC$ihnum = "; if ($cmdline eq "") { $cmd = ; chop $cmd; } else { $cmd = $cmdline; $cmdline = ""; print "$cmd\n"; } if ($inunix) { $cmd = &fixerase($cmd); } print HC_LOG "I$ihnum = $cmd\n"; $cmd =~ s|^ +||; $cmd =~ s|^\?|help |; $cmdre = $cmd; $cmdre =~ s|[+*?\(\)\|\[\]\{\}\\]||g; if ($cmdre eq "") { # nothing } elsif ($cmd =~ m/^ *help +gpl/i) { open(HELPOUT, "| more"); print HELPOUT $GPL_VERSION_2; close HELPOUT; } elsif ($cmd =~ m/^ *help +ops/i) { print " Hypercalc operators, listed in precedence order: op assoc description ----- ----- --------------------------------- () parentheses around subexpression ! factorial or gamma(x+1) - negation, as in 'sin(-x)' ^ right exponent x * / left multiply (x = *) and divide + - left add and subtract = assign value to variable "; } elsif ($cmd =~ m/^ *help +fun/i) { print " Hypercalc functions: exp(x) e to the power of x arctan(x) arctangent of x in radians sqrt(x) square root of x ln(x) log base e of x log(x) log base 10 of x logn(b,x) log base b of x root(r,x) r-th root of x "; } elsif ($cmd =~ m/^ *help +var/i) { print " Hypercalc has three predefined constants: pi 3.141592653589... e 2.718281828459... phi (1+sqrt(5))/2, the 'golden ratio' To define your own constants, type something like this: c = 2.99 * 10^8 Constants can be changed, and used as variables in iterated commands (like the Newton's method example, see 'help syn') and in BASIC programs (see 'help basic') The following special variables can be set: scale Number of digits to use in calculations (normally 14) Example: scale=100 format Printout format (several options, 'help format' for a list) Example: format = 2 The effect of setting a special variable inside a BASIC program is undefined. "; } elsif ($cmd =~ m/^ *help +format/i) { print " Values of 'format' special variable: 1 human-readable ASCII (the default at startup) 2 RHTF 3 HTML 4 Old-style numbers.rhtf label 5 Current numbers.rhtf label and link 6 Internal PT-VAL representation 7 Base 60 (HH:MM:SS.FRAC) "; } elsif ($cmd =~ m/^ *help +debug/i) { print " Debug bitmask fields: 1 function calls; eval_1; expand; docmd 2 I and R history; eval_2 4 open bc "; } elsif ($cmd =~ m/^ *help +basic/i) { print " BASIC programming commands: 10 Any input line starting with a number and a word, with no intervening symbol, is taken to be a BASIC input line and is added to the current program. A line can be any valid command, and if it ends in ';', the result will not be printed. for..next loops may be nested in a program. Put a ';' at the end of the for statement to prevent it from printing its value each time through. rem Comment end Stop program execution There are currently no conditional statements (if/then/else), subroutines or user-defined functions. list List the current program (will automatically indent loop bodies) cat List saved files new Erase the current program run Run the program save Save the current file old Load a saved program "; } elsif ($cmd =~ m/^ *help +syn/i) { print " Hypercalc syntax := [ ';' ] := '!!' | 'c' | [ '=' ] | 'quit' Here's an example. This computes the solution of x^x=10 using Newton's Method: c1: x = 2 c2: t = 10 c3: y = x^x c4: dy = y * (1 + ln(x)) c5: x = x + (t-y)/dy c6: c3 c7: c4;c5 c9: c3;c4;c5 c12: !! c15: !!;!! "; } elsif ($cmd =~ m|^help|) { print " Example Description 3 * 4 + 2 Multiply and add (answer 14) 2+3x4 Same answer, spaces are optional, x = * % + 2 Add 2 to previous answer 4-3-2 Gives -1 (left-to-right) 2^2^3 Gives 256 (right-to-left) 1.2 x 10 ^ 3456 Entering a large number 1.2e3456 Same number in scientific notation 2P100 Googolplex (P = 'powers of 10', '2P' -> 10^10^) 27! Factorial of 27 27!!! Same as ((27!)!)! c5 Repeat command C5 (as if you retyped it) c6;c7;c8 Repeat three typed commands in a row !! Repeat last typed line q quit debug 17 set debug bitmask For a complete list of operators, type 'help ops' For functions, 'help fun' and variables, 'help var' For a more on command syntax type 'help syn' For BASIC commands, 'help basic' For debugging features, 'help debug' "; } elsif ($cmd =~ m|^debug ([0-9]+)|) { $dbg1flag = $1; print "Debug flag set to $dbg1flag.\n"; } elsif ($cmd =~ m|^format *\= *([0-9]+)|) { print "Setting format.\n"; &format1($1); print "\n"; } elsif ($cmd =~ m/^ *([0-9]+) +(x *\=.*)$/) { # A number followed by "x =..." : This is a BASIC program line $ln = $1; $text = $2; $BASIC{$ln} = $text; } elsif ($cmd =~ m/^ *([0-9]+) +x[^a-z]/) { # A number followed by 'x': This is a calculation &docmd6($cmd); } elsif ($cmd =~ m/^ *([0-9]+) +([a-z].*)$/i) { # A number followed by a word: This is a BASIC program line $ln = $1; $text = $2; $BASIC{$ln} = $text; } elsif ($cmd =~ m/^ *([0-9]+) $/) { # A number only: Remove a line from the BASIC program $ln = $1; delete $BASIC{$ln}; } elsif ($cmd =~ m/^ *[0-9]+/) { # This is a calculation &docmd6($cmd); } elsif ($cmd =~ m/^ *([0-9]+) +(.*)$/) { # Input (or modify) a line for a BASIC program $ln = $1; $text = $2; $BASIC{$ln} = $text; } elsif ($cmd =~ m/^ *list *$/i) { # List the BASIC program $err = &BASIC_parse(); foreach $ln (sort {$a <=> $b} (keys %BASIC)) { print ("$ln " . (" " x $BASIC_indent{$ln}) . " $BASIC{$ln}\n"); } } elsif ($cmd =~ m/^ *cat *$/i) { # Catalog &BASIC_cat(); } elsif ($cmd =~ m/^ *new *$/i) { # NEW without name $PROGNAME = ""; undef %BASIC; } elsif ($cmd =~ m/^ *new +($p_var) *$/i) { # NEW with name $PROGNAME = $1; undef %BASIC; } elsif ($cmd =~ m/^ *new/i) { # NEW with illegal parameter print "*** NEW: Program name must be alphanumeric\n"; } elsif ($cmd =~ m/^ *run *$/i) { # Run the BASIC program &BASIC_run(); } elsif ($cmd =~ m/^ *save (.*)$/i) { &BASIC_save($1); } elsif ($cmd =~ m/^ *save *$/i) { &BASIC_save(""); } elsif ($cmd =~ m/^ *old (.*)$/i) { &BASIC_old($1); } elsif ($cmd =~ m/^ *old *$/i) { &BASIC_old(""); } else { # assignment and/or expression &docmd6($cmd); } } } $| = 1; $0 = "Hypercalc"; $cmdline = "@ARGV"; $inunix = (-d "/usr") + 0; $inwindows = ($ENV{"WINDIR"} ne "") + 0; $incygwin = (`uname` =~ m/CYGWIN/) + 0; if ($incygwin) { $inwindows = 0; } # Cygwin is UNIX-y enough for our needs. if ($inunix == $inwindows) { print ("Warning: Cannot determine if this is Windows or " . "MacOS X/Linux/UNIX/Cygwin (U=$inunix,W=$inwindows,C=$incygwin).\n" . "Please email the author, 'mrob at mrob.com' and he'll try to fix" . " this.\n\n"); } if ($inunix) { &seterase(); } chdir; # If the user has a ~/tmp directory we put our log there; otherwise no log. if (-d "tmp") { open(HC_LOG, "> tmp/hc_log.txt"); } else { open(HC_LOG, "> /dev/null"); } $bc_running = 0; # We start with a scale of 14 and use native floating-point primitives. $curscale = 14; $doround = 1; $curprim = "nat"; &f64_init(); &hi_init(); # Init the printing/formatting routines &format1(1); &mainloop(); if ($curprim eq "bc") { &fbc_close(); } else { &f64_close(); } close HC_LOG;