/*      $NetBSD: n_argred.S,v 1.10 2024/05/07 15:15:09 riastradh Exp $  */
/*
 * Copyright (c) 1985, 1993
 *      The Regents of the University of California.  All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of the University nor the names of its contributors
 *    may be used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 *      @(#)argred.s    8.1 (Berkeley) 6/4/93
 */

#include <machine/asm.h>

/*
 *  libm$argred implements Bob Corbett's argument reduction and
 *  libm$sincos implements Peter Tang's double precision sin/cos.
 *
 *  Note: The two entry points libm$argred and libm$sincos are meant
 *        to be used only by _sin, _cos and _tan.
 *
 * method: true range reduction to [-pi/4,pi/4], P. Tang  &  B. Corbett
 * S. McDonald, April 4,  1985
 */

        .hidden __libm_argred
ENTRY(__libm_argred, 0)
/*
 *  Compare the argument with the largest possible that can
 *  be reduced by table lookup.  %r3 := |x|  will be used in  table_lookup .
 */
        movd    %r0,%r3
        bgeq    abs1
        mnegd   %r3,%r3
abs1:
        cmpd    %r3,$0d+4.55530934770520019583e+01
        blss    small_arg
        jsb     trigred
        rsb
small_arg:
        jsb     table_lookup
        rsb
/*
 *  At this point,
 *         %r0  contains the quadrant number, 0, 1, 2, or 3;
 *      %r2/%r1  contains the reduced argument as a D-format number;
 *         %r3  contains a F-format extension to the reduced argument;
 *          %r4  contains a  0 or 1  corresponding to a  sin or cos  entry.
 */
END(__libm_argred)

        .hidden __libm_sincos
ENTRY(__libm_sincos, 0)
/*
 *  Compensate for a cosine entry by adding one to the quadrant number.
 */
        addl2   %r4,%r0
/*
 *  Polyd clobbers  %r5-%r0 ;  save  X  in  %r7/%r6 .
 *  This can be avoided by rewriting  trigred .
 */
        movd    %r1,%r6
/*
 *  Likewise, save  alpha  in  %r8 .
 *  This can be avoided by rewriting  trigred .
 */
        movf    %r3,%r8
/*
 *  Odd or even quadrant?  cosine if odd, sine otherwise.
 *  Save  floor(quadrant/2) in  %r9  ; it determines the final sign.
 */
        rotl    $-1,%r0,%r9
        blss    cosine
sine:
        muld2   %r1,%r1         # Xsq = X * X
        cmpw    $0x2480,%r1     # [zl] Xsq > 2^-56?
        blss    1f              # [zl] yes, go ahead and do polyd
        clrq    %r1             # [zl] work around 11/780 FPA polyd bug
1:
        polyd   %r1,$7,sin_coef # Q = P(Xsq) , of deg 7
        mulf3   $0f3.0,%r8,%r4  # beta = 3 * alpha
        mulf2   %r0,%r4         # beta = Q * beta
        addf2   %r8,%r4         # beta = alpha + beta
        muld2   %r6,%r0         # S(X) = X * Q
/*      cvtfd   %r4,%r4         ... %r5 = 0 after a polyd. */
        addd2   %r4,%r0         # S(X) = beta + S(X)
        addd2   %r6,%r0         # S(X) = X + S(X)
        jbr     done
cosine:
        muld2   %r6,%r6         # Xsq = X * X
        beql    zero_arg
        mulf2   %r1,%r8         # beta = X * alpha
        polyd   %r6,$7,cos_coef /* Q = P'(Xsq) , of deg 7 */
        subd3   %r0,%r8,%r0     # beta = beta - Q
        subw2   $0x80,%r6       # Xsq = Xsq / 2
        addd2   %r0,%r6         # Xsq = Xsq + beta
zero_arg:
        subd3   %r6,$0d1.0,%r0  # C(X) = 1 - Xsq
done:
        blbc    %r9,even
        mnegd   %r0,%r0
even:
        rsb
END(__libm_sincos)

#ifdef __ELF__
        .section .rodata
#else
        .text
#endif
        _ALIGN_TEXT

