Jets.s 206 KB
Newer Older
Hynek Baran's avatar
Hynek Baran committed
1
# J e t s 
Hynek Baran's avatar
Hynek Baran committed
2
#
Hynek Baran's avatar
Hynek Baran committed
3
# Differential calculus on jet spaces and diffieties
Hynek Baran's avatar
Hynek Baran committed
4 5
#

Hynek Baran's avatar
Hynek Baran committed
6 7 8
#
# See http://jets.math.slu.cz/ for further details
#
Hynek Baran's avatar
Hynek Baran committed
9

Hynek Baran's avatar
Hynek Baran committed
10 11 12 13 14 15
#
# A u t h o r s
#
# (c) 1999-2004 by Michal Marvan
# (c) 2005-2010 by Hynek Baran and Michal Marvan
#
Hynek Baran's avatar
Hynek Baran committed
16 17 18 19 20 21 22 23
#
#  L i c e n s e
#
# JETS is a free software distributed under the terms of the GNU General 
# Public License <http://www.gnu.org/copyleft/gpl.html> as published by 
# the Free Software Foundation.
#
# In particular, JETS comes with ABSOLUTELY NO WARRANTY.
Hynek Baran's avatar
Hynek Baran committed
24 25 26 27 28
#
# 
# A c k n o w l e d g m e n t s
#
# The first-named author was supported by GACR under project 201/07/P224
Hynek Baran's avatar
Hynek Baran committed
29
# Developement of Jets was gratefully supported by project CZ.1.07/2.3.00/20.0002. 
Hynek Baran's avatar
Hynek Baran committed
30

Hynek Baran's avatar
Hynek Baran committed
31 32

#
Hynek Baran's avatar
Hynek Baran committed
33 34
# History of changes
#
Hynek Baran's avatar
Hynek Baran committed
35 36
# v. 5.7  rel. 8 Dec 2010
# * Based on Jets code v. 5.6 as of 20 Jan 2010
Hynek Baran's avatar
Hynek Baran committed
37
#
Hynek Baran's avatar
Hynek Baran committed
38
# v. 5.71 rel. 4 Jan 2011
Hynek Baran's avatar
Hynek Baran committed
39
# * `divideout/unks/1`(): approved, better handling of polynomials
Hynek Baran's avatar
Hynek Baran committed
40
#
Hynek Baran's avatar
Hynek Baran committed
41 42 43 44
# v 5.72
# * `index/traceless` matrix index function added
# * store(): formatting of output of assignments changed (each on separate line now)
# * new `jet/gen/all` and `count/gen/all` functions 
Hynek Baran's avatar
Hynek Baran committed
45 46 47
#
# v 5.73 alpha
# * first attempt to fix a critical bug in `cc/new` (some cc's was omitted by run)
Hynek Baran's avatar
Hynek Baran committed
48 49 50 51 52 53 54 55 56 57
#
# v 5.73 beta
# * `divideout/unks`, `type/divideout`: fixed (it was ok in very old jets but broken later)
#   and little improved ((exp(F)+1) is divided out now)
# * Varordering: new ordering revdeg implemented;
#   [degree, revdeg] is equivalent to the graded reverse lexicographic order (tdeg or grevlex)
# * `resolve/normalizer` added
# * Reporter implemented (quite primitive, resolve only)  
#
# v 5.74 alpha
Hynek Baran's avatar
Hynek Baran committed
58 59
# * resolve completly rewritten
# * started use of parallel MaP version of map (in resolve and derive) but no still not supported by Maple 16
Hynek Baran's avatar
Hynek Baran committed
60 61 62
# * derive redesigned
# * new functions Varl, VarL, LVar
# * minor evalTD bug fixed (wrong result on matrices) 
Hynek Baran's avatar
Hynek Baran committed
63 64 65 66 67 68 69 70 71 72 73 74 75
# * nonzero() bug (missing factor@) fixed (probably big issiue)
# * minor bug in Varordering fixed (if none fiber vars)
# * `resolve/nonresrat` limit introduced (price overlap of resolvable epr. over nonresolvable enforces FAIL of resolve)
# * minor improvements in reporting of assembled cc's
# * TestUnkOrd() introduced
# * sizefactor
# * `resolve/lin/price` back in the simplest form
# * `divideout/normalizer` added
# * Maple 14 compatibilty
#
#
# v 5.75 
# * old resolve is default
Hynek Baran's avatar
Hynek Baran committed
76 77
#
# v. 5.76
Hynek Baran's avatar
Hynek Baran committed
78 79 80 81
# * apply fixed (last ERROR replaced with unevaluated apply)
#
# v 5.77
# * `derive/1` : treating of polynoms and sums disabled (a bug found)
Hynek Baran's avatar
Hynek Baran committed
82 83 84 85 86 87 88
#
# v 5.78
# package JetMachine introduced
# JetMachine:-Consequences introduced
# * `jet/gen/all` fixed
# * `convert/crack` introduced
# * reenabled Vars to accept sequential arguments
89 90 91 92 93 94 95 96 97
# * uncomment the line bellow to enable the new resolve implementation
# jets_new_resolve_enable := true;
#
#
# v 5.79
# * fixed missing local declaration in size() 
# * `resolve/nonresrat` reimplemented, resolve() FAILs when size condition is triggered
# * minor reporting improvements in `run/l`()
# * JetsProfiler introduced
Hynek Baran's avatar
Hynek Baran committed
98

99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
# * BasisExtractor function to symmetries/laws basis extraction introduced
# * `divideout/unks` improved (treating exponents)
# * `resolve/lin` unifyed
# * reporting in `resolve/lin` (lin resolvable/nonres. / nonlin)
# * CC:-markFF assertion commented out as a first aid, this issuie is to be resolved carefully
# * `TestUnkOrd/vars` weakened to report strong unknown ordering suspicients only
# * `run/l/extraders` introduced allowing to derive also subsets of cc in run (but none implementations besides empty)
# * fixed `AmICons/dependence` bug
# * fixed revdeg
# * putsize back
#
#
# v 5.80 beta
# * improved derive (noderives changed)
#
# 
# v 5.81
# * version mismatch cleanup
Hynek Baran's avatar
Hynek Baran committed
117 118 119 120 121 122
# * apply -> vfapply M.M. 20 Feb. 2014
# * run() HB bugfix #if Bytes() > Blimit then reduce() fi
# * `divideout/unks`() bugfix: nonnegative numeric exponents only are divided out
# * jetorder() uses max()
# * `size/*/<`() at lest one result is always returned (low ressize and putsize settings do not cause error)
# * derive reporting slightly changed
Hynek Baran's avatar
Hynek Baran committed
123 124 125
#
#
# v 5.82
Hynek Baran's avatar
Hynek Baran committed
126
# * coordinates(), newfibre() and `jet/aliases`() has new option separator:="_" ()
Hynek Baran's avatar
Hynek Baran committed
127
# * using coordinates(..., separator="__"), jets are prettyprinted as subscripts (prior Maple 17)
Hynek Baran's avatar
Hynek Baran committed
128
# * `jet/aliases/mainseparator` is a global variant of the above option
Hynek Baran's avatar
Hynek Baran committed
129
# * TestUnkOrd() functionality limited to simple dependecy test (memory leak in second part of the test found)
Hynek Baran's avatar
Hynek Baran committed
130
# * BasisExtractor: unknown U collected in the result
131 132 133 134
#
#
# v 5.83
# * removed reduce() when store() is called (inside `store/pds`())
135 136 137
#
#
# v 5.84
Hynek Baran's avatar
Hynek Baran committed
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
# * fixed `resolve/lin` new implementation reportfail bug (jets_new_resolve_enable only)
# * linderive() introduced (but not used yet) 
#
#
# v 5.85
# * clear() rewritten:
#   + clear(pds) is also clearing cc's (unless `clear/pds/ccclearing`:='suppress' is assigned)
#   + clear(puts) and clear(assignments) introduced
#   + clear(..., output=eq) output format modifier
# * `puts/assignments`() and `put/name/tab` defined (as an generalisation of `unks/assignments`() resp. `unk/tab`)
# * store() is saving:
#   + all put symbols (form `put/name/tab`) instead of unknowns only (`unk/tab`)
#   + unknowns order (restored from `unk/tab`)
#   + varordering and Varordering
#   + terminal issue fixed  
#
#
# v 5.86
# * version mismatch (5.85) hopefully fixed 
# * `put/sizelimit` and `size/1/DD` introduced
# * run() refactored (cc's are derived)
# * `put/1` distinguishes redundand and contradictory attempt of putting already assigned expressions
Hynek Baran's avatar
Hynek Baran committed
160 161 162 163 164 165 166 167
#
# v 5.87
# * not deriving cc's by default (`run/l/extraders` used in slightly different mode)
# * `run/1` reporting cleanup
# * Report Reportf macros rewritten, rt, rb no more needed
# * `size/1/DD` abandoned for a hidden bug
# * `put/limit/length`, `put/limit/size` instead `put/sizelimit` 
# * JetMachine[Consequences]:-AmICons: Using `AmICons/ignore`() for nonzero(), Varordering() and unknowns()
Hynek Baran's avatar
Hynek Baran committed
168
#
Hynek Baran's avatar
Hynek Baran committed
169 170 171 172 173
# v 5.88
# * fixed `resolve/lin` new implementation reportfail bug (jets_new_resolve_enable only)
# * linderive() introduced (but not used yet) 
# * `unks/TD` uses forceError=true in `vars/1` calls
# * JetMachine[Consequences] testing script changed (currentdir)
Hynek Baran's avatar
Hynek Baran committed
174 175
# * `clear/assignments` is clearing `put/name/tab` properly# current
# * JetMachine[Consequences] testing script changed (currentdir)
Hynek Baran's avatar
Hynek Baran committed
176
# * `clear/assignments` is clearing `put/name/tab` properly
Hynek Baran's avatar
Hynek Baran committed
177 178 179
#
# v 5.89
# * run: linear eqs are ALL passed to resolve
Hynek Baran's avatar
Hynek Baran committed
180 181 182 183
#
# v 5.90
# * jet: option remember bug (multiple eqs problem as reported by JSK) fixed by introducing `eqn/table`
# * lengthselect introduced
Hynek Baran's avatar
Hynek Baran committed
184 185 186 187 188 189 190
#
# v. 5.91
# * Vars has option cache`
# * resolve FAILs collecting introduced
# * NewIntSeq() introduced
# * Report macro changed to inlined function
# * New resolve implementation (jets_new_resolve_enable) removed from the Jets.s source file  to the Jets.newresolve.s file
Hynek Baran's avatar
Hynek Baran committed
191
#
Hynek Baran's avatar
Hynek Baran committed
192
# v 5.93
Hynek Baran's avatar
Hynek Baran committed
193
# * testing new nonlinear resolve implementation in file (cat(jets_new_resolve_enable, "/Jets.newresolve.s"))
194
# * LVar returns NULL on expressions without Vars
195
# * numerous bug fixes in the computation of coverings (MM 23.4.2019)
Hynek Baran's avatar
Hynek Baran committed
196

