nl-math.c

   1 /* nl-math.c
   2
   3     Copyright (C) 2008 Lutz Mueller
   4
   5     This program is free software: you can redistribute it and/or modify
   6     it under the terms of the GNU General Public License as published by
   7     the Free Software Foundation, either version 3 of the License, or
   8     (at your option) any later version.
   9
  10     This program is distributed in the hope that it will be useful,
  11     but WITHOUT ANY WARRANTY; without even the implied warranty of
  12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13     GNU General Public License for more details.
  14
  15     You should have received a copy of the GNU General Public License
  16     along with this program.  If not, see <http://www.gnu.org/licenses/>.
  17
  18 */
  19
  20 #include "newlisp.h"
  21 #include "protos.h"
  22
  23 #ifdef WINCE
  24 #define _exception exception
  25 #endif
  26
  27 #define OP_ADD 1
  28 #define OP_SUBTRACT 2
  29 #define OP_MULTIPLY 3
  30 #define OP_DIVIDE 4
  31 #define OP_BIT_OR 5
  32 #define OP_BIT_AND 6
  33 #define OP_BIT_XOR 7
  34 #define OP_SHIFTL 8
  35 #define OP_SHIFTR 9
  36 #define OP_POW 10
  37 #define OP_MODULO 11
  38 #define OP_SIN 12
  39 #define OP_COS 13
  40 #define OP_TAN 14
  41 #define OP_ASIN 15
  42 #define OP_ACOS 16
  43 #define OP_ATAN 17
  44 #define OP_SINH 18
  45 #define OP_COSH 19
  46 #define OP_TANH 20
  47 #define OP_ASINH 21
  48 #define OP_ACOSH 22
  49 #define OP_ATANH 23
  50 #define OP_SQRT 24
  51 #define OP_LOG 25
  52 #define OP_EXP 26
  53 #define OP_MIN 27
  54 #define OP_MAX 28
  55 #define OP_ABS 29
  56 #define OP_CEIL 30
  57 #define OP_FLOOR 31
  58 #define OP_NAN 32
  59 #define OP_ERRORFUNC 33
  60 #define OP_SIGNUM 34
  61 #define OP_ISNAN 35
  62
  63 #ifdef WIN_32
  64 int _matherr(struct _exception *e) {return 1;}
  65 #endif
  66
  67 CELL * incDec(CELL * params, int type)
  68 {
  69 SYMBOL * sPtr;
  70 INT64 lValue;
  71 double fValue;
  72 double adjust;
  73 CELL * cell;
  74
  75 if(params->next != nilCell)
  76         {
  77         getFloat(params->next, &adjust);
  78         adjust = adjust * type;
  79         type = 0;
  80         }
  81
  82 cell = evaluateExpression(params);
  83 if(cell->type != CELL_SYMBOL)
  84         {
  85         if(cell->type == CELL_DYN_SYMBOL)
  86                 sPtr = getDynamicSymbol(cell);
  87         else
  88                 return(errorProcExt(ERR_SYMBOL_EXPECTED, params));
  89         }
  90 else
  91         sPtr = (SYMBOL *)cell->contents;
  92
  93 if(isProtected(sPtr->flags))
  94         return(errorProcExt2(ERR_SYMBOL_PROTECTED, stuffSymbol(sPtr)));
  95
  96
  97 if(type == 0)
  98         {
  99         getFloat(cell, &fValue);
 100         deleteList((CELL *)sPtr->contents);
 101         fValue = fValue + adjust;
 102         cell = getCell(CELL_FLOAT);
 103 #ifndef NEWLISP64
 104         *(double *)&cell->aux = fValue;
 105 #else
 106         *(double *)&cell->contents = fValue;
 107 #endif
 108         sPtr->contents = (UINT)cell;
 109         }
 110 else
 111         {
 112         getInteger64(cell, &lValue);
 113         deleteList((CELL *)sPtr->contents);
 114         sPtr->contents = (UINT)stuffInteger64(lValue + type);
 115         }
 116
 117 pushResultFlag = FALSE;
 118 return((CELL *)sPtr->contents);
 119 }
 120
 121
 122 CELL * p_increment(CELL * params) { return(incDec(params, 1)); }
 123 CELL * p_decrement(CELL * params) { return(incDec(params, -1)); }
 124
 125
 126 CELL * arithmetikOp(CELL * params, int op)
 127 {
 128 INT64 number;
 129 INT64 result;
 130
 131 params = getInteger64(params, &number);
 132 result = number;
 133
 134 if(params == nilCell)
 135         {
 136         switch(op)
 137                 {
 138                 case OP_SUBTRACT:
 139                         result = - result;
 140                         break;
 141                 case OP_SHIFTL:
 142                         result <<= 1;
 143                         break;
 144                 case OP_SHIFTR:
 145                         result >>= 1;
 146                 default:
 147                         break;
 148                 }
 149         }
 150
 151 while(params != nilCell)
 152         {
 153         params = getInteger64(params, &number);
 154         switch(op)
 155                 {
 156                 case OP_ADD:            result += number; break;
 157                 case OP_SUBTRACT:       result -= number; break;
 158                 case OP_MULTIPLY:       result *= number; break;
 159                 case OP_DIVIDE:
 160                         if(number == 0) return(errorProc(ERR_MATH));
 161                         result /= number; break;
 162                 case OP_BIT_OR:         result |= number; break;
 163                 case OP_BIT_AND:        result &= number; break;
 164                 case OP_BIT_XOR:        result ^= number; break;
 165                 case OP_SHIFTL:         result <<= number; break;
 166                 case OP_SHIFTR:         result >>= number; break;
 167                 case OP_MODULO:
 168                         if(number == 0) return(errorProc(ERR_MATH));
 169                         result %= number; break;
 170                 default:
 171                         break;
 172                 }
 173         }
 174
 175 #ifndef NEWLISP64
 176 params = getCell(CELL_INT64);
 177 *(INT64 *)&params->aux = result;
 178 return(params);
 179 #else
 180 return(stuffInteger(result));
 181 #endif
 182 }
 183
 184
 185 CELL * p_add(CELL * params) { return(arithmetikOp(params, OP_ADD)); }
 186 CELL * p_subtract(CELL * params) { return(arithmetikOp(params, OP_SUBTRACT)); }
 187 CELL * p_multiply(CELL * params) { return(arithmetikOp(params, OP_MULTIPLY)); }
 188 CELL * p_divide(CELL * params) { return(arithmetikOp(params, OP_DIVIDE)); }
 189
 190 CELL * p_bitOr(CELL * params) { return(arithmetikOp(params, OP_BIT_OR)); }
 191 CELL * p_bitAnd(CELL * params) { return(arithmetikOp(params, OP_BIT_AND)); }
 192 CELL * p_bitXor(CELL * params) { return(arithmetikOp(params, OP_BIT_XOR)); }
 193 CELL * p_shiftLeft(CELL * params) { return(arithmetikOp(params, OP_SHIFTL)); }
 194 CELL * p_shiftRight(CELL * params) { return(arithmetikOp(params, OP_SHIFTR)); }
 195 CELL * p_modulo(CELL * params) { return(arithmetikOp(params, OP_MODULO)); }
 196
 197
 198 CELL * p_bitNot(CELL * params)
 199 {
 200 INT64 number;
 201 getInteger64(params, &number);
 202 return(stuffInteger64(~number));
 203 }
 204
 205
 206 /* ----------------------- float arithmetik ----------------------------- */
 207
 208 CELL * getFloat(CELL * params, double *);
 209 CELL * floatOp(CELL * params, int op);
 210 int compareFloats(CELL * left, CELL * right);
 211 int compareInts(CELL * left, CELL * right);
 212
 213
 214 CELL * p_addFloat(CELL * params) { return(floatOp(params, OP_ADD)); }
 215 CELL * p_subFloat(CELL * params) { return(floatOp(params, OP_SUBTRACT)); }
 216 CELL * p_mulFloat(CELL * params) { return(floatOp(params, OP_MULTIPLY)); }
 217 CELL * p_divFloat(CELL * params) { return(floatOp(params, OP_DIVIDE)); }
 218 CELL * p_minFloat(CELL * params) { return(floatOp(params, OP_MIN)); }
 219 CELL * p_maxFloat(CELL * params) { return(floatOp(params, OP_MAX)); }
 220 CELL * p_powFloat(CELL * params) { return(floatOp(params, OP_POW)); }
 221 CELL * p_modFloat(CELL * params) { return(floatOp(params, OP_MODULO)); }
 222
 223 CELL * floatOp(CELL * params, int op)
 224 {
 225 double number;
 226 double result;
 227
 228 params = getFloat(params, &result);
 229 if(params == nilCell)
 230         {
 231         if(op == OP_SUBTRACT)
 232                 result = - result;
 233         else if(op == OP_DIVIDE)
 234                 result = 1.0 / result;
 235         else if(op == OP_POW)
 236                 result = result * result;
 237         }
 238 while(params != nilCell)
 239         {
 240         params = getFloat(params, &number);
 241 #ifdef WIN_32
 242         if(isnan(number))
 243                 {
 244                 result = number;
 245                 break;
 246                 }
 247 #endif
 248         switch(op)
 249                 {
 250                 case OP_ADD:            result += number; break;
 251                 case OP_SUBTRACT:       result -= number; break;
 252                 case OP_MULTIPLY:       result *= number; break;
 253                 case OP_DIVIDE:
 254                         if(number == 0.0)
 255                                 return(errorProc(ERR_MATH));
 256                         result /= number;
 257                         break;
 258                 case OP_MIN: if(number < result) result = number; break;
 259                 case OP_MAX: if(number > result) result = number; break;
 260                 case OP_POW: result = pow(result, number); break;
 261                 case OP_MODULO:
 262                         if(number == 0.0)
 263                                 return(errorProc(ERR_MATH));
 264                         result = fmod(result, number);
 265                 default: break;
 266                 }
 267         }
 268
 269 params = getCell(CELL_FLOAT);
 270 #ifndef NEWLISP64
 271 memcpy((void *)&params->aux, (void *)&result, sizeof(double));
 272 #else
 273 *(double *)&params->contents = result;
 274 #endif
 275 return(params);
 276 }
 277
 278
 279 int compareFloats(CELL * left, CELL * right)
 280 {
 281 double leftFloat, rightFloat;
 282
 283 leftFloat = getDirectFloat(left);
 284 rightFloat = getDirectFloat(right);
 285
 286 if(isnan(leftFloat) && isnan(rightFloat)) return(0);
 287 if(isnan(leftFloat) || isnan(rightFloat)) return(9);
 288
 289 if(leftFloat < rightFloat) return(-1);
 290 if(leftFloat > rightFloat) return(1);
 291
 292 return(0);
 293 }
 294
 295 int compareInts(CELL * left, CELL * right)
 296 {
 297 INT64 leftnum;
 298 INT64 rightnum;
 299
 300 #ifndef NEWLISP64
 301 if(left->type == CELL_LONG)
 302         leftnum = (int)left->contents;
 303 else
 304         leftnum = *(INT64 *)&left->aux;
 305
 306 if(right->type == CELL_LONG)
 307         rightnum = (int)right->contents;
 308 else
 309         rightnum = *(INT64 *)&right->aux;
 310 #else
 311 leftnum = (UINT)left->contents;
 312 rightnum = (UINT)right->contents;
 313 #endif
 314
 315 if(leftnum < rightnum) return(-1);
 316 if(leftnum > rightnum) return(1);
 317
 318 return(0);
 319 }
 320
 321
 322 double getDirectFloat(CELL * param)
 323 {
 324 double floatNum = 0.0;
 325
 326 #ifndef NEWLISP64
 327 if(param->type == CELL_FLOAT)
 328         return(*(double *)&param->aux);
 329 else if(param->type == CELL_LONG)
 330         floatNum = (long)param->contents;
 331 else if(param->type == CELL_INT64)
 332         floatNum = *(INT64 *)&param->aux;
 333 #else
 334 if(param->type == CELL_FLOAT)
 335         return(*(double *)&param->contents);
 336 else if(param->type == CELL_LONG)
 337         floatNum = (long)param->contents;
 338 #endif
 339
 340 return(floatNum);
 341 }
 342
 343 /* ----------------------------- float functions ----------------------- */
 344
 345
 346 CELL * functionFloat(CELL *, int);
 347
 348 CELL * p_sin(CELL * params) { return(functionFloat(params, OP_SIN)); }
 349 CELL * p_cos(CELL * params) { return(functionFloat(params, OP_COS)); }
 350 CELL * p_tan(CELL * params) { return(functionFloat(params, OP_TAN)); }
 351 CELL * p_asin(CELL * params) { return(functionFloat(params, OP_ASIN)); }
 352 CELL * p_acos(CELL * params) { return(functionFloat(params, OP_ACOS)); }
 353 CELL * p_atan(CELL * params) { return(functionFloat(params, OP_ATAN)); }
 354 CELL * p_sinh(CELL * params) { return(functionFloat(params, OP_SINH)); }
 355 CELL * p_cosh(CELL * params) { return(functionFloat(params, OP_COSH)); }
 356 CELL * p_tanh(CELL * params) { return(functionFloat(params, OP_TANH)); }
 357 CELL * p_asinh(CELL * params) { return(functionFloat(params, OP_ASINH)); }
 358 CELL * p_acosh(CELL * params) { return(functionFloat(params, OP_ACOSH)); }
 359 CELL * p_atanh(CELL * params) { return(functionFloat(params, OP_ATANH)); }
 360 CELL * p_sqrt(CELL * params) { return(functionFloat(params, OP_SQRT)); }
 361 CELL * p_log(CELL * params) { return(functionFloat(params, OP_LOG)); }
 362 CELL * p_exp(CELL * params) { return(functionFloat(params, OP_EXP)); }
 363 CELL * p_abs(CELL * params) { return(functionFloat(params, OP_ABS)); }
 364 CELL * p_ceil(CELL * params) { return(functionFloat(params, OP_CEIL)); }
 365 CELL * p_floor(CELL * params) { return(functionFloat(params, OP_FLOOR)); }
 366 CELL * p_erf(CELL * params) { return(functionFloat(params, OP_ERRORFUNC)); }
 367 CELL * p_sgn(CELL * params) { return(functionFloat(params, OP_SIGNUM)); }
 368 CELL * p_isnan(CELL * params) { return(functionFloat(params, OP_ISNAN)); }
 369
 370 CELL * functionFloat(CELL * params, int op)
 371 {
 372 double floatN;
 373 double base;
 374 CELL * cell;
 375
 376 getFloat(params, &floatN);
 377 #ifdef WIN_32
 378 if(isnan(floatN))
 379         return( stuffFloat(&floatN) );
 380 #endif
 381
 382 switch(op)
 383   {
 384   case OP_SIN: floatN = sin(floatN); break;
 385   case OP_COS: floatN = cos(floatN); break;
 386   case OP_TAN: floatN = tan(floatN); break;
 387
 388   case OP_ASIN: floatN = asin(floatN); break;
 389   case OP_ACOS: floatN = acos(floatN); break;
 390   case OP_ATAN: floatN = atan(floatN); break;
 391
 392   case OP_SINH: floatN = sinh(floatN); break;
 393   case OP_COSH: floatN = cosh(floatN); break;
 394   case OP_TANH: floatN = tanh(floatN); break;
 395
 396   case OP_ASINH: floatN = asinh(floatN); break;
 397   case OP_ACOSH: floatN = acosh(floatN); break;
 398   case OP_ATANH: floatN = atanh(floatN); break;
 399
 400   case OP_SQRT: floatN = sqrt(floatN); break;
 401   case OP_LOG:
 402     if(params->next != nilCell)
 403       {
 404       getFloat(params->next, &base);
 405       floatN = log(floatN) / log(base);
 406       }
 407     else
 408       floatN = log(floatN);
 409     break;
 410   case OP_EXP: floatN = exp(floatN); break;
 411   case OP_ABS: floatN = (floatN < 0.0) ? -floatN : floatN; break;
 412   case OP_CEIL: floatN = ceil(floatN); break;
 413   case OP_FLOOR: floatN = floor(floatN); break;
 414   case OP_ERRORFUNC: floatN = erf(floatN); break;
 415   case OP_SIGNUM:
 416     if(params->next == nilCell)
 417         {
 418         if(floatN > 0.0) floatN = 1.0;
 419         else if(floatN < 0.0) floatN = -1.0;
 420         else floatN = 0.0;
 421         break;
 422         }
 423     else
 424         {
 425         if(floatN < 0.0) cell = params->next;
 426         else
 427             {
 428             params = params->next;
 429             if(floatN == 0.0) cell = params->next;
 430             else {
 431                  params = params->next;
 432                  cell = params->next;
 433                  }
 434             }
 435         }
 436
 437         return(copyCell(evaluateExpression(cell)));
 438     case OP_ISNAN:
 439                 return (isnan(floatN) ? trueCell : nilCell);
 440         default: break;
 441 }
 442
 443 cell = getCell(CELL_FLOAT);
 444 #ifndef NEWLISP64
 445 *(double *)&cell->aux = floatN;
 446 #else
 447 *(double *)&cell->contents = floatN;
 448 #endif
 449 return(cell);
 450 }
 451
 452
 453 CELL * p_atan2(CELL * params)
 454 {
 455 double floatX, floatY;
 456 CELL * cell;
 457
 458 params = getFloat(params, &floatX);
 459 getFloat(params, &floatY);
 460
 461 cell = getCell(CELL_FLOAT);
 462 #ifndef NEWLISP64
 463 *(double *)&cell->aux = atan2(floatX, floatY);
 464 #else
 465 *(double *)&cell->contents = atan2(floatX, floatY);
 466 #endif
 467 return(cell);
 468 }
 469
 470
 471 CELL * p_round(CELL * params)
 472 {
 473 double fNum;
 474 double precision = 1.0;
 475 long digits = 0;
 476 char * fmt;
 477 char * result;
 478
 479
 480 params = getFloat(params, &fNum);
 481 if(params != nilCell)
 482         getInteger(params, (UINT*)&digits);
 483
 484 if(digits > 0)
 485         {
 486         precision = pow(10.0, (double)(digits > 20 ? 20 : digits));
 487         fNum = precision * floor(fNum/precision + 0.5);
 488         }
 489 else
 490         {
 491         fmt = alloca(16);
 492         result = alloca(32);
 493         snprintf(fmt, 15, "%%.%df", (int)((digits < -16) ? 16 : -digits));
 494         snprintf(result, 31, fmt, fNum);
 495         fNum = atof(result);
 496         }
 497
 498 return(stuffFloat(&fNum));
 499 }
 500
 501
 502 CELL * p_rand(CELL * params)
 503 {
 504 long range;
 505 long n;
 506 CELL * dist, * cell;
 507 long rnum;
 508 double scale;
 509
 510 params = getInteger(params, (UINT *)&range);
 511 scale = range;
 512 if(range == 0)
 513         {
 514         srandom((unsigned)time(NULL));
 515         return(trueCell);
 516         }
 517
 518 rnum = ((scale * random())/(MY_RAND_MAX - 1));
 519 if(params->type == CELL_NIL)
 520         return(stuffInteger((UINT)rnum));
 521
 522 getInteger(params, (UINT *)&n);
 523
 524 dist = getCell(CELL_EXPRESSION);
 525
 526 cell = stuffInteger((UINT)rnum);
 527 dist->contents = (UINT)cell;
 528
 529 --n;
 530 while(n-- > 0)
 531         {
 532         rnum = ((scale * random())/(MY_RAND_MAX - 1));
 533         cell->next = stuffInteger((UINT)rnum);
 534         cell = cell->next;
 535         }
 536
 537 return(dist);
 538 }
 539
 540
 541 CELL * p_amb(CELL * params)
 542 {
 543 long len;
 544
 545 if((len = listlen(params)) == 0) return(nilCell);
 546
 547 len = random() % len;
 548
 549 while(len--) params = params->next;
 550
 551 return(copyCell(evaluateExpression(params)));
 552 }
 553
 554
 555 CELL *  p_seed(CELL * params)
 556 {
 557 INT64 seedNum;
 558
 559 getInteger64(params, &seedNum);
 560
 561 srandom((UINT)(seedNum & 0xFFFFFFFF));
 562
 563 return(stuffInteger(seedNum));
 564 }
 565
 566
 567 /* -----------------------  compare ops --------------------------- */
 568
 569 #define OP_LESS 1
 570 #define OP_GREATER 2
 571 #define OP_LESS_EQUAL 3
 572 #define OP_GREATER_EQUAL 4
 573 #define OP_EQUAL 5
 574 #define OP_NOTEQUAL 6
 575
 576 CELL * p_less(CELL * params) { return(compareOp(params, OP_LESS)); }
 577 CELL * p_greater(CELL * params) { return(compareOp(params, OP_GREATER)); }
 578 CELL * p_lessEqual(CELL * params) { return(compareOp(params, OP_LESS_EQUAL)); }
 579 CELL * p_greaterEqual(CELL * params) { return(compareOp(params, OP_GREATER_EQUAL)); }
 580 CELL * p_equal(CELL * params) { return(compareOp(params, OP_EQUAL)); }
 581 CELL * p_notEqual(CELL * params) { return(compareOp(params, OP_NOTEQUAL)); }
 582
 583 CELL * compareOp(CELL * params, int op)
 584 {
 585 CELL * left;
 586 CELL * right;
 587 long cnt = 0;
 588 int comp;
 589
 590 left = evaluateExpression(params);
 591
 592 while(TRUE)
 593         {
 594         if((params = params->next) == nilCell && cnt == 0)
 595                 {
 596                 if(isNumber(left->type))
 597                         right = stuffInteger64(0);
 598                 else if(left->type == CELL_STRING)
 599                         right = stuffString("");
 600                 else if(isList(left->type))
 601                         right = getCell(CELL_EXPRESSION);
 602                 else break;
 603                 pushResult(right);
 604                 }
 605         else
 606                 right = evaluateExpression(params);
 607         ++cnt;
 608         if((comp = compareCells(left, right)) == 9) return(nilCell);
 609         switch(op)
 610                 {
 611                 case OP_LESS:
 612                         if(comp >= 0) return(nilCell);
 613                         break;
 614                 case OP_GREATER:
 615                         if(comp <= 0) return(nilCell);
 616                         break;
 617                 case OP_LESS_EQUAL:
 618                         if(comp > 0) return(nilCell);
 619                         break;
 620                 case OP_GREATER_EQUAL:
 621                         if(comp < 0) return(nilCell);
 622                         break;
 623                 case OP_EQUAL:
 624                         if(comp != 0) return(nilCell);
 625                         break;
 626                 case OP_NOTEQUAL:
 627                         if(comp == 0) return(nilCell);
 628                 default:
 629                         break;
 630                 }
 631
 632         if(params->next == nilCell) break;
 633         left = right;
 634         }
 635
 636 return(trueCell);
 637 }
 638
 639 int compareSymbols(CELL * left, CELL * right)
 640 {
 641 SYMBOL * leftS;
 642 SYMBOL * rightS;
 643 int comp;
 644
 645 if(left->type == CELL_SYMBOL)
 646         leftS = (SYMBOL *)left->contents;
 647 else
 648         leftS = getDynamicSymbol(left);
 649
 650 if(right->type == CELL_SYMBOL)
 651         rightS = (SYMBOL *)right->contents;
 652 else
 653         rightS = getDynamicSymbol(right);
 654
 655 if(leftS->context != rightS->context)
 656         {
 657         leftS = leftS->context;
 658         rightS = rightS->context;
 659         }
 660
 661 if((comp = strcmp((char *)leftS->name, (char *)rightS->name)) == 0) return(0);
 662 return (comp > 0 ? 1 : -1);
 663 }
 664
 665
 666 /* returns equal: 0 less: -1 greater: 1 or 9 if NaN or Inf equal comparison */
 667 int compareCells(CELL * left, CELL * right)
 668 {
 669 int comp;
 670
 671 if(left->type != right->type)
 672         {
 673         if(left->type == CELL_FLOAT && ((right->type & COMPARE_TYPE_MASK) == CELL_INT))
 674                 return(compareFloats(left, right));
 675         if(((left->type & COMPARE_TYPE_MASK) == CELL_INT) && right->type == CELL_FLOAT)
 676                 return(compareFloats(left, right));
 677
 678         if((left->type & COMPARE_TYPE_MASK) == CELL_INT && (right->type & COMPARE_TYPE_MASK) == CELL_INT)
 679                 return(compareInts(left, right));
 680
 681         if(isSymbol(left->type) && isSymbol(right->type))
 682                 return(compareSymbols(left, right));
 683
 684         if(isNil(left))
 685                 {
 686                 if(isNil(right)) return(0);
 687                 else return(-1);
 688                 }
 689
 690         if(isTrue(left))
 691                 {
 692                 if(isTrue(right)) return(0);
 693                 if(isNil(right)) return(1);
 694                 return(-1);
 695                 }
 696
 697         comp = (left->type & COMPARE_TYPE_MASK) - (right->type & COMPARE_TYPE_MASK);
 698         if(comp == 0) return(0);
 699
 700         return( comp > 0 ? 1 : -1);
 701         }
 702
 703 switch(left->type)
 704         {
 705         case CELL_STRING:
 706                 comp = left->aux - right->aux;
 707                 if(comp == 0)
 708                         {
 709                         if((comp = memcmp((char *)left->contents,
 710                                         (char *)right->contents,
 711                                         left->aux)) == 0) return(0);
 712                         }
 713                 else if(comp > 0)
 714                         {
 715                         if((comp = memcmp((char *)left->contents,
 716                                           (char *)right->contents,
 717                                           right->aux - 1)) == 0)
 718                                 return(1);
 719                         }
 720                 else if(comp < 0)
 721                         {
 722                         if((comp = memcmp((char *)left->contents,
 723                                           (char *)right->contents,
 724                                           left->aux - 1)) == 0)
 725                                 return(-1);
 726                         }
 727
 728                 return (comp > 0 ? 1 : -1);
 729
 730         case CELL_SYMBOL:
 731         case CELL_DYN_SYMBOL:
 732                 return(compareSymbols(left, right));
 733
 734         case CELL_QUOTE:
 735         case CELL_EXPRESSION:
 736         case CELL_LAMBDA:
 737         case CELL_MACRO:
 738                 return(compareLists((CELL*)left->contents, (CELL*)right->contents));
 739         case CELL_ARRAY:
 740                 return(compareArrays((CELL*)left, (CELL*)right));
 741         case CELL_FLOAT:
 742                 return(compareFloats(left, right));
 743 #ifndef NEWLISP64
 744         case CELL_INT64:
 745                 if(*(INT64 *)&left->aux > *(INT64 *)&right->aux) return(1);
 746                 if(*(INT64 *)&left->aux < *(INT64 *)&right->aux) return(-1);
 747                 break;
 748 #endif
 749         case CELL_LONG:
 750         default:
 751                 if((long)left->contents > (long)right->contents) return(1);
 752                 if((long)left->contents < (long)right->contents) return(-1);
 753                 break;
 754         }
 755
 756 return(0);
 757 }
 758
 759 int compareLists(CELL * left, CELL * right)
 760 {
 761 int result;
 762
 763 while(left != nilCell)
 764         {
 765         if( (result = compareCells(left, right)) != 0)
 766                 return(result);
 767         left = left->next;
 768         right = right->next;
 769         }
 770 if(right == nilCell) return(0);
 771 return(-1);
 772 }
 773
 774
 775 /* ---------------------------- encryption -----------------------------
 776
 777   XOR one-pad enryption
 778
 779 */
 780
 781
 782 void encryptPad(char *encrypted, char *data, char * key, size_t dataLen, size_t keyLen)
 783 {
 784 size_t i;
 785 for(i = 0; i < dataLen; i++)
 786         *(encrypted + i) = *(data + i) ^ *(key + i % keyLen);
 787 }
 788
 789
 790 CELL * p_encrypt(CELL * params)
 791 {
 792 char * data;
 793 char * key;
 794 char * dataEncrypted;
 795 CELL * encrypted;
 796 size_t dataLen, keyLen;
 797
 798 getStringSize(params, &data, &dataLen, TRUE);
 799 getStringSize(params->next, &key, &keyLen, TRUE);
 800
 801 if(!keyLen) return(errorProcExt(ERR_STRING_EXPECTED, params->next));
 802
 803 dataEncrypted = (char *)allocMemory(dataLen + 1);
 804 *(dataEncrypted + dataLen) = 0;
 805
 806 encryptPad(dataEncrypted, data, key, dataLen, keyLen);
 807
 808 encrypted = getCell(CELL_STRING);
 809 encrypted->contents = (UINT)dataEncrypted;
 810 encrypted->aux = dataLen + 1;
 811
 812 return(encrypted);
 813 }
 814
 815
 816
 817 /* =============================  Fast Fourier Transform ======================== */
 818
 819 CELL * fft(CELL * params, int isign);
 820 CELL * makeListFromFloats(double num1, double num2);
 821 double getFloatFromCell(CELL *);
 822 void fastFourierTransform(double data[], unsigned int  nn, int isign);
 823
 824 CELL * p_fft(CELL * params)
 825 {
 826 return fft(params, 1);
 827 }
 828
 829
 830 CELL * p_ifft(CELL * params)
 831 {
 832 return fft(params, -1);
 833 }
 834
 835 CELL * fft(CELL * params, int isign)
 836 {
 837 CELL * listData;
 838 CELL * list;
 839 CELL * cell;
 840 double * data;
 841 unsigned int i, n, N;
 842
 843 getListHead(params, &listData);
 844
 845 n = listlen(listData);
 846
 847 N = 1;
 848 while(N < n) N <<= 1;
 849
 850 data = (double *)allocMemory(N * 2 * sizeof(double));
 851 list = listData;
 852 for(i = 0; i < n*2; i+=2)
 853         {
 854         if(isNumber(list->type))
 855                 {
 856                 data[i] = getFloatFromCell(list);
 857                 data[i+1] = (double)0.0;
 858                 list = list->next;
 859                 continue;
 860                 }
 861
 862         if(list->type != CELL_EXPRESSION)
 863                 {
 864                 freeMemory(data);
 865                 return(errorProcExt(ERR_LIST_OR_NUMBER_EXPECTED, list));
 866                 }
 867
 868         cell = (CELL *)list->contents;
 869         data[i] = getFloatFromCell(cell);
 870         data[i+1] = getFloatFromCell(cell->next);
 871         list= list->next;
 872         }
 873
 874
 875 for(i = n * 2; i < N * 2; i++)
 876         data[i] = 0.0;
 877
 878 fastFourierTransform(data, N, isign);
 879
 880 list = listData = getCell(CELL_EXPRESSION);
 881 if(isign == -1)
 882         for(i = 0; i < n * 2; i++) data[i] = data[i]/N;
 883
 884 list->contents = (UINT)makeListFromFloats(data[0], data[1]);
 885 list = (CELL*)list->contents;
 886 for(i = 2; i < N * 2; i += 2)
 887         {
 888         list->next = makeListFromFloats(data[i], data[i+1]);
 889         list = list->next;
 890         }
 891
 892 freeMemory(data);
 893
 894 return(listData);
 895 }
 896
 897
 898 CELL * makeListFromFloats(double num1, double num2)
 899         {
 900         CELL * cell;
 901         CELL * list;
 902
 903         list = getCell(CELL_EXPRESSION);
 904         list->contents = (UINT)stuffFloat(&num1);
 905         cell = (CELL *)list->contents;
 906         cell->next = stuffFloat(&num2);
 907         return(list);
 908         }
 909
 910
 911 double getFloatFromCell(CELL * cell)
 912         {
 913         double number;
 914 #ifndef NEWLISP64
 915         if(cell->type == CELL_FLOAT)
 916                 memcpy((void *)&number, (void *)&cell->aux, sizeof(double));
 917         else if(cell->type == CELL_LONG)
 918                 number = (int)cell->contents;
 919         else number = (double)*(INT64 *)&cell->aux;
 920 #else
 921         if(cell->type == CELL_FLOAT)
 922                 number = *(double *)&cell->contents;
 923         else
 924                 number = (long)cell->contents;
 925 #endif
 926         return(number);
 927         }
 928
 929
 930 /* Fast Fourier Transform
 931 // algorithm modified for zero-offset data[] and double precision from
 932 //
 933 // Numerical Recipes in 'C', 2nd Edition
 934 // W.H. Press, S.A. Teukolsky
 935 // W.T. Vettering, B.P. Flannery
 936 // Cambridge University Press, 1992
 937 */
 938
 939 #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
 940
 941 void fastFourierTransform(double data[], unsigned int nn, int isign)
 942 {
 943 unsigned int n, mmax, m, j, istep, i;
 944 double wtemp, wr, wpr, wpi, wi, theta;
 945 double tempr, tempi;
 946
 947 n = nn << 1;
 948 j = 1;
 949
 950 for(i = 1; i < n; i += 2)
 951         {
 952         if(j > i)
 953                 {
 954                 SWAP(data[j-1], data[i-1]);
 955                 SWAP(data[j], data[i]);
 956                 }
 957         m = n >> 1;
 958         while(m >= 2 && j > m)
 959                 {
 960                 j -= m;
 961                 m >>= 1;
 962                 }
 963         j += m;
 964         }
 965
 966
 967 mmax = 2;
 968 while(n > mmax)
 969         {
 970         istep = mmax << 1;
 971         theta = isign * (6.28318530717959/mmax);
 972         wtemp = sin(0.5 * theta);
 973         wpr = -2.0 * wtemp * wtemp;
 974         wpi = sin(theta);
 975         wr = 1.0;
 976         wi = 0.0;
 977         for(m = 1; m < mmax; m += 2)
 978                 {
 979                 for(i = m; i <= n; i += istep)
 980                         {
 981                         j = i + mmax;
 982                         tempr = wr * data[j-1] - wi * data[j];
 983                         tempi = wr * data[j] + wi * data[j-1];
 984                         data[j-1] = data[i-1] - tempr;
 985                         data[j] = data[i] - tempi;
 986                         data[i-1] += tempr;
 987                         data[i] += tempi;
 988                         }
 989                 wtemp = wr;
 990                 wr = wtemp * wpr - wi * wpi + wr;
 991                 wi = wi * wpr + wtemp * wpi + wi;
 992                 }
 993         mmax = istep;
 994         }
 995 }
 996
 997
 998 /* --------------------------- probability distributions ----------------- */
 999