sin_coef:
        .double 0d-7.53080332264191085773e-13   # s7 = 2^-29 -1.a7f2504ffc49f8..
        .double 0d+1.60573519267703489121e-10   # s6 = 2^-21  1.611adaede473c8..
        .double 0d-2.50520965150706067211e-08   # s5 = 2^-1a -1.ae644921ed8382..
        .double 0d+2.75573191800593885716e-06   # s4 = 2^-13  1.71de3a4b884278..
        .double 0d-1.98412698411850507950e-04   # s3 = 2^-0d -1.a01a01a0125e7d..
        .double 0d+8.33333333333325688985e-03   # s2 = 2^-07  1.11111111110e50
        .double 0d-1.66666666666666664354e-01   # s1 = 2^-03 -1.55555555555554
        .double 0d+0.00000000000000000000e+00   # s0 = 0

cos_coef:
        .double 0d-1.13006966202629430300e-11   # s7 = 2^-25 -1.8D9BA04D1374BE..
        .double 0d+2.08746646574796004700e-09   # s6 = 2^-1D  1.1EE632650350BA..
        .double 0d-2.75573073031284417300e-07   # s5 = 2^-16 -1.27E4F31411719E..
        .double 0d+2.48015872682668025200e-05   # s4 = 2^-10  1.A01A0196B902E8..
        .double 0d-1.38888888888464709200e-03   # s3 = 2^-0A -1.6C16C16C11FACE..
        .double 0d+4.16666666666664761400e-02   # s2 = 2^-05  1.5555555555539E
        .double 0d+0.00000000000000000000e+00   # s1 = 0
        .double 0d+0.00000000000000000000e+00   # s0 = 0

/*
 *  Multiples of  pi/2  expressed as the sum of three doubles,
 *
 *  trailing:   n * pi/2 ,  n = 0, 1, 2, ..., 29
 *                      trailing[n] ,
 *
 *  middle:     n * pi/2 ,  n = 0, 1, 2, ..., 29
 *                      middle[n]   ,
 *
 *  leading:    n * pi/2 ,  n = 0, 1, 2, ..., 29
 *                      leading[n]  ,
 *
 *      where
 *              leading[n]  := (n * pi/2)  rounded,
 *              middle[n]   := (n * pi/2  -  leading[n])  rounded,
 *              trailing[n] := (( n * pi/2 - leading[n]) - middle[n])  rounded .
 */
trailing:
        .double 0d+0.00000000000000000000e+00   #  0 * pi/2  trailing
        .double 0d+4.33590506506189049611e-35   #  1 * pi/2  trailing
        .double 0d+8.67181013012378099223e-35   #  2 * pi/2  trailing
        .double 0d+1.30077151951856714215e-34   #  3 * pi/2  trailing
        .double 0d+1.73436202602475619845e-34   #  4 * pi/2  trailing
        .double 0d-1.68390735624352669192e-34   #  5 * pi/2  trailing
        .double 0d+2.60154303903713428430e-34   #  6 * pi/2  trailing
        .double 0d-8.16726343231148352150e-35   #  7 * pi/2  trailing
        .double 0d+3.46872405204951239689e-34   #  8 * pi/2  trailing
        .double 0d+3.90231455855570147991e-34   #  9 * pi/2  trailing
        .double 0d-3.36781471248705338384e-34   # 10 * pi/2  trailing
        .double 0d-1.06379439835298071785e-33   # 11 * pi/2  trailing
        .double 0d+5.20308607807426856861e-34   # 12 * pi/2  trailing
        .double 0d+5.63667658458045770509e-34   # 13 * pi/2  trailing
        .double 0d-1.63345268646229670430e-34   # 14 * pi/2  trailing
        .double 0d-1.19986217995610764801e-34   # 15 * pi/2  trailing
        .double 0d+6.93744810409902479378e-34   # 16 * pi/2  trailing
        .double 0d-8.03640094449267300110e-34   # 17 * pi/2  trailing
        .double 0d+7.80462911711140295982e-34   # 18 * pi/2  trailing
        .double 0d-7.16921993148029483506e-34   # 19 * pi/2  trailing
        .double 0d-6.73562942497410676769e-34   # 20 * pi/2  trailing
        .double 0d-6.30203891846791677593e-34   # 21 * pi/2  trailing
        .double 0d-2.12758879670596143570e-33   # 22 * pi/2  trailing
        .double 0d+2.53800212047402350390e-33   # 23 * pi/2  trailing
        .double 0d+1.04061721561485371372e-33   # 24 * pi/2  trailing
        .double 0d+6.11729905311472319056e-32   # 25 * pi/2  trailing
        .double 0d+1.12733531691609154102e-33   # 26 * pi/2  trailing
        .double 0d-3.70049587943078297272e-34   # 27 * pi/2  trailing
        .double 0d-3.26690537292459340860e-34   # 28 * pi/2  trailing
        .double 0d-1.14812616507957271361e-34   # 29 * pi/2  trailing

