/**************************************************************************/
/*                                                                        */
/*                                 Sources                                */
/*                                from  the                               */
/*                           XSC SOFTWARE PACKAGE                         */
/*                                                                        */
/*                Copyright (c) 1997 University of Karlsruhe              */
/*                     Institute for Appl. Mathematics                    */
/*                                                                        */
/* This program is free software for non-commercial use. For details, see */
/* the file LICENSE.XSC, distributed with this software package.          */
/*                                                                        */
/* This program is distributed WITHOUT ANY WARRANTY. For details, see the */
/* file LICENSE.XSC, distributed with this software package.              */
/*                                                                        */
/**************************************************************************/

/* RCS $Id: i_sin.c,v 1.2 1997/09/25 17:31:08 ae59 Exp $ */
/*
$Log: i_sin.c,v $
 * Revision 1.2  1997/09/25  17:31:08  ae59
 * Copyright-message
 *
 * Revision 1.1  1994/01/27  13:23:04  ae59
 * Put under RCS control
 *
 */
/****************************************************************/
/*                                                              */
/*      Filename        : i_sin.c                               */
/*                                                              */
/*      Entries         : a_intv i_sin(a)                       */
/*                        a_intv a;                             */
/*                                                              */
/*      Arguments       : a = interval argument of function     */
/*                                                              */
/*      Description     : Interval sine function                */
/*                                                              */
/*      Date            : 1993-08-23                            */
/*      1993-08-23 : use B_TYPINI                        =b=    */
/*      1992-07-22 : cast to unsigned of jub, jlb        =b=    */
/*      1992-03-05 : chg unsigned to a_btyp in LhI decl. =b=    */
/*      1991-11-14 : save and restore flags                     */
/*      1991-09-06 : reduce maximum line length                 */
/*      1991-04-22 : real assignment replaced by R_ASSIGN()     */
/*      1991-02-21 : extern L...() declaration removed          */
/*      1991-01-22 : compare with r_zero replaced by r_sign()   */
/*      1990-10-18 : real operations                            */
/****************************************************************/

#ifndef ALL_IN_ONE
#ifdef AIX
#include "/u/p88c/runtime/o_defs.h"
#else
#include "o_defs.h"
#endif
#define local
extern a_real *r_mone,*r_one_;
#endif

extern a_btyp LhI;    /* entier(x/(pi/4)) mod 16 after call sin(x) */

#ifdef LINT_ARGS
local a_intv i_sin( a_intv a )
#else
local a_intv i_sin( a )
a_intv a;
#endif

