src/SyrComponent/SyrComponent_CheckOfUndefined.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 //  SuperVisionTest SyrComponent : example of component performing some mathinatical operations
  24 //  File   : SyrComponent_CheckOfUndefined.cxx
  25 //  Module : SuperVisionTest
  26 //
  27 #include <stdio.h>
  28 #include <unistd.h>
  29 #include <iostream>
  30 #include <fstream>
  31 #include <strstream>
  32 #include <string>
  33 #include <math.h>
  34 #include <cstring>
  35
  36 using namespace std;
  37
  38 // ------------------------------------------------------------------
  39 // NextPrime : Compute the first prime number greater or equal than an integer
  40 // ------------------------------------------------------------------
  41
  42 #define VALUESNBR 4
  43
  44 long NextPrime (const long me )
  45 {
  46
  47  struct svalue {int signiaib ;
  48                 int nbr ;} ;
  49
  50  struct svalue values[VALUESNBR] ;
  51  long ia ;
  52  long maxia ;
  53  long ib[4] ;
  54  int n[4] ;
  55 // int signiaib[4] = { -1 , +1 , +1 , -1 } ;
  56  int signiaib[4];
  57  signiaib[0] = -1;
  58  signiaib[1] = 1;
  59  signiaib[2] = 1;
  60  signiaib[3] = -1;
  61  long remain ;
  62
  63  int nbvalues ;
  64  int loop ;
  65  int nindd ;
  66  long minn ;
  67  long maxvn ;
  68  long premret = 0 ;
  69
  70    if (me < 0 || me >
  71 #ifdef DECOSF1
  72                       127149130704178201
  73 #else
  74                       2147483647
  75 #endif
  76                                ){
  77 //      Standard_RangeError::
  78 //       Raise("Try to apply NextPrime method with negative, null or too large value.");
  79         return 0 ;
  80    }
  81
  82    if ( me <= 7 ) {
  83      if ( me <= 1 )
  84        return 1 ;
  85      else if ( me <= 2 )
  86        return 2 ;
  87      else if ( me <= 3 )
  88        return 3 ;
  89      else if ( me <= 5 )
  90        return 5 ;
  91      else if ( me <= 7 )
  92        return 7 ;
  93    }
  94
  95    minn = ( me - 1 ) / 6 ;              // n minimum
  96    while ( 6*minn+1 < me ) {
  97         minn += 1 ;
  98       }
  99
 100    maxia = long( sqrt( (double)me ) / 6 + 1 ) ;
 101
 102    maxvn = minn + VALUESNBR ;
 103
 104    nbvalues = 0 ;
 105    for ( nindd = 0 ; nindd < VALUESNBR ; nindd++ ) {
 106       if ( 6*(nindd+minn)-1 < me ) {
 107         values[nindd].nbr = 1 ;
 108         values[nindd].signiaib = -1 ;
 109         nbvalues += 1 ;
 110       }
 111       else {
 112         values[nindd].nbr = 0 ;
 113         values[nindd].signiaib = 0 ;
 114       }
 115     }
 116
 117    for ( ia = 1 ; ia <= maxia ; ia++ ) {
 118       if ( nbvalues == VALUESNBR*2 ) {
 119         break ;
 120       }
 121       remain = -VALUESNBR ;
 122       ib[0] = ( minn + ia - remain ) / (6*ia - 1) ;
 123       n[0] = int ( 6*ia*ib[0] - ia - ib[0] - minn ) ;
 124       ib[1] = ( minn - ia - remain ) / (6*ia - 1) ;
 125       n[1] = int ( 6*ia*ib[1] + ia - ib[1] - minn ) ;
 126       ib[2] = ( minn + ia - remain ) / (6*ia + 1) ;
 127       n[2] = int ( 6*ia*ib[2] - ia + ib[2] - minn ) ;
 128       ib[3] = ( minn - ia - remain ) / (6*ia + 1) ;
 129       n[3] = int ( 6*ia*ib[3] + ia + ib[3] - minn ) ;
 130       for ( loop = 0 ; loop < 4 ; loop++ ) {
 131          if ( n[loop] >= 0 && n[loop] < VALUESNBR ) {
 132            if ( ( values[n[loop]].nbr == 0 ) ||
 133                 ( values[n[loop]].signiaib == signiaib[loop] ) ) {
 134              values[n[loop]].signiaib = -signiaib[loop] ;
 135              values[n[loop]].nbr += 1 ;
 136              if ( values[n[loop]].nbr <= 2 )
 137                nbvalues += 1 ;
 138            }
 139          }
 140        }
 141     }
 142    for ( nindd = 0 ; nindd < VALUESNBR ; nindd++ ) {
 143      if ( values[nindd].nbr == 0 ) {
 144        if ( me <= 6*(nindd+minn)-1 ) {
 145           premret = 6*(nindd+minn)-1 ;
 146           break ;
 147         }
 148        else if ( me <= 6*(nindd+minn)+1 ) {
 149          premret = 6*(nindd+minn)+1 ;
 150          break ;
 151        }
 152      }
 153      else if ( values[nindd].nbr == 1 ) {
 154         if ( values[nindd].signiaib > 0 ) {
 155           if ( me <= 6*(nindd+minn)-1 ) {
 156             premret = 6*(nindd+minn)-1 ;
 157             break ;
 158           }
 159         }
 160         else {
 161           if ( me <= 6*(nindd+minn)+1 ) {
 162             premret = 6*(nindd+minn)+1 ;
 163             break ;
 164           }
 165         }
 166       }
 167    }
 168
 169    if ( premret != 0 ) {
 170      return premret ;
 171    }
 172
 173  return NextPrime ( 6*(maxvn-1)+2) ;
 174
 175 }
 176
 177 static char * base( int b , long n ) {
 178   static char c[40] ;
 179   char ccc[2] ;
 180   int i = 39 ;
 181   c[39] = '\0' ;
 182   c[38] = '0' ;
 183   if ( n ) {
 184     while ( n ) {
 185       i-- ;
 186       sprintf( ccc ,"%d" , ( n % b ) ) ;
 187       c[i] = ccc[0] ;
 188       n = n / b ;
 189     }
 190   }
 191   else
 192     i = 38 ;
 193   return &c[i] ;
 194 }
 195
 196 static char * B2( long n ) {
 197   return base( 2 , n ) ;
 198 }
 199
 200 static char * B3( long n ) {
 201   return base( 3 , n ) ;
 202 }
 203
 204 static char * B4( long n ) {
 205   return base( 4 , n ) ;
 206 }
 207
 208 static char * B5( long n ) {
 209   return base( 5 , n ) ;
 210 }
 211
 212 static char * B7( long n ) {
 213   return base( 7 , n ) ;
 214 }
 215
 216 static char * B8( long n ) {
 217   return base( 8 , n ) ;
 218 }
 219
 220 static long S( long f , long n ) {
 221   int i ;
 222   long res = 0 ;
 223   for ( i = 0 ; i <= n ; i++ ) {
 224     res += 1 << (f*i) ;
 225   }
 226   return res ;
 227 }
 228
 229 static long S2( long n ) {
 230   return S( 2 , n ) ;
 231 }
 232
 233 static long S6( long n ) {
 234   return S( 6 , n ) ;
 235 }
 236
 237 static long rP( long odd ) {
 238   if ( (( odd - 1 ) % 6 ) == 0 )
 239     return 1 ;
 240   if ( (( odd - 3 ) % 6 ) == 0 )
 241     return 3 ;
 242   if ( (( odd - 5 ) % 6 ) == 0 )
 243     return 5 ;
 244   return 0 ;
 245 }
 246 static long nP( long odd ) {
 247   return ( odd / 6 ) ;
 248 }
 249
 250 static long rQ( long N ) {
 251   if ( ( N & 7 ) == 5 ) {
 252     return 5 ;
 253   }
 254   return ( N & 3 ) ;
 255 }
 256 static long nQ( long N ) {
 257   if ( ( N & 3 ) == 3 ) {
 258     return ( N / 4 ) ;
 259   }
 260   else {
 261     return ( N / 8 ) ;
 262   }
 263 }
 264
 265 static long rT( long N ) {
 266   return ( nP( N ) % 3 ) ;
 267 }
 268 static long nT( long N ) {
 269   return ( nP( N ) / 3 ) ;
 270 }
 271
 272 static char * fact( int n ) {
 273   static char chn[132] ;
 274   if ( n == 0 ) {
 275     strcpy( chn , "0" ) ;
 276   }
 277   else if ( n == 1 ) {
 278     strcpy( chn , "1" ) ;
 279   }
 280   else {
 281     int pow2 = 0 ;
 282     int pow3 = 0 ;
 283     while ( (n & 1 ) == 0 ) {
 284       pow2 += 1 ;
 285       n = (n >> 1 ) ;
 286     }
 287     while ( ( n % 3 ) == 0 ) {
 288       pow3 += 1 ;
 289       n = n / 3 ;
 290     }
 291     int pos = 0 ;
 292     if ( pow2 ) {
 293       if ( pow2 == 1 ) {
 294         sprintf( &chn[pos] , "%d" , 2 ) ;
 295       }
 296       else {
 297         sprintf( &chn[pos] , "2**%d" , pow2 ) ;
 298       }
 299       pos = strlen( chn ) ;
 300     }
 301     if ( pow3 ) {
 302       if ( pow2 ) {
 303         strcat( chn , "*" ) ;
 304         pos = strlen( chn ) ;
 305       }
 306       if ( pow3 == 1 ) {
 307         sprintf( &chn[pos] , "%d" , 3 ) ;
 308       }
 309       else {
 310         sprintf( &chn[pos] , "3**%d" , pow3 ) ;
 311       }
 312       pos = strlen( chn ) ;
 313     }
 314     if ( n != 1 ) {
 315       if ( pow2 || pow3 ) {
 316         sprintf( &chn[pos] , "*%d" , n ) ;
 317       }
 318       else {
 319         sprintf( &chn[pos] , "%d" , n ) ;
 320       }
 321     }
 322   }
 323   return chn ;
 324 }
 325
 326 static long Next( const long prev , long &pow2 ) {
 327   long next = ( 3*prev + 1 )/2 ;
 328   pow2 = 1 ;
 329   while ( ( next & 1 ) == 0 ) {
 330     next = ( next >> 1 ) ;
 331     pow2 += 1 ;
 332   }
 333   return next ;
 334 }
 335
 336 int T1a( int n0 , int n1 ) {
 337   return (1<<(6*n1+2))*(3*n0+0) + (1<<(6*n1+1)) + S2(3*n1+0) ;
 338 }
 339 int T1b( int n0 , int n1 ) {
 340   return (1<<(6*n1+6))*(3*n0+1) + (1<<(6*n1+5)) + S2(3*n1+2) ;
 341 }
 342 int T1c( int n0 , int n1 ) {
 343   return (1<<(6*n1+4))*(3*n0+2) + (1<<(6*n1+3)) + S2(3*n1+1) ;
 344 }
 345 int T2a( int n0 , int n1 ) {
 346   return (1<<(6*n1+7))*(3*n0+0) + S2(3*n1+2) ;
 347 }
 348 int T2b( int n0 , int n1 ) {
 349   return (1<<(6*n1+3))*(3*n0+1) + S2(3*n1+0) ;
 350 }
 351 int T2c( int n0 , int n1 ) {
 352   return (1<<(6*n1+5))*(3*n0+2) + S2(3*n1+1) ;
 353 }
 354
 355 int P1a1( int n0 , int n1 ) {
 356   return (1<<(6*n1+4))*(3*n0+0) + (1<<(6*n1+3)) + S2(3*n1+1) ;
 357 }
 358 int P1a2( int n0 , int n1 ) {
 359   return (1<<(6*n1+6))*(3*n0+0) + (1<<(6*n1+5)) + S2(3*n1+2) ;
 360 }
 361 int P1b1( int n0 , int n1 ) {
 362   return (1<<(6*n1+2))*(3*n0+1) + (1<<(6*n1+1)) + S2(3*n1+0) ;
 363 }
 364 int P1b2( int n0 , int n1 ) {
 365   return (1<<(6*n1+4))*(3*n0+1) + (1<<(6*n1+3)) + S2(3*n1+1) ;
 366 }
 367 int P1c1( int n0 , int n1 ) {
 368   return (1<<(6*n1+2))*(3*n0+2) + (1<<(6*n1+1)) + S2(3*n1+0) ;
 369 }
 370 int P1c2( int n0 , int n1 ) {
 371   return (1<<(6*n1+6))*(3*n0+2) + (1<<(6*n1+5)) + S2(3*n1+2) ;
 372 }
 373
 374 int P2a1( int n0 , int n1 ) {
 375   return (1<<(6*n1+3))*(3*n0+0) + S2(3*n1+0) ;
 376 }
 377 int P2a2( int n0 , int n1 ) {
 378   return (1<<(6*n1+5))*(3*n0+0) + S2(3*n1+1) ;
 379 }
 380 int P2b1( int n0 , int n1 ) {
 381   return (1<<(6*n1+5))*(3*n0+1) + S2(3*n1+1) ;
 382 }
 383 int P2b2( int n0 , int n1 ) {
 384   return (1<<(6*n1+7))*(3*n0+1) + S2(3*n1+2) ;
 385 }
 386 int P2c1( int n0 , int n1 ) {
 387   return (1<<(6*n1+3))*(3*n0+2) + S2(3*n1+0) ;
 388 }
 389 int P2c2( int n0 , int n1 ) {
 390   return (1<<(6*n1+7))*(3*n0+2) + S2(3*n1+2) ;
 391 }
 392
 393 typedef int (*fct)(int,int);
 394
 395 ostream * _Trace ;
 396
 397 void Compute( char * fctname , fct afct , int PT , int PP , int PPNext ) {
 398   int n0 ;
 399   int n1 ;
 400   long value = 0 ;
 401   long next = 0 ;
 402   for ( n1 = 0 ; n1 <= 3 ; n1++ ) {
 403     bool pn1 = true ;
 404     for ( n0 = 0 ; n0 <= 500 && value < 3000 ; n0++ ) {
 405       value = afct( n0 , n1 ) ;
 406       long pow2 ;
 407       next = Next( value , pow2 ) ;
 408       long kp2 ;
 409       if ( ( rP( value ) == 3 && rP( Next( next , kp2 ) ) == PPNext ) ||
 410              ( rT( value ) == PT && rP( value ) == PP &&
 411                rP( next ) == PPNext ) ) {
 412         if ( rQ( value ) == 3 ||
 413              ( rQ( value ) == 1 && ( nQ( value ) & 1 ) == 0 ) ) {
 414 //          if ( pn1 ) {
 415 //            *_Trace << "  n1 " << n1 ;
 416 //            *_Trace << " n0 " << n0 << " " ;
 417 //            pn1 = false ;
 418 //          }
 419 //          else {
 420             if ( value < 10 )
 421               *_Trace << "000" ;
 422             else if ( value < 100 )
 423               *_Trace << "00" ;
 424             else if ( value < 1000 )
 425               *_Trace << "0" ;
 426             *_Trace << value << " " ;
 427             long lownext = (next-1)/4 ;
 428             if ( 4*lownext+1 == next && ( lownext & 1 ) == 1 ) {
 429               *_Trace << "---" ;
 430             }
 431             *_Trace << fctname << " n0 " << n0 << " " ;
 432 //          }
 433           *_Trace << value << " = 6[3*" << nT( value ) << "+"
 434                   << rT( value ) << "]+" << rP( value ) << " --> "
 435                   << next << " (>>" << pow2 << ") = 6[3*"
 436                   << nT( next ) << "+" << rT( next )
 437                   << "] + " << rP( next ) << endl ;
 438         }
 439 #if 0
 440         else {
 441           if ( pn1 ) {
 442             *_Trace << "  n1 " << n1 ;
 443             *_Trace << " n0 " << n0 << " " ;
 444             pn1 = false ;
 445           }
 446           else {
 447             *_Trace << "       n0 " << n0 << " " ;
 448           }
 449           *_Trace << "---" << value << " = 6[3*" << nT( value ) << "+"
 450                   << rT( value ) << "]+" << rP( value ) << " = 4*["
 451                   << "2*" << nQ( value )/2 << "+1]+1" ;
 452           if ( Next( (value-1)/4 , kp2 ) != next )
 453             *_Trace << " " << Next( (value-1)/4 , kp2 ) << " != " << next ;
 454           *_Trace << endl ;
 455         }
 456 #endif
 457       }
 458     }
 459     value = 0 ;
 460   }
 461 }
 462
 463 int SeqSyr( int n ) {
 464   long pow2 ;
 465   int nstep = 0 ;
 466   int EQ ;
 467   int r ;
 468   int q ;
 469   int Sr ;
 470   int Sq ;
 471   int nn = n ;
 472   *_Trace << "   " << n ;
 473   while ( n != 1 ) {
 474     n = Next( n , pow2 ) ;
 475     nstep += 1 ;
 476     int EQ ;
 477     r = rQ( n ) ;
 478     q = nQ( n ) ;
 479     if ( r == 3 ) {
 480       EQ = 4 ;
 481     }
 482     else {
 483       EQ = 8 ;
 484     }
 485     if ( ( nstep % 2 ) == 0 ) {
 486       *_Trace << endl << "   " ;
 487     }
 488     Sr = rP( n ) ;
 489     Sq = nP( n ) ;
 490     *_Trace << " -> " << n << " = " << EQ << "*" << q << "+" << r
 491             << " = 6*" << Sq << "+" << Sr ;
 492     if ( r == 5 ) {
 493       break ;
 494     }
 495   }
 496   *_Trace << endl ;
 497   if ( r == 5 && q > nn/8 ) {
 498     return n ;
 499   }
 500   return 0 ;
 501 }
 502
 503 void little( int depth , int prt , int n , int r ) {
 504
 505   int n0 ;
 506   int r0 ;
 507   int n1 ;
 508   int r1 ;
 509   int r2 ;
 510   depth += 1 ;
 511
 512   if ( depth == prt )
 513     *_Trace << endl << depth << " little2 n = " << n << " r = " << r << endl ;
 514
 515   n0 = n*2 ;
 516   r0 = r*2 - 1 ;
 517   if ( ( r0 % 3 ) == 0 ) {
 518     r1 = r0/3 ;
 519     if ( depth == prt ) {
 520       *_Trace << "n->3n " << n0 << "n + " << r1 << " = 4[ "
 521               << n0/4 << "n + " << r1/4 << " ]+3" << endl ;
 522       SeqSyr( r1 ) ;
 523     }
 524   }
 525   else if ( ( ( n0 + r0 ) % 3 ) == 0 ) {
 526     r1 = ( n0 + r0 ) / 3 ;
 527     if ( depth == prt ) {
 528       *_Trace << "n->3n+1 " << n0 << "n + " << r1 << " = 4[ "
 529               << n0/4 << "n + " << r1/4 << " ]+3" << endl ;
 530       SeqSyr( r1 ) ;
 531     }
 532   }
 533   else if ( ( ( 2*n0 + r0 ) % 3 ) == 0 ) {
 534     r1 = ( 2*n0 + r0 ) / 3 ;
 535     if ( depth == prt ) {
 536       *_Trace << "n->3n+2 " << n0 << "n + " << r1 << " = 4[ "
 537               << n0/4 << "n + " << r1/4 << " ]+3" << endl ;
 538       SeqSyr( r1 ) ;
 539     }
 540   }
 541
 542   if ( depth == prt )
 543     *_Trace << endl << depth << " little4 n = " << n << " r = " << r << endl ;
 544
 545   n1 = n0*2 ;
 546   r2 = r0*2 + 1 ;
 547   if ( ( r2 % 3 ) == 0 ) {
 548     r2 = r2/3 ;
 549     if ( depth == prt ) {
 550       *_Trace << "n->3n " << n1 << "n + " << r2 << " = 8[ "
 551               << n1/8 << "n + " << r2/8 << " ]+1" << endl ;
 552       SeqSyr( r2 ) ;
 553     }
 554   }
 555   else if ( ( ( n1 + r2 ) % 3 ) == 0 ) {
 556     r2 = ( n1 + r2 ) / 3 ;
 557     if ( depth == prt ) {
 558       *_Trace << "n->3n+1 " << n1 << "n + " << r2 << " = 8[ "
 559               << n1/8 << "n + " << r2/8 << " ]+1" << endl ;
 560       SeqSyr( r2 ) ;
 561     }
 562   }
 563   else if ( ( ( 2*n1 + r2 ) % 3 ) == 0 ) {
 564     r2 = ( 2*n1 + r2 ) / 3 ;
 565     if ( depth == prt ) {
 566       *_Trace << "n->3n+2 " << n1 << "n + " << r2 << " = 8[ "
 567               << n1/8 << "n + " << r2/8 << " ]+1" << endl ;
 568       SeqSyr( r2 ) ;
 569     }
 570   }
 571
 572   if ( depth < prt ) {
 573     little( depth , prt , n0 , r1 ) ;
 574     little( depth , prt , n1 , r2 ) ;
 575   }
 576
 577 }
 578
 579 int main(int argc, char **argv) {
 580
 581   int k ;
 582   int i ;
 583   int j ;
 584   int l ;
 585   int l1 ;
 586   int n = 0 ;
 587   int p ;
 588   int P ;
 589   int Prem = 37 ;
 590
 591   _Trace = new ofstream( "/cas_01/jr/linux/SyrP.log" );
 592
 593   for ( l = 0 ; l <= 17 ; l++ ) {
 594     for ( l1 = -1 ; l1 <= 1 ; l1 += 2 ) {
 595       P = 6*l + l1 ;
 596       if ( P > 0 && NextPrime( P-1 ) == P ) {
 597         for ( k = -1 ; k <= 1 ; k += 2 ) {
 598           for ( i = 1 ; i <= 30 ; i++ ) {
 599             for ( j = 1 ; j <= 30 ; j++ ) {
 600               n = int( pow( (double)2 , (double)i ) * pow( (double)3 , (double)j )*P ) + k ;
 601               if ( n < 0 || n > pow( (double)2 , (double)30 ) ) {
 602                 break ;
 603               }
 604               if ( NextPrime( n-1 ) == n ) {
 605                 *_Trace << n << " = 2**" << i << " * 3**" << j << " * " << P << " + "
 606                         << k << endl ;
 607               }
 608 #if 0
 609               if ( NextPrime( n-1 ) != n && ( n % P ) == 0 ) {
 610                 *_Trace << n << " = 2**" << i << " * 3**" << j << " * " << P << " + "
 611                         << k << " = " ;
 612                 p = 5 ;
 613                 while ( p*p <= n ) {
 614                   while ( ( n % p ) == 0 ) {
 615                     *_Trace << p << "*" ;
 616                     n = n / p ;
 617                   }
 618                   p = NextPrime( p + 1 ) ;
 619                 }
 620                 *_Trace << n << endl ;
 621               }
 622 #endif
 623             }
 624             *_Trace << endl ;
 625           }
 626         }
 627       }
 628       *_Trace << endl ;
 629     }
 630   }
 631
 632 #if 0
 633   _Trace = new ofstream( "/cas_01/jr/linux/Syr.log" );
 634
 635   int n = 1 ;
 636   int prt ;
 637
 638 #if 0
 639   for ( prt = 1 ; prt <= 5 ; prt++ ) {
 640     *_Trace << endl << prt << ". 8n + 1 :" << endl ;
 641     little( 0 , prt , 8 , 1 ) ;
 642   }
 643
 644   for ( prt = 1 ; prt <= 5 ; prt++ ) {
 645     *_Trace << endl << prt << ". 4n + 3 :" << endl ;
 646     little( 0 , prt , 4 , 3 ) ;
 647   }
 648
 649   for ( prt = 1 ; prt <= 6 ; prt++ ) {
 650     *_Trace << endl << prt << ". 8n + 5 :" << endl ;
 651     little( 0 , prt , 8 , 5 ) ;
 652   }
 653 #endif
 654
 655   int i ;
 656   int max = 256 ;
 657   for ( i = 0 ; i <= max ; i++ ) {
 658     *_Trace << endl << endl << "8*" << i << "+5 = 4[2*" << i << "+1]+1 = " ;
 659     int ii = 2*i+1 ;
 660     while ( ((ii-1) & 7) == 4 ) {
 661       *_Trace << endl << "        8*" << ii/8 << "+5 = 4[2*" << ii/4
 662               << "+1]+1 = " ;
 663       ii = 2*(ii/4)+1 ;
 664     }
 665     *_Trace << "6*" << nP(8*i + 5) << "+" << rP(8*i + 5) << " :" << endl ;
 666     n = SeqSyr( 8*i + 5 ) ;
 667     n = n/8 ;
 668     while ( n && n > max ) {
 669       *_Trace << endl << "   ==> 8*" << n << "+5 = 4[2*" << n << "+1]+1 = " ;
 670       int ii = 2*n+1 ;
 671       while ( ((ii-1) & 7) == 4 ) {
 672         *_Trace << endl << "8*" << ii/8 << "+5 = 4[2*" << ii/4 << "+1]+1 = " ;
 673         ii = 2*(ii/4)+1 ;
 674       }
 675       *_Trace << "6*" << nP(8*n + 5) << "+" << rP(8*n + 5) << " :" << endl ;
 676       n = SeqSyr( 8*n + 5 ) ;
 677       n = n/8 ;
 678     }
 679   }
 680
 681 #if 0
 682   int x ;
 683   int rx ;
 684   int a ;
 685   int q ;
 686   int r ;
 687   int X ;
 688   int n ;
 689   int rn ;
 690   int xx ;
 691   int y ;
 692   int p ;
 693   int pow ;
 694
 695   *_Trace << "2**[x(n+1)] = (2**x-1)*S[x](n)+1" << endl << endl ;
 696
 697   *_Trace << "2**[(2x+1)(2n+rn+1)] = (2**(2x+1)-1)*S[2x+1](2n+rn)+1" << endl ;
 698   *_Trace << "                     = 2**(2x+1){(2**(2x+1)-1)*S[2x+1](2n+rn)-1}"
 699           << endl << endl ;
 700   *_Trace << "2**[2x(3n+rn+1)] = (2**2x-1)*S[2x](3n+rn)+1" << endl ;
 701   *_Trace << "                 = 2**2x{(2**2x-1)*S[2x](3n+rn-1)-1}"
 702           << endl << endl ;
 703   *_Trace << "2**(2x+1)-1 = 6S2(x-1)+1" << endl << endl ;
 704
 705   *_Trace << "2**(2(3x+1))-1 = 3(6*2*7*S6(x-1)+1) *4 :" << endl ;
 706   *_Trace << "2**(2(3x+2))-1 = 3(6*2**3*7*S6(x-1)+5) *4 :" << endl ;
 707   *_Trace << "2**(2(3x+3))-1 = 3(6*2**5*7*S6(x-1)+21)" << endl ;
 708   *_Trace << "               = 3**2(2**6*7*S6(x-1)+7)" << endl ;
 709   *_Trace << "               = 3**2*7(2**6*S6(x-1)+1)" << endl ;
 710   *_Trace << "               = 3**2*7*S6(x)" << endl ;
 711   *_Trace << "2**(2(3x))-1 = 3**2*7*S6(x-1)" << endl ;
 712   *_Trace << "2**(2(3x+4))-1 = 3(6*2*7*S6(x)+1) *4 :" << endl ;
 713   *_Trace << "2**(2(3x+5))-1 = 3(6*2**3*7*S6(x)+5)" << endl << endl ;
 714   *_Trace << "2**(2(3x+rx))-1 = 3**2*2**(2rx)*7*S6(x-1)+2**(2rx)-1" << endl
 715           << endl ;
 716
 717   *_Trace << "2**[(2x+1)(2n+rn+1)] = (6S2(x-1)+1)*S[2x+1](2n+rn)+1" << endl ;
 718   *_Trace << "2**[2(3x+rx)(3n+rn+1)] = {3**2*2**(2rx)*7*S6(x-1)+2**(2rx)-1}*"
 719           << endl << "                 S[2(3x+rx)](3n+rn)+1" << endl
 720           << endl ;
 721
 722   *_Trace << "(6x+5)y+x = 2**3*7*S6(x-1)" << endl ;
 723   *_Trace << "x = 0 ; y = 0" << endl ;
 724   *_Trace << "x = 1 ; y = 5" << endl ;
 725   *_Trace << "x = 2 ; y = 214 = 2*107 = 2*O153" << endl ;
 726   *_Trace << "x = 3 ; y = 10131 = 3*3377 = 3*O6461" << endl ;
 727   *_Trace << "x = 4 ; y = 514244 = 4*128561 = 4*O373061" << endl ;
 728   *_Trace << "x = 5 ; y = ?" << endl ;
 729   *_Trace << "x = 6 ; y = 1489853154 = 6*248308859 = 6*O1663162173" << endl
 730           << endl ;
 731
 732   for ( pow = 1 ; pow <= 31 ; pow++ ) {
 733     *_Trace << endl ;
 734     for ( x = 0 ; x <= 30 ; x++ ) {
 735       for ( n = 0 ; n <= 30 ; n++ ) {
 736         for ( rn = 0 ; rn <= 1 ; rn++ ) {
 737           if ( (2*x+1)*(2*n+rn+1) == pow && (2*x+1)*(2*n+rn+1) <= 31 ) {
 738             int val = (1 << (2*x+1)*(2*n+rn+1)) ;
 739             *_Trace << "2**(" << 2*x+1 << "*(" << 2*n+rn << "+1)) = " ;
 740             *_Trace << "2**[(2*" << x << "+1)(2*" << n << "+" << rn
 741                     << "+1)] = (6*" ;
 742             int SIG = S(2,x-1) ;
 743             int pow3 = 0 ;
 744             while ( SIG && (SIG % 3) == 0 ) {
 745               pow3 += 1 ;
 746               SIG = SIG/3 ;
 747             }
 748             *_Trace << "3**" << pow3 << "*" << SIG
 749                     << "+1) * " << fact( S(2*x+1,2*n+rn) )
 750                     << " + 1 = " << endl
 751                     << "      (6*S2(" << x-1 << ")+1)*S" << 2*x+1 << "("
 752                     << 2*n+rn << ")+1 = "
 753                     << "(6*S2(" << x << "-1)+1)*S[2*" << x << "+1](2*"
 754                     << n << "+" << rn << ")+1 = " << val << endl ;
 755             if ( val != (6*S(2,x-1)+1)*S(2*x+1,2*n+rn)+1 ) {
 756               *_Trace << val << " != " << (6*S(2,x-1)+1)*S(2*x+1,2*n+rn)+1
 757                       << endl ;
 758             }
 759           }
 760         }
 761       }
 762     }
 763     for ( x = 0 ; x <= 30 ; x++ ) {
 764       for ( rx = 0 ; rx <= 2 ; rx++ ) {
 765         for ( n = 0 ; n <= 30 ; n++ ) {
 766           for ( rn = 0 ; rn <= 2 ; rn++ ) {
 767             if ( 3*x+rx != 0 && 2*(3*x+rx)*(3*n+rn+1) == pow
 768                  && 2*(3*x+rx)*(3*n+rn+1) <= 31 ) {
 769               *_Trace << "2**[" << 2*(3*x+rx) << "*(" << 3*n+rn << "+1)] = " ;
 770               *_Trace << "2**[2(3*" << x << "+" << rx << ")(3*" << n << "+"
 771                       << rn << "+1)] = [9*" << (1 << (2*rx)) << "*7*"
 772                       << S(6,x-1) << "+" ;
 773               int DRX = (1 << (2*rx))-1 ;
 774               int pow3 = 0 ;
 775               while ( DRX && (DRX % 3) == 0 ) {
 776                 pow3 += 1 ;
 777                 DRX = DRX/3 ;
 778               }
 779               *_Trace << "3**" << pow3 << "*" << DRX << "] * " ;
 780               int SIG = S(2*(3*x+rx),3*n+rn) ;
 781               pow3 = 0 ;
 782               while ( SIG && (SIG % 3) == 0 ) {
 783                 pow3 += 1 ;
 784                 SIG = SIG/3 ;
 785               }
 786               *_Trace << "3**" << pow3 << "*" << fact( SIG )
 787                       << " + 1 = " << endl
 788                       << "   [3**2*2**" << 2*rx << "*7*S6("
 789                       << x << "-1)+" << (1 << (2*rx))-1 << "] * S"
 790                       << 2*(3*x+rx) << "(" << 3*n+rn << ")+1 = " << endl
 791                       << "      {3**2*2**(2*" << rx << ")*7*S6("
 792                       << x << "-1)+2**(2*" << rx << ")-1} * S[2(3*"
 793                       << x << "+" << rx << ")](3*" << n << "+" << rn << ")+1"
 794                       << endl ;
 795               int val = (1 << (2*(3*x+rx)*(3*n+rn+1)) ) ;
 796               if ( val != (3*3*(1<<(2*rx))*7*S(6,x-1)+(1<<(2*rx))-1)*S(2*(3*x+rx),3*n+rn)+1 ) {
 797                 *_Trace << val << " != " << (3*3*(1<<(2*rx))*7*S(6,x-1)+(1<<(2*rx))-1) * S(2*(3*x+rx),3*n+rn)+1 << endl ;
 798               }
 799             }
 800           }
 801         }
 802       }
 803     }
 804   }
 805   return 0;
 806
 807
 808   *_Trace << endl ;
 809   for ( x = 1 ; x <= 30 ; x++ ) {
 810     for ( n = 0 ; n <= 30 ; n++ ) {
 811       if ( x*(n+1) <= 30 ) {
 812         X = ((1 << x ) - 1)*S(x,n) + 1 ;
 813         if ( (1 << (x*(n+1)) ) == X ) {
 814           *_Trace << "2**[" << x << "(" << n << "+1)] = 2**[[" ;
 815           int pow2 = 0 ;
 816           int pow3 = 0 ;
 817           int p ;
 818           p = x ;
 819           while ( p && (p & 1) == 0 ) {
 820             pow2 += 1 ;
 821             p = (p >> 1 ) ;
 822           }
 823           while ( p && (p % 3) == 0 ) {
 824             pow3 += 1 ;
 825             p = p/3 ;
 826           }
 827           if ( pow2 ) {
 828             *_Trace << "2" ;
 829             if ( pow2 != 1 ) {
 830               *_Trace << "**" << pow2 ;
 831             }
 832             if ( pow3 || p != 1 ) {
 833               *_Trace << "*" ;
 834             }
 835           }
 836           if ( pow3 ) {
 837             *_Trace << "3" ;
 838             if ( pow3 != 1 ) {
 839               *_Trace << "**" << pow3 ;
 840             }
 841             if ( p != 1 ) {
 842               *_Trace << "*" ;
 843             }
 844           }
 845           if ( !((pow2 || pow3) && p == 1) ) {
 846             *_Trace << p ;
 847           }
 848           *_Trace << "][" ;
 849           p = n+1 ;
 850           pow2 = 0 ;
 851           pow3 = 0 ;
 852           while ( p && (p & 1) == 0 ) {
 853             pow2 += 1 ;
 854             p = (p >> 1 ) ;
 855           }
 856           while ( p && (p % 3) == 0 ) {
 857             pow3 += 1 ;
 858             p = p/3 ;
 859           }
 860           if ( pow2 ) {
 861             *_Trace << "2" ;
 862             if ( pow2 != 1 ) {
 863               *_Trace << "**" << pow2 ;
 864             }
 865             if ( pow3 || p != 1 ) {
 866               *_Trace << "*" ;
 867             }
 868           }
 869           if ( pow3 ) {
 870             *_Trace << "3" ;
 871             if ( pow3 != 1 ) {
 872               *_Trace << "**" << pow3 ;
 873             }
 874             if ( p != 1 ) {
 875               *_Trace << "*" ;
 876             }
 877           }
 878           if ( !((pow2 || pow3) && p == 1) ) {
 879             *_Trace << p ;
 880           }
 881           *_Trace << "]]" ;
 882           int d = (1 << x) - 1 ;
 883           pow3 = 0 ;
 884           while ( (d % 3) == 0 ) {
 885             pow3 += 1 ;
 886             d = d / 3 ;
 887           }
 888           *_Trace << " = (2**" << x << "-1)*S" ;
 889           if ( x & 1 ) {
 890             *_Trace << x << "[2*" << n/2 ;
 891             if ( n & 1 ) {
 892               *_Trace << "+1" ;
 893             }
 894             *_Trace << "]+1 = " ;
 895           }
 896           else {
 897             *_Trace << x << "[3*" << n/3 << "+" << (n % 3) << "]+1 = " ;
 898           }
 899           *_Trace << "{" ;
 900           if ( pow3 ) {
 901             *_Trace << "3" ;
 902             if ( pow3 != 1 ) {
 903               *_Trace << "**" << pow3 ;
 904             }
 905             *_Trace << "*" ;
 906           }
 907           if ( d == 1 ) {
 908           }
 909           else if ( d == 5 ) {
 910             *_Trace << "[5]" ;
 911           }
 912           else {
 913             if ( nP(d) ) {
 914               *_Trace << "[6*" << nP(d) << "+" << rP(d) << "]" ;
 915             }
 916             else {
 917               *_Trace << "[" << rP(d) << "]" ;
 918             }
 919           }
 920           *_Trace << "}*{" ;
 921           int Y = S(x,n) ;
 922           pow2 = 0 ;
 923           while ( (Y % 2) == 0 ) {
 924             pow2 += 1 ;
 925             Y = Y / 2 ;
 926           }
 927           pow3 = 0 ;
 928           while ( (Y % 3) == 0 ) {
 929             pow3 += 1 ;
 930             Y = Y / 3 ;
 931           }
 932           if ( pow2 ) {
 933             *_Trace << "2" ;
 934             if ( pow2 != 1 ) {
 935               *_Trace << "**" << pow2 ;
 936             }
 937             *_Trace << "*" ;
 938           }
 939           if ( pow3 ) {
 940             *_Trace << "3" ;
 941             if ( pow3 != 1 ) {
 942               *_Trace << "**" << pow3 ;
 943             }
 944             *_Trace << "*" ;
 945           }
 946           if ( nP(Y) && ( (nP(Y) & 1) == 0 || (nP(Y) % 3) == 0 ) ) {
 947             *_Trace << "[6*" ;
 948             int YY = nP(Y) ;
 949             pow2 = 0 ;
 950             while ( (YY % 2) == 0 ) {
 951               pow2 += 1 ;
 952               YY = YY / 2 ;
 953             }
 954             pow3 = 0 ;
 955             while ( (YY % 3) == 0 ) {
 956               pow3 += 1 ;
 957               YY = YY / 3 ;
 958             }
 959             if ( pow2 ) {
 960               *_Trace << "2" ;
 961               if ( pow2 != 1 ) {
 962                 *_Trace << "**" << pow2 ;
 963               }
 964               *_Trace << "*" ;
 965             }
 966             if ( pow3 ) {
 967               *_Trace << "3" ;
 968               if ( pow3 != 1 ) {
 969                 *_Trace << "**" << pow3 ;
 970               }
 971               *_Trace << "*" ;
 972             }
 973             *_Trace << YY << "+" << rP(Y) << "]}+1" << endl ;
 974           }
 975           else {
 976             if ( nP(Y) ) {
 977               *_Trace << "[6*" << nP(Y) << "+" << rP(Y) << "]}+1" << endl ;
 978             }
 979             else {
 980               *_Trace << "[" << rP(Y) << "]}+1" << endl ;
 981             }
 982           }
 983         }
 984         else {
 985           *_Trace << endl << (1 << (x*(n+1)) ) << " = 2**[" << x << "(" << n
 986                   << "+1)] != " << X << "(2**" << x
 987                   << "-1)*S" << x << "[" << n << "]+1" << endl ;
 988         }
 989       }
 990     }
 991   }
 992   return 0;
 993
 994   *_Trace << "S[x]( n ) = Sigma((2**x)*i) i = 0 , n" << endl ;
 995   *_Trace << "S[x]( aq + r ) = X*S[a*x]( q - 1 ) + Sx( r )" << endl ;
 996   for ( x = 2 ; x <= 10 ; x++ ) {
 997     *_Trace << endl << endl ;
 998     for ( a = 2 ; a <= 12 ; a++ ) {
 999       *_Trace << endl ;
1000       for ( q = 1 ; q <= 12 ; q++ ) {
1001         for ( r = 0 ; r < a ; r++ ) {
1002           if ( a*q + r <= 24 ) {
1003             int i = S( x , a*q+r ) - S( x , r ) ;
1004             int X = i / S( a*x , q-1 ) ;
1005             int p2 = 0 ;
1006             if ( X*S( a*x , q-1 ) == i && (i % S( a*x , q-1 ) ) == 0 ) {
1007               while ( ( X & 1 ) == 0 ) {
1008                 p2 += 1 ;
1009                 X = ( X >> 1 ) ;
1010               }
1011               int p3 = 0 ;
1012               while ( X > 0 && ( X % 3 ) == 0 ) {
1013                 p3 += 1 ;
1014                 X = ( X / 3 ) ;
1015               }
1016               *_Trace << "S( " << x << " , " << a << "*" << q << " + " << r
1017                       << " ) = 2**" << p2 << "*3**" << p3 << "*"<< X
1018                       << "*S[" << a << "*" << x << "]( " << q << "-1 ) + S( "
1019                       << x << " , " << r << " )" << endl ;
1020             }
1021           }
1022         }
1023       }
1024     }
1025   }
1026   return 0;
1027
1028   int PT ;
1029   int PP ;
1030   int PPNext ;
1031   for ( PP = 1 ; PP <= 5 ; PP+=4 ) {
1032     for ( PPNext = 1 ; PPNext <= 5 ; PPNext+=4 ) {
1033       for ( PT = 0 ; PT <= 2 ; PT++ ) {
1034         if ( PT == 0 && PP == 5 ) {
1035           *_Trace << endl << "T1a 2**(6n1+2)(3n0) + 2**(6n1+1) + S2(3n1+0) --> 6(3n0+0)+5"
1036                   << endl ;
1037           Compute( "T1a" , T1a , PT , PP , PPNext ) ;
1038         }
1039         if ( PT == 1 && PP == 1 ) {
1040           *_Trace << endl << "T2b 2**(6n1+3)(3n0+1) + S2(3n1+0) --> 6(3n0+1)+1"
1041                   << endl ;
1042           Compute( "T2b" , T2b , PT , PP , PPNext ) ;
1043         }
1044         if ( PT == 0 && PP == 1 ) {
1045           *_Trace << endl << "T2a 2**(6n1+7)(3n0) + S2(3n1+2) --> 6(3n0+0)+1"
1046                   << endl ;
1047           Compute( "T2a" , T2a , PT , PP , PPNext ) ;
1048         }
1049         if ( PT == 2 && PP == 5 ) {
1050           *_Trace << endl << "T1c 2**(6n1+4)(3n0+2) + 2**(6n1+3) + S2(3n1+1) --> 6(3n0+2)+5"
1051                   << endl ;
1052           Compute( "T1c" , T1c , PT , PP , PPNext ) ;
1053         }
1054         if ( PT == 2 && PP == 1 ) {
1055           *_Trace << endl << "T2c 2**(6n1+5)(3n0+2) + S2(3n1+1) --> 6(3n0+2)+1"
1056                   << endl ;
1057           Compute( "T2c" , T2c , PT , PP , PPNext ) ;
1058         }
1059         if ( PT == 1 && PP == 5 ) {
1060           *_Trace << endl << "T1b 2**(6n1+6)(3n0+1) + 2**(6n1+5) + S2(3n1+2) --> 6(3n0+1)+5"
1061                   << endl ;
1062           Compute( "T1b" , T1b , PT , PP , PPNext ) ;
1063         }
1064
1065 //      for ( PPNext = 1 ; PPNext <= 5 ; PPNext+=4 ) {
1066         *_Trace << endl << "P2a1 (1<<(6*n1+3))*(3*n0+0) + S2(3*n1+0) --> 6(3n0+0)+1"
1067                 << endl ;
1068         Compute( "P2a1" , P2a1 , PT , PP , PPNext ) ;
1069
1070         *_Trace << endl << "P2a2 (1<<(6*n1+5))*(3*n0+0) + S2(3*n1+1) --> 6(3n0+0)+1"
1071                 << endl ;
1072         Compute( "P2a2" , P2a2 , PT , PP , PPNext ) ;
1073
1074         *_Trace << endl << "P1b1 (1<<(6*n1+2))*(3*n0+1) + (1<<(6*n1+1)) + S2(3*n1+0) --> 6(3n0+1)+5"
1075                 << endl ;
1076         Compute( "P1b1" , P1b1 , PT , PP , PPNext ) ;
1077
1078         *_Trace << endl << "P1c1 (1<<(6*n1+2))*(3*n0+2) + (1<<(6*n1+1)) + S2(3*n1+0) --> 6(3n0+2)+5"
1079                 << endl ;
1080         Compute( "P1c1" , P1c1 , PT , PP , PPNext ) ;
1081
1082         *_Trace << endl << "P1a1 (1<<(6*n1+4))*(3*n0+0) + (1<<(6*n1+3)) + S2(3*n1+1) --> 6(3n0+0)+5"
1083                 << endl ;
1084         Compute( "P1a1" , P1a1 , PT , PP , PPNext ) ;
1085
1086         *_Trace << endl << "P2c1 (1<<(6*n1+3))*(3*n0+2) + S2(3*n1+0) --> 6(3n0+2)+1"
1087                 << endl ;
1088         Compute( "P2c1" , P2c1 , PT , PP , PPNext ) ;
1089
1090         *_Trace << endl << "P1b2 (1<<(6*n1+4))*(3*n0+1) + (1<<(6*n1+3)) + S2(3*n1+1) --> 6(3n0+1)+5"
1091                 << endl ;
1092         Compute( "P1b2" , P1b2 , PT , PP , PPNext ) ;
1093
1094         *_Trace << endl << "P2b1 (1<<(6*n1+5))*(3*n0+1) + S2(3*n1+1) --> 6(3n0+1)+1"
1095                 << endl ;
1096         Compute( "P2b1" , P2b1 , PT , PP , PPNext ) ;
1097
1098         *_Trace << endl << "P1a2 (1<<(6*n1+6))*(3*n0+0) + (1<<(6*n1+5)) + S2(3*n1+2) --> 6(3n0)+5"
1099                 << endl ;
1100         Compute( "P1a2" , P1a2 , PT , PP , PPNext ) ;
1101
1102         *_Trace << endl << "P2b2 (1<<(6*n1+7))*(3*n0+1) + S2(3*n1+2) --> 6(3n0+1)+1"
1103                 << endl ;
1104         Compute( "P2b2" , P2b2 , PT , PP , PPNext ) ;
1105
1106         *_Trace << endl << "P1c2 (1<<(6*n1+6))*(3*n0+2) + (1<<(6*n1+5)) + S2(3*n1+2) --> 6(3n0+2)+5"
1107                 << endl ;
1108         Compute( "P1c2" , P1c2 , PT , PP , PPNext ) ;
1109
1110         *_Trace << endl << "P2c2 (1<<(6*n1+7))*(3*n0+2) + S2(3*n1+2) --> 6(3n0+2)+1"
1111                 << endl ;
1112         Compute( "P2c2" , P2c2 , PT , PP , PPNext ) ;
1113       }
1114
1115     }
1116   }
1117 #endif
1118
1119 #endif
1120
1121   return 0;
1122 }
1123