sin.cpp

   1 // Sine function of numerical and symbolic arguments
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 void
   7 eval_sin(void)
   8 {
   9         push(cadr(p1));
  10         eval();
  11         sine();
  12 }
  13
  14 void
  15 sine(void)
  16 {
  17         save();
  18         p1 = pop();
  19         if (car(p1) == symbol(ADD))
  20                 sine_of_angle_sum();
  21         else
  22                 sine_of_angle();
  23         restore();
  24 }
  25
  26 // Use angle sum formula for special angles.
  27
  28 #define A p3
  29 #define B p4
  30
  31 void
  32 sine_of_angle_sum(void)
  33 {
  34         p2 = cdr(p1);
  35         while (iscons(p2)) {
  36                 B = car(p2);
  37                 if (isnpi(B)) {
  38                         push(p1);
  39                         push(B);
  40                         subtract();
  41                         A = pop();
  42                         push(A);
  43                         sine();
  44                         push(B);
  45                         cosine();
  46                         multiply();
  47                         push(A);
  48                         cosine();
  49                         push(B);
  50                         sine();
  51                         multiply();
  52                         add();
  53                         return;
  54                 }
  55                 p2 = cdr(p2);
  56         }
  57         sine_of_angle();
  58 }
  59
  60 void
  61 sine_of_angle(void)
  62 {
  63         int n;
  64         double d;
  65
  66         if (car(p1) == symbol(ARCSIN)) {
  67                 push(cadr(p1));
  68                 return;
  69         }
  70
  71         if (isdouble(p1)) {
  72                 d = sin(p1->u.d);
  73                 if (fabs(d) < 1e-10)
  74                         d = 0.0;
  75                 push_double(d);
  76                 return;
  77         }
  78
  79         // sine function is antisymmetric, sin(-x) = -sin(x)
  80
  81         if (isnegative(p1)) {
  82                 push(p1);
  83                 negate();
  84                 sine();
  85                 negate();
  86                 return;
  87         }
  88
  89         // sin(arctan(x)) = x / sqrt(1 + x^2)
  90
  91         // see p. 173 of the CRC Handbook of Mathematical Sciences
  92
  93         if (car(p1) == symbol(ARCTAN)) {
  94                 push(cadr(p1));
  95                 push_integer(1);
  96                 push(cadr(p1));
  97                 push_integer(2);
  98                 power();
  99                 add();
 100                 push_rational(-1, 2);
 101                 power();
 102                 multiply();
 103                 return;
 104         }
 105
 106         // multiply by 180/pi
 107
 108         push(p1);
 109         push_integer(180);
 110         multiply();
 111         push_symbol(PI);
 112         divide();
 113
 114         n = pop_integer();
 115
 116         if (n < 0) {
 117                 push(symbol(SIN));
 118                 push(p1);
 119                 list(2);
 120                 return;
 121         }
 122
 123         switch (n % 360) {
 124         case 0:
 125         case 180:
 126                 push_integer(0);
 127                 break;
 128         case 30:
 129         case 150:
 130                 push_rational(1, 2);
 131                 break;
 132         case 210:
 133         case 330:
 134                 push_rational(-1, 2);
 135                 break;
 136         case 45:
 137         case 135:
 138                 push_rational(1, 2);
 139                 push_integer(2);
 140                 push_rational(1, 2);
 141                 power();
 142                 multiply();
 143                 break;
 144         case 225:
 145         case 315:
 146                 push_rational(-1, 2);
 147                 push_integer(2);
 148                 push_rational(1, 2);
 149                 power();
 150                 multiply();
 151                 break;
 152         case 60:
 153         case 120:
 154                 push_rational(1, 2);
 155                 push_integer(3);
 156                 push_rational(1, 2);
 157                 power();
 158                 multiply();
 159                 break;
 160         case 240:
 161         case 300:
 162                 push_rational(-1, 2);
 163                 push_integer(3);
 164                 push_rational(1, 2);
 165                 power();
 166                 multiply();
 167                 break;
 168         case 90:
 169                 push_integer(1);
 170                 break;
 171         case 270:
 172                 push_integer(-1);
 173                 break;
 174         default:
 175                 push(symbol(SIN));
 176                 push(p1);
 177                 list(2);
 178                 break;
 179         }
 180 }
 181
 182 #if SELFTEST
 183
 184 static char *s[] = {
 185
 186         "sin(x)",
 187         "sin(x)",
 188
 189         "sin(-x)",
 190         "-sin(x)",
 191
 192         "sin(b-a)",
 193         "-sin(a-b)",
 194
 195         // check against the floating point math library
 196
 197         "f(a,x)=1+sin(float(a/360*2*pi))-float(x)+sin(a/360*2*pi)-x",
 198         "",
 199
 200         "f(0,0)",                       // 0
 201         "1",
 202
 203         "f(90,1)",                      // 90
 204         "1",
 205
 206         "f(180,0)",                     // 180
 207         "1",
 208
 209         "f(270,-1)",                    // 270
 210         "1",
 211
 212         "f(360,0)",                     // 360
 213         "1",
 214
 215         "f(-90,-1)",                    // -90
 216         "1",
 217
 218         "f(-180,0)",                    // -180
 219         "1",
 220
 221         "f(-270,1)",                    // -270
 222         "1",
 223
 224         "f(-360,0)",                    // -360
 225         "1",
 226
 227         "f(45,sqrt(2)/2)",              // 45
 228         "1",
 229
 230         "f(135,sqrt(2)/2)",             // 135
 231         "1",
 232
 233         "f(225,-sqrt(2)/2)",            // 225
 234         "1",
 235
 236         "f(315,-sqrt(2)/2)",            // 315
 237         "1",
 238
 239         "f(-45,-sqrt(2)/2)",            // -45
 240         "1",
 241
 242         "f(-135,-sqrt(2)/2)",           // -135
 243         "1",
 244
 245         "f(-225,sqrt(2)/2)",            // -225
 246         "1",
 247
 248         "f(-315,sqrt(2)/2)",            // -315
 249         "1",
 250
 251         "f(30,1/2)",                    // 30
 252         "1",
 253
 254         "f(150,1/2)",                   // 150
 255         "1",
 256
 257         "f(210,-1/2)",                  // 210
 258         "1",
 259
 260         "f(330,-1/2)",                  // 330
 261         "1",
 262
 263         "f(-30,-1/2)",                  // -30
 264         "1",
 265
 266         "f(-150,-1/2)",                 // -150
 267         "1",
 268
 269         "f(-210,1/2)",                  // -210
 270         "1",
 271
 272         "f(-330,1/2)",                  // -330
 273         "1",
 274
 275         "f(60,sqrt(3)/2)",              // 60
 276         "1",
 277
 278         "f(120,sqrt(3)/2)",             // 120
 279         "1",
 280
 281         "f(240,-sqrt(3)/2)",            // 240
 282         "1",
 283
 284         "f(300,-sqrt(3)/2)",            // 300
 285         "1",
 286
 287         "f(-60,-sqrt(3)/2)",            // -60
 288         "1",
 289
 290         "f(-120,-sqrt(3)/2)",           // -120
 291         "1",
 292
 293         "f(-240,sqrt(3)/2)",            // -240
 294         "1",
 295
 296         "f(-300,sqrt(3)/2)",            // -300
 297         "1",
 298
 299         "f=quote(f)",
 300         "",
 301
 302         "sin(arcsin(x))",
 303         "x",
 304
 305         // check the default case
 306
 307         "sin(1/12*pi)",
 308         "sin(1/12*pi)",
 309
 310         "sin(arctan(4/3))",
 311         "4/5",
 312
 313         "sin(-arctan(4/3))",
 314         "-4/5",
 315
 316         // phase
 317
 318         "sin(x-8/2*pi)",
 319         "sin(x)",
 320
 321         "sin(x-7/2*pi)",
 322         "cos(x)",
 323
 324         "sin(x-6/2*pi)",
 325         "-sin(x)",
 326
 327         "sin(x-5/2*pi)",
 328         "-cos(x)",
 329
 330         "sin(x-4/2*pi)",
 331         "sin(x)",
 332
 333         "sin(x-3/2*pi)",
 334         "cos(x)",
 335
 336         "sin(x-2/2*pi)",
 337         "-sin(x)",
 338
 339         "sin(x-1/2*pi)",
 340         "-cos(x)",
 341
 342         "sin(x+0/2*pi)",
 343         "sin(x)",
 344
 345         "sin(x+1/2*pi)",
 346         "cos(x)",
 347
 348         "sin(x+2/2*pi)",
 349         "-sin(x)",
 350
 351         "sin(x+3/2*pi)",
 352         "-cos(x)",
 353
 354         "sin(x+4/2*pi)",
 355         "sin(x)",
 356
 357         "sin(x+5/2*pi)",
 358         "cos(x)",
 359
 360         "sin(x+6/2*pi)",
 361         "-sin(x)",
 362
 363         "sin(x+7/2*pi)",
 364         "-cos(x)",
 365
 366         "sin(x+8/2*pi)",
 367         "sin(x)",
 368 };
 369
 370 void
 371 test_sin(void)
 372 {
 373         test(__FILE__, s, sizeof s / sizeof (char *));
 374 }
 375
 376 #endif