diff options
Diffstat (limited to 'Build/source/texk/web2c/mplibdir/mpmath.w')
-rw-r--r-- | Build/source/texk/web2c/mplibdir/mpmath.w | 1226 |
1 files changed, 1226 insertions, 0 deletions
diff --git a/Build/source/texk/web2c/mplibdir/mpmath.w b/Build/source/texk/web2c/mplibdir/mpmath.w new file mode 100644 index 00000000000..92f6f4a75a0 --- /dev/null +++ b/Build/source/texk/web2c/mplibdir/mpmath.w @@ -0,0 +1,1226 @@ +% $Id $ +% +% Copyright 2008-2010 Taco Hoekwater. +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU Lesser General Public License as published by +% the Free Software Foundation, either version 3 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 Lesser General Public License for more details. +% +% You should have received a copy of the GNU Lesser General Public License +% along with this program. If not, see <http://www.gnu.org/licenses/>. +% +% TeX is a trademark of the American Mathematical Society. +% METAFONT is a trademark of Addison-Wesley Publishing Company. +% PostScript is a trademark of Adobe Systems Incorporated. + +% Here is TeX material that gets inserted after \input webmac + +\font\tenlogo=logo10 % font used for the METAFONT logo +\font\logos=logosl10 +\def\MF{{\tenlogo META}\-{\tenlogo FONT}} +\def\MP{{\tenlogo META}\-{\tenlogo POST}} + +\def\title{Reading TEX metrics files} +\pdfoutput=1 + +@ Introduction. + +@ +@d hlp1(A) mp->help_line[0]=A; } +@d hlp2(A,B) mp->help_line[1]=A; hlp1(B) +@d help1 { mp->help_ptr=1; hlp1 /* use this with one help line */ +@d help2 { mp->help_ptr=2; hlp2 /* use this with two help lines */ + +@ @c +#include <w2c/config.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include "mplib.h" +#include "mpmp.h" /* internal header */ +#include "mpmath.h" /* internal header */ +@h + +@ @c +@<Declarations@>; + +@ @(mpmath.h@>= +@<Types@>; +@<Internal library declarations@>; + +@ Currently empty +@<Declarations@>= + +@* Math initialization. + +@<Types@>= +typedef struct math_data { + scaled max_scaled_; + scaled one_third_max_scaled_; + scaled unity_; + scaled two_; + scaled three_; + scaled half_unit_; + scaled three_quarter_unit_; + fraction fraction_one_; + fraction fraction_half_; + fraction fraction_two_; + fraction fraction_three_; + fraction fraction_four_; + angle ninety_deg_; + angle one_eighty_deg_; + angle three_sixty_deg_; +} math_data; + +@ @<Internal library declarations@>= +void * mp_initialize_math (MP mp); +void mp_free_math (MP mp); + +@ @c +void * mp_initialize_math (MP mp) { + math_data *math = (math_data *)mp_xmalloc(mp,1,sizeof(math_data)); + /* here are the constants for |scaled| objects */ + math->max_scaled_ = EL_GORDO; + math->one_third_max_scaled_ = one_third_EL_GORDO; + math->unity_ = unity; + math->two_ = two; + math->three_ = three; + math->half_unit_ = half_unit; + math->three_quarter_unit_ = three_quarter_unit; + /* |fractions| */ + math->fraction_one_ = fraction_one; + math->fraction_half_ = fraction_half; + math->fraction_two_ = fraction_two; + math->fraction_three_ = fraction_three; + math->fraction_four_ = fraction_four; + /* |angles| */ + math->ninety_deg_ = ninety_deg; + math->one_eighty_deg_ = one_eighty_deg; + math->three_sixty_deg_ = three_sixty_deg; + return (void *)math; +} + +void mp_free_math (MP mp) { + free(mp->math); +} + +@ Fixed-point arithmetic is done on {\sl scaled integers\/} that are multiples +of $2^{-16}$. In other words, a binary point is assumed to be sixteen bit +positions from the right end of a binary computer word. + +@d unity 0x10000 /* $2^{16}$, represents 1.00000 */ +@d two (2*unity) /* $2^{17}$, represents 2.00000 */ +@d three (3*unity) /* $2^{17}+2^{16}$, represents 3.00000 */ +@d half_unit (unity/2) /* $2^{15}$, represents 0.50000 */ +@d three_quarter_unit (3*(unity/4)) /* $3\cdot2^{14}$, represents 0.75000 */ + +@d EL_GORDO 0x7fffffff /* $2^{31}-1$, the largest value that \MP\ likes */ +@d one_third_EL_GORDO 05252525252 + +@ One of \MP's most common operations is the calculation of +$\lfloor{a+b\over2}\rfloor$, +the midpoint of two given integers |a| and~|b|. The most decent way to do +this is to write `|(a+b)/2|'; but on many machines it is more efficient +to calculate `|(a+b)>>1|'. + +Therefore the midpoint operation will always be denoted by `|half(a+b)|' +in this program. If \MP\ is being implemented with languages that permit +binary shifting, the |half| macro should be changed to make this operation +as efficient as possible. Since some systems have shift operators that can +only be trusted to work on positive numbers, there is also a macro |halfp| +that is used only when the quantity being halved is known to be positive +or zero. + +@<Internal library declarations@>= +#define half(A) ((A) / 2) +#define halfp(A) (integer)((unsigned)(A) >> 1) + +@ Todo: Here are some compilation tricks for problems to be sorted out later + +@<Internal library declarations@>= +#define integer_as_fraction(A) (fraction)(A) + + +@ Here is a procedure analogous to |print_int|. If the output +of this procedure is subsequently read by \MP\ and converted by the +|round_decimals| routine above, it turns out that the original value will +be reproduced exactly. A decimal point is printed only if the value is +not an integer. If there is more than one way to print the result with +the optimum number of digits following the decimal point, the closest +possible value is given. + +The invariant relation in the \&{repeat} loop is that a sequence of +decimal digits yet to be printed will yield the original number if and only if +they form a fraction~$f$ in the range $s-\delta\L10\cdot2^{16}f<s$. +We can stop if and only if $f=0$ satisfies this condition; the loop will +terminate before $s$ can possibly become zero. + +@<Internal library declarations@>= +void mp_print_scaled (MP mp, scaled s); + +@ @c +void mp_print_scaled (MP mp, scaled s) { /* prints scaled real, rounded to five digits */ + scaled delta; /* amount of allowable inaccuracy */ + if (s < 0) { + mp_print_char (mp, xord ('-')); + s = -s; /* print the sign, if negative */ + } + mp_print_int (mp, s / unity); /* print the integer part */ + s = 10 * (s % unity) + 5; + if (s != 5) { + delta = 10; + mp_print_char (mp, xord ('.')); + do { + if (delta > unity) + s = s + 0100000 - (delta / 2); /* round the final digit */ + mp_print_char (mp, xord ('0' + (s / unity))); + s = 10 * (s % unity); + delta = delta * 10; + } while (s > delta); + } +} + +@ Addition is not always checked to make sure that it doesn't overflow, +but in places where overflow isn't too unlikely the |slow_add| routine +is used. + +@<Internal library declarations@>= +integer mp_slow_add (MP mp, integer x, integer y); + +@ @c +integer mp_slow_add (MP mp, integer x, integer y) { + if (x >= 0) { + if (y <= EL_GORDO - x) { + return x + y; + } else { + mp->arith_error = true; + return EL_GORDO; + } + } else if (-y <= EL_GORDO + x) { + return x + y; + } else { + mp->arith_error = true; + return -EL_GORDO; + } +} + +@ The |make_fraction| routine produces the |fraction| equivalent of +|p/q|, given integers |p| and~|q|; it computes the integer +$f=\lfloor2^{28}p/q+{1\over2}\rfloor$, when $p$ and $q$ are +positive. If |p| and |q| are both of the same scaled type |t|, +the ``type relation'' |make_fraction(t,t)=fraction| is valid; +and it's also possible to use the subroutine ``backwards,'' using +the relation |make_fraction(t,fraction)=t| between scaled types. + +If the result would have magnitude $2^{31}$ or more, |make_fraction| +sets |arith_error:=true|. Most of \MP's internal computations have +been designed to avoid this sort of error. + +If this subroutine were programmed in assembly language on a typical +machine, we could simply compute |(@t$2^{28}$@>*p)div q|, since a +double-precision product can often be input to a fixed-point division +instruction. But when we are restricted to int-eger arithmetic it +is necessary either to resort to multiple-precision maneuvering +or to use a simple but slow iteration. The multiple-precision technique +would be about three times faster than the code adopted here, but it +would be comparatively long and tricky, involving about sixteen +additional multiplications and divisions. + +This operation is part of \MP's ``inner loop''; indeed, it will +consume nearly 10\pct! of the running time (exclusive of input and output) +if the code below is left unchanged. A machine-dependent recoding +will therefore make \MP\ run faster. The present implementation +is highly portable, but slow; it avoids multiplication and division +except in the initial stage. System wizards should be careful to +replace it with a routine that is guaranteed to produce identical +results in all cases. +@^system dependencies@> + +As noted below, a few more routines should also be replaced by machine-dependent +code, for efficiency. But when a procedure is not part of the ``inner loop,'' +such changes aren't advisable; simplicity and robustness are +preferable to trickery, unless the cost is too high. +@^inner loop@> + +@<Internal library declarations@>= +fraction mp_make_fraction (MP mp, integer p, integer q); + +@ We need these preprocessor values + +@d TWEXP31 2147483648.0 +@d TWEXP28 268435456.0 +@d TWEXP16 65536.0 +@d TWEXP_16 (1.0/65536.0) +@d TWEXP_28 (1.0/268435456.0) + + +@c +fraction mp_make_fraction (MP mp, integer p, integer q) { + fraction i; + if (q == 0) + mp_confusion (mp, "/"); +@:this can't happen /}{\quad \./@> { + register double d; + d = TWEXP28 * (double) p / (double) q; + if ((p ^ q) >= 0) { + d += 0.5; + if (d >= TWEXP31) { + mp->arith_error = true; + return EL_GORDO; + } + i = (integer) d; + if (d == (double) i && (((q > 0 ? -q : q) & 077777) + * (((i & 037777) << 1) - 1) & 04000) != 0) + --i; + } else { + d -= 0.5; + if (d <= -TWEXP31) { + mp->arith_error = true; + return -EL_GORDO; + } + i = (integer) d; + if (d == (double) i && (((q > 0 ? q : -q) & 077777) + * (((i & 037777) << 1) + 1) & 04000) != 0) + ++i; + } + } + return i; +} + + +@ The dual of |make_fraction| is |take_fraction|, which multiplies a +given integer~|q| by a fraction~|f|. When the operands are positive, it +computes $p=\lfloor qf/2^{28}+{1\over2}\rfloor$, a symmetric function +of |q| and~|f|. + +This routine is even more ``inner loopy'' than |make_fraction|; +the present implementation consumes almost 20\pct! of \MP's computation +time during typical jobs, so a machine-language substitute is advisable. +@^inner loop@> @^system dependencies@> + +@<Internal library declarations@>= +integer mp_take_fraction (MP mp, integer q, fraction f); + +@ @c +integer mp_take_fraction (MP mp, integer p, fraction q) { + register double d; + register integer i; + d = (double) p *(double) q *TWEXP_28; + if ((p ^ q) >= 0) { + d += 0.5; + if (d >= TWEXP31) { + if (d != TWEXP31 || (((p & 077777) * (q & 077777)) & 040000) == 0) + mp->arith_error = true; + return EL_GORDO; + } + i = (integer) d; + if (d == (double) i && (((p & 077777) * (q & 077777)) & 040000) != 0) + --i; + } else { + d -= 0.5; + if (d <= -TWEXP31) { + if (d != -TWEXP31 || ((-(p & 077777) * (q & 077777)) & 040000) == 0) + mp->arith_error = true; + return -EL_GORDO; + } + i = (integer) d; + if (d == (double) i && ((-(p & 077777) * (q & 077777)) & 040000) != 0) + ++i; + } + return i; +} + + +@ When we want to multiply something by a |scaled| quantity, we use a scheme +analogous to |take_fraction| but with a different scaling. +Given positive operands, |take_scaled| +computes the quantity $p=\lfloor qf/2^{16}+{1\over2}\rfloor$. + +Once again it is a good idea to use a machine-language replacement if +possible; otherwise |take_scaled| will use more than 2\pct! of the running time +when the Computer Modern fonts are being generated. +@^inner loop@> + +@<Internal library declarations@>= +integer mp_take_scaled (MP mp, integer q, scaled f); + +@ @c +integer mp_take_scaled (MP mp, integer p, scaled q) { + register double d; + register integer i; + d = (double) p *(double) q *TWEXP_16; + if ((p ^ q) >= 0) { + d += 0.5; + if (d >= TWEXP31) { + if (d != TWEXP31 || (((p & 077777) * (q & 077777)) & 040000) == 0) + mp->arith_error = true; + return EL_GORDO; + } + i = (integer) d; + if (d == (double) i && (((p & 077777) * (q & 077777)) & 040000) != 0) + --i; + } else { + d -= 0.5; + if (d <= -TWEXP31) { + if (d != -TWEXP31 || ((-(p & 077777) * (q & 077777)) & 040000) == 0) + mp->arith_error = true; + return -EL_GORDO; + } + i = (integer) d; + if (d == (double) i && ((-(p & 077777) * (q & 077777)) & 040000) != 0) + ++i; + } + return i; +} + + +@ For completeness, there's also |make_scaled|, which computes a +quotient as a |scaled| number instead of as a |fraction|. +In other words, the result is $\lfloor2^{16}p/q+{1\over2}\rfloor$, if the +operands are positive. \ (This procedure is not used especially often, +so it is not part of \MP's inner loop.) + +@<Internal library ...@>= +scaled mp_make_scaled (MP mp, integer p, integer q); + +@ @c +scaled mp_make_scaled (MP mp, integer p, integer q) { + register integer i; + if (q == 0) + mp_confusion (mp, "/"); +@:this can't happen /}{\quad \./@> { + register double d; + d = TWEXP16 * (double) p / (double) q; + if ((p ^ q) >= 0) { + d += 0.5; + if (d >= TWEXP31) { + mp->arith_error = true; + return EL_GORDO; + } + i = (integer) d; + if (d == (double) i && (((q > 0 ? -q : q) & 077777) + * (((i & 037777) << 1) - 1) & 04000) != 0) + --i; + } else { + d -= 0.5; + if (d <= -TWEXP31) { + mp->arith_error = true; + return -EL_GORDO; + } + i = (integer) d; + if (d == (double) i && (((q > 0 ? q : -q) & 077777) + * (((i & 037777) << 1) + 1) & 04000) != 0) + ++i; + } + } + return i; +} + +@ The following function divides |s| by |m|. |dd| is number of decimal digits. + +@<Internal library ...@>= +scaled mp_divide_scaled (MP mp, scaled s, scaled m, integer dd); + +@ @c +scaled mp_divide_scaled (MP mp, scaled s, scaled m, integer dd) { + scaled q, r; + integer sign, i; + sign = 1; + if (s < 0) { + sign = -sign; + s = -s; + } + if (m < 0) { + sign = -sign; + m = -m; + } + if (m == 0) + mp_confusion (mp, "arithmetic: divided by zero"); + else if (m >= (EL_GORDO / 10)) + mp_confusion (mp, "arithmetic: number too big"); + q = s / m; + r = s % m; + for (i = 1; i <= dd; i++) { + q = 10 * q + (10 * r) / m; + r = (10 * r) % m; + } + if (2 * r >= m) { + q++; + r = r - m; + } + mp->scaled_out = sign * (s - (r / mp->ten_pow[dd])); + return (sign * q); +} + + +@ The following function is used to create a scaled integer from a given decimal +fraction $(.d_0d_1\ldots d_{k-1})$, where |0<=k<=17|. + +@<Internal library declarations@>= +scaled mp_round_decimals (MP mp, unsigned char *b, quarterword k); + +@ @c +scaled mp_round_decimals (MP mp, unsigned char *b, quarterword k) { + /* converts a decimal fraction */ + unsigned a = 0; /* the accumulator */ + int l = 0; + for ( l = k-1; l >= 0; l-- ) { + if (l<16) /* digits for |k>=17| cannot affect the result */ + a = (a + (unsigned) (*(b+l) - '0') * two) / 10; + } + return (scaled) halfp (a + 1); +} + +@ The |scaled| quantities in \MP\ programs are generally supposed to be +less than $2^{12}$ in absolute value, so \MP\ does much of its internal +arithmetic with 28~significant bits of precision. A |fraction| denotes +a scaled integer whose binary point is assumed to be 28 bit positions +from the right. + +@d fraction_half 01000000000 /* $2^{27}$, represents 0.50000000 */ +@d fraction_one 02000000000 /* $2^{28}$, represents 1.00000000 */ +@d fraction_two 04000000000 /* $2^{29}$, represents 2.00000000 */ +@d fraction_three 06000000000 /* $3\cdot2^{28}$, represents 3.00000000 */ +@d fraction_four 010000000000 /* $2^{30}$, represents 4.00000000 */ + +@ Here is a typical example of how the routines above can be used. +It computes the function +$${1\over3\tau}f(\theta,\phi)= +{\tau^{-1}\bigl(2+\sqrt2\,(\sin\theta-{1\over16}\sin\phi) + (\sin\phi-{1\over16}\sin\theta)(\cos\theta-\cos\phi)\bigr)\over +3\,\bigl(1+{1\over2}(\sqrt5-1)\cos\theta+{1\over2}(3-\sqrt5\,)\cos\phi\bigr)},$$ +where $\tau$ is a |scaled| ``tension'' parameter. This is \MP's magic +fudge factor for placing the first control point of a curve that starts +at an angle $\theta$ and ends at an angle $\phi$ from the straight path. +(Actually, if the stated quantity exceeds 4, \MP\ reduces it to~4.) + +The trigonometric quantity to be multiplied by $\sqrt2$ is less than $\sqrt2$. +(It's a sum of eight terms whose absolute values can be bounded using +relations such as $\sin\theta\cos\theta\L{1\over2}$.) Thus the numerator +is positive; and since the tension $\tau$ is constrained to be at least +$3\over4$, the numerator is less than $16\over3$. The denominator is +nonnegative and at most~6. Hence the fixed-point calculations below +are guaranteed to stay within the bounds of a 32-bit computer word. + +The angles $\theta$ and $\phi$ are given implicitly in terms of |fraction| +arguments |st|, |ct|, |sf|, and |cf|, representing $\sin\theta$, $\cos\theta$, +$\sin\phi$, and $\cos\phi$, respectively. + +@<Internal library declarations@>= +fraction mp_velocity (MP mp, fraction st, fraction ct, fraction sf, + fraction cf, scaled t); + +@ @c +fraction mp_velocity (MP mp, fraction st, fraction ct, fraction sf, + fraction cf, scaled t) { + integer acc, num, denom; /* registers for intermediate calculations */ + acc = mp_take_fraction (mp, st - (sf / 16), sf - (st / 16)); + acc = mp_take_fraction (mp, acc, ct - cf); + num = fraction_two + mp_take_fraction (mp, acc, 379625062); + /* $2^{28}\sqrt2\approx379625062.497$ */ + denom = + fraction_three + mp_take_fraction (mp, ct, + 497706707) + mp_take_fraction (mp, cf, + 307599661); + /* $3\cdot2^{27}\cdot(\sqrt5-1)\approx497706706.78$ and + $3\cdot2^{27}\cdot(3-\sqrt5\,)\approx307599661.22$ */ + if (t != unity) + num = mp_make_scaled (mp, num, t); + /* |make_scaled(fraction,scaled)=fraction| */ + if (num / 4 >= denom) + return fraction_four; + else + return mp_make_fraction (mp, num, denom); +} + + +@ The following somewhat different subroutine tests rigorously if $ab$ is +greater than, equal to, or less than~$cd$, +given integers $(a,b,c,d)$. In most cases a quick decision is reached. +The result is $+1$, 0, or~$-1$ in the three respective cases. + +@<Internal library declarations@>= +integer mp_ab_vs_cd (MP mp, integer a, integer b, integer c, integer d); + +@ @c +integer mp_ab_vs_cd (MP mp, integer a, integer b, integer c, integer d) { + integer q, r; /* temporary registers */ + (void)mp; + @<Reduce to the case that |a,c>=0|, |b,d>0|@>; + while (1) { + q = a / d; + r = c / b; + if (q != r) + return (q > r ? 1 : -1); + q = a % d; + r = c % b; + if (r == 0) + return (q ? 1 : 0); + if (q == 0) + return -1; + a = b; + b = q; + c = d; + d = r; + } /* now |a>d>0| and |c>b>0| */ +} + + +@ @<Reduce to the case that |a...@>= +if (a < 0) { + a = -a; + b = -b; +}; +if (c < 0) { + c = -c; + d = -d; +}; +if (d <= 0) { + if (b >= 0) { + if ((a == 0 || b == 0) && (c == 0 || d == 0)) + return 0; + else + return 1; + } + if (d == 0) + return (a == 0 ? 0 : -1); + q = a; + a = c; + c = q; + q = -b; + b = -d; + d = q; +} else if (b <= 0) { + if (b < 0) + if (a > 0) + return -1; + return (c == 0 ? 0 : -1); +} + +@ We conclude this set of elementary routines with some simple rounding +and truncation operations. + +@ |floor_scaled| floors a |scaled| +@<Internal library declarations@>= +#define mp_floor_scaled(M,i) ((i)&(-65536)) + +@ |round_unscaled| rounds a |scaled| and converts it to |int| +@<Internal library declarations@>= +#define mp_round_unscaled(M,x) (x>=0100000 ? 1+((x-0100000) / 0200000) \ + : ( x>=-0100000 ? 0 : -(1+((-(x+1)-0100000) / 0200000)))) + +@ |round_fraction| rounds a |fraction| and converts it to |scaled| +@<Internal library declarations@>= +#define mp_round_fraction(M,x) (x>=2048 ? 1+((x-2048) / 4096) \ + : ( x>=-2048 ? 0 : -(1+((-(x+1)-2048) / 4096)))) + + + +@* Algebraic and transcendental functions. +\MP\ computes all of the necessary special functions from scratch, without +relying on |real| arithmetic or system subroutines for sines, cosines, etc. + +@ To get the square root of a |scaled| number |x|, we want to calculate +$s=\lfloor 2^8\!\sqrt x +{1\over2}\rfloor$. If $x>0$, this is the unique +integer such that $2^{16}x-s\L s^2<2^{16}x+s$. The following subroutine +determines $s$ by an iterative method that maintains the invariant +relations $x=2^{46-2k}x_0\bmod 2^{30}$, $0<y=\lfloor 2^{16-2k}x_0\rfloor +-s^2+s\L q=2s$, where $x_0$ is the initial value of $x$. The value of~$y$ +might, however, be zero at the start of the first iteration. + +@<Internal library declarations@>= +scaled mp_square_rt (MP mp, scaled x); + +@ @c +scaled mp_square_rt (MP mp, scaled x) { + quarterword k; /* iteration control counter */ + integer y; /* register for intermediate calculations */ + integer q; /* register for intermediate calculations */ + if (x <= 0) { + @<Handle square root of zero or negative argument@>; + } else { + k = 23; + q = 2; + while (x < fraction_two) { /* i.e., |while x<@t$2^{29}$@>|\unskip */ + k--; + x = x + x + x + x; + } + if (x < fraction_four) + y = 0; + else { + x = x - fraction_four; + y = 1; + }; + do { + @<Decrease |k| by 1, maintaining the invariant + relations between |x|, |y|, and~|q|@>; + } while (k != 0); + return (scaled) (halfp (q)); + } +} + + +@ @<Handle square root of zero...@>= +{ + if (x < 0) { + mp_print_err (mp, "Square root of "); +@.Square root...replaced by 0@>; + mp_print_scaled (mp, x); + mp_print (mp, " has been replaced by 0"); + help2 ("Since I don't take square roots of negative numbers,", + "I'm zeroing this one. Proceed, with fingers crossed."); + mp_error (mp); + }; + return 0; +} + + +@ @<Decrease |k| by 1, maintaining...@>= +x += x; +y += y; +if (x >= fraction_four) { /* note that |fraction_four=@t$2^{30}$@>| */ + x = x - fraction_four; + y++; +}; +x += x; +y = y + y - q; +q += q; +if (x >= fraction_four) { + x = x - fraction_four; + y++; +}; +if (y > (int) q) { + y -= q; + q += 2; +} else if (y <= 0) { + q -= 2; + y += q; +}; +k-- + +@ Pythagorean addition $\psqrt{a^2+b^2}$ is implemented by an elegant +iterative scheme due to Cleve Moler and Donald Morrison [{\sl IBM Journal +@^Moler, Cleve Barry@> +@^Morrison, Donald Ross@> +of Research and Development\/ \bf27} (1983), 577--581]. It modifies |a| and~|b| +in such a way that their Pythagorean sum remains invariant, while the +smaller argument decreases. + +@<Internal library ...@>= +integer mp_pyth_add (MP mp, integer a, integer b); + + +@ @c +integer mp_pyth_add (MP mp, integer a, integer b) { + fraction r; /* register used to transform |a| and |b| */ + boolean big; /* is the result dangerously near $2^{31}$? */ + a = abs (a); + b = abs (b); + if (a < b) { + r = b; + b = a; + a = r; + }; /* now |0<=b<=a| */ + if (b > 0) { + if (a < fraction_two) { + big = false; + } else { + a = a / 4; + b = b / 4; + big = true; + }; /* we reduced the precision to avoid arithmetic overflow */ + @<Replace |a| by an approximation to $\psqrt{a^2+b^2}$@>; + if (big) { + if (a < fraction_two) { + a = a + a + a + a; + } else { + mp->arith_error = true; + a = EL_GORDO; + }; + } + } + return a; +} + + +@ The key idea here is to reflect the vector $(a,b)$ about the +line through $(a,b/2)$. + +@<Replace |a| by an approximation to $\psqrt{a^2+b^2}$@>= +while (1) { + r = mp_make_fraction (mp, b, a); + r = mp_take_fraction (mp, r, r); /* now $r\approx b^2/a^2$ */ + if (r == 0) + break; + r = mp_make_fraction (mp, r, fraction_four + r); + a = a + mp_take_fraction (mp, a + a, r); + b = mp_take_fraction (mp, b, r); +} + + +@ Here is a similar algorithm for $\psqrt{a^2-b^2}$. +It converges slowly when $b$ is near $a$, but otherwise it works fine. + +@<Internal library declarations@>= +integer mp_pyth_sub (MP mp, integer a, integer b); + +@ @c +integer mp_pyth_sub (MP mp, integer a, integer b) { + fraction r; /* register used to transform |a| and |b| */ + boolean big; /* is the input dangerously near $2^{31}$? */ + a = abs (a); + b = abs (b); + if (a <= b) { + @<Handle erroneous |pyth_sub| and set |a:=0|@>; + } else { + if (a < fraction_four) { + big = false; + } else { + a = (integer) halfp (a); + b = (integer) halfp (b); + big = true; + } + @<Replace |a| by an approximation to $\psqrt{a^2-b^2}$@>; + if (big) + a *= 2; + } + return a; +} + + +@ @<Replace |a| by an approximation to $\psqrt{a^2-b^2}$@>= +while (1) { + r = mp_make_fraction (mp, b, a); + r = mp_take_fraction (mp, r, r); /* now $r\approx b^2/a^2$ */ + if (r == 0) + break; + r = mp_make_fraction (mp, r, fraction_four - r); + a = a - mp_take_fraction (mp, a + a, r); + b = mp_take_fraction (mp, b, r); +} + + +@ @<Handle erroneous |pyth_sub| and set |a:=0|@>= +{ + if (a < b) { + mp_print_err (mp, "Pythagorean subtraction "); + mp_print_scaled (mp, a); + mp_print (mp, "+-+"); + mp_print_scaled (mp, b); + mp_print (mp, " has been replaced by 0"); +@.Pythagorean...@>; + help2 ("Since I don't take square roots of negative numbers,", + "I'm zeroing this one. Proceed, with fingers crossed."); + mp_error (mp); + } + a = 0; +} + + +@ The subroutines for logarithm and exponential involve two tables. +The first is simple: |two_to_the[k]| equals $2^k$. The second involves +a bit more calculation, which the author claims to have done correctly: +|spec_log[k]| is $2^{27}$ times $\ln\bigl(1/(1-2^{-k})\bigr)= +2^{-k}+{1\over2}2^{-2k}+{1\over3}2^{-3k}+\cdots\,$, rounded to the +nearest integer. + +@d two_to_the(A) (1<<(unsigned)(A)) + +@<Declarations@>= +static const integer spec_log[29] = { 0, /* special logarithms */ + 93032640, 38612034, 17922280, 8662214, 4261238, 2113709, + 1052693, 525315, 262400, 131136, 65552, 32772, 16385, + 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1, 1 +}; + + +@ Here is the routine that calculates $2^8$ times the natural logarithm +of a |scaled| quantity; it is an integer approximation to $2^{24}\ln(x/2^{16})$, +when |x| is a given positive integer. + +The method is based on exercise 1.2.2--25 in {\sl The Art of Computer +Programming\/}: During the main iteration we have $1\L 2^{-30}x<1/(1-2^{1-k})$, +and the logarithm of $2^{30}x$ remains to be added to an accumulator +register called~$y$. Three auxiliary bits of accuracy are retained in~$y$ +during the calculation, and sixteen auxiliary bits to extend |y| are +kept in~|z| during the initial argument reduction. (We add +$100\cdot2^{16}=6553600$ to~|z| and subtract 100 from~|y| so that |z| will +not become negative; also, the actual amount subtracted from~|y| is~96, +not~100, because we want to add~4 for rounding before the final division by~8.) + +@<Internal library declarations@>= +scaled mp_m_log (MP mp, scaled x); + +@ @c +scaled mp_m_log (MP mp, scaled x) { + integer y, z; /* auxiliary registers */ + integer k; /* iteration counter */ + if (x <= 0) { + @<Handle non-positive logarithm@>; + } else { + y = 1302456956 + 4 - 100; /* $14\times2^{27}\ln2\approx1302456956.421063$ */ + z = 27595 + 6553600; /* and $2^{16}\times .421063\approx 27595$ */ + while (x < fraction_four) { + x = 2*x; + y -= 93032639; + z -= 48782; + } /* $2^{27}\ln2\approx 93032639.74436163$ and $2^{16}\times.74436163\approx 48782$ */ + y = y + (z / unity); + k = 2; + while (x > fraction_four + 4) { + @<Increase |k| until |x| can be multiplied by a + factor of $2^{-k}$, and adjust $y$ accordingly@>; + } + return (y / 8); + } +} + + +@ @<Increase |k| until |x| can...@>= +{ + z = ((x - 1) / two_to_the (k)) + 1; /* $z=\lceil x/2^k\rceil$ */ + while (x < fraction_four + z) { + z = halfp (z + 1); + k++; + }; + y += spec_log[k]; + x -= z; +} + + +@ @<Handle non-positive logarithm@>= +{ + mp_print_err (mp, "Logarithm of "); +@.Logarithm...replaced by 0@>; + mp_print_scaled (mp, x); + mp_print (mp, " has been replaced by 0"); + help2 ("Since I don't take logs of non-positive numbers,", + "I'm zeroing this one. Proceed, with fingers crossed."); + mp_error (mp); + return 0; +} + + +@ Conversely, the exponential routine calculates $\exp(x/2^8)$, +when |x| is |scaled|. The result is an integer approximation to +$2^{16}\exp(x/2^{24})$, when |x| is regarded as an integer. + +@<Internal library declarations@>= +scaled mp_m_exp (MP mp, scaled x); + +@ @c +scaled mp_m_exp (MP mp, scaled x) { + quarterword k; /* loop control index */ + integer y, z; /* auxiliary registers */ + if (x > 174436200) { + /* $2^{24}\ln((2^{31}-1)/2^{16})\approx 174436199.51$ */ + mp->arith_error = true; + return EL_GORDO; + } else if (x < -197694359) { + /* $2^{24}\ln(2^{-1}/2^{16})\approx-197694359.45$ */ + return 0; + } else { + if (x <= 0) { + z = -8 * x; + y = 04000000; /* $y=2^{20}$ */ + } else { + if (x <= 127919879) { + z = 1023359037 - 8 * x; + /* $2^{27}\ln((2^{31}-1)/2^{20})\approx 1023359037.125$ */ + } else { + z = 8 * (174436200 - x); /* |z| is always nonnegative */ + } + y = EL_GORDO; + }; + @<Multiply |y| by $\exp(-z/2^{27})$@>; + if (x <= 127919879) + return ((y + 8) / 16); + else + return y; + } +} + + +@ The idea here is that subtracting |spec_log[k]| from |z| corresponds +to multiplying |y| by $1-2^{-k}$. + +A subtle point (which had to be checked) was that if $x=127919879$, the +value of~|y| will decrease so that |y+8| doesn't overflow. In fact, +$z$ will be 5 in this case, and |y| will decrease by~64 when |k=25| +and by~16 when |k=27|. + +@<Multiply |y| by...@>= +k = 1; +while (z > 0) { + while (z >= spec_log[k]) { + z -= spec_log[k]; + y = y - 1 - ((y - two_to_the (k - 1)) / two_to_the (k)); + } + k++; +} + +@ The trigonometric subroutines use an auxiliary table such that +|spec_atan[k]| contains an approximation to the |angle| whose tangent +is~$1/2^k$. $\arctan2^{-k}$ times $2^{20}\cdot180/\pi$ + +@<Declarations@>= +static const angle spec_atan[27] = { 0, 27855475, 14718068, 7471121, 3750058, + 1876857, 938658, 469357, 234682, 117342, 58671, 29335, 14668, 7334, 3667, + 1833, 917, 458, 229, 115, 57, 29, 14, 7, 4, 2, 1 +}; + + +@ Given integers |x| and |y|, not both zero, the |n_arg| function +returns the |angle| whose tangent points in the direction $(x,y)$. +This subroutine first determines the correct octant, then solves the +problem for |0<=y<=x|, then converts the result appropriately to +return an answer in the range |-one_eighty_deg<=@t$\theta$@><=one_eighty_deg|. +(The answer is |+one_eighty_deg| if |y=0| and |x<0|, but an answer of +|-one_eighty_deg| is possible if, for example, |y=-1| and $x=-2^{30}$.) + +The octants are represented in a ``Gray code,'' since that turns out +to be computationally simplest. + +@d negate_x 1 +@d negate_y 2 +@d switch_x_and_y 4 +@d first_octant 1 +@d second_octant (first_octant+switch_x_and_y) +@d third_octant (first_octant+switch_x_and_y+negate_x) +@d fourth_octant (first_octant+negate_x) +@d fifth_octant (first_octant+negate_x+negate_y) +@d sixth_octant (first_octant+switch_x_and_y+negate_x+negate_y) +@d seventh_octant (first_octant+switch_x_and_y+negate_y) +@d eighth_octant (first_octant+negate_y) + +@<Internal library declarations@>= +angle mp_n_arg (MP mp, integer x, integer y); + +@ @c +angle mp_n_arg (MP mp, integer x, integer y) { + angle z; /* auxiliary register */ + integer t; /* temporary storage */ + quarterword k; /* loop counter */ + int octant; /* octant code */ + if (x >= 0) { + octant = first_octant; + } else { + x = -x; + octant = first_octant + negate_x; + } + if (y < 0) { + y = -y; + octant = octant + negate_y; + } + if (x < y) { + t = y; + y = x; + x = t; + octant = octant + switch_x_and_y; + } + if (x == 0) { + @<Handle undefined arg@>; + } else { + @<Set variable |z| to the arg of $(x,y)$@>; + @<Return an appropriate answer based on |z| and |octant|@>; + } +} + + +@ @<Handle undefined arg@>= +{ + mp_print_err (mp, "angle(0,0) is taken as zero"); +@.angle(0,0)...zero@>; + help2 ("The `angle' between two identical points is undefined.", + "I'm zeroing this one. Proceed, with fingers crossed."); + mp_error (mp); + return 0; +} + + +@ @<Return an appropriate answer...@>= +switch (octant) { +case first_octant: + return z; +case second_octant: + return (ninety_deg - z); +case third_octant: + return (ninety_deg + z); +case fourth_octant: + return (one_eighty_deg - z); +case fifth_octant: + return (z - one_eighty_deg); +case sixth_octant: + return (-z - ninety_deg); +case seventh_octant: + return (z - ninety_deg); +case eighth_octant: + return (-z); +}; /* there are no other cases */ +return 0 + +@ At this point we have |x>=y>=0|, and |x>0|. The numbers are scaled up +or down until $2^{28}\L x<2^{29}$, so that accurate fixed-point calculations +will be made. + +@<Set variable |z| to the arg...@>= +while (x >= fraction_two) { + x = halfp (x); + y = halfp (y); +} +z = 0; +if (y > 0) { + while (x < fraction_one) { + x += x; + y += y; + }; + @<Increase |z| to the arg of $(x,y)$@>; +} + +@ During the calculations of this section, variables |x| and~|y| +represent actual coordinates $(x,2^{-k}y)$. We will maintain the +condition |x>=y|, so that the tangent will be at most $2^{-k}$. +If $x<2y$, the tangent is greater than $2^{-k-1}$. The transformation +$(a,b)\mapsto(a+b\tan\phi,b-a\tan\phi)$ replaces $(a,b)$ by +coordinates whose angle has decreased by~$\phi$; in the special case +$a=x$, $b=2^{-k}y$, and $\tan\phi=2^{-k-1}$, this operation reduces +to the particularly simple iteration shown here. [Cf.~John E. Meggitt, +@^Meggitt, John E.@> +{\sl IBM Journal of Research and Development\/ \bf6} (1962), 210--226.] + +The initial value of |x| will be multiplied by at most +$(1+{1\over2})(1+{1\over8})(1+{1\over32})\cdots\approx 1.7584$; hence +there is no chance of integer overflow. + +@<Increase |z|...@>= +k = 0; +do { + y += y; + k++; + if (y > x) { + z = z + spec_atan[k]; + t = x; + x = x + (y / two_to_the (k + k)); + y = y - t; + }; +} while (k != 15); +do { + y += y; + k++; + if (y > x) { + z = z + spec_atan[k]; + y = y - x; + }; +} while (k != 26) + +@ Conversely, the |n_sin_cos| routine takes an |angle| and produces the sine +and cosine of that angle. The results of this routine are +stored in global integer variables |n_sin| and |n_cos|. + +@ Given an integer |z| that is $2^{20}$ times an angle $\theta$ in degrees, +the purpose of |n_sin_cos(z)| is to set +|x=@t$r\cos\theta$@>| and |y=@t$r\sin\theta$@>| (approximately), +for some rather large number~|r|. The maximum of |x| and |y| +will be between $2^{28}$ and $2^{30}$, so that there will be hardly +any loss of accuracy. Then |x| and~|y| are divided by~|r|. + +@d forty_five_deg 0264000000 /* $45\cdot2^{20}$, represents $45^\circ$ */ +@d ninety_deg 0550000000 /* $90\cdot2^{20}$, represents $90^\circ$ */ +@d one_eighty_deg 01320000000 /* $180\cdot2^{20}$, represents $180^\circ$ */ +@d three_sixty_deg 02640000000 /* $360\cdot2^{20}$, represents $360^\circ$ */ + +@d odd(A) ((A)%2==1) + +@<Internal library declarations@>= +void mp_n_sin_cos (MP mp, angle z, fraction *n_cos, fraction *n_sin); + +@ Compute a multiple of the sine and cosine + +@c +void mp_n_sin_cos (MP mp, angle z, fraction *n_cos, fraction *n_sin) { + quarterword k; /* loop control variable */ + int q; /* specifies the quadrant */ + fraction r; /* magnitude of |(x,y)| */ + integer x, y, t; /* temporary registers */ + while (z < 0) + z = z + three_sixty_deg; + z = z % three_sixty_deg; /* now |0<=z<three_sixty_deg| */ + q = z / forty_five_deg; + z = z % forty_five_deg; + x = fraction_one; + y = x; + if (!odd (q)) + z = forty_five_deg - z; + @<Subtract angle |z| from |(x,y)|@>; + @<Convert |(x,y)| to the octant determined by~|q|@>; + r = mp_pyth_add (mp, x, y); + *n_cos = mp_make_fraction (mp, x, r); + *n_sin = mp_make_fraction (mp, y, r); +} + + +@ In this case the octants are numbered sequentially. + +@<Convert |(x,...@>= +switch (q) { +case 0: + break; +case 1: + t = x; + x = y; + y = t; + break; +case 2: + t = x; + x = -y; + y = t; + break; +case 3: + x = -x; + break; +case 4: + x = -x; + y = -y; + break; +case 5: + t = x; + x = -y; + y = -t; + break; +case 6: + t = x; + x = y; + y = -t; + break; +case 7: + y = -y; + break; +} /* there are no other cases */ + + +@ The main iteration of |n_sin_cos| is similar to that of |n_arg| but +applied in reverse. The values of |spec_atan[k]| decrease slowly enough +that this loop is guaranteed to terminate before the (nonexistent) value +|spec_atan[27]| would be required. + +@<Subtract angle |z|...@>= +k = 1; +while (z > 0) { + if (z >= spec_atan[k]) { + z = z - spec_atan[k]; + t = x; + x = t + y / two_to_the (k); + y = y - t / two_to_the (k); + } + k++; +} +if (y < 0) + y = 0 /* this precaution may never be needed */ + + |