Hynek Baran's avatar
Hynek Baran committed
197 198 199 200 201 202 203
###########################################################################################
###########################################################################################
# Initialization and configuration
###########################################################################################
###########################################################################################

interface(screenwidth=120):
Hynek Baran's avatar
Hynek Baran committed
204
printf("Jets 5.93  as of May 2, 2019 (%s)\n",
Hynek Baran's avatar
Hynek Baran committed
205
      StringTools[FormatTime]("%Y-%m-%d %T", timestamp=FileTools[Status](__FILE__)[5]));
Hynek Baran's avatar
Hynek Baran committed
206 207 208 209 210 211 212 213 214 215 216

#
# Source code configuration, options and parameters
#

### optional features
# Some parts of Jets source code is disbbled by default,
# since it is not proprely finished or tested yet.
# Use assignments bellow (before reading Jets.s in) to enable such additional code.

#`Jets/opts`["Optionals"]["Multiord"] := true:
Hynek Baran's avatar
Hynek Baran committed
217
#jets_new_resolve_enable := ".";
Hynek Baran's avatar
Hynek Baran committed
218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265

#
# Debugging
#

### reporting
#
# Reporting about the progression of calculation is available independently of debug mode bellow.
# Syntax: reporting(derive = 0, resolve = 0, run = 0, cc = 0, pd = 0); 

### debug mode switch
#
# Some debugging functionality (as fuction arguments type checks, assumptions or logging) 
# is disabled by default but may be turned on by the macro bellow. 
# A lot of computing time is saved when debug mode is disabled, 
# but the AssertLevel() function and some more debugging functionality is not available.
# Turning debugging mode on may dramatically harm the time and memory usage.
#
# In the case of any problem (as an error raised or wrong result returned),
# please turn the debugging mode ON by removing the very first character # in the line bellow:
#
#$define jets_mode_debugging # <<<- TURN ON DEBUG MODE HERE by UNcommenting this line #  <- DEBUG ON/OFF


### debugging code
$ifdef jets_mode_debugging
# the case when debugging is ON:
print(`debugging mode enabled`);
print(kernelopts(version));
kernelopts(opaquemodules=false): 
$define ATC(a,t) a::t # Arguments Type Control enabled
$define DOE(a) a      # Debug Only Expressions enabled
$else
# the case when debugging is OFF:
$define ATC(a,t) a # Arguments Type Control disabled
$define DOE(a)     # Debug Only Expressions ignored
$endif

#
# Initialization tests
#

### single load test
# make sure that this source file is not already loaded:
if assigned(cat(`jets/read/flag/`,__FILE__)) then error "source file %1 already loaded.", __FILE__; fi;
if assigned(cat(`jets/read/flag/`,__FILE__)) then `quit`(255); fi; # force immediate exit from read() statement
assign(cat(`jets/read/flag/`,__FILE__), true):

Hynek Baran's avatar
Hynek Baran committed
266
### Reporting macros
267 268 269 270 271

ProcNameN := proc(n) option inline; debugopts('callstack')[1+1+3*n] end: # The list will contain 3n+4 elements. The first element is the name DEBUGSTACK. The remaining elements are in groups of 3.
ProcBaseNameN  := proc(n) option inline; StringTools:-Split(convert((ProcNameN(n)), string), "/")[1] end:
#ProcBaseName   := proc() option inline; remove(`=`, map(ProcBaseNameN, [$ 1..((nops(debugopts('callstack'))-1)/3)]), "unknown")[1] end: ### VERY slow
ProcBaseName   := proc() option inline; ProcBaseNameN(0) end:
Hynek Baran's avatar
Hynek Baran committed
272 273
ProcBaseSymbol := proc() option inline; convert(ProcBaseName(), symbol) end:
ProcBaseSYMBOL := proc() option inline; convert(StringTools:-UpperCase(ProcBaseName()), symbol) end:
Hynek Baran's avatar
Hynek Baran committed
274 275 276 277

Report := proc(l, m)
  option inline;   
  `if` (`report/tab`[ProcBaseSymbol()] > l, # l<=0 or ... for some magic reason wont work
Hynek Baran's avatar
Hynek Baran committed
278
        report(cat(ProcBaseSYMBOL(), '`:`'), [cat(ProcNameN(0),'`[`',l,'`]`',":"), `if`(type(m,list),op(m),m)]),
Hynek Baran's avatar
Hynek Baran committed
279 280 281 282 283 284
        NULL);
end:

Reportf := proc(l, m)
  option inline;  
  `if`(`report/tab`[ProcBaseSymbol()] > l,
Hynek Baran's avatar
Hynek Baran committed
285
       report(cat(ProcBaseSYMBOL(), '`:`'), sprintf(cat("%a[%a]: ", op(1,m)), ProcNameN(0), l, op(2..-1, m))),
Hynek Baran's avatar
Hynek Baran committed
286 287 288 289
       NULL);
end:


Hynek Baran's avatar
Hynek Baran committed
290 291 292

##############################
#MaP := evalf(Threads[Map]);
Hynek Baran's avatar
Hynek Baran committed
293
MaP := evalf(map):
Hynek Baran's avatar
Hynek Baran committed
294

Hynek Baran's avatar
Hynek Baran committed
295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318
#
# old Maple version backward compatibility
#
# the routines defined in newer versions only are called as compat[routine]
# in the code bellow are these routines 
# - assigned to this auxiliary table (if available in the current version)
# - implemented by Hynek for backward compatibility 

### Maple 14 is probably the minimum
if not(assigned(DifferentialAlgebra)) then # DifferentialAlgebra available since Maple 14
  lprint(kernelopts(version));
  error "Maple older than 14 detected, Jets cannot be used.";
fi;
if not(assigned(DifferentialAlgebra)) then
  `quit`(255);
fi;
### Maple 15 features:
if not(assigned(numelems)) then # numelems available since Maple 15
  printf("Warning: Maple of version older than 15 detected. Please upgrade to Maple 15.\n%s\n", kernelopts(version));
  compat[Record[packed]]:= Record: # Record[packed] not available, use Record instead
else
  # we have Maple 15 (or newer)
  compat[Record[packed]] := Record[packed]:
fi:
Hynek Baran's avatar
Hynek Baran committed
319

Hynek Baran's avatar
Hynek Baran committed
320 321 322 323 324 325 326 327 328 329 330 331 332 333 334
#
# Aux macros
#

# create generator of integer sequence
NewIntSeq := proc(initvalue::integer:=0) 
  local i := initvalue;
  proc({`set`::{integer,identical(NULL)}:=NULL, `get`::truefalse:=false})
    if type(`set`, integer) then i := `set`;
    elif `get` then return i;
    else i := i+1
    fi;
  end
end:

Hynek Baran's avatar
Hynek Baran committed
335 336 337 338 339 340 341
###########################################################################################
###########################################################################################
# Jets - Main procedures
###########################################################################################
###########################################################################################

# This is a set of main routines based on the original Jets.s by MM
Hynek Baran's avatar
Hynek Baran committed
342 343 344 345 346 347 348 349 350 351 352 353 354

#
#  D e c l a r a t i o n s
#

# A procedure whose purpose is to change settings of global variables 
# is called a declaration.
# If named in plural, it resets the previous declaration.
# If named in singular, it has a cummulative effect.

#
#  A u x i l i a r y   p r o c e d u r e s
#
Hynek Baran's avatar
Hynek Baran committed
355 356 357

Existing := proc(c,a) assigned(c||a) end:
Call := proc(c,a) c||a end: # || replaces . since 6.0
Hynek Baran's avatar
Hynek Baran committed
358
    
Hynek Baran's avatar
Hynek Baran committed
359 360
# sorting by size

Hynek Baran's avatar
Hynek Baran committed
361
size := proc() # HB generalized size of any number of arguments
362
  local a;
Hynek Baran's avatar
Hynek Baran committed
363
  #try 
Hynek Baran's avatar
Hynek Baran committed
364
    add(`size/1`(a), a in [args]);
Hynek Baran's avatar
Hynek Baran committed
365 366 367 368 369
  #catch:
  #  printf("\nsize: failed and returned infinity: %q\n", 
  #         StringTools[FormatMessage](lastexception[2..-1]));
  #  infinity;
  #end:
Hynek Baran's avatar
Hynek Baran committed
370 371
end:

Hynek Baran's avatar
Hynek Baran committed
372
sizefactor := 1000000000:
Hynek Baran's avatar
Hynek Baran committed
373

Hynek Baran's avatar
Hynek Baran committed
374
`size/1` := proc(a) # MM size evaluator
Hynek Baran's avatar
Hynek Baran committed
375
  global sizefactor;