middle:
        .double 0d+0.00000000000000000000e+00   #  0 * pi/2  middle
        .double 0d+5.72118872610983179676e-18   #  1 * pi/2  middle
        .double 0d+1.14423774522196635935e-17   #  2 * pi/2  middle
        .double 0d-3.83475850529283316309e-17   #  3 * pi/2  middle
        .double 0d+2.28847549044393271871e-17   #  4 * pi/2  middle
        .double 0d-2.69052076007086676522e-17   #  5 * pi/2  middle
        .double 0d-7.66951701058566632618e-17   #  6 * pi/2  middle
        .double 0d-1.54628301484890040587e-17   #  7 * pi/2  middle
        .double 0d+4.57695098088786543741e-17   #  8 * pi/2  middle
        .double 0d+1.07001849766246313192e-16   #  9 * pi/2  middle
        .double 0d-5.38104152014173353044e-17   # 10 * pi/2  middle
        .double 0d-2.14622680169080983801e-16   # 11 * pi/2  middle
        .double 0d-1.53390340211713326524e-16   # 12 * pi/2  middle
        .double 0d-9.21580002543456677056e-17   # 13 * pi/2  middle
        .double 0d-3.09256602969780081173e-17   # 14 * pi/2  middle
        .double 0d+3.03066796603896507006e-17   # 15 * pi/2  middle
        .double 0d+9.15390196177573087482e-17   # 16 * pi/2  middle
        .double 0d+1.52771359575124969107e-16   # 17 * pi/2  middle
        .double 0d+2.14003699532492626384e-16   # 18 * pi/2  middle
        .double 0d-1.68853170360202329427e-16   # 19 * pi/2  middle
        .double 0d-1.07620830402834670609e-16   # 20 * pi/2  middle
        .double 0d+3.97700719404595604379e-16   # 21 * pi/2  middle
        .double 0d-4.29245360338161967602e-16   # 22 * pi/2  middle
        .double 0d-3.68013020380794313406e-16   # 23 * pi/2  middle
        .double 0d-3.06780680423426653047e-16   # 24 * pi/2  middle
        .double 0d-2.45548340466059054318e-16   # 25 * pi/2  middle
        .double 0d-1.84316000508691335411e-16   # 26 * pi/2  middle
        .double 0d-1.23083660551323675053e-16   # 27 * pi/2  middle
        .double 0d-6.18513205939560162346e-17   # 28 * pi/2  middle
        .double 0d-6.18980636588357585202e-19   # 29 * pi/2  middle

leading:
        .double 0d+0.00000000000000000000e+00   #  0 * pi/2  leading
        .double 0d+1.57079632679489661351e+00   #  1 * pi/2  leading
        .double 0d+3.14159265358979322702e+00   #  2 * pi/2  leading
        .double 0d+4.71238898038468989604e+00   #  3 * pi/2  leading
        .double 0d+6.28318530717958645404e+00   #  4 * pi/2  leading
        .double 0d+7.85398163397448312306e+00   #  5 * pi/2  leading
        .double 0d+9.42477796076937979208e+00   #  6 * pi/2  leading
        .double 0d+1.09955742875642763501e+01   #  7 * pi/2  leading
        .double 0d+1.25663706143591729081e+01   #  8 * pi/2  leading
        .double 0d+1.41371669411540694661e+01   #  9 * pi/2  leading
        .double 0d+1.57079632679489662461e+01   # 10 * pi/2  leading
        .double 0d+1.72787595947438630262e+01   # 11 * pi/2  leading
        .double 0d+1.88495559215387595842e+01   # 12 * pi/2  leading
        .double 0d+2.04203522483336561422e+01   # 13 * pi/2  leading
        .double 0d+2.19911485751285527002e+01   # 14 * pi/2  leading
        .double 0d+2.35619449019234492582e+01   # 15 * pi/2  leading
        .double 0d+2.51327412287183458162e+01   # 16 * pi/2  leading
        .double 0d+2.67035375555132423742e+01   # 17 * pi/2  leading
        .double 0d+2.82743338823081389322e+01   # 18 * pi/2  leading
        .double 0d+2.98451302091030359342e+01   # 19 * pi/2  leading
        .double 0d+3.14159265358979324922e+01   # 20 * pi/2  leading
        .double 0d+3.29867228626928286062e+01   # 21 * pi/2  leading
        .double 0d+3.45575191894877260523e+01   # 22 * pi/2  leading
        .double 0d+3.61283155162826226103e+01   # 23 * pi/2  leading
        .double 0d+3.76991118430775191683e+01   # 24 * pi/2  leading
        .double 0d+3.92699081698724157263e+01   # 25 * pi/2  leading
        .double 0d+4.08407044966673122843e+01   # 26 * pi/2  leading
        .double 0d+4.24115008234622088423e+01   # 27 * pi/2  leading
        .double 0d+4.39822971502571054003e+01   # 28 * pi/2  leading
        .double 0d+4.55530934770520019583e+01   # 29 * pi/2  leading