1000 #define DIST_RANDOM 0
1001 #define DIST_NORMAL 1
1002
1003 double getRandom(double offset, double scale, int type)
1004 {
1005 int i;
1006 double randnum;
1007
1008 if(type == DIST_RANDOM)
1009         {
1010         randnum = random();
1011         return(scale * randnum/MY_RAND_MAX + offset);
1012         }
1013
1014 if(type == DIST_NORMAL)
1015         {
1016         randnum = 0.0;
1017         for(i = 0; i < 12; i++)
1018                 randnum += random() % 1024;
1019         return(scale * (randnum - 6144)/1024 + offset);
1020         }
1021
1022 return(0.0);
1023 }
1024
1025 CELL * randomDist(CELL * params, int type)
1026 {
1027 double randnum;
1028 size_t i;
1029 size_t n = 0;
1030 double scale = 1.0;
1031 double offset = 0.0;
1032 CELL * dist;
1033 CELL * cell;
1034
1035 if(params->type != CELL_NIL)
1036         {
1037         params = getFloat(params, &offset);
1038         params = getFloat(params, &scale);
1039         if(params->type != CELL_NIL)
1040                 getInteger(params, (UINT*)&n);
1041         }
1042
1043 if( n == 0)
1044         {
1045         randnum = getRandom(offset, scale, type);
1046         return(stuffFloat(&randnum));
1047         }
1048
1049 dist = getCell(CELL_EXPRESSION);
1050 randnum = getRandom(offset, scale, type);
1051 dist->contents = (UINT)stuffFloat(&randnum);
1052 cell = (CELL*)dist->contents;
1053 for(i = 1; i < n; i++)
1054         {
1055         randnum = getRandom(offset, scale, type);
1056         cell->next = stuffFloat(&randnum);
1057         cell = cell->next;
1058         }
1059
1060 return(dist);
1061 }
1062
1063
1064 CELL * p_normal(CELL * params)
1065 {
1066 return(randomDist(params, DIST_NORMAL));
1067 }
1068
1069 CELL * p_random(CELL * params)
1070 {
1071 return(randomDist(params, DIST_RANDOM));
1072 }
1073
1074 CELL * p_randomize(CELL * params)
1075 {
1076 CELL * list;
1077 CELL * cell;
1078 size_t length, i, j;
1079 CELL * * vector;
1080 int repetition = 0;
1081
1082 getListHead(params, &list);
1083
1084 if((length = listlen(list)) <= 1)
1085   {
1086   cell = getCell(CELL_EXPRESSION);
1087   cell ->contents = (UINT)copyList(list);
1088   return(cell);
1089   }
1090
1091 repetition = getFlag(params->next);
1092
1093 /* build index vector */
1094 cell = list;
1095 vector = allocMemory(length * sizeof(UINT));
1096 for(i = 0; i < length; i++)
1097         {
1098         vector[i] = copyCell(cell);
1099         vector[i]->next = (void *)i;
1100         cell = cell->next;
1101         }
1102
1103 /* reorganize randomly */
1104 RANDOMIZE:
1105 for(i = 0; i < (length - 1); i++)
1106   {
1107   j = random() % length;
1108   cell = vector[i];
1109   vector[i] = vector[j];
1110   vector[j] = cell;
1111   }
1112
1113 /* check that new sequence is different */
1114 if(!repetition)
1115     {
1116     for(i = 0; i < length; i++)
1117         if(vector[i]->next != (void *)i) break;
1118
1119     if(i == length)  goto RANDOMIZE;
1120     }
1121
1122 /* relink the list */
1123 cell = list = getCell(CELL_EXPRESSION);
1124 cell->contents  = (UINT)vector[0];
1125 cell = (CELL *)cell->contents;
1126 for(i = 1; i < length; i++)
1127         {
1128         cell->next = vector[i];
1129         cell = cell->next;
1130         }
1131 cell->next = nilCell;
1132 free(vector);
1133
1134 return(list);
1135 }
1136
1137 /* ---------------------------------------------------------------------
1138    probZ  - probability of normal z value
1139
1140    Adapted from a polynomial approximation in:
1141       Ibbetson D, Algorithm 209
1142       Collected Algorithms of the CACM 1963 p. 616
1143
1144    Note:
1145       This routine has six digit accuracy, so it is only useful for absolute
1146       z values < 6.  For z values >= to 6.0, poz() returns 0.0.
1147
1148    propChi2 - popbablitiy of CHI-2 for df degrees of freedom
1149
1150    Adapted from:
1151       Hill, I. D. and Pike, M. C.  Algorithm 299
1152       Collected Algorithms for the CACM 1967 p. 243
1153       Updated for rounding errors based on remark in
1154       ACM TOMS June 1985, page 185
1155
1156    critChi2 - Compute critical chi-square value for p and df
1157    critZ - Compute critical Z-value from p
1158
1159 */
1160
1161 double probChi2(double chi2, int df);
1162 double critChi2(double p, int df);
1163 double probZ(double z);
1164 double critZ(double p);
1165 double gammaln(double xx);
1166 double betai(double a, double b, double x);
1167 double gammap(double a, double x);
1168 double betacf(double a, double b, double x);
1169 static double gser(double a, double x, double gln);
1170 double gcf(double a, double x, double gln);
1171
1172 CELL * p_probabilityZ(CELL * params)
1173 {
1174 double z;
1175 double p;
1176
1177 getFloat(params, &z);
1178
1179 p = probZ(z);
1180 return(stuffFloat((double *)&p));
1181 }
1182
1183
1184 double probChi2(double chi2, int df)
1185 {
1186 return(1.0 - gammap(df/2, chi2/2));
1187 }
1188
1189 CELL * p_probabilityChi2(CELL * params)
1190 {
1191 double chi2;
1192 double df;
1193 double q;
1194
1195 params = getFloat(params, &chi2);
1196 getFloat(params, &df);
1197
1198 q = probChi2(chi2, df);
1199
1200 return(stuffFloat((double *)&q));
1201 }
1202
1203
1204 CELL * p_criticalChi2(CELL  * params)
1205 {
1206 double p;
1207 long df;
1208 double chi;
1209
1210 params = getFloat(params, &p);
1211 getInteger(params, (UINT *)&df);
1212
1213 chi = critChi2(p, df);
1214
1215 return(stuffFloat((double *)&chi));
1216 }
1217
1218
1219 CELL * p_criticalZ(CELL * params)
1220 {
1221 double p;
1222 double Z;
1223
1224 getFloat(params, &p);
1225 Z = critZ(p);
1226
1227 return(stuffFloat((double *)&Z));
1228 }
1229
1230
1231 #define Z_MAX 6.0  /* Maximum meaningful z value */
1232
1233 double probZ(double z)
1234 {
1235 double y, x, w;
1236
1237
1238 if (z == 0.0) x = 0.0;
1239 else
1240         {
1241         y = 0.5 * (z < 0.0 ? -z : z);
1242         if (y >= (Z_MAX * 0.5)) x = 1.0;
1243         else
1244                 {
1245                 if (y < 1.0)
1246                         {
1247                         w = y * y;
1248                         x = ((((((((0.000124818987 * w
1249                                 - 0.001075204047) * w + 0.005198775019) * w
1250                                 - 0.019198292004) * w + 0.059054035642) * w
1251                                 - 0.151968751364) * w + 0.319152932694) * w
1252                                 - 0.531923007300) * w + 0.797884560593) * y * 2.0;
1253                         }
1254                 else {
1255                         y -= 2.0;
1256                         x = (((((((((((((-0.000045255659 * y
1257                                 + 0.000152529290) * y - 0.000019538132) * y
1258                                 - 0.000676904986) * y + 0.001390604284) * y
1259                                 - 0.000794620820) * y - 0.002034254874) * y
1260                                 + 0.006549791214) * y - 0.010557625006) * y
1261                                 + 0.011630447319) * y - 0.009279453341) * y
1262                                 + 0.005353579108) * y - 0.002141268741) * y
1263                                 + 0.000535310849) * y + 0.999936657524;
1264                         }
1265                 }
1266         }
1267 return z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5);
1268 }
1269
1270 #define BIGX 20.0        /* max value to represent exp(x) */
1271
1272 double ex(double x)
1273 {
1274 return (x < -BIGX) ? 0.0 : exp(x);
1275 }
1276
1277
1278 double critChi2(double p, int df)
1279 {
1280 #define CHI_EPSILON 0.000001   /* Accuracy of critchi approximation */
1281 #define CHI_MAX 99999.0        /* Maximum chi-square value */
1282 double minchisq = 0.0;
1283 double maxchisq = CHI_MAX;
1284 double chisqval;
1285
1286 if (p <= 0.0) return maxchisq;
1287 else if (p >= 1.0) return 0.0;
1288
1289 chisqval = df / sqrt(p);    /* fair first value */
1290 while ((maxchisq - minchisq) > CHI_EPSILON)
1291         {
1292         if (gammap(df/2, chisqval/2) < p) minchisq = chisqval;
1293         else maxchisq = chisqval;
1294         chisqval = (maxchisq + minchisq) * 0.5;
1295         }
1296 return chisqval;
1297 }
1298
1299
1300 double critZ(double p)
1301 {
1302 #define Z_EPSILON 0.000001   /* Accuracy of critchi approximation */
1303 double minZ = 0.0;
1304 double maxZ = Z_MAX;
1305 double Zval;
1306
1307 if (p <= 0.0) return maxZ;
1308 else if (p >= 1.0) return 0.0;
1309
1310 Zval = 2.0;    /* fair first value */
1311 while ((maxZ - minZ) > Z_EPSILON)
1312         {
1313         if (probZ(Zval) < p) minZ = Zval;
1314         else maxZ = Zval;
1315         Zval = (maxZ + minZ) * 0.5;
1316         }
1317 return Zval;
1318 }
1319
1320
1321 /* ----------------------- betai and gammaln fucntions ----------------------*/
1322
1323
1324 static int paramError;
1325
1326 CELL * p_beta(CELL * params)
1327 {
1328 double a, b, beta;
1329
1330 params = getFloat(params, &a);
1331 getFloat(params, &b);
1332
1333 beta = exp(gammaln(a) + gammaln(b) - gammaln(a+b));
1334
1335 return(stuffFloat(&beta));
1336 }
1337
1338 CELL * p_betai(CELL * params)
1339 {
1340 double a, b, x, result;
1341
1342 params = getFloat(params, &x);
1343 params = getFloat(params, &a);
1344 getFloat(params, &b);
1345
1346 paramError = 0;
1347
1348 result = betai(a,b,x);
1349
1350 if(paramError)
1351         return(nilCell);
1352
1353 return(stuffFloat(&result));
1354 }
1355
1356
1357 CELL * p_gammaln(CELL * params)
1358 {
1359 double x, result;
1360
1361 getFloat(params, &x);
1362
1363 paramError = 0;
1364
1365 result = gammaln(x);
1366
1367 if(paramError)
1368         return(nilCell);
1369
1370 return(stuffFloat(&result));
1371 }
1372
1373
1374 CELL * p_gammai(CELL * params)
1375 {
1376 double a, x, result;
1377
1378 params = getFloat(params, &a);
1379 getFloat(params, &x);
1380
1381 result = gammap(a, x);
1382
1383 return(stuffFloat(&result));
1384 }
1385
1386
1387 /* these functions are adapted from:
1388    Numerical Recipes in C
1389    W.H.Press, S.A. Teukolsky, W.T. Vettering, B.P. Flannery
1390    Cambridge University Press (c) 1992
1391 */
1392
1393 #define ITMAX 100
1394 #define EPS7 3.0e-7
1395
1396 double betai(double a, double b, double x)
1397 {
1398 double bt;
1399
1400 if (x < 0.0 || x > 1.0)
1401         {
1402         paramError = 1;
1403         return(0.0);
1404         }
1405
1406 if (x == 0.0 || x == 1.0)
1407         bt = 0.0;
1408 else
1409         bt = exp(gammaln(a+b)-gammaln(a)-gammaln(b)+a*log(x)+b*log(1.0-x));
1410
1411 if (x < (a+1.0)/(a+b+2.0))
1412         return (bt * betacf(a,b,x) / a);
1413 else
1414         return (1.0 - bt * betacf(b,a,1.0-x) / b);
1415 }
1416
1417
1418 double betacf(double a, double b, double x)
1419 {
1420 double qap,qam,qab,em,tem,d;
1421 double bz,bm,bp,bpp;
1422 double az,am,ap,app,aold;
1423 int m;
1424
1425 bm = az = am = 1.0;
1426 qab=a+b;
1427 qap=a+1.0;
1428 qam=a-1.0;
1429 bz=1.0-qab*x/qap;
1430 for (m=1;m<=ITMAX;m++)
1431         {
1432         em=(double) m;
1433         tem=em+em;
1434         d=em*(b-em)*x/((qam+tem)*(a+tem));
1435         ap=az+d*am;
1436         bp=bz+d*bm;
1437         d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem));
1438         app=ap+d*az;
1439         bpp=bp+d*bz;
1440         aold=az;
1441         am=ap/bpp;
1442         bm=bp/bpp;
1443         az=app/bpp;
1444         bz=1.0;
1445         if (fabs(az-aold) < (EPS7 * fabs(az))) return az;
1446         }
1447
1448 return(paramError = 1);
1449 }
1450
1451
1452 double gammaln(double xx)
1453 {
1454 double x,y,tmp,ser;
1455 static double cof[6] = {
1456  76.18009172947146,
1457 -86.50532032941677,
1458  24.01409824083091,
1459  -1.231739572450155,
1460   0.1208650973866179e-2,
1461  -0.5395239384953e-5};
1462
1463 int j;
1464
1465 if(xx == 0.0)
1466         fatalError(ERR_INVALID_PARAMETER_0, NULL, 0);
1467
1468 y=x=xx;
1469 tmp=x+5.5;
1470 tmp -= (x+0.5)*log(tmp);
1471 ser=1.000000000190015;
1472
1473 for (j=0;j<=5;j++) ser += cof[j]/++y;
1474
1475 return -tmp+log(2.5066282746310005*ser/x);
1476 }
1477
1478
1479 #define EPS9 3.0e-9
1480 #define FPMIN 1.0e-307
1481
1482 /* incomplete Gamma function
1483 //
1484 // gammap = P(a,x)
1485 // gammaq = Q(a,x) = 1.0 - P(a,x)
1486 //
1487 // prob-chi2 = gammap(df/2, chi2/2)
1488 // of beeing as good as observed (>=)
1489 */
1490 double gammap(double a, double x)
1491 {
1492 double gln;
1493 gln = gammaln(a);
1494 return (x < a + 1) ? gser(a, x, gln) : (1.0 - gcf(a ,x , gln));
1495 }
1496
1497
1498 static double gser(double a, double x, double gln)
1499 {
1500 double ap, del, sum;
1501 int n;
1502
1503 ap = a;
1504 del = 1.0/a;
1505 sum = del;
1506
1507 for (n = 1 ; n <= ITMAX ; n++)
1508         {
1509         ++ap;
1510         del *= x / ap;
1511         sum += del;
1512         if (fabs(del) < fabs(sum) * EPS9)
1513         return sum * exp(-x + a * log(x) - gln);
1514         }
1515
1516 return sum * exp(-x + a * log(x) - gln);
1517 }
1518
1519
1520 double gcf(double a, double x, double gln)
1521 {
1522 double b, c, d, h, an, del;
1523 int i;
1524
1525 b = x + 1 - a;
1526 c = 1.0/FPMIN;
1527 d = 1.0/b;
1528 h = d;
1529
1530 for (i = 1 ; i <= ITMAX ; i++)
1531         {
1532         an = -i * (i - a);
1533         b += 2;
1534
1535         d = an * d + b;
1536         if (fabs(d) < FPMIN) d = FPMIN;
1537
1538         c = b + an/c;
1539         if (fabs(c) < FPMIN) c = FPMIN;
1540
1541         d = 1.0/d;
1542         del = d * c;
1543         h *= del;
1544         if (fabs(del-1.0) < EPS9) break;
1545         }
1546
1547 return exp(-x + a * log(x) - gln) * h;
1548 }
1549
1550
1551 /* ------------------------------------- Binomial distribution -------------------- */
1552
1553 CELL * p_binomial(CELL * params)
1554 {
1555 long n, k;
1556 double bico, p, binomial;
1557
1558 params = getInteger(params, (UINT *)&n);
1559 params = getInteger(params, (UINT *)&k);
1560 getFloat(params, &p);
1561
1562 bico = exp(gammaln(n + 1.0) - gammaln(k + 1.0) - gammaln(n - k + 1.0));
1563
1564 #ifdef TRU64
1565 binomial = bico * pow(p, (double)k) * pow(1.0 - p, (double)(n - k));
1566 #else
1567 binomial = bico * pow(p, k) * pow(1.0 - p, (n - k));
1568 #endif
1569
1570 return(stuffFloat(&binomial));
1571 }
1572
1573 /* -------------------------------------------------------------------------------- */
1574
1575 CELL * p_series(CELL * params)
1576 {
1577 double fromFlt, factFlt, current;
1578 ssize_t i;
1579 ssize_t count;
1580 CELL * sequence;
1581 CELL * cell;
1582
1583 params = getFloat(params, &fromFlt);
1584 params = getFloat(params, &factFlt);
1585 if(isnan(fromFlt) || isnan(factFlt))
1586         return(errorProc(ERR_INVALID_PARAMETER_NAN));
1587 getInteger(params, (UINT *)&count);
1588 sequence = getCell(CELL_EXPRESSION);
1589 if(count > 0)
1590         {
1591         cell = stuffFloat(&fromFlt);
1592         sequence->contents = (UINT)cell;
1593         current = fromFlt;
1594         for(i = 1; i < count; i++)
1595                 {
1596                 current *= factFlt;
1597                 cell->next = stuffFloat(&current);
1598                 cell = cell->next;
1599                 }
1600         }
1601
1602 return(sequence);
1603 }
1604
1605 /* ------------------------------- prime numbers ---------------------------- */
1606
1607 /*
1608 * adapted for newLISP from the following code:
1609 *
1610 * factor.c -- print prime factorization of a number
1611 * Ray Gardner -- 1985 -- public domain
1612 * Modified Feb. 1989 by Thad Smith > public domain
1613 *
1614 * This version takes numbers up to the limits of double precision.
1615 */
1616
1617
1618 CELL * pushFactor (INT64 d, int k, INT64 * prevFact, CELL * factList)
1619 {
1620 if (! *prevFact)
1621         {
1622         factList->contents = (UINT)stuffInteger64(d);
1623         factList = (CELL*)factList->contents;
1624         k--;
1625         }
1626
1627 while(k--)
1628         {
1629         factList->next = stuffInteger64(d);
1630         factList = factList->next;
1631         }
1632
1633 (*prevFact)++;
1634
1635 return(factList);
1636 }
1637
1638
1639 CELL * p_factor (CELL * params)
1640 {
1641 INT64 d, n;
1642 long k;
1643 INT64 prevFact = 0;
1644 CELL * factList;
1645 CELL * next;
1646
1647 getInteger64(params, &n);
1648
1649 d = n + 1;     /* test for roundoff error */
1650
1651 if ( (n + 3 != d + 2)  || (n < 2) )
1652       return(nilCell);
1653
1654 next = factList = getCell(CELL_EXPRESSION);
1655
1656 if ( n > 2 )
1657       {
1658       d = 2;
1659       for ( k = 0; (n % d) == 0; k++ )
1660               n /= d;
1661       if ( k )  next = pushFactor(d, k, &prevFact, next);
1662
1663       for ( d = 3; d * d <= n; d += 2 )
1664             {
1665             for ( k = 0; (n % d) == 0; k++ )
1666                 n /= d;
1667                 if ( k ) next = pushFactor(d, k, &prevFact, next);
1668             }
1669       }
1670
1671 if ( n > 1 )
1672       {
1673       if ( ! prevFact )
1674         factList->contents = (UINT)stuffInteger64(n);
1675       else  next = pushFactor(n, 1, &prevFact, next);
1676       }
1677
1678 return(factList);
1679 }
1680
1681 CELL * p_gcd(CELL * params)
1682 {
1683 INT64 m, n;
1684 INT64 t, r;
1685
1686 params = getInteger64(params, &m);
1687 if(m < 0) m = -m;
1688 n = m;
1689
1690 while(params != nilCell)
1691         {
1692         params = getInteger64(params, &n);
1693         if(n < 0) n = -n;
1694
1695         ITERATE_GCD:
1696         if (m < n) { t = m; m = n; n = t; }
1697         if(n == 0) {n = m; continue; } /* for (gcd x 0) => x */
1698         r = m % n;
1699         if (r == 0) { m = n; continue; }
1700         m = n; n = r;
1701         goto ITERATE_GCD;
1702         }
1703
1704 return(stuffInteger64(n));
1705 }
1706
1707
1708 /* ------------------------------- financial functions ---------------------- */
1709
1710
1711 CELL * p_pmt(CELL * params)
1712 {
1713 long nper, type;
1714 double rate, pv;
1715 double fv = 0.0;
1716 double pmt = 0.0;
1717 double inc;
1718
1719 params = getFloat(params, &rate);
1720 params = getInteger(params, (UINT *)&nper);
1721 params = getFloat(params, &pv);
1722 if(params != nilCell)
1723         {
1724         params = getFloat(params, &fv);
1725         getInteger(params, (UINT *)&type);
1726         }
1727 else type = 0;
1728
1729
1730 if(rate == 0)
1731     pmt = (-pv - fv) / nper;
1732 else
1733     {
1734     inc = pow(1 + rate, (double)nper);
1735     pmt = - (pv * inc + fv) / ((1 + rate * type) * ((inc - 1) / rate));
1736     }
1737
1738 return stuffFloat(&pmt);
1739 }
1740
1741
1742 CELL * p_pv(CELL * params)
1743 {
1744 long nper, type;
1745 double rate, pmt, pv;
1746 double fv = 0.0;
1747 double inc;
1748
1749 params = getFloat(params, &rate);
1750 params = getInteger(params, (UINT *)&nper);
1751 params = getFloat(params, &pmt);
1752
1753 if(params != nilCell)
1754         {
1755         params = getFloat(params, &fv);
1756         getInteger(params, (UINT *)&type);
1757         }
1758 else type = 0;
1759
1760 if(rate == 0)
1761     pv = - pmt * nper - fv;
1762 else
1763     {
1764     inc = pow(1 + rate, (double)nper);
1765     pv = (-pmt * (1 + rate * type) * ((inc - 1) / rate) - fv) / inc;
1766     }
1767
1768 return stuffFloat(&pv);
1769 }
1770
1771
1772 CELL * p_fv(CELL * params)
1773 {
1774 double rate, pmt, pv;
1775 long nper, type;
1776 double inc, fv;
1777
1778 params = getFloat(params, &rate);
1779 params = getInteger(params, (UINT *)&nper);
1780 params = getFloat(params, &pmt);
1781 params = getFloat(params, &pv);
1782 if(params != nilCell)
1783     getInteger(params, (UINT *)&type);
1784 else type = 0;
1785
1786
1787 if(rate == 0)
1788     fv = - pmt * nper - pv;
1789 else
1790     {
1791     inc = pow(1 + rate, (double)nper);
1792     fv = -pmt * (1 + rate * type) * ((inc - 1) / rate) - pv * inc;
1793     }
1794
1795 return stuffFloat(&fv);
1796 }
1797
1798
1799 CELL * p_nper(CELL * params)
1800 {
1801 double rate, pmt, pv;
1802 double fv = 0.0;
1803 long type;
1804 double R, c, nper;
1805
1806 params = getFloat(params, &rate);
1807 params = getFloat(params, &pmt);
1808 params = getFloat(params, &pv);
1809
1810 if(params != nilCell)
1811         {
1812         params = getFloat(params, &fv);
1813         getInteger(params, (UINT *)&type);
1814         }
1815 else type = 0;
1816
1817 if(rate == 0)
1818     nper = (-pv - fv) / pmt;
1819 else if(rate > -1.0)
1820     {
1821     R = 1.0 + rate;
1822     c = pmt * (type ? R : 1.0) / rate;
1823     nper = log((c - fv) / (c + pv)) / log(R);
1824     }
1825 else nper = sqrt(-1.0);   /* NaN */
1826
1827
1828 return stuffFloat(&nper);
1829 }
1830
1831
1832 CELL * p_npv(CELL * params)
1833 {
1834 CELL * list;
1835 double rate, cashFlow, fNum;
1836 int count;
1837
1838 params = getFloat(params, &rate);
1839 list = evaluateExpression(params);
1840 if(!isList(list->type))
1841     return(errorProcExt(ERR_LIST_EXPECTED, params));
1842 list = (CELL *)list->contents;
1843
1844 cashFlow = 0.0;
1845 count = 0;
1846 while(list !=  nilCell)
1847     {
1848 /*
1849     if(list->type == CELL_FLOAT)
1850         memcpy((void *)&fNum, (void *)&list->aux, sizeof(double));
1851     else if(list->type == CELL_LONG)
1852         fNum = (int)list->contents;
1853 */
1854         if(isNumber(list->type))
1855                 fNum = getFloatFromCell(list);
1856     else
1857         fNum = 0.0;
1858
1859     cashFlow += fNum / pow((1.0 + rate), (double)++count);
1860     list = list->next;
1861     }
1862
1863 return stuffFloat(&cashFlow);
1864 }
1865
1866
1867 /* Internal Rate of Return - irr
1868
1869    adapted from a C++ program by:
1870
1871    Bernt Arne Odegaard, 1999 http://finance.bi.no/~bernt
1872
1873    Note, that some data have multiple solutions, this search
1874    algorithm returns only one
1875
1876 */
1877
1878
1879 double cfPv(int N, int times[], double amounts[], double rate)
1880 {
1881 double PV = 0.0;
1882 int t;
1883
1884 for(t = 0; t < N; t++)
1885     PV += amounts[t] / pow((1.0 + rate), (double)times[t]);
1886
1887 return(PV);
1888 }
1889
1890
1891 #define PRECISION 1.0e-5
1892 #define MAX_ITERATIONS 50
1893 #define IRR_ERROR -1.0
1894
1895 double irr(int N, int times[], double amounts[], double x2)
1896 {
1897 double est1, est2;
1898 double x1 = 0.0;
1899 double f, rtb, dx;
1900 double x_mid, f_mid;
1901 int i;
1902
1903 est1 = cfPv(N, times, amounts,  x1);
1904 est2 = cfPv(N, times, amounts,  x2);
1905
1906 for(i=0; i<MAX_ITERATIONS; i++)
1907     {
1908     if(est1 * est2 < 0.0) break;
1909
1910     if(fabs(est1) < fabs(est2))
1911         est1 = cfPv(N, times, amounts,  x1 += 1.6 * (x1 - x2));
1912     else
1913         est2 = cfPv(N, times, amounts,  x2 += 1.6 * (x2 - x1));
1914     }
1915
1916 if (est2 * est1 > 0.0) return (IRR_ERROR);
1917
1918 f = cfPv(N, times, amounts, x1);
1919 dx=0;
1920
1921 if(f < 0.0)
1922     {
1923     rtb = x1;
1924     dx=x2-x1;
1925     }
1926 else
1927     {
1928     rtb = x2;
1929     dx = x1-x2;
1930     };
1931
1932 for(i = 0; i< MAX_ITERATIONS; i++)
1933     {
1934     dx *= 0.5;
1935     x_mid = rtb + dx;
1936     f_mid = cfPv(N, times, amounts, x_mid);
1937
1938     if(f_mid <= 0.0)
1939         rtb = x_mid;
1940
1941     if((fabs(f_mid) < PRECISION) || (fabs(dx) < PRECISION))
1942         return x_mid;
1943     };
1944
1945 return(IRR_ERROR);
1946 }
1947
1948 CELL * p_irr(CELL * params)
1949 {
1950 CELL * amounts;
1951 CELL * times;
1952 CELL * list;
1953 double result;
1954 double guess = 0.5;
1955 double * amountsVec;
1956 int * timesVec;
1957 int i, N = 0;
1958 long number;
1959
1960 params = getListHead(params, &amounts);
1961 list = amounts;
1962 if(params != nilCell)
1963         {
1964         params = getListHead(params, &times);
1965         if(params != nilCell)
1966                 getFloat(params, &guess);
1967         }
1968
1969 else
1970         times = NULL;
1971
1972 while(list != nilCell)
1973         {
1974         ++N;
1975         list = list->next;
1976         }
1977
1978 amountsVec = (double *)allocMemory(N * sizeof(double));
1979 timesVec = (int *)allocMemory(N * sizeof(int));
1980
1981 for(i = 0; i < N; i++)
1982         {
1983         if(isNumber(amounts->type))
1984                 {
1985                 amountsVec[i] = getFloatFromCell(amounts);
1986                 amounts = amounts->next;
1987                 }
1988         else
1989                 {
1990                 free(amountsVec);
1991                 free(timesVec);
1992                 return(errorProcExt(ERR_NUMBER_EXPECTED, amounts));
1993                 }
1994
1995         if(times == NULL)
1996                 timesVec[i] = i + 1;
1997         else
1998                 {
1999                 times = getIntegerExt(times, (UINT *)&number, FALSE);
2000                 timesVec[i] = number;
2001                 }
2002         }
2003
2004 result = irr(N, timesVec, amountsVec, guess);
2005 if(result < 0.00001)
2006         result = 0.0;
2007 else
2008         {
2009         result = floor((10000 * result + 0.5))/10000;
2010         }
2011
2012 free(amountsVec);
2013 free(timesVec);
2014
2015 if(result == IRR_ERROR)
2016         return(nilCell);
2017
2018 return(stuffFloat(&result));
2019 }
2020
2021 /* ----------------------------------- CRC32 ----------------------- */
2022
2023 /* Algorithm from: http://www.w3.org/TR/PNG-CRCAppendix.html */
2024
2025 unsigned int update_crc(unsigned int crc, unsigned char *buf, int len);
2026
2027 CELL * p_crc32(CELL * params)
2028 {
2029 char * data;
2030 size_t len;
2031
2032 params = getStringSize(params, &data, &len, TRUE);
2033 return(stuffInteger(update_crc(0xffffffffL, (unsigned char *)data, (int)len) ^ 0xffffffffL));
2034 }
2035
2036 /* Update a running CRC with the bytes buf[0..len-1]--the CRC
2037    should be initialized to all 1's, and the transmitted value
2038    is the 1's complement of the final running CRC (see the
2039    crc() routine below)).
2040 */
2041
2042 unsigned int update_crc(unsigned int crc, unsigned char *buf, int len)
2043 {
2044 unsigned int crc_table[256];
2045 unsigned int c;
2046 int n, k;
2047
2048 /* make crc table */
2049 for (n = 0; n < 256; n++)
2050     {
2051     c = (unsigned int) n;
2052     for (k = 0; k < 8; k++)
2053        {
2054        if (c & 1)
2055            c = 0xedb88320L ^ (c >> 1);
2056        else
2057            c = c >> 1;
2058        }
2059    crc_table[n] = c;
2060    }
2061
2062 c = crc;
2063
2064 /* update crc */
2065 for (n = 0; n < len; n++)
2066     c = crc_table[(c ^ buf[n]) & 0xff] ^ (c >> 8);
2067
2068 return c;
2069 }
2070
2071 /* ----------------------------------- Bayesian statistics -------------------------- */
2072
2073 /*
2074
2075    (bayes-train M1 M2 [M3 ... Mn] D)
2076
2077    takes two or more lists M1, M2 ... with tokens from a joint set of tokens T.
2078    Token can be symbols or strings other datatypes are ignored.
2079    Tokerns are placed in a common dictionary in context D and the frquency
2080    is counted  how often they occur in each category Mi.
2081
2082    The M categories represent data models for which a sequence of tokens can be
2083    classified see (bayes-query ...)
2084
2085    Each token in D is a content addressable symbol in D containing a
2086    list of ferquencies of that token how often it occurs in each category.
2087    String tokens are prepended with an undersocre _ before convering them
2088    to a symbol in D.
2089
2090    The function returns a list of token numbers in the different category models:
2091    (N-of-tokens-in-M1 N-of-tokens-in-M2)
2092
2093    (bayes-train M1 M2 [M3 ... Mn] D)
2094
2095 */
2096
2097 CELL * p_bayesTrain(CELL * params)
2098 {
2099 CELL * * category;
2100 CELL * list;
2101 CELL * count;
2102 int * total;
2103 int i, idx, maxIdx = 0;
2104 SYMBOL * ctx;
2105 SYMBOL * sPtr;
2106 SYMBOL * totalPtr;
2107 char * token;
2108
2109 list = params;
2110 while(list != nilCell) list = list->next, maxIdx++;
2111 --maxIdx; /* last is a context not a category */
2112 if(maxIdx < 1) errorProc(ERR_MISSING_ARGUMENT);
2113
2114 category = alloca(maxIdx * sizeof(CELL *));
2115 total = alloca(sizeof(int));
2116 token = alloca(MAX_STRING + 1);
2117
2118 for(idx = 0; idx < maxIdx; idx++)
2119   {
2120   params = getListHead(params, &list);
2121   category[idx] = list;
2122   }
2123
2124 if((ctx = getCreateContext(params, TRUE)) == NULL)
2125   return(errorProcExt(ERR_SYMBOL_OR_CONTEXT_EXPECTED, params));
2126
2127 for(idx = 0; idx < maxIdx; idx++)
2128   {
2129   list = category[idx];
2130   total[idx] = 0;
2131   while(list != nilCell)
2132     {
2133
2134     switch(list->type)
2135       {
2136       case CELL_STRING:
2137         if(list->aux > MAX_STRING) continue;
2138         *token = '_';
2139         memcpy(token + 1, (char *)list->contents, list->aux);
2140         break;
2141       case CELL_SYMBOL:
2142         strncpy(token, ((SYMBOL *)list->contents)->name, MAX_STRING);
2143         break;
2144       }
2145
2146     sPtr = translateCreateSymbol(token, CELL_EXPRESSION, ctx, TRUE);
2147     if(((CELL*)sPtr->contents)->type == CELL_NIL)
2148       {
2149       /* create list with maxIdx 0s */
2150       count = getCell(CELL_EXPRESSION);
2151       sPtr->contents = (UINT)count;
2152       count->contents = (UINT)stuffInteger(0);
2153       count = (CELL *)count->contents;
2154       for(i = 1; i < maxIdx; i++)
2155         {
2156         count->next = stuffInteger(0);
2157         count = count->next;
2158         }
2159       }
2160
2161     /* increment value at idx */
2162     count = (CELL *)sPtr->contents;
2163     if(count->type != CELL_EXPRESSION)
2164       return(errorProcExt(ERR_LIST_EXPECTED, count));
2165     count = (CELL *)count->contents;
2166     for(i = 0; i < idx; i++)
2167       count = count->next;
2168     if(count->type == CELL_LONG)
2169       count->contents++;
2170 #ifndef NEWLISP64
2171         else if(count->type == CELL_INT64)
2172           *(INT64 *)&count->aux += 1;
2173 #endif
2174         else
2175       return(errorProcExt(ERR_NUMBER_EXPECTED, count));
2176     total[idx]++;
2177     list = list->next;
2178     } /* done category idx */
2179   }
2180
2181 totalPtr = translateCreateSymbol("total", CELL_EXPRESSION, ctx, TRUE);
2182 if(((CELL *)totalPtr->contents)->type == CELL_NIL)
2183   {
2184   count = getCell(CELL_EXPRESSION);
2185   totalPtr->contents = (UINT)count;
2186   count->contents = (UINT)stuffInteger(0);
2187   count = (CELL *)count->contents;
2188   for(i = 1; i < maxIdx; i++)
2189     {
2190     count->next = stuffInteger(0);
2191     count = count->next;
2192     }
2193   }
2194
2195 count = (CELL *)totalPtr->contents;
2196 if(count->type != CELL_EXPRESSION)
2197     return(errorProcExt(ERR_LIST_EXPECTED, count));
2198 count = (CELL *)count->contents;
2199 for(i = 0; i < maxIdx; i++)
2200   {
2201   if(count->type == CELL_LONG)
2202         count->contents += total[i];
2203 #ifndef NEWLISP64
2204   else if(count->type == CELL_INT64)
2205         *(INT64 *)&count->aux += total[i];
2206 #endif
2207   count = count->next;
2208   }
2209
2210 return(copyCell((CELL *)totalPtr->contents));
2211 }
2212
2213 /*
2214    (bayes-query L D [trueChainedBayes])
2215    (bayes-query L D [trueChainedBayes] [trueFloatProbs)
2216
2217    Takes a list of tokens L and a trained dictionary D and returns a list of the
2218    combined probabilities of the tokens to be in one category A versus the other
2219    categories B. All token in L should occur in D, for tokens which are not in D
2220    equal probability is asssumed over categories.
2221
2222    Token can be strings or symbols. If token are strings or symbols they are
2223    prepended with an underscore when looked up in D. When frequencies in D wehere
2224    learned using bayes-train, the underscore was automatically prepended during
2225    learning.
2226
2227    for 2 categories:
2228
2229                      p(tkn|A) * p(A)
2230    p(A|tkn) = ---------------------------------
2231               p(tkn|A) * p(A) + p(tkn|B) * p(B)
2232
2233    p(A|tkn) is the posterior for A which gets prior p(A) for the next token
2234    the priors p(A) and p(B) = p(not A) = 1 are substituted after every token
2235
2236    for N categories:
2237
2238                     p(tkn|Mc) * p(Mc)
2239    p(Mc|tkn = -------------------------------    ; Mc is one of N categories
2240                 sum-i-N( p(tkn|Mi) * p(Mi) )
2241
2242
2243    When in chain Bayes mode p(Mc|tkn) is the posterior for Mc and replaces
2244    the prior p(Mc) for the next token. In chain Bayes mode tokens with 0 frequency
2245    in one category will effectiviely put the probability of this category to 0.
2246    This causes queries resulting in 0 probabilities for all categories to yield NaN values.
2247
2248    The pure chain Bayes mode is the more sensitive one for changes, when a token count of 0
2249    occurs the resulting probability goes to a complete 0 in this category, while other
2250    categories gain higher probabilities. In the Fisher's  Chi2 mode the the zero-category is
2251    still assigned a probability resulting from other tokens with non-zero counts.
2252
2253    Use same or simlar total in training categries when using Fisher Chi2.
2254
2255 */
2256
2257 CELL * p_bayesQuery(CELL * params)
2258 {
2259 CELL * tokens;
2260 CELL * tkn;
2261 char * token;
2262 SYMBOL * ctx;
2263 SYMBOL * sPtr;
2264 SYMBOL * totPtr;
2265 CELL * count;
2266 CELL * total;
2267 int i, idx;
2268 int nTkn = 0;   /* number of token in query */
2269 int maxIdx = 0; /* number of catagories */
2270 int N = 0;      /* total of tokens in all categories */
2271 int chainBayes, probFlag;
2272 double * priorP;
2273 double * postP;
2274 double * p_tkn_and_cat;
2275 double p_tkn_and_all_cat;
2276 double * Pchi2 = NULL;
2277 double * Qchi2 = NULL;
2278 double sumS = 1.0;
2279 long countNum, totalNum;
2280
2281 CELL * result = NULL;
2282 CELL * cell = NULL;
2283
2284 params = getListHead(params, &tokens);
2285 params = getContext(params, &ctx);
2286
2287 chainBayes = getFlag(params);
2288 probFlag = getFlag(params->next);
2289
2290 totPtr = translateCreateSymbol("total", CELL_EXPRESSION, ctx, FALSE);
2291 if(totPtr == NULL)  return(nilCell);
2292
2293 total = (CELL *)totPtr->contents;
2294 if(total->type != CELL_EXPRESSION)
2295     return(errorProcExt(ERR_LIST_EXPECTED, (CELL *)totPtr->contents));
2296
2297 /* get number of categories maxIdx and total N for when
2298    no probabilities are specified but frequencies/counts  */
2299 total = (CELL *)total->contents;
2300 while(total != nilCell)
2301     {
2302     if(probFlag == FALSE)
2303                 {
2304                 if(total->type == CELL_LONG)
2305                 N += total->contents;
2306 #ifndef NEWLISP64
2307                 else if(total->type == CELL_INT64)
2308                         N += *(INT64 *)&total->aux;
2309 #endif
2310                 }
2311     total = total->next;
2312     maxIdx++;
2313     }
2314
2315 priorP = alloca(maxIdx * sizeof(double));
2316 postP = alloca(maxIdx * sizeof(double));
2317 p_tkn_and_cat = alloca(maxIdx * sizeof(double));
2318 if(!chainBayes)
2319     {
2320     Pchi2 = alloca(maxIdx * sizeof(double));
2321     Qchi2 = alloca(maxIdx * sizeof(double));
2322     }
2323
2324 total = (CELL *)((CELL *)totPtr->contents)->contents;
2325 for(idx = 0; idx < maxIdx; idx++)
2326     {
2327         if(probFlag == TRUE)
2328 #ifndef NEWLISP64
2329         priorP[idx] = *(double *)&total->aux;
2330 #else
2331                 priorP[idx] = (double)total->contents;
2332 #endif
2333         else
2334                 {
2335                 if(total->type == CELL_LONG)
2336                         priorP[idx] = (double)total->contents / N;
2337 #ifndef NEWLISP64
2338                 else /* INT64 */
2339                         priorP[idx] = (double)*(INT64 *)&total->aux / N;
2340 #endif
2341                 }
2342
2343 /*    printf("priorP[%d]=%f\n", idx, priorP[idx]);  */
2344
2345     total = total->next;
2346     if(!chainBayes) Pchi2[idx] = Qchi2[idx] = 0.0;
2347     }
2348
2349 token = alloca(MAX_STRING + 1);
2350 tkn = tokens;
2351 while(tkn != nilCell) tkn = tkn->next, nTkn++;
2352
2353 tkn = tokens;
2354 for(i = 0; i < nTkn; i++)
2355     {
2356     switch(tkn->type)
2357       {
2358       case CELL_STRING:
2359         if(tkn->aux > MAX_STRING) continue;
2360         *token = '_';
2361         memcpy(token + 1, (char *)tkn->contents, tkn->aux);
2362         break;
2363       case CELL_SYMBOL:
2364         strncpy(token, ((SYMBOL *)tkn->contents)->name, MAX_STRING);
2365         break;
2366       }
2367
2368     if((sPtr = lookupSymbol((char *)token, ctx)) == NULL) continue;
2369
2370     count = (CELL *)sPtr->contents;
2371     if(count->type != CELL_EXPRESSION) continue;
2372
2373     count = (CELL *)(CELL *)count->contents;
2374     total = (CELL *)((CELL *)totPtr->contents)->contents;
2375
2376     p_tkn_and_all_cat = 0.0;
2377     for(idx = 0; idx < maxIdx; idx++)
2378         {
2379         if(probFlag)
2380 #ifndef NEWLISP64
2381                         p_tkn_and_cat[idx] = *(double *)&count->aux * priorP[idx];
2382 #else
2383                         p_tkn_and_cat[idx] = *(double *)&count->contents * priorP[idx];
2384 #endif
2385                 else /*   p(M) * p(tkn|M) */
2386                         {
2387                         getInteger(count, (UINT *)&countNum);
2388                         getInteger(total, (UINT *)&totalNum);
2389                 p_tkn_and_cat[idx] = priorP[idx] * countNum / totalNum;
2390                         }
2391
2392         p_tkn_and_all_cat += p_tkn_and_cat[idx];
2393 /*
2394         printf("token[%d] p(tkn(M) * p(tkn|M)[%d]=%lf prior[%d]=%lf\n", i, idx,
2395             (double)p_tkn_and_cat[idx], idx, priorP[idx]);
2396 */
2397
2398         count = count->next;
2399         total = total->next;
2400         }
2401
2402     for(idx = 0; idx < maxIdx; idx++)
2403         {
2404         if(!chainBayes)
2405             {
2406             postP[idx] = p_tkn_and_cat[idx] / p_tkn_and_all_cat;
2407
2408             priorP[idx] = 1.0 / maxIdx; /* will cancel out */
2409
2410             if(postP[idx] == 0.0)
2411                 Qchi2[idx] += log(3e-308) * -2.0;
2412             else
2413                 Qchi2[idx] += log(postP[idx]) * -2.0;
2414
2415             if(postP[idx] == 1.0)
2416                 Pchi2[idx] += log(3e-308) * -2.0;
2417             else
2418                 Pchi2[idx] += log(1.0 - postP[idx]) * -2.0;
2419             }
2420         else
2421             {
2422             postP[idx] = p_tkn_and_cat[idx] / p_tkn_and_all_cat;
2423             priorP[idx] = postP[idx];
2424             }
2425 /*
2426         printf("p_tkn_and_cat[%d] / p_tkn_and_all_cat = %lf / %lf = %lf\n",
2427           idx, p_tkn_and_cat[idx], p_tkn_and_all_cat, postP[idx]);
2428 */
2429         }
2430
2431     tkn = tkn->next;
2432     }
2433
2434 if(!chainBayes)
2435     {
2436     sumS = 0.0;
2437     for(idx = 0; idx < maxIdx; idx++)
2438       {
2439       postP[idx] = (probChi2(Qchi2[idx],  2 * nTkn) - probChi2(Pchi2[idx],  2 * nTkn) + 1.0) / 2.0;
2440       sumS += postP[idx];
2441       }
2442     }
2443
2444
2445 for(idx = 0; idx < maxIdx; idx++)
2446     {
2447     /* normalize probs from fisher's Chi2 */
2448     if(!chainBayes)
2449         postP[idx] /= sumS;
2450
2451     if(idx == 0)
2452         {
2453         result = getCell(CELL_EXPRESSION);
2454         result->contents = (UINT)stuffFloat(postP + idx);
2455         cell = (CELL *)result->contents;
2456         }
2457     else
2458         {
2459         cell->next = stuffFloat(postP + idx);
2460         cell = cell->next;
2461         }
2462     }
2463
2464 return(result);
2465 }
2466
2467 /*
2468 //
2469 // Copyright (C) 1992-2007 Lutz Mueller <lutz@nuevatec.com>
2470 //
2471 // This program is free software; you can redistribute it and/or modify
2472 // it under the terms of the GNU General Public License version 2, 1991,
2473 // as published by the Free Software Foundation.
2474 //
2475 // This program is distributed in the hope that it will be useful,
2476 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2477 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
2478 // GNU General Public License for more details.
2479 //
2480 // You should have received a copy of the GNU General Public License
2481 // along with this program; if not, write to the Free Software
2482 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
2483 //
2484 */
2485
2486 /*
2487  'unify' for Prolog like unification of s-expressions:
2488   (unify '(f (g A) A) '(f B xyz)) => binds A to xyz and B to (g xyz)
2489
2490   variable symbols must start with an upper-case letter, variables
2491   containing nil are considered unbound
2492
2493   after linear left to right unification each variable symbol is expanded
2494   by all other variable symbols
2495 */
2496
2497 #ifdef SUPPORT_UTF8
2498 #include <wctype.h>
2499 #endif
2500
2501 typedef struct
2502         {
2503         CELL * left;
2504         CELL * right;
2505         void * next;
2506         } TERMSET;
2507
2508 void pushSet(TERMSET * * root, CELL * left, CELL * right);
2509 CELL * unify(CELL * left, CELL * right);
2510 int unifyGetType(CELL * cell);
2511 void substitute(CELL * expr, CELL * sym, TERMSET * tset);
2512 CELL * subsym(CELL * expr, CELL * sym, CELL * cell);
2513 CELL * finishUnify(CELL * result);
2514 int occurCheck(CELL * symCell, CELL * expr);
2515 void printStack(TERMSET * tset);
2516 void freeTermSet(TERMSET * * tset);
2517
2518 TERMSET * mgu = NULL;
2519 TERMSET * ws = NULL;
2520
2521 int bindFlag;
2522
2523 #ifdef DEBUG
2524 int debugFlag;
2525 #endif
2526
2527 CELL * p_unify(CELL * params)
2528 {
2529 CELL * left, * right;
2530 CELL * envHead;
2531 left = evaluateExpression(params);
2532 params = params->next;
2533 right = evaluateExpression(params);
2534 params = params->next;
2535
2536 if(params != nilCell)
2537         {
2538         params = getListHead(params, &envHead);
2539         while(envHead != nilCell)
2540                 {
2541                 if(envHead->type != CELL_EXPRESSION)
2542                         return(errorProcExt(ERR_LIST_EXPECTED, envHead));
2543                 pushSet(&ws, copyCell((CELL*)envHead->contents),
2544                                         copyCell(((CELL*)envHead->contents)->next));
2545                 envHead = envHead->next;
2546                 }
2547         }
2548
2549 bindFlag =  getFlag(params);
2550
2551 #ifdef DEBUG
2552 debugFlag =  getFlag(params->next);
2553 if(debugFlag) printStack(ws);
2554 #endif
2555
2556 return(unify(left, right));
2557 }
2558
2559 #define UNIFY_ATOM 0
2560 #define UNIFY_LIST 1
2561 #define UNIFY_VAR 2
2562
2563 void pushSet(TERMSET * * root, CELL * left, CELL * right)
2564 {
2565 TERMSET * set;
2566
2567 set = (TERMSET *)callocMemory(sizeof(TERMSET));
2568
2569 set->left = left;
2570 set->right = right;
2571 set->next = *root;
2572 *root = set;
2573 }
2574
2575 void  popSet(TERMSET * * root, CELL * * left, CELL * * right)
2576 {
2577 TERMSET * set;
2578
2579 set = *root;
2580 *root = set->next;
2581
2582 *left = set->left;
2583 *right = set->right;
2584
2585 free(set);
2586 }
2587
2588
2589 CELL * unify(CELL * left, CELL * right)
2590 {
2591 int leftType, rightType;
2592 CELL * lCell, * rCell;
2593
2594 pushSet(&ws, copyCell(left), copyCell(right));
2595
2596 while(ws != NULL)
2597         {
2598         popSet(&ws, &left, &right);
2599
2600 #ifdef DEBUG
2601         if(debugFlag)
2602                 {
2603                 printf("unify:");
2604                 printCell(left, FALSE, OUT_CONSOLE);
2605                 printf(" ");
2606                 printCell(right, FALSE, OUT_CONSOLE);
2607                 printf("\n");
2608                 }
2609 #endif
2610
2611         leftType = unifyGetType(left);
2612         rightType = unifyGetType(right);
2613
2614         if( (leftType == UNIFY_ATOM && rightType == UNIFY_ATOM) ||
2615         (left->contents == right->contents))
2616                 {
2617                 if(compareCells(left, right) == 0)
2618                         {
2619                         deleteList(left);
2620                         deleteList(right);
2621                         continue;
2622                         }
2623                 else
2624                         {
2625                         deleteList(left);
2626                         deleteList(right);
2627                         return(finishUnify(nilCell));
2628                         }
2629                 }
2630
2631         if(leftType == UNIFY_VAR && !occurCheck(left, right))
2632                 {
2633                 substitute(right, left, mgu); /* expand(right-expr, left-sym) in mgu set */
2634                 substitute(right, left, ws); /* expand(right-expr, left-sym) in  ws set */
2635
2636 #ifdef DEBUG
2637                 if(debugFlag)
2638                         {
2639                         printf("ws stack\n");
2640                         printStack(ws);
2641                         }
2642 #endif
2643                 pushSet(&mgu, left, right);
2644 #ifdef DEBUG
2645                 if(debugFlag)
2646                         {
2647                         printf("mgu stack\n");
2648                         printStack(mgu);
2649                         }
2650 #endif
2651                 continue;
2652                 }
2653
2654         if(rightType == UNIFY_VAR && !occurCheck(right, left))
2655                 {
2656                 substitute(left, right, mgu); /* expand(left-expr, right-sym) in mgu set */
2657                 substitute(left, right, ws); /* expand(left-expr, right-sym) in ws set */
2658 #ifdef DEBUG
2659                 if(debugFlag)
2660                         {
2661                         printf("ws stack\n");
2662                         printStack(ws);
2663                         }
2664 #endif
2665                 pushSet(&mgu, right, left);
2666 #ifdef DEBUG
2667                 if(debugFlag)
2668                         {
2669                         printf("mgu stack\n");
2670                         printStack(mgu);
2671                         }
2672 #endif
2673                 continue;
2674                 }
2675
2676         if(leftType == UNIFY_LIST && rightType == UNIFY_LIST)
2677                 {
2678 #ifdef DEBUG
2679                 if(debugFlag)
2680                         {
2681                         printf("lists:");
2682                         printCell(left, FALSE, OUT_CONSOLE);
2683                         printf(" ");
2684                         printCell(right, FALSE, OUT_CONSOLE);
2685                         printf("\n");
2686                         }
2687 #endif
2688                 if(listlen((CELL*)left->contents) != listlen((CELL*)right->contents))
2689                         {
2690                         deleteList(left);
2691                         deleteList(right);
2692                         return(finishUnify(nilCell));
2693                         }
2694
2695                 lCell = (CELL*)left->contents;
2696                 rCell = (CELL*)right->contents;
2697                 while(lCell != nilCell)
2698                         {
2699                         pushSet(&ws, copyCell(lCell), copyCell(rCell));
2700                         lCell = lCell->next;
2701                         rCell = rCell->next;
2702                         }
2703                 deleteList(left);
2704                 deleteList(right);
2705                 continue;
2706                 }
2707
2708         deleteList(left);
2709         deleteList(right);
2710         return(finishUnify(nilCell));
2711         }
2712
2713 return(finishUnify(trueCell));
2714 }
2715
2716 int unifyGetType(CELL * cell)
2717 {
2718 SYMBOL * sPtr;
2719 #ifdef SUPPORT_UTF8
2720 int wchar;
2721 #endif
2722
2723 if(isSymbol(cell->type))
2724         {
2725     if(cell->type == CELL_SYMBOL)
2726         sPtr = (SYMBOL *)cell->contents;
2727     else
2728         sPtr = getDynamicSymbol(cell);
2729
2730 #ifndef SUPPORT_UTF8
2731         return((toupper(*sPtr->name) == *sPtr->name) ? UNIFY_VAR : UNIFY_ATOM);
2732 #else
2733     utf8_wchar(sPtr->name, &wchar);
2734         return((towupper(wchar) == wchar) ? UNIFY_VAR : UNIFY_ATOM);
2735 #endif
2736         }
2737 else if(isList(cell->type))
2738         return(UNIFY_LIST);
2739
2740 return(UNIFY_ATOM);
2741 }
2742
2743
2744 int occurCheck(CELL * symCell, CELL * expr)
2745 {
2746 CELL * cell;
2747
2748 if(!isEnvelope(expr->type))
2749         return(symCell->contents == expr->contents);
2750
2751 cell = (CELL*)expr->contents;
2752
2753 while(cell != nilCell)
2754         {
2755         if(symCell->contents == cell->contents) return(1);
2756         if(isEnvelope(cell->type) && occurCheck(symCell, cell)) return(1);
2757         cell = cell->next;
2758         }
2759
2760 return(0);
2761 }
2762
2763 void substitute(CELL * expr, CELL * sym, TERMSET * tset)
2764 {
2765 if(tset == NULL)
2766         {
2767 #ifdef DEBUG
2768         if(debugFlag)
2769                 {
2770                 printf("empty set in substitute %s for ", ((SYMBOL*)sym->contents)->name);
2771                 printCell(expr, FALSE, OUT_CONSOLE);
2772                 printf("\n");
2773                 }
2774 #endif
2775         return;
2776         }
2777
2778 while(tset != NULL)
2779     {
2780 #ifdef DEBUG
2781         if(debugFlag)
2782                 {
2783                 printf("substitute %s for ", ((SYMBOL*)sym->contents)->name);
2784                 printCell(expr, FALSE, OUT_CONSOLE);
2785                 printf("\n");
2786                 printf("left:");
2787                 printCell(tset->left, FALSE, OUT_CONSOLE);
2788                 }
2789 #endif
2790         tset->left = subsym(expr, sym, tset->left);
2791 #ifdef DEBUG
2792         if(debugFlag)
2793                 {
2794                 printf("->");
2795                 printCell(tset->left, FALSE, OUT_CONSOLE);
2796                 printf(" right:");
2797                 printCell(tset->right, FALSE, OUT_CONSOLE);
2798                 }
2799 #endif
2800         tset->right = subsym(expr, sym, tset->right);
2801 #ifdef DEBUG
2802         if(debugFlag)
2803                 {
2804                 printf("->");
2805                 printCell(tset->right, FALSE, OUT_CONSOLE);
2806                 printf("\n");
2807                 }
2808 #endif
2809     tset = tset->next;
2810     }
2811 }
2812
2813
2814 CELL * subsym(CELL * expr, CELL * sym, CELL * cell)
2815 {
2816 SYMBOL * sPtr;
2817 CELL * sCell;
2818
2819 sPtr = (SYMBOL *)sym->contents;
2820 sCell = (CELL *)sPtr->contents;
2821 sPtr->contents = (UINT)expr;
2822
2823 if(isEnvelope(cell->type))
2824     {
2825     expr =  expand(copyCell(cell), sPtr);
2826         deleteList(cell);
2827     }
2828 else
2829     {
2830     if(cell->contents == (UINT)sPtr)
2831                 {
2832                 expr = copyCell(expr);
2833                 deleteList(cell);
2834                 }
2835         else
2836                 expr = cell;
2837     }
2838
2839 sPtr->contents = (UINT)sCell;
2840 return(expr);
2841 }
2842
2843 CELL * finishUnify(CELL * result)
2844 {
2845 CELL * left, * right;
2846 CELL * list = NULL;
2847 CELL * assoc,  * cell = NULL;
2848 SYMBOL * sPtr;
2849
2850 if(result == nilCell)
2851         {
2852         freeTermSet(&mgu);
2853         freeTermSet(&ws);
2854         return(nilCell);
2855         }
2856
2857 /* make an association list of all bound variables */
2858 while(mgu != NULL)
2859         {
2860         popSet(&mgu, &left, &right);
2861
2862         if(!bindFlag)
2863                 {
2864                 assoc = getCell(CELL_EXPRESSION);
2865                 assoc->contents = (UINT)left;
2866                 left->next = right;
2867                 if(list == NULL)
2868                         {
2869                         list = getCell(CELL_EXPRESSION);
2870                         list->contents = (UINT)assoc;
2871                         }
2872                 else
2873                         cell->next = assoc;
2874                 cell = assoc;
2875                 }
2876         else
2877                 {
2878                 sPtr = (SYMBOL*)left->contents;
2879                 if(isProtected(sPtr->flags))
2880                         return(errorProcExt2(ERR_SYMBOL_PROTECTED, stuffSymbol(sPtr)));
2881                 sPtr->contents = (UINT)right;
2882                 }
2883
2884         }
2885
2886 freeTermSet(&ws);
2887
2888 if(!list || bindFlag) return(getCell(CELL_EXPRESSION));
2889 return(list);
2890 }
2891
2892
2893 void freeTermSet(TERMSET * * tset)
2894 {
2895 TERMSET * set;
2896
2897 while(*tset != NULL)
2898         {
2899         set = *tset;
2900         deleteList(set->left);
2901         deleteList(set->right);
2902         *tset = set->next;
2903         free(set);
2904         }
2905 }
2906
2907 #ifdef DEBUG
2908 void printStack(TERMSET * tset)
2909 {
2910 while(tset != NULL)
2911         {
2912         printCell(tset->left, FALSE, OUT_CONSOLE);
2913         printf("<->");
2914         printCell(tset->right, FALSE, OUT_CONSOLE);
2915         printf("\n");
2916         tset = tset->next;
2917         }
2918 }
2919 #endif
2920
2921 /* eof */