Hynek Baran's avatar
Hynek Baran committed
376
  local b;
Hynek Baran's avatar
Hynek Baran committed
377 378
  #b := eval(eval(a,  TD=`size/1/DD`), pd=`size/1/DD`);
  evalf(nops(Vars(a)) + length(a)/sizefactor);
Hynek Baran's avatar
Hynek Baran committed
379 380
end:

Hynek Baran's avatar
Hynek Baran committed
381 382 383 384 385 386 387 388
## there was a hidden bug, consider
## 1/(pd(F,u)-pd(F,v)) --> 1/(DD(F,d)-DD(F,d)) = 1/0
##
#`size/1/DD` := proc(f, p)
#  # do not count too much of TD(f,p), pd(f,p) derivating monomial p
#	local d; # dummy 
#	'DD'(f, cat(d $ degree(p))) # slight penalty for higher derivatives
#end:
Hynek Baran's avatar
Hynek Baran committed
389

Hynek Baran's avatar
Hynek Baran committed
390 391 392 393 394 395 396
#sizesort := proc(ql,pr)  # ql = list of expressions, pr = sizeing proc
#  local qls;
#  qls := map(proc(q,pr) [q,pr(q)] end, ql, pr);
#  qls := sort(qls, proc(a,b) evalb(op(2,a) < op(2,b)) end);
#  map (proc(a) op(1,a) end, qls)
#end:
sizesort := proc(L::list, sizef)
Hynek Baran's avatar
Hynek Baran committed
397 398
  local l;
  return map(attributes,sort([seq](setattribute(SFloat(sizef(l),0),l), l=L),`<`));
Hynek Baran's avatar
Hynek Baran committed
399 400
end proc:
# see http://www.mapleprimes.com/blog/joe-riel/sorting-with-attributes
Hynek Baran's avatar
Hynek Baran committed
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420

# reducing products

reduceprod :=  proc(a)   
  if type (a,`^`) then
    if type (op(2,a), posint) then procname(op(1,a))
    elif type (op(2,a), negint) then 1
    else a
    fi
  elif type (a,`*`) then convert(map(procname,[op(a)]),`*`)
  else a
  fi
end:

#
#  V a r i a b l e s 
#

# Declaration of coordinates (only plural form):

Hynek Baran's avatar
Hynek Baran committed
421
coordinates := proc (blist,flist, MaxAliasDegree::integer:=0, {separator::string:=`jet/aliases/mainseparator`})
Hynek Baran's avatar
Hynek Baran committed
422 423
  global `b/var/list`,`b/var/s`,`b/dim`,`b/<</list`,
    `f/var/list`,`f/var/s`,`f/dim`,`f/<</list`;
Hynek Baran's avatar
Hynek Baran committed
424
  local res;
Hynek Baran's avatar
Hynek Baran committed
425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446
  if nargs = 0 then
    ERROR(`arguments should be:\n
[list of base variables], [list of fibre variables], optional maxorder`) 
  fi;
  if not type([op(blist)], list(name)) then
    ERROR(`base coordinates must be unassigned names`)
  fi;
  if not type([op(flist)], list(name)) then
    ERROR(`fibre coordinates must be unassigned names`)
  fi;
  if nops([op(blist),op(flist)]) <> nops({op(blist),op(flist)}) then
    ERROR(`coordinates must be mutually different`)
  fi;
  `b/var/list` := blist;
  `b/var/s` := {op(blist)};
  `b/dim` := nops(`b/var/s`);
  `b/<</list` := `b/var/list`:
  `f/var/list` := flist;
  `f/var/s` := {op(flist)};
  `f/dim` := nops(`f/var/s`);
  `f/<</list` := `f/var/list`;
  noderives():
447
  doderives();
Hynek Baran's avatar
Hynek Baran committed
448 449
  refresh ();
# Result:
Hynek Baran's avatar
Hynek Baran committed
450 451 452 453 454 455
  if MaxAliasDegree > 0 then 
    res := `jet/aliases`(`f/var/list`, MaxAliasDegree, ':-separator'=separator);
    if separator <> "_" then # for backward compatibility, allways alias u_x
       `jet/aliases`(`f/var/list`, MaxAliasDegree, ':-separator'="_") 
    fi; 
    res;
Hynek Baran's avatar
Hynek Baran committed
456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522
  else `b/var/list`, `f/var/list`
  fi
end:

`b/var/list` := []: `b/var/s` := {}: `b/dim` := 0:
`f/var/list` := []: `f/var/s` := {}: `f/dim` := 0:

# Type checking procedures:

`type/b/var` := proc (x) member (x, `b/var/s`) end:
`type/f/var` := proc (f) member (f, `f/var/s`) end:

# Jet variables are unassigned calls to 'jet'.

`type/j/var` := specfunc(anything,jet):

# Variables are of type 'var'.

`type/var` := proc(x)
  if type (x,'name') then type (x,{`b/var`, `f/var`})
  else type (x,specfunc(anything,jet))
  fi
end:

#
#  C o u n t s 
#

# A count is a product of nonnegative integer powers, but 1 is not a count.

`type/count` := proc (x) 
  local x1;  
  if type (x,`^`) then type (op(2,x), posint) 
  elif type (x,`*`) then
    for x1 in x do if not type (x1,'count') then RETURN(false) fi od;
    true
  elif x = 1 then false
  else true
  fi
end:

# A count is proper when it is a product or a power

`type/count/proper` := {`*`,`^`}:

# A count is of first order when it is neither a product nor a power

`type/count/1order` := proc(x) not type(x,`count/proper`) end:

# A b/count is a product of nonnegative integer powers of base variables,
# but 1 is not a b/count.

`type/b/count` := proc (x) 
  local x1;  
  if type (x,`^`) then type (op(1,x),`b/var`) and type (op(2,x), posint) 
  elif type (x,`*`) then
    for x1 in x do if not type (x1,`b/count`) then RETURN(false) fi od;
    true
  else type (x,`b/var`)
  fi
end:

# Dealing with counts: 
# `count/f`(x) := the first ``factor'' of the count x;
# `count/r`(x) := x/`count/f`(x);
# `count/length`(x) := the sum of the powers in x;
# `count/sgn`(x) := the parity of `count/length`x;
Hynek Baran's avatar
Hynek Baran committed
523
# `transform/count`(x) :=  the product of non-negative powers in x # = numerator 
Hynek Baran's avatar
Hynek Baran committed
524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549
# (transforms x into 1 or count).

`count/f` := proc (x)
  if type (x,`^`) then op(1,x)
  elif type (x,`*`) then `count/f`(op(1,x))
  else x
  fi
end:

`count/r` := proc (x)
  if type (x,`^`) then x/op(1,x)
  elif type (x,`*`) then x/`count/f`(op(1,x))
  else NULL
  fi
end:

`count/length` := proc (x)
  if x = 1 then 0 
  elif type (x,`^`) then op(2,x)
  elif type (x,`*`) then convert (map (`count/length`, [op(x)]), `+`)
  else 1
  fi
end:

`count/sgn` := proc (x) 1 - 2*(`count/length`(x) mod 2) end:

Hynek Baran's avatar
Hynek Baran committed
550 551 552 553 554 555 556 557 558
#`transform/count` := proc(x)
#  local aux;
#  if type (x,`^`) then if op(2,x) < 0 then 1 else x fi
#  elif type (x,`*`) then aux := map(`transform/count`, x);
#  else x
#  fi
#end:

`transform/count` := op(numer):
Hynek Baran's avatar
Hynek Baran committed
559

Hynek Baran's avatar
Hynek Baran committed
560
`count/gen/all` := proc(b::list, m::{posint,0}) # HB
Hynek Baran's avatar
Hynek Baran committed
561 562 563 564 565 566 567 568 569 570 571 572 573
  description "Generates all counts in variables b up to m-th order";
  `count/gen/all/1`(b,m)[2..-1]; # 1 is not a count
end:

`count/gen/all/1` := proc(b::list, m::{posint,0}) # HB
  if nops(b) = 0 or m<=0 then 
    1
  else
    seq(seq(i*j, j in [procname(b[2..-1], m-`count/length`(i))]), i in seq(b[1]^i, i=0..m));
  fi
end:


Hynek Baran's avatar
Hynek Baran committed
574 575 576 577 578 579 580 581 582
# Jet order of a variable q

varorder := proc(q)
  if type(q,'name') then 0
  else `count/length`(op(2,q))
  fi
end:

jetorder := proc(f)
Hynek Baran's avatar
Hynek Baran committed
583
  max(jetorders(f)) # HB
Hynek Baran's avatar
Hynek Baran committed
584 585
end:

Hynek Baran's avatar
Hynek Baran committed
586
jetorders := proc(f) # HB
Hynek Baran's avatar
Hynek Baran committed
587
  map(varorder,`vars/1`(f))
Hynek Baran's avatar
Hynek Baran committed
588 589
end:

Hynek Baran's avatar
Hynek Baran committed
590 591 592 593 594 595 596 597 598 599 600 601 602
# If you assign
# `jet/aliases/mainseparator` := "__" ;
# BEFORE reading-in jets.s
# jets will be pretyprinted using subscripts, i. e. aliased using "__" instead of "_"
# jet(u,x) = u__x
# and old fashion jet aliases will be still available for backward compatibility
# jet(u,x)=u__x=u_x
if not(assigned(`jet/aliases/mainseparator`)) then `jet/aliases/mainseparator` := "_" fi:

# Creating aliases for jet variables. Arguments are:
# f1,...,fn,r, where r = jet order and 
# f's are fibre variables (optional)

Hynek Baran's avatar
Hynek Baran committed
603 604 605 606
# Creating aliases for jet variables. Arguments are:
# f1,...,fn,r, where r = jet order and 
# f's are fibre variables (optional)