twoOverPi:
        .double 0d+6.36619772367581343076e-01

        .text
        _ALIGN_TEXT

table_lookup:
        muld3   %r3,twoOverPi,%r0
        cvtrdl  %r0,%r0                 # n = nearest int to ((2/pi)*|x|) rnded
        subd2   leading[%r0],%r3                # p = (|x| - leading n*pi/2) exactly
        subd3   middle[%r0],%r3,%r1     # q = (p - middle  n*pi/2) rounded
        subd2   %r1,%r3                 # r = (p - q)
        subd2   middle[%r0],%r3         # r =  r - middle  n*pi/2
        subd2   trailing[%r0],%r3               # r =  r - trailing n*pi/2  rounded
/*
 *  If the original argument was negative,
 *  negate the reduce argument and
 *  adjust the octant/quadrant number.
 */
        tstw    4(%ap)
        bgeq    abs2
        mnegf   %r1,%r1
        mnegf   %r3,%r3
/*      subb3   %r0,$8,%r0      ...used for  pi/4  reduction -S.McD */
        subb3   %r0,$4,%r0
abs2:
/*
 *  Clear all unneeded octant/quadrant bits.
 */
/*      bicb2   $0xf8,%r0       ...used for  pi/4  reduction -S.McD */
        bicb2   $0xfc,%r0
        rsb
/*
 *                                              p.0
 */
#ifdef __ELF__
        .section .rodata
#else
        .text
#endif
        _ALIGN_TEXT
/*
 * Only 256 (actually 225) bits of 2/pi are needed for VAX double
 * precision; this was determined by enumerating all the nearest
 * machine integer multiples of pi/2 using continued fractions.
 * (8a8d3673775b7ff7 required the most bits.)           -S.McD
 */
        .long   0
        .long   0
        .long   0xaef1586d
        .long   0x9458eaf7
        .long   0x10e4107f
        .long   0xd8a5664f
        .long   0x4d377036
        .long   0x09d5f47d
        .long   0x91054a7f
        .long   0xbe60db93
bits2opi:
        .long   0x00000028
        .long   0
/*
 *  Note: wherever you see the word `octant', read `quadrant'.
 *  Currently this code is set up for  pi/2  argument reduction.
 *  By uncommenting/commenting the appropriate lines, it will
 *  also serve as a  pi/4  argument reduction code.
 */
        .text

/*                                              p.1
 *  Trigred  preforms argument reduction
 *  for the trigonometric functions.  It
 *  takes one input argument, a D-format
 *  number in  %r1/%r0 .  The magnitude of
 *  the input argument must be greater
 *  than or equal to  1/2 .  Trigred produces
 *  three results:  the number of the octant
 *  occupied by the argument, the reduced
 *  argument, and an extension of the
 *  reduced argument.  The octant number is
 *  returned in  %r0 .  The reduced argument
 *  is returned as a D-format number in
 *  %r2/%r1 .  An 8 bit extension of the
 *  reduced argument is returned as an
 *  F-format number in %r3.
 *                                              p.2
 */
trigred:
/*
 *  Save the sign of the input argument.
 */
        movw    %r0,-(%sp)
