src/SyrComponent/SyrComponent_CheckOfUndefined.cxx

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