Hex Artifact Content
Not logged in

Artifact 1aa5b814e32c4af29d22930a08d0fdd57a4bb32e:


0000: 2f 2a 20 6d 70 66 72 5f 73 69 6e 20 2d 2d 20 73  /* mpfr_sin -- s
0010: 69 6e 65 20 6f 66 20 61 20 66 6c 6f 61 74 69 6e  ine of a floatin
0020: 67 2d 70 6f 69 6e 74 20 6e 75 6d 62 65 72 0a 0a  g-point number..
0030: 43 6f 70 79 72 69 67 68 74 20 32 30 30 31 2c 20  Copyright 2001, 
0040: 32 30 30 32 20 46 72 65 65 20 53 6f 66 74 77 61  2002 Free Softwa
0050: 72 65 20 46 6f 75 6e 64 61 74 69 6f 6e 2c 20 49  re Foundation, I
0060: 6e 63 2e 0a 0a 54 68 69 73 20 66 69 6c 65 20 69  nc...This file i
0070: 73 20 70 61 72 74 20 6f 66 20 74 68 65 20 4d 50  s part of the MP
0080: 46 52 20 4c 69 62 72 61 72 79 2e 0a 0a 54 68 65  FR Library...The
0090: 20 4d 50 46 52 20 4c 69 62 72 61 72 79 20 69 73   MPFR Library is
00a0: 20 66 72 65 65 20 73 6f 66 74 77 61 72 65 3b 20   free software; 
00b0: 79 6f 75 20 63 61 6e 20 72 65 64 69 73 74 72 69  you can redistri
00c0: 62 75 74 65 20 69 74 20 61 6e 64 2f 6f 72 20 6d  bute it and/or m
00d0: 6f 64 69 66 79 0a 69 74 20 75 6e 64 65 72 20 74  odify.it under t
00e0: 68 65 20 74 65 72 6d 73 20 6f 66 20 74 68 65 20  he terms of the 
00f0: 47 4e 55 20 4c 65 73 73 65 72 20 47 65 6e 65 72  GNU Lesser Gener
0100: 61 6c 20 50 75 62 6c 69 63 20 4c 69 63 65 6e 73  al Public Licens
0110: 65 20 61 73 20 70 75 62 6c 69 73 68 65 64 20 62  e as published b
0120: 79 0a 74 68 65 20 46 72 65 65 20 53 6f 66 74 77  y.the Free Softw
0130: 61 72 65 20 46 6f 75 6e 64 61 74 69 6f 6e 3b 20  are Foundation; 
0140: 65 69 74 68 65 72 20 76 65 72 73 69 6f 6e 20 32  either version 2
0150: 2e 31 20 6f 66 20 74 68 65 20 4c 69 63 65 6e 73  .1 of the Licens
0160: 65 2c 20 6f 72 20 28 61 74 20 79 6f 75 72 0a 6f  e, or (at your.o
0170: 70 74 69 6f 6e 29 20 61 6e 79 20 6c 61 74 65 72  ption) any later
0180: 20 76 65 72 73 69 6f 6e 2e 0a 0a 54 68 65 20 4d   version...The M
0190: 50 46 52 20 4c 69 62 72 61 72 79 20 69 73 20 64  PFR Library is d
01a0: 69 73 74 72 69 62 75 74 65 64 20 69 6e 20 74 68  istributed in th
01b0: 65 20 68 6f 70 65 20 74 68 61 74 20 69 74 20 77  e hope that it w
01c0: 69 6c 6c 20 62 65 20 75 73 65 66 75 6c 2c 20 62  ill be useful, b
01d0: 75 74 0a 57 49 54 48 4f 55 54 20 41 4e 59 20 57  ut.WITHOUT ANY W
01e0: 41 52 52 41 4e 54 59 3b 20 77 69 74 68 6f 75 74  ARRANTY; without
01f0: 20 65 76 65 6e 20 74 68 65 20 69 6d 70 6c 69 65   even the implie
0200: 64 20 77 61 72 72 61 6e 74 79 20 6f 66 20 4d 45  d warranty of ME
0210: 52 43 48 41 4e 54 41 42 49 4c 49 54 59 0a 6f 72  RCHANTABILITY.or
0220: 20 46 49 54 4e 45 53 53 20 46 4f 52 20 41 20 50   FITNESS FOR A P
0230: 41 52 54 49 43 55 4c 41 52 20 50 55 52 50 4f 53  ARTICULAR PURPOS
0240: 45 2e 20 20 53 65 65 20 74 68 65 20 47 4e 55 20  E.  See the GNU 
0250: 4c 65 73 73 65 72 20 47 65 6e 65 72 61 6c 20 50  Lesser General P
0260: 75 62 6c 69 63 0a 4c 69 63 65 6e 73 65 20 66 6f  ublic.License fo
0270: 72 20 6d 6f 72 65 20 64 65 74 61 69 6c 73 2e 0a  r more details..
0280: 0a 59 6f 75 20 73 68 6f 75 6c 64 20 68 61 76 65  .You should have
0290: 20 72 65 63 65 69 76 65 64 20 61 20 63 6f 70 79   received a copy
02a0: 20 6f 66 20 74 68 65 20 47 4e 55 20 4c 65 73 73   of the GNU Less
02b0: 65 72 20 47 65 6e 65 72 61 6c 20 50 75 62 6c 69  er General Publi
02c0: 63 20 4c 69 63 65 6e 73 65 0a 61 6c 6f 6e 67 20  c License.along 
02d0: 77 69 74 68 20 74 68 65 20 4d 50 46 52 20 4c 69  with the MPFR Li
02e0: 62 72 61 72 79 3b 20 73 65 65 20 74 68 65 20 66  brary; see the f
02f0: 69 6c 65 20 43 4f 50 59 49 4e 47 2e 4c 49 42 2e  ile COPYING.LIB.
0300: 20 20 49 66 20 6e 6f 74 2c 20 77 72 69 74 65 20    If not, write 
0310: 74 6f 0a 74 68 65 20 46 72 65 65 20 53 6f 66 74  to.the Free Soft
0320: 77 61 72 65 20 46 6f 75 6e 64 61 74 69 6f 6e 2c  ware Foundation,
0330: 20 49 6e 63 2e 2c 20 35 39 20 54 65 6d 70 6c 65   Inc., 59 Temple
0340: 20 50 6c 61 63 65 20 2d 20 53 75 69 74 65 20 33   Place - Suite 3
0350: 33 30 2c 20 42 6f 73 74 6f 6e 2c 0a 4d 41 20 30  30, Boston,.MA 0
0360: 32 31 31 31 2d 31 33 30 37 2c 20 55 53 41 2e 20  2111-1307, USA. 
0370: 2a 2f 0a 0a 23 69 6e 63 6c 75 64 65 20 22 67 6d  */..#include "gm
0380: 70 2e 68 22 0a 23 69 6e 63 6c 75 64 65 20 22 67  p.h".#include "g
0390: 6d 70 2d 69 6d 70 6c 2e 68 22 0a 23 69 6e 63 6c  mp-impl.h".#incl
03a0: 75 64 65 20 22 6d 70 66 72 2e 68 22 0a 23 69 6e  ude "mpfr.h".#in
03b0: 63 6c 75 64 65 20 22 6d 70 66 72 2d 69 6d 70 6c  clude "mpfr-impl
03c0: 2e 68 22 0a 0a 69 6e 74 20 0a 6d 70 66 72 5f 73  .h"..int .mpfr_s
03d0: 69 6e 20 28 6d 70 66 72 5f 70 74 72 20 79 2c 20  in (mpfr_ptr y, 
03e0: 6d 70 66 72 5f 73 72 63 70 74 72 20 78 2c 20 6d  mpfr_srcptr x, m
03f0: 70 5f 72 6e 64 5f 74 20 72 6e 64 5f 6d 6f 64 65  p_rnd_t rnd_mode
0400: 29 20 0a 7b 0a 20 20 69 6e 74 20 70 72 65 63 79  ) .{.  int precy
0410: 2c 20 6d 2c 20 6f 6b 2c 20 65 2c 20 69 6e 65 78  , m, ok, e, inex
0420: 61 63 74 2c 20 6e 65 67 3b 0a 20 20 6d 70 66 72  act, neg;.  mpfr
0430: 5f 74 20 63 2c 20 6b 3b 0a 0a 20 20 69 66 20 28  _t c, k;..  if (
0440: 4d 50 46 52 5f 49 53 5f 4e 41 4e 28 78 29 20 7c  MPFR_IS_NAN(x) |
0450: 7c 20 4d 50 46 52 5f 49 53 5f 49 4e 46 28 78 29  | MPFR_IS_INF(x)
0460: 29 0a 20 20 20 20 7b 0a 20 20 20 20 20 20 4d 50  ).    {.      MP
0470: 46 52 5f 53 45 54 5f 4e 41 4e 28 79 29 3b 0a 20  FR_SET_NAN(y);. 
0480: 20 20 20 20 20 4d 50 46 52 5f 52 45 54 5f 4e 41       MPFR_RET_NA
0490: 4e 3b 0a 20 20 20 20 7d 0a 0a 20 20 69 66 20 28  N;.    }..  if (
04a0: 4d 50 46 52 5f 49 53 5f 5a 45 52 4f 28 78 29 29  MPFR_IS_ZERO(x))
04b0: 0a 20 20 20 20 7b 0a 20 20 20 20 20 20 4d 50 46  .    {.      MPF
04c0: 52 5f 43 4c 45 41 52 5f 46 4c 41 47 53 28 79 29  R_CLEAR_FLAGS(y)
04d0: 3b 0a 20 20 20 20 20 20 4d 50 46 52 5f 53 45 54  ;.      MPFR_SET
04e0: 5f 5a 45 52 4f 28 79 29 3b 0a 20 20 20 20 20 20  _ZERO(y);.      
04f0: 4d 50 46 52 5f 53 45 54 5f 53 41 4d 45 5f 53 49  MPFR_SET_SAME_SI
0500: 47 4e 28 79 2c 20 78 29 3b 0a 20 20 20 20 20 20  GN(y, x);.      
0510: 4d 50 46 52 5f 52 45 54 28 30 29 3b 0a 20 20 20  MPFR_RET(0);.   
0520: 20 7d 0a 0a 20 20 70 72 65 63 79 20 3d 20 4d 50   }..  precy = MP
0530: 46 52 5f 50 52 45 43 28 79 29 3b 0a 20 20 6d 20  FR_PREC(y);.  m 
0540: 3d 20 70 72 65 63 79 20 2b 20 5f 6d 70 66 72 5f  = precy + _mpfr_
0550: 63 65 69 6c 5f 6c 6f 67 32 20 28 28 64 6f 75 62  ceil_log2 ((doub
0560: 6c 65 29 20 70 72 65 63 79 29 20 2b 20 41 42 53  le) precy) + ABS
0570: 28 4d 50 46 52 5f 45 58 50 28 78 29 29 20 2b 20  (MPFR_EXP(x)) + 
0580: 31 33 3b 0a 0a 20 20 6d 70 66 72 5f 69 6e 69 74  13;..  mpfr_init
0590: 32 20 28 63 2c 20 6d 29 3b 0a 20 20 6d 70 66 72  2 (c, m);.  mpfr
05a0: 5f 69 6e 69 74 32 20 28 6b 2c 20 6d 29 3b 0a 0a  _init2 (k, m);..
05b0: 20 20 2f 2a 20 66 69 72 73 74 20 64 65 74 65 72    /* first deter
05c0: 6d 69 6e 65 20 73 69 67 6e 20 2a 2f 0a 20 20 6d  mine sign */.  m
05d0: 70 66 72 5f 63 6f 6e 73 74 5f 70 69 20 28 63 2c  pfr_const_pi (c,
05e0: 20 47 4d 50 5f 52 4e 44 4e 29 3b 0a 20 20 6d 70   GMP_RNDN);.  mp
05f0: 66 72 5f 6d 75 6c 5f 32 75 69 20 28 63 2c 20 63  fr_mul_2ui (c, c
0600: 2c 20 31 2c 20 47 4d 50 5f 52 4e 44 4e 29 3b 20  , 1, GMP_RNDN); 
0610: 2f 2a 20 32 2a 50 69 20 2a 2f 0a 20 20 6d 70 66  /* 2*Pi */.  mpf
0620: 72 5f 64 69 76 20 28 6b 2c 20 78 2c 20 63 2c 20  r_div (k, x, c, 
0630: 47 4d 50 5f 52 4e 44 4e 29 3b 20 20 20 20 20 20  GMP_RNDN);      
0640: 2f 2a 20 78 2f 28 32 2a 50 69 29 20 2a 2f 0a 20  /* x/(2*Pi) */. 
0650: 20 6d 70 66 72 5f 66 6c 6f 6f 72 20 28 6b 2c 20   mpfr_floor (k, 
0660: 6b 29 3b 20 20 20 20 20 20 20 20 20 20 20 20 20  k);             
0670: 20 20 20 20 2f 2a 20 66 6c 6f 6f 72 28 78 2f 28      /* floor(x/(
0680: 32 2a 50 69 29 29 20 2a 2f 0a 20 20 6d 70 66 72  2*Pi)) */.  mpfr
0690: 5f 6d 75 6c 20 28 63 2c 20 6b 2c 20 63 2c 20 47  _mul (c, k, c, G
06a0: 4d 50 5f 52 4e 44 4e 29 3b 0a 20 20 6d 70 66 72  MP_RNDN);.  mpfr
06b0: 5f 73 75 62 20 28 6b 2c 20 78 2c 20 63 2c 20 47  _sub (k, x, c, G
06c0: 4d 50 5f 52 4e 44 4e 29 3b 20 20 20 20 20 20 2f  MP_RNDN);      /
06d0: 2a 20 30 20 3c 3d 20 6b 20 3c 20 32 2a 50 69 20  * 0 <= k < 2*Pi 
06e0: 2a 2f 0a 20 20 6d 70 66 72 5f 63 6f 6e 73 74 5f  */.  mpfr_const_
06f0: 70 69 20 28 63 2c 20 47 4d 50 5f 52 4e 44 4e 29  pi (c, GMP_RNDN)
0700: 3b 20 2f 2a 20 63 61 63 68 65 64 20 2a 2f 0a 20  ; /* cached */. 
0710: 20 6e 65 67 20 3d 20 6d 70 66 72 5f 63 6d 70 20   neg = mpfr_cmp 
0720: 28 6b 2c 20 63 29 20 3e 20 30 3b 0a 20 20 6d 70  (k, c) > 0;.  mp
0730: 66 72 5f 63 6c 65 61 72 20 28 6b 29 3b 0a 0a 20  fr_clear (k);.. 
0740: 20 64 6f 0a 20 20 20 20 7b 0a 20 20 20 20 20 20   do.    {.      
0750: 6d 70 66 72 5f 63 6f 73 20 28 63 2c 20 78 2c 20  mpfr_cos (c, x, 
0760: 47 4d 50 5f 52 4e 44 5a 29 3b 0a 20 20 20 20 20  GMP_RNDZ);.     
0770: 20 6d 70 66 72 5f 6d 75 6c 20 28 63 2c 20 63 2c   mpfr_mul (c, c,
0780: 20 63 2c 20 47 4d 50 5f 52 4e 44 55 29 3b 0a 20   c, GMP_RNDU);. 
0790: 20 20 20 20 20 6d 70 66 72 5f 75 69 5f 73 75 62       mpfr_ui_sub
07a0: 20 28 63 2c 20 31 2c 20 63 2c 20 47 4d 50 5f 52   (c, 1, c, GMP_R
07b0: 4e 44 4e 29 3b 0a 20 20 20 20 20 20 65 20 3d 20  NDN);.      e = 
07c0: 32 20 2b 20 28 2d 4d 50 46 52 5f 45 58 50 28 63  2 + (-MPFR_EXP(c
07d0: 29 29 20 2f 20 32 3b 0a 20 20 20 20 20 20 6d 70  )) / 2;.      mp
07e0: 66 72 5f 73 71 72 74 20 28 63 2c 20 63 2c 20 47  fr_sqrt (c, c, G
07f0: 4d 50 5f 52 4e 44 4e 29 3b 0a 20 20 20 20 20 20  MP_RNDN);.      
0800: 69 66 20 28 6e 65 67 29 0a 09 6d 70 66 72 5f 6e  if (neg)..mpfr_n
0810: 65 67 20 28 63 2c 20 63 2c 20 47 4d 50 5f 52 4e  eg (c, c, GMP_RN
0820: 44 4e 29 3b 0a 0a 20 20 20 20 20 20 2f 2a 20 74  DN);..      /* t
0830: 68 65 20 61 62 73 6f 6c 75 74 65 20 65 72 72 6f  he absolute erro
0840: 72 20 6f 6e 20 63 20 69 73 20 61 74 20 6d 6f 73  r on c is at mos
0850: 74 20 32 5e 28 65 2d 6d 29 20 3d 20 32 5e 28 45  t 2^(e-m) = 2^(E
0860: 58 50 28 63 29 2d 65 72 72 29 20 2a 2f 0a 20 20  XP(c)-err) */.  
0870: 20 20 20 20 65 20 3d 20 4d 50 46 52 5f 45 58 50      e = MPFR_EXP
0880: 28 63 29 20 2b 20 6d 20 2d 20 65 3b 0a 20 20 20  (c) + m - e;.   
0890: 20 20 20 6f 6b 20 3d 20 28 65 20 3e 3d 20 30 29     ok = (e >= 0)
08a0: 20 26 26 20 6d 70 66 72 5f 63 61 6e 5f 72 6f 75   && mpfr_can_rou
08b0: 6e 64 20 28 63 2c 20 65 2c 20 47 4d 50 5f 52 4e  nd (c, e, GMP_RN
08c0: 44 4e 2c 20 72 6e 64 5f 6d 6f 64 65 2c 20 70 72  DN, rnd_mode, pr
08d0: 65 63 79 29 3b 0a 0a 20 20 20 20 20 20 69 66 20  ecy);..      if 
08e0: 28 6f 6b 20 3d 3d 20 30 29 0a 09 7b 0a 09 20 20  (ok == 0)..{..  
08f0: 6d 20 2b 3d 20 42 49 54 53 5f 50 45 52 5f 4d 50  m += BITS_PER_MP
0900: 5f 4c 49 4d 42 3b 0a 09 20 20 6d 70 66 72 5f 73  _LIMB;..  mpfr_s
0910: 65 74 5f 70 72 65 63 20 28 63 2c 20 6d 29 3b 0a  et_prec (c, m);.
0920: 09 7d 0a 20 20 20 20 7d 0a 20 20 77 68 69 6c 65  .}.    }.  while
0930: 20 28 21 6f 6b 29 3b 0a 0a 20 20 69 6e 65 78 61   (!ok);..  inexa
0940: 63 74 20 3d 20 6d 70 66 72 5f 73 65 74 20 28 79  ct = mpfr_set (y
0950: 2c 20 63 2c 20 72 6e 64 5f 6d 6f 64 65 29 3b 0a  , c, rnd_mode);.
0960: 0a 20 20 6d 70 66 72 5f 63 6c 65 61 72 20 28 63  .  mpfr_clear (c
0970: 29 3b 0a 0a 20 20 72 65 74 75 72 6e 20 69 6e 65  );..  return ine
0980: 78 61 63 74 3b 20 2f 2a 20 69 6e 65 78 61 63 74  xact; /* inexact
0990: 20 2a 2f 0a 7d 0a                                 */.}.