/*
 *  Extract the exponent field.
 */
        extzv   $7,$7,%r0,%r2
/*
 *  Convert the fraction part of the input
 *  argument into a quadword integer.
 */
        bicw2   $0xff80,%r0
        bisb2   $0x80,%r0       # -S.McD
        rotl    $16,%r0,%r0
        rotl    $16,%r1,%r1
/*
 *  If  %r1  is negative, add  1  to  %r0 .  This
 *  adjustment is made so that the two's
 *  complement multiplications done later
 *  will produce unsigned results.
 */
        bgeq    posmid
        incl    %r0
posmid:
/*                                              p.3
 *
 *  Set  %r3  to the address of the first quadword
 *  used to obtain the needed portion of  2/pi .
 *  The address is longword aligned to ensure
 *  efficient access.
 */
        ashl    $-3,%r2,%r3
        bicb2   $3,%r3
        mnegl   %r3,%r3
        movab   bits2opi[%r3],%r3
/*
 *  Set  %r2  to the size of the shift needed to
 *  obtain the correct portion of  2/pi .
 */
        bicb2   $0xe0,%r2
/*                                              p.4
 *
 *  Move the needed  128  bits of  2/pi  into
 *  %r11 - %r8 .  Adjust the numbers to allow
 *  for unsigned multiplication.
 */
        ashq    %r2,(%r3),%r10

        subl2   $4,%r3
        ashq    %r2,(%r3),%r9
        bgeq    signoff1
        incl    %r11
signoff1:
        subl2   $4,%r3
        ashq    %r2,(%r3),%r8
        bgeq    signoff2
        incl    %r10
signoff2:
        subl2   $4,%r3
        ashq    %r2,(%r3),%r7
        bgeq    signoff3
        incl    %r9
signoff3:
/*                                              p.5
 *
 *  Multiply the contents of  %r0/%r1  by the
 *  slice of  2/pi  in  %r11 - %r8 .
 */
        emul    %r0,%r8,$0,%r4
        emul    %r0,%r9,%r5,%r5
        emul    %r0,%r10,%r6,%r6

        emul    %r1,%r8,$0,%r7
        emul    %r1,%r9,%r8,%r8
        emul    %r1,%r10,%r9,%r9
        emul    %r1,%r11,%r10,%r10

        addl2   %r4,%r8
        adwc    %r5,%r9
        adwc    %r6,%r10
/*                                              p.6
 *
 *  If there are more than five leading zeros
 *  after the first two quotient bits or if there
 *  are more than five leading ones after the first
 *  two quotient bits, generate more fraction bits.
 *  Otherwise, branch to code to produce the result.
 */
        bicl3   $0xc1ffffff,%r10,%r4
        beql    more1
        cmpl    $0x3e000000,%r4
        bneq    result
more1:
/*                                              p.7
 *
 *  generate another  32  result bits.
 */
        subl2   $4,%r3
        ashq    %r2,(%r3),%r5
        bgeq    signoff4

        emul    %r1,%r6,$0,%r4
        addl2   %r1,%r5
        emul    %r0,%r6,%r5,%r5
        addl2   %r0,%r6
        jbr     addbits1

signoff4:
        emul    %r1,%r6,$0,%r4
        emul    %r0,%r6,%r5,%r5

addbits1:
        addl2   %r5,%r7
        adwc    %r6,%r8
        adwc    $0,%r9
        adwc    $0,%r10
/*                                              p.8
 *
 *  Check for massive cancellation.
 */
        bicl3   $0xc0000000,%r10,%r6
/*      bneq    more2                   -S.McD  Test was backwards */
        beql    more2
        cmpl    $0x3fffffff,%r6
        bneq    result
more2:
/*                                              p.9
 *
 *  If massive cancellation has occurred,
 *  generate another  24  result bits.
 *  Testing has shown there will always be
 *  enough bits after this point.
 */
        subl2   $4,%r3
        ashq    %r2,(%r3),%r5
        bgeq    signoff5

        emul    %r0,%r6,%r4,%r5
        addl2   %r0,%r6
        jbr     addbits2

signoff5:
        emul    %r0,%r6,%r4,%r5

addbits2:
        addl2   %r6,%r7
        adwc    $0,%r8
        adwc    $0,%r9
        adwc    $0,%r10