Hynek Baran's avatar
Hynek Baran committed
607
`jet/aliases` := proc ({separator::string:="_"})
Hynek Baran's avatar
Hynek Baran committed
608 609
  global `b/var/list`;
  local aux,flist,r;
Hynek Baran's avatar
Hynek Baran committed
610 611 612
  r := _rest[-1];
  if _nrest = 1 then flist := `f/var/list` 
  else flist := _rest[1..-2] 
Hynek Baran's avatar
Hynek Baran committed
613 614 615 616
  fi;
  if not type (`b/var/list`,'list'(symbol)) then
    ERROR (`base variables must be symbols`) 
  fi;
Hynek Baran's avatar
Hynek Baran committed
617 618 619
  #aux := expand((1 + convert(`b/var/s`,`+`))^r - 1);  
  #aux := map(proc(a) a/coeffs(a) end, [op(sort(aux, `b/var/list`))]);
  aux := sort([`count/gen/all`(`b/var/list`, r)]); # HB
Hynek Baran's avatar
Hynek Baran committed
620
  alias(op(map(proc(u,aux) if type(u,symbol) then
Hynek Baran's avatar
Hynek Baran committed
621
        op(map(proc(x,u) cat(u, separator, `bcount//str`(x, separator)) = 'jet'(u,x) end,
Hynek Baran's avatar
Hynek Baran committed
622 623 624
          aux, u))
      else
        op(map(proc(x,u)
Hynek Baran's avatar
Hynek Baran committed
625
            cat(op(0,u), separator, `bcount//str`(x, separator))[op(u)] = 'jet'(u,x)
Hynek Baran's avatar
Hynek Baran committed
626 627 628 629 630
          end, aux, u))
      fi
    end, flist, aux)));
end:

Hynek Baran's avatar
Hynek Baran committed
631
`bcount//str` := proc (x, separator)
Hynek Baran's avatar
Hynek Baran committed
632 633 634
  if type (x,'name') then x
  elif type (x,`^`) then 
    if op(2,x) > 3 then cat(op(2,x),op(1,x)) else op(1,x) $ op(2,x) fi
Hynek Baran's avatar
Hynek Baran committed
635
  elif type (x,`*`) then `bcount//str/_`(op(map(`bcount//str`, [op(x)], separator)), ':-separator'=separator)
Hynek Baran's avatar
Hynek Baran committed
636 637 638 639
  else ERROR (`not a count`, x)
  fi 
end:

Hynek Baran's avatar
Hynek Baran committed
640
`bcount//str/_` := proc ({separator::string:="_"})
Hynek Baran's avatar
Hynek Baran committed
641
  local i,a,b,ans;
Hynek Baran's avatar
Hynek Baran committed
642 643 644 645
  ans := _rest[1];
  for i from 2 to _nrest do
    a := _rest[i-1]; b := _rest[i];
    if length(a) + length(b) > 2 then ans := ans,separator,b else ans := ans,b fi
Hynek Baran's avatar
Hynek Baran committed
646 647 648 649 650 651 652
  od 
end:

`count//seq` := proc (x)
  if type (x,`^`) then `count//seq`(op(1,x)) $ op(2,x)
  elif type (x,`*`) then op(map(`count//seq`,[op(x)]))
#  elif type (x,symbol) then x
Hynek Baran's avatar
Hynek Baran committed
653
#  else cat(op(1,x),'_',`bcount//str`(op(2,x), separator)) 
Hynek Baran's avatar
Hynek Baran committed
654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689
  else x
  fi
end:

#
#  N a m e s
#

`name/tab` := table([]):

# Unassigned names may acquire a meaning through a declaration.
# Such names are said to be registered.
# For example, variables are registered names since they are declared
# through the coordinates command.

# Meanings other than 'variable' are stored in `name/tab`.
# registered() prints meanings of all names;
# registered(a) prints all names of meaning a;

# WARNING: After  a s s u m e,  names lose their registration and must be
# registered again.

registered := proc() 
  local a;
  if nargs = 0 then
    op(select(proc(e) op(1,e) = eval(op(1,e)) end, op(2,op(`name/tab`))))
  else a := args[1];
    op(map(proc(e,a) op(1,e) end, select(
      proc(e,a)
        op(1,e) = eval(op(1,e)) and type(op(1,e),a) 
      end, op(2,op(`name/tab`)), a), a))
  fi
end:

# To clear a declaration from all names:

Hynek Baran's avatar
Hynek Baran committed
690 691
clear := proc(opts::seq(name), {output::identical(expr, eq):=expr})
  local a,aux,ans,nex;
Hynek Baran's avatar
Hynek Baran committed
692
  forget(Vars);
Hynek Baran's avatar
Hynek Baran committed
693 694 695 696 697 698 699 700 701 702 703 704 705 706 707
  aux := {opts};
  #lprint(aux);
  #if aux minus {} <> {args} then
  #  ERROR(`not a name`, op({args} minus aux)) 
  #fi;
  
  if has(aux, puts) then aux := aux union {cc, pds, assignments}
  elif has(aux, pds) and `clear/pds/ccclearing`<>'suppress' then aux := aux union {cc} fi;
  
  aux, nex := selectremove(proc(a) Existing(`clear/`,a) end, aux);
  
  #lprint(aux,nex);
  #if aux <> {args} then
  #  ERROR(`invalid argument`, op({args} minus aux)) 
  #fi;
Hynek Baran's avatar
Hynek Baran committed
708
  ans := NULL;
Hynek Baran's avatar
Hynek Baran committed
709 710 711
  aux := sort(convert(aux,list)); # assignments is the first in the alphabet :)
  for a in aux do ans := ans, Call(`clear/`,a)(_options) od;
  return [ans];
Hynek Baran's avatar
Hynek Baran committed
712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742
end:

#
#   P a r a m e t e r s
#

parameter := proc()
  global `name/tab`;
  local a,aux;
  aux := select (type,{args},name);
  if aux <> {args} then ERROR(`not a name`, op({args} minus aux)) fi;
  aux := select(
    proc(a)
      type (a,'dep')
      or type (a,{`b/var`,`f/var`,'constant'})
      or (type (`name/tab`[a],'set')
        and `name/tab`[a] minus {'parameter','nonlocal'} <> {})
    end, {args});
  if aux <> {} then ERROR (`name(s) already used`, op(aux)) fi;
  for a in [args] do 
    if type (`name/tab`[a],'set') then
      `name/tab`[a] := `name/tab`[a] union {'parameter'}
    else `name/tab`[a] := {'parameter'}
    fi
  od;
  registered('parameter')
end:

`clear/parameter` := proc()
  local a,aux;
  global `name/tab`;
Hynek Baran's avatar
Hynek Baran committed
743
  print(`Clearing parameter's.`);
Hynek Baran's avatar
Hynek Baran committed
744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786
  aux := {registered('parameter')};
  for a in aux do
    `name/tab`[a] := `name/tab`[a] minus {'parameter'} 
  od;
  op(aux)
end:

parameters := proc()
  `clear/parameter`();
  parameter(args)
end:

# Parameters are of type 'parameter':

`type/parameter` := proc (x)
  if type (`name/tab`[x],'set') then
    member ('parameter',`name/tab`[x])
  else false
  fi
end: 

# Type 'ar' is either 'parameter' or 'var':

`type/ar` := {'var','parameter'}:

`type/ar/count` := proc (x) 
  local x1;  
  if x = 1 then false 
  elif type (x,`^`) then type (op(1,x),'ar') and type (op(2,x), posint) 
  elif type (x,`*`) then
    for x1 in x do if not type (x1,`ar/count`) then RETURN(false) fi od;
    true
  else type (x,'ar')
  fi
end:

#
#   D e p e n d e n c e
#

# Names may depend on variables:

dependence := proc ()
Hynek Baran's avatar
Hynek Baran committed
787
  local f,fs,es,e, opts; 
Hynek Baran's avatar
Hynek Baran committed
788
  global `dep/tab`;
Hynek Baran's avatar
Hynek Baran committed
789
  forget(Vars);
Hynek Baran's avatar
Hynek Baran committed
790 791
  fs := select(type, {args}, function);
  es := select(type, {args}, `=`);
Hynek Baran's avatar
Hynek Baran committed
792 793
  opts := select(type, {args}, symbol); # HB
  if fs union es union opts <> {args} then # HB
Hynek Baran's avatar
Hynek Baran committed
794 795 796 797 798 799 800 801 802 803
    ERROR (`wrong arguments`, op({args} minus fs minus es)) 
  fi; 
  es := es union map(proc(f) op(0,f) = {op(f)} end, fs);
  fs := select(
    proc(e)
      type (op(1,e),{`b/var`,`f/var`,'constant','parameter'})
    end, es);
  if fs <> {} then
    ERROR (`name(s) already used`, op(fs))
  fi; 
Hynek Baran's avatar
Hynek Baran committed
804

Hynek Baran's avatar
Hynek Baran committed
805 806 807 808 809 810 811 812
  for e in es do 
    if type (op(2,e), set('ar')) then `dep/tab`[op(1,e)] := op(2,e)
    else ERROR (`forbidden type of dependence`)
    fi
  od;
  refresh(); 
  es := {op(op(2,op(`dep/tab`)))};
  es := select(proc(e) evalb(op(1,e) = eval(op(1,e))) end, es);
Hynek Baran's avatar
Hynek Baran committed
813 814 815 816 817 818 819
  # HB:
  if `subs` in opts then
    seq((op(1,e)=op(1,e)(op(op(2,e)))), e = es) 
  else  
  # :HB
    seq(op(1,e) = op(2,e), e = es)
  fi   
Hynek Baran's avatar
Hynek Baran committed
820 821 822 823 824 825 826 827 828
end:

`dep/tab` := table ():