{
        /* --- 2*Pi rounded downwards --- */
static a_btyp o_2pil[] = B_TYPINI( 0x401921fbL, 0x54442d18L );  /* chopped */

        a_real *twopil = (a_real *)&o_2pil[0];
/*      static unsigned long c2pi[2] = {0x401921fb, 0x54442d18};
        a_real            *twopil= (a_real*)&c2pi[0];    */

        a_btyp        rc;
        a_intv          res;
        a_intv          clb, cub;
        unsigned        jub, jlb, d;

        E_SPUSH("i_sin")

        /* --- point argument --- */
        if (i_point(a)) {
            /* printf("point argument!\n"); */
            if (r_sign(a.INF)==0) {
              res = a;
              rc = 0;
            }
            else {
              rc = i_invp(Lsin,a.INF,&res.INF,&res.SUP);
              /* printf("entier part Lhi: %d\n", LhI); */
            }
        }

        /* --- interval argument --- */
        else if (i_iv(a)) {
/*          printf("interval argument!\n");
*/
            rc = i_invp(Lsin,a.INF,&clb.INF,&clb.SUP);
/* korrigiert 6.8.90
            jlb = (LhI/2)%4;
*/
            jlb = (unsigned) (LhI/2);      /* entier(lb/(pi/2)) */

/*          printf("entier part Lhi lb mod 16: %u\n", LhI);
*/
            rc += i_invp(Lsin,a.SUP,&cub.INF,&cub.SUP);
            jub = (unsigned) (LhI/2);      /* entier(lb/(pi/2)) mod 4 */

/*          printf("entier part Lhi ub mod 16: %u\n", LhI);
*/
            if (rc==0) {
            /* --- d --- */
            if (jlb>jub) jub += 8;
            d = jub - jlb;        /* d in {0,..,7} */

            if (r_lt(r_addd(a.INF,*twopil),a.SUP))  d = 4;

            if (d>4) d=4;

            jlb%=4;
            jub%=4;

/*          printf("d = jub - jlb            : %u\n", d);
            printf("entier part Lhi lb mod 4 : %u\n", jlb);
            printf("entier part Lhi ub mod 4 : %u\n", jub);
*/
            switch (d) {
            case 0: /* -------------------------- d = 0 --- */
                /* jub = jlb */
                switch (jlb) {
                case 0:
                case 3: /* increasing */
                    R_ASSIGN(res.INF,clb.INF);
                    R_ASSIGN(res.SUP,cub.SUP);
                    break;
                case 1:
                case 2: /* decreasing */
                    R_ASSIGN(res.INF,cub.INF);
                    R_ASSIGN(res.SUP,clb.SUP);
                    break;
                    }
                break; /* d == 0 */

            case 1: /* -------------------------- d = 1 --- */
                /* jub = jlb + 1 */
                switch (jlb) {
                case 0:
                    R_ASSIGN(res.INF,
                             (r_lt(cub.INF,clb.INF))?cub.INF:clb.INF);
                    R_ASSIGN(res.SUP,*r_one_);
                    break;
                case 1: /* decreasing */
                    R_ASSIGN(res.INF,cub.INF);
                    R_ASSIGN(res.SUP,clb.SUP);
                    break;
                case 2:
                    R_ASSIGN(res.INF,*r_mone);
                    R_ASSIGN(res.SUP,
                             (r_gt(cub.SUP,clb.SUP))?cub.SUP:clb.SUP);
                    break;
                case 3: /* increasing */
                    R_ASSIGN(res.INF,clb.INF);
                    R_ASSIGN(res.SUP,cub.SUP);
                    break;
                }
                break; /* d == 1 */

            case 2: /* -------------------------- d = 2 --- */
                /* jub = jlb + 2 */
                switch (jlb) {
                case 0:
                    /* jub = 2 */
                    R_ASSIGN(res.INF,cub.INF);
                    R_ASSIGN(res.SUP,*r_one_);
                    break;
                case 1:
                    /* jub = 3 */
                    R_ASSIGN(res.INF,*r_mone);
                    R_ASSIGN(res.SUP,clb.SUP);
                    break;
                case 2:
                    /* jub = 0 */
                    R_ASSIGN(res.INF,*r_mone);
                    R_ASSIGN(res.SUP,cub.SUP);
                    break;
                case 3:
                    /* jub = 1 */
                    R_ASSIGN(res.INF,clb.INF);
                    R_ASSIGN(res.SUP,*r_one_);
                    break;
                }
                break; /* d == 2 */

            case 3: /* -------------------------- d = 3 --- */
                /* jub = jlb + 3 */
                switch (jlb) {
                case 0:
                case 2:
                    R_ASSIGN(res.INF,*r_mone);
                    R_ASSIGN(res.SUP,*r_one_);
                    break;
                case 1:
                    R_ASSIGN(res.INF,*r_mone);
                    R_ASSIGN(res.SUP,
                             (r_gt(cub.SUP,clb.SUP))?cub.SUP:clb.SUP);
                    break;
                case 3:
                    R_ASSIGN(res.INF,
                             (r_lt(cub.INF,clb.INF))?cub.INF:clb.INF);
                    R_ASSIGN(res.SUP,*r_one_);
                    break;
                }
                break; /* d == 3 */

            case 4: /* -------------------------- d = 4 --- */
                /* jub = jlb + 4 */
                R_ASSIGN(res.INF,*r_mone);
                R_ASSIGN(res.SUP,*r_one_);
                break; /* d == 4 */
            } /* switch(jdiff) */

          }
        if (r_lt(res.INF,*r_mone)) R_ASSIGN(res.INF,*r_mone);
        if (r_gt(res.SUP,*r_one_)) R_ASSIGN(res.SUP,*r_one_);
        }/* interval argument */
        else
            /* --- no correct interval --- */
            rc = 1;

        /* --- error --- */
        if (rc)
            e_trap(INV_ARG,4,E_TDBL+E_TEXT(5),&a.INF,
                             E_TDBL+E_TEXT(6),&a.SUP);

        E_SPOPP("i_sin")
        return(res);
}