/*                                              p.10
 *
 *  The following code produces the reduced
 *  argument from the product bits contained
 *  in  %r10 - %r7 .
 */
result:
/*
 *  Extract the octant number from  %r10 .
 */
/*      extzv   $29,$3,%r10,%r0 ...used for  pi/4  reduction -S.McD */
        extzv   $30,$2,%r10,%r0
/*
 *  Clear the octant bits in  %r10 .
 */
/*      bicl2   $0xe0000000,%r10        ...used for  pi/4  reduction -S.McD */
        bicl2   $0xc0000000,%r10
/*
 *  Zero the sign flag.
 */
        clrl    %r5
/*                                              p.11
 *
 *  Check to see if the fraction is greater than
 *  or equal to one-half.  If it is, add one
 *  to the octant number, set the sign flag
 *  on, and replace the fraction with  1 minus
 *  the fraction.
 */
/*      bitl    $0x10000000,%r10                ...used for  pi/4  reduction -S.McD */
        bitl    $0x20000000,%r10
        beql    small
        incl    %r0
        incl    %r5
/*      subl3   %r10,$0x1fffffff,%r10   ...used for  pi/4  reduction -S.McD */
        subl3   %r10,$0x3fffffff,%r10
        mcoml   %r9,%r9
        mcoml   %r8,%r8
        mcoml   %r7,%r7
small:
/*                                              p.12
 *
 *  Test whether the first  29  bits of the ...used for  pi/4  reduction -S.McD
 *  Test whether the first  30  bits of the
 *  fraction are zero.
 */
        tstl    %r10
        beql    tiny
/*
 *  Find the position of the first one bit in  %r10 .
 */
        cvtld   %r10,%r1
        extzv   $7,$7,%r1,%r1
/*
 *  Compute the size of the shift needed.
 */
        subl3   %r1,$32,%r6
/*
 *  Shift up the high order  64  bits of the
 *  product.
 */
        ashq    %r6,%r9,%r10
        ashq    %r6,%r8,%r9
        jbr     mult
/*                                              p.13
 *
 *  Test to see if the sign bit of  %r9  is on.
 */
tiny:
        tstl    %r9
        bgeq    tinier
/*
 *  If it is, shift the product bits up  32  bits.
 */
        movl    $32,%r6
        movq    %r8,%r10
        tstl    %r10
        jbr     mult
/*                                              p.14
 *
 *  Test whether  %r9  is zero.  It is probably
 *  impossible for both  %r10  and  %r9  to be
 *  zero, but until proven to be so, the test
 *  must be made.
 */
tinier:
        beql    zero
/*
 *  Find the position of the first one bit in  %r9 .
 */
        cvtld   %r9,%r1
        extzv   $7,$7,%r1,%r1
/*
 *  Compute the size of the shift needed.
 */
        subl3   %r1,$32,%r1
        addl3   $32,%r1,%r6
/*
 *  Shift up the high order  64  bits of the
 *  product.
 */
        ashq    %r1,%r8,%r10
        ashq    %r1,%r7,%r9
        jbr     mult
/*                                              p.15
 *
 *  The following code sets the reduced
 *  argument to zero.
 */
zero:
        clrl    %r1
        clrl    %r2
        clrl    %r3
        jbr     return
/*                                              p.16
 *
 *  At this point,  %r0  contains the octant number,
 *  %r6  indicates the number of bits the fraction
 *  has been shifted,  %r5  indicates the sign of
 *  the fraction,  %r11/%r10  contain the high order
 *  64  bits of the fraction, and the condition
 *  codes indicate where the sign bit of  %r10
 *  is on.  The following code multiplies the
 *  fraction by  pi/2 .
 */
mult:
/*
 *  Save  %r11/%r10  in  %r4/%r1 .              -S.McD
 */
        movl    %r11,%r4
        movl    %r10,%r1
/*
 *  If the sign bit of  %r10  is on, add  1  to  %r11 .
 */
        bgeq    signoff6
        incl    %r11
signoff6:
/*                                              p.17
 *
 *  Move  pi/2  into  %r3/%r2 .
 */
        movq    $0xc90fdaa22168c235,%r2
/*
 *  Multiply the fraction by the portion of  pi/2
 *  in  %r2 .
 */
        emul    %r2,%r10,$0,%r7
        emul    %r2,%r11,%r8,%r7