`type/dep` := proc(a) type(`dep/tab`[a],'set') end:

`clear/dependence` := proc()
  local e,el,ans;
  global `dep/tab`;
Hynek Baran's avatar
Hynek Baran committed
829
  print(`Clearing dependence's.`);
Hynek Baran's avatar
Hynek Baran committed
830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848
  el := op(2,op(`dep/tab`));
  ans := {seq((op(1,e))(op(op(2,e))), e = el)};
  `dep/tab` := table([]);
  ans 
end:

dependences := proc()
  `clear/dependence`();
  dependence(args)
end:

#
#   T o t a l   d e r i v a t i v e
#

TD := proc (f, x::`b/count`)
  global `EVALTD/s`,`EVALTD/b`,`EVALTD/n`;
  if nargs = 1 then RETURN (`TD~`(procname,f)) fi;
  if nargs > 2 then RETURN (TD(TD(f, args[2..nargs-1]), args[nargs])) fi;
Hynek Baran's avatar
Hynek Baran committed
849
  if type (f,{'numeric', 'complex'}) then 0 # HB 5. 5. 2008
Hynek Baran's avatar
Hynek Baran committed
850 851 852 853 854 855
  elif type (f,'name') then
    if member (f,[constants]) then 0
    elif type (f,`b/var`) then if f = x then 1 else 0 fi
    elif type (f,`f/var`) then jet(f,x)
    elif type (f,'parameter') then 0
    elif type (f,'dep') then
Hynek Baran's avatar
Hynek Baran committed
856
      if `vars/1`(f) minus `b/var/s` = {} then pd(f,x)
Hynek Baran's avatar
Hynek Baran committed
857 858 859 860 861 862 863 864 865 866 867 868 869 870 871
      elif not `EVALTD/b` and not member(f,`EVALTD/s`)
        and `count/length`(x) > `EVALTD/n` then 'TD'(f,x) 
      elif type (x,'name') then 
        if type (f,'dep') then `TD/dep/s`(f,x,`dep/tab`[f]) fi
      else TD (TD(f,`count/f`(x)), `count/r`(x)) 
      fi
    elif Existing(`TD/indexed/`,op(0,f)) then
      Call(`TD/indexed/`,op(0,f))(op(f),x)
    else 'TD'(f,x)
    fi
  elif type (f,`+`) then map (procname, f, x)  
  elif type (f,'function') then 
    if type (f, `j/var`) then jet(op(1,f),x*op(2,f)) 
    elif type (f, specfunc(anything,pd)) then
      if type (op(1,f),'dep') then 
Hynek Baran's avatar
Hynek Baran committed
872
        if `vars/1`(op(1,f)) minus `b/var/s` = {} then pd(op(1,f), x*op(2,f))
Hynek Baran's avatar
Hynek Baran committed
873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915
        elif not `EVALTD/b` and not member (op(1,f),`EVALTD/s`)
          and `count/length`(x) > `EVALTD/n` then 'TD'(f,x)
        elif type (x,'name') then `TD/dep/s`(f,x,`dep/tab`[op(1,f)]) 
        else TD (TD(f,`count/f`(x)), `count/r`(x)) 
        fi
      else 'TD'(f,x)
      fi
    elif type (f, specfunc(anything,TD)) then 'TD'(op(1,f), x*op(2,f))
    elif type (x,'name') then
      if Existing(`der/`,op(0,f)) then
        Call(`der/`,op(0,f))(procname,op(f),x)
      else ERROR (`unexpected function`, f) 
      fi
    else TD (TD(f,`count/f`(x)), `count/r`(x))
    fi
  elif type (x,'name') then
    if type (f,`*`) then `der/*` (procname, f, x)
    elif type (f,`^`) then `der/^` (procname, op(1,f), op(2,f), x)
    else ERROR (`unexpected object`, f) 
    fi
  else TD (TD(f,`count/f`(x)), `count/r`(x))
  fi
end:

`TD~` := proc (pn,f)
  if type (pn,'indexed') then 
    if nops(pn) <> 1 then ERROR (`too many indices`, op(pn)) 
    else RETURN (TD(f,op(1,pn)))
    fi
  else eval(f)
  fi
end:

`TD/dep/s` := proc (f,x,s)
  convert(map(proc(p,f,x) pd(f,p)*TD(p,x) end, [op(s)], f,x), `+`)
end:

# Evaluating total derivatives

evalTD := proc (f) 
  local a,b,`EVALTD/b/old`,`EVALTD/s/old`;
  global `EVALTD/b`, `EVALTD/s`;
  `EVALTD/b/old` := `EVALTD/b`; `EVALTD/s/old` := `EVALTD/s`;
Hynek Baran's avatar
Hynek Baran committed
916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933
  if type(f,sequential) or type(f, matrix) then # HB
    map('procname', f, args[2..-1])
  else
    if nargs = 1 then `EVALTD/b` := true
		elif type (args[2], set(name)) then `EVALTD/b` := false;
		 `EVALTD/s` := args[2]
		else ERROR (`second argument must be a set of names`)
		fi; 
		b := eval(f); 
		a := {}; # a either different from a or evalTD(a) = a
		while hasfun(b,TD) and b <> a do
			a := b;
			b := traperror(eval(a))
		od;
		`EVALTD/b` := `EVALTD/b/old`; `EVALTD/s` := `EVALTD/s/old`;
		if b = lasterror then ERROR (lasterror) fi;
		b
  fi;
Hynek Baran's avatar
Hynek Baran committed
934 935
end: 

Hynek Baran's avatar
Hynek Baran committed
936 937


Hynek Baran's avatar
Hynek Baran committed
938 939 940 941 942 943 944 945 946 947 948
`EVALTD/s` := {}:
`EVALTD/b` := false:
`EVALTD/n` := 0:

#
#   P a r t i a l   d e r i v a t i v e s
#

pd := proc (f, p::`ar/count`)
  if nargs = 1 then RETURN (`pd~`(procname,f)) fi;
  if nargs > 2 then RETURN(pd(pd(f, args[2..nargs-1]), args[nargs])) fi;
Hynek Baran's avatar
Hynek Baran committed
949
  if type (f,{'numeric', 'complex'}) then 0 # HB 5. 5. 2008
Hynek Baran's avatar
Hynek Baran committed
950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996
  elif type (f,'name') then
    if member (f,[constants]) then 0
    elif type (f,'ar') then if f = p then 1 else 0 fi
    elif Existing (`pd/indexed/`,op(0,f)) then
      Call(`pd/indexed/`,op(0,f))(op(f),p)
    else `pd//tab`(f,p)
    fi
  elif type (f,`+`) then map (procname, f, p)  
  elif type (f,'function') then 
    if type (f, `j/var`) then if f = p then 1 else 0 fi
    elif type (f, specfunc(anything,pd)) then `pd//tab`(op(1,f), op(2,f)*p) 
    elif type (p,'ar') then
      if type (f, specfunc(anything,TD)) then `pd/TD`(op(f), p)     
      elif Existing(`der/`,op(0,f)) then
        Call(`der/`,op(0,f)) (procname,op(f),p)
      else `der//unknown`(procname,op(0,f),op(f),p)  # M.M. 26.11.2004
      fi
    else pd (pd(f,`count/f`(p)), `count/r`(p))
    fi
  elif type (p,'ar') then
    if type (f,`*`) then `der/*` (procname, f, p)
    elif type (f,`^`) then `der/^` (procname, op(1,f), op(2,f), p)
    else ERROR (`unexpected object`, f) 
    fi
  else pd (pd(f,`count/f`(p)), `count/r`(p))
  fi
end:

`pd~` := proc (pn,f)
  if type (pn,'indexed') then 
    if nops(pn) <> 1 then ERROR (`too many indices`, op(pn)) 
    else RETURN (pd(f,op(1,pn)))
    fi
  else eval(f)
  fi
end:

`pd/TD` := proc (f,x,p)
  option remember; 
  if type(x, 'name') then `pd/TD/1`(f,x,p)
  else `pd/TD/1`(evaln(TD(f,`count/r`(x))), `count/f`(x), p)
  fi
end:

`pd/TD/1` := proc (f,x,p)
  option remember; 
  local qs;
Hynek Baran's avatar
Hynek Baran committed
997
  qs := [op(`vars/1`(f))];
Hynek Baran's avatar
Hynek Baran committed
998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057
  TD(pd(f,p),x) + convert(map(
    proc(q,f,x,p)
      local qp; qp := pd(TD(q,x),p); 
      if qp <> 0 then pd(f,q)*qp else NULL fi end,
    qs, f, x, p), `+`)
end:

#
#   D e r i v a t i v e s   o f   f u n c t i o n s 
#

# For every concrete nary function  a  there should exist a function
# `der/f` of arguments  der,f1,...,fn,x  which returns the derivative
# der of a(f1,...,fn) with respect to x. The actual parameter der will
# eventually be one of pd, TD.

`der/` := proc (der,f,x) der(f,x) end:

`der/Diff` := proc (der,f,x) der(f,x) end:

`der/*` := proc (der,f,x) local i,s,di;
  s := 0;
  for i to nops(f) do di := der(op(i,f),x);
    if di <> 0 then s := s + subsop (i=di, f) fi
  od; s
end:

`der/&*` := proc ()
  local d,f,x,i,s,di;
  d := args[1];
  f := &*(args[2..nargs-1]); x := args[nargs];
  s := 0;
  for i to nops(f) do di := d(op(i,f),x);
    if di <> 0 then s := s + subsop (i=di, f) fi
  od; s
end:

`der/^` := proc (der,f,g,x)
  if type (g,'integer') then g*f^(g-1)*der(f,x)
  else g*f^(g-1)*der(f,x) + ln(f)*f^g*der(g,x)
  fi
end:

`der//unknown` := proc (der,f)  # M.M. 26.11.2004
  local as,x,i,ind;
  global DER;
  if type(f,'indexed') and op(0,f) = DER then 
    ind := op(f);
    if nargs = 4 then
      as := args[3]
    else ERROR(`wrong number of arguments`)
    fi;
    x := args[nargs];
    add(DER[op(sort([i,ind]))](as)*der(op(i,as), x), i = 1..(nargs - 3))
  else
    as := [args[3..(nargs - 1)]];
    x := args[nargs];
    add(DER[i](f(op(as)))*der(as[i], x), i = 1..(nargs - 3))
    fi
end:
Hynek Baran's avatar
Hynek Baran committed
1058 1059 1060 1061 1062 1063 1064
# HB:
# there is an example how to specify derivatives of general function:
`der/MyFun`    := proc (der,f,x) 'D_MyFun'(f)*der(f,x) end:
`der/D_MyFun`  := proc (der,f,x) 'D2_MyFun'(f)*der(f,x) end:
`der/D2_MyFun` := proc (der,f,x) 'D3_MyFun'(f)*der(f,x) end:
`der/D3_MyFun` := proc (der,f,x) 'D4_MyFun'(f)*der(f,x) end:
# :HB
Hynek Baran's avatar
Hynek Baran committed
1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092

`der/exp` := proc (der,f,x) exp(f)*der(f,x) end:
`der/ln`  := proc (der,f,x) der(f,x)/f end:

`der/sin` := proc (der,f,x) cos(f)*der(f,x) end:
`der/cos` := proc (der,f,x) -sin(f)*der(f,x) end:
`der/tan` := proc (der,f,x) der(f,x)/cos(f)^2 end:
`der/cot` := proc (der,f,x) -der(f,x)/sin(f)^2 end:
`der/csc` := proc (der,f,x) -cos(f)*der(f,x)/sin(f)^2 end:
`der/sec` := proc (der,f,x) sin(f)*der(f,x)/cos(f)^2 end:

`der/arcsin` := proc (der,f,x) der(f,x)/sqrt(1-f^2) end:
`der/arccos` := proc (der,f,x) -der(f,x)/sqrt(1-f^2) end:
`der/arccot` := proc (der,f,x) -der(f,x)/(1+f^2) end:
`der/arctan` := proc (der,f,x) der(f,x)/(1+f^2) end:

`der/sinh` := proc (der,f,x) cosh(f)*der(f,x) end:
`der/cosh` := proc (der,f,x) sinh(f)*der(f,x) end:
`der/tanh` := proc (der,f,x) der(f,x)/cosh(f)^2 end:
`der/coth` := proc (der,f,x) -der(f,x)/sinh(f)^2 end:
`der/csch` := proc (der,f,x) -cosh(f)*der(f,x)/sinh(f)^2 end:
`der/sech` := proc (der,f,x) -sinh(f)*der(f,x)/cosh(f)^2 end:

`der/arcsinh` := proc (der,f,x) der(f,x)/sqrt(1+f^2) end:
`der/arccosh` := proc (der,f,x) der(f,x)/sqrt(f^2-1) end:
`der/arccoth` := proc (der,f,x) der(f,x)/(f^2-1) end:
`der/arctanh` := proc (der,f,x) der(f,x)/(1-f^2) end:

Hynek Baran's avatar
Hynek Baran committed
1093 1094 1095 1096
`der/dilog` := proc (der,f,x) # AS
  (ln(f)/(1-f))*der(f,x)
end:

Hynek Baran's avatar
Hynek Baran committed
1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160
`der/LambertW` := proc (der,f,x)
  LambertW(f)/(1 + LambertW(f))*der(f,x)/f
end:

`der/AiryAi` := proc(der,n,f,x) 
  if nargs = 3 then AiryAi(1,n)*der(n,f) 
  elif n = 1 then AiryAi(f)*f*der(f,x)
  elif pd(n,x) = 0 then AiryAi(n + 1, f)*der(f,x)
  else 'der'(AiryAi(n,f),x)
  fi 
end:

`der/AiryBi` := proc(der,n,f,x) 
  if nargs = 3 then AiryBi(1,n)*der(n,f) 
  elif n = 1 then AiryBi(f)*f*der(f,x)
  elif pd(n,x) = 0 then AiryBi(n + 1, f)*der(f,x)
  else 'der'(AiryBi(n,f),x)
  fi 
end:

`der/RootOf` := proc(der,f,x) 
  if has(f,x) then ERROR(`Sorry, this version of Jets cannot handle`) 
  else 0
  fi 
end:

# Printing partial derivatives

`print/pd` := proc (q,x) 'Diff'(q,`count//seq`(x)) end: 

# Diff is mapped to pd

unprotect(Diff):

Diff := proc(f) 
  if nargs = 1 then
    'Diff(f)' # Maple procedures may use Diff for other purposes
  else pd(f,convert([op(2..nargs, [args])], `*`))
  fi 
end:

#
#   A s s i g n i n g   p a r t i a l   d e r i v a t i v e s 
#

# Values assigned to partial derivatives are stored in a table.

`pd/tab` := table([]):

# Extracting values from the pd table:
# `pd//tab`(f,p) returns the value of pd(f,p).
# If ambiguous, chooses that of minimal size

`pd//tab` := proc (f,p)
  option remember; 
  local el,el1,e,ans,m,q,ql,rt,lb; 
  global `pd/tab`,`dep/tab`; 
  if type (f,'dep') and ars(p) minus `dep/tab`[f] <> {} then RETURN(0) fi;
  lb := `PD:`; rt := `report/tab`[pd];
  if assigned(`pd/tab`[f]) then
    if assigned(`pd/tab`[f][p]) then RETURN(`pd/tab`[f][p]) fi;
    el := op(2,op(`pd/tab`)[f]); # list of already assigned derivatives
    el := map(proc(e,p) p/op(1,e) = op(2,e) end, el, p); # reciprocals
    el := select(proc(e) type (op(1,e),'count') end, el); # select counts 
Hynek Baran's avatar
Hynek Baran committed
1161
    if el = [] then RETURN('pd'(f,p)) fi; # if no counts at all
Hynek Baran's avatar
Hynek Baran committed
1162 1163
    el1 := select(proc(e) not type (op(1,e),{`*`,`^`}) end, el); # 1st order 
    if el1 <> [] then
Hynek Baran's avatar
Hynek Baran committed
1164
      if select(type, el1, anything='numeric') <> [] then ans := 0  
Hynek Baran's avatar
Hynek Baran committed
1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207
      else
        if rt > 3 then report(lb,[`extending table 1 step: `, nops(el1)]) fi;
        el1 := map(proc(e) Simpl(eval(pd(op(2,e),op(1,e)))) end, el1); 
        if nops(el1) > 1 then
          el1 := sizesort(el1,size);
          ans := op(1,el1); # the smallest size derivative
        else
          ans := el1[1]
        fi
      fi;
      if rt > 3 then report(lb,[`storing a value for`, 'pd'(f,p)]) fi;
      `put/1`('pd'(f,p), ans);
      RETURN (ans)
    else # no 1st order counts
      if select(type, el, anything='numeric') <> [] then
        ans := 0; `put/1`('pd'(f,p), ans); RETURN (ans)
      else
        el := map(proc(e,p) [op(e),`count/length`(op(1,e))] end, el);
        el := sort(el, proc(e1,e2) evalb(op(3,e1) < op(3,e2)) end);
        m := op(3,el[1]); # length of the shortest count
        el := select(proc(e,m) evalb(op(3,e) = m) end, el, m); 
        for e in el do # running through lowest (= m-th) order derivatives
          ql := `pd//tab/1`(op(1,e)); # list of indeterminates
          for q in ql do ans := pd(op(2,e),q);
            if rt > 3 then report(lb,[`storing `, 'pd'(f,p/op(1,e)*q)]) fi;
            `put/1`('pd'(f,p/op(1,e)*q), ans)
          od
        od;
        RETURN(`pd//tab`(f,p))
      fi
    fi 
  fi;
  'pd'(f,p)
end:

`pd//tab/1` := proc(q) # count -> list of indeterminates
  if type (q,`*`) then 
    map(proc(z) if type (z,`^`) then op(1,z) else z fi end, [op(q)]);
  elif type (q,`^`) then [op(1,q)]
  else [q]
  fi
end:

Hynek Baran's avatar
Hynek Baran committed
1208 1209 1210 1211 1212 1213 1214 1215 1216 1217
`pd/tab/print/unk` := proc(u) # print all assigments of unknown u stored in pd/tab
  global `pd/tab`;
  if assigned(`pd/tab`[u]) then
    map(i -> lprint('pd'(u,op(i))=pd(u,op(i))), [indices(`pd/tab`[u])])
  else
    printf("%a not found in pd/tab\n", u);
  fi
end:


Hynek Baran's avatar
Hynek Baran committed
1218 1219 1220 1221 1222 1223 1224 1225
#`pd/decomp` := proc(e)
#	 # decomposes pd(f,p) to (f,p) and f to (f,1)
#   if type (e,'name') then (e,1)
#   elif type (e,specfunc(anything,'pd')) then	(op(1,e), op(2,e))
#  else error e, "cannot be decomposed, its not symbol or pd" fi;
#end:


Hynek Baran's avatar
Hynek Baran committed
1226

Hynek Baran's avatar
Hynek Baran committed
1227 1228

# put assigns values to pd's via `pd/tab` or `dep/tab`
Hynek Baran's avatar
Hynek Baran committed
1229 1230
put := proc()
  local e;
Hynek Baran's avatar
Hynek Baran committed
1231
  forget(Vars);
Hynek Baran's avatar
Hynek Baran committed
1232
  if type ([args], 'list'(`=`)) then
Hynek Baran's avatar
Hynek Baran committed
1233 1234
    for e in args do 
    	#`put/items/add`(lhs(e));