/*
 *  Multiply the fraction by the portion of  pi/2
 *  in  %r3 .
 */
        emul    %r3,%r10,$0,%r9
        emul    %r3,%r11,%r10,%r10
/*
 *  Add the product bits together.
 */
        addl2   %r7,%r9
        adwc    %r8,%r10
        adwc    $0,%r11
/*
 *  Compensate for not sign extending  %r8  above.-S.McD
 */
        tstl    %r8
        bgeq    signoff6a
        decl    %r11
signoff6a:
/*
 *  Compensate for  %r11/%r10  being unsigned.  -S.McD
 */
        addl2   %r2,%r10
        adwc    %r3,%r11
/*
 *  Compensate for  %r3/%r2  being unsigned.    -S.McD
 */
        addl2   %r1,%r10
        adwc    %r4,%r11
/*                                              p.18
 *
 *  If the sign bit of  %r11  is zero, shift the
 *  product bits up one bit and increment  %r6 .
 */
        blss    signon
        incl    %r6
        ashq    $1,%r10,%r10
        tstl    %r9
        bgeq    signoff7
        incl    %r10
signoff7:
signon:
/*                                              p.19
 *
 *  Shift the  56  most significant product
 *  bits into  %r9/%r8 .  The sign extension
 *  will be handled later.
 */
        ashq    $-8,%r10,%r8
/*
 *  Convert the low order  8  bits of  %r10
 *  into an F-format number.
 */
        cvtbf   %r10,%r3
/*
 *  If the result of the conversion was
 *  negative, add  1  to  %r9/%r8 .
 */
        bgeq    chop
        incl    %r8
        adwc    $0,%r9
/*
 *  If  %r9  is now zero, branch to special
 *  code to handle that possibility.
 */
        beql    carryout
chop:
/*                                              p.20
 *
 *  Convert the number in  %r9/%r8  into
 *  D-format number in  %r2/%r1 .
 */
        rotl    $16,%r8,%r2
        rotl    $16,%r9,%r1
/*
 *  Set the exponent field to the appropriate
 *  value.  Note that the extra bits created by
 *  sign extension are now eliminated.
 */
        subw3   %r6,$131,%r6
        insv    %r6,$7,$9,%r1
/*
 *  Set the exponent field of the F-format
 *  number in  %r3  to the appropriate value.
 */
        tstf    %r3
        beql    return
/*      extzv   $7,$8,%r3,%r4   -S.McD */
        extzv   $7,$7,%r3,%r4
        addw2   %r4,%r6
/*      subw2   $217,%r6                -S.McD */
        subw2   $64,%r6
        insv    %r6,$7,$8,%r3
        jbr     return
/*                                              p.21
 *
 *  The following code generates the appropriate
 *  result for the unlikely possibility that
 *  rounding the number in  %r9/%r8  resulted in
 *  a carry out.
 */
carryout:
        clrl    %r1
        clrl    %r2
        subw3   %r6,$132,%r6
        insv    %r6,$7,$9,%r1
        tstf    %r3
        beql    return
        extzv   $7,$8,%r3,%r4
        addw2   %r4,%r6
        subw2   $218,%r6
        insv    %r6,$7,$8,%r3
/*                                              p.22
 *
 *  The following code makes an needed
 *  adjustments to the signs of the
 *  results or to the octant number, and
 *  then returns.
 */
return:
/*
 *  Test if the fraction was greater than or
 *  equal to  1/2 .  If so, negate the reduced
 *  argument.
 */
        blbc    %r5,signoff8
        mnegf   %r1,%r1
        mnegf   %r3,%r3
signoff8:
/*                                              p.23
 *
 *  If the original argument was negative,
 *  negate the reduce argument and
 *  adjust the octant number.
 */
        tstw    (%sp)+
        bgeq    signoff9
        mnegf   %r1,%r1
        mnegf   %r3,%r3
/*      subb3   %r0,$8,%r0      ...used for  pi/4  reduction -S.McD */
        subb3   %r0,$4,%r0
signoff9:
/*
 *  Clear all unneeded octant bits.
 *
 *      bicb2   $0xf8,%r0       ...used for  pi/4  reduction -S.McD */
        bicb2   $0xfc,%r0
/*
 *  Return.
 */
        rsb