Hynek Baran's avatar
Hynek Baran committed
1235
    	`put/1`(op(e), _rest) 
Hynek Baran's avatar
Hynek Baran committed
1236
    od
Hynek Baran's avatar
Hynek Baran committed
1237 1238 1239 1240 1241 1242
  else ERROR (`wrong arguments`)
  fi;
  `put/evaluate`();
  NULL
end:

Hynek Baran's avatar
Hynek Baran committed
1243

Hynek Baran's avatar
Hynek Baran committed
1244
unassign('`put/limit/size`', '`put/limit/length`'): # use with extreme care!
Hynek Baran's avatar
Hynek Baran committed
1245

Hynek Baran's avatar
Hynek Baran committed
1246 1247 1248 1249 1250 1251
`put/1` := proc(p,a) # {`put/limit/size`::{numeric,infinity}:=infinity, `put/limit/length`::{numeric,infinity}:=infinity}
  global `report/tab`, `unk/s`,`unk/<</list`,`dep/tab`,`pd/tab`, `put/name/tab`, `put/limit/size`, `put/limit/length`;
  Report(0, [p]);
  Report(1, [p=size(a)]);
  Report(2, [p=[LVar(a),size(a)]]);
  Report(4, [p=a]);
Hynek Baran's avatar
Hynek Baran committed
1252 1253 1254
  if type (p,'name') then
    `unk/s` := map(`put/1/remove`, `unk/s`, p);
    `unk/<</list` := map(`put/1/remove`, `unk/<</list`, p);
Hynek Baran's avatar
Hynek Baran committed
1255
     `put/name/tab`[p] := a;
Hynek Baran's avatar
Hynek Baran committed
1256
    assign(p = a);
Hynek Baran's avatar
Hynek Baran committed
1257 1258 1259
  elif type (p,specfunc(anything,'pd')) then
    if a = 0 and type (op(1,p),'dep') and type(op(2,p),'var') then 
      `dep/tab`[op(1,p)] := `dep/tab`[op(1,p)] minus {op(2,p)}
Hynek Baran's avatar
Hynek Baran committed
1260
    else     	  
Hynek Baran's avatar
Hynek Baran committed
1261 1262 1263 1264
       if type(`put/limit/length`, numeric) and length(a) > `put/limit/length` then
          printf("put: Ignoring too long (size=%a, length=%a>%a) argument %a=%a\n", size(a), length(a), `put/limit/length`, p, a);   
       elif type(`put/limit/size`, numeric) and size(a) > `put/limit/size` then
          printf("put: Ignoring too large (size=%a>%a, length=%a) argument %a=%a\n", size(a),`put/limit/size`, length(a), p, a);             	  
Hynek Baran's avatar
Hynek Baran committed
1265
       else
Hynek Baran's avatar
Hynek Baran committed
1266 1267 1268
      	 if type(`put/limit/size`, numeric) or type(`put/limit/length`, numeric) then 
      	   printf("put: Putting (size=%a<=%a, length=%a<=%a) %a=%a\n", size(a),`put/limit/size`, length(a), `put/limit/length`, p, a); 
      	 fi;
Hynek Baran's avatar
Hynek Baran committed
1269 1270 1271 1272 1273
         `pd/tab`[op(1,p)][op(2,p)] := a;
       fi;	
    fi;
  else 
    if simpl(p-a) = 0 then 
Hynek Baran's avatar
Hynek Baran committed
1274
      printf ("ignoring redundant put of %q\n", p);
Hynek Baran's avatar
Hynek Baran committed
1275
    else
Hynek Baran's avatar
Hynek Baran committed
1276 1277
      lprint(p);
      lprint(a);
Hynek Baran's avatar
Hynek Baran committed
1278
      error "Contradiction in put( %1 = %2 )", p, a;
Hynek Baran's avatar
Hynek Baran committed
1279
    fi;
Hynek Baran's avatar
Hynek Baran committed
1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312
  fi
end:

`put/1/remove` := proc(a,b) if a = b then NULL else a fi end:

`put/evaluate` := proc()
  local e,el,f,ft,q,qs,eq,rpt; 
  global `dep/tab`,`pd/tab`,`nonzero/s`;
# `refresh/1`(`pd//tab`);
  refresh();
  el := op(2,op(`pd/tab`));
  rpt := table([]);
  for e in el do
    f := op(1,e); ft := op(2,e);
    qs := map(proc(x) op(1,x) end, op(2,op(ft)));
    rpt[f] := table([]);
    for q in qs do eq := Simpl(eval(ft[q]));
      if type(q,'var') and eq = 0 and type (f,'dep') then
        `dep/tab`[f] := `dep/tab`[f] minus {q}
      else rpt[f][q] := eq
      fi
    od
  od;
  `pd/tab` := op(rpt);
  nonzero();
  map(
    proc(a)
      if Simpl(eval(a)) = 0 then
        ERROR(`declared nonzero became zero`, a)
      fi
    end, `nonzero/s`);
end:

Hynek Baran's avatar
Hynek Baran committed
1313 1314 1315 1316 1317 1318 1319 1320 1321
`puts/assignments` := proc()
  global `put/name/tab`;
  op(op(op(`put/name/tab`)));  # 'unk_i'=value_i, 'unk_j'=value_j, ... (unordered)
end:

`clear/assignments` := proc({output::identical(expr, eq):=expr})
  description "unassign all symbols assigned by put()";
  global `put/name/tab`;
  local as;
Hynek Baran's avatar
Hynek Baran committed
1322
  forget(Vars);
Hynek Baran's avatar
Hynek Baran committed
1323 1324
  print(`Clearing assignments.`);
  as := `puts/assignments`();
Hynek Baran's avatar
Hynek Baran committed
1325 1326
  map(unassign@op, [indices(`put/name/tab`)]); # unassign assigned symbols
  `put/name/tab` := table(); # clear the table of assigned symbols
Hynek Baran's avatar
Hynek Baran committed
1327 1328 1329 1330 1331 1332 1333 1334 1335
  if output=eq then 
    return as;
  elif output=expr then
    return op(map(lhs-rhs, [as])) ;
  else 
    error "Unknown output type %1", output; 
  fi;
end:

Hynek Baran's avatar
Hynek Baran committed
1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360

#`put/items/add` := proc(e)
#  global `put/items/table`;
#  local f, p;
#  f, p := `pd/decomp`(e);
#  `put/items/table`[f] := `put/items/table`[f] union {p};
#end:
#
#`index/put/items/table` := proc(Idx::list,Tbl::table,Entry::list)
#    if (nargs = 2) then # get returns {} if unassigned
#        if assigned(Tbl[op(Idx)]) then Tbl[op(Idx)];
#        else {};
#        end if;
#    else # assign
#        Tbl[op(Idx)] := op(Entry);
#    end if;
#end proc:
#
#`clear/put/items` := proc()
#	 global `put/items/table`;
#	 `put/items/table` := table(`put/items/table`):
#end:
#
#`clear/put/items`():

Hynek Baran's avatar
Hynek Baran committed
1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373
# To convert the pd table into a list:

pds := proc()
  op(map(`pds/1`, op(2,op(`pd/tab`))))
end:

`pds/1` :=  proc(p)
   op(map(
     proc(e,f) 
       'pd'(f,op(1,e)) = op(2,e) 
     end, op(2,op(2,p)), op(1,p))) 
end:

Hynek Baran's avatar
Hynek Baran committed
1374
# To convert the pd table into a list of expressions while clearing all
Hynek Baran's avatar
Hynek Baran committed
1375 1376
# assignments to pd's:

Hynek Baran's avatar
Hynek Baran committed
1377 1378
`clear/pds` := proc({output::identical(expr, eq):=expr})
  description "unassign all pd's and cc's assigned by put()";
Hynek Baran's avatar
Hynek Baran committed
1379
  local t,aux;
Hynek Baran's avatar
Hynek Baran committed
1380
  global `pd/tab`, `clear/pds/ccclearing`;
Hynek Baran's avatar
Hynek Baran committed
1381
  forget(Vars);
Hynek Baran's avatar
Hynek Baran committed
1382
  print(`Clearing pds.`);
Hynek Baran's avatar
Hynek Baran committed
1383 1384 1385 1386
  reduce();
  t := copy(`pd/tab`);
  `pd/tab` := table([]);
  refresh();
Hynek Baran's avatar
Hynek Baran committed
1387
  #`clear/put/items`():
Hynek Baran's avatar
Hynek Baran committed
1388 1389 1390 1391 1392 1393 1394 1395 1396
  aux := map(`pds/1`, op(2,op(t)));
  aux := sizesort(aux, size);
  if output=eq then 
    return op(aux);
  elif output=expr then
    return op(map(lhs-rhs, aux)) ;
  else 
    error "Unknown output type %1", output; 
  fi;
Hynek Baran's avatar
Hynek Baran committed
1397 1398 1399 1400
end:

reduce := proc()
  local ans1,ans2;
1401
  print(`Reducing...`);
Hynek Baran's avatar
Hynek Baran committed
1402 1403 1404 1405 1406 1407 1408 1409 1410 1411
  ans1 := `reduce/pd`();
  ans2 := `unks/assignments`();
  refresh();
  ans1,ans2 
end:

`unks/assignments` := proc()
  global `unk/tab`;
  local aux;
  aux := map(op,[entries(`unk/tab`)]);
Hynek Baran's avatar
Hynek Baran committed
1412
  op(map(proc(f) if f <> eval(f) then f = eval(f) fi end, aux)) # 1 = 'unk_1', 2='unk_2', ... (in the order of unknowns())
Hynek Baran's avatar
Hynek Baran committed
1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426
end:

`reduce/pd` := proc()
  local e,el,f,ft,p,pl,q,qs,qe,eq,rpt; 
  global `dep/tab`,`pd/tab`;
  el := op(2,op(`pd/tab`));
  rpt := table([]);
  for e in el do
    f := op(1,e); ft := op(2,e);
    pl := map(proc(x) op(1,x) end, op(2,op(ft)));
    if f = eval(f) then
      # remove by independence
      if type(f,'dep') then
        pl := select(proc(p,f) global `dep/tab`;
Hynek Baran's avatar
Hynek Baran committed
1427
            evalb(`vars/1`(p) minus `dep/tab`[f] = {}) 
Hynek Baran's avatar
Hynek Baran committed
1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456
          end, pl, f)
      fi;
      qs := {op(pl)}; qe := {};
      # remove by pd consequence
      for p in pl do
        for q in qs do
          if type (q/p,'count') then qe := qe union {q} fi;
        od;
        qs := qs minus qe
      od
    else 
      # remove by assignment
      qs := {}; eq := eval(f)
    fi;
    # evaluate
    rpt[f] := table([]);
    for q in qs do eq := Simpl(eval(ft[q]));
      if type(q,'var') and eq = 0 and type (f,'dep') then
        `dep/tab`[f] := `dep/tab`[f] minus {q}
      else rpt[f][q] := eq
      fi
    od
  od;
  `pd/tab` := op(rpt);
  pds()
end:

# Checking still unresolved compatibility conditions

Hynek Baran's avatar
Hynek Baran committed
1457 1458 1459 1460
`cc/time` := 0.0 :

# The old cc implementation
# This cc implementation is deprecated but still available, 
Hynek Baran's avatar
Hynek Baran committed
1461 1462
# if you need to use it globally, set 
# cc := proc() `cc/old`(args, _options) end:
Hynek Baran's avatar
Hynek Baran committed
1463 1464
# or you may call `cc/old`() at any moment without any side-effects

Hynek Baran's avatar
Hynek Baran committed
1465 1466
cc := proc() `cc/new`(args, _options) end:

Hynek Baran's avatar
Hynek Baran committed
1467 1468 1469 1470
`cc/old` := proc()
  local e,el,f,ft,zl,p,pl,q,i,j,p1,q1,z,ans,sl,aux,rt,lb,time0; 
  global `dep/tab`,`pd/tab`,`cc/s`,`cc/tab`,`cc/time`;
  time0 := time();
Hynek Baran's avatar
Hynek Baran committed
1471 1472 1473
  ans := NULL; 
  zl := [];
  lb := `CC:`; rt := `report/tab`[cc];
Hynek Baran's avatar
Hynek Baran committed
1474
  
Hynek Baran's avatar
Hynek Baran committed
1475
  el := op(2,op(`dep/tab`)); # list of stored dependences
Hynek Baran's avatar
Hynek Baran committed
1476
  if rt > 3 then report(lb,["looking for cc - dep/tab"]) fi;
Hynek Baran's avatar
Hynek Baran committed
1477 1478 1479 1480 1481 1482
  for e in el do 
    f := op(1,e); # an unknown
    ft := op(2,e); # its dependence set
    if rt > 3 then report(lb,[`dependence set:`, 'f' = f, 'ft'= ft]) fi;
    aux := NULL;
    if f <> eval(f) then # f assigned
Hynek Baran's avatar
Hynek Baran committed
1483
      for p in `vars/1`(eval(f)) minus ft do
Hynek Baran's avatar
Hynek Baran committed
1484 1485 1486 1487 1488 1489 1490
        # pd(f,p) should be zero
        aux := aux, [f, eval(f), p, 0, 1, 1] 
      od;
      zl := [op(zl), aux];
      if rt > 2 then report(lb,[f, cat(`ass/dep: `,nops(zl),` c.c.`)]) fi
    fi
  od;  
Hynek Baran's avatar
Hynek Baran committed
1491 1492
  
  if rt > 2 then report(lb,["looking for cc - pd/tab"]) fi;
Hynek Baran's avatar
Hynek Baran committed
1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508
  el := op(2,op(`pd/tab`)); # list of stored partial derivatives
  for e in el do
    f := op(1,e); ft := op(2,e);
    if rt > 2 then report(lb,['f' = f]) fi;
    pl := map(proc(x) op(1,x) end, op(2,op(ft)));
    aux := NULL;
    if f <> eval(f) then # f assigned
      for p in pl do 
        # pd(f,p) should have the expected value
        aux := aux, [f, f, p, ft[p], 1, `count/length`(p)]; 
      od;
      zl := [op(zl), aux];
      if rt > 2 then report(lb,[f, cat(`ass/pd: `,nops([aux]),` c.c.`)]) fi
    else 
      if type(f,'dep') then # dependence
        for p in pl do
Hynek Baran's avatar
Hynek Baran committed
1509
          if `vars/1`(p) minus `dep/tab`[f] <> {} then
Hynek Baran's avatar
Hynek Baran committed
1510 1511 1512 1513 1514 1515
              # pd(f,p) should be zero 
            if ft[p] <> 0 then
              aux := aux, [f, ft[p], 1, 0, 1, `count/length`(p)]
            fi
          fi; 
          if ft[p] <> 0 then
Hynek Baran's avatar
Hynek Baran committed
1516
            for q in `vars/1`(ft[p]) minus `dep/tab`[f] do
Hynek Baran's avatar
Hynek Baran committed
1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530
              # pd(f,p*q) should vanish for each q from pd(f,p)
              aux := aux, [f, ft[p], q, 0, 1, `count/length`(p*q)]
            od
          fi
        od
      fi;
      zl := [op(zl), aux];
      if rt > 2 then report(lb,[f, cat(`dep/pd: `,nops([aux]),` c.c.`)]) fi;
      # cross differentiation
      aux := NULL;
      for i to nops(pl) do
        for j to i - 1 do p := op(i,pl); q := op(j,pl);
          if ft[p] <> 0 or ft[q] <> 0 then 
            p1 := `transform/count`(p/q); q1 := `transform/count`(q/p);
Hynek Baran's avatar
Hynek Baran committed
1531 1532
            if p*q1 <> q*p1 then ERROR(`this should not happen`) fi; # NejSpolNas
            aux := aux, [f, ft[p], q1, ft[q], p1, `count/length`(p*q1)] # ii kind
Hynek Baran's avatar
Hynek Baran committed
1533 1534 1535 1536 1537 1538 1539
          fi
        od
      od;
      if rt > 2 then report(lb,[f, cat(`cross-derivatives: `,nops([aux]))]) fi;
      zl := [op(zl), aux];
    fi
  od;  
Hynek Baran's avatar
Hynek Baran committed
1540
  
Hynek Baran's avatar
Hynek Baran committed
1541 1542
  if rt > 1 then report(lb,[cat(`total number: `,nops(zl),` c.c.`)]) fi;
  zl := [op({op(zl)} minus `cc/s`)];
Hynek Baran's avatar
Hynek Baran committed
1543
  inc(`cc/count/total`, nops(zl));
Hynek Baran's avatar
Hynek Baran committed
1544 1545
  if rt > 1 then report(lb,[cat(`of them new: `,nops(zl),` c.c.`)]) fi;
  zl := select(proc(z,m) evalb(op(6,z) = m) end,
Hynek Baran's avatar
Hynek Baran committed
1546
    zl, min(op(map(proc(z) op(6,z) end, zl)))); ## !
Hynek Baran's avatar
Hynek Baran committed
1547 1548 1549 1550 1551 1552 1553
#   zl := select(
#     proc(z,zl) local z1;
#       for z1 in zl do
#         if type (op(6,z1)/op(6,z),'count') then RETURN(true) fi
#       od;
#       false
#     end, zl, zl);
Hynek Baran's avatar
Hynek Baran committed
1554
   inc(`cc/count/computed`, nops(zl));
Hynek Baran's avatar
Hynek Baran committed
1555 1556 1557 1558 1559 1560 1561 1562 1563 1564
  if rt > 1 then report(lb,[cat(`of them minimal: `,nops(zl),` c.c.`)]) fi;
  if rt > 1 then report(lb,[`to be computed:`, nops(zl),` c.c.`]) fi; 
  for z in zl do 
    aux := `cc/1`(op(1..5,z)); 
    if aux = 0 then `cc/s` := `cc/s` union {z} 
    else `cc/tab`[f,op(z)] := aux;
      ans := ans, aux
    fi;
  od;
  if rt > 0 then report(lb,[`number of c.c.:`, nops({ans})]) fi;
Hynek Baran's avatar
Hynek Baran committed
1565
  inc(`cc/time`, time()-time0);
Hynek Baran's avatar
Hynek Baran committed
1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584
  {ans}
end:

`cc/s` := {}:

`cc/1` := proc(f,g,p,h,q)
  local ans,pt;
  global `cc/tab`;
  if assigned(`cc/tab`[f][g,p,h,q]) then
    Simpl(eval(`cc/tab`[f][g,p,h,q]));
  elif p = 1 then
    if q = 1 then Simpl(eval(g - h))
    else Simpl(eval(g - pd(h,q)))
    fi
  elif q = 1 then Simpl(eval(pd(g,p) - h))
  else Simpl(eval(pd(g,p) - pd(h,q)))
  fi
end:

Hynek Baran's avatar
Hynek Baran committed
1585 1586
# The common initialization

Hynek Baran's avatar
Hynek Baran committed
1587 1588
`clear/cc` := proc()
  global `cc/s`,`cc/tab`;
Hynek Baran's avatar
Hynek Baran committed
1589
  local T;
Hynek Baran's avatar
Hynek Baran committed
1590
  print(`Clearing cc's.`);
Hynek Baran's avatar
Hynek Baran committed
1591
  `cc/s` := {};
Hynek Baran's avatar
Hynek Baran committed
1592
  `cc/tab` := table([]);  
Hynek Baran's avatar
Hynek Baran committed
1593 1594
  CC:-init();
  return NULL